Study of cardiovascular function using a coupled left ventricle and systemic circulation model

To gain insight into cardio-arterial interactions, a coupled left ventricle-systemic artery (LV–SA) model is developed that incorporates a three-dimensional finite-strain left ventricle (LV), and a physiologically-based one-dimensional model for the systemic arteries (SA). The coupling of the LV model and the SA model is achieved by matching the pressure and the flow rate at the aortic root, i.e. the SA model feeds back the pressure as a boundary condition to the LV model, and the aortic flow rate from the LV model is used as the input for the SA model. The governing equations of the coupled system are solved using a combined immersed-boundary finite-element (IB/FE) method and a Lax–Wendroff scheme. A baseline case using physiological measurements of healthy subjects, and four exemplar cases based on different physiological and pathological scenarios are studied using the LV–SA model. The results of the baseline case agree well with published experimental data. The four exemplar cases predict varied pathological responses of the cardiovascular system, which are also supported by clinical observations. The new model can be used to gain insight into cardio-arterial interactions across a range of clinical applications.

As changes in arterial properties can alter the heart function and vice versa (Noguchi et al., 2011), in this paper, we focus on the coupling of the heart and the arteries, by combining a models of a 3D left ventricle (LV) with a systemic arteries (SA) model that uses a structural tree description of the vascular beds containing the smaller arteries. The 3D LV is reconstructed from a dataset of in vivo magnetic resonance imaging (MRI) of a healthy volunteer (Gao et al., 2015a), which includes fluid-structure interaction (FSI). The SA model is based on the development by Olufsen (1999) and Olufsen et al. (2000), which includes both large arteries and remote vascular beds.

