A Computationally Efficient Formal Optimization of Regional Myocardial Contractility in a Sheep with Left Ventricular Aneurysm

doi:10.1115/1.3148464

J Biomech Eng. Author manuscript; available in PMC 2009 December 15.

Published in final edited form as:

J Biomech Eng. 2009 November 1; 131(11): 111001.

doi: 10.1115/1.3148464

PMCID: PMC2793686

NIHMSID: NIHMS159053

A Computationally Efficient Formal Optimization of Regional Myocardial Contractility in a Sheep with Left Ventricular Aneurysm

Kay Sun,^1,⁴ Nielen Stander,⁵ Choon-Sik Jhun,^1,⁴ Zhihong Zhang,^1,⁴ Takamaro Suzuki,^1,⁴ Guan-Ying Wang,^3,⁴ Maythem Saeed,³ Arthur W. Wallace,^2,⁴ Elaine E. Tseng,^1,⁴ Anthony J. Baker,^3,⁴ David Saloner,^3,⁴ Daniel R. Einstein,⁶ Mark B. Ratcliffe,^1,⁴ and Julius M. Guccione^1,⁴^*

¹ Department of Surgery, University of California, San Francisco, CA

² Department of Anesthesia, University of California, San Francisco, CA

³ Department of Radiology, University of California, San Francisco, CA

⁴ Department of Veterans Affairs Medical Center, San Francisco, CA

⁵ Livermore Software Technology Corporation, Livermore, CA

⁶ Biological Monitoring and Modeling, Pacific Northwest National Laboratory, Olympia, WA

^* Corresponding Author: Julius M. Guccione, Ph.D., UCSF/VA Medical Center (112D), 4150 Clement Street, San Francisco, CA 94121, Email: GuccioneJ/at/surgery.ucsf.edu

See other articles in PMC that cite the published article.

Abstract

A non-invasive method for estimating regional myocardial contractility in vivo would be of great value in the design and evaluation of new surgical and medical strategies to treat and/or prevent infarction-induced heart failure. As a first step towards developing such a method, an explicit finite element (FE) model-based formal optimization of regional myocardial contractility in a sheep with left ventricular (LV) aneurysm was performed using tagged magnetic resonance (MR) images and cardiac catheterization pressures. From the tagged MR images, 3-dimensional (3D) myocardial strains, LV volumes and geometry for the animal-specific 3D FE model of the LV were calculated, while the LV pressures provided physiological loading conditions. Active material parameters (T_{max_B} and T_{max_R}) in the non-infarcted myocardium adjacent to the aneurysm (borderzone) and in myocardium remote from the aneurysm were estimated by minimizing the errors between FE model-predicted and measured systolic strains and LV volumes using the successive response surface method for optimization. The significant depression in optimized T_{max_B} relative to T_{max_R} was confirmed by direct ex vivo force measurements from skinned fiber preparations. The optimized values of T_{max_B} and T_{max_R} were not overly sensitive to the passive material parameters specified. The computation time of less than 5 hours associated with our proposed method for estimating regional myocardial contractility in vivo makes it a potentially very useful clinical tool.

Keywords: tagged magnetic resonance imaging, finite element modeling, numerical optimization, cardiac mechanics

INTRODUCTION

A non-invasive method for estimating regional myocardial contractility is a “holy grail” of cardiology. Such a method would be of great value in the design and evaluation of new surgical and medical strategies to treat and/or prevent infarction-induced heart failure. Once the constitutive equation for the myocardium is established the effect of therapeutic changes on regional geometry (i.e., surgical remodeling) and/or material properties (i.e., medicine, gene therapy, cell therapy) can be evaluated and the success or failure of a proposed therapy predicted. With clinical experience, such a method could be used as a diagnostic modality to risk stratify patients early after a myocardial infarction who are at risk for adverse remodeling and the development of heart failure.

It is impossible to determine myocardial material properties (in the form of three-dimensional constitutive equations) from ventricular pressure-volume relations alone. Instead, our laboratory recently used cardiac catheterization, magnetic resonance (MR) imaging with myocardial tissue tagging [1], MR diffusion tensor imaging [2], and a finite element (FE) method (developed specifically for cardiac mechanics [3]) to measure regional systolic myocardial material properties in the beating hearts of four sheep with left ventricular (LV) aneurysm [4] and six sheep with LV aneurysm repaired surgically [5]. With knowledge of these myocardial material properties, we were able to quantify the effect of aneurysm plication on regional myocardial stress distributions. The 3D stress distributions in the myocardium are important to regional ventricular function because both regional coronary blood flow [6] and myocardial oxygen consumption [7] are influenced by ventricular wall stress [8]. Changes in ventricular wall stress are believed to be stimuli for hypertrophy and remodeling [9]. There have been no successful methods developed to measure stress in the intact heart wall—primarily because of its large deformations and the tissue injury caused by implanted transducers [10].

Although our previous studies [4–5] represent significant advancements in FE modeling of hearts with myocardial infarction, because of long computation times, they both employed a manually directed pseudo-optimization. In other words, a formal nonlinear optimization of material constants was not feasible. The objective of the present study is to describe a method to formally optimize regional myocardial contractility in vivo that is at least an order of magnitude faster than that used in our previous studies.

