MODELING CARDIOVASCULAR AND RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE LAURA M. ELLWEIN1 , SCOTT R. POPE2 , ALIANG XIE3 , JERRY J. BATZEL4 , C.T. KELLEY5 , METTE S. OLUFSEN5

1) Olin Engineering Center, Marquette University, 1515 West Wisconsin Av, Milwaukee, WI 2) SAS, 100 SAS Campus Drive, Cary, NC 3) Pulmonary Physiology Laboratory, William S. Middleton Veterans Hospital, 2500 Overlook Tarace, Madison, WI 4) Inst. for Math. and Scientific Computing, Univ. of Graz, Heinrichstrasse 36, Graz, Austria 5) Dept. of Mathematics, North Carolina State University, Campus Box 8205, Raleigh, NC

Abstract. This study develops a coupled cardiovascular-respiratory model that predicts cerebral blood flow velocity (CBFV), arterial blood pressure, end-tidal CO2 , and ejection fraction for a patient with congestive heart failure. The model is a lumped parameter model giving rise to a system of ordinary differential equations. We use sensitivity analysis and subset selection to identify a set of model parameters that can be estimated given the patient data. Gradient based nonlinear optimization is used to estimate the subset of parameters. Optimization was caried out first for the cardiovascular submodel and subsequently for the respiratory model. Once a set of optimal parameters were found, the coupled model was computed to confirm that the model is still able to predict the observed data. Results showed that with the approach and methods presented in this paper it is possible to examine and quantify identifiability of model parameters. Using this approach we identified 5 key cardiovascular parameters and 4 key respiratory parameters. Nonlinear optimization techniques was used to estimate these parameters and we tested that values for all parameters were physiologically reasonable for a patient with congestive heart failure.