Methodology
The coupled left ventricle and the systemic artery (LV-SA) model is shown in Fig. 1. The methodologies for the 3D LV and the 1D SA models have been published elsewhere (Gao et al., 2014a;Olufsen et al., 2000Olufsen et al., , 2012, but are briefly described here to explain the coupling procedure.

The LV model
The LV model consists of the valvular and inflow/outflow tracts (assumed passive), and the active LV region. The model is solved using a combined immersed boundary finite element (IB/FE) method. Let Ω & R 3 denote the physical domain occupied by the fluid-structure interaction (FSI) system, in which x ¼ ðx 1 ; x 2 ; x 3 Þ A Ω are fixed Eulerian coordinates. Let U & R 3 denote the reference configuration of the immersed solid, in which X ¼ ðX 1 ; X 2 ; X 3 Þ A U are Lagrangian coordinates. χ ðX; tÞ describes the physical position of in which μ is the fluid viscosity, uðx; tÞ is the fluid velocity of the blood, and pðx; tÞ is the pressure, NðXÞ is the exterior unit normal to U, and δðxÞ is the three-dimensional Dirac delta function. We assume that the fluid and the solid have the same density ρ. P s ðX; tÞ ¼ detðFÞσ s F À T is the first Piola-Kirchhoff stress tensor, and σ s is the structure Cauchy stress tensor, where σ a ¼ T 0 Tðf fÞ is the active stress, while the contractile tension T is described by the myofilament model of Niederer et al. (2006), which is triggered by a prescribed intracellular calcium transit (Hunter et al., 1998), as shown in Fig. 3. T 0 is introduced to make the contraction patient-specific. The passive stress σ p is determined through the Holzapfel-Ogden strain energy function (Holzapfel and Ogden, 2009), as detailed in the Appendix A.
Equations (1) and (2) are discretized using a finite-difference method, and Eqs. (3) and (4) are discretized with a finite-element method. The material parameters in (12) are determined inversely by fitting the measured end-diastolic volume and myocardial strains using a multi-step optimization procedure (Gao et al., 2015b). T 0 is determined by matching the measured stroke volume (Gao et al., 2014a). The valvular region is modelled as a neo-Hookean material, with the shear modulus adjusted so that its maximum displacement agrees with MRI measurements. The inflow/outflow tracts are both assumed to be rigid, with the inlet and outlet annuli fixed in space.

The systemic arterial model
The SA model consists of 24 large arteries modelled as a onedimensional cross-sectional-area-averaged flow and pressure. Each terminal vessel in the network of the large arteries is coupled with a group of small arteries (the vascular bed), which are modelled as an asymmetric structured-tree to provide outflow boundary conditions (Olufsen, 1999;Olufsen et al., 2000Olufsen et al., , 2012.
The governing equations for the SA model are Pðx; tÞÀP 0 ¼ 4 3 where Q is the volumetric flow rate, A is the cross-sectional area, P is the averaged cross-sectional pressure, P 0 is the constant external pressure, ρ is the density, ν is the kinematic viscosity, R is the radius of the vessel. A 0 and r 0 are the cross sectional area when P ¼ P 0 , δ n is the width of the boundary layer (δ n ⪡R), h is the wall thickness, and E is the arterial Young's modulus, computed as, where k i (i¼ 1-3) are material constants (Olufsen et al., 2000). These equations are solved numerically using the Lax-Wendroff scheme (Lax and Wendroff, 1960). With each vascular bed, each parent artery with radius r p divides into two daughter arteries with radius r d 1 ¼ αr p and r d 2 ¼ βr p , (0oβ oα o 1). The bifurcation process continues until the radius of the daughter vessel reaches the minimum radius r min ¼ 100 μm. The radius exponent ξ, the asymmetry ratio η, and the area ratio γ (only two of which are independent) are given by

Coupling of the SA model and the LV model
The coupling is achieved by matching the pressure (P À ) and the flow rate (Q À ) at the outflow tract (corresponding to the beginning of the ascending aorta) of the 3D LV model to these of 1D SA model (P þ and Q þ ) at interface plane Γ a (Fig. 2). Subscripts ' þ' and ' À ' are used for representing variables in the SA and the LV models, respectively. Before coupling, the SA model is initialized for four periods using a prescribed cardiac output from a decoupled LV model (Gao et al., 2015a). The parameters of the baseline case are summarized in Appendix B.
Since the LV and SA models are coupled only in systole when the aortic valve (AV) is open, a set of simplified boundary conditions are used in our coupled model as detailed in Appendix C. In the SA model, the time step is Δt SA ¼ t period =N ¼ 0:9=8192 % 1:10 Â 10 À 4 s, in which t period is the length of period and N is the number of time steps during one period. In the LV model, a basic time step Δt LV ¼ 1:22 Â 10 À 4 s is used in the diastolic and relaxation phases, but smaller time steps are used in systole (0:125Δt LV ) as much higher structural stress is generated from the myocardial contraction. Interpolation is used to ensure the pressure and the flow rate are exchanged at the right time points when the two models are coupled. Simulation of the 3D LV model is implemented using the open-source IBAMR framework (https:// github.com/IBAMR/IBAMR); the 1D SA model and the interactions are implemented using C þ þ within the same framework. The LV-SA model converges to periodic solutions in the second period after coupling, and the maximum pressure difference between the 2nd and the 3rd periods is less than 1%. Hence the results from the second period are presented. The simulation time for a cardiac cycle is around 168 h (7 days) on a local Linux workstation with eight Intel(R) Xeon(R) CPU cores (2.65 GHz) and 32 GB RAM.

Cases different to the baseline
We also consider the following pathological scenarios.
Case 1 -Stiffening of the arterial wall: Increased arterial stiffness is associated with an increased risk of cardiovascular events such as myocardial infarction and stroke (Lee and Oh, 2010;Avolio et al., 1983). Here we increase k 3 in (9) by 100% in the large arteries to investigate how arterial stiffness affects the LV function.
Case 2 -Vascular bed rarefaction: Rarefaction, a reduction in the density of small arterioles, is often associated with Type II diabetes (Hopkins and McLoughlin, 2002), and has a significant impact on cardiovascular system (Olufsen et al., 2012;Safar et al., 2003). Rarefaction is simulated by decreasing the radius exponent ξ in (10) from 2.76 to 2.4, thus reducing the total number of small vessels by 71% compared to the baseline case.
Case 3 -Decreased compliance: LV hypertrophy or fibrosis is usually linked with a stiffer myocardium (Beckett et al., 2008). We increase the myocardium stiffness by doubling the values of a, a f , a s and a fs in Eq. (12) in Appendix A.
Case 4 -Increasing the myocardial contractility: Myocardial contractility, or isotropy, is one of the three primary mechanisms that regulate the LV stroke volume. To model the effects of increased LV contractility, we increase T 0 in (5) by 33% from the baseline case.
Cases 1_iso and 2_iso: We also apply Cases 1 and 2 to an isolated SA model as Case 1_iso and Case 2_iso, respectively, to identify coupling effects.

The baseline case
The results of the coupled model in systole for the baseline case are plotted in Fig. 3; the intracellular calcium transient taken from Hunter et al. (1998) is also shown for comparison. The corresponding clinical indices are summarized in Table 1. The ejection fraction (EF) is 51%. The pressure curve closely follows the active stress, with a small delay behind the intracellular calcium transient.
The deformed LV geometry (in red) at end-diastole and endsystole are shown in Fig. 4(a-d), superimposed on the corresponding cine MR images. Qualitative assessment by visual inspection suggests that the computed LV geometries agree well   (1), peak value ¼10 À 7 μmol. All results are normalized with their respective peak values, except for the pressure at outflow tract, which is normalized with the peak LV center cavity pressure. The average active stress (2) has peak value ¼ 68.9 kPa, the LV center cavity pressure (3) has peak value ¼113.1 mmHg, the aortic flow rate (4) has peak value ¼532 mL s À 1 , and the pressure at the outflow tract (5) has peak value¼ 111.7 mmHg. Time t¼ 0 indicates the beginning of the isovolumetric contraction.
with the MR images in most parts of the LV, with exception of the apex, where our model overshoots slightly. This is presumably due to the fact that our model does not include the pericardium, which would constrain the apical motion. The myofibre strain field at end-systole is shown in Fig. 5(a). The mean value of the strain is À 0.17, which is comparable to previously reported systolic strain values (Moore et al., 2000). Figure 5(b) shows the active stress (σ a ) during mid-systole. All regions appear to contract equally hard with the exception of the apical region. Flow patterns at mid-diastole and mid-systole are shown in Fig. 5(c) and (d) respectively. Since there is no mitral valve, the flow vortex in Fig. 5(c) is not as distinct as when the mitral valve is present (Yin et al., 2010). This inaccuracy in the flow fields may not significantly affect the overall results as it is commonly agreed that blood flow affects the LV deformation mostly through the pressure. The pressure and flow waveforms in selected arteries are shown in Fig. 6, with the indices summarized in Table 2. A clear trend in the time delay of the peak pressure during systole with distance away from the heart can be seen from Fig. 6(a). At 0.25 s, there is a pronounced dicrotic notch in the ascending aorta due to the reflected flow at end-systole, which fades away in the distal arteries. A similar time delay is present in the flow rate waveform in the distal arteries ( Fig. 6(b)); the peak flow rate arrives approximately 64 ms earlier in the ascending aorta than that in the abdominal artery. Figure 6(c) and (d) plots the pressure and flow rate waveforms in the long arteries such as the carotid, brachial and femoral arteries. The systolic (126 mmHg) and diastolic Table 1 Comparison of systolic LV pump functions of the baseline and Cases 1-4. All changes for Cases 1-4 are in percentage compared to the baseline case.

Indices
T peak P peak À ao P peak À LV Q peak À ao SV SW t ej  (81 mmHg) pressures in the mid-brachial artery are within the normal range (Sesso et al., 2000). In summary, the results from both the LV and the systemic circulation in the baseline case lie within a normal physiological range (Starling, 1993;Moore et al., 2000;Sesso et al., 2000).

Cases 1-4
The results of the active stress, LV pressure and the aortic flow rates from the coupled model in systole are compared for all the cases in Fig. 7. We can see that stiffened arteries (Case 1) result in a higher ventricular pressure, and a reduced aortic flow rate. Rarefaction (Case 2) causes a similar flow rate reduction as in Case 1, but has the highest increase in pressure and active stress. Decreased LV compliance (Case 3) leads to a much smaller enddiastolic volume. As the active stress is a function of strain, a smaller end-diastolic volume indicates a small expansion of the LV, hence both the active stress and LV pressure decrease, leading to a lower aortic flow rate and shorter systolic ejection duration. Increased contractility (Case 4) leads to a moderate pressure rise compared to the baseline case, but causes the highest flow rate. Interestingly, Case 4 does not yield the highest active stress ( Table 1).
The pressure-volume loops are shown in Fig. 7(d). Cases 1-3 all lead to a reduced stroke volume, with Case 3 (stiffener LV) being the most severe. The decreased compliance (Case 3) has the most notable reduction in the stroke work (indicated by the enclosed area of the pressure-volume loop), while increased contractility (Case 4) is associated with the largest stroke work. The indices of the LV pump function are summarized in Table 1. Case 3 (stiffer LV) has the worst pump function indices, with stroke volume reduced by 31%, and significant drops in the LV pressure and active tension. Compared to the baseline case, there is an increase in the peak active tension and LV pressure, and a smaller decrease in the stroke volume in Cases 1 and 2, although the stroke work is higher in Case 2, and lower in Case 1. Both the stroke volume and the stroke work are significantly increased in Case 4 due to the faster ejection speed.
The changes in the LV are carried through to the systemic arteries, as shown in Fig. 8. The key indices of the SA system are summarized in Tables 2 and 3 for selected arteries.

Case 1_iso and Case 2_iso
To identify LV effects on the systemic circulation in pathological situations, we also compare the results from the isolated systemic circulation model with these of the coupled model for Cases 1 and 2. As expected, the general trends in the pressure and flow rate waveforms are similar in both cases. However the isolated SA model tends to over-estimate the results by up to 7% for pressure, and up to 20% for flow rate, as shown in Table 4.

Discussion
Our coupled model makes it possible to study more detailed interactions between the LV and the systemic arteries under various normal and pathological conditions during systole. Although the LV and the systemic models were initially derived from different healthy human subjects, through inverse parameter estimation, we are able to match the measurements of the LV subject according to the end-diastolic and stroke volumes for the baseline case. All of the computed indices are in line with the previously reported values (Sesso et al., 2000;Moore et al., 2000;Mahadevan et al., 2008).
The coupled model is used to examine four different scenarios. Cases 1 and 2 show how the LV reacts to an increased LV afterload due to changes in the systemic circulation. Cases 3 and 4 show the impact of the changed LV function on the systemic circulation. For both Cases 1 and 2, the LV needs to generate a higher active tension and pressure, yet still suffers from a reduced stroke volume. This could promote adverse remodelling of the LV in the longer term. Rarefaction also increases both the peak and mean pressures in the large arteries, which agrees with clinical observations that rarefaction in vascular beds can lead to hypertension (Noon et al., 1997;He et al., 2010;Antonios, 2006). In Case 3, the stiffer LV represents a potential heart disease, post-myocardial fibrosis. This scenario gives the most markedly reduced pumping performance: the stroke volume, active tension, and LV pressure are all significantly decreased. With increased T 0 (Case 4), an immediate

Table 2
Summary of the indices of the pressure waveform: the peak arterial blood pressure (P peak ), the time when pressure peaks (t peak ), and the arterial trough pressure (P trough ), recorded at mid-ascending aorta, mid-femoral and mid-brachial arteries.

Position
Mid-ascending aorta Mid-femoral Mid-brachial Indices t peak P peak P trough t peak P peak P trough t þ2.9 þ 15.9 þ 34.6 þ 0.5 þ19.3 þ 37.5 þ 2.9 þ 11.9 þ 34.6 Case 3 (%) À 4.7 À 11.5 À 8.6 À 1.9 À 14.3 À 8.8 effect is the increased stroke volume (by 13%). This suggests that during certain conditions, myocardium could remodel itself to enhance contractility in order to meet the flow demand (Zhang et al., 2010). However increased myocardial stress in this case may further induce myocyte hypertrophy and a consequent increase in myocyte death rate (Gomez et al., 2001;Abbate et al., 2003). Cuff (brachial) pressure is routinely measured by clinicians as an indication of the central (ascending aorta) pressure and the pressure inside LV. However, the difference between the cuff and the central pressures in vivo is unclear. Our model shows that the pressure drop between the LV and the ascending aorta is minor, see Fig. 9(a), and the pressure difference between the brachial artery and the center of the LV for the baseline case is about 12.41 mmHg. This agrees with Kroeker and Wood (1955) who observed that, although mean and diastolic blood pressure are relatively constant throughout the arterial tree, there is a gradual increase in systolic pressure moving from the aorta to the peripheral arteries. Indeed, the pressure measured at the brachial artery can be 5-20 mmHg higher than that in the ascending aorta. This pressure amplification arises principally because of increased vessel stiffness and changes in vessel geometries. Our results also suggest that this pressure difference varies, as shown in Fig. 9(b). Stiffer arteries (Case 1) and increased contractility (Case 4) cause the greater pressure amplifications, while the difference is less in rarefaction (Case 2). The most prominant pressure amplification is seen in Case 3 (stiffer LV).
We can show how the arterial elastance (E A , the ratio of endsystolic LV pressure to stroke volume) changes the shape of the P-V loop and the end-of-systolic P-V relationship (ESPVR) in Fig. 10. E A can be increased either by stiffer arteries (Case 1) or rarefaction (Case 2). The interaction between the LV and the systemic circulation can also be assessed by the LV chamber elastance (E S ), the slope of ESPVR (Sunagawa et al., 1983;Chirinos, 2013), and by the ratio E A =E S . Clinically, E S can be obtained by producing step-wise pharmacological afterload variations without inducing inotropic changes. In the simulated Cases 1 and 2, the afterload for the coupled LV model (aortic pressure) is increased, but the contractility is not, therefore we can obtain E S by fitting the three endsystolic points (Cases 1, 2 and baseline) at the P-V curves (Fig. 10). The ratio E A =E S in the baseline case is 0.72, which lies in the reported range for normal humans (0.62-0.82) (Redfield et al., 2005). In fact, it is close to the value when the myocardial energetic efficiency is maximized (0.7) (De Tombe et al., 1993). E A =E S is 0.87 and 0.95 for Cases 1 and 2, suggesting that the coupling between the LV and systemic circulation is sub-optimal. Our model predictions are consistent with clinical interpretations of how E A and E s change the P-V curve, but provide more quantitative information on exactly what causes such changes. Such a coupled model could be used together with measurable values (such as cuff pressure), to give greater clinical insight.
Comparison between the results of the coupled model and that of the isolated circulation model suggests that the isolated circulation model overestimates the peak pressure and the flow rate. This is an example of how the circulation system is affected by the LV. When either increase in the arterial stiffness or rarefaction, the LV tends to increase the active tension and decrease the cardiac output, which in turn affects the systemic circulation.
We now mention the limitations of the model. The LV and systemic arteries are from two different volunteers. To make the coupled model work physiologically and agree with the measured systolic and diastolic volumes of the LV subject, a number of parameters are tuned, such as the external pressure p 0 in the SA model, the myocardial passive stiffness, and T 0 , which controls the magnitude of the active tension. All other parameters in the SA model are kept the same as in Olufsen et al. (2000). Therefore the coupled model is not patient-specific per se. In our model of the systemic arteries, the pulse wave velocity (PWV) is determined by the pressure-area relationship for the arteries (the tube law), specifically cðPÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi ðA=ρÞð∂P=∂AÞ p . We have studied the clinical estimates of the PWV in a model of the pulmonary circulation (Qureshi and Hill, 2015). It would be possible in principle to use the clinical estimates of the PWV to assess the tube law although errors in the measurements may make this difficult. However, it is probably even more interesting to use wave intensity analysis to identify waves reflected from bifurcations in the arterial system and to study their contribution to the pressure loading on the LV during systole in health and disease. This work is underway.
Throughout the study, we have kept the end-diastole pressure at 8 mmHg, although this value is known to change in pathological conditions. In addition, we use velocity boundary conditions to mimic the action of the AV, thus the pressure drop across the AV is Table 3 Summary of indices of the flow waveform: the peak arterial flow (Q peak ), the time when flow peaks (t peak ), and the maximum backflow (Q max À back ), recorded at midascending aorta, mid-femoral and mid-brachial arteries.

Position
Mid-aorta Mid-femoral Mid-brachial Indices t peak Q peak Q max À back t peak Q peak Q max À back t peak Q peak Q max À back Unit s mL s À 1 mL s À 1 s m L s À 1 mL s À 1 s m L s À 1 mL s À 1  not accounted for. The LV and the systemic arteries coupling only occurs during systole (from end-diastole to end-systole) because when the AV is closed, the flow in the arterial system is separated from the LV. Moreover, the LV diastolic filling is simplified by not including the mitral valve, the left atrium and the pulmonary circulation. While our coupled model can address a range of clinical questions as we have illustrated here, more involved pathological problems may require the development of more sophisticated models of cardiovascular system that include the four chambers of the heart, as well as the pulmonary circulation.

Conclusions
We have developed a coupled model for studying the interaction between a three-dimensional contracting left ventricle model and a structured-tree based cross-sectional-area-averaged systematic circulation model. The coupled model can predict details of the left ventricular dynamics, pressure and flow rate profiles at any position in the systemic arteries throughout the cardiac cycles, and thus it provides a powerful in silico tool for exploring and understanding cardio-arterial interactions. This new model is used to study a number of pathological changes in the left ventricle and the systemic circulation. The results from the coupled model are consistent with clinical observations. We find that stiffening of the arterial wall and functional rarefaction in the remote vascular beds cause higher blood pressure along with higher LV active tension, but with reduced stroke volume. A stiffer LV leads to severely impaired pump function with low active tension, stroke volume and low blood pressure. Increased contractility can help the heart to maintain a higher stroke volume, but gives rise to an elevated pressure in the circulation. Furthermore, the model can be combined with clinical measurements, such as cuff pressure, to infer pressure profiles inside the LV. With further development towards patient-specific individualization, the model can be applied to a range of clinical studies for exploring the causes and development of cardiovascular diseases and their potential treatments.  9. (a) The predicted left ventricular/aortic pressures vs cuff pressure for the baseline case, and (b) comparison of the pressure difference from the LV to station "X" at the brachial artery as marked in Fig. 1, for all the simulated cases. Fig. 10. Derivation of the effective arterial elastance (E A ) and LV end-systolic elastance (E S , the slope of the ESPVR or end-systolic pressure-volume relation) based on the P-V loops for the Baseline case, and Cases 1 and 2.

Table 4
Comparison of the peak pressure (P peak ) and peak flow rate (Q peak ) between the coupled model and the isolated SA model for Cases 1 and 2, estimated at the midascending aorta, mid-femoral and mid-brachial arteries.