METHODS

The sheep used in this study was treated in compliance with the “Guide for the Care and Use of Laboratory Animals” prepared by the Institute of Laboratory Animal Resources, National Research Council, and published by the National Academy Press (revised 1996).

Experimental Measurements

Data collected from one male adult sheep [11] was used to demonstrate methodology and accuracy of the FE optimization tool. Briefly, the sheep underwent anteroapical myocardial infarct following the procedures described in Markovitz et al. [12]. At 14 weeks post-myocardial infarction, a series of orthogonal short- and long-axis tagged MR images were acquired as described in detail previously [11]. The tags were laid down at end of diastole by synchronization to the R-wave of the electrocardiographic (ECG) signal and tagged images were captured as the heart continues through systole. LV pressure was measured with a non-ferromagnetic transducer-tipped pressure catheter (model SPC-320; Millar Instruments, Houston, TX) inserted into the LV via sterile neck incisions as described previously [11]. The end-diastolic and end-systolic LV pressures (P_ED, P_SD) recorded were used to define the endocardial boundary conditions of the FE model.

A customized version of the MR image tagging post-processing software, FindTags (Laboratory of Cardiac Energetics, National Institutes of Health, Bethesda, MD) was used to contour the endocardial and epicardial LV surfaces and also to segment the systolic tags for each image slice [13]. Systolic myocardial strains (6 Lagrangian Green’s strain tensor components in cylindrical coordinates, circumferential, longitudinal and radial) at midwall and around the circumference in each short-axis slice were calculated from tag-line deformation using the four dimensional B-spline-based motion tracking technique [14] (Fig. 1).

Fig. 1

3D cardiac strain analysis from in vivo tagged MR images. Endocardial and epicardial contours as well as segmented tag-lines were traced from (a) short and (b) long-axis MR images to create (c) a 3D geometry. (d) Each short axis slice was divided into (more ...)

FE Model

An FE model was created using early diastole as the initial unloaded reference state since the LV pressure is lowest at this point and therefore stress is at a minimum. From the LV contours at early diastole, aneurysm, remote and borderzone regions were determined based on the ventricular wall thickness. Specifically, the borderzone region is defined as the steep transition in wall thickness between remote and aneurysm regions [15]. Surface meshes were then created from the LV contours to replicate the in vivo geometry (Rapidform, INUS Technology, Inc., Sunnyvale, CA). The spaces between the endocardium and epicardium surfaces were filled with 8-noded brick elements with a single integration point for computational efficiency to generate a volumetric mesh that is refined into three elements transmurally (Truegrid, XYZ Scientific Applications, Inc., Livermore, CA). Each zone, remote, borderzone and infarct, were assigned different material properties. The inner endocardial layer was lined with a layer of shell elements that extends to LV base to form an enclosed volume for LV volume measurements. The shell elements were modeled as an extremely soft linearly elastic material (Young’s modulus of 1×10⁻¹⁰ kPa and Poisson’s ratio of 0.3) that offered no mechanical response. The reliability of the model predictions was tested with a mesh convergence study to find the minimum number of elements needed to produce accurate results within the fastest computation time. The mesh convergence study determined that 2496 elements is required and further mesh refinement only results in a 1% change in strain predictions (Fig. 2).

Fig. 2

Creation of the FE model of the LV using geometry from in vivo tagged MR images. Endocardial and epicardial contours extracted from (Fig. 1a) short and (Fig. 1b) long-axis MR images were used to generate (a) a surface mesh with three distinct LV regions (more ...)

Cardiac myofiber angles of −37°, 23° and 83° were assigned at the epicardium, midwall and endocardium, respectively, in the remote and borderzone regions [16]. At the aneurysm region, fiber angles were set to 0° in order to use experimentally determined aneurysm material parameters with respect to this direction [17]. In other words, the constitutive equation for the aneurysm is in terms of strain components referred to cardiac (i.e., circumferential and longitudinal) coordinates instead of fiber coordinates. Nodes at the LV base were restricted to displace horizontally and circumferential displacements were constrained at the basal epicardial nodes. The inner endocardium wall was loaded to the measured in vivo end-diastolic and end-systolic LV pressures. Surface meshes and subsequent volumetric meshes were also created from the LV contours at end-diastole and end-systole to provide the end-diastolic and end-systolic LV volumes.

Constitutive Model

Nearly incompressible, transversely isotropic, hyperelastic constitutive laws for passive [18] and active myocardium [19] were modeled in a user-defined material subroutine in the explicit FE solver, LS-DYNA (Livermore Software Technology Corporation, Livermore, CA). The passive myocardium mechanics is described by the strain energy function, W, that is transversely isotropic with respect to the local fiber direction,

(1)

where C, b_f, b_t and b_fs are diastolic myocardial material parameters. E₁₁ is strain in fiber direction, E₂₂ is cross-fiber in-plane strain, E₃₃ is radial strain transverse to the fiber direction, and the rest are shear strains.

Systolic contraction was modeled as the sum of the passive stress derived from the strain energy function and an active fiber directional component, T₀, which is a function of time, t, peak intracellular calcium concentration, Ca₀, sarcomere length, l, and maximum isometric tension achieved at the longest sarcomere length, T_max [19],

