Anisotropic Elastography for Local Passive Properties and Active Contractility of Myocardium from Dynamic Heart Imaging Sequence

Major heart diseases such as ischemia and hypertrophic myocardiopathy are accompanied with significant changes in the passive mechanical properties and active contractility of myocardium. Identification of these changes helps diagnose heart diseases, monitor therapy, and design surgery. A dynamic cardiac elastography (DCE) framework is developed to assess the anisotropic viscoelastic passive properties and active contractility of myocardial tissues, based on the chamber pressure and dynamic displacement measured with cardiac imaging techniques. A dynamic adjoint method is derived to enhance the numerical efficiency and stability of DCE. Model-based simulations are conducted using a numerical left ventricle (LV) phantom with an ischemic region. The passive material parameters of normal and ischemic tissues are identified during LV rapid/reduced filling and artery contraction, and those of active contractility are quantified during isovolumetric contraction and rapid/reduced ejection. It is found that quasistatic simplification in the previous cardiac elastography studies may yield inaccurate material parameters.


INTRODUCTION
Major heart diseases including ischemia, hypertrophic myocardiopathy, and myocardial dilatation are growingly pervasive with high morbidity and mortality at tremendous social and healthcare costs [1]. It is commonly observed that heart diseases are accompanied with impaired myocardium, which typically shows increase in the passive viscoelastic stiffness and/or decrease in active contractile ability. For example, in ischemia, there is loss of blood supply with the subsequent development of localized or diffuse fibrosis. The ischemic muscle thus becomes stiff with passive properties [2][3][4], and its active contractility is commonly weakened [5]. In the dilatation/hypokinesis states, the myocardial contractility is weakened [6]. Such impairments reduce the heart pumping ability at normal ventricular wall stress, and require increase in stress to maintain the need of circulation. Chronic high myocardial stress may lead to further damages, and may eventually cause complete heart failure. Thus, modality development for detecting, localizing, and evaluating the impaired myocardium should have far-reaching impacts on diagnosis of heart diseases, monitoring therapy, and design of surgery.
Cardiac biomechanics [7][8][9][10][11][12][13][14] has been developed to simulate anatomically accurate heart model undergoing dynamic deformation, taking into consideration the sheet-like myofiber architecture, finite-strain deformation, nonlinearity and passive/active behaviors of myocardium. With prescribed material properties of the myocardial tissue, cardiac biomechanics calculates myocardial motion and intramural strain/stress in response to ventricular blood pressure and bioelectric stimulus. Importance of the wall stress has been well recognized, as it influences major physiological factors including the ventricular pumping performance, the energy consumption in the myocardial tissue, the coronary blood flow, the vulnerability of regions to ischemia and infarction, and the remodeling of myocardium and risk of arrhythmias [15][16][17][18]. However, stress can only be calculated from the strain with knowledge of the material properties. Inaccurate constitutive description of myocardium will necessarily lead International Journal of Biomedical Imaging to inaccurate wall stress, and may mislead the clinical evaluations.
There have been in vitro experiments to characterize the passive elastic properties of heart muscle [19][20][21][22][23][24][25]. While these experiments provide important knowledge of the mechanical behaviors of myocardium, the main disadvantages are that the cut specimens are not in intact condition, and that the varieties due to the diversity of location, species, sex, age and health condition of the heart, and so forth, cannot be fully addressed. Furthermore, such tests are inappropriate in clinics. Therefore, the fundamental importance of accurate constitutive description of an individual living heart urges a cost-efficient, in vivo and noninvasive modality that reveals the passive viscoelastic properties and active contractility of normal and abnormal myocardium.
Elastography is a biomechanical technique that quantifies the material properties of an object based on in vivo measurements of deformation and force, instead of in vitro tests with cut specimens. A comprehensive elastography framework for heart tissues should include the following key elements: material anisotropy due to the sheet-like myofiber architecture, nonlinearity and viscosity of the passive mechanical properties, bioelectricity-stimulated active contraction stress, inhomogeneity and dynamic deformation, and so forth. The pioneering cardiac elastography studies include [26][27][28][29][30][31]. Han et al. [28] and Gotteiner et al. [31] conducted a series of cardiac elastography simulations, in which the heart muscle is modeled as a uniform (homogeneous) and isotropic nonlinear material, and the heart beating is approached as being quasistatic. The 2D simulations of Creswell et al. [29] and Moulton et al. [30] further took into consideration the myocardial anisotropy. In these works, the effects of time-dependent deformation have yet to be investigated, and no attempt has been made to characterize the active contraction capability of normal and impaired heart muscles. Therefore, the resulting myocardial material parameters may not be adequate and sufficiently accurate for clinical applications.
There have been tremendous developments in highspeed, high-resolution biomedical imaging technologies. The 3D anatomic structures of a living heart can now be accurately depicted [32][33][34][35]. The laminar orthotropic myofiber architecture can be reconstructed from diffusion-tensor MRI data [36][37][38]. The dynamic deformation has also become measurable, for instance, Lin and Robb [34] and Eusemann et al. [35] used a CT-based dynamic spatial reconstructor to track the endocardial and epicardial displacement, Kanai et al. [39] developed ultrasonic phased-tracking method to obtain the velocity field in the heart wall. Therefore, it is time to develop a dynamic cardiac elastography methodology that takes advantages of these information-rich measurements and gives accurate material parameters for both the passive behaviors and active contractibility of normal and impaired myocardium.
The proposed DCE is to identify the material parameters of myocardium, with which the computed time-dependent displacement matches optimally to the experimentally measured displacement, as described in Sections 2.1 and 2.2. In Section 2.3, a dynamic adjoint method is derived for cost-efficient computation of the gradients of the objective function with respect to the material parameters, and is compared to the previous direct method. The overall DCE framework is described in Section 2.4. Simulations are presented in Section 3 using a numerical left ventricle (LV) phantom containing an ischemic region, whose passive stiffness is much higher than the normal tissue and the active contractility is much weaker. An orthotropic viscoelastic model is employed for the passive behavior, according to the myofiber architecture. A simplified active contraction stress model is applied. The displacement on endocardial and epicardial LV surfaces is extracted from the forward simulation at discrete time steps, and is used as measurement for DCE reconstruction. Material parameters for the passive viscoelastic behaviors are identified with the LV passive (filling) phase, as in Section 3.5, and those for the active contraction are identified with the LV active (contraction) phase, as in Section 3.6. In Section 4, simulations are conducted to investigate the effects of measurement errors with the myofiber architecture and the displacement, respectively. To demonstrate the necessity of using dynamic deformation for cardiac elastography, reconstruction assuming quasistatic deformation is conducted, and is compared with the present results. Conclusion is given in Section 5.

