Learn more: PMC Disclaimer | PMC Copyright Notice
Left Ventricular Finite Element Model Bounded by a Systemic Circulation Model
Short abstract
A series of models were developed in which a circulatory system model was coupled to an existing series of finite element (FE) models of the left ventricle (LV). The circulatory models were used to provide realistic boundary conditions for the LV models. This was developed for the JSim analysis package and was composed of a systemic arterial, capillary, and venous system in a closed loop with a varying elastance LV and left atria to provide the driving pressures and flows matching those of the FE model. Three coupled models were developed, a normal LV under normotensive aortic loading (116/80 mm Hg), a mild hypertension (137/89 mm Hg) model, and a moderate hypertension model (165/100 mm Hg). The initial step in the modeling analysis was that the circulation was optimized to the end-diastolic pressure and volume values of the LV model. The cardiac FE models were then optimized to the systolic pressure/volume characteristics of the steady-state JSim circulatory model solution. Comparison of the stress predictions for the three models indicated that the mild hypertensive case produced a 21% increase in the average fiber stress levels, and the moderate hypertension case had a 36% increase in average stress. The circulatory work increased by 18% and 43% over that of the control for the mild and moderate hypertensive cases, respectively.
Introduction
Finite element models have been applied to study the effects of ischemia [1,2] and infarction [3–5] on the left ventricle as well as used to evaluate the effectiveness of cardiac interventions such as the repair of infarction induced aneurysms [6–9]. One of the limitations of these types of FE models is that they have not been connected to a model of the circulatory system so that the reflections of the pressure waves back to the heart from the arterial system and their influence on the contractile process are not represented. Instead, a priori boundary conditions are assigned to the models based on typical vascular pressure loading to predict the cardiac strain/stress.
In order to overcome this problem, Kerckhoffs et al. [10] embedded a full zero-dimensional (0D) circulatory system into the Continuity simulation package1, thus providing the appropriate boundary conditions to a 3D FE biventricular mechanical model [11,12]. The tight integration of the circulation with the solid model allows for direct exchange of LV pressure/volume information between the two parts of the simulation package, but this system is computationally expensive.
The primary goal of the following work was to create a framework that allows for the communication of pressure and volume values back and forth between two distinct and dissimilar modeling systems, thus linking a 0D circulatory model to a 3D finite-element LV mechanical model. The circulatory model was run under the JSim [13] analysis package2 while the LV models were developed for and analyzed using the general purpose, nonlinear, large deformation FE package NIKE3D [14]. Due to limitations in the communication pathway between the circulatory and FE models, this type of coupling can be considered “weak” or “loosely” coupled rather than the direct coupling illustrated by the work of Kerckhoffs et al. [10].
Materials and Methods
The coupled JSim/FE models were composed of three parts, each of which will be described in detail below. The first component was a circulatory system run using JSim. The second component was the FE LV mechanical model that runs under NIKE3D. The third component, the JSim/FE interface, controls the communication of information between these two modeling systems as well as providing the means to optimize the FE LV active contraction parameters of the FE model so that it reproduces the pressure/volume characteristics of the JSim simulation.
JSim Circulatory Model.
The JSim model of the systemic circulation contains the left atria, left ventricle, the aorta, and the rest of the systemic circulation (arteries, capillaries, and venous return) (Fig. (Fig.1)1) and is a reduced form of the Neal and Bassingthwaighte [15,16] model. The systemic circulations were modeled using three lumped Windkessel compartments in series, one compartment for the aorta, another for the arteries and capillary blood, and another for the venous blood. The right heart and the pulmonary systems were lumped into the venous return. The left atria and left ventricle are represented using time-varying elastance models. All the equations are listed in Appendix A.
Schematic of the systemic JSim model of FE heart model and model of circulatory system. The labels are P for the pressure values, F for the flow, R for the resistance, C for the compliance, and L for the inertance at the given locations shown above. The left atria and ventricle are represented as a time varying elastance units.
Left Ventricle Mechanical Model.
A complete description of the LV model can be found in Veress et al. [1,5] (Fig. (Fig.2).2). Briefly, the passive myocardium was represented as transversely isotropic [17] with a time varying elastance active contraction component [18,19].
The total Cauchy stress in the fiber direction a is given as and is the sum of the active stress T(a) and the passive stress tensor T(p) component due to the transversely isotropic material model in the fiber direction as follows:
The active fiber stress T(a) is defined as
where T max = 135.7 KPa is the isometric tension under maximal activation at the peak intracellular calcium concentration of Ca0 = 4.35 μM [19]. Ct governs the shape of the activation curve [19], which is only defined during contraction. The product of the constant Tmax and Ct defines the primary boundary condition in the mechanical model input file and was zero during diastole and controls active contraction during systole. The length dependent calcium sensitivity ECa 50 is governed by the following equation:
where (Ca 0)max = 4.35 μM is the peak intracellular calcium concentration, B = 4.75 μm−1 governs the shape of the peak isometric tension-sarcomere length relation, l 0 = 1.58 μm is the sarcomere length at which no active tension develops, and l is the sarcomere length, which is the product of the fiber stretch ratio , and the unloaded length lr = 2.04 μm [18,19].
Coupled Analysis.
The coupled analysis begins with the analysis of the baseline FE LV model in order to provide the JSim circulatory model with diastolic pressure/volume characteristics based on literature values [19,20]. The end-diastolic pressure/volumes were passed to JSim where a SENSOP optimization routine [21,22] was utilized to tailor the circulatory model so that the diastolic pressure and volume characteristics of the FE model were reproduced in the circulatory model. This was accomplished through the optimization of the systemic resistance (Rar) and capacitance (Car) values and the maximum elastance of the left ventricle (Emaxlv).
Upon completion of the diastolic optimization, the circulatory model was run until equilibrium was achieved (five cardiac cycles). JSim then initiated the interface program by issuing a command to the Linux operating system. The interface program provided the communication pathway between the JSim circulatory model and the ventricular FE mechanical model.
JSim Interface Program.
The interface program was composed of two parts. The first part reads in pressure and volume values at four points during the ejection phase of systole exported from JSim as a text file. These were 408, 413, 417, and 420 ms from the beginning of diastole (reference configuration t = 0 with an 800 ms cardiac cycle time). Peak systolic pressure was achieved at 417 ms. These time points represent the ejection phase of systole when the aortic valve is open. The NIKE3D LV model file was read by the interface program, and the JSim-derived LV pressure values were substituted for those of the original Nike model. The interface program initiated the FE analysis of the model, which was the starting point for the subsequent systolic optimizations. For this initial run, the FE analysis of the entire cardiac cycle was completed with the analysis starting with the model at reference, beginning diastolic configuration, proceeding through passive filling followed by contraction, and finally relaxation from the contracted state.
The active contraction stress values of the LV model were then optimized until the difference in the NIKE3D volumes (Appendix B) and the JSim volumes reached a user defined tolerance. This was accomplished using a secant [24] type iteration scheme (Eq. (4)), which uses an estimate of the derivative based upon the current and previous iterations (Eq. (5)). The derivative of the active contraction stress/error function was used to determine the next active contractile stress values via the setting of Ct as given below:
where n is the current iteration, n – 1 is the previous iteration, and n + 1 is the next iteration. α is a damping factor having a value of 0.5. The damping factor was necessary due to the highly nonlinear nature of the stress/volume relationship. The derivative f′ was defined as:
with f being the error in volumes
The iteration scheme ran until the user defined tolerance was reached, (1.5 ml). Then the analysis for this time point was completed, and the analysis would proceed to the next time point. The optimization process is shown schematically in Fig. Fig.3.3. In order to save computational time, the FE analysis used for each optimization iteration was completed up to and including the time point being optimized and then was subsequently terminated.
Schematic of the coupled system analysis protocol. The initial step is the optimization of the circulation to the end-diastolic pressure and volume values produced by an initial run of the FE model. JSim was then run until equilibrium was achieved in the model. This was followed by optimization of the FE model to reproduce the JSim pressure and volume values at each of the four time points. Once the correct values were achieved for a given time point, the process was repeated for the next time point until all of the time points had been optimized.
Application: Normal and Hypertensive Hearts.
Three cases were developed and analyzed with the first case being a normal LV model coupled to the circulation model that represented normotensive loading (116/80 mm Hg) measured at the aorta. Two pathologic cases were modeled with the first being mild hypertension (137/89 mm Hg) and the second being moderate hypertension (160/100 mm Hg). The hypertension circulatory models were created by systematically increasing the peripheral resistance until the hypertension values were achieved without an increase in the ejection fraction. The geometries of all of the FE models were identical, thus representing cases of hypertension without subsequent remodeling. The starting parameters for the normotensive systolic optimizations were based upon a previous LV model [1]. Additionally, the hypertension cases used the optimized normal FE model as the starting point for the systolic optimizations. The entire cardiac cycle was analyzed following the optimization of the final systolic time point.
FE Model Data Analysis.
Upon completion of the analysis of the three cases, the average fiber stresses predicted by each of the FE models were compared. The circulatory pressure and volume characteristics were used to estimate the circulatory work.
Results
The diastolic optimizations within JSim were completed in 1 min for 150 iterations using the SENSOP optimizer run on a four core 3.2 GHz Linux computational server. The optimized diastolic LV pressure value fell within 1.0 mm Hg of the NIKE3D _target value, and the end-diastolic LV volume came within 1.0 ml of the _target value for all cases.
The systolic optimizations run within the interface program (Figs. 4(a) and 4(b)) required four NIKE3D iterations to ensure proper LV volume output (Fig. (Fig.5).5). The total analysis time was approximately 4 h. The curves labeled “starting values” in Fig. Fig.55 are the NIKE3D volume values on the initial run following the substitution of the JSim systolic pressure values into the LV model.
(a) Aortic pressure curves (dashed black) and the left ventricular pressure curve (solid black) defined in the JSim circulatory system are given. (b) The JSim derived LV volume curve (black) used in the mild hypertensive systolic optimizations. The vertical gray lines indicate the four systolic time points for both graphs.
Optimization to the JSim systolic volume values (gray) was achieved after four iterations all of the time points in each of the cases: normotensive (top), mild hypertension (middle), and moderate hypertension (bottom)
The circulatory model was able to produce realistic hemodynamic values (Table (Table1).1). The cardiac output (stroke volume) increased in the mild and moderate hypertension cases compared with the normal model. The hypertension cases showed unchanged ejection fraction values from that of the normal case.
Table 1
Circulatory hemodynamic data
| Normal | Mild Hypertension | Moderate Hypertension | |
|---|---|---|---|
| EDV | 128 | 129 | 131 |
| ESV | 60 | 60 | 60 |
| SV | 69 | 70 | 71 |
| EF | 53% | 54% | 54% |
Note: EDV = end-diastolic volume, ml; ESV = end-systolic volume, ml; SV = stoke volume, ml; EF = ejection fraction = (EVD – ESV)/EVD, dimensionless.
The FE analysis prediction for the average LV fiber stress was 17.9 KPa for the normotensive case. The mild hypertensive case had an average fiber stress of 21.6 KPa, and the moderate hypertensive case had a value of 24.3 KPa (Fig. (Fig.6).6). These represented increases of 21% and 36% in end-systolic fiber stress.
Hypertension increases the fiber stress as seen in these ventricular short axis slices. These are fiber strain distributions for normal, mild, and moderate hypertension (left). The location of the slices is shown on the right.
The circulatory work indicated that the increase in peripheral resistance substantially increased the work of the heart. The circulatory work for the baseline normal case was 0.857 J, 1.015 J for mild hypertension, and 1.223 J for the moderate hypertension case. These represent 18% and 43% increases in workload, respectively.
Discussion
The use of the 0D JSim model to define the boundary conditions for the LV models resulted in a dual success, namely the development and parameterization of an FE model with constrained time-dependent physiological outflow impedance and a 0D equivalent circulatory model defining systemic arterial pressures and flows. The system in its present form is capable of modeling and having the FE LV model respond to changes in the circulatory system through altered loading of the FE model.
The analysis of the normotensive and the hypertensive cases indicate that even mild hypertension can cause a marked increase in total LV wall stress as demonstrated by the fiber stress results. The circulatory work results were consistent with these findings indicating that relatively mild increases in afterload resulted in substantially increased circulatory work values.
The optimization routines utilized in this study proved to be efficient. The SENSOP method used in this study was robust and converged quickly for the diastolic optimization of the circulatory system. The systolic optimization (secant method) showed reasonable convergence characteristics provided the starting pressure values in the FE model were within 5 mm Hg of the JSim-defined pressures. Outside this range at the end-systolic time point, the optimization had difficulty converging as the iteration solutions tended to oscillate between positive and negative errors or simply did not converge.
The presented work is similar in scope with the work of Wenk et al. [24] in that a weak coupling was sought between a circulatory model and a FE based model. Wenk et al. initially tune the FE models to pressure/volume relationships found in their animal models. The models were then used to generate a family of pressure/volume curves through the definition of end-systolic pressure/volume and the end-diastolic pressure/volume relationships. These functions are passed to the circulatory model in order to tune the resistance and capacitance components. The strength of the Wenk model was that the circulatory system was tuned using a more complete relationship rather than the end-diastolic volumes and pressures as was done in our methodology. In systole, we simply let the circulatory system define the entire systolic pressure/volume relationship and then matched the FE model volumes using the pressures defined by the circulatory system. The use of discrete points to characterize the systolic contractile phase allows for optimization of experimental systolic pressure/volume data directly rather than fitting to a pressure/volume relationship.
Limitations.
The primary limitation of the current implementation of this type of coupled modeling is that the transfer of information at the diastolic and four systolic time points is unidirectional instead of being bidirectional at all times. The communication necessary for heart failure and ischemia/infarction needs to be bidirectional at both diastole and systole so that the deficiencies in LV function will alter the response of the circulatory system.
Another drawback was that changes cannot be made to the circulatory parameters within a single cardiac cycle for either analysis system due to limitations inherent to both JSim and NIKE3D. The changing of parameters during the runs within JSim causes numerical instabilities that greatly alter the pressure/volume curves leading to erroneous results. While NIKE3D does have a mechanism by which the run may be interrupted and restarted, the boundary conditions cannot be changed for the continuation of the simulation.
The pressure wave forms produced by the circulatory system model tended to be overly simple. For example, the dicrotic notch does not appear in the aortic pressure waveform (Fig. (Fig.4).4). However, overall the waveforms do capture the relevant features of the LV pressure/volume relationship.
This work demonstrates a methodology used to couple two completely separate and dissimilar simulation packages, providing results without changing the underlying source code of either system. The system contains relatively generic components such that the interface program could be modified and applied to other dissimilar systems. With continued development, the system will be able to be used to model the boundary conditions for models that are used to define circulatory based pathologies.
Acknowledgment
We would like to acknowledge these sources of funding: NIH Grant Nos. R01-EB08407 for the JSim development, T15-HL088516 for the model archiving at www.physiome.org, and R03 EB008450, R01 EB07219, and R01 EB000121, and the Director, Office of Science, Office of Biological and Environmental Research, Biological Systems Science Division of the US Department of Energy under Contract No. DE-AC02-05CH11231 for FE model development.
Appendix A: Equations Defining the JSim Model
The equations used to describe the function of the JSim 0D circulatory system are given below.
Pressure
Pressure = elastance (volume minus resting volume) (in mm Hg)
Flow
Flow = pressure drop/resistance (in l/min)
ΔFlow equals pressure drop/inertance
Parameter Values for Normotensive Case
Emaxla = 2.07 mm Hg/ml, maximum left atrium elastance
Eminla = 0.15 mm Hg/ml, minimum left atrium elastance
Emaxlv = 5.01 mm Hg/ml, maximum left ventricle elastance
Eminlv = 0.16 mm Hg/ml, minimum left ventricle elastance
Cao = 0.14 ml/mm Hg, compliance or capacitance of aorta
Car = 1.32 ml/mm Hg, compliance of systemic arteries
Cvein = 42.1 ml/mm Hg, compliance of systemic veins
restVlad = 43.9 ml, end-diastolic rest volume left atrium
restVlas = 54.4 ml, end-systolic rest volume right atrium
restVlvd = 56.2 ml, end-diastolic rest volume left ventricle
restVlvs = 37.6 ml, end-systolic rest volume right ventricle
restVao = 22.8 ml, rest volume of aorta
restVar = 184.1 ml, rest volume of systemic arteries
restVvein = 593.3 ml, rest volume of systemic veins
Rmitvalve = 0.003 mm Hg*s/ml, resistance of mitral valve
Raorvalve = 0.001 mm Hg*s/ml, resistance of aortic valve
Lao = 1.e-8 mm Hg*s2/ml, inertance of aorta
Rar = 0.78 mm Hg*sec/ml, resistance of systemic arteries
Rvein = 0.26 mm Hg*sec/ml, resistance of systemic veins
HeartRate = 1.25 Hz, heart rate (75 beats/min)
PRint = 0.16 s, P-R interval
ylaoffset = 0.5 s, offset for left atrium activation
ylvoffset = 0.0 s, offset for left ventricle activation
Variable Definitions
Ela(t) mm Hg/ml, elastance of left atrium
Elv(t) mm Hg/ml, elastance of left ventricle
yla(t) dimensionless, activation function for artrial elastance
ylv(t) dimensionless, activation function for ventricle elastance
restVla(t) ml, rest volume of left atrium
restVlv(t) ml, rest volume of left ventricle
Pla(t) mm Hg, pressure in left atrium
Plv(t) mm Hg, pressure in left ventricle
Pao(t) mm Hg, pressure in aorta
Par(t) mm Hg, pressure in systemic arteries
Pvein(t) mm Hg, pressure in systemic veins
Fmit(t) liters/min, flow in left atrium
Faov(t) liters/min, flow in left ventricle
Far(t) liters/min, flow in systemic arteries
Fvein(t) liters/min, flow in systemic veins
Initial Conditions of State Variables (Normotensive Case), t = 0 s
Fao(0) = 0.30 liters/min, flow in aorta
Vla(0) = 59.1 ml, volume of left atrium
Vlv(0) = 125.9 ml, volume of left ventricle
Vao(0) = 34.5 ml, volume of aorta
Var(0) = 327.5 ml, volume of systemic arteries
Vvein(0) = 1910.1 ml, volume of systemic veins
State Variables for Normotensive Case
Fao(t) liters/min, flow of aorta
Vla(t) ml, volume of left atrium
Vlv(t) ml, volume of left ventricle
Vao(t) ml, volume of aorta
Var(t) ml, volume of systemic arteries
Vlv(t) ml, volume of systemic veins
Appendix B: LV Volume Calculation
A variation of the “stacked coin” methodology used clinically to determine LV volumes from echocardiographic images was utilized [25]. The endocardial nodes in the FE LV model corresponding to short axis slices were used to determine an effective radius through fitting a circle to the endocardial node points using a least squares approximation [26]. The volume between each of these layers was approximated as a truncated circular cone so that each incremental volume was determined by
where R 1 is the radius of a layer of endocardial nodes, and R 2 is the radius of an adjacent set of endocardial nodes (next layer). h is the axis distance between the layers as determined by taking the difference between the average of the long axis coordinates for a given layer and the average of the adjacent layer. The incremental volumes were summed to determine the total LV lumen volume. However, given the modular nature of the interface program, other measures of LV volume could easily be substituted for this one.
Footnotes
1For more information, see www.continuity.ucsd.edu.
2For more information, see www.physiome.org.
Contributor Information
A. I. Veress, Department of Mechanical Engineering,
Department of Bioengineering,
University of Washington,
Seattle, WA 98195.
G. M. Raymond, Department of Bioengineering,
University of Washington,
Seattle, WA 98195.
G. T. Gullberg, Life Science Division,
E. O. Lawrence Berkeley National Laboratory,
Berkeley, CA 94720;
Department of Radiology,
University of California San Francisco,
San Francisco, CA 94122.
J. B. Bassingthwaighte, Department of Bioengineering,
University of Washington,
Seattle, WA 98195.