(2)

where S is the second Piola-Kirchoff stress tensor, p is the hydrostatic pressure introduced as the Lagrange multiplier needed to ensure incompressibility and was calculated from the bulk modulus of water, J is the Jacobian of the deformation gradient tensor, C is the right Cauchy-Green deformation tensor and the Dev is the deviatoric projection operator,

(3)

is the deviatoric contribution of the strain energy function, W (Eq. 1). The assumption of near incompressibility of the myocardium requires the decoupling of the strain energy function into dilational and deviatoric components,

(4)

where U is the volumetric contribution.

The active fiber directional stress component is defined by a time-varying elastance model, which at end-systole, is reduced to [20],

(5)

with m and b as constants, and the length-dependent calcium sensitivity, ECa₅₀, is given by,

(6)

where B is a constant, (Ca₀)_max is the maximum peak intracellular calcium concentration, l₀ is the sarcomere length at which no active tension develops and l_R is the stress-free sarcomere length. The material constants for active contraction were found to be [21]: Ca₀ = 4.35 μmol/L, (Ca₀)_max = 4.35 μmol/L, B = 4.75 μm⁻¹, l₀ = 1.58 μm, m = 1.0489 secμm⁻¹, b = −1.429 sec, and l_R was set at 1.85 μm, the sarcomere length in the unloaded configuration. Based on the biaxial stretching experiments [22] and FE analyses [4, 23], cross-fiber, in-plane stress equivalent to 40% of that along the myocardial fiber direction was added.

Since tagged MR images were acquired during systole, only systolic myocardial strains could be determined. This meant only systolic material parameters, T_max in remote (T_{max_R}) and borderzone region (T_{max_B}), could be optimized. T_max in the aneurysm region was set to zero as the coronary ligations were permanently put in place to create dyskinetic infarcts. Diastolic materials parameters, b_f, b_t and b_fs, were set to the average optimized values for sheep obtained previously from Walker et al. [4]: b_f = 49.25, b_t = 19.25, b_fs = 17.44. C in the aneurysm region (C_I) was defined as 10 times stiffer than that in the remote (C_R) [4]. C_R was determined by calibrating it such that the predicted end-diastolic LV volume matched the measured value. Initial ranges for T_{max_R} were between 0.1 and 1000.0 kPa, and between 0.1 and 500.0 kPa for T_{max_B}.

Material Parameter Optimization

The commercial FE optimization software, LS-OPT (Livermore Software Technology Corporation, Livermore, CA), uses a systematic search methodology to automatically explore the parameter space and find an optimum design [24]. In this application, the optimum design is the ideal myocardial material parameters that will produce 3D myocardial strains that come closest to matching those measured in vivo. LS-OPT is based on the successive response surface method (SRSM). It works by first selecting experimental design points which consist of systolic myocardial material parameters, T_{max_R} and T_{max_B} as components, within their specified ranges (Fig. 3a). The selection process is based on the Design of Experiments approach using the D-Optimality criterion [25]. This method uses a subset of all the possible design points as a basis to solve

Fig. 3

A graphical representation of the RSM optimization approach employed by LS-OPT. (a) For the optimization of 2 parameters, 5 experimental points including the initial guess were selected using a D-Optimal method bounded by their respective upper and lower (more ...)

where X is the coefficient matrix of the normal equations:

The vector â represents the coefficients of the basis functions (in this case linear monomials) and vector y represents the computed values (e.g. strain or volume values) at all the experimental design points.

The subset is usually selected from an ⁿ -factorial design where is chosen a priori as the number of grid points in any particular dimension and n is the number of variables. A genetic algorithm is used to solve the resulting discrete maximization problem. The proper selection of experimental points results in the most accurate response surface and therefore a faster convergence rate and greater accuracy of the optimization solution [24].

Using the selected experimental points, FE simulations were performed and the strains and LV volumes were calculated (Fig. 3b) for all points. A linear response surface was then fitted to each strain and volume by means of a least square fitting method. The approximate MSE is defined as the difference between the response surface predicted and experimental results (end-diastolic and end-systolic volumes and strains).

(7)

where n is the in vivo strain point, N is the total number of in vivo strain points, E_{i j,S} are the predicted strains at each strain point, V_ED and V_SD are the predicted end-diastolic and end-systolic LV volumes, respectively. The overbar represents experimental in vivo measurements. Strain in the radial direction, E₃₃, is excluded as it cannot be measured with sufficiency accuracy with tagged MR images [26, 27]. The goal of the optimization is to minimize the MSE.

The leap-frog dynamic trajectory method was applied towards the minimization of the MSE constructed from these response surfaces to determine the optimum parameter set [28–30] (Fig. 3c). This optimum then becomes the initial experimental point for the next iteration. The range for each parameter was adjusted by either contraction or translation or both depending on how close the current optimum is to the previous one and the degree of oscillation of the optimization solution across the iterations [24]. For the next sub-region of parameter space, the D-Optimality criterion was applied again to find a new set of experimental points for the Response Surface Method (RSM) to be applied once more [31]. By repeating these steps for each successive iteration, subregions of the parameter space were reduced and shifted until a final optimum design was found (Fig. 4). LS-OPT is also able to compute confidence intervals of the optimized material parameters in order to assess their reliability [24]. Fig. 5 summarizes the methodology.