DYNAMIC CARDIAC ELASTOGRAPHY FRAMEWORK
In this section, an optimization-based algorithm is derived for DCE to identify the passive properties and active contractility of heart muscle. The myocardium is assumed undergoing viscoelastic deformation, involving active contraction stress. An objective function is proposed to measure the discrepancy between measured and computed displacements in a cardiac cycle. A dynamic adjoint method is developed for analytical gradients of the objective function with respect to the to-be-identified material parameters, for the purpose of improving the numerical efficiency for elastography reconstruction. A flowchart is given of the optimization-based DCE reconstruction procedure.

Material model and motion equation
A biological object Ω 0 is described as undergoing viscoelastic deformation, with active contraction stress. Its boundary consists of Γ 0 T , where surface force is applied, and Γ 0 u on which displacement and velocity are prescribed. The displacement at location x and time t is u(x, t), the corresponding experimental measurement is U m (x, t). The displacement-strain relation is assumed linear, that is, the strain is ε = (Öu + (Öu) T )/2. We use "¡" to denote time derivatives, for example, strain rate isε and acceleration isü. The stress σ(x, t) of myocardium is calculated with a Voigt-type viscoelastic constitutive relation: Yi Liu et al.

3
The first term describes passive elastic behavior with a strain energy W that depends on material parameters p p . The second term is a viscous stress proportional to the strain ratė ε by a forth-order tensor C v (p v ), where p v denotes viscosity parameters. The third term is the active contraction stress σ f due to the calcium concentration Ca 2+ . It is described by active parameters p a , and depends on the local strain ε. The mechanical properties for myocardium is very complex, and there may not be an accurate yet neat constitutive model [26][27][28][29][30][31]. In general, the viscous stress depends on the whole deformation history. For a heart, however, the effects of history have been "shaken down" after many beating cycles, and can be equivalently described with the current strain ε and strain rateε, and so forth. Thus the Voigt-type viscoelastic model (1) is considered a first-order approximation. The proposed DEC itself can adapt to more complex material model. Assuming no body force, the viscoelastic motion equation and the corresponding boundary and initial conditions are where ρ 0 is the mass density, T is the surface force and n is the unit outer normal to the boundary Γ 0 T . For a heart, T is well approximated as zero on the epicardial surfaces and is P(x, t)n on the endocardial surfaces.
Here P(x, t) denotes the blood pressure. It changes with time and varies in the four heart chambers. The dynamic displacement u and its time derivativeu are continuous in Ω 0 and prescribed on Γ 0 u . The variational weak form of (2) can be shown as