1. Introduction. The strong pathophysiological links between sleep-related breathing disorders and cardiovascular diseases have recently received raised attention among clinicians and researchers. On one hand, sleep-related respiratory disorders, especially obstructive apneas, have been identified as an independent risk factor for hypertension, myocardial infarction and stroke due to the associated recurrent upper airway vibration, intermittent hypoxia, surges of sympathetic neural activities, and repetitive arousals and resulting sleep fragmention a (see for example [23, 16]). Central sleep apnea has been considered a marker for congestive heart failure (CHF) and also been associated with the deterioration of heart function in CHF [21]. On other hand, epidemiological data reveal that about one-half of patients with CHF develop periodic breathing [17, 53]. CHF is characterized by a reduction in cardiac output, elevation of left atrial pressure and pulmonary capillary wedge pressure, accumulation of fluid in the lung, reduced functional residual capacity (FRC), and prolongation of circulation time. All of these factors may destabilize breathing. In fact, a greater central and peripheral chemosensitivity to CO2 (i.e., controller gain) has been observed in CHF patients [54, 9] and CHF animal models [57]. Respiratory changes duiring sleep combined with these cardiovascular and control sensitivity changes can contribute to the occurrence of sleep apnea ([12, 37, 4]. Furthermore, the low cardiac output also reduces brain perfusion [9, 36] and attenuates cerebrovascular chemosensitivity [22, 63]. The combination of hypoperfusion and low cerebrovascular reactivity impairs the CBF protective mechanism for breathing stability [64]. Finally, CHF and its concomitant cardiovascular insufficiency and potential for reduced cerebral perfusion has been associated with impaired cognitive function [5], while syncope is an issue in patients with CHF [26]. In general, although the Key words and phrases. Cardiovascular modeling, Respiratory modeling, Parameter estimation, Model validation, Sensitivity analysis, Subset selection. 1

2

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

underlying mechanisms for the interaction between the respiratory disorders and cardiaovascular disorders have not fully understood, the two disorders eventually exacerbate each other, creating a vicious cycle. Modeling and especially models that can be applied in the clinical setting can contribute to the development of methods for assessment of respiratory function such as estimating central and peripheral control gain and CO2 reactivity. Models can also provide insight into the control mechanisms that govern respiratory function, respiratory interaction with cardiovascular control, and mutual regulation between ventilation and cerebral blood flow (CBF) through partial pressure of arterial CO2 . Recent models that address aspects of the above topics include work by Dong and Langford [13], which developed a model to study factors affecting stable behavior of the cardiovascular-respiratory system in heart failure and Bidini et al. [38, 39] who applied a combined cardiovascular-respiratory model to study the Valsalva maneuver and factors affecting cerebral blood flow. Recent experimental studies include work investigating the differential role of partial pressures of CO2 and sympathetic response on vasoconstriction [1], and the role of cerebral vascular function in setting ventilatory response via influence on CO2 [63]. Thus, mathematical models that integrate ventilatory and cardiovascular system dynamics and control functions have the potential to be used as a diagnostic tool as well as to obtain more detailed knowledge of the dynamics of the underlying complex systems including aspects related to cardiovascular diseases. Furthermore, models that can be used to investigate the cardiovascularrespiratory coupling both in health and disease has potential to help understanding and diagnosing clinical conditions. To understand this complex dynamics models should be applicable for analyzing dynamics during typical orthostatic and ventilatory tests such postural change, step changes in inspiratory CO2 concentrations, or rhythmic breathing at various rates. These requirements indicate that models should include a fair degree of complexity. On the other hand, for a given model to have an impact in clinical settings, it is essential that the model can predict patient specific outcomes such as systemic and cerebrovascular resistances. The latter is often studied by solving an inverse problem, using a model combined with experimental data to predict patient specific parameters that minimize the difference between observed and computed quantities such as arterial blood pressure, cerebral blood flow velocity, or end-tidal CO2 . However, estimating parameters from complex nonlinear models is difficult, in particular since parameters may be correlated or not sensitive to given observations. Developing better techniques for analyzing such complex data is the main focus in this manuscript. To achieve an effective compromise between model complexity and clinical application requires a modular approach to model analysis. This entails several considerations. First, model reduction should be applied to reflect only the key functional elements under investigation. This simplifies the parameter estimation problem but care needs to be taken to avoid loss of the models’ predictive power. Second, a given model must be analyzed with regard to available data so that key model parameters sensitive to the observed data can be effectively and reliably estimated within a reasonable time frame. Other less sensitive parameters should be estimated using information from literature, statistics, or allometric scaling. Conversely, sensitivity analysis can provide suggestions for experimental design. However, sensitivity analysis alone does not fully enable estimation of all model parameters. Therefore, this analysis must be coupled with a more advanced analysis investigating correlations between model parameters, e.g., two correlated sensitive parameters should not both be estimated. This study develops a coupled pulsatile cardiovascular-respiratory model applicable to clinical problems related to CHF. We use sensitivity analysis to rank model parameters from the most to the least sensitive with respect to the available data. Furthermore, we study correlation of model parameters using subset selection.

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE Vcv Ccv

3

Cca Vca Rc

qc

Cerebral Tissue

cB,g

pcv

pca

MB,g

Head Rcv

Vpa Cpa

Cpv Vpv ppa

Lungs

qpa

qcv

Rpv

Vsv Csv

qsv

psv Rtv

pa,g

qpv

VD

ca,g

pD2,g pD1,g

prv Erv Vrv Right Ventricle

qsa

plv

pI,g

Vlv Elv

pexp,CO2

Rs

qca

Rmv

pD3,g

csv,g

Rca

ppv

qp

Rp

Left Ventricle

Csa Vsa psa

Rav

Thorax

qs

Systemic Tissue cS,g

Body

MS,g

Figure 1. Compartmental model of systemic and pulmonary circulations. Systemic (subscript s), cerebral (subscript c), and pulmonary (subscript p) arteries (subscript a) and veins (subscript v); Vessels carrying oxygenated blood (systemic arteries and pulmonary veins) are red, while vessels carrying deoxygenated blood (systemic veins and pulmonary arteries) are blue. All vascular compartments represent a group of vessels with similar pressure p [mmHg]. Each vascular compartment is characterized by its volume V [ml] and compliance C [ml/mmHg]. The left (dark red, subscript lv) and right (dark blue, subscript rv) ventricles generate pulsatile pressure plv , prv [mmHg]. These two compartments are defined using time-varying elastance E [mmHg/ml]. Flow q [ml/sec] between compartments are opposed by constant resistances R [mmHg sec/ml]. Tissue compartments (magenta) account for exchange of gases (O2 and CO2 ). Each vascular bed is characterized by a metabolic rate M [mlST P D /sec], and the gas concentrations in the tissues are denoted by c [mlST P D /ml]. The lungs are represented by three dead-space compartments each predicting the partial pressure of the gases pDi,g [mmHg]. The partial pressure of the gases in the inspired air pI,g , and the end tidal partial pressure of expiratory CO2 pexp,CO2 are marked separately.

2. Methods. The cardiorespiratory model depicted in Figure 1 is a lumped parameter compartmental model designed to predict systemic arterial blood pressure psa , cerebral blood flow velocity vc , and end tidal partial pressure of expiratory CO2 pexp,CO2 for a CHF patient resting in supine position. For the purposes of this paper, we are considering dynamic oscillatory variations in cardiovascular and respiratory system quantities around a resting steady state. Hence, control mechanisms are not responding to system stress and control responses are not directly invoked. For this reason we do not consider control submodels in this study. Note also that H and V˙ IE are inputs to the model so that variations in the controls of these quantities are implicitly included in the model.

4

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

2.1. Cardiovascular system. The cardiovascular system (see Figure 1) is represented by a closed circuit with 3 arterial compartments, 3 venous compartments, and 2 ventricular compartments. Vascular compartments represent systemic arteries and veins in the body and the brain, as well as pulmonary arteries and veins. Each compartment lumps vessels with similar transmural pressure p(t) [mmHg], volume V (t) [ml], and compliance C [ml/mmHg] (constant). For each compartment i, the stressed volume Vi,str (t) is given by Vi,str (t) = Vi (t) − Vi,unstr = Ci pi (t),

(1)

where Vi (t) [ml] is the time-dependent total volume and Vi,unstr [ml] (constant) is the unstressed volume, at zero transmural pressure. Flow between adjacent compartments are characterized by a constant resistance R [mmHg sec/ml] and volumetric flow rate q(t) [ml/sec]. Incorporating Ohm’s law, the net change in volume for each compartment is given by dVi (t) = qin (t) − qout (t), dt

qin (t) =

pi−1 (t) − pi (t) , Rin

qout (t) =

pi (t) − pi+1 (t) . Rout

(2)

The subscripts i − 1 and i + 1 refer to upstream and downstream compartments in relation to compartment i, respectively. A system of differential equations is obtained by differentiating (1) and equating with (2), dpi (t) qin (t) − qout (t) . (3) = dt Ci An equation of this form is associated with each vascular compartment. The heart is represented by the left and right ventricles, modeled to generate the driving pressures for the systemic and pulmonary systems. This is realized by imposing a time-varying pressure predicted as a function of ventricular volume. Thus, for the ventricular compartments we obtain a differential equation by imposing volume conservation, i.e., dViv (t) = qin (t) − qout (t), i = l, r (4) dt where the flows q are determined similarly to (2). We use a time-varying elastance model to predict ventricular pressures piv (t) [48, 51, 55, 60], piv (t) = Eiv (t)[Viv (t) − Vid ],

i = l, r,

(5)

where Vd [ml] denotes the volume at zero end-systolic pressure [48, 52] and Eiv (t) [mmHg/ml] denotes the time-varying elastance. In the above equations subscript i denote the left and right ventricles, respectively. For each ventricle, time-varying elastance E(t) is modeled using a piecewise sinusoidal function of the form     πt   /2 + Em 0 ≤ t ≤ TM (EM − Em ) 1 − cos   TM      (6) E(t) = (E − E ) cos π(t − TM ) + 1 /2 + E m TM ≤ t ≤ TM + TR M m   T  R   Em TM + TR ≤ t ≤ T,

where TM and TR denote the time for maximum (systolic) elastance (TM ) and the remaining time to relaxation (TR ). To account for varying heart rate, we express these times as fractions of the heart period T , i.e., TM,f rac = TM /T and TR,f rac = TR /T . We assume that peak elastance occurs at the same time for both left and right ventricles; thus we assume that the values of the two parameters TM and TR remain the same between the two sides of the heart. Em and EM denote that minimum and maximum elastance of each ventricle. The maximum and minimum elastance differ significantly between the right and the left ventricle, in particular for the CHF patient studied in this manuscript.

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

5

Thus we include four elastance parameters Eml , Emr and EMl , EMr and two timing parameters TM and TR . Similar to previous studies [45, 14, 50], the ventricular valves are modeled using time-varying resistances {Rmv (t), Rav (t), Rtv (t), Rpv (t)} defined as a function of the pressure drop across the valves. A small baseline resistance is used to define an “open” valve (subscript “o”) and a resistance that is several orders of magnitude larger is used to define the “closed” valve (subscript “c”). An exponential function is applied to describe the degree of openness as a function of the pressure gradient. The effective resistance of a valve is then defined by i h (7) Rvalve = min Rvalve,o + e−k(pin −pout ) , Rvalve,c .

The parameter Rvalve,o is the small resistance allowing flow out of the ventricle, k describes the speed of the transition from open to closed, and Rvalve,c is a value large enough to effectively shut off the flow through the valve. Since this function is non-smooth at the junctions of the exponential and Rvalve,c , a smoothing function [11] ! X exp(xi /ε) , min(x) = −ε log ε

i

where 0 < ε < 1 denotes the degree of smoothness (in this study we used ε = 0.5) and x denote the vector to be minimized. All equations and initial conditions for the cardiovascular model are given in the Appendix. 2.2. Respiratory system. The respiratory system is modeled using components allowing for gas exchange, transport, and metabolism. We model the dynamics of the two metabolite gases, O2 and CO2 , while we do not account for inert gases such as N2 . As shown in Figure 1, a tissue compartment connects each peripheral arterial compartment with each peripheral venous compartment. Dynamic gas concentrations in each tissue compartment reflect O2 consumption and CO2 production during metabolism as well as convection via blood flow into and out of each compartment. The model also includes two pulmonary components. The alveolar space has a dynamic volume, in which O2 and CO2 are exchanged between the lungs and the pulmonary vasculature. The lungs have been modeled using rigid dead space compartments connecting the alveolar space with the atmosphere. The rigid-walled dead space is divided into three compartments of the same size to account for some of the effects of the pulmonary branching. Gas concentrations and partial pressures are predicted in all compartments. The inputs to the model include HR and volumetric airflow V˙ IE , which provides information about tidal volume VT and ventilation frequency fR . Standard material balance equations describe respiration in the tissue compartments. The following symbol convention will be used: c represents concentration, T represents a generic tissue compartment which can be chosen to be systemic tissue S or brain tissue B. Generic gas quantity is denoted by g and gas fractional amount by F . Venous outflow is denoted by v and arterial inflow by a. Thus we denote the total amount of a gas in a tissue compartment by AT,g [ml] which is given by AT,g = VT,g cT,g . This equation describes AT,g by the product of the effective tissue volume V [ml] (constant) available for the gas and the concentration c [mlgas /mlblood ] of the gas in the volume. The change in amount of gas in a compartment is equal to the sum of the amount of gas produced or consumed by metabolism M [ml/sec] and the amount added or removed by the bloodstream qT [ml/sec]. dcT,g dVT,g dAT,g = cT,g + VT,g = MT,g + qT (ca,g − cT,g ). dt dt dt

6

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

Note, the concentration in the systemic arteries (i.e., in the body and the brain) are the same, since no metabolism has been accounted for on this side. Hence we denote the arterial concentration for ca . Thus, Assuming dVT,g /dt = 0, this equation reduces to dcT,g = MT,g + qT (ca,g − cT,g ). (8) dt The tissue compartments are considered well-mixed and equilibrated with exiting venous concentration. Note that the gas concentration csv,g is calculated as a mixture of the concentrations from the incoming systemic venous stream qs and cerebral stream qcv so that cS,g qs + cB,g qcv , csv,g = qs + qcv where cS,g and cB,g are the systemic and cerebral concentrations of each gas, respectively. Let A denote the alveolar compartment; then the quantity of alveolar gas is given by VT,g

VA,g = VA FA,g , which is the product of alveolar volume VA and gas fraction FA,g . The change in quantity of alveolar gas is represented by a mass balance relation parallel to the one developed for the tissue compartment, taking into account this time a time-varying alveolar volume, and with gas transport via fluid flowing to and from the lungs via pulmonary capillaries, and gas exchange with environment via inspiration and expiration. Let subscript p represent the pulmonary compartment. Thus we have, dVA,g dVA FA,g dVA = FA,g + VA = Fi,g + qp (csv,g − ca,g ). dt dt dt dt In this representation Fi,g denotes the fraction of gas in the air that is either being inspired or expired into the alveolar compartment. Hence, i = D3 during inspiration since inspired air is coming from the adjacent dead space region and i = A during expiration since air leaving the alveoli is alveolar air equilibrated with the pulmonary capillaries. Rearranging to express the rate of change of alveolar gas fraction gives   dVA dFA,g (FD3 ,g − FA,g ) + qp (csv,g − ca,g ), inspiration = VA dt q (c dt expiration. p sv,g − ca,g ), Gas fractions are converted to partial pressures via the relationship

FA,g = pA,g /(pamb − pwater ) = pA,g /713,

where pamb is the ambient air pressure of 760 mmHg and water vapor partial pressure pwater equals 47 mmHg at body temperature of 37◦ , i.e.,   dVA (pD3 ,g − pA,g ) VA dpA,g + qp (csv,g − ca,g ), inspiration = dt 713 q (c 713 dt expiration. p sv,g − ca,g ),

Because blood gas concentrations are reported in [mlST P D /ml] (STPD is the Standard Temperature (0◦ C), barometric Pressure at sea level and Dry gas), but alveolar volumes are in BTPS (standardized to Body Temperature, barometric Pressure at sea level, Saturated with water vapor: body temperature and pressure, saturated), we convert tissue concentrations to BTPS. Incoming air is immediately humidified once it enters the nasal passages [7, 8] and expired air is a composition of alveolar air and dead space air at BTPS, therefore terms with those quantities are not converted. The final equations are obtained by using the conversions and multiplying through by 713, i.e.,   dVA dpA,g (pD3 ,g − pA,g ) + 863 · qp (csv,g − ca,g ), inspiration (9) VA = dt 863 · q (c dt − c ), expiration. p

sv,g

a,g

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

7

Additional modeling considerations concern the relationship between the pulmonary capillaries and the systemic arteries. First, an anatomical shunt is present in the pulmonary circulation. This consists of O2 -deficient blood that bypasses the alveoli for physiological or pathological reasons. Typical shunt values range from 2-5%, with a larger shunt percentage possibly indicating a pathological condition [8, 7]. The subject studied has no known pulmonary health issues so we assume a 2% shunt. Therefore only 98% of the cardiac output becomes oxygenated, altering the alveolar equations to read as   dVA dpA,g (pD3 ,g − pA,g ) + 0.98 · 863qp (cv,g − ca,g ), inspiration = (10) VA dt 0.98 · 863q (c − c ), dt expiration. p v,g a,g

Second, we note that the thin alveolar wall allows for almost immediate equilibration of gases between the alveoli denoted by A and the pulmonary capillaries; thus we assume that the concentrations of blood gases is the same in the pulmonary capillaries and in the systemic arteries. This concentration is denoted by psa,g , i.e., we assume pA,g = psa,g as an auxiliary equation. Connected to the alveolar compartment are three compartments of equal volume with a total volume representing anatomical dead space. Each dead space is considered a well-mixed compartment with units VBT P S , in ml. Material balance equations for the dead space compartments reflect change in gas levels due to airflow, with opposite directions of flow for inspiration versus expiration. The relation FA,g = pA,g /713 holds, and equation units are all in BTPS. Thus following the derivation of (10), dpD1 ,g dt dpDi ,g VDi dt dpDi ,g Expiration: VDi dt dpD3 ,g VD3 dt

Inspiration: VD1

= = = =

dVA (pI,g − pD1 ,g ), dt dVA (pDi−1 ,g − pDi ,g ), i=2,3. dt dVA (pDi ,g − pDi+1 ,g ), i=1,2, dt dVA (pD3 ,g − psa,g ). dt

(11)

(12)

Pressure pI,g is the partial pressure of the gas in the inspired air. Note that dVA /dt is positive during inspiration and negative during expiration. The rate of change of alveolar volume dVA /dt is equivalent to the ventilation airflow V˙ IE . Thus, the alveolar volume can be predicted as Z (13) VA = V˙ IE dt. To give a correct value for the alveolar volume, additional information is needed to offset the integration constant. We have chosen to adjust this assuming that the minimum alveolar volume should match the FRC [ml], which we determine as a function of height H [cm] and weight W [kg], i.e., we assume min(VA ) = FRC,

FRC = (3.8 · H/100 − 3.41 · W/H − 2.74) · 1000

The above estimation for FRC is based on work by Stocks and Quanjer [56]. Furthermore, it should be noted that the use of the ”dot” notation in V˙ IE is a standard way to indicate a volumetric flowrate as opposed to a rate of change of volume. The airflow V˙ IE is measured during experimentation and is used as an input to the model to drive the gas concentration dynamics. O2 and CO2 have different affinities for hemoglobin, therefore behave differently in the gas versus liquid phases. Gas dissociation laws are used to convert alveolar gas pressures to blood gas concentrations.

8

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

In this study we use equations also used by Batzel et al. [4], cT,CO2

=

KCO2 pT,CO2 + kCO2 ,

cT,O2

=

K1 (1 − e−K2 pT ,O2 )2 .

Note that the law for CO2 is linear while the law for O2 is exponential. This reflects general behavior, but does not account for the the dependencies of each gas on the other (Bohr and Haldane effects) as well as factors such as pH and temperature. All respiratory equations and initial values are listed in the Appendix. 3. Experimental Methods. The data analyzed in this study include continuous (sampled at 128 Hz) measurements of systemic arterial blood pressure (psa ), cerebral blood flow velocity vc measured from the middle cerebral artery (MCA), and end tidal partial pressure of CO2 in expiratory air. In addition we have a measure for left ventricular ejection fraction, as well as anthropometric measurements of height, weight, and gender. These data are obtained from Dr. Skatrud’s group, Department of Medicine, University of Wisconsin. The subject is a male age 55, height 178 cm, weight 82.3 kg with CHF but no known sleep apnea syndrome. Data used for this study were recorded while the subject was in semirecumbent position during normal breathing (for 10 min). Experiments were started between 8.00 and 9.00 am, to standardize the effect of diurnal variability of cerebral vasomotor reactivity [2]. A 2-MHz pulsed Doppler ultrasound system (Neurovision 500 M; Multigon Industries, Yonkers, NY) was used to continuously measure cerebral blood flow velocity in the proximal segment of the middle cerebral artery (MCA). The MCA was insonated through the right temporal bone window using search techniques described in Otis and Ringelstein [46]. After detection and optimization of the Doppler signal, the probe was mechanically secured using a headband device and probe holder to provide a fixed angle of insonation for the duration of the experiment. The subject, in the semirecumbent position, was asked to keep his head still and eyes open throughout the experiment. Heart rate was obtained from the electrocardiogram, and arterial pressure (psa ) was measured beat by breat in the middle finger of the left hand by photoelectric plethysmography (Finapres, Ohmeda, Louisville CO). Tidal volume (VT ) and breathing frequency were measured with a pneumotachograph Model 3700, Hans Rudolph, Kansas City, MO) that was attached to a leak-free nasal mask. End tidal PETCO2 tensions was sampled from the mask and measured by gas analyzers (#S-3A/I & CD-3; Ametek, Pittsburgh PA). Detailed descriptions of the experimental protocol can be found in [63]. 4. Parameterization. 4.1. Cardiovascular parameters. Nominal parameter values and initial conditions for the cardiovascular model were estimated from the subject’s anthropometric measurements. A complete list of all initial parameter values can be found in Table 1. We start by setting up predictions for nominal parameter values used to predict time-varying elastance. Parameters needed include maximum and minimum elastance of the left ventricle (EMi , Emi , i = l, r), as well as fractions predicting the timing of the cardiac cycle (TM,f rac , TR,f rac ). Nominal parameter values for end diastolic elastance are obtained from estimated diastolic ventricular pressure pdia , end diastolic volume EDV, and zero pressure volume Vd as piv,dia , i = l, r. Emi = EDViv − Vid Diastolic ventricular pressures are set using literature values [29]. The subject studied has CHF and is thus expected to have an enlarged left ventricle, while the size of the right ventricle is closer to normal. We estimate these volumes using results for CHF patients without sleep-apnea reported by Tkacova et al. [59].

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

9

End systolic elastances are estimated using a similar relation, but as functions of systolic ventricular pressure psys , end systolic volume ESV, and zero pressure volume Vd . Systolic ventricular pressure is obtained from the maximum measured arterial pressure. End systolic volume is obtained by subtracting stroke volume (SV) from the end diastolic volume: ESV = EDV − SV. Stroke volume is estimated from end left ventricular diastolic volume and ejection fraction (EF, measured 26% for the subject studied): SV = EF · EDVlv . The timing fractions TM,f rac and TR,f rac were estimated from literature values suggested by Ottesen and and Danielsen [47] and Heldt [30]. We assume that the ejection of the left and the right sides of the heart are synchronous, thus we use the same values for both ventricles. Initial values for blood pressures were obtained partly from the data and partly from literature estimates. All arterial pressures were scaled relative to the measured arterial pressure, whereas venous pressures were set using standard literature values. Initial blood pressures were used both in computation of initial resistances and compliances and as initial conditions for the differential equations. Initial flows were scaled to cardiac output, which was computed from stroke volume SV and mean heart rate HR, i.e., CO = SV · HR. Mean heart rate was obtained from the measurements. Flows in the circuit was distributed to let 20% cardiac output go to the brain, while 80% was directed to the systemic arteries. Resistors were predicted using Ohm’s law as R=

pi − pi+1 qi

and compliances were predicted using the pressure volume relation (1). Estimates of total blood volume was computed as a function of body surface area (BSA, m2 ), estimated from Mosteller’s formula [42] p BSA = (W · H)/3600,

where W is the subjects weight in kg and H is the subjects height in cm. Total blood volume (ml) was computed using Baker’s formula [3] as Vt = (23.9 · BSA − 1.229) · 1000.

Distributions of volume and prediction of unstressed volumes were obtained using values from Beneken and DeWitt [6]. 4.2. Respiratory parameters. Nominal values for metabolic rates and tissue volumes are given in Table 2. As is standard practice, all tissue gas volumes and blood gas concentrations are given in STPD, thus units are consistent. Metabolic rates were set using standard allometric scaling proportional to body mass by the power of 3/4 [61, 62] using values in Table 2. We used metabolic rates given in [4, 34, 27] for total systemic (including both the body and the brain) CO2 (MCO2 ) and O2 (MO2 ) combined with a metabolic rate for the brain for CO2 (MB,CO2 ). To compute metabolic rates for oxygen in the brain MB,O2 we assumed that the ratio of metabolism between oxygen and carbon dioxide in the brain was approximately equal to one given that brain tissue burns primarily carbohydrate (glucose) ([27]). Finally, we used these four metabolic rates to adjust metabolic rates in the body. Similarly, for volumes of the gases we used total systemic volumes and brain volumes to get the volumes in the body, based on values given by Batzel [4]. Alveolar volume was a dynamic quantity as shown in (13). Dead space volume was set proportional to the body weight in pounds (the subject weight 82.3 kg or 181.3 lbs). Finally, the coefficients of the dissociation equations are assumed to be independent of the size of the subject.

10

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

Table 1. Nominal values for all cardiovascular parameters. Parameter EDVlv [ml] EDVrv [ml] EF SV [ml] ESVl ESVr HR [bpm] CO [l/min] pdsa [mmHg] plv,sys [mmHg] prv,sys [mmHg] plv,dia [mmHg] plv,dia [mmHg] ppa [mmHg] ppv [mmHg] psa [mmHg] psv [mmHg] pca [mmHg] pcv [mmHg] qs [ml/sec] qca , qc , qcv [ml/sec] qpa , qpv [ml/sec] Rmv , Rav [mmHg sec/ml] Rtv , Rpv [mmHg sec/ml] Rs [mmHg sec/ml] Rca [mmHg sec/ml] Rc [mmHg sec/ml] Rcv [mmHg sec/ml] Rp [mmHg sec/ml] Vsa [ml] Vsv [ml] Vca [ml] Vcv [ml] Vpa [ml] Vpv [ml] Csa [ml/mmHg] Csv [ml/mmHg] Cca [ml/mmHg] Ccv [ml/mmHg] Cpa [ml/mmHg] Cpv [ml/mmHg] Ac [cm2 ]

Physiologic description End diastolic volume (left ventricle) End diastolic volume (right ventricle) Ejection fraction Stroke volume End systolic volume (left ventricle) End systolic volume (right ventricle) mean heart rate Cardiac output Measured arterial pressure Systolic left ventricular pressure Systolic right ventricular pressure Diastolic left ventricular pressure Diastolic right ventricular pressure Pulmonary arterial pressure Pulmonary venous pressure Systemic arterial pressure Systemic venous pressure Cerebral arterial pressure Cerebral venous pressure Systemic flow Flow in cerebral circulation Flow in pulmonary circulation Mitral and aortic valve resistances Tricuspid and pulmonary valve resistances Systemic resistance Cerebral arterial resistance Cerebral resistance Cerebral venous resistance Pulmonary resistance Systemic artrial blood volume Systemic venous blood volume cerebral arterial blood volume Cerebral venous blood volume Pulmonary arterial blood volume Pulmonary arterial blood volume Systemic arterial compliance Systemic venous compliance Cerebral arterial compliance Cerebral venous compliance Pulmonary arterial compliance Pulmonary venous compliance Cerebral scaling factor

Value 312 100 0.26 EF · EDVl EDVl − SV EDVr − SV

Reference [59] [8, 29] measured

measured SV · HR max(pdsa ) 30 3 6 20 3.3 pdsa 6.6 psa · 0.99 7 0.8 CO 0.2 CO [8, 29] 0.001 0.001 (psa − psv )/qs (psa − pca )/qca (pca − pcv )/qc (pcv − psv )/qcv (ppa − ppv )/qp 0.1178 · Vt 0.6091 · Vt 0.237 · Vt 0.936 · Vt 0.0288 · Vt 0.1243 · Vt 0.3 · Vsa /psa 0.08 · Vsv /psv 0.22 · Vca /pca 0.08 · Vcv /pcv 0.58 · Vpa /ppa 0.11 · Vpv /ppv 0.3

measured measured [8, 29] [8, 29] [8, 29] [8, 29] [8, 29] measured [8, 29] measured [8, 29] [8, 29] [8, 29] estimated estimated

[6] [6] [6] [6] [6] [6] [6] [6] [6] [6] [6] [6] estimated

4.3. Parameter estimation. Model parameters discussed in the parameter above are obtained using physiological considerations and allometric scaling, so even though they to some extend have been adapted to the patient with CHF studied here, significant variations from these standard values is to be expected. Below we discuss how to estimate model parameters though the solution to the inverse problem: given data for cerebral blood flow velocity vc , arterial blood pressure psa , end tidal partial pressure of expiratory CO2 pexp,CO2 , and cardiac ejection fraction EF find the set of parameters that minimize the difference between computed and measured values of the observed quantities. The coupled model uses heart rate and airflow as inputs to predict these quantities. However, the model contains a large number of parameters (37) including 24 cardiovascular parameters and 13 respiratory parameters. The cardiovascular parameters include 9 resistances

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

11

Table 2. Nominal values for metabolic rates, tissue volumes, and gas dissociation constants. All values are adapted from [4, 34, 27]. Parameter MCO2 [mlST P D /sec] MO2 [mlST P D /sec] MB,CO2 [mlST P D /sec] MB,O2 [mlST P D /sec] MS,CO2 [mlST P D /sec] MS,O2 [mlST P D /sec] VT,CO2 [mlST P D ] VT,O2 [mlST P D ] VB,CO2 [mlST P D ] VB,O2 [mlST P D ] VS,CO2 [mlST P D ] VS,O2 [mlST P D ] VA,CO2 [mlBT P S ] VA,O2 [mlBT P S ] VD [mlBT P S ] K1 [mlST P D /ml] K2 [mmHg−1 ] KCO2 [mlST P D mmHg/ml] kCO2 [mlST P D /ml]

Physiologic description Systemic tissue metabolic rate of CO2 Systemic tissue metabolic rate of O2 Cerebral tissue metabolic rate of CO2 Cerebral tissue metabolic rate of O2 Systemic (body) tissue metabolic rate of CO2 Systemic (body) tissue metabolic rate of O2 Systemic tissue volume of CO2 Systemic tissue volume of O2 Cerebral tissue volume of CO2 Cerebral tissue volume of O2 Systemic tissue volume of CO2 Systemic tissue volume of O2 Alveolar tissue volume of CO2 Alveolar tissue volume of O2 Total dead space volume Dissociation coefficient for O2 Dissociation coefficient for O2 Dissociation coefficient for CO2 Dissociation coefficient for CO2

Value 4.333 5.167 0.875 MB,CO2 MCO2 − MB,CO2 MO2 − MB,O2 15000 6000 900 1000 VT,CO2 − VB,CO2 VT,O2 − VB,O2 3200 2500 181.1 0.200 0.046 0.0065 0.244

R = {Rs ,Rca ,Rc ,Rcv ,Rp ,Rav,o ,Rmv,o ,Rpv,o ,Rtv,o }, 6 compliances C = {Csa ,Cca ,Ccv ,Csv , Cpa ,Cpv }, 8 heart parameters θheart = {Vrd ,Vld ,Emr ,Eml ,EMr ,EMl ,TM,f rac , TR,f rac }, and a scaling factor Ac . The respiratory parameters include 3 metabolic rates M = {MCO2 ,MO2 ,MB,CO2 }, 4 gas tissue volumes V = {VT,CO2 , VT,O2 , VB,CO2 , VB,O2 }, the lung dead space volume VD , 4 dissociation constants K = {K1 ,K2 ,KCO2 ,kCO2 , }, and a fraction indicating alveolar air in the exiting air stream falv . Solving this inverse problem uniquely is not possible using this complex model, since it is likely that many combinations of these parameters can give rise to the same solution. Thus, the goal here is to identify a limited set of uncorrelated parameters that can be estimated reliably. To do so, we use sensitivity analysis to rank parameters from the most to the least sensitive, and we use subset selection to select a set of uncorrelated parameters. Subsequently, we use gradient based nonlinear optimization to estimate this reduced set of parameters, while the remaining parameters will be fixed at their nominal values. The data analyzed in this study comprise baseline values, thus it is essential that the solution of the system allows all states to reach “steady state” (in this study steady state means that the model should display steady oscillatory behavior). Timescales in the cardiovascular model are short, and this model will reach steady state within a few cycles. However, the respiratory model operates on much slower timescales. To analyze how long the model should run to reach steady state we analyzed timescales for each of the two systems. Time constants in the cardiovascular model are all of the same order of magnitude, while time constants in the respiratory model vary significantly. The largest time constant in the respiratory model is associated with the effective systemic tissue CO2 volume, which is approximately 14.1 liter. To compare the two models, we look at time constants related to cerebral arterial pressure pca and concentration of systemic tissue CO2 (cS,CO2 ). These time constants appear in the equations   pca − pcv 1 psa − pca dpca , − = dt Cca Rca Rc MS,CO2 + qs (ca,CO2 − cS,CO2 ) dcS,CO2 . = dt VS,CO2

12

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN 0.57

cS,CO2 [mlSTPD/ml]

0.565

0.56

0.555

0.55

0.545

0.54 0

1000

2000

3000 4000 time [sec]

5000

6000

7000

Figure 2. Dynamics of concentration of CO2 in the systemic tissue computed using nominal parameter value, sequence between the two vertical lines denotes the last sequence of the data. The time constant for these two equations are approximately (Cca Rca Rc )/(Rca + Rc ) ≈ 0.0167 seconds and VS,CO2 /qs ≈ 257 seconds, respectively. These estimates indicate that we need long sequences of baseline data to reach steady state dynamics. In this study we analyze 172 seconds of data, thus to ensure steady state dynamics we repeat the data-segment 40 times to reach an end time of approximately 7000 seconds. Figure 2 show dynamics of the tissue CO2 concentration. Another important issue is that each respiratory cycle contain approximately 4-5 cardiac cycles, requiring significantly smaller time steps to resolve cardiovascular dynamics than respiratory dynamics. Consequently, computation time for one simulation with the coupled cardiovascular-respiratory model is approximately 23 min on a 3 GHz Mac Pro with 4 Dual-Core Intel processors and 12 GB memory, rendering it computationally intensive to solve the parameter estimation problem using the coupled model. The cardiovascular and respiratory systems are weakly coupled. The respiratory model depends on blood flow through the pulmonary, systemic and brain tissues, while as stated the cardiovascular model does not depend on quantities from the respiratory model. Physiologically, the thoracic arteries and veins depend on respiration since change of volume in the lungs imposes oscillation of the external pressure in the tissue surrounding the thoracic arteries and veins. Such oscillations can be observed in measurements of arterial pressure, even on the finger pressure measurements used for analysis in this study. However, this slow oscillation is a secondary effect not included in the present study. Furthermore, when control of the system is included respiration will impact cerebral vascular resistance and compliance through autoregulation. Based on the above observations, we decoupled the models. The parameter estimation problem was solved by first estimating parameters for the cardiovascular model. Using the optimized cardiovascular parameters, systemic, pulmonary and cerebral mean blood flows were computed and used as inputs to the respiratory model. Finally, respiratory model parameters were estimated. A solution to the original coupled model was found using the optimized cardiovascular and respiratory parameter values and this solution was compared with the experimental data. This was done to ensure that the optimized parameters from the decoupled models could predict the dynamics of the complete model. 4.4. Least squares cost. As mentioned earlier, model parameters are estimated by minimizing the least squares error between computed and measured values of arterial pressure, cerebral blood flow velocity, and left ventricular ejection fraction (for the cardiovascular model) and end tidal partial

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

13

pressure of expiratory CO2 for the respiratory model. Similar to previous studies [50, 15] we note that to accurately predict cardiovascular dynamics it is important to account for changes in systolic and diastolic values of arterial pressure and cerebral blood flow velocity. For each timeseries, the data are measured at a fixed frequency (128 Hz); thus the quantities (pressure, velocity, partial pressure of CO2 ) are observed at N equally spaced times, ti . The timeseries analyzed contain a given number of cardiac cycles (M ), thus for arterial pressure and cerebral blood flow velocity we have M observations of systolic and diastolic values. In addition, we have an average value for ejection fraction (EF). Since the study is done during steady state, we assume that the ejection fraction should be the same for each of the M cardiac cycles. In summary, we have the following observations: • Arterial blood pressure psa (ti ), i = 1 : N [mmHg] • Cerebral blood flow velocity vc (ti ), i = 1 : N [cm/sec] • Systolic arterial blood pressure psa,sys,j , j = 1 : M [mmHg] • Diastolic arterial blood pressure psa,dia,j , j = 1 : M [mmHg] • Systolic cerebral blood flow velocity vc,sys,j , j = 1 : M [cm/sec] • Diastolic cerebral blood flow velocity vc,dia,j , j = 1 : M [cm/sec] • Left ventricular ejection fraction EFj , j = 1 : M = 0.26 • End tidal partial pressure of expiratory CO2 pexp,CO2 (ti ), i = 1 : N [mmHg]. Each of these data vectors has an associated predicted vector, which we denote by a superscript (p). p p The predicted vectors are ppsa , vcp , ppsa,sys , ppsa,dia , vc,sys , vc,dia , EF p , and ppexp,CO2 . Since the parameter estimation problem is solved separately for the cardiovascular and the respiratory models, we define two residual vectors Rcar and Rresp , as T  Rcar = Rpsa , Rvc , Rpsa,sys , Rpsa,dia , Rva,sys , Rva,dia , REF , T  Rresp = Rpexp,CO2 .

where each component is scaled to account for the number of elements in the vector and relative to the data, i.e., T  p ppi (tN ) − pi (tN ) pi (t1 ) − pi (t1 ) 1 , ,··· , Ri = √ pi (t1 ) pi (tN ) K where i = psa , vc , psa,sys , psa,dia , . . . , EF , and K = N, M , respectively. Thus, the vector Rcar has 2N + 5M entries, while Rresp has N entries. For each of the two models, the least squares cost J is defined by J = kRi k22 = RiT Ri ,

i = car, resp.

(14)

Minimizing the cost function is a nonlinear least squares problem. Solutions were found using the Levenberg-Marquardt variant of the Gauss-Newton optimization method [32]. Before solving the least squares problem we used sensitivity analysis and subset selection to prune the set of parameters. 4.5. Sensitivity analysis. Sensitivity analysis gives a measure of how much the output of a model is affected by changes in the model parameters. For the cardiovascular part of the model, sensitivities were computed for arterial blood pressure, ppsa , cerebral blood flow velocity vcp , and ejection fraction EF p . For each quantity x, the sensitivity to a parameter θ is a time-varying quantity, given by θ0 dxi for i = 1, . . . , N. dθ xi θ=θ0

To remove the effects of poorly scaled parameters and output, each component of the derivative is weighted by the parameter value and the reciprocal of predicted pressure. This gives sensitivities as dimensionless quantities, allowing them to be compared amongst each other.

14

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

To rank parameters from the most to the least sensitive we use a weighted 2-norm to obtain a sensitivity for each quantity. For the cardiovascular model, we compute sensitivities with respect to three quantities (psa , vc , and EF), then take the mean over the three quantities, i.e, √ √ √ Scar,θ = mean{kSppsa ,θ k2 / N, kSvcp ,θ k2 / N kSEF p ,θ k2 / M }. However, the respiratory model is only validated against one quantity (end tidal partial pressure of expiratory CO2 , i.e., the sensitivity is given by √ Sresp,θ = kSpexp,CO2 k2 / N . If the model is insensitive to a parameter, it should be not be considered during an optimization. 4.6. Subset Selection. Let θ be the vector of P model parameters. Let R be the model residual. Subset selection looks at the Jacobian, dR/dθ, and attempts to find a subset of columns that are “maximally independent”. We employ a method from [24] which also appears in [25]. The method was also used in [50, 15]. The method includes the following steps • Compute a Jacobian matrix at nominal parameter values, R′ (θ0 ) = dR/dθ|θ=θ0 , with error estimate ε. • Find a singular value decomposition of R′ (θ0 ) = U ΣV T with Σ = diag(σ1 , . . . , σP ) and singular values σ1 ≥ σ2 ≥ · · · ≥ σP . • Determine the number of columns, k, by counting the number of singular values larger than ε. • Use QR factorization with column pivoting [25] on VkT where Vk is the matrix consisting of the first k columns of V . The factorization is VkT P = QR. • The first k components of θT P = [θn1 · · · θnP ] should be used in the optimization.

Label the vector of k chosen parameters θk and the vector of remaining parameters θP −k . These steps ensure that the condition number of the reduced Jacobian, κ(dR/dθk ) = kdR/dθk k2 k(dR/dθk )† k2 remains small, and that dR/dθk does a good job of representing dR/dθP −k [25]. 5. Results. We first present results from the sensitivity analysis, then we discuss the subset selection algorithm and present simulation results obtained from the nonlinear least squares optimization. 5.1. Sensitivity analysis. For both the cardiovascular and the respiratory models we computed and ranked relative sensitives as described above. Figure 3 shows that both models contain both sensitive and insensitive parameters, but no clear jump separates the two groups. For both models we defined all parameters with a sensitivity greater than 10−2 as sensitive. This choice allows us to include all parameters that we can account for using physiological arguments. For the cardiovascular model this leaves 16 parameters including xcar,sens = {Vld ,Emr ,Rc ,Eml ,TM,f rac,Rp ,Csv ,Csa ,Cpv ,EMl ,EMr , Rs ,Vrd ,Rca ,Rcv ,Cca }. Note that all heart valve resistances are insensitive. From previous work [50, 15] we know that Ac is correlated with Rc , thus we keep this parameter constant and do not include this parameter in our list of parameters that we seek to estimate using nonlinear optimization. For the respiratory model, including parameters with sensitivity larger than 10−2 allows us to pick 6 parameters including xresp,sens = {MCO2 ,VD ,KCO2 ,kCO2 ,VT,CO2 ,VB,CO2 ,MB,CO2 }. Note, we do not include any data predicting partial pressure of O2 , thus all parameters related to prediction of O2 partial pressures/concentrations are insensitive.

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

15

Figure 3. Ranked sensitivities, top graph shows ranking for the cardiovascular model and bottom graph shows ranking for the respiratory model. For both graphs, blue squares denote sensitivities computed using the initial (nominal) parameter values, and the red diamonds denote sensitivities computed using the optimized parameters. Note, that the y-axis is a log scale. 5.2. Subset selection. For each subsystem, we investigated correlation between parameters using subset selection. For the cardiovascular model, subset selection allowed us to extract five uncorrelated parameters including xcar,sub = {Rc , Rs , Ca ,TM,f rac, EMl }. It should be noted that these are the same cardiovascular parameters found to be uncorrelated in our previous studies [50, 15]. Estimation of these five parameters using nonlinear optimization did not allow us to accurately predict the observed data. An analysis of the remaining parameters (not picked by subset selection) revealed that including a sixth parameter Cca allowed us to better predict the data. We remind the reader that subset selection was based on initial parameter values. We repeated the subset selection with optimized parameter values, and at this stage Cca was included in the subset. Furthermore, it should be noted, that if we ran the models using data from a young subject, the same cardiovascular parameters were identified, but for the young subject it was not necessary to also estimate arterial cerebral compliance. Before predicting sensitivities for the respiratory model we used optimized cardiovascular parameters to compute the mean flow in the systemic q s , cerebral q c , and pulmonary q p circulations. Along with airflow data these flows were used as inputs to the model. These quantities were kept fixed during the subset selection. Results from subset selection allowed us to estimate 4 parameters including xresp,sub = {KCO2 , kCO2 , MS,CO2 , VD }. It should be noted that none of the selected parameters were insensitive. 5.3. Nonlinear optimization. For both models parameters were estimated using the LevenbergMarquardt variant of the Gauss-Newton optimization method [32]. Optimized parameters are given in Table 3. For the cardiovascular model data were subsampled at 64 Hz to speed up simulations, and computations settled at steady state within the data analyzed. For the respiratory model data were subsampled at 4.2667 Hz to speed up simulations, furthermore to ensure that steady state was reached respiratory data were repeated 40 times. The comparison used to compute the least squares error were done over the last repeat of the data. Results of optimizations are shown in Figures 4

16

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

Table 3. Optimized model parameters, note the value for TM,f rac is the fraction relative to the length of the cardiac cycle. Averaged model quantities include ejection fraction, cardiac output, mean heart rate, as well as systemic, cerebral, and pulmonary flow. Parameter Rs Rc Csa Cca EM l TM,f rac EF CO HR qs qc qp KCO2 kCO2 MS,CO2 VD

Description Systemic resistance Cerebral resistance Systemic arterial compliance Cerebral arterial compliance Maximum left ventricular elastance Time for maximum ventricular elastance Ejection fraction Cardiac output mean heart rate Mean systemic flow Mean cerebral flow Mean pulmonary flow Dissociation constant Dissociation constant Metabolic rate systemic tissue Total deadspace volume

Optimized 1.28 2.54 0.536 2.96 0.721 0.139 0.262 4.92 62.5 55.0 26.9 82.0 0.00688 0.261 4.20 151

Unit mmHg sec/ml mmHg sec/ml mmHg/ml mmHg/ml ml/mmHg

l/min beats/min ml/sec ml/sec ml/sec mlST P D mmHg/ml mlST P D mmHg/ml mlST P D mmHg/sec ml

and 5. The cardiovascular model used heart rate as an input, and the respiratory model used airflow and mean blood flows computed from the cardiovascular model with optimized parameters. 6. Discussion. Results from the cardiovascular and respiratory models (Figures 4, 5, and 6) showed that the modular model approach used in this study enable prediction of both cardiovascular and respiratory quantities. Results from the coupled model were indistinguishable from the results obtained with each of the two models indicating that using mean values for flows do not impact the respiratory model significantly, which is to be expected since respiration occurs over a much slower timescale as the cardiovascular oscillations. In addition to estimating observed quantities, we also computed all internal states and most of these were within physiological bounds as discussed below. It should be noted that techniques used for parameter estimation do not guarantee that internal states cannot fluctuate away from physiological values, even though initial parameter estimates were within physiological bounds. This is one of the main problems with the proposed method. Besides including bounds on all model parameters, we did not impose additional constraints limiting the states within certain bounds. Consequently some states did drift away from expected values. In regards to cardiovascular characteristics as reflected by the model for this subject with CHF, we note that approximately 30% of cardiac output was utilized as cerebral blood flow. This is large, however, cerebral blood flow velocity measured for this subject was significantly larger than mean velocities reported in the literature [28]. This could be due to poor circulation in the body or exceptionally high values from the transcranial Doppler measurements. This subject exhibits a severely enlarged left ventricle (consistent with left ventricular failure) and reduced pumping effectiveness, with an ejection fraction of 26%. Neverless, the cardiac output (modeled) is close to normal. This type of behavior is seen in approximately half the patients suffering from CHF [43]. With optimized parameters the maximum left ventricular volume was approximately 317 ml, which is close to initial volume chosen as 312 ml. The minimum ventricular volume was approximately 238 ml, yielding a stroke volume of 79 ml and with the mean heart rate of 62.5 beats/min this gives an approximate CO of 4.9 computed from the flows in the model (see Table 3). Systemic venous pressures were approximately 7.2 mmHg somewhat higher than normal, which is expected for a patient with CHF (see Table 4). Pulmonary arterial pressures were approximately 23 mmHg

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

100 80 100

vc [cm/sec]

100

Data Model

80 60 40 0

50

100 time [sec]

150

psa [mmHg]

120

Data Model

100 80 60 150

152

154

156

158

vc [cm/sec]

100 80 60 40 150

152

154 156 time [sec]

0.32

0.28 0.27

Data Model

80 50

100 time [sec]

150

Data Model

100 80 152

154

156

158

160 Data Model

100 90 80 152

154 156 time [sec]

0.28

158

160

Model Data

0.275 0.27 0.265 0.26 0.255

0.26 0.25

150

90

150

Ejection Fraction (EF)

0.29

100

100

60 150

160

0.3

50

120

Model Data

0.31 Ejection Fraction (EF)

158

80

0

160 Data Model

Data Model

100

60 0

150

vc [cm/sec]

50

psa [mmHg]

60 0

120 psa [mmHg]

Data Model

vc [cm/sec]

psa [mmHg]

120

17

50

100 Heartbeat

150

0.25

50

100 Heartbeat

150

Figure 4. Simulation results cardiovascular model. Top panel shows arterial pressure and cerebral blood flow velocity for the full dataset, second panel shows a zoom for 150 ≤ t ≤ 160 seconds. Bottom panel shows computed and expected ejection fraction (set at 26%). Left graphs give results with initial (nominal) parameter values, while the right graphs show results with the optimized values.

18

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

35

35

30

30 pexp,CO2 [mmHg]

40

pexp,CO2 [mmHg]

40

25

25

20

20

15

15

10

10

5

5

0 0

50

100 time [sec]

0 0

150

35

35

30

30

25

150

25

20

20

15

15

10

10

5

5

0 100

100 time [sec]

pexp,CO2 [mmHg]

40

pexp,CO2 [mmHg]

40

50

110

120 130 time [sec]

140

150

0 100

110

120 130 time [sec]

140

150

Figure 5. Estimation of respiratory parameters minimizing difference between measured and computed values of partial pressure of CO2 . Left graphs show computations with initial parameters and right graphs show results with optimized parameters. arterial exp

40 35

pCO2 [mmHg]

30 25 20 15 10 5 0 100

110

120 130 time [sec]

140

150

Figure 6. Model simulation of psa,CO2 and pexp,CO2 showing relation between these quantities during inhalation and exhalation

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

19

and pulmonary venous pressures were approximately 5.5 mmHg. This is also to be expected for a patient with left ventricular CHF (see Table 4). Heart rate (not modeled) was within normal range. For this model both cerebral and systemic vascular resistance and compliance were kept constant. The total resistance RT = (1/RB + 1/Rs ) = 0.86, where RB = Rca + Rc + Rcv is fairly low, again this is to be expected in CHF where CO is lower than normal, and for this subject in particular, where a large quantity of the blood goes to the brain. Moreover, the somewhat large proportion of flow to the brain could also be a result of the somewhat simplified modeling included for the cerebral arteries. The blood flow velocity measured is from the middle cerebral arteries, this has then been scaled (assuming constant area of the major cerebral arteries) to compute a total flow to the brain. Without information about anatomical quantities for the size of the cerebral arteries, this scaling factor cannot be computed accurately, and the magnitude of the scaling factor Ac would impact the optimized values for the cerebrovascular resistance Rc since these two parameters are correlated as stated earlier. However, if e.g., the area of all major cerebral vessels (PCA, MCA and ACA) were measured e.g., with an angiogram, then this scaling parameter could be computed from data. For the respiratory quantities values for dead space pD,CO2 and pD,O2 , arterial values pa,CO2 and pa,O2 , and systemic values pS,CO2 and pS,O2 (calculated from dissociation laws) are reasonable. Dynamics of pA,CO2 (equivalent with pa,CO2 ) are shown in Figures 5 and 6. Dead space values for CO2 oscillate between zero (inhaled air) and end-expiratory air approximately the same as alveolar air. Oxygen values follow similar reasonable patterns. It should be noted that brain partial pressures pB,CO2 and pB,O2 did not match typical values so well in the optimized case. Estimates of the dissociation constants varied from the initial guesses which resulted in pB,CO2 approximately 43 mmHg (too low) and pB,O2 approximately 48 mmHg (too high). There are likely several reasons for this. Firstly, the data for cerebral blood velocity for this subject was rather high, which would tend to depress brain pB,CO2 and raise pB,O2 . Secondly, we did not have data on O2 which might impact the results. Further, research indicates that normal cerebral metabolism is distorted in CHF and cerebral energy deficits can exist in certain cerebral tissue [35]. This complicates the picture for the estimation process. These observations lead us to conclude that the proposed methodology and model can be used for quantitative prediction of changes within several subjects with CHF. To use this model to compare quantities from a larger population it should be emphasized that the model should be used with care. It is important that nominal parameters reflect known characteristics from a subject with CHF, in particular since only a subset of the parameters will be estimated. Conclusions based on analysis of parameters obtained from several subjects should all be interpreted related to the nominal parameter values used. Furthermore, additional measurements that would increase predictability of the model parameters include measurements of cardiac output and estimation of diameters of major cerebral vessels. Knowing these quantities would allow better prediction of systemic and arterial flow as well as help to get better estimation of the scaling factor needed, which would help in prediction of the scaling factor relating MCA CBFV to CBF. Furthermore, if constants used in dissociation laws (optimized with our respiratory model) are fixed, subset selection allows estimation of the tissue and brain CO2 metabolic rates as well as the deadspace volume. The parameter values for deadspace volume were identical to simulation reported here, while parameters for tissue metabolic rate went down slightly (from 4.20 to 4.16) and the parameter for brain CO2 metabolic rate was reduced from its initial value (from 1.04 to 0.87). Furthermore, simulations with this set of parameters leads to similar dynamics that were observed in results presented above. Again, it should be noted that it is essential that caution is used when choosing parameters fed to the subset selection. Subset selection will give a subset of uncorrelated parameters that can be estimated given particular nominal parameter values and initial conditions for the differential equations.

20

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

Table 4. Generic state changes in left systolic CHF. State Mean pap [mmHg] Mean pvp [mmHg] Mean pas [mmHg] Mean pvs [mmHg] Mean plv [mmHg] Mean pla [mmHg] LVEDP [mmHg] LVESP [mmHg] LVEDV [ml] LVESV [ml] LVEF [ml] CO [mlmin] SV [ml] HR [beat/min] Systemic resistance [mmHg · s /ml] Fluid retention [ml]

Value increase increase decrease or stable increase steady to increase increase increase steady or increase increase increase decrease decrease decrease increase increase increase

Source [44] [40, 31] [49] [31] [41] [49] [41] [41] [10, 49] [10, 49] [10, 20] [49] [49] [49, 20] [31] [31, 19]

Acknowledgments. Olufsen and Ellwein were partially supported by National Science Foundation under grant DMS-0616597 and OISE-524249, Batzel was suported in part by FWF (Austrian Research Funds) project P18778-N13, Pope and Kelley was supported in part by National Science Foundation under grant DMS-0707220, and Xie was supported in part by the American Lung Association of Wisconsin.

Appendix.

6.1. Cardiovascular model. Blood pressures, ventricular volumes, and intrathoracic pressure. dppa dt dppv dt dpsa dt dpsv dt dpca dt dpcv dt dVlv dt dVrv dt

= = = = = = = =

 ppa − ppv prv − ppa /Cpa − Rpv Rp   ppv − plv ppa − ppv /Cpv − Rp Rmv   psa − psv psa − pca plv − psa /Csa − − Rav Rs Rca   pcv − psv psv − prv psa − psv /Csv + − Rs Rcv Rtv   psa − pca pca − pcv − /Cca Rca Rc   pcv − psv pca − pcv /Ccv − Rc Rcv plv − psa ppv − plv − Rmv Rav prv − ppa psv − prv − . Rtv Rpv 

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

21

Table 5. Initial values for the cardiovascular model. Note initial values for the cardiovascular system are set using the corresponding nominal parameter values given in Table 1. Variable ppa [mmHg] ppv [mmHg] pdsa [mmHg] psv [mmHg] pca [mmHg] pcv [mmHg]

Physiologic description Pulmonary arterial pressure Pulmonary venous pressure Systemic arterial pressure Systemic venous pressure Cerebral arterial pressure Cerebral venous pressure

Value 20 3.3 pdsa 6.6 psa · 0.99 7

Reference [8, 29] [8, 29] measured [8, 29] measured [8, 29]

Vlv [ml] Vrv [ml]

Left ventricular volume Right ventricular volume

EDVlv = 312 EDVrv = 100

[59] [8, 29]

Note the airflow V˙ IE is defined below. Heart valves Rav Rmv Rpv Rtv

i h = min Rav,o + e−2(plv −psa ) , Rav,c i h = min Rmv,o + e−2(ppv −plv ) , Rmv,c h i = min Rpv,o + e−2(prv −ppa ) , Rpv,c i h = min Rtv,o + e−2(psv −prv ) , Rtv,c .

Ventricular pressures plv (t) = prv (t) =

Elv (t)[Vlv (t) − Vld ]

Erv (t)[Vrv (t) − Vrd ].

where

Elv (t) =

Erv (t) =

    πt   Eml + (EMl − Eml ) 1 − cos /2   TM      π(t − TM ) + 1 /2 Eml + (EMl − Eml ) cos   TR    Eml     πt   /2 E + (E − E ) 1 − cos mr Mr mr   TM      π(t − TM ) + 1 /2 E + (E − E ) cos mr Mr mr   TR    Emr

0 ≤ t ≤ TM TM ≤ t ≤ TM + TR TM + TR ≤ t ≤ T, 0 ≤ t ≤ TM TM ≤ t ≤ TM + TR TM + TR ≤ t ≤ T.

Note, the timing parameters TM and TR are the same for both the left and right ventricles. Input to this model is heart rate HR. Output from the model is arterial blood pressure psa (t) and cerebral blood flow velocity vc = (pca − pcv )/Rc . Initial values for these differential equations are listed in Table 5 (below), while nominal values for all model parameters are listed in Table 1 and optimized values for estimated parameters are given in Table 3.

22

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

6.2. Respiratory model. dcS,CO2 dt dcS,O2 dt dcB,CO2 dt dcB,O2 dt

=

MS,CO2 + qs (ca,CO2 − cS,CO2 )/VS,CO2

=

−MS,O2 + qs (ca,O2 − cS,O2 )/VS,O2

=

MB,CO2 + qc (ca,CO2 − cB,CO2 )/VB,CO2

=

−MB,O2 + qc (ca,O2 − cB,O2 )/VB,O2 .

Inspiration dpD1,CO2 dt dpD1,O2 dt dpD2,CO2 dt dpD2,O2 dt dpD3,CO2 dt dpD3,O2 dt dpa,CO2 dt dpa,O2 dt

= V˙ ie (pi,CO2 − pD1,CO2 )/VD1 = V˙ ie (pi,O2 − pD1,O2 )/VD1 = V˙ ie (pD1,CO2 − pD2,CO2 )/VD2 = V˙ ie (pD1,O2 − pD2,O2 )/VD2 = V˙ ie (pD2,CO2 − pD3,CO2 )/VD3 = V˙ ie (pD2,O2 − pD3,O2 )/VD3 = 863 · 0.98 · qp (cv,CO2 − ca,CO2 ) + V˙ ie (pD3,CO2 − pa,CO2 )/VA = 863 · 0.98 · qp (cv,O2 − ca,O2 ) + V˙ ie (pD3,O2 − pa,O2 )/VA .

Expiration dpD1,CO2 = V˙ ie (pD1,CO2 − pD2,CO2 )/VD1 dt dpD1,O2 = V˙ ie (pD1,O2 − pD2,O2 )/VD1 dt dpD2,CO2 = V˙ ie (pD2,CO2 − pD3,CO2 )/VD2 dt dpD2,O2 = V˙ ie (pD2,O2 − pD3,O2 )/VD2 dt dpD3,CO2 = V˙ ie (pD3,CO2 − pa,CO2 )/VD3 dt dpD3,O2 = V˙ ie (pD3,O2 − pa,O2 )/VD3 dt dpa,CO2 = 863 · 0.98 · qp (cv,CO2 − ca,CO2 )/VA dt dpa,O2 = 863 · 0.98 · qp (cv,O2 − ca,O2 )/VA . dt Note, instantaneous airflow is obtained from integrating the air flow velocity data as described in equation (13). Inputs to the stand alone respiratory model are airflow velocity, average systemic, pulmonary, and cerebral blood flow obtained by solving the cardiovascular model. For the coupled model, these flows

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

23

are computed at any instant in time. Initial values for the respiratory model are given in Table 6, while nominal parameters for all respiratory parameters are listed in Table 2, and optimized values for estimated parameters are given in Table 3. Table 6. Initial values for the respiratory model. Variable pD1,CO2 [mmHg] pD1,O2 [mmHg] pD2,CO2 [mmHg] pD2,O2 [mmHg] pD3,CO2 [mmHg] pD3,O2 [mmHg] pa,CO2 [mmHg] pa,O2 [mmHg] cS,CO2 [mlSTPD /ml] cS,O2 [mlSTPD /ml] cB,CO2 [mlSTPD /ml] cB,O2 [mlSTPD /ml]

Physiologic description CO2 partial pressure lung deadspace 1 O2 partial pressure lung deadspace 1 CO2 partial pressure lung deadspace 2 O2 partial pressure lung deadspace 2 CO2 partial pressure lung deadspace 3 O2 partial pressure lung deadspace 3 Systemic arterial CO2 partial pressure Systemic arterial O2 partial pressure Systemic tissue CO2 concentration Systemic tissue O2 concentration Cerebral tissue CO2 concentration Cerebral tissue O2 concentration

Value 5 159 6 158 7 157 40 100 0.543 0.128 0.569 0.112

Reference [27, 33] [27, 33] [27, 33] [27, 33] [27, 33] [27, 33] [58] [18] [18] [18] [27, 33] [27, 33]

REFERENCES [1] P.N. Ainslie, J.C. Ashmead, K. Ide, B.J. Morgan, and M. J. Poulin, Differential responses to CO2 and sympathetic stimulation in the cerebral and femoral circulations in humans, J Physiol, 566 (2005), 613–624. [2] S. F. Ameriso, J. G. Mohler, M. Suarez, and M. Fisher, Morning reduction of cerebral vasomotor reactivity, Nerology, 44 (1994), 1907–1909. [3] R. J. Baker, D. D. Kozoll, and K. A. Meyer, The use of surface area as a basis for establishing normal blood volume, Surg Gynecol Obstet, 104 (1957), 183–189. [4] J. J. Batzel, F. Kappel, and S. Timischl-Teschl, A cardiovascular-respiratory control system model including state delay with application to congestive heart failure in humans, J Math Biol, 50 (2005), 293-335. [5] S. J. Bennett, M. J. Sauv¨ e, and R. M. Shaw, A conceptual model of cognitive deficits in chronic heart failure, J Nurs Scholarsh, 37 (2005), 222–228. [6] J. Beneken and B. DeWitt, A physical approach to hemodynamic aspects of the human cardiovascular system, in “Physical Bases of Circulatory Transport: Regulation and Exchange” (eds. E. Reeve and A. Guyton), W. B. Saunders, Philadelphia, PA, 1967, 1–45. [7] R. M. Berne and M. N. Levy, Physiology, C. V. Mosby Company, Second edition, 1988. [8] W. Boron and E. Boulpaep, Medical Physiology, W. B. Saunders, Philadelphia, PA, 2003. [9] S. N. Caparas, M.J Clair, R. S. Krombach, J. W. Hendrick, W. V. Houck, S. B. Kribbs, R. Mukherjee, G. E. Tempel, F. G. Spinale, Brain blood flow patterns after the development of congestive heart failure: effects of treadmill exercise, Crit Care Med, 28 (2000), 209–214. [10] K. Chatterjee and B. Massie. Systolic and diastolic heart failure: differences and similarities, J Card Fail, 13 (2007), 569–576. [11] X. Chen, L. Qi, and K-L. Teo, Smooth convex approximation to the maximum eigenvalue function, J Global Opt, 30 (2004), 253–270. [12] N. S. Cherniack, Apnea and periodic breathing during sleep, New England J Med, 341 (1999), 985–987. [13] F. Dong and W. F. Langford, Models of Cheyne-Stokes respiration with cardiovascular pathologies, J Math Biol, 57 (2008), 497–519. [14] L. Ellwein, H. T. Tran, C. Zapata, V. Novak, and M. S. Olufsen, Sensitivity analysis and model assessment: mathematical models for arterial blood flow and blood pressure, J Cardiovasc Eng, 8 (2008), 94–108. [15] L. M. Ellwein, Cardiovascular and respiratory regulation, modeling and parameter estimation, Ph.D. Thesis, Applied Mathematics, North Carolina State University, Raleigh, NC, 2008. [16] G. E. Foster, M. J. Poulin, and P. J. Hanly, Intermittent hypoxia and vascular function: implications for obstructive sleep apnea, Exp Physiol, 92 (2007), 51–65. [17] S. Javaheri, T. J. Parker, J. D. Liming, W. S. Corbett, H. Nishiyama, L. Wexler, G. A. Roselle, Sleep apnea in 81 ambulatory male patients with stable heart failure. Types and their prevalences, consequences, and presentations, Circulation, 97 (1998), 2154–2159, 1998. [18] S. I. Fox, Human Physiology, Wm C Brown Publishers, Dubuque, IA, 1993.

24

L.M. ELLWEIN, S.R. POPE, A. XIE, J. J. BATZEL, C.T. KELLEY, M.S. OLUFSEN

[19] G. S. Francis. Pathophysiology of chronic heart failure, Am J Med, 110 (2001), 37S–46S. [20] N. Fukuda, H. Shinohara, K. Sakabe, S. Morishita, Y. Fukuda, and Y. Tamura, Evaluation of right ventricular function in patients with congestive left heart failure by the Doppler derived total cardiac performance index (tei index), J Echocard, 5 (2007), 11–17. [21] A. Garcia-Touchard, V. K. Somers, L. J. Olson, and S. M. Caples, Central sleep apnea: implications for congestive heart failure, Chest, 133 (2008), 1495–1504. [22] D. Georgiadis, M. Sievert, S. Cencetti Cerebrovascular reactivity is impaired in patients with cardiac failure, Eur Heart J, 21 (2000), 407–413. [23] J. M. Golbin, V. K. Somers, and S. M. Caples, Obstructive Sleep Apnea, cardiovascular disease, and pulmonary hypertension, Proc Am Thorac Soc, 5 (2008), 200–206. [24] G. H. Golub, V. Klema, and G. W. Stewart, Rank Degeneracy and the Least Squares Problem, Technical Report TR-456, Department of Computer Science, University of Maryland, College Park, MD, 1976. [25] G. Golub and C. Van Loan, Matrix Computations, The Johns Hopkins Univ Press, Baltimore, MD, 1989. [26] R. Gopinathannair, A. Mazur, and B. Olshansky, Syncope in congestive heart failure, J Cardiol, 15 (2008), 303–312. [27] F. S. Grodins, J. Buell, and A. J. Bart, Mathematical analysis and digital simulation of the respiratory control system, J Appl Physiol, 22 (1967), 260–276. [28] N. Gruhn, F. S. Larsen, S. Boesgaard, G. M. Knudsen, S. A. Mortensen, G. Thomsen, J. Aldershvile. Cerebral blood flow in patients with chronic heart failure before and after heart transplantation, Stroke, 32 (2001), 2530– 2533. [29] A. Guyton and J. Hall, Textbook of Medical Physiology, W. B. Saunders, Philadelphia, PA, 1996. [30] T. Heldt, Computational Models of Cardiovascular Response to Orthostatic Stress, Ph.D. Thesis, Harvard/MIT Division of Health Science and Technology, MIT, 2004. [31] A. M. Katz, Heart failure: a hemodynamic disorder complicated by maladaptive proliferative responses, J Cell Mol Med, 7 (2003), 1–10. [32] C. T. Kelley, Iterative methods for optimization, SIAM, Philadelphia, PA, 1999. [33] M. C. K. Khoo, A model-based evaluation of the single breath CO2 test, J Appl Physiol, 68 (1990), 393–399. [34] M. C. K. Khoo, A. Gottschalk, and A. I. Pack, Sleep-induced periodic breathing and apnea: A theoretical study, J Appl Physiol, 70 (1991), 2014–2024. [35] C. W. Lee, J. H. Lee, J. J. Kim, S. W. Park, M. K. Hong, S. T. Kim, T. H. Lim, and S. J. Park, Cerebral metabolic abnormalities in congestive heart failure detected by proton magnetic resonance spectroscopy, J Am Coll Cardiol, 33 (1999), 1196 –1202. [36] C. W. Lee, J. H. Lee, T. H. Lim, H. S. Yang, M. K. Hong, J. K. Song, S. W. Park, S. J. Park, J. J. Kim, Prognostic significance of cerebral metabolic abnormalities in patients with congestive heart failure, Circulation, 103 (2001), 2784–2787. [37] G. S. Longobardo, N. S. Cherniack, and B. Gothe, Factors affecting respiratory system stability, Ann Biomed Eng, 17 (1989), 377–396. [38] K. Lu, J. Clark, F. Ghorbel, D. Ware, and A. Bidani, A human cardiopulmonary system model applied to the analysis of the Valsalva maneuver, Am J Physiol, 281 (2001), H2661–H2679. [39] K. Lu, W. Clark, F. H. Ghorbel, C. S. Robertson, D. L. Ware, J. B. Zwischenberger, and A. Bidani, Cerebral autoregulation and gas exchange studied using a human cardiopulmonary model, Am J Physiol Heart Circ Physiol, 286 (2004), 584–601. [40] S. Magder, S. Veerassamy, and J. H. Bates. A further analysis of why pulmonary venous pressure rises after the onset of lv dysfunction, J Appl Physiol, online first Oct. 9, 2008. [41] S. Mak, E. R. Azevedo, P. P. Liu, and G. E. Newton. Effect of hyperoxia on left ventricular function and filling pressures in patients with and without congestive heart failure, Chest, 120 (2001), 467–473. [42] R. Mosteller, Simplified calculation of body surface area, N Engl J Med, 317 (1987), 1098. [43] J. O. Mudd and D. A. Kass. Tackling heart failure in the twenty-first century, Nature, 451 (2008), 919–928. [44] K. M. O’Dell, J. S. Kalus, S. Kucukarslan, and B. Czerska. Nesiritide for secondary pulmonary hypertension in patients with end-stage heart failure, Am J Health Syst Pharm, 62 (2005), 606–609. [45] M. S. Olufsen, J. T. Ottesen, H. T. Tran, L. M. Ellwein, L. A. Lipsitz and V. Novak. Blood pressure and blood flow variation during postural change from sitting to standing: model development and validation, J Appl Physiol, 99 (2005), 1523–1537. [46] S. M. Otis and E. B. Ringelstein, The transcranial Doppler examination: principles and applications of transcranial Doppler sonography, in Neurosonology, by C. H. Tegeler, V. L. Babikian, and C. R. Gomez (eds), Mosby, St. Louis, MO, [47] J. T. Ottesen and M. Danielsen, Modeling ventricular contraction with heart rate changes, J Theo Biol, 22 (2003), 337–346.

CARDIOVASCULAR & RESPIRATORY DYNAMICS IN CONGESTIVE HEART FAILURE

25

[48] J. T. Ottesen, M. S. Olufsen and J. Larsen, Applied Mathematical Models in Human Physiology, SIAM, Philadelphia, PA, 2004. [49] W. W. Parmley. Pathophysiology of congestive heart failure, Clin Cardiol, 15 (Suppl) (1992), I5–I12. [50] S. R. Pope, L. M. Ellwein, C. L. Zapata, V. Novak, C. T. Kelley, M. S. Olufsen, Estimation and identification of parameters in a lumped cerebrovascular model, Math Biosci Eng, 6 (2009). [51] V. Rideout Mathematical and computer modeling of physiological systems, Prentice Hall, Englewood Cliffs, NJ, 1991. [52] G. S. Schroff, J. S. Janicki, and K. T. Weber, Left ventricular systolic dynamics in terms of its chamber mechanical properties, Am J Physiol Heart Circ Physiol, 14 (1983), H110–H124. [53] D. D. Sin, F. Fitzgerald, J. D. Parker, G. Newton, J.S. Floras, T.D Bradley, Risk factors for central and obstructive sleep apnea in 450 men and women with congestive heart failure, Am J Respir Crit Care Med, 160 (1999), 1077-1078. [54] P. Solin, T. Roebuck, D. P.Johns, E.H. Walters, M.T. Naughton, Peripheral and central ventilatory responses in central sleep apnea with and without congestive heart failure, Am J Respir Crit Care Med, 162 (2000), ; 162(6):2194-2200. [55] N. Stergiopoulos, J. Meister, and N. Westerhof, Determinants of stroke volume and systolic and diastolic aortic pressure, Am J Physiol, 270 (1996), H2050–H2059. [56] J. Stocks and P. H. Quanjer, Reference values for residual volume, functional residual capacity and total lung capacity, J Biomech Eng, 124 (2002), 9–20. [57] S. Y. Sun, W. Wang, I. H. Zucker, H. D. Schultz, Enhanced peripheral chemoreflex function in conscious rabbits with pacing-induced heart failure, J Appl Physiol, 86 (1999), 1264–1272. [58] V. Thomas, F. Costes, T. Busso, Estimation of arterial PCO2 from a lung model during ramp exercise in healthy young subjects, Resp Physiol Neurobiol, 156 (2007), 259–265. [59] R. Tkacova, M. J. Hall, P. P. Liu, F. S. Fitzgerald, and T. D. Bradley, Left ventricular volume in patients with heart failure and cheyne-stokes respiration during sleep, Am J Respir Crit Care Med, 156 (1997), 1549–1555. [60] M. Ursino Interaction between carotid baroregulation and the pulsating heart: a mathematical model, Am J Physiol, 44 (1998), H1733–H1747. [61] G. B. West, J. H. Brown, and B. J. Enquist, A general model for the origin of allometric scaling laws in biology, Science, 276 (1997), 122–126. [62] C. R. White and R. S. Seymour, Allometric scaling of mammalian metabolism, J Exp Biol, 208 (2008), 1611–1619. [63] A. Xie, J. B. Skatrud, R. Khayat, J. Dempsey, B. Morgan, and D. Russell, Cerebrovascular response to carbon dioxide in patients with congestive heart failure, Am J Resp Crit Care Med, 172 (2005), 371–378. [64] A. Xie, J. V. Skatrud, B. Morgan, B. Chenuel, R. Khayat, K. Reichmuth, J. Lin, J. A. Dempsey, Influence of cerebrovascular function on the hypercapnic ventilatory response in healthy humans, J Physiol, 577 (Pt 1) (2006), 319–29.