Global sensitivity analysis informed model reduction and selection applied to a Valsalva maneuver model
E. Benjamin Randall, Nicholas Z. Randolph, Alen Alexanderian, Mette S. Olufsen
GGlobal sensitivity analysis informed model reduction and selection appliedto a Valsalva maneuver model
E. Benjamin Randall a,b, ∗ , Nicholas Z. Randolph a , Alen Alexanderian a , Mette S. Olufsen a, ∗∗ a Department of Mathematics, North Carolina State University, Raleigh, NC b Department of Molecular and Integrative Physiology, University of Michigan, Ann Arbor, MI
Abstract
In this study, we develop a methodology for model reduction and selection informed by globalsensitivity analysis (GSA) methods. We apply these techniques to a control model taking systolicblood pressure and thoracic tissue pressure data as inputs, predicting heart rate in response to theValsalva maneuver (VM). The study compares four GSA methods based on Sobol’ indices (SIs)quantifying parameter influence on the difference between the model output and the heart ratedata. The GSA methods include standard scalar SIs determining the average parameter influenceover the time interval studied and three time-varying methods analyzing how parameter influencechanges over time. The time-varying methods include a new technique, termed limited-memorySIs, predicting parameter influence using a moving window approach. Using the limited-memorySIs, we perform model reduction and selection to analyze the necessity of modeling both the aorticand carotid baroreceptor regions in response to the VM. We compare the original model to threesystematically reduced models including (i) the aortic and carotid regions, (ii) the aortic regiononly, and (iii) the carotid region only. Model selection is done quantitatively using the Akaike andBayesian Information Criteria and qualitatively by comparing the neurological predictions. Resultsshow that it is necessary to incorporate both the aortic and carotid regions to model the VM.
Keywords:
Mathematical modeling, Sobol’ indices, Time-dependent processes, AkaikeInformation Criterion, Bayesian Information Criterion ∗ Corresponding author ∗∗ Principal corresponding author
Email addresses: [email protected] (E. Benjamin Randall), [email protected] (Nicholas Z. Randolph), [email protected] (Alen Alexanderian), [email protected] (Mette S. Olufsen)
Preprint submitted to May 27, 2020 a r X i v : . [ q - b i o . Q M ] M a y bbreviations In order of appearance:Systolic blood pressure SBPIntrathoracic pressure ITPValsalva maneuver VMLocal sensitivity analysis LSAGlobal sensitivity analysis GSAQuantity of interest QoISobol’ indicies SIsPointwise-in-time Sobol’ indices PTSIsGeneralized Sobol’ indices GSIsLimited-memory Sobol’ indices LMSIsAkaike information criterion with correction AICcBayesian information criterion BICDelay differential equations DDEElectrocardiogram ECGPiecewise cubic Hermite interpolating polynomial PCHIPRespiratory sinus arrhythmia RSAOrdinary differential equations ODEStandard deviation SDOne-at-a-time analysis OTAInitial conditions ICs2 . Introduction
Mathematical models describing cardiovascular processes are typically complex with many nonlin-ear interactions. The level of detail necessary to explain these interactions varies with the type ofsystem studied and the questions investigated. Such models with many interrelated states can bedifficult to analyze and parameter values difficult to assign, especially with limited data to informthe model. Simplifying or reducing such models may facilitate model-based data analysis providedthat the reduced model retains similar behavior to the original. Though formal model reductionmethods exist [3], they typically focus on analyzing model input-output relationships, ignoringpredicted quantities for which there is no available data. Moreover, these techniques identify onereduced model when in fact many reduced models can be generated. In this study, we develop asystematic approach for model reduction and selection via global sensitivity analysis. We applythis protocol to our model [35] that uses systolic blood pressure (SBP) and intrathoracic tissuepressure (ITP) data as inputs to predict heart rate responses to the Valsalva maneuver (VM). Weuse model selection to investigate if it is necessary to differentiate between the aortic and carotidbaroreceptor afferent signals in the model and determine which compartments should be modeled.The VM is conducted by forcefully exhaling against a resistance [46]. The VM stimulates thebaroreceptor reflex ( baroreflex ) triggered by the deformation of baroreceptors in the aortic archand carotid sinus. These receptors sense changes in blood pressure initiating a signal cascadetransmitted through parasympathetic (vagal) and sympathetic neurons controlling heart rate, vas-cular resistance, vascular compliance, and cardiac contractility [46]. In healthy adults, these signalsmaintain homeostatic heart rate and blood pressure; however, in disease they can become impaired.To understand how to treat patients with impaired autonomic function, it is essential to know howthe signals are impacted. To trust model predictions, it is prudent to use a well-validated model inwhich parameters can be estimated uniquely when comparing the model output to data.To effectively perform optimization, we estimate parameters that are most influential to themodel output. We employ local (LSA) or global (GSA) sensitivity analysis techniques determin-ing the influence of the model parameters on a quantity of interest (QoI). LSA (e.g., [7, 27, 30])compute partial derivatives of the QoI with respect to the parameters, perturbing them one at atime about a known “local” value, quantifying the effect of each parameter without accounting forparameter interactions. GSA methods [38] compute parameter influence by analyzing parametersand their interactions over the entire prescribed parameter space. Popular GSA methods includeSobol’ indices (SIs) [41] (used in this study), Morris screening [28, 29], generalized sensitivity func-tions [16], and moment-independent importance measures [14]. Morris screening is computationallyinexpensive, but results are first-order approximations that do not account for higher order inter-actions [28]. Generalized sensitivity functions characterize the sensitivity of a model to nominalparameter values and take into account parameter interactions; however, they make more stringentassumptions on local parameter identifiability [16]. Moment-independent importance measures fo-cus on constructing probability density functions for the QoI that are computationally expensivebut intractable for larger systems [14].We compute GSA using SIs [37, 41] due to their ease of calculation and wide usage [6, 18, 22, 44].SIs apportion relative contributions of the overall effect on the QoI to each parameter. They wereoriginally developed to analyze parameter influence in models with scalar outputs [41] (referredto here as scalar
SIs) but have recently been extended to analyze time-varying QoIs [1]. Oneway to compute time-varying SIs is to calculate the index at every time point [2, 18] (referred tohere as pointwise-in-time
Sobol’ indices (PTSIs)). This approach is simple but neglects the timecorrelation structure and history of the model output. To solve this problem, Alexanderian et al.[1] proposed a method incorporating time dependence, termed generalized
Sobol’ indices (GSIs),3hat have many applications for models with time-dependent outputs. However, for simulationsthat involve significant, short-lived disturbances to the steady-state behavior within the time courseof interest, integration over the entire time interval can average the effects of model parameters,diminishing features that only play a role during short time intervals. We illustrate this pointusing our VM model, where the breath hold causes a substantial yet transient model response. Tomitigate these effects, we propose a new method referred to here as limited-memory
Sobol’ indices(LMSIs). These indices account for parameter interactions, time correlations, and a snapshot ofthe history of the variance by applying weights from a moving window of width ∆ < T for T thelength of the time interval.In addition to sensitivity analysis, we propose a methodology for model analysis, reduction, andselection. For model reduction, we employ the LMSIs to select parameters that are noninfluentialover the entire time interval. These parameters are “removed” using one of two processes: (i)by excision of the mathematical equation(s) associated with the noninfluential parameter [7, 26];or (ii) by fixing the noninfluential parameter at its nominal value [41]. The former approachchanges the structure of the model. This process generates a set of reduced models, upon whichwe perform model selection using the Akaike Information Criterion with correction (AICc) andBayesian Information Criterion (BIC) [47]. Finally, we determine if the reduced models can predictheart rate and autonomic responses to the VM.In summary, new contributions of this article include: (i) the development of LMSIs for time-dependent processes incorporating a moving window to resolve the transient nature of the VM; and(ii) a model reduction and selection protocol informed by the LMSIs. These methods are illustratedusing our model predicting the heart rate response to VM [35]. Finally, we use model selection totest the need for including afferent signaling from the carotid region only, the aortic region only, orboth regions simultaneously when predicting heart rate dynamics.This study is organized as follows: Section 2 describes the VM; Section 3 presents the mathemat-ical methods, including the LMSIs, GSA-informed model reduction, and model selection; Section 4compiles the results of our analyses; Section 5 discusses key findings; and, finally, Section 6 statesour conclusions.
2. Valsalva maneuver
The VM is a clinical test involving forced expiration against a closed airway while maintainingan open glottis [11]. It stimulates the parasympathetic and sympathetic nervous systems viathe baroreflex in response to blood pressure changes in the aortic arch and the carotid sinus [4].In an effort to maintain homeostasis, the body responds to a sudden decrease in blood pressure,causing compensatory responses in heart rate, vascular resistance, vascular compliance, and cardiaccontractility, which in turn restore blood pressure to baseline. The VM is divided into four phases(illustrated in Figure 1):I. Initiation of the breath hold causes a sharp increase in blood pressure and a decrease in heartrate.II. Phase II is split into two parts:i. Early phase II - blood pressure drops significantly, triggering parasympathetic withdrawaland accelerating heart rate.ii. Late phase II - delayed sympathetic activation accelerates heart rate further and increasesperipheral vascular resistance, which in turn, increases blood pressure.III. Release of the breath hold results in a sharp decrease in blood pressure, triggering additionalparasympathetic withdrawal. 4
DataOptimized Model Predictions
I II III IV I II III IV E C G ( m V ) P ( mm H g ) E ff e r en t F i r i ng ( % ) H ( bp m ) P t h ( mm H g ) H ( bp m ) T p,b T s Data Model
Figure 1: A Valsalva maneuver (VM) for a representative control subject. The VM phases are represented asalternating gray (I and III) and light gray (II and IV) boxes. Phase II is divided into early and late parts (verticaldotted black line).
Data : (a) Electrocardiogram (ECG, mV). (b) Heart rate (H, bpm). (c) Blood pressure (P, mmHg)with systolic blood pressure (bold) indicated. (d) Thoracic pressure ( P th , mmHg) given in equation (2). Optimizedmodel predictions : (e) Efferent baroreflex-mediated parasympathetic ( T p,b , dimensionless, purple) and sympathetic( T s , dimensionless, green) signals. (f) Heart rate model output ( H , bpm, red) fitted to data (blue). odel with nominal parameter values ForwardModelAnalysis GSA ParameterEstimation ModelSelection
Original ModelModel Model M RankparameterinfluenceRecalibrate parameters ReduceModel
Model Reduction
Figure 2: Workflow diagram illustrating the steps in the procedure outlined in Section 3. From left to right: A modelis developed and nominal parameter values are determined a priori . Forward model analysis produces outputs,quantities of interest, which undergo global sensitivity analysis (GSA). The results from the GSA are used to reducethe model iteratively (orange arrow), producing M reduced models. For each model, a subset of parameters isestimated to fit data. The reduced model that captures the original model best, both qualitatively and quantitatively,is selected. Iterative Model Reduction Step (insert): Using GSA, parameters are ranked from the most to the leastinfluential, and parameters below a preset threshold are termed noninfluental. The noninfluential parameters are“removed”, i.e., they can be fixed at their nominal values or the model components associated with these parameterscan be analytically removed. Each iteration of this step produces a new reduced model. For each reduced model,nominal parameter values are recalibrated and estimated. Forward model evaluations are conducted and GSA isperformed again.
IV. Increased sympathetic activation causes blood pressure to overshoot and return to baselinewithin 30 s, while normalization of parasympathetic activity gives rise to a sharp drop in heartrate, which subsequently returns to baseline.
3. Methods
We apply the GSA techniques to analyze a system of delay differential equations (DDEs) with 6states and 25 parameters that takes electrocardiogram (ECG), SBP, and ITP data as inputs andpredicts heart rate as the output. Figure 2 depicts the model reduction and selection workflow,which includes the following components:1.
Forward model analysis computes the QoI by solving the system with nominal parametervalues. In this study, the QoI is the time-varying residual vector r ( t ; θ ) = (cid:20) H m ( t ; θ ) − H d ( t ) H d ( t ) , ..., H m ( t N ; θ ) − H d ( t N ) H d ( t N ) (cid:21) , (1)where H m ( t j ; θ ) denotes the heart rate model output at time t j for j = 1 , . . . , N and H d ( t j )denotes the corresponding observed heart rate data. Model predictions depend on the pa-rameter vector θ ∈ Ω p ⊆ R p , where Ω p is the prescribed parameter space of dimension p .2. Global sensitivity analysis determines the influence of the model parameters on the varianceof the QoI. We compute scalar SIs with respect to the Euclidean norm of r ( || r ( t ; θ ) || ) aswell as the time-varying PTSIs, GSIs, and LMSIs with respect to r .(a) Parameter influence is computed with the chosen index. We compute LMSIs and de-termine “noninfluential” parameters as those that have an index that does not exceed aprescribed threshold.(b)
Model reduction is performed by “removing” noninfluential parameters using one of twoapproaches: 6. by fixing noninfluential parameters at their nominal value; or,ii. by analytically removing equations affected by noninfluential parameters.(c)
Model recalibration is required if equations containing noninfluential parameters are mod-ified or removed. This ensures that the reduced model produces similar predictions tothe original model.(d)
GSA on the reduced model is conducted (again) to test that all parameter influenceindices are above the threshold, i.e., to test if all parameters are influential. This stepis necessary since model reduction can cause new parameters to be pushed below thenoninfluential threshold.Steps (b)-(d) are repeated iteratively until the system has no more noninfluential parameters.3.
Reduced models : We create a set M = { m , m , . . . , m M } , where m is the original model and m k for k = 1 , . . . , M are the reduced models.4. Parameter estimation : To determine a subset of identifiable parameters (ˆ θ ), we use subset se-lection analyzing the Fisher Information Matrix as suggested in [26, 29, 30, 33]. A limitationof this method is that it is local in nature, since the model is evaluated at specific parametervalues. Each of the reduced models in M are fitted to data estimating the identifiable influ-ential parameters by minimizing the least squares error between the model predictions anddata (the QoI).5. Model selection : To ensure that the reduced models produce outputs within physiologicallyacceptable ranges, we perform model selection both quantitatively and qualitatively. Quan-titatively, we select the reduced model that best fits the data (heart rate) by calculatingrelative AICc and BIC values. Qualitatively, we assume the original model is the true sig-nal and compare the outputs (parasympathetic and sympathetic activity) from the reducedmodels to the original.
The data used in this study is measured in a 21-year-old female healthy volunteer following aprotocol described in our previous work [35]. Figure 1 displays the collected ECG, blood pressure,and heart rate signals. To determine SBP (Figure 1c) and respiration, we interpolate the localmaxima of consecutive blood pressure waveforms in the pulse pressure data and consecutive QRS-complex amplitudes in the ECG data, respectively, via a piecewise cubic Hermite interpolatingpolynomial (PCHIP) in MATLAB R (cid:13) P th ) combines the measured ITP and the ECG-derived respiration, giving P th j = ITP j t s ≤ t j ≤ t e R M − R m ¯ R I − ¯ R E R j + ( R m − ¯ R E ) otherwise , (2)where R denotes the ECG-derived respiration, R M = 6 and R m = 3 . R I and ¯ R E are the mean values for the end of inspiration and expirationcalculated from the data, and t s and t e are the start and end times of the VM extracted directlyfrom the ITP data [35]. Heart rate is computed from the ECG signal using a built-in LabChart R (cid:13) function that detects RR intervals. Taking SBP and P th as inputs, the model predicts heart rate ( H ) as an output. Figure 3 displaysa schematic of the model. The upper section, denoted with solid black arrows, represents the7 edulla Output Inputs
Afferent
Efferent
SBP
Baroreflex-mediatedParasympatheticBaroreflex-mediatedSympathetic Heart RateRespiratory-mediatedParasympatheticCarotid StrainAortic Strain NeuralIntegration
Baroreflex MechanismRespiratory Sinus Arrhythmia
Figure 3: Model diagram reproduced from [35]. The model takes systolic blood pressure (SBP, mmHg, Figure 1c) andthe thoracic pressure ( P th , mmHg, Figure 1d, equation (2)) as inputs. The baroreflex model (solid arrows) predictsthe baroreceptor strains in the carotid sinus ( ε b,c , dimensionless, equation (5)) as a function of SBP and the aorticarch ( ε b,a , dimensionless, equation (5)) as a function of SBP and P th . These signals are integrated in the medulla( n , s − , equation (6)). A signaling cascade initiates a dynamic change in baroreceptor-mediated parasympathetic( T p,b , dimensionless, equation (8)) and sympathetic ( T s , dimensionless, equation (9)) pathways. The respiratory sinusarrhythmia model (dotted arrows) depends solely on P th , as the thorax oscillates with normal breathing. Fluctuationsin respiration modulate the respiratory-mediated parasympathetic outflow ( T p,r , dimensionless, equation (11)). Allefferent outflows are combined to control the heart rate ( H , bpm, equation (12)). baroreflex pathway, and the lower section denoted with dotted black arrows encodes heart ratecontrol via respiratory sinus arrhythmia (RSA).The model assumes that the pressure in the carotid sinus is the same as the measured SBP andthat the pressure in the aortic arch can be predicted from the SBP and P th by P c = SBP and P a = SBP − P th . (3)The baroreceptors are embedded in the arterial wall, which deforms nonlinearly as arterial pressureincreases [45]. Similarly to Mahdi et al. [25], we compute the strain of the arterial wall ε w,j for j = c or a for carotid or aortic, respectively, as ε w,j = 1 − (cid:115) e − q w ( P j − s w ) A + e − q w ( P j − s w ) , (4)where q w (mmHg − ) and s w (mmHg) are the steepness and half-saturation value of the sigmoidalrelationship and A (dimensionless) is an offset parameter. The strain of the baroreceptors as theydeform along with the arterial wall is predicted using a linear ordinary differential equation (ODE)of the form d ε b,j d t = − ε b,j + K b ε w,j τ b , (5)where K b (dimensionless) and τ b (s) are the gain and time-scale, respectively. The medulla receivessignals from the baroreceptors, which sense the relative strain from rest [4]. We compute anintegrated neural signal via a convex combination of the relative strains of both the aortic andcarotid regions as n = B ( ε w,c − ε b,c ) + (1 − B )( ε w,a − ε b,a ) , B ∈ [0 , . (6)8ote, if B = 0 (s − ), the model depends solely on the aortic baroreceptors, and if B = 1, it dependssolely on the carotid baroreceptors.After the medulla receives signals from the baroreceptors, it inhibits the parasympathetic andstimulates the sympathetic nervous systems. To model these interactions, we use the sigmoidalfunctions G p,b = 11 + e − q p,b ( n − s p,b ) and G s = 11 + e q s ( n − s s ) , (7)where q k (s) and s k (s − ) are the steepness and half-saturation values, respectively, for k = p, b or s denoting baroreflex-mediated parasympathetic and sympathetic, respectively. The efferent signals,transmitted from the medulla to the heart, are the solutions to the differential equationsd T p,b d t = − T p,b + K p,b G p,b τ p,b and (8)d T s d t = − T s ( t − D s ) + K s G s τ s , (9)where K k (dimensionless) and τ k (s) are the gains and time-scale, respectively, and D s (s) is thediscrete delay in the sympathetic outflow. Figure 1e displays the time series for T p,b and T s for thissubject with optimized parameter values.The respiratory-mediated parasympathetic outflow is modulated by the medulla as well. Thissignal has many inputs which are lumped into the thoracic pressure signal P th . The integration of P th into the medulla is determined by solving a decreasing sigmoid function G p,r = 11 + e q p,r ( P th − s p,r ) , (10)where q p,r (mmHg − ) and s p,r (mmHg) are the steepness and half-saturation value. The efferentsignal is modeled with the ODE d T p,r d t = − T p,r + K p,r G p,r τ p,r , (11)where K p,r (dimensionless) and τ p,r (s) denote the gain and time-scale.Inspired by previous work [23, 31, 32], the heart rate is the solution to the ODEd H d t = − H + ˜ Hτ H , (12)where τ H (s) is a time-scale and˜ H = H I (1 − H p,b T p,b + H p,r T p,r + H s T s ) . (13) H I (bpm) is the age-dependent intrinsic heart rate (when the sinoatrial node is denervated [15]),while H p,b (dimensionless), H p,r (dimensionless), and H s (dimensionless) scaling parameters. Figure1f shows the optimized model fit to the heart rate data.In summary, heart rate is predicted as the solution to a system of stiff ODEs and DDEs of theform d x d t = f ( t, x ( t ) , x ( t − D s ); θ ) , x ( t ) = x for t ∈ [ − D s , , (14)9here x = [ ε b,c , ε b,a , T p,b , T p,r , T s , H ] T ∈ R is the vector of states, f : R → R is the righthand side of the system, t denotes time in seconds, D s (s) is the discrete delay, and θ ∈ R is avector of parameters, including θ = [ A, B, K b , K p,b , K p,r , K s , τ b , τ p,b , τ p,r , τ s , τ H , q w , q p,b , q p,r , q s ,s w , s p,b , s p,r , s s , H I , H p,b , H p,r , H s , t s , t e ] T . (15)Table 1 lists the nominal parameter values. Calculations of the initial conditions for the states aregiven in Appendix A. The model is solved using the Fortran implementation of the stiff DDE solverRADAR5 [9] with relative and absolute tolerances of 10 − . The time-varying H model output wasinterpolated at the discretized time points from the heart rate data. As described in Appendix A,the model was fitted to the heart rate data by estimating a parameter subset that minimizes theleast squares error between the model output and the data.Upper and lower parameter bounds are given in Table 1. Parameters not calculated a priori are varied by ±
50% of their nominal values. Parameters calculated from the data have bounds setto the mean value ± τ s , which interactswith the delay D s and can cause instability [36]. To remain in the stable region, τ s is varied ±
50% of its nominal value. The parameter B denoting the convex combination of aortic and carotidbaroreceptors varies from 0 to 1. The GSA and the parameter estimation methods are performedon the logarithm of the parameters. Therefore, we enforce positivity by setting the lower boundsthat became negative to 0.01. We use variance-based SIs [37, 41] for the GSA. This section discusses four different methods forcomputing SIs: scalar [37], pointwise-in-time [2], generalized [1], and limited-memory (new). Areference table of symbols is compiled in Table 2.
Consider a mathematical model f with a scalar output y dependent on θ ∈ Ω p ⊆ R p , a vector of p uncertain model parameters with a prescribed parameter space Ω p ; that is, y = f ( θ ) . (16)For each parameter θ i , we compute its contribution to the variance of y [40, 41]. Assuming thatthe parameters are independent, the main effect on f by varying θ i is given by S i ( f ) = V θ i ( E θ ∼ i [ f | θ i ]) V ( f ) , (17)where V ( · ) and E [ · ] denote the variance and expectation operators and θ ∼ i is the vector θ withoutparameter θ i , that is θ ∼ i = [ θ , θ , . . . , θ i − , θ i +1 , . . . , θ p ] T . (18)The total-effect scalar SI on f includes both the main and higher order effects of θ i and is given by T i ( f ) = E θ ∼ i [ V θ i ( f | θ ∼ i )] V ( f ) = 1 − V θ ∼ i ( E θ i [ f | θ ∼ i ]) V ( f ) . (19)10 able 1: Parameters. Parameter Nom Physiological Lower Upper BoundDescription Units value Range bound bound Explanationand Symbol (Mean ± SD)
Offset A ± B s − , K b ± K p,b ± K p,r ± K s ± τ b s 0.9 0.45 1.8 Nom ± τ p,b s 1.8 6.5 ± ∗ ± τ p,r s 6 9.6 ± ∗ ± τ s s 10 14 ± ∗ ± τ H s 0.5 0.25 0.75 Nom ± q w mmHg − ± q p,b s 10 5 15 Nom ± q p,r mmHg − ± q s s 10 5 15 Nom ± s w mmHg 123 ± ∗∗
83 163 Mean ± s p,b s − ± ∗∗ ± s p,r mmHg 4.88 ± ∗∗ ± s s s − ± ∗∗ ± H I bpm 100 ± ∗∗
86 114 Mean ± H p,b ± ∗ ± H p,r ± ∗ ± H s ± ∗ ± D s s 3VM start/end t s s data t e s data Nom: nominal. SD: standard deviation. VM: Valsalva maneuver. Absence of units denotesdimensionless parameter. ∗ denotes parameter range is from optimized parameter valuesfrom [35]. ∗∗ denotes values calculated from data. able 2: Summary of Sobol’ index (SI) symbols. SI Type Symbol Reference
Main TotalScalar S i ( f ) T i ( f ) [37]Time-varying Pointwise-in-time (PTSIs) S i ( f ; t ) T i ( f ; t ) [2]Generalized (GSIs) S i ( f ; I t ) T i ( f ; I t ) [1]Limited-memory (LMSIs) S i ( f ; I ∆ ( t )) T i ( f ; I ∆ ( t )) This study i = 1 , . . . , p denotes the parameter index I t = [0 , t ] and I ∆ ( t ) = [ t − ∆ , t ] for integration window ∆. The following formulations attempt to account for changes in parameter influence over time.
Pointwise-in-time Sobol’ indices (PTSIs) : Consider a model f with time-varying output y ( t )on an interval I T = [0 , T ] for end time point T >
0, i.e., y ( t ) = f ( t ; θ ) , t ∈ I T . (20)For this model, the main effect PTSI of f by varying θ i at time t is given by S i ( f ; t ) = V θ i ( E θ ∼ i [ f ( t ; · ) | θ i ]) V ( f ( t ; · )) , t ∈ I T , (21)and the total effect PTSI of f by varying θ i at time t is given by T i ( f ; t ) = E θ ∼ i [ V θ i ( f ( t ; · ) | θ ∼ i )] V ( f ( t ; · )) , t ∈ I T . (22) Generalized Sobol’ indices (GSIs) : The indices in equations (21) and (22) ignore time correla-tions and the history of the time-dependent signal. The method proposed by Alexanderian et al. [1]takes into account time correlations by integrating the numerators and denominators of equations(21) and (22) over time. Let I t = [0 , t ] be a time-varying interval. The the main effect GSI of f byvarying θ i over the interval I t is given by S i ( f ; I t ) = (cid:90) t V θ i ( E θ ∼ i [ f ( τ ; · ) | θ i ]) d τ (cid:90) t V ( f ( τ ; · )) d τ , t ∈ I T . (23)Similarly, the total effect GSI on f by varying θ i over the interval I t is given by T i ( f ; I t ) = (cid:90) t E θ ∼ i [ V θ i ( f ( τ ; · ) | θ ∼ i )] d τ (cid:90) t V ( f ( τ ; · )) d τ , t ∈ I T . (24)The integrals in equations (23) and (24) are computed following the method in Appendix B. Limited-memory Sobol’ indices (LMSIs) : By integrating over I t , equations (23) and (24)compute the parameter influence up to time t , which can be interpreted as an average across I t . However, this method might miss the contribution of the parameters during fast, transient12isturbances, especially for problems where an extended baseline is necessary to obtain steady-stateconditions. To remedy this, we propose a new technique, limited-memory Sobol’ indices, to analyzeparameter influence as the QoI responds to transient disturbances in its steady-state behavior. Todo so, we introduce a moving integration window I ∆ ( t ) of width ∆. The shape of the window andmagnitude of ∆ is necessarily problem-dependent. By implementing this window, the integrationinterval I T necessarily decreases by ∆; that is, the new integration interval is I T − ∆ = [∆ , T ]. Themain effect LMSI on f by varying θ i over the interval I ∆ ( t ) is given by S i ( f ; I ∆ ( t )) = (cid:90) tt − ∆ V θ i ( E θ ∼ i [ f ( τ ; · ) | θ i ]) w ( τ ) d τ (cid:90) tt − ∆ V ( f ( τ ; · )) w ( τ ) d τ , t ∈ I T − ∆ , (25)where w ( t ) denotes the weights at each time point t j determined by the window of choice. Similarly,the total effect LMSI on f by varying θ i over the interval I ∆ ( t ) is given by T i ( f ; I ∆ ( t )) = (cid:90) tt − ∆ E θ ∼ i [ V θ i ( f ( τ ; · ) | θ ∼ i )] w ( τ ) d τ (cid:90) tt − ∆ V ( f ( τ ; · )) w ( τ ) d τ , t ∈ I T − ∆ . (26)By integrating the variances within a moving window, we can capture the transient changes inparameter influence (as in the PTSIs) but avoid missing short-lived responses (as in the GSIs).There are many factors contributing to choosing the appropriate type, shape, and width ofthe window. In this study, we let the physiology guide window selection. Since both the bloodpressure and the change in blood pressure affect the baroreflex modulation of heart rate, recentblood pressure and heart rate values contribute more to the current heart rate. Moreover, futuretime points have no bearing on the current heart rate. To accommodate these features, we used aGaussian trailing window. Since the VM is typically 15 s in length, a window of about 15 s shouldsuffice. We chose to use a trailing half-Gaussian integration window for ∆ = 15 s before time point t j to coincide with the typical length of the VM. Further discussion of the moving window can befound in Section 5.3. The GSA was conducted using Monte Carlo integration including all model parameters ( p = 23)except the Valsalva start and end times, t s and t e , which are extracted from the data. We computed L = 10 ( p + 2) resulting in 25,000 function evaluations [12, 22] and tested for convergence using L = 10 ( p + 2) with 250,000 function evaluations and L = 10 ( p + 2) with 2,500,000 functionevaluations, which produced similar results. For the reduced models, we use L evaluations, sincethe results did not change with larger sample sizes. To compute the variance and expectationgiven in equations (21)-(26), we used the method proposed by Saltelli et al. [37] (described indetail in Appendix B). To approximate the integrals in equations (23)-(26), we used the trapezoidquadrature scheme. We compare the performance of each of the four methods discussed in Section 3.3. The scalar SIsare calculated with respect to the Euclidean norm of the residual r , that is, S i ( || r || ) and T i ( || r || ),where r is given in equation (1). Using || r || as the scalar model output gives a decent indication of13he sensitivity of r to the parameters at steady-state, since the Euclidean norm can be consideredan average of the signal over time. However, this disregards the changes in parameter influenceas fast disturbances occur in the model output. Hence, in response to the time-varying r , wesimultaneously compute the PTSIs, GSIs, and LMSIs. To identify parameters that do not influence the model output, we rank the total effect scalar SIs T i ( || r || ) and group them into three categories: most, moderately, and least influential. η = 10 − is the threshold for moderately influential parameters assigned based on the clear separation inparameter influence values shown in Figure 4. η = 10 − is assigned in accordance with theLSA conducted in our previous work [35]. The parameter groups designated by the scalar indicesmotivated the grouping used for the time-varying indices (PTSIs, GSIs, and LMSIs). We define aparameter θ i as noninfluential over the time series if T i ( r ; t ) < η ∀ t ∈ I T . (27)Parameters that satisfy this criterion will either be fixed at their nominal value or removed fromthe analysis by analytically excising the model component in which they appear. Using the sensitivity ranking, we develop steps for an iterative model reduction methodology in-formed by the GSA. We use the results from the LMSIs to identify components of the model toremove. The steps are as follows:1. Compute the time-varying indices T i for all parameters to be considered.2. Determine if each T i satisfies the criterion in equation (27) for all time points. If so, this setof noninfluential parameters, θ NI , is flagged for removal.3. Analyze the parameters in θ NI and determine if it is possible to remove the equations asso-ciated with each parameter. This step is problem-specific and inherently changes the mathe-matical structure of the model. Choose a parameter θ k that has the least influence. Removethe components of the model associated with that parameter and restructure the model. Wewould like to stress that this is an iterative and model-specific process, as there are manyinstances where removing model components can be detrimental. Some changes to the modelequations must occur simultaneously, which we will exemplify in the next section.4. Recalibrate parameters to obtain new parameter vector θ (cid:63) for newly generated model f (cid:63) andrun the model. f (cid:63) joins the set of models, M .5. Repeat GSA and then steps 2-4 on f (cid:63) . Iterate until the parameters in θ NI are not algebraicallyremovable. Fix all remaining parameters in θ NI at their nominal values.This process generates a set of models M = { m , m , . . . , m M } for m the full model and M thenumber of m i reduced models. To compare model performance between the full and reduced models, we use a statistical measurethat calculates a trade off between how well the model fits the data (goodness of fit) with howcomplex the model is (number of estimated parameters) [34]. To quantify this comparison, weuse the AICc and BIC to calculate a regression between the model output and the data [47]. Weassume that the residual errors are independent and identically distributed with mean zero and14nite variance. By predicting the maximum likelihood estimate, which is equivalent to minimizingthe least squares error J (equation (A.9)), we computeAICc = N log (cid:16) JN (cid:17) + 2(ˆ p + 2) (cid:16) NN − (ˆ p + 2) − (cid:17) and (28)BIC = N log (cid:16) JN (cid:17) + (ˆ p + 2) log( N ) , (29)where N is the number of data points and ˆ p the length of the optimized parameter subset vectorˆ θ [5]. To compare the models, we report the relative index of each model from the minimal model,that is AICc rel = e (AICc min − AICc i ) / and (30)BIC rel = e (BIC min − BIC i ) / . (31)This statistical technique is useful when determining the goodness of fit to data. However, there areother predicted model outputs, such as T p,b and T s , that are of clinical importance since they cannotbe measured without costly and invasive procedures that blunt the signals with anesthetization.Therefore, we must also assess the model performance qualitatively. We do so by comparing thereduced model to the full model predictions, assuming that the full model produces the true signal.Since the traces for T p,b are very similar, we employ the metric Q = | max( T s,m ) − max( T s,m i ) | (32)determining the reduced model that minimizes the distance between the maxima of their respective T s predictions.In summary, these metrics account for both the model fits to data and the predicted quantitiesfor which there is no data. Both aspects are crucial for an effective reduced model. The data andrun-time environment for the model code can be found at DOI: 10.5281/zenodo.3856447.
4. Results
This section discusses the outcomes from the scalar Sobol’ indices (SIs) (Figure 4) and the time-varying SIs: pointwise-in-time (PTSIs), generalized (GSIs), and limited-memory (LMSIs) (Figure6). The latter is computed using a moving integration window of width ∆. We also compare severalwindow widths (Figure 8). Using outcomes from the LMSIs, we inform a model reduction. Lastly,we select the best model by computing the relative Akaike Information Criterion with correction(AICc rel , equation (30)) and relative Bayesian Information Criterion (BIC rel , equation (31)) inTable 5 and by examining the predicted model outputs T p,b and T s to give physiological predictions(Figure 7). The scalar SIs are shown in Figure 4 and the time-varying SIs are shown in Figure 6. For eachanalysis, the parameters are divided into three influence groups: most ( > η ), moderately (between η and η ), and least ( < η ) influential. 15 arametersRelative Sensitivity Ranking Figure 4: Relative scalar sensitivity analysis ranking with respect to || r || . Influence thresholds η = 10 − and η = 10 − are indicated with horizontal dashed lines. Global sensitivity analysis results (blue bars) are computedusing the scalar total effect Sobol’ indices ( T i ( || r || )) given in equation (19), scaled by the maximum sensitivity valueand plotted on a logarithmic scale for the y -axis. Local sensitivity analysis results (black x’s) are reproduced from[35]. Figure 4 shows the ranking of the total effect scalar SIs, ranked with respect to the Euclidean normof the model output residual || r || , that is, T i ( || r || ). These are plotted with the local sensitivityanalysis (LSA) results reproduced from Randall et al. [35]. The five most influential parameters,i.e., { H s , K s , K b , H p,r , K p,r } > η , are associated with RSA and the sympathetic tone during andafter the VM. These parameters influence the controllers of heart rate both at rest (during thebaseline) and during the transient impact of the VM. In comparison, the LSA determined that theparameters s p,b , s w , A , H I were the most influential. Note that the LSA is evaluated at particularinstant within the parameter space Ω p and each of the parameters is changed one at a time, whichcould lead to a different parameter ranking.The subset of least influential parameters, i.e., T i ( || r || ) < η , is θ NI,sca = { q w , s s , B, s p,b , τ s , τ H , τ p,b , τ b } . (33)Figure 5 displaying 10 samples from the least sensitive parameters varied between their upper andlower bounds (given in Table 1) provides useful insight into how each parameter affects the output.Parameters s s , s p,b , τ H , and τ b have the least effect on the model output. The small effect ofchanging parameters s s and s p,b is most likely due to their narrow distributions. In Randall et al.[35], these parameters are calculated a priori and do not vary much among the subjects studied.Parameters τ H and τ b have the least effect on the residual across all analyses. Figure 6 shows results for the most (Figure 6a), moderately (Figure 6d), and least (Figure 6g)influential parameters for the PTSIs. In the calculation of the PTSIs, the relative importance ofeach parameter is determined with respect to each time point. Comparisons between parameterinfluence measures at different time points can be difficult to interpret within one influence trace.The total effect SIs ( T i ( r ; t )) display rapid fluctuations. Since our goal is to determine regionswhere the influence of certain parameters increase or decrease substantially, the highly oscillatorynature of these solutions makes it difficult to compare the effect of each parameter, especially forthe the least influential parameters. 16 i m e ( s ) H ( bp m ) H ( bp m ) I II III IV I II III IV I II III IV I II III IV
Figure 5: Noninfluential parameters as determined by the total effect scalar Sobol’ indices in Figure 4. The modelwas evaluated within the physiological range of each parameter given in Table 1 (cyan) plotted with the heart ratedata (blue). The VM phases are represented as alternating gray (I and III) and light gray (II and IV) boxes. PhaseII is divided into early and late parts (vertical dotted black line).
The most influential parameters are H s , K s , K b , H p,r , and K p,r (Figure 6a). Our results showthat the sensitivity for K b does not change significantly during the time interval, which is expectedsince the baroreceptors are always active. Parameters K s and H s are most influential, due to theactivation of the sympathetic response after the onset of the VM. Finally, H p,r and K p,r decreasein influence below η after the onset of the VM, due to parasympathetic withdrawal that occurs inphase II of the VM.Parameters that remain under the threshold η at all time include s s , s p,b , and τ b . The half-saturation values are most likely small due to their narrow parameter distributions. It is surprising,though, that τ b remains below η , since all the other time-scales rise above η at least one timepoint. These insensitive parameters are prime candidates for removal via the methods discussed inSection 3.The pointwise-in-time technique designates almost all parameters as at least moderately influ-ential at one point in time, which coincides with the results from the LSA. The oscillatory natureof the results complicates the determination of parameter influence level at different regions in thetime series [1]. The PTSIs show that some of the parameters in θ NI,sca become moderately influ-ential over the time series, namely q w , B , τ s , τ H , and τ p,b . However, these results do not providemuch additional information that is not already given in both the scalar SI analysis and the LSA. As shown in Figures 6b, e, and h, the differences between the PTSIs and GSIs are evident. Theformer provides one case where the parameter influence oscillates significantly both at rest andduring dynamic changes. It does not take into account the history of the variance, and eventhe baseline results are difficult to interpret. The latter places too much emphasis on the time-dependence, and hence, averages the signal over extended periods of baseline activity, missing any ofthe potential transient changes in parameter influence on the model output.These methods provideextremes for the analysis of these time-varying signals, as shown in Figure 8b where an examplePTSI is in blue and GSI in red. These results suggest that there is a need for a method that17
Pointwise-in-time Generalized Limited-Memory
Most influentialModerately influentialLeast influential (a) (b) (c)(d) (e) (f)(g) (h) (i)
I II III IV I II III IV I II III IV T o t a l t i m e - v a r y i ng S obo l i nd i c e s Figure 6: Time-varying total effect Sobol’ indices (SIs) for the most (top row), moderately (middle row), and least(bottom row) influential parameters as determined by the scalar SIs in Figure 4 for the pointwise-in-time (a, d,g), generalized (b, e, h), and limited-memory (c,f,i) SIs. Thresholds η = 10 − and η = 10 − are indicated withhorizontal dashed black lines. The VM phases are represented as alternating gray (I and III) and light gray (II andIV) boxes. Phase II is divided into early and late parts (vertical dotted black line). incorporates both the transient nature of the PTSIs and the smoothing capabilities of the GSIs.As shown in Figure 6h, the least influential subset of parameters determined by the GSIs such that T ( r ; I t ) < η for all time t is θ NI,G = { q w , s s , B, s p,b , τ s , τ p,b , τ b } . (34)Parameters τ p,r (Figure 6e) and τ H (Figure 6h) are notable exceptions to parameters withrelatively constant traces. τ p,r clearly increases in influence after the onset of the VM, which isto be expected as the influence of RSA decreases substantially during the VM. The time-scale forheart rate, τ H , is moderately influential at the beginning of the time series but decreases in influenceafter the onset of the VM. This is surprising and may imply that τ H plays an important role inestablishing steady-state behavior from the initial conditions and then decreases in influence oncesteady-state is achieved. Fluctuations in parameter influence correspond to different control mechanisms that activate anddeactivate during the VM. We expect to observe the change in importance of parameters associated18 able 3: Noninfluential parameters from Sobol’ index methods.
Parameter Sobol’ indicesScalar Generalized Limited-memory θ NI,sca (33) θ NI,G (34) θ NI,LM (35) q p,r (cid:88) q w (cid:88) (cid:88) s s (cid:88) (cid:88) (cid:88) B (cid:88) (cid:88) (cid:88) s p,b (cid:88) (cid:88) (cid:88) τ s (cid:88) (cid:88) (cid:88) τ H (cid:88) τ p,b (cid:88) (cid:88) τ b (cid:88) (cid:88) (cid:88) with activated control mechanisms that affect the heart rate dynamically. As shown in Figures 6c,f, and i, the LMSIs remain relatively constant before the VM. The LMSIs plotted in the figurehave a trailing half-Gaussian moving window with ∆ = 15 s. It is important to choose a windowthat is short enough such that it can retain important features in the signal but long enough to beeasily interpretable. Since the VM typically lasts 15 s, we compare the output from several windowwidths for ∆ = 5 , , , and 20 s (Figure 8d). For smaller values of ∆, the LMSI converges to thePTSI, i.e., as ∆ → S ( f ; I ∆ ( t )) → S i ( t ). For larger values of ∆, the LMSI converges to the GSI,i.e., as ∆ → T , S i ( f ; I ∆ ( t )) → S i ( f ; I t ). Excluding ∆ = 5 and 20, the traces for ∆ = 10 and 15 aresimilar. To remain consistent with the VM itself, we chose a window of 15 s.The LMSIs provide a balance between the PTSIs and GSIs by smoothing out the highly oscil-latory PTSIs and incorporating the time-correlation structure preserved by the GSIs, if only fora small window of time. The LMSIs retain some modulation due to the transient nature of theVM, providing distinct changes in parameter influence rankings before, during, and after the VM.Therefore, we conclude that if the QoI is operating in steady-state, computing the scalar SIs wouldsuffice for the analysis. However, during a transient disturbance from baseline, LMSIs provide amore informative parameter ranking.The LMSIs show clear changes in parameter influence as time evolves. In particular, we observeincreases in influence of parameters associated with T s during the VM that subsequently decreasetowards the end of phase IV. This behavior coincides with the activation of sympathetic tone duringthe VM and then deactivation after the release of the breath hold.For the least influential parameters, q p,r , q w , τ H , and τ p,b become moderately influential withinthe interval I T , while the subset θ NI,LM = { s s , B, s p,b , τ s , τ b } (35)remains below η . The half-saturation values s p,b and s s are also in this subset due to their narrowdistribution. τ b has consistently been the least influential throughout all of the analyses and isexpected to be in this subset. It is surprising that B and τ s remain below η for all t , since B determines the the contribution of the aortic versus carotid baroreceptors. The necessity of thisparameter has been supported in previous studies [20, 35], but the aim of the model reduction inthe next section explores whether it is structurally necessary to include both baroreceptor regions.The sympathetic time-scale τ s has been found in our previous work to be nonlinearly related tothe delay D s [36]. In this study, we have held D s constant at 3 s and have only varied τ s withinparameter bounds that maintain stable model behavior (Table 1). Though our previous work showsthe importance of τ s , it may not be as influential to the heart rate residual in this range at D s = 3.In summary, the scalar SIs provide a clear parameter influence ranking that is similar to theranking produced by the LSA. The time-varying techniques elucidate when a particular parameter19 able 4: Estimated parameter values. Symbol Model m m m m B † † τ p,b τ p,r H p,b H p,r H s † The B parameter was held constant atthis value and not included as a part ofthe subset. becomes important. Our results show that the GSIs have a clear benefit over the PTSIs, as theyincorporate the time correlation structure of the signal. However, GSIs miss the effects of transientdisturbances to steady-state behavior. Hence, we conclude that the LMSIs proposed in this studyillustrate parameter influence traces that correspond to the modeled physiological phenomena.Table 3 summarizes the subsets of least parameters from the scalar SIs (equation (33)), GSIs(equation (34)), and LMSIs (equation (35)). The PTSIs were not included as their results wereinconclusive. The LMSIs determined the smallest subset of noninfluential parameters for the entiretime interval with 5, which are considered for removal in the model reduction in Section 4.2. As stated in the protocol given in the Methods section (Section 3), we can “remove” parametersfrom consideration in two ways: (i) fix the parameter at its nominal value and rerun the analysisor (ii) analytically excise equations associated with the parameter in this subset. We first considerthe parameters in this subset that we will fix at their nominal values based on the model structure.Secondly, we remove equations associated with noninfluential parameters and develop a suite ofmodels upon which we can perform statistics.
Parameters s p,b and s s are the half-saturation values appearing in equation (7). Their nominal val-ues are computed a priori based on the assumption that 80% of the baseline heart rate is controlledby the baroreflex-mediated parasympathetic tone and 20% the baroreflex-mediated sympathetictone [19, 35]. The coefficient of variation (SD/mean) for s p,b is 0.1%, which implies a very narrowdispersion.As shown in our previous work [36], τ s , appearing in equation (9), is nonlinearly related tothe delay D s and this relationship can lead to instability. Therefore, we restricted the upper andlower bounds of τ s in this study to ±
50% rather than use the entire physiological range (mean ± p did not intersect an oscillatory regime.Therefore, we choose to fix τ s at its nominal value.After fixing these three parameters, we rerun the GSA calculating LMSIs. This analysis pro-duced a similar plot to Figure 6i and τ b and B are still in θ NI,LM (results not shown). These resultswere computed with 25,000 samples, mimicking the strategy used with the original model. Sinceholding these parameters fixed produced very similar results, we will consider this analogous to theoriginal model and refer to it as model m . From θ NI,LM , the time-scale for the baroreceptor strain τ b and the convex combination parameter B are flagged for removal by changing the equations associated with these parameters.20 emoving τ b : The parameter τ b appears in equation (5), the baroreceptor strains for the carotidand aortic baroreceptors. Since this value can be small and not impact the result, we can rearrangeequation (5) as d ε b,j d t = − ε b,j + K b ε w,j τ b ⇒ τ b d ε b,j d t + ε b,j = K b ε w,j . (36)By letting τ b = 0, equation (5) is replaced with the algebraic relation ε ∗ b,j = K b ε w,j . (37)From here on, we will let the asterisk denote the new equations of the reduced system determinedby the removal of τ b . With this substitution, we have reduced the number of ODEs by two.Substituting equation (37) into equation (6) gives n = B ( ε w,c − ε ∗ b,c ) + (1 − B )( ε w,a − ε ∗ b,a )= B ( ε w,c − K b ε w,c ) + (1 − B )( ε w,a − K b ε w,a )= (1 − K b ) (cid:16) Bε w,c + (1 − B ) ε w,a (cid:17) = (1 − K b ) n ∗ (38)and we define n ∗ = Bε w,c + (1 − B ) ε w,a . (39)By substituting equation (38) into equation (7), we obtain G ∗ p,b = 11 + e − q p,b ( n − s p,b ) (40)= 11 + e − q p,b ((1 − K b ) n ∗ − s p,b ) (41)= 11 + e − q ∗ p,b ( n ∗ − s ∗ p,b ) , (42)and we define new parameter values q ∗ p,b = q p,b (1 − K b ) and (43) s ∗ p,b = n ∗ + 1 q ∗ p,b ln (cid:16) K p,b ¯ T p,b − (cid:17) , (44)where ¯ T p,b is the steady-state value. New parameters q ∗ s and s ∗ s are determined similarly. The new21ystem of equations without the baroreceptor strain differential equations is ε w,j = 1 − (cid:115) e − q w ( ¯ P j − s w ) A + e − q w ( ¯ P j − s w ) for j = c or a, (45) n ∗ = Bε w,c + (1 − B ) ε w,a , (46)d T p,b d t = 1 τ p,b ( − T p,b + K p,b G p,b ) , G p,b = 11 + e − q ∗ p,b ( n ∗ − s ∗ p,b ) (47)d T s d t = 1 τ s ( − T s ( t − D s ) + K s G s ) , G s = 11 + e q ∗ s ( n ∗ − s ∗ s ) (48)d T p,r d t = 1 τ p,r ( − T p,r + K p,r G p,r ) , G p,r = 11 + e − q p,r ( P th − s p,r ) (49)d H d t = 1 τ H ( − H + ˜ H ) , ˜ H = H I (1 − H p,b T p,b + H p,r T p,r + H s T s ) . (50)This new system of ODEs and DDEs consists of 4 states and 23 parameters. Let · ∗ denote thereduced system d x ∗ d t = f ∗ ( t, x ∗ ( t ) , x ∗ ( t − D s ); θ ∗ ) , (51)where f ∗ is given by equations (45)–(50), x ∗ = [ T p,b , T p,r , T s , H ] T ∈ R , D s is the delay, and θ ∗ = [ A, B, K p,b , K p,r , K s , τ p,b , τ p,r , τ s , τ H , q w , q ∗ p,b , q p,r , q ∗ s ,s w , s ∗ p,b , s p,r , s ∗ s , H I , H p,b , H p,r , H s , t s , t e ] T ∈ R . (52)We refer to this reduced model as m .To determine a subset of parameters to optimize for m , we perform subset selection using thestructured correlation analysis described in Appendix A and obtainˆ θ m = [ B, τ p,b , τ p,r , H p,b , H p,r , H s ] T . (53)Table 4 summarizes the optimized parameters for m and m . We perform GSA on m , holdingparameters s p,b , s s , and τ s fixed and observe that B is the only remaining parameter below η forall time t (results not shown). Removing B : The convex combination parameter B in equation (39) determines the relativecontribution of the afferent signals from the aortic and carotid baroreceptors. The necessity ofdelineating these baroreceptor regions has been explored in two previous studies [20, 35] and it issurprising that this parameter has been flagged for removal using the protocol given in the Methodssection (Section 3). It is unknown how the aortic and carotid signals are integrated in the medullaand not clear whether it is sufficient to model only one of these regions. From Figure 6i, we observethat the influence of B is close to zero, but as the VM progresses, B increases, corresponding tothe stimulation of the baroreceptors. After the breath hold ends, B returns to zero. To remove theequation associated with B , we observe in equation (39) that B vanishes when it is zero, indicatingthat the aortic baroreceptor strain solely influences the efferent response, or one, indicating carotideffects. This produces the following models: m where B = 0 and n ∗ = ε w,a , (54)producing a system of ODEs and DDEs with 4 states and 22 parametersd x ∗ d t = f ∗ ( t, x ∗ ( t ) , x ∗ ( t − D s ); θ ∗ ) , (55)22 HR ( bp m )
012 012 T pb T s m m m m (a) (b) (c) I II III IV I II III IV I II III IV
Figure 7: Plots of full ( m , solid) and reduced models. Reduced model 1 ( m , dotted) - removed 2 states and twoparameters with parameter B estimated. Reduced model 2 ( m , dashed) - same as m except with B = 0. Reducedmodel 3 ( m , dash-dotted) - same as m except with B = 1. The Valsalva maneuver (VM) phases are represented asalternating gray (I and III) and light gray (II and IV) boxes. Phase II is divided into early and late parts (verticaldotted black line). (a) Model fits (red) to the heart rate data (blue). Insert shows a zoom of the heart rate fit in latephase II and phase III of the VM. (b) Model predictions of the efferent baroreflex-mediated parasympathetic ( T p,b ,purple) for the full and reduced models. (c) Model predictions of the efferent baroreflex-mediated sympathetic ( T s ,green) for the full and reduced models. where f ∗ is the right hand side, x ∗ = [ T p,b , T p,r , T s , H ] T ∈ R and θ ∗ = [ A, K p,b , K p,r , K s , τ p,b , τ p,r , τ s , τ H , q w , q ∗ p,b , q p,r , q ∗ s ,s w , s ∗ p,b , s p,r , s ∗ s , H I , H p,b , H p,r , H s , t s , t e ] T ∈ R ; (56)and m where B = 1 and n ∗ = ε w,c (57)producing a system of ODEs and DDEs with 4 states and 22 parametersd x ∗ d t = f ∗ ( t, x ∗ ( t ) , x ∗ ( t − D s ); θ ∗ ) , (58)where f ∗ is the right hand side, x ∗ = x ∗ and θ ∗ = θ ∗ . Using the structured correlation analysismethod given in Appendix A, we obtainˆ θ m = ˆ θ m = [ τ p,b , τ p,r , H p,b , H p,r , H s ] T . (59)Table 4 lists the estimated parameters for models m and m . We perform GSA on each model,which results in all parameters above threshold η at at least one time point within the interval I T (results not shown). Hence, there are no parameters remaining in the noninfluential subset.In summary, the model reduction methodology has produced 4 possible models, summarized inTable 4: m where parameter s p,b , s s , and τ s are held constant; m where the equations associatedwith τ b are removed; m where τ b is removed and B = 0; and m where τ b is removed and B = 1.Figure 7 plots each of the model fits to the data along with their respective T p,b and T s traces. To determine which model best fit the data, we perform model selection protocol both quantitativelyand qualitatively. The quantitative approach involves computing the relative AICc and BIC givenin equations (30) and (31), respectively, assessing the model fit to the data, while the qualitativeanalysis investigates the importance for prediction of T p,b and T s as they evolve in time. The latter23 able 5: Statistical analysis for model selection. Model Cost Parameters Quantitative Qualitative J (10 − ) p AICc rel
BIC rel Q (A.9) (30) (31) (32) m m m m AICc - Akaike Information Criterion with correction in equation (28).BIC - Bayesian Information Criterion in equation (29) . is important, since the goal of our model is not only to fit the data but also to predict the neuraltones that cannot be measured in vivo without blunting the parasympathetic and sympatheticresponses with anesthetization. Hence, we test the reduced models by comparing the T p,b and T s predictions to those obtained with the original model m . We computed the relative AICc and BIC values, AICc rel and BIC rel , respectively, for each of themodels, assuming that the residuals are independent and identically distributed. From Table 5, weobserve that m has the greatest AICc rel and BIC rel values. This is surprising, since m estimatesmore model parameters than both m and m , which typically yield greater AICc rel and BIC rel values [47]. Therefore, this analysis suggests that modeling both the aortic and carotid baroreceptorregions is necessary to predict heart rate. Figure 7 shows the H , T p,b , and T s traces for all of the models considered: m (solid curve), m (dotted curve), m (dashed curve) and m (dash-dotted curve). The model fits to heart rate data(Figure 7a) are all similar, which is expected since the model is calibrated to this representativedata set. The least squares costs of the fits are of the same magnitude (Table 5). Interestingly, m is the only model that is able to fit phase III of the data (Figure 7a insert), suggesting from aqualitative standpoint that m captures the most features of the heart rate data. Figure 7b displaysthe predicted T p,b trace for the models, which are all similar and difficult to compare. However,Figure 7c shows the trajectories for T s , exhibiting the greatest deviation from m . Since traces forheart rate and T p,b for all four models are similar, we use T s to compare the reduced models to thefull model via the Q value (equation (32)). Table 5 compiles these metrics and m has the lowest Q value. This result gives credence to the assertion that both the aortic and carotid baroreceptorsmust be included to predict heart rate and the neurological signals.In summary, we conclude that m , the reduced model after removing the equations associatedwith τ b , is the best model for the biological questions investigated here. The AICc rel and BIC rel values are the greatest while the qualitative metric Q was the lowest. This model predicts theheart rate the “best” with respect to complexity and number of parameters estimated, and thesympathetic output is closest to the original model.
5. Discussion
This study performs model reduction and selection using global sensitivity analysis (GSA) on acardiovascular model predicting parasympathetic ( T p,b ) and sympathetic ( T s ) nervous outflow andheart rate ( H ) in response to the Valsalva maneuver (VM). For the GSA, we used Sobol’ indices(SIs), conducting our analysis with respect to a scalar quantity of interest (QoI), computed using24he Euclidean norm of the heart rate residual || r || , and a time-varying QoI, the vector r given inequation (1). Computation for the scalar SIs [37] was done with respect to || r || while the time-varying pointwise-in-time (PTSIs) [2], generalized (GSIs) [1], and limited-memory (LMSIs) Sobol’indices are with respect to r . A novel component is the introduction of the LMSIs using a movingintegration window ∆ motivated by the transient VM process. The scalar SIs determined a rankingof parameter influence on the model output averaged over the entire time interval. We were ableto categorize parameters based on their influence on the model output into three groups: most,moderately, and least influential. Additionally, the LMSIs informed a model reduction protocolthat generated four models, m , m , m , and m . We analyzed the performance of these modelsqualitatively, comparing the model predicted signals to the original model m , and quantitatively,calculating the relative Akaike Information Criterion with correction (AICc rel , equation (30)) andthe Bayesian Information Criterion (BIC rel , equation (31)) for each model. m proved to be thebest model, implying the necessity of modeling both the aortic and carotid baroreceptor regions topredict heart rate accurately. This study focused on the use of a GSA to analyze the sensitivity of the model output to theinput parameters. These methods are computationally expensive, whereas a LSA method maysuffice. In Randall et al. [35], we performed a LSA on the neural model presented here, reproducedin Figure 4. LSA is useful in its relative ease in computation, especially using approximationsof derivatives with forward or central differences [17] and in its ability to calculate time-varyingsensitivities. In steady-state with nominal parameters close to their optimal value, LSA is veryuseful in ranking parameter influence as shown in our previous work [7, 24, 27, 30]. The resultsfrom the LSA are similar to the ranking generated with scalar SIs, especially in regard to the leastinfluential parameters τ H and τ b . Since every parameter in the LSA was above η = 10 − , everyparameter is influential. However, this is just a snapshot of the model sensitivity at one instancein the parameter space. The benefit of GSA is that it explores the entire parameter space Ω p andincorporates parameter interactions. Though there are some differences, the parameter rankingsfor both the LSA and GSA are similar, which could be due to the a priori calculation of nominalparameter values. Others have also found agreement in the calculation of the local and globalparameter influences [22]. In this study, we used SIs for the GSA. Though there are many other methods that explore theparameter space, including Morris screening [29] and derivative-based GSA methods [42], we usethis approach due to its broad applications. We developed LMSIs including a moving windowof width ∆, since the VM induces fast, transient changes in the steady-state behavior. Movingwindows have been implemented in graphical sensitivity analysis, but this method has difficultycapturing the effects between parameters [14]. The LMSIs mitigate this issue by incorporatingboth parameter interactions and the history of the parameter variance. Another advantage of thisapproach is that it can be used for virtually any modeling effort analyzing parameter sensitivityover time.Surely, not every problem requires time-varying GSA, and choosing to use this approach dependsgreatly on the questions asked and the QoI. Scalar SIs are appropriate if the QoI is a scalar, insteady-state, at a particular time point in the interval, or periodic with a constant frequency.However, if the objective is to understand parameter influence during a transition of states, thistime-varying approach could provide rich information. One such study is that of Calvo et al. [6],which calculated SIs for the parameters of a cardiovascular model studying head-up tilt at rest and25
Uniform Gaussian HR ( bp m ) Centered Trailing T o t a l t i m e - v a r y i ng S obo l i nd i c e s Figure 8: Moving window of integration, ∆, for a generic parameter θ i . The Valsalva maneuver (VM) phases arerepresented as alternating gray (I and III) and light gray (II and IV) boxes. Phase II is divided into early and lateparts (vertical dotted black line). (a) Common moving windows include uniform (top) for both centered and trailingwindows, full Gaussian for centered windows, and half-Gaussian for trailing windows. (b) Comparison of uniform(black) versus (Gaussian) windows of width ∆ = 15 s. The Gaussian window shows a steeper decline in parameterinfluence in phase IV. (c) Comparison of centered (black) versus trailing (red) windows of width ∆ = 15 s overlaid theheart rate (HR, bpm) data (blue, right axis). The centered window shows an increase in parameter influence beforethe VM occurs. (d) Comparison of varying window lengths for increasing ∆ = [5 , , ,
20] s. Pointwise-in-timeSobol’ indices (blue) and generalized Sobol’ indices (red) are also plotted. during the tilt. The transition between these states as the tilt occurs is of great clinical interest,and hence, this study can benefit from using the LMSIs proposed here to measure the changesin parameter influence. Another is the study by Sumner et al. [44] that claims to analyze time-dependent parameter influence with SIs. However, the QoI for this model is the state predictingglycogen synthase kinase evaluated when t = 60, which is a scalar value, and only time-dependentin the sense that the QoI is a time point. If other time points are of interest, especially involvingpeaks of various compounds, either GSIs or LMSIs can be used to quantify the parameter influenceover the entire time span up to 60 s. Moving windows have been used in signal processing for decades and relatively recently in graphicalsensitivity analysis [14, 43]. A moving window, ∆, typically results in smoothing of the signal todetermine mean behavior and enhance interpretability. However, there are several considerationswhen choosing an appropriate ∆ for analysis of the biological problem: shape, type, and width.Shape refers to the functional form providing the weights for the window. Type refers to whetherthe window is calculated as trailing the current time point t j or centered about t j . Width refers tothe time interval over which the window applies. We strongly suggest allowing the features of thesystem studied to dictate choice of window. In our study, much effort went into determining themost appropriate window for our problem. 26hape and type of window are interrelated, as choosing one can have bearing on the other.There are many potential window shapes, but here we discuss the most common types: uniformand Gaussian. Figure 8a displays uniform, half-Gaussian, and full Gaussian window. A uniformwindow will apply equal weight to each time point in the interval and, therefore, is independent oftype of window. However, a Gaussian shape can be influenced by whether the window is trailing orcentered. Choosing a half-Gaussian shape may be appropriate for trailing windows, placing moreemphasis on the most recent time points, whereas a full Gaussian would be more applicable for acentered window. Figure 8b compares the traces from both a uniform (black) and a Gaussian (red)shape, though the differences are negligible. We allowed the physiology to dictate the appropriatewindow, so we chose a trailing window with the shape of a half-Gaussian, diminishing importance ofthe signal going further back in time. The half-Gaussian window defines a continuous function alongthe interval I t , as the weights of the points outside the interval are zero. Figure 8c compares theresults for trailing half-Gaussian and a centered full Gaussian overlaid the heart rate trace. Clearly,the centered window places the parameter influence before the VM occurs, showing a change inparameter influence before the change to the steady-state even begins. Since the trailing windowcoincides with the VM itself, especially the increase in parameter influence during late phase II, wechose to allow the window to trail t j , as future time points do not impact the heart rate. Due to the overall model complexity, understanding the biological implications of the results andparameter interactions can be difficult. Therefore, model reduction can simplify these interactionsand still retain its predictive power. Many model reduction techniques exist from engineering andcontrol theory [3], aiming to reduce large numbers of state variables with many nonlinearities byattempting to mitigate the same inherent problem we address here: to what extent do the inputparameters affect the output. Our method using GSA to inform an analytical model reductionuses this idea to make appropriate choices for the exclusion of certain model components, unlikemethods that solely approximate the input-output relationship without considering the other pre-dicted model quantities [39]. To our knowledge, there are no previous studies that inform a modelreduction based on time-varying GSA methods for physiological models. For biological systemswith scalar QoIs, there has been a methodology proposed for model reduction via GSA by Marinoet al. [26]; however, a statistical analysis of a group of reduced models was not conducted.In this study, we have developed the LMSIs to specifically address this issue of substantial,transient disturbances in the steady-state behavior of the model and how the parameters impactthe variance of the model output. An advantage of this approach is that it takes into account higherorder parameter interactions and also the time correlation structure of the problem. Moreover, theimplementation of the moving window increased the interpretability of the PTSIs, which are toooscillatory to decipher. Hence, the LMSIs can provide a clearer ranking while still retaining thetransient fluctuations during the VM.Though we acknowledge that SIs may be impractical for very large differential equation systemswith hundreds of state variables and parameters, such as pharmacokinetics models, we suggest amulti-level approach. One can use simpler to compute but possibly less informative GSA measuresfirst to identify a set of unimportant parameters and then perform a more comprehensive analysison the remaining parameters [13]. We propose our GSA-informed model reduction methodologyas an alternative approach to the balanced truncation method in the model reduction formulationproposed by Snowden et al. [39]. 27 .5. Model selection
To our knowledge, no previous studies have performed a model selection protocol for cardiovascularmodels in response to the VM. We employ relative AICc and BIC scores to compare the modelfits to data. To consider the effect of reducing the full model on the other predicted quantitiesfor which data is difficult to acquire, we combine quantitative and qualitative approaches to selectwhether the delineation between the baroreceptors of the carotid sinus and aortic arch is necessaryto model, and, if not, which pathway should be modeled. We perform our analysis on our set ofmodels M = { m , m , m , m } where m is the original model with parameters s p,b , s s , and τ s held fixed at their nominal values, m is a reduced model of four states with τ b removed, m is thesame as m with convex combination parameter B = 0 emphasizing the aortic contribution only,and m is the same as m with B = 1 considering the carotid contribution only. The statisticalanalysis shows that m fit the data best with the greatest relative AICc and BIC scores comparedto the other models. This implies that it is necessary to model both the aortic and the carotidbaroreceptor regions to predict heart rate. We compared each of the reduced models to m via themetric Q given in equation (32), which measures the distance of the peak of T s,m . m has thelowest Q value is the closest to m . Ultimately, we conclude that m is the best model both tofit the heart data and predict the autonomic tones, and therefore, both regions are necessary tomodel the VM. Previous studies have modeled baroreceptor stimulation of only the carotid region[21, 23], though to our knowledge there are no studies solely modeling the aortic baroreceptors.Our previous work and the works of others [20, 35] support that modeling both regions is importantand the analysis of the present study also agrees with this assertion. The GSA method of choice is highly dependent on the model formulation and the QoI. Choosing acomputationally expensive GSA may not be feasible for models with long evaluation times, and theGSA results will differ based on the choice of QoI. For our choice of GSA, SIs assume that the modelparameters are statistically independent and may not be reasonable for some models. Since modelreduction was the focus of this study, we chose to take the analytical reduction approach, thoughwe could have fixed τ b at a constant value. Furthermore, analytical model reduction may also beimpractical for very large systems with many parameters. In this case, setting the parameters totheir nominal values may be more reasonable. Lastly, the results of the statistical analysis usingAICc and BIC scores is highly dependent on the available data. If a different heart data set hadbeen used, there is a possibility that the outcome of the model selection protocol could have beendifferent. However, we conducted this analysis on three representative control subjects and achievedsimilar results as those produced in Section 4 (not shown).
6. Conclusions
In this study, we introduced a model reduction and selection process informed by global sensitivityanalysis on a neurological model predicting heart rate and autonomic nervous function in responseto the Valsalva maneuver. We have established a systematic framework for parameter dimensionreduction and model simplification using global sensitivity analysis and developed the novel limited-memory Sobol’ indices that are aware of time history as well as transient dynamics. Furthermore, weapplied this methodology to a complex physiological model of the Valsalva maneuver response withextensive numerical experiments, simplifying the model while retaining important characteristics.From this methodology, we conclude that modeling both the aortic and carotid baroreceptor regionsare necessary to achieve the appropriate dynamics of the Valsalva maneuver.28 . Acknowledgments
We would like to thank Dr. Pierre Gremaud for his help in the development and critical analysisof the limited-memory Sobol’ indices.
8. Funding
This study was supported in part by the National Science Foundation under awards NSF/DMS1246991 and NSF/DMS 1557761 and in part by the National Institute of Health HL139813.
Appendix A.
This appendix contains calculations of the initial conditions (ICs) of the system in equation (14)and details concerning the parameter subset selection methods.
Appendix A.1. Initial conditions
The ICs are computed to ensure the model is in steady-state by taking into account the baselinesystolic blood pressure ¯ P , thoracic pressure ¯ P th , and heart rate ¯ H values, that is,¯ P c = ¯ P and ¯ P a = ¯ P − ¯ P th . (A.1)From here on, ¯ · denotes the steady-state condition. The ICs for the carotid ( ε b,c ) and aortic ( ε b,a )baroreceptor strains (equation (5)) are computed as¯ ε j = K b ¯ ε w,j = K b (cid:32) − (cid:115) e − q w ( ¯ P j − s w ) A + e − q w ( ¯ P j − s w ) (cid:33) , (A.2)where j = c or a for carotid and aortic, respectively. The steady-state neural integration equationis the algebraic relation ¯ n = B (cid:0) ¯ ε w,c − ¯ ε b,c (cid:1) + (1 − B ) (cid:0) ¯ ε w,a − ¯ ε b,a (cid:1) . (A.3)The efferent states have ICs ¯ T p,b = 0 . T s = 0 . T p,r = K p,r ¯ G p,r = K p,r e − q p,r ( ¯ P th − s p,r ) . (A.4)The IC for the model output heart rate is ¯ H . Appendix A.2. Subset selection and parameter estimation
Motivated by our previous work [30, 33, 45], we determine a subset of parameters ˆ θ to optimizeusing the structured correlation analysis approach. We compute the i th column of a sensitivitymatrix ( S ) of the model residual r with respect to the logarithm of the parameters (which vary inmagnitude) at time t j as S ij = ∂ r ( t j ) ∂ log θ i = ∂∂θ i H m ( t j ; θ ) − H d ( t j ) H d ( t j ) θ i = ∂H m ( t j ; θ ) ∂θ i θ i H d ( t j ) (A.5)29or H m the model output heart rate and H d the observed data. This analysis is local in thatthe derivatives computed are evaluated at the nominal parameter values. To determine possiblepairwise correlations between influential parameters as shown in our previous work [30], we calculatea correlation matrix c from a covariance matrix C = ( ˜S T ˜S ) − (A.6)as c ij = C ij (cid:112) C ii C jj , (A.7)where ˜ S are the columns of S corresponding to the parameters in the subset ˆ θ . The matrix c issymmetric with | c ii | = 1 and | c ij | ≤
1. For this study, we assigned a threshold of 0.9 for correlatedparameters.In the structured correlation analysis, we removed parameters from consideration a priori , suchas the sigmoidal steepness parameters ( q i ) and the half-saturation values ( s i ). Optimizing theseparameters can force the model to produce linear relations where nonlinearity occurs, producingresults that are not physiological [35]. We also excluded the time parameters t s and t e fromconsideration as they were determined directly from the data and are naturally very influential tothe model output.Applying these methods to the given model produced the identifiable parameter subsetˆ θ = [ B, τ p,b , τ p,r , H p,b , H p,r , H s ] T . (A.8)More details can be found in [35]. We fit H to data by minimizing the least squares cost ( J )between the model output and the data, i.e., J (ˆ θ ) = r T r + (cid:32) max j H m ( t j ; θ ) − max j H d ( t j )max j H d ( t j ) (cid:33) , (A.9)where r is the time-dependent residual vector given in equation 1. The last term in J is included toensure accurate estimates during the Valsalva maneuver. Figure 1f displays the model fit to heartrate data. Appendix B.
This appendix describes the implementation of approximating the integrands for equations (23)-(26). We employed estimators proposed by Saltelli et al. [37]. Let i = 1 , . . . , p for p the numberof parameters, t j be a time point where j = 1 , . . . , N for N the total number of discretized timepoints, and l = 1 , . . . , L be the index of L model function evaluations. Also, t = 0 and t N = T for T the end time point. We approximate the numerator for equation (21) and the integrands for thenumerators of equations (23) and (25) as V θ i ( E θ ∼ i [ f ( t j ; · ) | θ i ]) ≈ L L (cid:88) l =1 f ( t j , B ) l ( f ( t j , A ( i ) B ) l − f ( t j , A ) l ) (A.1)and the numerator for equation (22) and the integrands for the numerators of equations (24) and(26) as E θ ∼ i ( V θ i ( f ( t j ; · ) | θ ∼ i )) ≈ L L (cid:88) l =1 f ( t j , A ) l ( f ( t j , A ) l − f ( t j , A ( i ) B ) l ) , (A.2)30here A and B are two independent sampling matrices determined quasi-randomly using Sobol’sets. The matrix A ( i ) B contains all of the columns of A except the i th column, which is swappedwith the i th column of B . The denominators for equations (21) and (22) and the integrands ofthe denominators of (23)-(26) were approximated using the variance function var in MATLAB R (cid:13) ReferencesReferences [1] A. Alexanderian, P. A. Gremaud, R. C. Smith. Variance-based sensitivity analysis for time-dependent processes.
Reliab Eng Syst Safe , 196: 1–17. 2020. DOI: 10.1016/j.ress.2019.106722.[2] A. Alexanderian, J. Winokur, I. Sraj, A. Srinivasa, M. Iskandarani, W. C. Thacker, O. M.Knio. Global sensitivity analysis in an ocean general circulation model: a sparse spectralprojection approach.
Comput Geosci , 16 (3): 757–778. 2012. DOI: 10.1007/s10596-012-9286-2.[3] B. Besselink, U. Tabak, A. Lutowska, N. van de Wouw, H. Nijmeijer, D. J. Rixen, M. E.Hochstenbach, W. H. A. Schilders. A comparison of model reduction techniques from structuraldynamics, numerical mathematics, and systems and control.
J Sound Vib , 332 (19): 4403–4422.2013. DOI: 10.1016/j.jsv.2013.03.025.[4] W. F. Boron, E. L. Boulpaep.
Medical Physiology: A Cellular and Molecular Approach . 3rdEdition, Elsevier Inc., Philadelphia, PA, 2017. ISBN: 978-1455743773.[5] K. P. Burnham, D. R. Anderson.
Model selection and multimodel inference: a practi-cal information-theoretic approach . 2nd Edition, Springer-Verlag, New York, 2002. DOI:10.1007/b97636.[6] M. Calvo, V. Le Rolle, D. Romero, N. Behar, P. Gomis, P. Mabo, A. I. Hernandez. Globalsensitivity analysis of a cardiovascular model for the study of the autonomic response tohead-up tilt testing. in:
Conf Proc IEEE Eng Med BiolSoc , 5458–5461. 2018. DOI:10.1109/EMBC.2018.8513536.[7] L. M. Ellwein, H. T. Tran, C. Zapata, V. Novak, M. S. Olufsen. Sensitivity analysis and modelassessment: mathematical models for arterial blood flow and blood pressure.
Cardiovasc Eng
SIAM J Numer Anal
17 (2): 238–246. 1980. DOI: 10.1137/0717021.[9] N. Guglielmi, E. Hairer. Implementing Radau IIA methods for stiff delay differential equations.
Computing
67 (1): 1–12. 2001. DOI: 10.1007/s006070170013.[10] J. E. Hall.
Guyton and Hall Textbook of Medical Physiology . 13th Edition, Elsevier, Inc,Philadelphia, PA, 2016. ISBN: 978-1-4557-7005-2.[11] W. F. Hamilton, R. A. Woodbury, H. T. Harper, Jr. Arterial, cerebrospinal and venous pres-sures in man during cough and strain.
Am J Physiol
141 (1): 42–50. 1944. DOI: 10.1152/aj-plegacy.1944.141.1.42. 3112] J. L. Hart.
Extensions of global sensitivity analysis: theory, computation, and applications .Ph.D. thesis, North Carolina State University. 2018.[13] J. L. Hart, P. A. Gremaud, T. David. Global sensitivity analysis of high-dimensional neuro-science models: an example of neurovascular coupling.
Bull Math Biol
81: 1805–1828. 2019.DOI: 10.1007/s11538-019-00578-0.[14] B. Iooss, P. Lemaitre. A review on global sensitivity analysis methods. in:
UncertaintyManagement in Simulation-Optimization of Complex Systems: Algorithms and Applications .Research/Computer Science Interfaces Series, Springer, Boston, MA, 59: 101–122. 2015. DOI:10.1007/978-1-4899-7547-8 5[15] A. D. Jose, D. Collison. The normal range and determinants of the intrinsic heart rate in man.
Cardiovasc Res
J InverseIll-Posed Probl
25 (4): 499–519. 2017. DOI: 10.1515/jiip-2016-0024.[17] C. T. Kelley.
Iterative Methods for Optimization . Frontiers in Applied Mathematics, Society forIndustrial and Applied Mathematics, Philadelphia, PA. 1996. DOI: 10.1137/1.9781611970920.[18] A. Kiparissides, S. S. Kucherenko, A. Mantalaris, E. N. Pistikopoulos. Global sensitivityanalysis challenges in biological systems modeling.
Ind Eng Chem Res
48 (15): 7168–7180.2009. DOI: 10.1021/ie900139x.[19] P. I. Korner, A. M. Tonkin, J. B. Uther. Reflex and mechanical circulatory effects ofgraded Valsalva maneuvers in normal man.
J Appl Physiol
40 (3): 434–440. 1976. DOI:10.1152/jappl.1976.40.3.434.[20] S. A. Kosinski, B. E. Carlson, S. L. Hummel, R. D. Brook, D. A. Beard. Computationalmodel-based assessment of baroreflex function from response to Valsalva maneuver.
J ApplPhysiol
125 (6): 1944–1967. 2018. DOI: 10.1152/japplphysiol.00095.2018.[21] V. Le Rolle, A. I. Hernandez, P. Y. Richard, G. Carrault. An autonomic nervous systemmodel applied to the analysis of orthostatic tests.
Model Simul Eng
PLoS One
13 (7): 1–38. 2018. DOI: 10.1371/journal.pone.0200917.[23] K. Lu, J. W. Clark Jr., F. H. Ghorbel, D. L. Ware, A. Bidani. A human cardiopulmonarysystem model applied to the analysis of the Valsalva maneuver.
Am J Physiol Heart CircPhysiol
281 (6): H2661–H2679. 2001. DOI: 10.1152/ajpheart.2001.281.6.H2661.[24] G. Mader, M. S. Olufsen, A. Mahdi. Modeling cerebral blood flow velocity during orthostaticstress.
Ann Biomed Eng
43 (8): 1748–1758. 2015. DOI: 10.1007/s10439-014-1220-4.[25] A. Mahdi, J. Sturdy, J. T. Ottesen, M. S. Olufsen. Modeling the afferent dynamics ofthe baroreflex control system.
PLoS Comput Biol
J Theor Biol
Math Biosci
Techno-metrics
33 (2): 161–174. 1991. DOI: 10.2307/1269043.[29] C. H. Olsen, J. T. Ottesen, R. C. Smith, M. S. Olufsen. Parameter subset selection techniquesfor problems in mathematical and biology.
Biol Cybern
113 (1-2): 121–138. 2019. DOI:10.1007/s00422-018-0784-8.[30] M. S. Olufsen, J. T. Ottesen. A practical approach to parameter estimation applied to modelpredicting heart rate regulation.
J Math Biol
67 (1): 39–68. 2013. DOI: 10.1007/s00285-012-0535-8.[31] M. S. Olufsen, J. T. Ottesen, H. T. Tran, L. M. Ellwein, L. A. Lipsitz, V. Novak. Blood pressureand blood flow variation during postural change from sitting to standing: model developmentand validation.
J Appl Physiol
99 (4): 1523–1537. 2005. DOI: 10.1152/japplphysiol.00177.2005.[32] M. S. Olufsen, H. T. Tran, J. T. Ottesen, L. A. Lipsitz, V. Novak. Modeling baroreflexregulation of heart rate during orthostatic stress.
Am J Physiol Regul Integr Comp Physiol
291 (5): R1355–R1368. 2006. DOI: 10.1152/ajpregu.00205.2006.[33] S. R. Pope, L. M. Ellwein, C. Zapata, V. Novak, C. T. Kelley, M. S. Olufsen. Estimationand identification of parameters in a lumped cerebrovascular model.
Math Biosci Eng
Biomech Model Mechanobiol
18 (1):219–243. 2019. DOI: 10.1007/s10237-018-1078-8.[35] E. B. Randall, A. Billeschou, L. S. Brinth, J. Mehlsen, M. S. Olufsen. A model-based analysisof autonomic nervous function in response to the Valsalva maneuver.
J Appl Physiol
127 (5):1386–1402. 2019. DOI: 10.1152/japplphysiol.00015.2019.[36] E. B. Randall, N. Z. Randolph, M. S. Olufsen. Persistent instability in a nonhomogeneousdelay differential equation system of the Valsalva maneuver.
Math Biosci
ComputPhys Commun
181 (2): 259–270. 2010. DOI: 10.1016/j.cpc.2009.09.018.[38] R. C. Smith.
Uncertainty quantification: theory, implementation, and applications . Society forIndustrial and Applied Mathematics (SIAM), Philadelphia, PA. 2014. ISBN: 978-1611973211.3339] T. J. Snowden, P. H. van der Graaf, M. J. Tindall. Model reduction in mathematical phar-macology.
J Pharmacokinet Pharmacodyn
45 (4): 537–555. 2018. DOI: 10.1007/s10928-018-9584-y.[40] I. M. Sobol. Sensitivity estimates for nonlinear mathematical models.
Math Mod Comp Exp
Math Comput Simul
55 (1-3): 271–280. 2001. DOI: 10.1016/S0378-4754(00)00270-6.[42] I. M. Sobol, S. S. Kucherenko. Derivative based global sensitivity measures.
Procedia SocBehav Sci
Reliab Eng Syst Safe
93 (1): 28–54. 2008. DOI: j.ress.2006.10.012.[44] T. Sumner, E. Shephard, I. D. L. Bogle. A methodology for global sensitivity analysis oftime-dependent outputs in systems biology modelling.
J R Soc Interface
Adv Appl Math Mech
Neurologist
16 (4): 215–222. 2010. DOI: 10.1097/NRL.0b013e3181cf86ab.[47] E. Wit, E. ven den Heuvel, J.-W. Romeijn. All models are wrong : an introduction to modeluncertainty.