Roskilde University Cardiovascular dynamics during head-up tilt assessed via pulsatile and non-pulsatile models Williams, Nakeya; Brady, Renee; Gilmore, Steven; Gremaud, Pierre; Tran, Hien; Ottesen, Johnny T.; Mehlsen, Jesper; Olufsen, Mette Published in: Journal of Mathematical Biology DOI: 10.1007/s00285-019-01386-9 Publication date: 2019 Document Version Peer reviewed version Citation for published version (APA): Williams, N., Brady, R., Gilmore, S., Gremaud, P., Tran, H., Ottesen, J. T., Mehlsen, J., & Olufsen, M. (2019). Cardiovascular dynamics during head-up tilt assessed via pulsatile and non-pulsatile models. Journal of Mathematical Biology, 79(3), 987–1014. https://doi.org/10.1007/s00285-019-01386-9 General rights Copyright and moral rights for the publications made accessible in the public portal are retained by the authors and/or other copyright owners and it is a condition of accessing publications that users recognise and abide by the legal requirements associated with these rights. • Users may download and print one copy of any publication from the public portal for the purpose of private study or research. • You may not further distribute the material or use it for any profit-making activity or commercial gain. • You may freely distribute the URL identifying the publication in the public portal. Take down policy If you believe that this document breaches copyright please contact rucforsk@ruc.dk providing details, and we will remove access to the work immediately and investigate your claim. Download date: 28. Nov. 2021 J Math Biol manuscript No. (will be inserted by the editor) Cardiovascular dynamics during head-up tilt assessed via pulsatile and non-pulsatile models Nakeya D. Williams · Renee Brady · Steven Gilmore · Pierre Gremaud · Hien T Tran · Johnny T Ottesen · Jesper Mehlsen · Mette S. Olufsen Received: date / Accepted: date Abstract This study develops non-pulsatile and pulsatile models for prediction of blood flow and pressure during head-up tilt (HUT). This test is used to diagnose potential pathologies within the autonomic control system, which acts to keep the cardiovascular system at homeostasis. We show that mathematical modeling can be used to predict changes in cardiac contractility, vascular resistance, and arterial compliance, quantities that cannot be measured, but are useful to assess the system’s state. These quantities are Nakeya D. Williams United States Military Academy, West Point, NY E-mail: nakeya.williams@usma.edu Rene Brady H Lee Moffitt Cancer Center and Research Institute, Tampa, FL E-mail: renee.brady@moffitt.org Steven Gilmore NC State University, Raleigh, NC E-mail: sjgilmo2@ncsu.edu Pierre Gremaud NC State University, Raleigh, NC E-mail: gremaud@ncsu.edu Hien Tran NC State University, Raleigh, NC E-mail: tran@ncsu.edu Johnny T Ottesen Roskilde University, Roskilde, Denmark E-mail: johnny@ruc.dk Jesper Mehlsen Frederiksberg and Bispebjerg Hospital, Frederiksberg, Denmark E-mail: Jesper.Mehlsen@regionh.dk Mette Olufsen NC State University, Raleigh, NC E-mail: msolufse@ncsu.edu 2 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen predicted as time-varying parameters modeled using piecewise linear splines. Having models with various levels of complexity formulated with a common set of parameters, allows us to combine long-term non-pulsatile simulations with pulsatile simulations on a shorter time-scale. We illustrate results for a representative subject tilted head-up from supine position to a 60 degree angle. The tilt is maintained for 5 minutes before the subject is tilted back down. Results show that if volume data is available for all vascular compartments three parameters can be identified, cardiovascular resistance, vascular compliance, and ventricular contractility, whereas if model predictions are made against arterial pressure and cardiac output data alone, only two parameters can be estimated either resistance and contractility or resistance and compliance. Keywords Cardiovascular dynamics modeling · head-up tilt · pulsatile vs. non-pulsatile modeling · parameter estimation · orthostatic intolerance Mathematics Subject Classification (2000) 92C30 · 92C35 · 92C42 · 92C50 · 76Z05 1 Introduction Emergency rooms and syncope clinics see a large number of patients who experience chronic episodes of dizziness and syncope. These symptoms may be associated with orthostatic intolerance: the inability to maintain blood pressure and flow in response to active standing or head-up tilt (HUT). Orthostatic intolerance is triggered by a number of factors, the most important being dysautonomia, a disorder associated with the autonomic nervous system [12]. This phenomenon is difficult to diagnose not only due to the ambiguity of the disorder, but also due to the inability to obtain clinical data for quantities pertinent to cardiovascular regulation. The HUT protocol typically includes continuous measurements of heart rate and blood pressure, and sometimes cardiac output. The protocol starts by measuring these quantities while the subject is resting in supine position. After steady oscillating values for heart rate and blood pressure are recorded, the subject is tilted head up to a 60 degree angle. After a period of 5-30 minutes, the subject is tilted back to supine position. Upon HUT, blood is pooled in the lower extremities causing a drop in blood pressure in the upper body, while blood pressure in the lower body is increased. When the subject is returned to supine position, the blood pressure returns to supine levels. In response to HUT, the autonomic nervous system elicits an increase in heart rate, cardiac contractility, and peripheral resistance, redistributing blood volume and thereby re-establishing homeostasis. It is known that the cardiovascular control system modulates the aforementioned parameters in response to an orthostatic stimulus. This modulation can be attenuated or accentuated in disease, yet these quantities cannot be measured. Therefore, diagnosis is made from system level measurements such as heart rate, blood pressure, and cardiac output. Several groups have successfully developed mathematical models of the cardiovascular system and its control (e.g. [8, 9, 14, 17, 22, 21, 27, 25]). Only some Pulsatile and non-pulsatile cardiovascular model 3 of these compare predictions to data and only a few have attempted to predict HUT dynamics [13, 14, 22]. These models vary in complexity; there is a clear trade-off between simplicity and the number of model parameters. For example, models by Olufsen [19] and Lim [13] include 100 parameters, while the model by Pope et al. [21] uses 15 parameters to predict cerebral blood flow and arterial blood pressure. Other contributions include the model by TenVoorder [27] predicting short-term blood pressure and heart rate variability for a healthy young male, and models by Ursino predicting heart rate regulation [28, 29, 25]. The latter compared the model output with experimental data but did not address parameter estimation. Moreover, to our knowledge no previous studies have developed a set of models (pulsatile and non-pulsatile) that can be simulated using a single parameter set. The same applies to non-pulsatile models. Several non-pulsatile models exist [5, 30, 1], yet most incorporate an additional equation (derived from Frank-Starling’s law [26, 6]) that relates stroke volume to arterial and venous pressure and volume. This study develops a lumped parameter model that predicts controlled quantities from patient-specific measurements of heart rate and blood pressure. To do so, we develop pulsatile and non-pulsatile models predicting dynamics over short (1-3 minutes) and long (> 5 minutes) time scales. The novelty is the development of a non-pulsatile model obtained by averaging the pulsatile model over each cardiac cycle. This methodology allows us to identify a parameter set that can represent the same physiological quantities across two models. The non-pulsatile and pulsatile models are compared, and results show that the two models elicit similar dynamics. The non-pulsatile model has multiple advantages. It is less complex making it easier to run simulations over long time intervals (15-45 minutes), e.g. to study the system collapse during fainting. Secondly, the non-pulsatile model can easily be coupled with models representing other parts of the physiological response, e.g. the respiratory system, with a time-scale of hours, the renal system, which acts on a time-scale of hours-days, or the hormonal system, which has a time-scale of hours [3, 6]. Finally, having both models allows us to speed up computations in studies for which it is not necessary to account for discrete events associated with opening and closing of the heart valves, and to identify short segments over which we can predict the more complex pulsatile dynamics. 2 Methods This study develops simple patient-specific pulsatile and non-pulsatile models that use heart rate as an input to estimate arterial blood pressure during HUT. We first describe the data used for model predictions, followed by a description of the model development. Next, we discuss gravitational effects and autonomic regulation associated with the HUT procedure. 4 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen 2.1 Data Pulsatile measurements include ECG signals recorded using standard precordial leads and blood pressure recorded using photoplethysmography (Finapres Medical Systems B.V.). Data was collected at the Coordinating Research Centre at Frederiksberg and Bispebjerg Hospital, Frederiksberg, Denmark from a healthy young male volunteer with no known heart or vascular diseases. The subject was left to rest in supine position for 10 minutes. Subsequently, he was tilted to an angle of 60 degrees at a speed of 15 degrees per second measured by way of an electronic marker. The subject remained tilted for five minutes, and was then tilted back to supine position at the same tilt speed. 2.2 Lumped parameter cardiovascular models Similar to our previous study [31], we use a lumped parameter model estimating volume (V ), flow (q), and pressure (p), in five compartments within the systemic circulation, see Figure 2 and Table 1. The model includes four compartments representing arteries and veins in the upper and lower body, as well as a compartment representing the heart. As in our previous studies, we omit the pulmonary circuit from our model since we do not have any pulmonary circulation data and our work focuses on aspects of the systemic circulation. However, with available data there is no reason why that the pulmonary circuit could not be included in future work. For each compartment, a pressure-volume relation can be defined as Vi,s = Vi − Vi,us = Ci (pi − pi,us ), (1) where Vi,s denotes the stressed volume in compartment i, Vi (ml) is the total volume, Vi,us (ml) the unstressed volume, Ci (ml/mmHg) the compartment compliance, pi (mmHg) the instantaneous blood pressure, and pi,us (mmHg) 110 110 110 100 100 100 90 90 90 80 80 80 70 70 70 60 60 60 50 50 p (mmHg) & H (bpm) Supine 40 0 100 Supine HUT 200 300 400 time (sec) (a) 500 600 40 50 120 140 160 180 40 150 200 time (sec) time (sec) (b) (c) 250 Fig. 1: (a) Pulsatile carotid blood pressure (light gray) and heart rate (dark gray) over the entire experiment (head-up tilt (HUT) followed by head-down tilt). (b) Dynamics during supine position (before HUT), and (c) Blood pressure and heart rate dynamics during the initial HUT phase. Pulsatile and non-pulsatile cardiovascular model 5 (assumed constant) is the associated unstressed pressure. For each compartment, the change in volume is given by dVi = qin − qout , dt (2) where qin and qout denote the flow in and out of the compartment. Using a linear relationship analogous to Ohm’s law, the volumetric flow q (ml/s) between compartments can be computed as q= (a) pin − pout , R (3) (b) Fig. 2: (a) Compartment model. For each compartment we denote blood pressure by p (mmHg), volume by V (ml), and compliance by C (ml/mmHg). The compartments represent the upper body arteries (subscript au), lower body arteries (subscript al), upper body veins (subscript vu), lower body veins (subscript vl), and the left heart (subscript lh). Resistances R (mmHg s/ml) are placed between all compartments: Ral denotes the resistance between arteries in the upper and lower body, Raup and Ralp denote resistance between arteries and veins in the upper and lower body, respectively, while Rmv and Rav denote the resistance to flow entering and leaving the heart. For the pulsatile model, the two heart valves, the mitral valve and the aortic valve, are modeled as diodes. (b) Differences between the pulsatile (top panel) and non-pulsatile (center panel) models in the left heart compartment, and a sketch of a generic model component (bottom panel). Table 1: Abbreviations (subscripts) used in the compartmental model. Abbreviation au al aup alp vu vl lh av mv Name upper body arteries lower body arteries upper body ”peripheral” arteries (vascular bed) lower body ”peripheral” arteries (vascular bed) upper body veins lower body veins the left heart (ventricle and atrium) aortic valve mitral valve 6 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen where pin and pout are the pressures on either side of the resistor R (mmHg s/ml). Similar to [31], the left ventricular pressure is given by plh = Elh (t)Vlh,s , Vlh,s = Vlh − Vlh,us , (4) where Elh (mmHg/ml) is the left heart elastance (the reciprocal of its compliance) and Vlh,s is the stressed left heart volume, Vlh the total volume and Vlh,us is the unstressed volume. Pumping is achieved by introducing a variable elastance function [4] of the form     EM −Em π t̃  , t̃ ≤ TS 1 − cos  2 TS + Em      Elh (t̃) = EM −Em cos π(t̃−TS ) + 1 + Em , TS ≤ t̃ ≤ TR (5) 2 TR −TS     Em , TR ≤ t̃ ≤ T, where t̃ is the time within a cardiac cycle T = 1/H, for the heart rate H. Parameters Em and EM denote the minimum and maximum elastance, respectively. TS = αS T and TR = αR T denote the time at which maximum elastance and minimum elastance occur, respectively. For each cardiac cycle, elastance is increased for 0 < t̃ < TS and decreased for TS < t̃ < TR , while for TR < t̃ < T elastance is kept constant at its minimum value. Values for T and αS are obtained from data, while αR is a model parameter. The time-varying elastance function is illustrated in Figure 3(d). Finally, heart valves are modeled as the electrical current in a diode, i.e. the flow through the valve is given by (3) when the valve is open and is zero otherwise. Figure 3(c) shows flow in and out of the valves over a cardiac cycle. 2.2.1 Non-pulsatile Model The non-pulsatile model is derived by integrating the pulsatile model during filling and ejection. The derivation presented here is specific to the elastance function (5), but can be generalized to any function that can be integrated during filling and ejection. Integration during the two phases allows us to eliminate the differential equation predicting volume in the heart, replacing it with an empirical equation relating venous and arterial pressures and volume. The derivation is carried out in two steps first determining relations characterizing filling followed by integration during ejection. They are all build on the assumption that the change in left ventricular volume (2), is given by dVlh = qmv − qav , dt where during filling (pvu > plh ) the mitral valve is open, i.e., dVlh pvu − plh , = qmv = dt Rmv (6) Pulsatile and non-pulsatile cardiovascular model 7 whereas during ejection (plh > pau ), the aortic valve is open, i.e., dVlh plh − pau . = −qav = − dt Rav (7) To obtain the non-pulsatile model, we integrate over one period computing the average flow through the heart Q, which during filling is given by Q = q mv = 1 T Z T TR dVlh dt. dt The latter holds since the mitral valve is closed for 0 < t < TR (as shown in Figure 3). Integration gives Q= i 1h VED − VES , T (8) where VED and VES denote the end-diastolic and end-systolic ventricular volumes, respectively. Similarly, during ejection the average flow through the heart, which by conservation of flow, also amounts to Q, can be computed from Q = q av 1 =− T Z TS 0 dVlh dt. dt Note that during ejection, the volume in the left ventricle is decreasing, indicated by a minus in front of the integral. Again, this holds since the aortic valve is closed for TS < t < T, i.e. Q= i 1h VED − VES . T Keeping in mind that heart rate H = 1/T , the stroke volume Vstr can be computed from Vstr = VED − VES . (9) Finally, using (8) cardiac output Q can be computed as Q = HVstr . (10) To approximate values for VED and VES , we integrate (6) and (7) during filling and ejection, respectively. 8 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen Filling: Inserting (4) in (6) gives V̇lh + Elh pvu + Elh Vlh,us Vlh = . Rmv Rmv Assuming pvu = pvu is constant and Elh (t) = Em , integration gives Vlh (t) ≈ pvu + Vlh,us + Ce−Em t/Rmv . Em (11) Assuming Vlh (0) = VES (start of filling), evaluation of (11) at end-diastolic volume V = VED gives   p (12) VED ≈ VES e−Em TD /Rmv + vu + Vlh,us 1 − e−Em TD /Rmv , Em where TD = T − TR is the time to fill, as shown in Figure 3. Ejection: Inserting (4) in (7) gives V̇lh + Elh Elh Vlh,us + pau Vlh = . Rav Rav Assuming pau = pau is constant and Elh (t) = EM , integrating gives Vlh (t) ≈ pau + Vlh,us + Ce−EM t/Rav . EM (13) Assuming Vlh (0) = VED (start of ejection), evaluation of (13) at end-systolic volume V = VES gives p   au VES ≈ e−EM TR /Rav VED + + Vlh,us 1 − e−EM TR /Rav , (14) EM where TR is the time at end-systole. Integrated model: Defining k1 = e−EM TR /Rav and k2 = e−Em TD /Rmv , allows us to rewrite (12) and (14) as  p  VED = k2 VES + vu + Vlh,us 1 − k2 Em   p au + Vlh,us 1 − k1 . VES = k1 VED + EM Solving (15) and (16) gives p  (1 − k )k  (1 − k ) p 1 2 2 au + Vlh,us + vu + Vlh,us . VED = EM (1 − k1 k2 ) Em (1 − k1 k2 ) p p  (1 − k )k  (1 − k ) 2 1 1 vu au VES = + Vlh,us + Vlh,us + Em (1 − k1 k2 ) EM (1 − k1 k2 ) (15) (16) (a) p (mmHg) Pulsatile and non-pulsatile cardiovascular model 100 Left Heart Arterial 50 0 0 0.1 0.2 0.3 0.4 0.5 SYSTOLE 0.6 0.7 0.8 0.9 1 DIASTOLE 150 V (ml) (b) VED 100 50 VES 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 1 q/q max (c) q av q 0.5 E (mmHg/ml) 0 (d) 9 0 0.1 0.2 0.3 0.4 mv 0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1 3 2 1 0 TS 0 0.1 TR 0.2 0.3 0.4 time (sec) Fig. 3: Systolic and diastolic phases of the cardiac cycle of length T (scaled to 1 sec). (a) Left ventricular (black) and aortic (blue dashed) pressures. (b) Left ventricular volume. The horizontal dashed lines represent end-diastolic (VED ) and end-systolic (VES ) volumes. (c) Flow through the mitral (qmv ) and aortic (qav ) valves, note these are zero when valves are closed. (d) Time-varying elastance during a cardiac cycle. The maximum elastance is found at t̃ = TS and the minimal elastance at t̃ = TR . for k1 k2 6= 1. Finally, using (17) stroke volume becomes Vstr = (1 − k1 )(1 − k2 )  pau p  − vu . k1 k2 − 1 EM Em Note, that for the exponentials we have k1 ≈ k2 ≈ 0 (since the resistances Rav and Rmv are small). Thus, assuming that (k1 , k2 ) → (0, 0) the stroke volume can be approximated by p p  Vstroke ≈ − au − vu . (17) EM Em Pulsatile and non-pulsatile models: Using these relations, the pulsatile and non-pulsatile differential equations can be summarized as Pulsatile dVau dt dVal dt dVvl dt dVvu dt dVlh dt = qav − qaup − qal = qal − qalp = qalp − qvl = qvl + qaup − qmv = qmv − qav Non-pulsatile dVau dt dVal dt dVvl dt dVvu dt = Q − qaup − qal (18) = qal − qalp (19) = qalp − qvl (20) = qvl + qaup − Q (21) 10 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen where pau − pvu Raup pau − pal = Ral pal − pvl = Ralp pvl − pvu = Rvl qaup = qal qalp qvl plh − pau Rav pvu − plh = Rmv qav = qmv Q = H Vstr = H  pau pvu − Em EM  where pi Ci = Vi − Vi,us by equation (1) with i referring to the compartments in question i = {au, vu, al, vl, aup, alp, lh} listed in Table 1 and illustrated in Figure 2. Note we use the same notation for the pulsatile and non-pulsatile models, yet for the non-pulsatile models flows q, volumes V and pressures p refer to quantities averaged over the cardiac cycle, fulfilling assumptions (that the aortic and venous pressure is constant during filling and ejection, respectively) for the derivation of the non-pulsatile model. The main differences between the two models are highlighted in blue (pulsatile) and red (non-pulsatile). For the pulsatile model, the left ventricular pressure (plh ) is predicted from (4), while the non-pulsatile model does not have a heart compartment. These equations were solved in Matlab using ODE15s. 2.2.2 Initial conditions Both pulsatile and non pulsatile models conserve volume, hence the systems of equations can be reduced. For the pulsatile model Vtot = Vau + Val + Vvl + Vvu + Vlh , whereas for the non-pulsatile model the total volume is Vtot = Vau + Val + Vvl + Vvu . For the non-pulsatile model, we let Val = Vtot − Vau − Vvl − Vvu allowing us to rewrite the system of equations as dx = Ax + b, dt (22) where x = {Vau , Vvl , Vvu }, A is a constant coefficient matrix, and b (a vector function of the total volume Vtot ). At nominal parameter values (assuming a constant heart rate), the eigenvalues of A are real and negative, i.e. the system reaches a steady state. We set initial conditions for the model to this equilibrium point. Note, in general, heart rate is an input that varies with time. 2.2.3 Nominal parameter values Literature values and subject specific information were integrated to identify nominal values for all model parameters (resistances, compliances, heart, and , Pulsatile and non-pulsatile cardiovascular model 11 HUT parameters). Similar to calculations by Hoppensteadt and Peskin [10], nominal parameter values were obtained by backward calculation of ( 1) and (3) using mean approximations for all pressures, flows, and volumes at rest (i.e. before HUT). Pressure. The mean pressure in the upper arteries, p̄au , was estimated from the averaged pressure data. The pressure drop along the large arteries is small. As a result, in supine position, the lower body arterial pressure is of the same magnitude as the upper body pressure [6, 3]. To enforce flow in the right direction, we let p̄al = 0.98p̄au . No data are available for the venous pressures, therefore these were approximated from literature [6, 3]. As suggested by [10], we set the upper body venous pressure p̄vu = 3.5 and the lower body venous pressure p̄vl = 3.75. Volume. The total blood volume was estimated from body surface area [24] as  (3.47 · BSA − 1.954) 1000, Female Vtot = (23) (3.29 · BSA − 1.229) 1000, Male p where BSA = h w/3600 denotes the body surface area, h (cm) the height and w (kg) the weight of the subject studied. We use common physiological approximations for the distribution of blood volume, assuming that 85% of the total blood volume is in the systemic circuit (the pulmonary circuit is not included in the model), with 20% in the arteries and 80% in the veins [6, 3]. We further assume that 85% of the blood volume is in the upper body and 15% is in the lower body. Finally, for each compartment, we distribute the stressed and unstressed blood volume as proposed by [2], given in Table 2. Table 2: For each compartment, the volume is estimated as fractions of the total volume Vtot , and the total compartment blood volume is separated between a stressed and an unstressed volume, i.e., Vtot,i = Vi,s + Vi,us for i ∈ {au, al, vu, vl}. Stressed volumes are calculated as a fraction of the total systemic volume. Values are computed using volume distributions reported by [2]. In this study, the upper body compartments contain arteries and veins in the head, thorax, and abdomen, while the lower body compartments contain arteries and veins in the legs. Volume Vtot,sys Vtot,a Vtot,v Vtot,au Vtot,al Vtot,vu Vtot,vl Vau,s Val,s Vvu,s Vvl,s Position Total systemic circuit Total systemic arteries Total systemis veins Total upper body arteries Total lower body arteries Total upper body veins Total lower body veins Stressed upper body arteries Stressed lower body arteries Stressed upper body veins Stressed lower body veins Value 0.85 Vtot 0.20 Vtot,sys 0.80 Vtot,sys 0.85 Vtot,a 0.15 Vtot,a 0.85 Vtot,v 0.15 Vtot,v 0.18 Vtot,au 0.18 Vtot,al 0.05 Vtot,vu 0.05 Vtot,vl 12 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen Cardiac output was estimated from the assumption that the entire volume is circulated in the body within one minute [4]. The average flow between the upper and lower body is determined under the assumption that in supine position, 80% of the blood flows between the upper body arteries and veins, while 20% of the blood flow supports the lower extremities [2]. Parameters. Using these assumptions combined with Ohm’s Law (3) and the pressure volume equation (1), we are able to estimate resistors and capacitors (in supine position) as p̄in − p̄out , q̄ Vs V̄ − V̄us = , C= p̄ p̄ R= where p̄, q̄, V̄ denote mean values for the blood pressures, flow, and volumes. Finally, for the heart model, parameters representing the minimum and maximum elastance, as well as timing of the pump function, must be estimated. The minimum left ventricular elastance can be calculated using the pressurevolume relation (1), where the left ventricular volume equals the end-diastolic volume (VED ), i.e. p̄vu = Em (VED − Vlh,us ). Similarly, the maximum left ventricular elastance can be predicted by assessing the same relation at the end-systolic phase. For this case, p̄au = EM (VES − Vlh,us ). For both parameters, we assume that the unstressed value of the ventricular volume Vlh,us = 10 ml, which was used in previous studies (e.g. [4]). The timing of the pump is achieved via parameters TS = αS T, TR = αR T and T , where αS (estimated from data) denotes the fraction of the cardiac cycle at which the elastance is maximal and αR (a parameter) determines the fraction of the cardiac cycle at which elastance is minimal. Table 3 specifies parameter values and units for all model parameters. 2.2.4 Modeling HUT As the subject is tilted head up (shown in Figure 4), blood is pooled in the lower extremities leading to an increase in pressure in the lower body, while pressure in the upper body is decreased. When the subject is tilted back to supine position, blood is returned to the upper body, decreasing pressure in the legs and increasing pressure in the upper body. Similar to [19, 18], the effect of gravity is included by adding hydrostatic pressure to qal , and subtracting Pulsatile and non-pulsatile cardiovascular model 13 Table 3: Model parameters: The nominal parameter values and equations. Note, LB and UB refer to lower body and upper body, respectively. Name Vtot l w CO p̄au p̄al p̄vu p̄vl q̄up q̄low Rav Rmv Definition total blood volume height weight cardiac output UB mean arterial bp LB mean arterial bp UP mean venous bp LB mean venous bp UB flow LB flow open aortic valve open mitral valve Raup UB peripheral resist Ralp LB peripheral resist Ral UB arterial resist Rvl LB venous resist Cau Cal Cvu Cvl VED VES Vlh,us UB arterial comp LB arterial comp UB venous comp LB venous comp end-diastolic vol end-systolic vol unst ventricular vol Em min elastance EM max elastance T cardiac cycle length ND (non dimensional), ∗ ratio Equation Eqn (23) – – Vtot /60 mean(bp data) 0.98 p̄au – – 0.8 CO 0.2 CO – – p̄au − p̄vu q̄up p̄al − p̄vl q̄low p̄au − p̄al q̄low p̄vl − p̄vu q̄low V̄au,s /p̄au V̄al,s /p̄al V̄vu,s /p̄vu V̄vl,s /p̄vl – – – p̄vu VED − Vlh,us p̄au VES − Vlh,us data of T Nom Value 5584 186 83 93.07 68.08 66.72 3.50 3.75 74.45 18.61 0.001 0.001 Opt Value – – – – – – – – – – – – Units ml cm kg ml/s mmHg mmHg mmHg mmHg ml/s ml/s mmHg s/ml mmHg s/ml 0.8673 0.9290 mmHg s/ml 3.3828 – mmHg s/ml 0.0731 – mmHg s/ml 0.0134 – mmHg s/ml 2.1334 0.3765 46.1081 7.5943 142 47 10 2.1366 – – – – – – ml/mmHg ml/mmHg ml/mmHg ml/mmHg ml ml ml 0.0265 0.0265 mmHg/ml 1.8399 – mmHg/ml 0.98 – s it from qvl in Ohm’s law (3), determining the flow q as the ratio of pressure p to resistance R as ρghtilt sin (φ(t)) + pin − pout , R  0 t < tus     vt (t − tus ) tus ≤ t < tus + tue  π 60 tus + tue ≤ t < tds φ(t) = 180   −v (t − t ) t  t ds ds ≤ t ≤ tds + tde   0 t > tds + tde , q= (24) where ρ (g/ml) is the blood density, g (cm/s2 ) is the gravitational acceleration constant, htilt (cm) is the absolute height between the upper body and lower body compartments, φ(t) is the tilt angle (in radians), vt = 15 degrees/s is the tilt speed, while tus , tue , tds , and tde denote the times at which HUT is started and ended for tilting up and back down, respectively. The combined 14 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen term ρghtilt sin (φ(t)) denotes the hydrostatic pressure between the upper and lower body compartments. 2.2.5 Modeling effects of cardiovascular regulation Upon HUT, the baroreceptor nerve firing is modulated by the aortic and carotid sinus baroreceptors sensing changes in the stretch of the arterial wall. Typically, HUT leads to a decrease in blood pressure mediating an increase in sympathetic outflow along with parasympathetic withdrawal. Sympathetic stimulation elicits changes in peripheral vascular resistance (Raup and Ralp ), vascular tone (represented by compliance C), and cardiac contractility (determined from the elastance parameters EM and Em ), while parasympathetic withdrawal primarily affects heart rate (H) and cardiac contractility (again EM and Em ), [3]. Tilting back to supine position elicits the opposite response. Heart rate is used as an input, consequently, parasympathetic regulation is implicitly accounted for in the model. Regulation of cardiac contractility was modeled by allowing the minimum elastance Em of the left heart to vary with time, while keeping the maximal elastance constant. All capacitors appear in series. As a result only one is identifiable. Given that data is available for the upper body arterial pressure, we chose to let Cau vary with time. Within the pulsatile model this allows us to accurately fit pulse pressure. Finally, vascular resistance in both the upper and lower body (Raup and Ralp , respectively) are also varied with time. The resistances associated with compartments representing the upper and lower body arteries appear in parallel, hence both resistances are not identifiable. Thus, Raup is controlled directly, while we let Ralp = kRaup , where k is the ratio of the optimized supine values of Raup and Ralp . Fig. 4: The HUT test: The subject depicted is tilted to an angle of 60 degrees at a constant speed of 15 degrees per second. Pulsatile and non-pulsatile cardiovascular model 15 Similar to [18] we predict Raup (and indirectly Ralp ), Cau and Em as timevarying quantities expressed using piecewise linear functions defined by X(t) = N X γi Ki (t), (25) i=1  t−t i−1  ,    ti − ti−1 ti+1 − t Ki (t) =     ti+1 − ti 0, ti−1 ≤ t ≤ ti ti ≤ t ≤ ti+1 otherwise where the unknown coefficients γi , i = 1, . . . , N are the new parameters to be estimated to predict the control. N is the number of nodes along the time span analyzed. As shown in Figure 5 gridpoints are equidistant during rest (in supine position) and after HUT (∆t = 5 s), but more dense during and immediately following the tilt up (and back down), where ∆t = 2 s. The actual tilt-time ( 14 s) is marked by a red marker on the Figure. 2.3 Model Summary The nonpulsatile model, summarized in equations (18-21) can be written on the form dx = f (x, t; θ), (26) dt where x = Vau , Val , Vvl , Vvu contains 11 parameters θ = {Raup , Ral , Rvl , Ralp , Cau , Cal , Cvl , Cvu , Em , EM , TR }. (27) This model is first solved during rest (before head-up tilt) with constant parameters and constant heart rate, i.e. the model will take the form (22) used 1.4 ∆t ti marker ti 1.2 1 0.8 0.6 0 20 40 60 80 100 120 140 160 180 200 time (s) Fig. 5: Position of the nodes before, during, and after head-up tilt. Nodes are distributed uniformly during sitting and standing (∆t15s), while node intervals are decreased to ∆t = 2 s immediately before and after the tilt. A similar pattern is used when the subject is tilted back down. The marker line indicates the duration of the tilt (∆t = 14 sec), noting how long it takes to tilt the subject up, a similar interval applies to the time when the subject is tilted back down. 16 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen to determine initial conditions. Subsequently, we solve the model during headup tilt, assuming that parameters controlled by the baroreflex system vary in time. For this study, we assume that the output data ydata can be described fully by the model plus an error term representing the measurement noise, i.e. ydata = g(x, t; θ) + ε, (28) where ymodel = g(x, t; θ) denotes the model output. For validation simulations we use surrogate data obtained from the pulsatile model for all states as well as arterial pressure, i.e. ydata = {Vau , Val , Vvl , Vvu , pau }. (29) For predictions we use the actual data, i.e. ydata = {pau , CO}. (30) 2.4 Model Analysis To analyze the proposed model we employed sensitivity analysis to determine the sensitivity of the model output to the parameters; subset selection to determine if sensitive parameters are correlated, and parameter estimation to minimize the least squares error between model predictions and data. Two approaches were used: a local approach analyzing the model in a small neighborhood around known parameter values and a global approach examining the model behavior over the entire parameter space. More detailed discussion of this type of approach can be found in our resent study [15]. The analysis was done using both for the surrogate model output defined in (29) and for actual model output defined in (30). Due to computational constraints, we only performed the analysis during steady state (before HUT) where the parameters are assumed constant. We justified this by noting that at any short time-interval (during the HUT) the parameters are approximately constant and therefore we expect that results found during rest can be transferred to HUT. Sensitivity analysis. For this study, we compute sensitivities of the model output ymodel to its parameters using a local approach based on the forward difference approximation ∂yk yk (t, θ̃ + δei ) − yk (t, θ̃) , = δ ∂ θ̃i where " i ei = 0 . . . 0 b 1 0 ... #T Pulsatile and non-pulsatile cardiovascular model Rel sens ranked 10 17 0 Sensitive parameters 10 -1 Insensitive parameters 10 -2 Pseudo data data E m C vu C au C al C vl R aup E M R alp R vl R al Fig. 6: Ranked sensitivities computed for surrogate and actual model output. Except for a few parameters, sensitivity ranking remain the same for the two models. Uncorrelated parameters are marked with black squares. Dashed squares indicate that the two parameters are correlated when tested with model output containing the actual data. √ is the unit vector in the i-th component direction, δ = ξ is the step-size, and ξ = 10−8 is the integration tolerance. Subsequently, we ranked the time-varying sensitivities by scaling the sensitivity matrix with the two-norm allowing us to separate parameters into two sets: sensitive and insensitive (Fig. 6). Subset selection. To determine potential parameter correlations, expected in a model derived from electrical circuit components with resistors and capacitors in series, we combined two approaches. First, we note that the model contains two parallel circuits predicting the flow in the upper and lower body. In steady state, the model can be reformulated as an equivalent circuit with one branch giving a priori information on parameter correlations. To eliminate these, we kept parameters associated with the lower extremities constant, while we analyzed parameters associated with the upper body θ = {Raup , Cau , Cvu , TR , Em , EM , xVlh,us }. Subsequently, we used singular value decomposition and QR factorization to identify a subset of parameters that can be identified given the model and available data [21]. For the surrogate data ydata = {Vau , Val , Vvl , Vvu , pau } three parameters are identifiable θ̂ρ = {Raup , Cau , Em }. For the actual data ydata = {pau , CO} the parameters Em and Cau are correlated, i.e. we can estimate two parameters either θ̂ρ = {Raup , Cau } or θ̂ρ = {Raup , Em }. Parameters, not included in the subset were kept constant at their nominal values. To verify this local result, only valid near the nominal parameter values, we tested the two subsets using a global analysis. For this analysis, we used the Delayed Rejection Adaptive Metropolis algorithm DRAM [7], an efficient Metropolis Hastings type Markov chain Monte Carlo (MCMC) algorithm that 18 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen estimates parameter distributions, which can be mapped pairwise to determine if two distributions are correlated. Results from this method (shown in Fig. 10), confirmed our local result that it is only possible to estimate two of the three parameters. Parameter estimation. To fit the model to data, we estimate identifiable parameters θ̃ ˆ θ̃ = arg min J(θ̃), (31) θ̃ where J = rt r denotes the least squares cost and r ( is the model residual given by T  1 yN − ỹ(tM ) y1 − ỹ(t1 ) y2 − ỹ(t2 ) , (32) r= √ , , ..., ỹ(t1 ) ỹ(t2 ) ỹ(tM ) M Here yi denotes the model output (ymodel ) while ỹ(ti ) denote the corresponding data (ydata ), and M is the length of the model output vector. To estimate the identifiable model parameters (θ̂ρ ), we used the LevenbergMarquart gradient method embedded in the Matlab function fmincon to minimize the least squares error, equivalent to maximizing the Gaussian likelihood function for the data. 3 Results The results are separated into three parts: validation of the pulsatile (P) and non-pulsatile (NP) models using data measured in supine position (before HUT), in response to HUT, and over the entire experiment (head-up tilt followed by head-down tilt). 3.1 Pulsatile and non pulsatile simulations (supine position) The first step involved estimating the subset of parameters θ̂ρ = {Raup , Cau , Em } minimizing the least squares error between pulsatile model predictions and measurements of arterial blood pressure pau and cardiac output (CO) at rest (before HUT). Second, we computed mean pressures and volumes from the pulsatile model, and estimated the same subset of parameters using the non-pulsatile model minimizing the least squares error between the non-pulsatile model predictions and the mean of predictions from the pulsatile model, i.e., the model residual r was given by " # P 1 xP i − xN i , x = {Vau , Val , Vvl , Vvu , pau } (33) r= √ M xP i where the superscript P, N P refers to the pulsatile and non-pulsatile models, respectively, see Figure 7. For these simulations 3 parameters were estimated. Pulsatile and non-pulsatile cardiovascular model 19 3.2 Pulsatile and non pulsatile predictions (HUT) The next set of simulations present HUT dynamics. First, we estimate the time-varying parameters for the pulsatile model. Nominal parameter values for this simulations are set using estimates from supine simulations ensuring that the models agree at rest. For these simulations, we estimate a total of 3 × N parameters where N is the number of nodes (see Figure 5). However, given that the nodes are distributed in time, at any time interval between nodes, only 3 parameters are estimated. For both pulsatile and non-pulsatile models, given that the locations of nodes are larger than δt = 0.02, the problem is well defined. To check convergence for both the pulsatile and non-pulsatile models, we added and removed every other node and checked that the solution remained invariant. As shown in Figure 8(a), upon HUT upper body arterial blood pressure drops due to pooling of blood in the legs. Results with estimated time-varying parameters are shown in Figure 8(b). From these results, we calculated the mean pressure and volumes used for predictions with the non-pulsatile model. It should be noted that the mean values obtained from model predictions are higher than the mean of the data (compare magenta and black lines on Figure 8(b). This can be explained from the fact that the curve computed from the model does not account for wave-reflection, and therefore has a larger area under the curve (see Figure 8(c)). For consistency between tests, all compar- 100 100 pau (mmHg) 90 80 70 60 50 80 70 60 0 20 40 50 60 P Pm NP 4.4 4.2 4 3.8 3.6 0 20 40 60 P Pm NP 4.4 4.2 pvl (mmHg) pvu (mmHg) P Pm NP 90 pal (mmHg) data P Pm NP 4 3.8 3.6 3.4 3.4 3.2 3.2 0 20 40 time (sec) 60 0 20 40 60 time (sec) Fig. 7: Prediction of blood pressure in supine position (before HUT). Upper body arterial blood pressure data (labeled data) (pau , black) is shown in the top left panel. Pulsatile model based predictions for the arterial (pau and pal ) and venous (pvu and pvl ) compartments are shown in light gray (labeled P ). The magenta lines (labeled Pm) depict the mean of the computed pulsatile pressure. Finally, the cyan lines (labeled NP ) show predictions from the non-pulsatile model minimizing the least squares cost (33). 20 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen isons between the pulsatile and non-pulsatile models were done against mean values of pulsatile predictions. Two sets of simulations were performed to compare pulsatile and nonpulsatile model predictions: first we show that the pulsatile and non-pulsatile models predict the same time-varying parameters when estimated against all states (the four arterial and venous volumes) as well as the upper body arterial pressure (pau ), minimizing (33). Second, we predict time-varying parameters minimizing the least squares error between predictions of pressure (pau ) and CO, mimicking data available. Results of these computations are shown in Figure 9. Time-varying parameters were determined from estimates of γRaup,i , γCau,i and γEm,i combined using the piecewise linear function (25). Results (Figure 9) show that the time-varying parameters agree for both formulations. However, this is not valid if the least squares cost function only includes arterial pressure pau and CO. For these, both formulations give similar predictions of Raup , while predictions of the other two variables do not agree, likely because with the reduced data, the three parameters can no longer be identified. To study potential parameter correlations, we conducted a DRAM simulation with 10,000 samples for the two models. Results showed that with a least squares cost including all volume states and pau , all parameters are identifiable, while simulations minimizing the pressure pau and CO, revealed a pairwise correlation between Cau and Em , see Figure 10. As a result of this analysis, we conducted pulsatile and non-pulsatile simulations estimating two parameters (Raup and Em ) and (Raup and Cau ), fixing the third (correlated parameter) at its nominal value. Results with the pulsatile model show that when excluding Cau , we were not able to capture pulse-pressure with the pul100 100 data P Pm pau (mmHg) 90 80 80 70 70 60 60 50 50 40 90 data P data mean Pm 90 70 60 40 0 50 100 time (sec) (a) 150 data P data mean Pm 80 0 50 100 time (sec) (b) 150 50 20 21 22 23 24 25 time (sec) (c) Fig. 8: Estimation of pulsatile dynamics during HUT. Panel (a) shows data (black) and model computations with nominal parameter values (gray) of the pulsatile upper body arterial pressure (pau ) during HUT. Panel (b) shows data (black) and model computations (gray) with time-varying parameter values of the pulsatile upper body arterial pressure (pau ) during HUT. The magenta lines (labeled P ) represent the mean of the pulsatile model computations in each panel, while the dashed black lines (labeled data mean) in panels (b) and (c) represent the mean of the pulsatile data. Note, the mean of the pulsatile model is higher than the mean of the data. This is because the model based predictions do not account for the reflected wave. Panel (c) shows a zoom over four cardiac cycles illustrating that the mean of the data is lower than the mean of the model based predictions. Pulsatile and non-pulsatile cardiovascular model 21 satile model (results not shown), while estimations with Raup and Cau allowed us to predict most features in the data, see Figure 11. 3.2.1 Pulsatile model with non pulsatile controlled quantities A goal of this study is to demonstrate that the simpler non-pulsatile cardiovascular model for HUT can be used in place of the more complex pulsatile model, 100 120 pau (mmHg) 80 pal (mmHg) P Pm NP NP PD 90 70 60 100 80 P Pm NP NP PD 60 50 40 0 50 100 40 150 0 50 100 (a) 25 5 pvu (mmHg) 4 3.5 3 2.5 20 pvl (mmHg) P Pm NP NP PD 4.5 2 1.5 150 (b) 15 10 P Pm NP NP PD 5 0 50 100 0 150 50 time (sec) 100 150 time (sec) (c) (d) 1.2 1.1 1 0.9 P NP NP PD 0.8 0 50 100 time (sec) (e) 150 Cau (ml/mmHg) Raup (mmHg s/ml) P NP NP PD 2.5 2 1.5 Em (mmHg/ml) 0.04 1.3 P NP NP PD 0.035 0.03 0.025 0.02 1 0 50 100 time (sec) (f) 150 0.015 0 50 100 150 time (sec) (g) Fig. 9: Pulsatile predictions during HUT with time-varying estimates of Raup , Cau and Em . Panels (a-d) include pulsatile model output (light gray labeled P ) and its mean (magenta curves labeled Pm), optimized non-pulsatile predictions (cyan curves labeled NP ) obtained by minimizing the least squares cost over all volumes and pau (33) as well as over pseudo data (pau and CO, (black curves labeled NP PD)). Panel (e) depicts computations of the upper body peripheral resistance for the pulsatile model computations (magenta curve labeled P ) and the non-pulsatile model predictions minimizing the cost function in (33) (cyan curve labeled NP ) and the cost function with pseudo data (30) (black curve labeled NP PD). Note that predictions of the peripheral vascular resistance Raup (e) are consistent over all simulations, while predictions of Cau and Em (f,g) vary, impacting predictions of pvu . 22 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen while still giving accurate dynamics and pertinent information to clinicians. The top row of Figure 12 shows model results obtained using the non pulsatile time-varying predictions of (Raup , Cau , Em ), estimated against all volumes and pau , in the pulsatile model. Results show that pulsatile predictions agree well for this set of parameters. Second, we used estimates of the time-varying parameters (obtained against the cost including pau and CO, mimicking data). Even though predictions of Cau and Em differ from those obtained with the pulsatile model, due to their inherent correlation, pulsatile predictions are 1.2 -0.02 Cau Raup -0.04 -0.06 0.8 -0.08 -0.1 2000 1 6000 10000 6000 10000 0.6 2000 6000 -0.1 -0.07 Raup -0.03 -4 -3.8 -3.6 Em -3.4 10000 0.6 0.8 1 Cau 1.2 Em -3.4 -3.6 -3.8 -4 2000 (a) (b) Raup Cau 1 0.8 0.6 Cau Em -3.5 -3.7 -3.9 -0.08 -0.06 -0.04 0.6 0.8 1 (c) Fig. 10: DRAM simulations with 10,000 samples estimating Raup , Em , and Cau with the non-pulsatile model (in supine position - before HUT). Dark gray predictions are obtained by minimizing the least squares error between the model predictions and the pseudo data (pau and CO), magenta estimates are obtained by minimizing over all volumes and pau . Note that for available data pau and CO parameters Em and Cau are correlated. This explains variation observed in predictions shown in Figure 9. (a) shows chains after burn-in predictions have been discarded, (b) shows parameter distributions, and (c) shows pair-wise correlations. Pulsatile and non-pulsatile cardiovascular model 23 still reasonable. Finally, we ran the pulsatile model with two non-pulsatile parameters (Raup and Cau ). Again, the close agreement between pulsatile and non-pulsatile predictions allowed us to predict pulsatile quantities. 3.3 Independent non pulsatile model: Estimation dynamics over the entire experiment To study dynamics during the entire experiment (over approximately 12 min), we only estimated parameters with the non-pulsatile model. To avoid problems with correlations between parameters, this simulation was done estimating the two independent parameters Raup and Cau . For this simulation, the objective was to capture ”long-term” dynamics. Results shown here were done using 181 nodes within the piecewise linear function (compared with 45 nodes for the HUT computations, estimating approximately 3 minutes of data). Results of this computation are shown in Figure 13. Note, during the initial HUT phase results agree well with those obtained earlier, see Figure 11. Again, note that predictions of the mean pressure data is excellent. 100 70 60 50 40 1.2 1.1 1 0.9 50 100 time (sec) (a) 150 2 1.5 1 0.8 0 P NP PD 2.5 1.3 Cau (ml/mmHg) pau (mmHg) 80 1.4 Raup (mmHg s/ml) data P Pm NP PD 90 0.7 P NP PD 0 50 100 time (sec) (b) 150 0.5 0 50 100 150 time (sec) (c) Fig. 11: Optimization of pulsatile (magenta curve labeled P ) and non-pulsatile (black curve labeled NP PD) models including time-varying predictions of the peripheral resistance Raup (b) and arterial compliance Cau (c). Panel (a) shows estimates with the pulsatile model (gray labeled P ) compared to data (black labeled data) and the non-pulsatile model (black curve labeled NP PD) compared to the mean of the pulsatile model output (magenta curve labeled Pm). Note since Em is kept constant the estimates immediately following HUT overestimates pulse pressure, yet predictions of time-varying parameters with the pulsatile and non-pulsatile models agree. pau (mmHg) 24 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen 100 100 100 90 90 90 80 80 80 70 70 70 60 60 50 50 60 50 40 40 0 40 0 50 100 150 0 50 pal (mmHg) (a) (b) 120 120 110 110 100 100 100 90 90 90 80 80 80 70 70 70 60 60 60 50 50 100 150 50 pvu (mmHg) 100 0 150 5 4 4 4 3 3 3 2 2 2 150 0 50 (g) 100 25 20 20 15 15 10 50 100 time (sec) (j) 50 150 0 100 150 100 150 25 20 15 10 5 5 150 (i) 10 0 0 150 (h) 25 100 (f) 5 100 50 (e) 5 50 150 50 0 (d) 0 100 (c) 110 50 50 150 120 0 pvl (mmHg) 100 5 0 50 100 time (sec) (k) 150 0 50 time (sec) (l) Fig. 12: Pulsatile predictions during HUT with time-varying estimates from the pulsatile (black) and non-pulsatile (gray) models. The rows show predictions for blood pressures pau , pal , pvu , and pvl . The first column (panels a, d, g and j) show predictions using non-pulsatile estimates of Raup , Cau , and Em obtained by minimizing over all volumes and pau (33). The second column (panels b, e, h, and k) show predictions minimizing over pseudo-data pau and CO, and the third column (panels c, f, i, and l) show predictions when only Raup and Em are estimated. 25 110 95 100 90 90 85 pal (mmHg) pau (mmHg) Pulsatile and non-pulsatile cardiovascular model 80 70 60 50 40 100 80 75 70 65 200 300 400 500 60 100 600 200 300 (a) 400 500 600 500 600 500 600 (b) 25 3.5 pvu (mmHg) pvl (mmHg) 20 15 10 3 2.5 5 0 100 200 300 400 500 2 100 600 200 300 (c) 400 (d) 100 1.1 Raup (ml/mmHg) CO (ml/s) 95 90 85 80 75 70 65 100 1 0.9 0.8 200 300 400 500 600 100 200 300 time (sec) 400 time (sec) (e) (f) Cau (ml/mmHg) 3.5 3 2.5 2 1.5 1 0.5 100 200 300 400 500 600 time (sec) (g) Fig. 13: Predictions over the full data with non-pulsatile model estimates of Raup and Em minimizing the least squares error between measured pressure data and CO. (a-d) shows predictions of pressure in the four compartments, (e) shows predictions of CO, (f-g) shows the time-varying parameter estimates. Stars indicate nodes in the piecewise function. These predictions were done using significantly fewer nodes as the validation computations during HUT, yet results follow observations with the more dense model (Figure 11). 26 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen 4 Discussion This study has shown that it is possible to develop interchangeable pulsatile and non-pulsatile models that can predict dynamics during HUT, and that time-varying parameters (Raup , Cau and Em ) can be predicted by both models when accounting for all volumes and pau . To our knowledge, this is the first study developing pulsatile and non-pulsatile models that can be represented by the same parameter set. Results here were obtained using a compartment model with heart rate as an input. HUT was imposed by accounting for gravitational pooling of blood in the lower extremities, and the autonomic response to HUT was included via time-varying parameters estimating vascular resistance, arterial compliance, and minimum elastance. The non-pulsatile model was compared with averages from the pulsatile model for a healthy young adult (see Figures 7 and 9). Finally, we showed that the two models are interchangeable, allowing parameter estimates obtained with the non-pulsatile model to be used within the pulsatile model (see Figure 12). The two models agree for parameters estimated by accounting for all volumes and pau . However, in a clinical setting measurements of blood volume are not available. As a result, we also estimated parameters minimizing the least squares cost only including quantities that can be measured clinically (upper body arterial pressure pau and cardiac output (CO)). Typically, HUT studies only measure heart rate and arterial blood pressure, but as discussed earlier [31] without measurements of CO predictions are not reliable. We recommend to include one measurement at baseline (in supine position) and one measurement after the HUT. For the results presented here, we did not have CO measurements. As a result, we estimated CO at rest using an allometric calculation combined with a 20% reduction during HUT as observed in literature [32, 33]. Results of these computations revealed that even though all three parameters can be distinguished with the pulsatile model, two of the three parameters are correlated. Therefore comparison over the entire experiment (tilt-up followed by tilt-down) was done estimating only two parameters. Physiologically, the autonomic control system does vary all three parameters even though cardiac contractility was kept constant in this study. One way to account for all variations would be to use information about correlation to built a model relating cardiac contractility and arterial compliance, or add additional measurements, e.g. of venous pressure. Results in this study were obtained using a piecewise linear function to estimate time-varying parameters. While this approach works [31, 16], it has the disadvantage that results critically depend on node placement. In a previous study [16], we used an ensemble Kalman filter to predict the time-varying parameters. The advantage of this method is that it in addition to the instantaneous value provide a measure of uncertainty, though uncertainty can also be calculated for the spline method, e.g. as discussed by Ramsay and Silverman [23]. This study did not use a filtering approach as it is difficult to set up a stable filter that can account for pulsatility within each period as well as variation in the response to the tilt. Future efforts will study if Kalman filtering Pulsatile and non-pulsatile cardiovascular model 27 can improve predictions presented here. An advantage of estimating the timevarying parameters using a piecewise linear function is that they encode total variation accounting for all mechanisms controlling the response to HUT. It is expected that the baroreflex system is engaged in regulating the controlled quantities, but it is likely that both cardiopulmonary receptors controlling effectors in response to changes in volume, and the respiratory sinus center also contribute to controlling the effectors [22]. While these controls are all encoded in the time-varying parameters, the estimated splines do not explain how each mechanism encode the response. One way to investigate the mechanistic relations is to build a model relating changes in hemodynamic quantities to changes in controlled parameters, a feature we plan to investigate in future studies. Finally, results presented here were obtained for a single healthy control subject for which we had measurement of pressure both at the level of the heart (where pressure increases during HUT) and at the level of the carotid arteries, where pressure drops in response to HUT. Typically, experimental HUT studies [11, 20] only measure blood pressure at the level of the heart. Due to the location with respect to center of gravity [31], to set up models predicting physiological response using this type of data requires further studies. One approach is to account for change the center of gravity to preprocess measured data, as suggested in our previous study [31]. With this type of preprocessing it should be possible to extend the methodology proposed here to analyze how parameters change within a large population allowing us to determine how controlled parameters change with disease. 5 Conclusion In summary, we have developed a non-pulsatile model interchangeable with a pulsatile model and shown that it can be used to predict HUT dynamics. We showed that two time-varying parameters can be identified by simulations analyzing typical experimental data including upper body arterial pressure and cardiac output (CO), while estimation of the third parameters require additional data, e.g. venous pressure, or development of a model explicitly describing the correlation between compliance and elastance. These models (the pulsatile and non-pulsatile models) have many potential benefits in applications that require analysis of data over long time-scales, they are simpler and faster to simulate than its pulsatile counterpart . Acknowledgements Results from this study was obtained with support from the Virtual Physiological Rat (VPR) project NIH-NIGMS#1P50GM094503, by the National Science Foundation under awards NSF-DMS #1022688, # 1557761, and the NCSU Research Training grant NSF-DMS#1246991. 28 Williams, Brady, Gilmore, Gremaud, Tran, Ottesen, Mehslen, Olufsen References 1. Batzel, J., Kappel, F., Schneditz, D., Tran, H.T.: Cardiovascular and respiratory systems: modeling, analysis, and control. SIAM, Philadelphia, PA (2007) 2. Beneken, J., Dewit, B.: A physical approach to hemodynamic aspects of the human cardiovascular system. In: E. Reeve, A. Guyton (eds.) Physical Bases of Circulatory Transport: Regulation and Exchange, pp. 1–45. W.B. Saunders, Philadelphia, PE (1967) 3. Boron, W., Boulpaep, E.: Medical Physiology, 3rd edn. Elsevier, Philadelphia, PA (2017) 4. Ellwein, L.: Cardiovascular and respiratory modeling. Ph.D. thesis, Department of Mathematics, NC State University, Raleigh, NC (2008) 5. Fink, M., Batzel, J., Kappel, F.: An optimal control approach to modeling the cardiovascular-respiratory system: an application to orthostatic stress. Cardiovsc Eng 4, 27–38 (2004) 6. Guyton, A., Hall, J.: Textbook of medical physiology, 13th edn. Elsevier, Philadelphia, PA (2013) 7. Haario, H., Laine, M., Mira, A., Saksman, E.: Dram: Efficient adaptive mcmc. Statistics and Computing 16, 339–354 (2006) 8. Heldt, T., Karmm, R., Mark, R.: Computational modeling of cardiovascular response to orthostatic stress. J. Appl. Phyisiol. 92, 1239–1254 (2002) 9. van Heusden, K., Gisolf, J., Stok, W., Dijkstra, S., Karemaker, J.: Mathematical modeling of gravitational effects on the circulation: importance of the time course of venous pooling and blood volume changes in the lungs. Am J Phsyiol 291, H2152–H2165 (2006) 10. Hoppensteadt, F., Peskin, C.: Modeling and simulation in medicine and the life sciences. Springer, New York, NY (2002) 11. Kenny, R., O’Shea, D., Parry, S.: The Newcastle protocols for head-up tilt table testing in the diagnosis of vasovagal syncope, carotid sinus hypersensitivity, and related disorders. Heart 83, 564–569 (2000) 12. Lanier, J., Mote, M., Clay, E.: Evaluation and management of orthostatic hypotension. Am Fam Physician 84, 527–536 (2011) 13. Lim, E., Chan, G., Dokos, S., Ng, S., Latif, L., Vandenberghe, S., Karunanithi, M., Lovell, N.: A cardiovascular mathematical model of graded head-up tilt. PLoS One 8, e77357 (2013) 14. Lim, E., Salamonsen, R., Mansouri, M., Gaddum, N., Mason, D., Timms, D., Stevens, M., Fraser, J., Akmeliawati, R., Lovell, N.: Hemodynamic response to exercise and headup tilt of patients implanted with a rotary blood pump: a computational modeling study. Artificial Organs 39, E24–E35 (2014) 15. Marquis, A., Arnold, A., Dean-Bernhoft, C., Carlson, B., Olufsen, M.: Practical identifiability and uncertainty quantification of a pulsatile cardiovascular model. Math Biosci 304, 9–24 (2018) 16. Matzuka, B., Mehlsen, J., Tran, H., Olufsen, M.: Using Kalman filtering to predict time-varying parameters in a model predicting baroreflex regulation during head-up tilt. IEEE Trans Biomed Eng 62, 1992–2000 (2015) 17. Melchior, F., Srinivasen, R., Charles, J.: Mathematical modeling of human cardiovascular system for simulation of orthostatic response. Am J Phsyiol 262, H1920–H1933 (1992) 18. Olufsen, M., Ottesen, J.: A practical approach to parameter estimation applied to model predicting heart rate regulation. J Math Biol 67, 39–68 (2013) 19. Olufsen, M., Ottesen, J., Tran, H., Ellwein, L., Lipsitz, L., Novak, V.: Blood pressure and blood flow variation during postural change from sitting to standing: model development and validation. J Appl Physiol 99, 1523–1537 (2005) 20. Parry, S., Reeve, P., Lawson, J., Shaw, F., Norton, M., Frearson, R., Kerr, S., Newton, J.: The newcastle protocols 2008: an update on head-up tilt table testing and the management of vasovagal syncope and related disorders. Heart 95, 416–420 (2009) 21. Pope, S., Ellwein, L., Zapata, C., Novak, V., Kelly, C., Olufsen, M.: Estimation and identification of parameters in a lumped cerebrovascular model. Math. Biosci. Eng. 6, 93–115 (2011) 22. Porta, A., Bassani, T., Bari, V., Tobaldini, E., Takahashi, A., Catai, A., Montano, N.: Model-based assessment of baroreflex and cardiopulmonary couplings during graded head-up tilt. Comp Biol Med 42, 298–305 (2011) Pulsatile and non-pulsatile cardiovascular model 29 23. Ramsey, J., Silverman, B.: Functional Data Analysis, 2nd edn. Springer, New York, NY (2002) 24. Shoemaker, W.: Fluids and electrolytes in the acutely ill adult. In: W. Shoemaker, S. Ayres, A. Greenvik, P. Holbrook (eds.) Textbook of critical care. W.B. Saunders, Philadelphia, PA (1989) 25. Silvani, A., Magosso, E., Bastianini, S., Lenzi, P., Ursino, M.: Mathematical modeling of cardiovascular coupling: Central autonomic commands and baroreflex control. Auton Neurosci 162, 66–71 (2011) 26. Starling, E., Visscher, M.: The regulation of the energy output of the heart. J Physiol 62, 243–261 (1926) 27. TenVoorde, B., Kingma, R.: A baroreflex model of short term blood pressure and heart rate variability. Stud Health Technol Inform 71, 179–200 (2000) 28. Ursino, M.: Interaction between carotid baroregulation and the pulsating heart: a mathematical model. Am J Physiol 275, H1733–H1747 (1998) 29. Ursino, M.: A mathematical model of the carotid baroregulation in pulsating conditions. IEEE Trans Biomed Eng 46, 382–392 (1999) 30. van de Vooren, H., Gademan, M., Swenne, C., TenVoorde, B., Schalij, M., der Wall, E.V.: Baroreflex sensitivity, blood pressure buffering, and resonance: What are the links? Computer simulation of healthy subjects and heart failure patients. J Appl Physiol 102, 1348–1356 (2007) 31. Williams, N., Wind-Willassen, O., Program, R., Mehlsen, J., Ottesen, J., Olufsen, M.: Patient specific modeling of head-up tilt. Math Med Biol 31, 365–392 (2014) 32. Zaidi, A., Benitez, D., Gaydecki, P., A.Vohra, Fitzpatrick, A.: Haemodynamic effects of increasing angle of head up tilt. Heart 83, 181–184 (2000) 33. Zhang, J., Critchley, L., Lee, D., Khaw, K.: The effect of head up tilting on bioreactance cardiac output and stroke volume readings using suprasternal transcutaneous doppler as a control in healthy young adults. J Clin Monit Comput 30, 519–526 (2015)