|Home | About | Journals | Submit | Contact Us | Français|
In this study we present a novel, robust method to couple finite element (FE) models of cardiac mechanics to systems models of the circulation (CIRC), independent of cardiac phase. For each time step through a cardiac cycle, left and right ventricular pressures were calculated using ventricular compliances from the FE and CIRC models. These pressures served as boundary conditions in the FE and CIRC models. In succeeding steps, pressures were updated to minimize cavity volume error (FE minus CIRC volume) using Newton iterations. Coupling was achieved when a predefined criterion for the volume error was satisfied. Initial conditions for the multi-scale model were obtained by replacing the FE model with a varying elastance model, which takes into account direct ventricular interactions. Applying the coupling, a novel multi-scale model of the canine cardiovascular system was developed. Global hemodynamics and regional mechanics were calculated for multiple beats in two separate simulations with a left ventricular ischemic region and pulmonary artery constriction, respectively. After the interventions, global hemodynamics changed due to direct and indirect ventricular interactions, in agreement with previously published experimental results. The coupling method allows for simulations of multiple cardiac cycles for normal and pathophysiology, encompassing levels from cell to system.
Three-dimensional (3D) finite element (FE) mechanics models of the whole heart have been used extensively to investigate regional myofiber strains and stresses during normal and abnormal function, including ischemia, ventricular pacing, myofiber disarray, and heart failure (for an overview see Ref. 21). To ensure reasonably realistic pressure waveforms and volume changes in the ventricles, FE heart models of the ventricles have been coupled to simple afterload models, i.e. two-40 or three-element Windkessel models,2,47 while filling of the ventricles2,20,40 has been typically prescribed by an unphysiologic (i.e. linear) increase in diastolic pressure. Furthermore, the timing of mitral valve closure is usually a parameter input to the models, not determined by the physics of the problem. To apply a more realistic preload, Watanabe et al.47 have coupled a model of pulmonary venous return to a left ventricular FE model, in which the left atrium was modeled as a time-varying elastance. In order to prescribe realistic boundary conditions, 3D fluid mechanics models of large arteries have been coupled to zero- or one-dimensional models.8,11,44
Three major improvements are needed in studies in which 3D FE ventricular solid mechanics models are coupled to preloads and afterloads: (1) A more realistic filling phase is needed in relaxing ventricles. In the normal heart, the mitral and tricuspid valves open while the ventricles are still generating some active force, an important aspect of diastolic heart failure. This has not been taken into account yet in existing models. (2) Conservation of blood mass: A higher cardiac output results in increased venous return. In other words, FE models of the heart need to be embedded in a closed circulation. (3) Steady-state solutions need to be assessed, independent of transient effects due to varying initial conditions, as well as transient alterations due to interventions (e.g. a change in impedance or contractility). This means that simulations of multiple cardiac cycles are needed. Furthermore, a multi-scale cardiovascular model, in which a 3D FE ventricular model is coupled to a model of the circulation, will allow the investigation of phenomena such as direct and indirect ventricular interactions,48 and the dependence of global hemodynamics on changes at the cellular level.
In this study we present a novel method to couple FE models of cardiac mechanics to closed-loop lumped models of the circulation. The method was set up in a modular manner that is independent of the FE solution algorithm, and allows coupling to different circulatory models of arbitrary complexity. Special attention was also paid to prescribing initial conditions. To reduce computation time, initial conditions were obtained by coupling the circulatory model to a simple time-varying elastance model of the heart which takes into account ventricular interactions. Steady-state values of the state variables from this simpler coupled model were then used as initial conditions for the multi-scale coupled model. To investigate the performance of the multi-scale cardiovascular model, we performed two interventions: one by applying acute ischemia in the left ventricle; and one by applying pulmonary artery constriction (PAC).
The complete cardiovascular model, which consists of a 3D finite element ventricular model coupled to a lumped parameter systems model of the circulation, is illustrated in Fig. 1a. Synchronous activation of the myofibers in the ventricular walls results in myofiber contraction, generating active forces, thus increasing the pressures in the ventricular cavities. After the aortic and pulmonary valves open, the finite element model drives the blood – characterized by pressures, volumes, and flows – through the systemic and pulmonic circulations. Blood returning from the systemic and pulmonic circulations fills the respective ventricles for the next cycle.
In order to determine the initial conditions for the working model, a time-varying elastance model of the ventricles that includes direct ventricular interactions was developed. The latter model was coupled to the same circulation model. Computation times of a varying elastance model are much shorter than those of a finite element model. Hence, the varying elastance model was used to obtain parameters for the circulation, and steady-state variable values from this simplified model served as initial conditions for the finite element model coupled to the circulation.
Continuity of flow was imposed to couple the finite element model of the heart to the circulation properly. LV and RV compliances were computed at each time step to calculate cavity pressures. These pressures served as hemodynamic boundary conditions in the finite element and circulation model. This calculation also yielded co-compliances, which characterize direct interactions between the ventricles. The individual components of the model are described below in more detail.
The anatomic model of the canine ventricles has been published32 and used40 before, and includes 3D left and right ventricular geometry and a 3D myofiber angle distribution. The unloaded resting LV and RV cavity volumes were 26.1 and 22.3 ml, respectively. Total wall volume was 131.8 ml. Passive stress was modeled using a strain energy function for an exponential transversely isotropic material. Parameter values for this model were determined14 and validated41 previously for canine myocardium.
Active stress was modeled using a time-varying material law46 in the myofiber direction. Active stress is also generated transverse to the myofibers, and was about 40% of active stress generated in the myofiber direction.42 Furthermore, an equation was added to make time to peak stress dependent on sarcomere length. Thus for each ventricle, larger cavity volumes result in a longer times to peak pressure.4 The parameters that control the length of a twitch were adjusted (m = 131.2 ms/µm and b = −94.34 ms46), such that isovolumic pressure curves matched those of dogs measured by Burkhoff et al.4
Ischemia was modeled by decreasing myofilament calcium sensitivity transmurally in the apical region (~30% of wall volume) that is normally perfused by the left anterior descending (LAD) coronary artery, as in the study of Mazhari et al.28 (Fig. 1B). This resulted in a decrease in active force in the ischemic region of about 80%.
A time-varying elastance model of the LV and RV was developed in order to simplify the determination of initial conditions for the finite element coupled model. Commonly used time-varying elastance models relate cavity pressure p to cavity volume V through a time-dependent elastance E(t) (the reciprocal of compliance). For a ventricle this relation is
Using both ventricular time-varying compliances and co-compliances, a relation can be written between LV and RV pressures and volumes, which also takes into account ventricular interaction:
where and are column vectors with LV and RV pressures and volumes, respectively. C̲ is the 2 by 2 compliance matrix with compliances and co-compliances on the diagonal and off-diagonal, respectively.
The parameters of the compliance matrix were obtained from the pressure–volume relationship of the finite element model. Figure 2 shows the LV and RV pressure and volume relations for the passive and fully activated canine FE model, decoupled from the circulation. Increasing LV pressure (in phase 1), while keeping RV pressure at zero, not only results in increasing LV volume, but also in a decreasing RV volume, as a result of direct ventricular interaction through the myocardium. Next, in phase 2, LV pressure is kept constant and RV pressure is increased: this also affects LV volume. In phase 3, when RV pressure is kept at a constant level, and LV pressure is decreased, RV volume increases. In the last phase, LV volume increases, while RV pressure is decreased. Thus, the PV-relation changes in one ventricle when the pressure is changed in the other ventricle, which has been demonstrated experimentally during diastole.48 Hence, besides compliances for the LV and the RV (dVL/dpL and dVR/dpR), there exist cross- or co-compliances (dVL/dpR and dVR/dpL). These co-compliances are a quantification of direct ventricular interaction, where a change in pressure in one ventricle affects the volume changes in the other ventricle through the myocardium. For more details see Appendix A.
The systemic and pulmonary circulations (Fig. 1a, Appendix B) were each modeled as two lumped Windkessel compartments in series, one compartment for arterial and capillary blood and one for venous blood. Each compartment was characterized by a resistance and compliance parameter. Given a pressure drop Δp along a circulatory segment, the resistance parameter R determines forward flow Q through the segment using the fluid analog of Ohm’s Law. For the valves, flow is zero when Δp becomes negative. Flow through the valves was modeled in an “all or nothing” fashion (rectifiers in Fig. 1a): they were either completely open or closed (see e.g. Eq. B11 in Appendix B). Compliance parameters C (fixed capacitors in Fig. 1a) determine the pressure–volume (PV) relationship for each segment through the fluid analog of the law of capacitance. All segments were therefore modeled with linear PV relationships. Using conservation of mass (for incompressible blood leading to conservation of volume), the volume change of a Windkessel segment was determined by inflow minus outflow.
The atria were represented with time-varying elastance models27 (Appendix B) in which maximum and minimum elastance parameters set the peak systolic and diastolic stiffness throughout the cardiac cycle, respectively. Unloaded atrial volumes were also dependent on the atrial activation level. The atria contracted 120 ms before ventricular contraction (PQ interval in the ECG).
In order to estimate appropriate initial volumes and circulatory parameter values for the FE-driven circulatory model, a simpler cardiovascular model in the style of Rideout33 was developed. In this simpler model the circulatory components described above were coupled to time-varying elastance models of the right and left ventricles that also take into account ventricular interaction (Appendix A). This closed-loop, varying elastance driven system runs faster than real-time and was tuned to produce physiologically typical canine target values for pressures, flows and volumes throughout the circulation. These target values were determined using literature and global hemodynamic properties of the FE model, as follows.
The target value for cardiac output was computed using a formula that relates body weight (estimated at 20 kg) to cardiac output.9 Heart rate was set at 100 beats per min. Stroke volume was then calculated by dividing cardiac output by heart rate. Target end-diastolic volume for the left ventricle was computed by multiplying the FE model’s diastolic unloaded LV volume by the ratio of reference human end-diastolic LV volume to human unloaded LV volume. The value for the reference human unloaded diastolic volume was determined using a formula from Klotz et al23 and a physiologically normal end-diastolic pressure and volume pair for human (0.7 kPa, 126 ml). End-diastolic volume of the RV was calculated in the same way.
End-diastolic pressures for LV and RV were calculated from the FE model’s passive compliance and end-diastolic volumes. Average right atrial pressure was estimated to be the same as end-diastolic RV pressure. Due to the steeper diastolic PV curve of the LV, average left atrial pressure was estimated to be about 0.2 kPa less than end-diastolic LV pressure, in accordance with the average (not end-) diastolic LV pressure.
End-systolic volume for each ventricle was calculated by subtracting the stroke volume from the end-diastolic volume. Target end-systolic ventricular pressures for maximally activated ventricles were calculated using the FE model’s compliance at maximum activation. By taking the end-systolic left ventricular pressure and then subtracting the product of the cardiac output and the estimated aortic valve resistance, the target pressure at the proximal arterial end of the systemic circulation was determined. This same method was used to compute the target pressure at the proximal end of the pulmonary circulation as well. To compute the target pressure at the proximal end of the systemic venous compartment, a reference literature value was used,9 whereas the target upstream venous pressure in the pulmonary circulation was calculated by adding half the pressure drop across the pulmonary circuit to the target left atrial pressure.
With the distribution of target pressures and the cardiac output defined, resistance parameters for each circulatory segment were calculated by dividing upstream compartment pressure by cardiac output. Aortic and pulmonary artery valve resistances were set a priori at 0.007 and 0.004 kPa.s/ml. Next, the distribution of target blood volumes throughout the circulation was determined. First the total blood volume was calculated from the estimated animal weight using a formula from Altman and Ditmer9 (86.2 ml of blood/kg). Target blood volumes for all circulatory segments other than the ventricles were computed by multiplying the total blood volume by reference human volume fractions.16 Compliance parameters for each circulatory component were then determined by dividing the compartment’s target blood volume by the target upstream pressure.
With the resistance and compliance parameters thus obtained, the varying elastance driven closed-loop model was run, and model outputs were compared with target values. By only adjusting the resistance of the systemic venous compartment, we were able to balance venous return and cardiac ejection such that steady-state model solutions were within 5% of all target pressures, flows and volumes.
This final parameter set and initial conditions, determined from variables in steady state, were then copied into the FE-driven circulatory model for simulation of normal cardiac dynamics. To simulate pulmonary artery constriction, the parameter for pulmonary valve resistance was increased to 0.030 kPa.s/ml. This produced a peak-systolic pressure gradient of about 5.5 kPa between the RV and the pulmonary arterial segment in the FE model coupled to the circulation. A gradient of this magnitude is classified as a “mild stenosis” in the clinical literature.3
The ventricular–vascular coupling method is based on the method of Bovendeerd et al.,2 in which ventricular cavity pressures are estimated and updated until convergence has been achieved. Here, throughout the whole cardiac cycle for each new time step, LV and RV cavity pressures were extrapolated in time using pressure history, just as in Bovendeerd et al.2 But next, in our method, these pressures served as hemodynamic boundary conditions both in the FE model and circulatory models (Fig. 3). When the difference in cavity volumes from the finite element model (VFE) and the circulatory model (Vcirc, see Eq. B16 and B26) was small enough, time was incremented. Otherwise, LV and RV pressures were perturbed in the next 2 updates to obtain ventricular compliances (ΔVLV/ΔpLV, ΔVRV/ΔpRV [ml/kPa]) and co-compliances (ΔVRV/ΔpLV, ΔVLV/ΔpRV [ml/kPa]) from the FE and circulatory models. The volume residual (FE cavity volume minus circulatory cavity volume) was minimized in the remaining updates using modified Newton iterations.
Figure 3 shows a detailed flowchart of the coupling between the FE and the circulatory model. The four major modules (pressure update, circulatory model, perturbation of circulatory model, and FE model) are independent of one another. Thus, any arbitrary circulatory model (including the simpler Windkessel models) can be “plugged in”, as well as any arbitrarily chosen ventricular model of cardiac (electro-)mechanics (with LV-only or a biventricular anatomic model).
A pressure updating algorithm (Fig. 3a) was needed in order to obtain the correct pressures that ensure continuity of flow between the FE and circulatory models. Ventricular cavity pressures were estimated and updated in several steps (in a variable number of update steps) within each time step:
If the resulting cavity volumes from the finite element and circulation model are not within a preset tolerance, then:
At the beginning (update 0) of each new time step n + 1, cavity pressures are estimated using a fourth order Adams–Bashforth scheme:24
is a vector containing left and right ventricular pressures, and
is the temporal derivative of those pressures for time step n.
For the first time step, the initial conditions in the circulatory model for left and right ventricular cavity volumes are taken from the FE model, and an initial estimate for ventricular pressures. In all following time steps, pressures serve as boundary conditions in the circulatory model, together with initial conditions that are solutions of the previous converged time step. Solving the circulatory model yields new cavity volumes . Next, the FE model is solved with pressures , which yields FE cavity volumes . To decrease calculation time, the Jacobian of the FE problem is updated only with each new time step.
The residual vector is defined as the difference in FE and circulatory cavity volumes:
When the residual, normalized to circulatory cavity volume:
has not satisfied a preset criterion δ = 0.01%, LV and RV pressures are perturbed in the next two pressure update steps.
In this update step, LV pressure is perturbed while keeping RV pressure constant:
and consequently the FE model is solved with these pressure values.
The solution of the previous update step – with perturbed LV pressure – yields new LV and RV cavity volumes. Together with pressure and volume values from update step 0, compliance ΔVLV,FE/ΔpLV and co-compliance ΔVRV,FE/ΔpLV are computed.
Next, RV pressure is perturbed while keeping LV pressure constant:
such that the compliance ΔVRV,FE/ΔpRV and co-compliance ΔVLV,FE/ΔpRV can be computed.
With the previously computed compliances and co-compliances, the FE compliance matrix is assembled.
In each following update step the FE model is solved with updated LV and RV pressures. The LV and RV pressures are updated using Newton’s method, where the residual is minimized:
is the derivative of the residual (Eq. 6) with respect to LV and RV pressures. This represents the system compliance matrix, and are the FE and circulatory compliance matrices, respectively. Note that co-compliances are included off-diagonal in these matrices.
As soon as the FE compliance matrix is fully assembled after update 2, it is kept constant for all following update steps within a time step, avoiding large computation times (updating (co)compliances necessitates the repetition of perturbation steps). Hence, a modified Newton approach is used here.
As soon as the residual is small enough (< 0.01%), a new time step is entered and the pressure update procedure starts over again at update 0.
Ventricular compliances, needed for the circulatory compliance matrix (in Eq. 11), are also obtained by perturbation. In each update step, the circulatory model is solved three times: two times with perturbed estimated pressures in order to assemble the circulatory compliance matrix, and one time with unperturbed pressures, from which the cavity volumes are used to calculate the residual in Eq. 6.
The circulatory compliance matrix is updated every update, because the computation time to solve the circulatory model is insignificant compared to solving the FE model.
First, initial conditions for the FE model were determined by solving the circulatory model when it is coupled to the time-varying elastance model. Parameters of the varying elastance model were fitted to the global (cavity PV) hemodynamic behavior of the FE model, as described in the Appendix.
The FE model was then linearly inflated in 25 steps (uncoupled from the circulatory model) to the initial cavity volumes, as determined from the circulatory model. At this point, the actual target FE simulation was started, fully coupled to the circulatory model. Start of active stress generation in the ventricular myocardium (120 ms after atrial stimulation) was synchronous at a basic cycle length of 600 ms.
Three simulations with the combined cardiovascular model were performed. In the NORM simulation, a normal heart was simulated to investigate when steady state was reached. Steady state was defined as the state where the LV and RV stroke volume difference, normalized to LV stroke volume, was less than 1%, keeping all parameters constant. In the ISCH simulation, ischemia was suddenly induced (Fig. 1) after multiple normal beats. In the PAC simulation the pulmonary artery was constricted.
For every time step, global hemodynamics such as pressures, volumes and flows were computed. Regional values of strains were referred to the end-diastolic state of the preceding beat. For every beat, stroke volume was also calculated. One beat was defined from end-diastole to end-diastole (at mitral valve closure).
The FE anatomic model of the canine heart was discretized into 48 tricubic Hermite elements, with 1968 degrees of freedom. A time step of 4 ms was used.
The non-linear FE model was solved with a modified Newton iteration scheme. Integration was performed with 3 × 3 × 3 Gaussian quadrature points. Convergence was reached when both the maximal incremental displacement solution and the maximal value of the residuals were lower than 10−5. The Jacobian was calculated and factorized in the first iteration of a new time step and when the solution was diverging. The system of linear equations was solved with SuperLU,26 a direct solver optimized for sparse matrices. The circulatory model was integrated in time with the Radau5 solver.15 Absolute and relative tolerances were set at 10−4. The initial guess for the step size was set at 10−3.
The combined multi-scale model was solved with the Continuity 6.3 package (http://www.continuity.ucsd.edu). The three models (NORM, ISCH and PAC) are publicly available in example “biomechanics11” after downloading and installing the latest release of Continuity 6.3. The ISCH simulation was solved on a Linux platform with a 3.2 GHz Intel® Pentium® 4 processor and 3.5 GB of RAM. On average, it took 2.1 min to solve one time step, resulting in a total simulation time of 156 h and 25 min for a total of 30 beats (this includes 12 normal beats initially). Client and server were on separate computers, and data was transferred over the network every time step. The PAC simulation was solved on a Windows XP™ platform (version 2002), on an AMD Sempron™ 3000+ processor at 1.81 GHz and 256 MB of RAM. Solving took 1.8 min per time step on average, resulting in a total simulation time of 84 h and 10 min for a total of 18 beats. Client and server were on the same computer.
Maximum total RAM needed was about 200 MB.
The initial conditions were determined with the ODE solver JSIM 1.6 (http://www.nsr.bioeng.washington.edu/PLN/Members/butterw/JSIMDOC1.6/JSim_Home.stx) on a Windows XP™ platform (version 2002), on an AMD Sempron™ 3000+ processor at 1.81 GHz.
To assess whether meshing was sufficiently fine, the number of transmural elements was doubled and the first normal beat was simulated. Compared with the coarser mesh, maximum absolute deviation in LV and RV cavity volumes was 0.2 and 0.4%, respectively. For LV and RV cavity pressures this was 0.4 and 1.0%, respectively. At the moment of peak LV pressure, the root mean squared difference for myofiber strains was 1.2%. Thus, mesh refinement was sufficient.
To assess the speed and stability of the simulations owing to the co-compliances, a normal beat was simulated without inclusion of co-compliances (as in Refs.2,18). On average, the number of FE iterations needed per time step was 59 (SD = 48) with co-compliances and 71 (SD = 60) without. An attempt to simulate a heart beat with ischemia without using co-compliances failed in the first time step (the FE model diverged) due to an over-estimation of cavity pressures.
Figure 4 shows systemic and pulmonic pressures, volumes, and flows for the first 12 and last 6 beats of the NORM simulation. Early rapid filling occurred with a peak of 0.25 and 0.21 l/s for LV and RV, respectively (Table 1). The atrial kick was more prominent in the RV than in the LV (a small negative flow value, just before it becomes zero and then positive). For every beat, the pulmonary valve opened 28 ms earlier than the aortic valve.
Steady state (with a relative stroke volume difference of < 1%) was reached at beat 21. At the 46th beat, the relative stroke volume difference was 0.06%. Initially, the change in LV stroke volume was larger than that of the RV (Fig. 4). From the PV loops (Fig. 5a) it can be seen that this change is due to a decreasing LV end-diastolic volume, while LV end-systolic volume remains relatively constant.
LV myofiber strains were about −0.07 to −0.10 at end-systole (Fig. 5c). At this time, the RV was contracted more (myofiber strains ~ −0.13). On average, 5.1 pressure updates per time step were needed, within a range from 1 to 7, for a preset value of the volume error convergence criterion δ = 0.01%.
First, 12 normal beats were calculated. After ischemia was induced during filling of the 12th beat (at t = 7200 ms), LV systolic pressures decreased (Table 1, Figs. 5a, ,6a)6a) and LV stroke volume dropped by 50% (Table 1, Fig. 6a) due to an increase of end-systolic LV volume (Fig. 5a). Both LV and RV PV loops shifted rightward to larger end-diastolic volumes and lower maximum pressures as compared to the loops from the normal heart.
In the center of the ischemic region, myofiber strains were less negative, changing from −0.07 (in the normal heart) to about zero (Fig. 5c) in the first beat with ischemia (beat 13). At this moment, RV contraction was larger than the LV, but smaller compared to the NORM simulation at end-systole. The pulmonary valve opened 44 ms before the aortic valve (Table 1). In the last beat, LV and RV myofiber strains were more negative compared to the ones in beat 13, while strains in the ischemic region were ~ −0.04. In the last beat, the pulmonary valve opened 24 ms before the aortic valve. The simulation with ischemia did not affect the number of pressure updates needed (average 5.1 per time step) for volume convergence.
After PAC was applied at t = 7200 ms, LV systolic pressure initially increased (Table 1, Figs. 5b, ,6b)6b) but returned to normal in the last beat. LV stroke volume also initially increased (Table 1, Fig. 6b) due to a decrease of end-systolic LV volume (Fig. 5b). Stroke volume then decreased to a value lower than normal in the last beat due to a decrease of LV end-diastolic volume.
Peak RV pressure increased by 185% in the first beat with PAC compared to the normal heart, and gradually increased further by 195% in the last beat (Table 1, Fig. 6b) compared to the normal heart. RV end-diastolic pressure and volume increased also. RV stroke volume initially dropped, but gradually increased again over time, settling at a value below normal. The LV PV loops shifted to the left (Fig. 5b), while RV PV loops shifted to the right, compared to the loops from the normal heart.
At end-systole in the first beat with PAC, RV myofiber strains were less negative than in the normal heart (Fig. 5c). RV contraction was again larger than the LV. Over time, change in end-systolic myofiber strains became more negative. Opening of the pulmonary valve occurred earlier up to 36 ms before the aortic valve in the last beat (Table 1), presumably due to the higher RV dp/dtmax, compared to normal. The septum was shifted towards the LV.
A new multi-scale method was developed that couples 3D finite element models of ventricular cardiac mechanics to lumped models of the circulation. This method enables realistic hemodynamic boundary conditions to be used in all four cardiac phases and multiple cardiac cycles to be computed. Transient alterations in global hemodynamics and regional mechanics were seen after applying two interventions, one in the FE model (ischemia) and one in the circulatory model (pulmonary artery constriction). The method is robust, enabling fast convergence, due to inclusion of co-compliances to update pressures.
The cardiovascular model was set up in a modular fashion. The FE model, pressure update algorithm, circulatory model and the perturbation of the circulatory model are all separate. For example, the FE model only needs left and right ventricular pressure as hemodynamic boundary conditions and does not need information about other pressures or volumes in the circulation. The circulatory model does not need information about material properties (e.g. a region with inhibited contraction due to ischemia) or shape (e.g. dog, rabbit, or pig) of the FE heart model. A modular setup enables the easy exchange of circulatory models and FE models.
The function of the model was demonstrated by coupling a FE model of a dog heart to a closed systemic and pulmonic circulation. The atria were also included as time-varying elastance chambers, relating atrial pressure to atrial volume and time.
The method for cardiovascular coupling presented in this paper has several advantages over previously published methods.2,18,19,40 First, the method presented here is independent of cardiac phase: the opening and closing of the valves is determined by the circulatory model. Phase-independence opens up many possibilities:
Second, in FE models with a left and right ventricle, fast convergence is achieved by taking into account compliances and ventricular interaction (through co-compliances) from the FE model and compliances from the circulatory model. In previous methods ventricular interaction was not taken into account, which resulted in many pressure update steps and/or instability.
Perturbation of the FE model is needed in order to compute the compliances. Fortunately, FE perturbation does not substantially increase computation time, since the Jacobian is only assembled and factorized at the beginning of a time step, not at every update. Assembly and factorization of the Jacobian require the most time in FE simulations. A simulation in which the Jacobian was updated for every pressure update step took 3.2 times longer to compute.
Normal global results (pressures, volumes, flows) were similar to typical measured values (e.g. Verbeek et al.43). Immediately after the ventricles were activated, the atrioventricular valves closed. The pulmonary valve opened before the aortic valve. RV systolic pressure increased and decreased before LV pressures, which is also in agreement with experiments.43 The opening of the pulmonary valve before the aortic valve is reflected in the magnitude of end-systolic myofiber strains, which was to be expected from the volume data. Hence, myofibers shorten more in the RV free wall than in the rest of the heart.
Immediately after inducing ischemia, LV and RV stroke volume reduced by 50 and 10%, respectively. The former was due to a loss of active force generation in the ischemic region. The latter is a manifestation of direct ventricular interaction, i.e. a change in LV function modulates RV function through the walls and septum.36 As can be seen in Fig. 2b, a lower LV systolic pressure resulted in a higher RV volume. Thus, end-systolic RV volume was increased, lowering RV stroke volume. LV PV loops shifted to the right, in agreement with experiments.17 Because outflow of the RV was larger than that of the LV (LV outflow was reduced), pulmonic pressures and volumes were increasing, and thus preload of the LA and LV increased. This in turn increased LV end-diastolic pressures and volumes (which results in the right-shift of the loops). Because aortic pressure was decreasing, the aortic valve closed at lower pressures, decreasing the end-systolic volume. Therefore, in the last beat, LV stroke volume was only reduced by 18% compared with normal.
The RV PV loops also shifted to the right. This phenomenon has also been shown experimentally.17 Owing to reduced LV outflow, systemic volumes were decreasing. But because RV outflow was also reduced, flow into the RA was still larger than outflow of the RA, increasing RA volume. This increased the RV preload and the RV PV loops shifted to the right.
Myofiber strain in the ischemic region at end-systole was larger than in the non-ischemic part of the LV and compared to the normal heart, which qualitatively agrees well with experimental data1,28,45,46 and previously performed simulations in the LV.1,28,46 However, in the ISCH simulation these strains were about zero, becoming more negative in the last beat while in the experiments positive values were reported. These quantitative differences may be explained by the difference in size of the ischemic region, material properties, and cavity pressures.
In the PAC simulation, the RV PV loops shifted to the right, while those for the LV shifted leftwards, which has been shown in experiments.10,22,31 The cause of the RV PV loop right-shift is similar to the situation with the ischemic LV. The higher pulmonary valve impedance causes a lower RV stroke volume. LV stroke volume is still close to normal, hence the venous return increases the RV preload.
Interestingly, the LV stroke volume initially increased. This is again explained by direct ventricular interaction. From Fig. 2b (phase 2) it appears that an increase in RV pressure decreases LV volume. Hence, end-systolic LV volume was lower and stroke volume larger, so the higher RV pressures assist the LV to eject blood. Later, however, owing to the decreased RV stroke volume, LV preload decreased and therefore LV stroke volume also decreased.
A change in the transseptal pressure gradient (due to a higher RV pressure than normal) resulted in a movement of the septal wall towards the LV, in agreement with experiments22,31 and 2D FE models.12,31
Special attention was paid to setting initial conditions with a simpler (and thus faster) zero-dimensional time-varying elastance model that replaced the computationally expensive 3D FE model. Even with these initial conditions, steady state was not immediately obtained when applied in the 3D FE model coupled to the circulatory model. The elastances of the 3D model are non-linear, especially in the passive state. Hence, the strategy of obtaining initial conditions could be improved by also taking these non-linearities into account in the varying elastance model.
The mechanical effects of ischemia were simplified by decreasing the peak of intracellular calcium bound to troponin C in the active stress model. This resulted in a ~80% decrease in active stress. Events in cardiac myocytes undergoing ischemia are much more complicated than assumed here,5 and more elaborate electrophysiological computational studies of ischemia have been developed.34,35,38 Dynamic whole heart FE models of cardiac mechanics have been used to investigate strains in more detail during regional ischemia.1,28 The results in these studies were acute in terms of one cardiac cycle after ischemic occurrence. Although the methods presented in our paper could be used to gain information on strains during ischemia in a steady state, the aim of the simulation of regional ischemia was to demonstrate the dynamic and modular nature of the method.
This cardiovascular model lacks neurohumoral stimulation. The in vivo heart is part of an integrated system that compensates for a fall in blood pressure (e.g. due to ischemia). When arterial blood pressure falls, neurohumoral stimulation increases cardiac output and blood supply is redistributed throughout the body. A gradual drop in aortic pressure would be compensated for through feedback from the baroreceptors, increasing contractility and beating frequency. Integrated numerical models containing feedback from the baroreceptors,39 signal transduction,37 cardiac (electro)mechanics, and the hemodynamics of the circulation will allow a more complete investigation of regulation of myocardial contractility on a multi-level scale.
This heart model also lacks a pericardium, which would further modulate ventricular interaction48 and stroke volumes, especially during abnormal conditions.13 Although ventricular interaction still exists in the absence of the pericardium,29,48 effects of ventricular interaction could be enhanced by inclusion of pericardial pressure.6,12
Unloaded volumes for the arterial/capillary and venous Windkessel segments in the circulatory model were assumed to be zero. Coupled with the fact that model parameters must account for significant volumes of blood in each segment, this required high compliances throughout the circulation. This is a limitation related to the low number of Windkessel segments in our simplified circulatory model. We changed the unloaded volume of the systemic arterial/capillaries segment from 0 to to 98 ml (0.28 of its normal volume7) to see how it would affect the circulatory model driven by the varying elastance model. To keep the correct distribution of blood volumes, compliance Cas had to be reduced to 23.4 ml/kPa. Pulse pressure only increased by 0.2 kPa when these changes were made. This illustrates the fact that in order to produce realistic pulse pressures and volume distributions, more detail is needed in the circulatory model – i.e. the arterial/capillary Windkessel needs to be divided into several more Windkessel segments.
The parameters controlling atrial contraction were hand-tuned so that when the target cardiac output was reached, average atrial volumes would fall within 20% of their target volumes in the steady state and the atrial “kick” would be physiologically realistic. We attempted to use elastance and unloaded volume values from Lau et al.25 for the atria. However, these values produced an unrealistically high atrial kick (40–50% of the stroke volume) and also produced average atrial volumes that were too low in the steady state. Current parameters have been chosen so that atrial volumes match the mean target values for volumes and pressures and provide a realistic contribution to ventricular filling. Furthermore, the effects of breathing were not modeled here. If a periodic, negative pleural pressure was applied external to the heart chambers, then elastance parameters in the atria would have to be increased to produce the target atrial pressures and volumes.
A novel method was developed that couples anatomically detailed 3D finite element models of ventricular cardiac mechanics to lumped parameter systems models of the systemic and pulmonary circulations. The proposed multi-scale modeling method is cardiac phase-independent and allows simulations of multiple cardiac cycles. The method is efficient and robust. Models of ischemia and pulmonary artery constriction agreed with experiments and enabled the investigation of direct and indirect ventricular interactions.
This work was supported by the National Biomedical Computation Resource (NIH Grant P41 RR08605) (to A.D.M), National Science Foundation Grants BES-0096492 and BES-0506252 (to A.D.M) and BES-0506477 (to M.L.N.), NIH Grant HL32583 (to J.H.O.), and NIH Grant EB001973 (to J.B.B.). This investigation was conducted in a facility constructed with support from Research Facilities Improvement Program Grant Number C06 RR-017588-01 from the National Center for Research Resources, National Institutes of Health. A.D.M. and J.H.O. are co-founders of Insilicomed Inc., a licensee of UCSD-owned software used in this research. Furthermore, we are grateful to our programmers Sherief Abdel-Rahman, Ryan Brown, and Fred Lionetti for their excellent work on improving and extending Continuity.
In the heart, the relation between ventricular volumes and pressures is written as:
In which Δ is a vector with LV and RV instantaneous ventricular volumes minus the rest volumes:
where yv is a ventricular activation function:
where Vx,rest,d and Vx,rest,s are diastolic and systolic unloaded volumes.
Matrix C̲ is the ventricular time-varying compliance matrix:
C̲max and C̲min are compliances for the fully active and passive state, respectively.
From the pressure and volume curves in Fig. 2 it can be seen that ventricular co-compliances are pressure-dependent: i.e. the RV volume change for a LV pressure change is different at a constant low and high RV pressure. Hence CLR and CRL in Eq. A1 are written as a function of pressure (see also Table A1).
The same procedure is performed for maximally activated ventricles (Fig. 2B).
Using the time-varying elastance model in the more common way (with input volume and output pressure), Eq. A1 is rewritten as:
The atrial elastances are driven by an activation function
Left atrial pressure is given by
where LA elastance and rest volume (volume at zero pressure) are given by
Right atrial pressure, elastance, and rest volume are given by
Notice that ventricular volume changes (Eqs. B16 and B26) are purely determined by ventricular in- and outflows. In the case of coupling between the FE and circulatory model, cavity pressures come from the update algorithm. In the case of the procedure for determining the initial conditions, cavity pressures are calculated by the ventricular time-varying elastance model.