Fig. 4

A graphical representation of the SRSM. From the initial experimental points and parameter space bounded by the ranges of the 2 parameters, T_{max_R} and T_{max_B}, the first implementation of the RSM produced an initial response surface and its optimum. This (more ...)

Fig. 5

A flowchart illustrating the process involved in determining the optimum myocardial material parameters from tagged MR images and LV pressures from cardiac catheterization.

It is important to note here that a two-step process was required to compute Green’s strain referred to the end-diastolic state. Using index notation, the Green’s strain to be compared with the experimentally measured values is defined as:

(8)

where F is the deformation gradient tensor, δ is the kronecker delta, Y refers to the end-diastolic (ED) state and x to the end-systolic (ES) state. Each run starts from a given undeformed state X, so that the deformation gradient referred to ED can be computed as:

(9)

where F_mK and F_iK are the deformation gradient components at the ED and ES phases respectively, both referred to the undeformed configuration X.

The executable ‘mri2lso’ (which reads a text file containing the measured 3D strain data in simple tabular format and converts it to an LS-OPT command file) is available at ftp://ftp.lstc.com/user/ls-opt as part of the LS-OPT distribution.

RESULTS

The pressure catheter recorded P_ED and P_SD at 12.89 mmHg and 101.59 mmHg, respectively. These recorded pressures were offset by the minimum recorded LV pressure of 2.04 mmHg to get P_ED of 10.85 mmHg and P_SD of 99.55 mmHg, which were applied to the inner endocardium layer in order to create an initial stress-free model. C_R was calibrated to be 0.95 kPa and C_I defined as 10 times stiffer at 9.5 kPa such that V_ED was accurately predicted at 123.6 mL with the measured value at 123.4 mL.

The minimum MSE of 4.17, consisting of 960 strain and 2 LV volume data points, was reached in 10 iterations. The optimized T_{max_R} and T_{max_B} for this sheep are 190.1 kPa and 60.3 kPa, respectively, with 90% confidence intervals at 14.9% and 16.9%, respectively. The precision of the optimized parameters brought on by the narrowing of the parameter space over the 10 iterations (Fig. 6). V_SD was also accurately predicted at 110.8 mL, only 4.9% higher than the measured value of 105.6 mL. The predicted systolic strains using the optimized material parameters were in generally decent agreement with the in vivo measured strains. The insertion points of the right ventricle (RV) to the LV showed the largest difference between the measured and predicted strains since the RV was not included in the model. Fig. 7 illustrates the similarities in the circumferential strains that were measured in vivo from tagged MR images and predicted from this FE model. The root mean square (RMS) error for the circumferential strain component between the 137 pairs of measured and predicted strains in the remote zone was 0.048 and in the borderzone, the RMS error was 0.070 with 55 pairs of strain points.

Fig. 6

The narrowing of the parameter space for the systolic parameters, T_{max_R} and T_{max_B}, with each iteration resulted in a precise final converged optimum.

Fig. 7

Circumferential strains predicted from the present FE model are generally in decent agreement with the values measured in vivo from tagged MR images. Slice 1 is the most basal while slice 16 is the most apical. I is the posterior right ventricular insertion, (more ...)

The optimization process was performed once using the upper bounds of the parameters as the initial guess and repeated using the lower bounds as the initial guess. Both optimizations with different starting points produced the same converged results (Table 1). Each forward simulation of the FE model takes about 10 minutes on a single 2.4 GHz processor. For the optimization of 2 material parameters, 5 experimental points were selected using a D-Optimal method for each iteration. The entire optimization process involving 10 iterations or 50 experimental points plus 10 leap-frog optimizations required about 4.5 hours on four 2.4 GHz processors.

Table 1

Initial and converged systolic material parameters.

DISCUSSION

A very efficient or fast method was developed in order to formally optimize regional myocardial contractility from tagged MR images and cardiac catheterization pressures. Our approach was demonstrated for data from a single sheep, 14 weeks after anteroapical myocardial infarction. The proposed method involves performing FE simulations using the customized commercial FE solver (LS-DYNA) that was programmed with the passive and active myocardial material laws. The forward FE solutions are fed into the optimization software (LS-OPT), which was customized to determine the systolic myocardial material parameters using the SRSM approach by targeting the in vivo systolic strains and LV volumes. The in vivo systolic strains and LV volumes were determined from tagged MR images, which also provided the LV endocardial and epicardial contours that were used to generate the FE model. Finite element model loading conditions were obtained from cardiac catheterization measurements of LV pressures.

Validation of the Method