Objective function
The elastography reconstruction is thus to find material parameters, denoted with p = p p , p v , p a , with which the displacement u calculated via (1) and (2) approaches experimental measurement U m optimally. This is to minimize the objective function in which the time interval [0, T 0 ] is typically a part of a cardiac cycle. Second-order tensors χ Ω (x, t) and χ Γ (x, t) are weight functions defined in Ω 0 and on Γ 0 T , respectively. In medical applications, displacement U m is typically measured at M discrete time steps t = t I (I = 1, 2, . . . , M). Therefore, χ Ω (x, t) and χ Γ (x, t) should be interpreted as δ-functions with respect to time, that is, In fact, these δ-functions lead to impulsive force on the adjoint system described below.

Dynamic adjoint method for the gradients
The simplest optimization methods only require usersupplied objective function Φ(p). The material parameters p are searched with trial-and-error method, finite-deference gradients [30], or approximate gradients [40]. Without usersupplied gradients, the optimization-based reconstruction procedure is very time-consuming, considering the large number of iterative steps needed for convergent results and the high computational expense for solving time-dependent equation (2) in each step. The difficulty for calculating the analytic gradients ∂Φ/∂p comes from the fact that the dependency of u on p is implicit and highly nonlinear. Here, we derive a dynamic adjoint method for efficient and accurate computation of ∂Φ/∂p.
Following the standard treatment for constrained optimization, the weak form (3) is introduced into the objective function (4), and a Lagrangian is formed as where u ¾ κ and w ¾ κ 0 . Noting that the weak form (3) is satisfied by u with any w ¾ κ 0 , it can be shown that Φ(p) = L(p; u, w) and δΦ = δL for arbitrary displacement fields w, δw, and δu that belong to κ 0 , where "δ" represents variation. This allows suitable choice for w, with which only simple terms are kept in δL and δΦ. To find the optimal choice, variational derivations are conducted. First, the variation of L(p; u, w) is expanded as for which the terms related to δw have been canceled out since δw ¾ κ 0 and u satisfies the weak form (3). The forthorder tensor L T = ∂ 2 W/∂ε∂ε + ∂σ f /∂ε is the tangent elastic modulus. The stress variation due to the change of material parameters p is denoted as δ p σ. By integrating by part, the 4 International Journal of Biomedical Imaging terms with δü and δε are simplified as noting that δu( (7) can be further simplified, and (6) becomes Compared to the variational weak form (3), it is readily shown that the following choice for w will eliminate the first two integrals in (8) for arbitrary δu ¾ κ 0 , that is, w ¾ κ 0 and Ö ¡ L T : Note that (9) gives end conditions for w, instead of initial conditions. Thus, only the last term remains in (8), that is, The above variational derivations reflect the original idea of the continuum-based optimality criteria technique for topology optimization [41], where w defined in (9) is called the adjoint displacement. Recent elastography implementations of the adjoint method have been made by Tardieu and Constantinescu [42] for determination of elastic coefficients from indentation experiments, and by Oberai et al. [43] and Liu et al. [44] for static isotropic and anisotropic problems, respectively.
Alternatively, the gradients ∂Φ/∂p can be calculated with a direct method. Denoting the gradient of u with respect to the ith material parameter as u (i) = ∂u/∂p i , and the associated strain as ε (i) = (Öu (i) + (Öu (i) ) T )/2, it is shown that where Compared to the direct method, the most favorable advantage of the adjoint method is that it calculates the gradients ∂Φ/∂p most efficiently. With suitably chosen time integral algorithm, for instance, the standard Newmark-β scheme employed in our simulations, solutions for w and u share the same finite-element effective stiffness matrix K eff = [β(Δt) 2 K T +γΔtC+M] (Δt is time step, β = 1/12 and γ = 1/2 are used) and its Cholesky factorization at any time step, which are most time consuming in finite-element computations. Note that matrix K T corresponds to the tangent elastic modulus L T , C corresponds to viscosity tensor C v , and M is the mass matrix. Therefore, solving (9) for w and then integrating (10) for ∂Φ/∂p just add a small fraction of cost to the solution for u (2), which is required in optimization-based elastography reconstruction. In contrast, the direct method (11) and (12) involves solving dynamic displacement u (i) for every ∂Φ/∂p i . When the number of material parameters increases, the computational expense with the direct gradient method increases proportionally, while that with the proposed adjoint method rises slightly only because more integrals are calculated (10). The tradeoff for using adjoint method is additional complexity in programming. Because the solution for w (9) starts at T 0 and ends at 0, it needs a backward numerical scheme, which must be consistent with the forward scheme for u. The method also needs additional memory to store the factorization of K eff at every time step.
Recently, Oberai et al. [45] gave an in-detail evaluation on the numerical efficiency of static elastography reconstructions using the adjoint method.

Optimization-based DCE reconstruction procedure
The DCE framework involves data acquisition, material modeling, and reconstruction of material parameters, shown Is Φ(p) small enough? with the flowchart in Figure 1. Data acquisition includes reconstruction of the anatomic structures of a heart and the myofiber architecture, extracting displacement data from dynamic cardiac imaging sequence, and measurement of the blood pressures and calcium kinetics. In this work, data acquisition is considered well ready (e.g., [33][34][35][36][37][38]). Formulas for the passive and active stresses are chosen, with material parameters p = p p , p v , p a to be reconstructed. The reconstruction procedure is optimization-based, making use of a limited-memory BFGS (L-BFGS) optimization subroutine [46], for which user-supplied gradients are required. When solving u with finite-element method, the factorization of K eff is stored at each time step, and is used for the adjoint displacement w.

SIMULATIONS
This section presents DCE simulations with a threedimensional left ventricle (LV) phantom containing an ischemic region, demonstrating the validity and utility of the proposed algorithm. The passive stress-strain relation is assumed linear viscoelastic, and is orthotropic according to the myofiber architecture. A simplified model is employed for the active stress. Forward computation is first carried on the phantom undergoing filling and contraction motions. Displacement on the endocardial and epicardial LV surfaces is then extracted from the forward results at discrete time steps, and is used as measurement for DCE reconstruction of the passive and active parameters of the normal and ischemic tissues.

LV phantom and myofiber architecture
Our DCE simulations use a 3D LV phantom, as shown in Figures 2(a) and 2(b). To approximate the shape of a left ventricle (e.g., [47]), a thick-wall hollow ellipsoid is employed, whose surfaces are defined by where λ = λ epi = 1 for the epicardial surface and λ = λ endo = 0.6875 for the endocardial surface. The length unit is cm. The base is z = z Base = 2.5 cm, and the apex is z = z Apex = 5.0 cm. Shown with Figure 2(a), a red region is embedded into the wall, simulating an ischemic region. Clinically, the location, size, and shape of an impaired region can be 3D depicted from medical imaging measurements. The myocardial tissue exhibits laminar histological structure [36][37][38], and is described with the local myofiber structure-based material axes (n f , n s , n n ) as schematically shown in Figures 2(b) and 2(c), where n f is the fiber direction, n s is the cross-fiber direction in the plane of myofiber layer, and n n = n f ¢ n s . These material axes vary continuously from epicardium to endocardium and from base to apex. The present simulations use a mathematical approximation to the myocardial architecture. At a point (x, y, z) in the LV wall, local coordinates (e h , e η , e ξ ) are defined (Figure 2(b)) as: e h is along the local thickness direction (x/a 2 , y/b 2 , z/c 2 ); e η is along direction ( y/b 2 , x/a 2 , 0), that is, it lies in the x y plane and is perpendicular to e h ; the third direction is e ξ = e h ¢ e η . In consistency with the descriptions of LeGrice et al. [36] and Costa et al. [48], and as schematically shown in Figure 2(c), the myocardial fiber is well approximated as lying in the tangential plane (e η , e ξ ), that is, n f = cos αe η + sin αe ξ , where α is the fiber angle. There is a sheet angle β between n s and e h , so that n s = cos βe h + sin β(e h ¢ n f ). Finally, n n = n f ¢ n s . Angles α and β change transmurally and from base to apex. According to the experiment observations [36,48], a linear transmural distribution is assumed for α and a quadratic distribution is assumed for β, that is, where t is the relative wall depth measured from the epicardial surface, the angles (α epi , α endo , β endo , β mid , β epi ) are assumed vary linearly from the basal site (z Base = 2.5 cm) to the apical site (z Apex = 5.0 cm), as where h = (z Base z)/(z Base z Apex ). Based on the measured distribution of angles α and β in canine ventricular myocardium [48], the following architectural parameters are assumed (in degree) for the simulations:

Linear anisotropic viscoelastic passive behavior and active contraction stress
In the present phantom simulations, a linear anisotropic viscoelastic model is employed for the passive behaviors of normal and ischemic myocardial tissues, and a directional active contraction stress is adopted, that is, (17) in which the strain energy W E is calculated on the local structure-based material axes (n f , n s , n n ) as where ε f f = n f ¡ ε ¡ n f , ε ss = n s ¡ ε ¡ n s , ε nn = n n ¡ ε ¡ n n , ε f s = n f ¡ε ¡n s , ε f n = n f ¡ε ¡n n , and ε ns = n n ¡ε ¡n s . This is a linear approximation to the exponential strain energy employed in the forward heart simulations of Guccione et al. [8] and Usyk et al. [12]. An isotropic Voigt-type viscosity is assumed. The active contraction stress is assumed directly related to the fiber direction n f . It depends on the current calcium concentration Ca 2+ and the relative sarcomere stretch λ, which is calculated with λ = n f ¡ε ¡n f + 1 for small-strain approach. There have been complex models for T f [49,50] in heart muscle, and three-dimensional description of active stress σ f [51,52]. For the purpose of the present DCE phantom simulations, a simplified model of Rachiev and Hayashi [53] is employed: in which λ 0 = 0.8 and λ m = 1.2 are the reference sarcomere stretches for myocardium [54], and T max is the reference contraction stress. The material parameters p = (E f f , E ss , E nn , E ns , E f n , E f s , γ, T max ) are given for normal and ischemic myocardial where E's and T max are in kPa, and γ is in kPa ¡ second. The ratios between E's reflect the material anisotropy, in consistency with the exponential model of Guccione et al. [8] and Usyk et al. [12]. The magnitudes of E's are chosen such that the end-filling sarcomere stretch lies mostly around 1.1 in the LV wall, and that the ischemic region are three times as stiff as the normal region [4]. For the normal tissue, a reference contractility T max = 800 kPa is chosen for most realistic LV contraction deformation. The ratio between T max and the passive stiffness E's is consistent with those found in the literature [18]. The contractility of the ischemic myocardium is typically much weaker than the normal tissue [5,55], and changes with the degree of damage and healing. Therefore, its T max is assumed at 200 kPa, that is, 25% of the normal tissue. Although the normal tissue is assigned with one set of material parameters, it is heterogeneous owning to the spatial variation of myofiber architecture. So is the ischemic tissue. It is noted that the objective of the present simulations is to investigate the feasibility of cardiac dynamic elastography, instead of the actual mechanical responses in a real heart. Therefore, the phantom, loading conditions, material models, and parameters are not exact.
Material parameters (E f f , E ss , E nn , E ns , E f n , E f s , γ, T max ) are to be identified in the following DCE reconstructions. The reference stretches λ 0 and λ m show relatively small variations in different hearts [54], and are assumed known in the simulations. The mass density of myocardium is also assumed known at ρ 0 = 1g/cm 3 . Kinetics of Ca 2+ is described in the next subsection.

LV pressure and calcium kinetics profiles
An empirical LV blood pressure profile P(t) (e.g., [16]) is used, as in Figure 4(a). The passive phase involves rapid/reduced filling (time t < t 0 = 0.34 second) and artery contraction (t 0 t < t 1 = 0.48). The active phase includes LV isovolumetric contraction (t 1 t < t 2 = 0.60) and rapid/reduced ejection (t 2 t < t 3 = 0.88). In the active phase, the calcium concentration Ca 2+ (t) is assumed uniform in the LV wall, and is modeled mathematically as [49] Ca where k = 10, as in Figure 4(b). In practice, LV blood pressure and electromechanical activity in an individual heart can be obtained noninvasively with established techniques [56,57], and are used for DCE reconstructions.

Dynamic deformation of LV
In the phantom simulations, the endocardial surface is subjected to uniform blood pressure P(t) (Figure 4(a)), and the epicardial surface is force free. The basal surface (z = 2.5 cm) is fixed, that is, u = 0 at z = 2.5 cm. With material model  made with a normal phantom without the ischemic region. The differences of epicardial/endocardial displacement components, u x , u y , and u z , between ischemic and normal phantoms are plotted in Figures 5(b), 5(c), and 5(d) for t = 0.4 second and Figures 6(b), 6(c), and 6(d) for t = 0.6 second, respectively. The comparison shows that the ischemic region, which is passively stiffer and actively weaker, changes the displacement characteristics. In other words, the material parameters of normal and ischemic tissues may be extracted from the displacement measurements.

DCE reconstruction for the passive viscoelastic parameters
The DCE reconstruction is divided into two steps. First, the viscoelastic parameters E's and γ are identified using the epicardial/endocardial displacement measured in the passive phase, where no active contraction occurs in the LV wall. Then, the active contraction stress T max is reconstructed using the displacement in the active phase. The reconstruction follows the optimization procedure (Figure 1), and uses the proposed adjoint gradients (9), (10). The epicardial/endocardial displacement is measured at time t = 0.1, 0.2, 0.3, 0.4 second (Figure 4(a)). The initial guess is picked near to the real parameters of the normal tissue, in order to reach the global minimum of the objective function Φ(p), which corresponds to the exact material parameters. In practice, the range of the normal tissue parameters may be obtained from available database. Therefore, such an initial guess is feasible. The results are presented in row "Ideal Input" of Table 1. The reconstruction converges after about 250 iteration steps. It yields the exact passive viscoelastic parameters p Normal for the normal tissue, and very accurate p Ischemic for the ischemic tissue. Very fast convergence of p Normal is observed. For instance, parameters (E f f , E ss , E nn , E ns , E f n , E f s ; γ) of normal tissue reach (188.8, 37.45, 37.56, 28.14, 26.82, 29.39; 0.844) at the 25th iterative step, and (185. 59, 36.24, 36.33, 28.35, 28.06, 28.85; 0.922) at the 50th step, then the main adjustment is for the viscosity coefficient γ. For the ischemic tissue, E f f climbs steadily from the initial guess 125.0, reaches the highest point of 680 at around the 200th step, and decreases gradually to the result of 559.3. The other parameters in p Ischemic , however, experiences fluctuations. There seems to be a combining effect of the elasticity E's and viscosity γ in the monotonous LV filling motion, that is, both contribute to the effective stiffness. The weakening effect owning to low elasticity E's may be balanced by a high viscosity γ. During the reconstruction, E's and γ both growth from the initial guess, with a much higher increasing rate for the latter. After about 150 steps, when p Normal have converged very closely to the exact values, the ischemic γ begins to decrease from the highest value of about 6.6, and the E's increase with faster speed toward the exact values. For the ischemic tissue, the faster convergence and better accuracy of E f f indicates its relatively stronger role in the LV deformation. The primary deformation pattern in the passive phase is LV dilatation, where the dominant strain-stress relation occurs mainly along the myofiber direction, and the shear strains on the material axes (n f , n s , n n ) are relatively weak. Therefore, the influences of E f f are most significant. As a result, E f f for ischemic region is much better identified than the other E's. There is also size effect, that is, small size of the ischemic region weakens the effects of p Ischemic on the overall LV motion. As shown in Figures 5 and 6, deformation variations due to the ischemic region are localized, that is, occur mainly around the ischemic region. This may explain the much faster convergence and higher accuracy for p Normal than p Ischemic . Furthermore, DCE reconstruction has been performed using displacement at nine time steps, that is, t = 0.05, 0.1, 0.15, . . . , 0.45 second (Figure 4(a), both solid and open dots). The improvement is not significant. Better results are expected if transmural displacement field can be measured, for instance, with MRI tagging technique [58]. Reconstructions have also been performed using the direct gradients (11), (12), and similar results are obtained. However, the computational expense is much higher, by a factor of about 5, noting that there are 14 displacement fields u (i) = ∂u/∂p i to solve (12).
In comparison with our previous static anisotropic elastography simulations [44], the present DCE needs more iterative steps and shows overall lower accuracy. This is due to the approximate numerical scheme for solving the dynamic motion equations (2), (9), that is, Newmark-β method in this work.

DCE reconstruction for the active contraction parameters
This section identifies the active contraction stresses T max (19) for the normal and ischemic tissues, denoted as T Normal and T Ischemic , respectively. The epicardial/endocardial displacement is measured at time t = 0.5, 0.6, 0.7, 0.8 second during the active phase (Figure 4(a)). Two approaches are made.
The first approach employs the exact passive material parameters given in (20). The initial guess is 500.00 for both T Normal and T Ischemic . It is near to the exact value for the normal tissue. Thirteen iterative steps are needed to get convergent results, T Normal = 800.0 and T Ischemic = 200.0, respectively, which are the exact values. Similar to the above observation on the convergence of p Normal and p Ischemic , T Normal converges quickly toward the exact value, then T Ischemic accelerates its convergent speed. More specifically, T Normal has reached 806.8 at the 8th iterative reconstruction step, while T Ischemic is at 445.2, about twice of the exact value. Then, T Normal gradually drops to the exact value, and T Ischemic approaches rapidly to 200.0. As discussed above for p Ischemic , two factors probably contribute to the late convergence of T Ischemic , that is, the small size of ischemic region and the low value of T Ischemic (25% of T Normal ).
The second approach uses the previously reconstructed passive viscoelastic material parameters (E f f , E ss , E nn , E ns , E f n , E f s ; γ) (row "Ideal Input" of Table 1). The pattern of convergence is similar to the above with exact passive parameters, while there are some reconstruction errors. At the 8th iterative step, T Normal = 805.6 and T Ischemic = 446.9. Then T Ischemic starts to converge quickly. After 13 steps, the gradients ∂Φ/∂p approach zero, and reconstruction results are T Normal = 799.5 and T Ischemic = 197.8.

DISCUSSIONS
The above DCE reconstructions use exact measurements, that is, exact myofiber architecture, exact epicardial and endocardial displacement, and exact blood pressure. The following simulations investigate the influences of the errors with myofiber architecture and displacement measurements, respectively. Furthermore, a comparison is made between the dynamic cardiac elastography and the previous quasistatic approach.

Errors with myofiber architecture
In vivo measurement of the myofiber architecture has been a challenge in cardiac imaging techniques. Recently, Scollan et al. [38] developed a diffusion-tensor MRI to estimate the fiber orientation in a living heart. Another feasible way is to approximate the myofiber architecture using available ex vivo experimental database [36,55,[59][60][61], based on the understanding that variation of myofiber architecture in individual hearts is relatively small after some anatomy-based mapping.
To simulate the errors with myofiber architecture, randomly picked values between 5 and 5 degrees are added to the "exact" fiber angle α and sheet angle β (14)- (16) in each element. Reconstruction is then conducted to identify the passive parameters (E f f , E ss , E nn , E ns , E f n , E f s ; γ). The results are given in the row "Errors-F" of Due to the small size and weak influences of the ischemic region on the LV deformation, p Ischemic are more sensitive, and show reconstruction errors as large as 10%. The overall patterns of convergence are similar to those of DCE reconstruction with exact myofiber architecture (Section 3.5).

Errors with displacement measurement
Epicardial/endocardial displacement of a beating heart can be computed from dynamic cardiac imaging sequences (Lin & Robb [34], Eusemann et al. [35]), and inevitably contains errors. The response of DCE reconstruction to inaccurate displacement measurement is investigated herein. Randomly generated 5% to +5% relative error is added to each "exact" epicardial/endocardial displacement component. The results are shown in the row "Errors-D" of Table 1.
The parameters p Normal converge stably, and reach accurate results. In comparison, p Ischemic are affected by the displacement errors more significantly. The ischemic parameters E ss , E nn , E ns , E f n , and E f s show error as high as about 15%. However, the material anisotropy, described by the ratios between E's is maintained very well, as well as the heterogeneity due to the ratios between p Normal and p Ischemic . The ischemic viscosity parameter γ shows the largest error, that is, the result is 2.51, compared to the exact value 1.00. This may probably be explained by the combining effect of E's and γ, as mentioned above, and the weak influence of ischemic γ on the deformation, which may be somehow concealed by the errors with displacement measurements. DCE reconstruction has also been conducted with both ¦5 degree myofiber architecture errors and ¦5% displacement errors. The errors are randomly picked, and are different from the above. The results are given in the row "Errors-FD" of Table 1. In overall, p Normal are robust against these errors, with maximum error 7% for γ. While p Ischemic are more sensitive to the measurement errors, the material anisotropy and heterogeneity are well maintained. In fact, the resulting ischemic E's are somehow more accurate than those with displacement errors only, probably due to the differences in the random errors. Based on our other DCE reconstructions with randomly picked measurement errors, it is shown that the influences of myofiber architecture errors are in general less significant than the displacement errors. This could be explained with concept of homogenization, that is, the random myofiber architecture errors at a material point may be compensated by the errors in its neighborhood. In fact, the perturbation on deformation due to material variation at a point is highly localized in a solid. In view of the less accurate reconstruction results for the viscosity parameters, we speculate that methods may be needed to bring out or isolate the viscosity effects, or instance, a modality to use ultrasonic impulse to stimulate a low-magnitude jump-type deformation in the myocardium.

On the quasistatic approach
The previous cardiac elastography studies [28][29][30][31] made quasistatic assumption for the LV passive (filling) deformation. Consequently, the effects of inertia and viscosity were omitted. In those studies, the displacement is typically taken near to the highest filling pressure, that is, around t = 0.4 or 0.45 second for the present blood pressure profile (Figure 4(a)). To investigate the quasistatic assumption, reconstruction is conducted by deactivating all the terms related to inertia and viscosity in our DCE algorithm and the codes. Exact myofiber architecture data and displacement at t = 0.4 and 0.45 second are used. The results are given in the last three rows of Table 1. "Static A" and "Static B" use displacement measured at 0.4 and 0.45 second, respectively, and "Static C" uses displacement at both 0.4 and 0.45 second. Note that the viscosity parameters γ cannot be identified with quasistatic assumption.
It is shown clearly that the quasistatic elastography does not yield accurate material parameters for either the normal or the ischemic myocardium. This may indicate the necessity of dynamic description for cardiac elastography. In a real heart, the left ventricle dilates to about the maximum chamber volume near to the end-filling stage, where the velocity field in the wall is overall very low and the effects of viscosity are about minimal. From this point of view, quasistatic elastography is probably a good approach for the elastic parameters E's, and failure of the present quasistatic reconstruction may only indicate that the "real" material parameters we chose do not give almost static deformation at the highest filling blood pressure (around t = 0.4 and 0.45 second). On the other hand, it is not clear if the real acceleration field at the end-filling stage is low enough for omitting the inertia effects. Another concern is that, even at the end-filling stage where the overall velocity could be very low, the quasistatic displacement may still be qualitatively  different from the dynamic displacement. Consider a onedimensional spring-mass-dashpot system, for example, the dynamic maximum displacement (at zero velocity) is higher than the quasistatic maximum displacement under the same force. For the present LV phantom, the differences between dynamic and static displacements on the epicardial surface, with the same material parameters and blood pressure, are given in Figure 7 at t = 0.45 second. It is shown that the difference are global and have clear distribution patterns, compared to the much more localized differences due to the ischemic region (Figures 5 and 6), and are much more significant in magnitude than the latter. This suggests that the errors due to quasistatic assumption may overwhelm the effects of the ischemic region. In fact, the quasistatic reconstruction results for p Ischemic (Table 1) are qualitatively incorrect. Therefore, quasistatic cardiac elastography needs careful evaluations, particularly using in vivo measurements. For the viscosity and active contraction parameters, like γ and T max in the present model, however, the proposed dynamic cardiac elastography is indispensable.

CONCLUSION
In this work, a dynamic cardiac elastography (DCE) framework has been developed for identification of the material parameters for the passive viscoelastic behaviors and active contractility of myocardium, using time-dependent displacement extracted from dynamic cardiac imaging sequence. The DCE method takes into consideration the myofiber architecture-induced orthotropy and heterogeneity of the material properties. The reconstruction algorithm employs an iterative optimization procedure, and uses usersupplied gradients of the objective function. A dynamic adjoint method is derived for cost-efficient calculation of the gradients. DCE simulations have been conducted to identify the passive and active parameters for normal and ischemic myocardium in a left ventricle phantom. For the normal myocardium, the passive viscoelastic parameters and reference active contraction stress are identified with high accuracy, and are not considerably sensitive to the errors with the myofiber architecture and displacement measurements. For the small ischemic region in the ventricular wall, which is stiffer passively and weaker actively than the normal tissue, the reconstruction results are very accurate with ideal measurements, and are satisfactory when errors are added to the myofiber architecture and displacement measurements. The previous quasistatic assumption for cardiac elastography is investigated. It is suggested that consideration of dynamic deformation may be necessary, particularly for the viscosity and active contraction parameters.