In order to really validate our method for estimating regional myocardial contractility in vivo, T_{max_R} and T_{max_B} need to be estimated independently using another method (preferably one in which forces can be directly measured). We have had great difficulty in making comparisons of T_{max_R} and T_{max_B} obtained in the same hearts with our biaxial stretcher because LV regions that are thick enough for tagged MRI are too thick for biaxial testing and LV regions that are thin enough for biaxial testing are too thin for tagged MRI. This forced us to seek another method which involves measurement of active stress developed in “skinned” (de-membraned using Triton-X) ovine LV myofibers over a range of calcium concentrations in the muscle bath. Direct measurements of maximal force (normalized by cross-sectional area) developed by skinned myofibers dissected from remote regions of three sheep hearts with LV aneurysm (38.1 ± 2.3 kPa) were significantly greater than that developed by skinned myofibers dissected from borderzone regions (22.2 ± 3.0 kPa; p < 0.03). Part of the reason why both of the optimized T_max values reported in the Results (T_{max_R} = 190.1 kPa; T_{max_B} = 60.3 kPa) are so much greater than the skinned fiber values is the latter are fit to a much simpler (modified Hill) equation than Eq. 5 in the Methods. Specifically, the modified Hill equation does not include the factor of 0.5 in front of T_max, nor the factor in the parentheses, the value of which depends on sarcomere length. This sarcomere length-dependence can be removed from Eq. 5 by setting the parameter m = 0 and b ≥ 2000. Rerunning the optimization with the modified Hill equation in place of Eq. 5 results in T_{max_R} = 100.0 kPa and T_{max_B} = 36.3 kPa. It is not difficult to explain the higher in vivo contractility than ex vivo contractility in terms of tissue damage caused by the dissection process. Nevertheless, the significant depression in optimized T_{max_B} relative to T_{max_R} was confirmed by direct ex vivo force measurements from skinned fiber preparations.

Sensitivity Analysis

Effect of passive and active parameters specified on Tmax

To estimate how sensitive the model formal optimization is to the specified passive material parameter values, simulations were rerun after increasing or decreasing by 10% (one-at-a-time) the values specified in the Methods for b_f, b_t, and b_fs. The greatest sensitivity of the optimized T_{max_R} or T_{max_B} to these three passive material parameters was only a 5.6% increase in T_{max_R} in response to a 10% decrease in b_t and a 2.1% decrease in T_{max_B} in response to a 10% decrease in b_fs. The model formal optimization was much more sensitive to some of the specified active material parameter values, especially l₀ and l_R. Specifically, a 5% increase in l₀ resulted in a 32% increase in T_{max_R} and 5% decrease in l_R resulted in a 195% increase in T_{max_R}. However, we are quite confident that the values specified in the Methods for l₀ and l_R are very well defined by experimental data in the literature. The active material parameter Ca₀ also “controls” model myocardial contractility (like T_max) but instead of simply scaling active stress development, it affects the shape of the active stress versus sarcomere relationship (i.e., linear at Ca₀ = 2 μmol/L, “concave up” at Ca₀ < 2 μmol/L; concave down at Ca₀ > 2 μmol/L; see Fig. 8 of Guccione and McCulloch 1993). Not surprisingly, 10% and 20% decreases in Ca₀ resulted in increases of 7.7% and 18.9% in T_{max_R} and increases of 4.6% and 10.5% in T_{max_B}, respectively. Unfortunately, to the best of our knowledge, there is no way to reliably measure peak intracellular calcium concentration in the intact beating LV myocardium.

Effect of errors in strain measurement on Tmax

To estimate how sensitive the model formal optimization is to errors in strain measurement, simulations were rerun after adding or subtracting fixed offsets of 0.015 and 0.03 in circumferential strain (E_cc) measurements. The greatest sensitivity of the optimized T_{max_R} or T_{max_B} to uncertainty in circumferential strain measurement was a 112% increase in T_{max_R} in response to a fixed decrease of 0.03 in E_cc and a 33% increase in T_{max_B} in response to fixed increase of 0.03 in E_cc. Not surprisingly, these regional “indices” of myocardial contractility are very sensitive to the degree of systolic LV circumferential shortening.

Relative Contributions of Terms in Objective Function (Eq. 7)

We performed a detailed study of the effect of the relative contributions of strains and volume on T_{max_R} and T_{max_B} by varying the weights over a range of two orders of magnitude (LS-OPT could not handle weights below 1%). Under all (weighting) conditions T_{max_R} = 190 kPa and T_{max_B} = 60 kPa. The only “measurable” effect was in the 3^rd significant figure of T_{max_B} (60.3 kPa for equal weighting versus 60.4 kPa for weighting of 98% for strain and 1% for each of the two volumes). Moreover, the confidence intervals were equally “tight” for these two (extreme weighting) cases.

Comparison to Previous Work

Previous attempts at estimating myocardial material parameters from an intact heart using strains have been oversimplified by assigning non-physiological material properties [32] and loading conditions [33]. Some studies with more realistic FE models however only used 2 strain measurements to fit for the optimum material parameters [34-34], while others that used 500 to 700 strain measurements for optimization required a significant amount of computing time and power [4, 36]. There have also been labor intensive, trial-and-error approaches towards optimization [4, 34] as well as automated optimization algorithms. The elegant study by Augenstein et al. [37] is an example of a nonlinear automated optimization procedure, which was independently validated using phantoms.

The optimization algorithms used by previous studies [32–33, 36, 38] were numerical gradient-based, which may lead to local minima rather than the global one due to local sensitivities of the nonlinear problem and require a good initial guess to the solution. This study employed the gradient-free optimization approach called the successive response surface method (SRSM), which is better able to reach the global minimum [39, 40]. It is based on the iterative construction and minimization of response surfaces that have been fitted onto the residuals between the measured and FE predicted strains. The approximations avoid local minima and the shrinkage of the parameter space with each iterative response surface ensures the detailed variations of the nonlinear problem near the global minimum are captured. The optimized myocardial contractility in the remote region was found to be 190.1 kPa while in the borderzone, the contractility was 3.15 times less at 60.3 kPa. The marked difference in contractility in the two regions is greater than the 50% reduction reported by Guccione et al. [34]. The reason for the discrepancy can be attributed to the greater accuracy and precision of the present FE model and optimization method. 960 in vivo systolic strain measurements (12 circumferentially, 16 longitudinally, 5 cylindrical strain components) and 2 LV volumes formed the objectives of optimization in this study while the other study only used 2 strain measurements (only circumferential strains in just the anterior and posterior borderzone regions). The more data points available for optimization, the more precise the optimized results will be. This relationship is reflected by the relatively small 90% confidence intervals for the optimized T_{max_R} and T_{max_B}, 14.9% and 16.9%, respectively. Guccione et al. [34] did not report any analysis on the confidence of the optimized results.

Advantages of SRSM

The SRSM optimization algorithm is a gradient-free approach that is suitable for highly nonlinear problems with multiple local minima. SRSM is based on the sequential optimization of the linear response surfaces constructed from the FE simulation results. The size of the successive sub-regions, within which the response surfaces are formed, is adapted based on contraction and panning parameters designed to prevent oscillations and premature convergence [41]. As SRSM starts with a large sub-region of the parameter space that is contracted automatically as the optimum is approached, local sensitivities do not cause large departures from the previous design. This inherent adjustment limit of the parameter space in the algorithm eliminates the step-size dilemma in gradient-based optimization methods [41]. In previous myocardial material optimization studies that used numerical gradients [32–33, 36, 38], finite differences were taken at 5 to 10% interval for each material parameter. If the interval was set too large, accuracy is lost. Conversely, if the interval was set too small, the gradient may be spurious and miss the true optimum.

Compared to another gradient-free optimization approach, like the genetic algorithm used by Nair et al. [35], SRSM is more efficient. Even though that previous study only used 2 strain measurements for optimization, 5400 FE simulations and 150 “chromosomal generations” or iterations were required for convergence and the total computation time took 25 to 40 days. Although this study optimized 2 fewer material parameters, a much larger number of data points, 960 strain measurements and 2 LV volumes, were used for optimization and the total computation time was only 4.5 hrs. The genetic method takes longer to compute probably because it adds randomness into each material parameters selection process for FE analysis in order to widen the parameter space to find the global minimum but doing so also slowed the rate of convergence.

Limitations

There are a few limitations in this study. First, since tagged MR images were acquired during the systolic phase, only systolic strains were calculated and therefore just the systolic myocardial material parameters can be optimized. Due to the lack of diastolic strains, the diastolic myocardial material parameters, b_f, b_t, b_fs, were taken from Walker et al. [4], which had optimized those parameters for four sheep. The other diastolic material parameters, C_R and C_I, had to be calibrated to the measured end-diastolic LV volume. For future studies, 2 sets of tagged MR images acquired during diastole and systole will be required. Strain analysis will be performed on both image sets to produce diastolic and systolic strains, which will then be used separately to optimize for diastolic and systolic myocardial material parameters, respectively.

The second limitation of this study is the assumption of an initial stress-free state, which is present in all previous FE simulations of the heart based on in vivo images. Not all the blood in the heart is ejected out during systole and it is this residual blood at early diastole that results in a small amount of stress on the heart. Our previous FE model study suggests that residual stress produced by surgical remodeling has little effect on ventricular function and regional mechanics [42].

The next limitation is with the use of transversely isotropic material properties to describe the myocardium, which is consistent with biaxial tests that observed stiffness in the muscle fiber direction is greater than in cross-fiber direction. However, histological sections of myocardium show myofibers arranged into branching laminae, implying stiffness is lower normal to the laminar plane than within it [43, 44]. Direct measurements of material orthotropy have yet been determined due to experimental limitations of triaxial testing of excised myocardium. Orthotropic material properties will be modeled in future studies and optimized diastolic material parameters will be obtained to determine the degree of anisotropy.

The last limitation is the absence of the RV in the FE model. Only the LV was modeled since it experiences the largest stress and strain as it has to eject blood from the heart with enough pressure to reach the rest of the body. Although the RV just needs to contract with sufficient force for the blood to reach the nearby lungs, the deformation of the RV seems to have a significant effect on the LV wall strains at the RV insertion points. This effect is evident by the largest deviations between the measured and predicted strains right at the RV insertion points. Future studies will need to include the RV in the model and strain analysis.

CONCLUSIONS

In summary, we have presented a method to formally optimize regional myocardial contractility in vivo that is at least an order of magnitude faster than that used in our previous studies [4–5]. The significant depression in optimized T_{max_B} relative to T_{max_R} was confirmed by direct ex vivo force measurements from skinned fiber preparations. The optimized values of T_{max_B} and T_{max_R} were not overly sensitive to the passive material parameters specified. The computation time of less than 5 hours associated with our proposed method for estimating regional myocardial contractility in vivo makes it a potentially very useful clinical tool.

Acknowledgments

The authors acknowledge financial support from NIH grants R01-HL-77921, R01-HL-63348 and R01-HL-84431.

References

1. Guccione JM, Walker JC, Beitler JR, Moonly SM, Zhang P, Guttman MA, Ozturk C, McVeigh ER, Wallace AW, Saloner DA, Ratcliffe MB. The effect of anteroapical aneurysm plication on end-systolic three-dimensional strain in the sheep: a magnetic resonance imaging tagging study. J Thorac Cardiovasc Surg. 2006;131(3):579–586 e3. [PMC free article] [PubMed]

2. Walker JC, Guccione JM, Jiang Y, Zhang P, Wallace AW, Hsu EW, Ratcliffe MB. Helical myofiber orientation after myocardial infarction and left ventricular surgical restoration in sheep. J Thorac Cardiovasc Surg. 2005;129(2):382–90. [PubMed]

3. Costa KD, Hunter PJ, Wayne JS, Waldman LK, Guccione JM, McCulloch AD. A three-dimensional finite element method for large elastic deformations of ventricular myocardium: II--Prolate spheroidal coordinates. J Biomech Eng. 1996;118(4):464–72. [PubMed]

4. Walker JC, Ratcliffe MB, Zhang P, Wallace AW, Fata B, Hsu EW, Saloner D, Guccione JM. MRI-based finite-element analysis of left ventricular aneurysm. Am J Physiol Heart Circ Physiol. 2005;289(2):H692–700. [PubMed]

5. Walker JC, Ratcliffe MB, Zhang P, Wallace AW, Hsu EW, Saloner DA, Guccione JM. Magnetic resonance imaging-based finite element stress analysis after linear repair of left ventricular aneurysm. J Thorac Cardiovasc Surg. 2008;135(5):1094–102. 1102, e1–2. [PubMed]

6. Jan KM. Distribution of myocardial stress and its influence on coronary blood flow. J Biomech Eng. 1985;18(11):815–20. [PubMed]

7. Sarnoff SJ, Braunwald E, Welch GH, Jr, Case RB, Stainsby WN, Macruz R. Hemodynamic determinants of oxygen consumption of the heart with special reference to the tension-time index. Am J Physiol. 1958;192(1):148–56. [PubMed]

8. Yin FC. Ventricular wall stress. Circ Res. 1981;49(4):829–42. [PubMed]

9. Grossman W. Cardiac hypertrophy: useful adaptation or pathologic process? Am J Med. 1980;69(4):576–84. [PubMed]

10. Huisman RM, Elzinga G, Westerhof N, Sipkema P. Measurement of left ventricular wall stress. Cardiovasc Res. 1980;14(3):142–53. [PubMed]

11. Zhang P, Guccione JM, Nicholas SI, Walker JC, Crawford PC, Shamal A, Acevedo-Bolton G, Guttman MA, Ozturk C, McVeigh ER, Saloner DA, Wallace AW, Ratcliffe MB. Endoventricular patch plasty for dyskinetic anteroapical left ventricular aneurysm increases systolic circumferential shortening in sheep. J Thorac Cardiovasc Surg. 2007;134(4):1017–24. [PMC free article] [PubMed]

12. Markovitz LJ, Savage EB, Ratcliffe MB, Bavaria JE, Kreiner G, Iozzo RV, Hargrove WC, 3rd, Bogen DK, Edmunds LH., Jr Large animal model of left ventricular aneurysm. Ann Thorac Surg. 1989;48(6):838–45. [PubMed]

13. Guttman MA, Zerhouni EA, McVeigh ER. Analysis and visualization of cardiac function from MR images. IEEE Comp Graph Appl. 1997;17(1):30–8. [PMC free article] [PubMed]

14. Ozturk C, McVeigh ER. Four-dimensional B-spline based motion analysis of tagged MR images: introduction and in vivo validation. Phys Med Biol. 2000;45(6):1683–702. [PMC free article] [PubMed]

15. Moustakidis P, Maniar HS, Cupps BP, Absi T, Zheng J, Guccione JM, Sundt TM, Pasque MK. Altered left ventricular geometry changes the border zone temporal distribution of stress in an experimental model of left ventricular aneurysm: a finite element model study. Circulation. 2002;106(12 Suppl 1):I168–75. [PubMed]

16. Omens JH, May KD, McCulloch AD. Transmural distribution of three-dimensional strain in the isolated arrested canine left ventricle. Am J Physiol. 1991;261(3 Pt 2):H918–28. [PubMed]

17. Moonly S. Experimental and computational analysis of left ventricular aneurysm mechanics. University of California, San Francisco with University of California; Berkeley, San Francisco, CA: 2003.

18. Guccione JM, McCulloch AD, Waldman LK. Passive material properties of intact ventricular myocardium determined from a cylindrical model. J Biomech Eng. 1991;113(1):42–55. [PubMed]

19. Guccione JM, Waldman LK, McCulloch AD. Mechanics of active contraction in cardiac muscle: Part II--Cylindrical models of the systolic left ventricle. J Biomech Eng. 1993;115(1):82–90. [PubMed]

20. Tozeren A. Continuum rheology of muscle contraction and its application to cardiac contractility. Biophys J. 1985;47(3):303–9. [PubMed]

21. Guccione JM, Costa KD, McCulloch AD. Finite element stress analysis of left ventricular mechanics in the beating dog heart. J Biomech. 1995;28(10):1167–77. [PubMed]

22. Lin DH, Yin FC. A multiaxial constitutive law for mammalian left ventricular myocardium in steady-state barium contracture or tetanus. J Biomech Eng. 1998;120(4):504–17. [PubMed]

23. Usyk TP, Mazhari R, McCulloch AD. Effect of laminar orthotropic myofiber architecture on regional stress and strain in the canine left ventricle. J Elasticity. 2000;61:143–64.

24. Stander N, Roux W, Eggleston T, Craig K. LS-OPT User’s Manual Version 3.3 2007

25. Myers RH, Montgomery DC. Response Surface Methodology - Process and Product Optimization Using Design Experiments. Wiley; New York: 1995.

26. Declerck J, Denney TS, Ozturk C, O’Dell W, McVeigh ER. Left ventricular motion reconstruction from planar tagged MR images: a comparison. Phys Med Biol. 2000;45(6):1611–32. [PMC free article] [PubMed]

27. Denney TS, Jr, Gerber BL, Yan L. Unsupervised reconstruction of a three-dimensional left ventricular strain from parallel tagged cardiac images. Magn Reson Med. 2003;49(4):743–54. [PubMed]

28. Snyman JA. A new and dynamic method for unconstrained minimization. Appl Math Modeling. 1982;6:449–462.

29. Snyman JA. An improved version of the original leap-frog dynamic method for unconstrained minimization. Appl Math Modeling. 1994;8:453–460.

30. Snyman JA. The LFOPC leap-frog algorithm for constrained optimization. Computers Math Appl. 2000;40:1085–96.

31. Roux W, Stander N, Haftka RT. Response surface approximations for structural optimization. Intl J Num Methods Eng. 1998;42:517–34.

32. Moulton MJ, Creswell LL, Downing SW, Actis RL, Szabo BA, Pasque MK. Myocardial material property determination in the in vivo heart using magnetic resonance imaging. Int J Card Imaging. 1996;12(3):153–67. [PubMed]

33. Okamoto RJ, Moulton MJ, Peterson SJ, Li D, Pasque MK, Guccione JM. Epicardial suction: a new approach to mechanical testing of the passive ventricular wall. J Biomech Eng. 2000;122(5):479–87. [PubMed]

34. Guccione JM, Moonly SM, Moustakidis P, Costa KD, Moulton MJ, Ratcliffe MB, Pasque MK. Mechanism underlying mechanical dysfunction in the border zone of left ventricular aneurysm: a finite element model study. Ann Thorac Surg. 2001;71(2):654–62. [PubMed]

35. Nair AU, Taggart DG, Vetter FJ. Optimizing cardiac material parameters with a genetic algorithm. J Biomech. 2007;40(7):1646–50. [PubMed]

36. Remme EW, Hunter PJ, Smiseth O, Stevens C, Rabben SI, Skulstad H, Angelsen BB. Development of an in vivo method for determining material properties of passive myocardium. J Biomech. 2004;37(5):669–78. [PubMed]

37. Augenstein KF, Cowan BR, LeGrice IJ, Nielsen PM, Young AA. Method and apparatus for soft tissue material parameter estimation using tissue tagged Magnetic Resonance Imaging. J Biomech Eng. 2005;127(1):148–57. [PubMed]

38. Augenstein KF, Cowan BR, LeGrice IJ, Young AA. Estimation of cardiac hyperelastic material properties from MRI tissue tagging and diffusion tensor imaging. Med Image Comput Comput Assist Interv Int Conf Med Image Comput Comput Assist Interv. 2006;9(Pt 1):628–35. [PubMed]

39. Einstein DR, Freed AD, Stander N, Fata B, Vesely I. Inverse parameter fitting of biological tissue: a response surface approach. Ann Biomed Eng. 2005;33(12):1819–1830. [PubMed]

40. Chen K, Fata B, Einstein DR. Characterization of the highly nonlinear and anisotropic vascular tissues from experimental inflation data: a validation study toward the use of clinical data for in-vivo modeling and analysis. Ann Biomed Eng. 2008;36(10):1668–80. [PMC free article] [PubMed]

41. Stander N, Craig K. On the robustness of a simple domain reduction scheme for simulation-based optimization. Eng Comput. 2002;19(4):431–50.

42. Guccione JM, Moonly SM, Wallace AW, Ratcliffe MB. Residual stress produced by ventricular volume reduction surgery has little effect on ventricular function and mechanics: a finite element model study. J Thorac Cardiovasc Surg. 2001;122(3):592–9. [PubMed]

43. LeGrice IJ, Smaill BH, Chai LZ, Edgar SG, Gavin JB, Hunter PJ. Laminar structure of the heart: ventricular myocyte arrangement and connective tissue architecture in the dog. Am J Physiol. 1995;269(2 Pt 2):H571–82. [PubMed]

44. LeGrice IJ, Hunter PJ, Young AA, Smaill BH. The architecture of the heart: a data-based model. Phil Trans R Soc Lond. 2001;A 359:1217–32.