INTRODUCTION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION RESULTS DISCUSSION REFERENCES

CARDIAC MECHANICAL PERFORMANCE depends largely on structural properties of the myocardium and the behavior of its cells. Moreover, myocardial properties may vary due to long-term changes in such factors as hemodynamic load and electrical activation sequence. For example, in eccentric and concentric hypertrophies, the ventricular thickness is adjusted to restore normal relation between the inner radius and the wall thickness (see Ref. 27 for a review). Volumetric growth can also be a local phenomenon. During asynchronous electrical activation, a reduction of the wall thickness in the early activated region was observed (23), which was probably induced by the heterogeneity of fiber strain distribution due to asynchronous electric activation (22). In some cases such an adaptation of the ventricle was associated with fiber reorientation, e.g., during healing of an infarcted region (32), whereas in other cases fiber orientation in the wall remained similar to the orientation in the normal heart (7, 20). Regardless to the many cases of adaptation described in the literature, it is still not clear which signal is used at the cellular level to regulate the growth and remodeling of the cardiac tissue. A growing body of evidence points to some mechanical factor as a stimulus (19). Some potential candidates for this mechanical factor are fiber stress, fiber strain, generated mechanical work, or ATP consumed. Because experimental assessment of these quantities with sufficient spatial resolution is difficult, mathematical models have been developed to predict them. In one of these models, measured wall geometry and fiber orientation were used to compute the stress and strain in the wall, and inhomogeneous distributions of stress and strain were obtained (13). In other model studies, it has been found that distributions of stress and strain are very sensitive to changes in fiber orientation. With the fiber orientation varied within the range of measured values, distributions of developed stress were found to be virtually homogeneous, or strongly inhomogeneous, with stress levels varying by more than a factor 2 (3).

Because the wide range of available experimental fiber orientation data excludes a reliable prediction of local mechanics in mathematical models, model fiber orientation has been optimized for homogeneous spatial distribution of local mechanics. The rationale behind this approach was the assumption that myofibers would strive for the same optimal mechanical load. In one study, active myofiber stress was chosen as the relevant aspect of mechanical load, and fiber orientation was optimized for optimal homogeneity of myofiber stress (3). In other studies the variation of fiber strain at the beginning of ejection (24) or during ejection (25) was minimized, and a fiber orientation close to the measured one was predicted. The latter studies suggest that fiber reorientation might be an attractive mechanism for the cell to adapt to changes in mechanical load (1).

There are several issues that were not addressed in the model studies mentioned above. First, it is not clear whether the fiber orientation obtained by optimizing strain or stress distribution will be different depending on which mechanical factor is optimized. Second, it is not clear whether another potential stimulus, e.g., consumed ATP, would yield the same fiber orientation. Third, it is not clear whether such an adaptation of fiber orientation would make the heart function better as a pump. For example, the question is whether such an adaptation would increase the ejection fraction or cardiac efficiency. In other words, Rijcken et al. (24, 25) did not address the main output and input of the left ventricle during optimization of fiber orientation.

The aim of this study is to investigate the influence of fiber orientation on the ejection fraction, efficiency, and the heterogeneity of the distributions of fiber stress, fiber strain, and ATP consumption. A finite element model similar to that of Ref. 3 was used with active properties described by the Huxley-type cross-bridge model (36).

There are several differences between this study and earlier studies (3, 24, 25). First, we used a dynamic model (computing the state of the ventricle during a cycle) to study the influence of fiber orientation on the distribution of fiber stress and strain. In Bovendeerd et al. (3) a dynamic model was used, but very few fiber orientation distributions were considered. In the study by Rijcken et al. (24, 25), only two systolic states were considered with measured ventricular pressure and volume used as a model input. Second, because we used a dynamic model, we were able to relate ejection fraction of the ventricle to the heterogeneity of fiber stress and strain-an aspect missing in the previous studies. Third, we used a Huxley-type cross-bridge model for description of the active properties of the muscle (36). Because Huxley-type cross-bridge models relate mechanical properties of the muscle to the consumption of ATP by cross-bridges, it was possible to study how the variation of the fiber orientation influences the energy consumption distribution within the ventricle.

	MODEL DESCRIPTION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION RESULTS DISCUSSION REFERENCES

Geometry. In the reference state of the model, defined as the state with zero transmural pressure, the endocardial and epicardial surfaces are represented by truncated focal ellipsoids (29), leaving a thick wall between them. The volume of left ventricular wall and cavity was set to 142 and 40 ml, respectively. On the basis of geometrical data presented by Streeter and Hanna (29), the papillary muscle volume was set to 4 ml and the common focal length (C) of the ellipsoids was set to 43 mm.

(1)

(2)

(3)

(4)

(5)

View larger version (23K): [in this window] [in a new window]   Fig. 1.   Fiber orientation in the left ventricle is quantified by the transverse angle (t) shown in the transverse cross section of the left ventricle (A) and by the helix fiber angle (h) (B). The fibers (bold lines) with helix angle of +45° (subendocardial), 45° (subepicardial), and 0° (middle) are shown in B (fine distribution of finite elements is not shown).

Constitutive behavior. In the model, myocardial tissue is assumed to consist of fluid, a connective tissue matrix, and muscle fibers. The total Cauchy stress () developed in myocardial tissue is divided into the following: 1) the uniaxial active stress (_a) generated by the contractile element parallel to the muscle fiber direction vector (e_f), 2) the three-dimensional passive stress (_p) resulting from the tissue deformation, and 3) hydrostatic pressure (pI) of the fluid trapped in the solid

(6)

where I is the unity tensor.

(7)

(8)

(9)

(10)

(11)

(12)

View larger version (12K): [in this window] [in a new window]   Fig. 2.   The stiffness factor in normal and stiffened simulations. Eq, equator.

Governing equations and initial and boundary conditions. Calculations are based on the law of conservation of momentum. Neglecting inertial (16) and gravitational effects, the conservation of momentum is given by

(13)

The mechanical activation of the left ventricle is assumed to be simultaneous in the most of the simulations. To check the sensitivity of the model to the activation sequence, the electrical activation sequence adopted by Bovendeerd et al. (2) was used in some simulations. Here, activation starts in the subendocardial apical region and reaches the epicardial surface at the base at ~40 ms, in agreement with experimental data (8).

Analysis. The regional differences of sarcomere strain, stress, and ATP consumption are quantified by average variation of these functions during systole and relaxation of the ventricle. The average variation of function f is given by

(14)

where r is a position vector,  is a domain in which the variance is found, V is the volume of , is the average of function f in  at time moment t, and t_bs and t_bd are time moments at the beginning of systole and beginning of diastole, respectively. The time moment t_bd is defined as a time moment at which the relaxing ventricle has a pressure of 1 kPa. Because of the definition of var(f) used here, the units of var(f) and f are the same. To exclude the influence of stress concentration in apex and basal regions to the computation of the variance, the domain  included the complete left ventricle except the rings with the thickness of two finite elements in the apex region and one element near the base, i.e.,   (1.16, 2.70). The variances of the following functions were computed in this work: half-sarcomere length (l_s strain in fiber direction), _a, and ATP consumption during a cycle V. When ATP consumption variance was computed, no integration by time was required

(15)

Numerical methods. In the model, the state of stress and strain is known, once the cross-bridge distribution functions ["internal variables" for our model (9)] and the position vector r are given for every material point in the tissue. The governing equations were discretized using finite element method in conjunction with Galerkin's method. As a result of discretization, the system of nonlinear equations was composed. The unknowns in every finite element node were as follows: displacement r_i, the cross-bridge distribution functions, and the sarcomere shortening velocity. Additional unknowns were the left ventricular pressure and the aortic pressure. To find the unknowns, time was discretized and it was assumed that the sarcomere shortening velocity and efflux of the blood from the left ventricle was constant during each time step. The composed system of nonlinear equations was solved using the Newton iterative method implemented with the use of NITSOL software (21). The change in the cross-bridge distribution functions within the time step was found by integrating the corresponding equations using DVODE software (5). Finite element discretization was performed with the use of Diffpack software (6). We used a 20-node mixed finite element with displacement approximated using all 20 nodes, stress approximated by trilinear functions using only 8 corner nodes, and elementwise constant pressure approximation. The accuracy of solution was tested by comparison of different spatial discretizations and varying the time-step duration. According to our tests, the variances of sarcomere length, _a, and ATP consumption changed <10% when the following grid was refined by doubling elements in all directions: 6 elements between endocardial and epicardial surfaces and 10 elements in apex-base direction. In the circumferential direction, we used two elements because only one-eighth of the left ventricle was simulated. The time step used in all simulations was 10 ms. The reduction of time step to 2.5 ms changed the computed variances <3%.

Simulations performed. First, optimization was performed with respect to var(l_s/2), var(_a), and var(V) at P₃ = 5° over a large (P₁,P₂) design space: P₁ and P₂ were varied from 40° to +60° and 100° to 0°, respectively. To prevent numerical instabilities, these simulations were performed in the stiffened left ventricle. Mechanical activation was chosen either simultaneous or the same as the electrical activation sequence to check the sensitivity of var(l_s/2), var(_a), and var(V) to the changes in the mechanical activation sequence.

	RESULTS

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION RESULTS DISCUSSION REFERENCES

Fiber orientation. A typical dependence of the sarcomere length, developed stress and ATP consumption variances on helix fiber angle transmural distribution is presented in Fig. 3. Angles P₁ and P₂ represent the _h in the midwall (P₁) and the slope (P₂) of transmural _h changes. There are several local minima of the variances in the (P₁,P₂) plane. In this simulation, the smallest variance of sarcomere length was found at P₁ = 20° and P₂ = 70°. It is clear that only one minimum provides required ejection fraction and high efficiency (mechanical work performed to pump blood divided by amount of consumed ATP) of the ventricle (Fig. 3, D and E). In the simulations presented in Fig. 3, we activated the left ventricle simultaneously. If the mechanical activation of the ventricle was assumed to be the same as electrical activation sequence, then similar results were obtained (Fig. 4).

View larger version (55K): [in this window] [in a new window]   Fig. 3.   The variances of half-sarcomere length (strain, A), developed stress (B), and ATP consumption (VATP) during a beat (C) are compared with the ejection fraction (D) and the relative efficiency of the left ventricle (E) at different helix fiber angle transmural distributions. The simulations were performed with a stiffened model, which resulted in smaller ejection fraction and smaller differences in the variances. In the top left corner of A-E, the computations did not converge. Note that there is only one region (at P1 and P2 angles approximately +20° and 70°, respectively) in the (P1, P2) plane with relatively small variances and high ejection fraction. In these simulations, angle P3 was equal to 5° and the left ventricle was activated simultaneously.

View larger version (41K): [in this window] [in a new window]   Fig. 4.   The computed variance of half-sarcomere length (strain, A) and ejection fraction (B) for the heterogeneously activated left ventricle. The mechanical activation was taken the same as electrical activation of the ventricle. Note that the absolute values of the half-sarcomere length variance increased if compared with the variance presented in Fig. 3. However, the overall pattern of the variances, ejection fraction, and the efficiency dependency on the helix fiber angle transmural distribution is very similar to the one obtained for the simultaneously activated left ventricle (see Fig. 3).

View larger version (16K): [in this window] [in a new window]   Fig. 5.   Helix fiber angle found by optimizing sarcomere length (strain), developed stress, and ATP consumption distribution within the left ventricular wall at different transverse (P3) fiber orientations. The h in the midwall (P1) and the slope (P2) of transmural h changes are shown in A and B, respectively. Error bars indicate the discretization level of the mesh in the (P1,P2) plane used to find the optimal location. Note the small difference between the fiber orientations predicted by optimizing strain, stress, or ATP consumption distributions.

View larger version (18K): [in this window] [in a new window]   Fig. 6.   The relative variances of sarcomere length (strain), developed stress, and ATP consumption computed as functions of transverse fiber orientation P3 using optimal values of P1 and P2 (see Fig. 5). Note the differences in optimal P3 angle predicted by minimum of strain, stress, or ATP consumption variances.

View larger version (14K): [in this window] [in a new window]   Fig. 7.   The normalized ejection fraction and efficiency of the left ventricle as functions of transverse fiber orientation P3 at P1 and P2 equal to 22.5° and 69°, respectively. The measured helix fiber distribution (P3 value around 10-15°) corresponds to the position of the maxima of the ejection fraction and efficiency of the left ventricle.

View larger version (14K): [in this window] [in a new window]   Fig. 8.   Transmural distribution of h at equator predicted by optimizing sarcomere strain, developed stress, and ATP consumption distributions. Predicted h is the same regardless to the distribution used in the optimization. The measurements are of the equatorial region of the human left ventricle () (28) and the dog left ventricular wall at anterior (), lateral (), and posterior () sites (17). Endo, endocardial; Epi, epicardial.

View larger version (12K): [in this window] [in a new window]   Fig. 9.   Distribution of transverse angle t at midwall region predicted by the model when either sarcomere strain, developed stress, and ATP consumption distributions were optimized. Prediction is compared with the measurements () performed by Bovendeerd et al. (4).

Model testing. In the test simulations, we used the fiber orientation angles predicted by minimizing sarcomere length variance: P₁ = 22.5°, P₂ = 69°, and P₃ = 5°. The computed hemodynamic properties and deformation of the ventricle were as follows. The peak systolic pressure was equal to 21.8 kPa. The relative increase in equatorial wall thickness (_d), outer equatorial ventricular radius (_R), and outer ventricular length (_L) during a systole were +23%, 7%, and 3%, respectively. According to the measurements of Olsen et al. (18), end-systolic values of _d, _R, and _L were approximately +17%, 10%, and 5%, respectively.

View larger version (22K): [in this window] [in a new window]   Fig. 10.   The apex rotation angle-pressure relationship computed by the model is compared with the measurements of Gibbons Kroeker et al. (10) (). The positive apex rotation angle denoted clockwise rotation of the apex with respect to the base (as viewed from the apex). EIVC, end of isovolumetric contraction period; ES, end of systole. The computations were performed with fiber orientation predicted by minimizing the variance of sarcomere length (, for simultaneous left ventricle activation and dashed line with  for nonsimultaneous activation) or developed stress (, for simultaneous activation). The symbols on the computed lines indicate EIVC and ES points. Note that the apex rotation angle-pressure relationship is sensitive to the fiber orientation and mechanical activation sequence.

1

View larger version (20K): [in this window] [in a new window]   Fig. 11.   The pressure-volume relationship computed for isovolumetric and ejecting contractions of the left ventricle. The volume-axis intercept of the end-systolic pressure-volume relationship (ESPVR) was found by approximating the ESPVR by a linear function. The linear ESPVR was fitted using end-systolic points of isovolumetric () and ejecting () contractions.

View larger version (12K): [in this window] [in a new window]   Fig. 12.   ATP consumption of the complete left ventricle as a function of the pressure-volume area (PVA) for isovolumetric () and ejecting () contractions. The linear function (solid line) was fitted using isovolumetric contraction data. Note that ATP consumption-PVA relationship is linear and is the same for isovolumetric and ejecting contractions, in correspondence with measurements (30).

View larger version (39K): [in this window] [in a new window]   Fig. 13.   Average ATP consumption within the left ventricular wall predicted by the model. The contour lines are shown at the following values: 1.0, 1.2, 1.3, and 1.4 ATP molecules per myosin head. Note that the highest ATP consumption is between the endocardium and middle of the ventricular wall.

	DISCUSSION

TOP ABSTRACT INTRODUCTION MODEL DESCRIPTION RESULTS DISCUSSION REFERENCES

According to our simulations, the variances of the sarcomere length, developed stress, and ATP consumption during a beat have very similar dependencies on transmural course of _h. The optimal transverse angle value is also similar if the variance of the sarcomere length or developed stress is minimized. The dependence of sarcomere length, developed stress, and ATP consumption variances on _h distribution is not simple: the variances have several minima at different _h distributions. However, we identified only one region in the (P₁,P₂) plane with high ejection fraction and high efficiency of the left ventricle and relatively homogeneous distributions of sarcomere strain, developed stress, and ATP consumption within the ventricular wall.

In this study, we tested whether the fibers may be oriented to minimize the heterogeneity in either stress or strain distributions from the theoretical point of view. If this hypothesis is correct and the strain or stress distribution is indeed regulating the fiber orientation in the ventricular wall, then two criteria have to be met. First, the fiber orientation predicted by minimizing the variance of the strain or stress distributions should be close to the measured orientation. Second, from an evolutionary point of view, the resulting fiber orientation should lead to better cardiac performance, i.e., the heart should function better as a pump. These criteria are required but not sufficient to prove the hypothesis, and we were only able to check whether we can exclude some stimuli if one of these conditions was unsatisfied.

The predicted fiber orientation angles are close to the measured ones (Figs. 8 and 9), in accordance with the earlier studies performed on the models where only two systolic states were considered (24, 25). The difference between predicted and measured _h in subendocardial and subepicardial regions is caused by a linear approximation of _h. The proportional approximation of _t may be also a reason of underestimation of this angle near the apex and the base (Fig. 9). When more complex functions were used to approximate _h and _t, a much better fit between predicted and measured fiber orientation was found (4). However, it would be very difficult (if at all possible) to perform the present analysis with 12 parameters describing _h and _t as in Ref. 4. Namely, the computation of the left ventricle deformation, energy consumption and ejection fraction in 12-dimensional parameter space would require many simulations and a very complicated analysis of obtained solutions.

From the results obtained in this study (Figs. 3 and 4), we conclude that the fiber orientation obtained by minimizing the variance of the strain or stress distributions leads to high left ventricular ejection fraction and stability in design. By comparing the variances of the sarcomere strain and developed stress with ejection fraction at different fiber orientations (Fig. 3), we have shown that a relatively homogeneous distribution of strain and stress leads to high ejection fraction and efficiency of the left ventricle. It is important to note that there are local minima of the variances in the high ejection fraction region in the (P₁,P₂) plane. Thus, if the strain or stress are used as stimulus, the fiber orientation would then be stabilized in the high ejection fraction region in the (P₁,P₂) plane. If small changes in the fiber orientation were to occur, then the mechanical stimulus should return the fiber orientation to the original one, i.e., the fiber orientation is stable. On these terms, the ATP consumption variance is rather small, thus indicating that relatively homogeneous coronary perfusion of the left ventricle is required.

In earlier studies (3, 24, 25), only fiber strain or developed stress were used as a stimuli to find the fiber orientation in the ventricle. Here we tested whether another potential stimulus, e.g., consumed ATP, would yield the same fiber orientation. However, the fiber orientation obtained by minimizing ATP consumption variance reproduced the measured _h only. The _t obtained from ATP consumption distribution was equal to zero, as opposed to the measured data (see Fig. 9). Partially, the misprediction of _t may be caused by the changes of the average ATP consumption when P₃ is varied. Namely, the amount of consumed ATP is maximal at P₃ = 20° (not shown) and is ~15% smaller at P₃ = 10° when only P₃ is varied at fixed P₁ and P₂ values close to the optimal (P₁,P₂), as in Fig. 7. Such changes in average ATP consumption are possibly reducing the absolute values of the differences in the consumption at the different left ventricular wall positions. To test this hypothesis, we normalized the ATP consumption variance (Fig. 6) by the average value of the consumption. This shifted the minima to P₃ = 2.5  5°. The similar procedure applied to the variances of sarcomere stress and developed stress did not affect the position of the minima (Fig. 6). The shift of the minima positions was insensitive to the selected (P₁,P₂) when either P₁ or P₂ was modified by ±2°. So the _t predicted by ATP consumption variance normalized by the average value is closer to the measured data, but is still underestimated. This may indicate that ATP consumption is not used as a signal orienting the fibers in the left ventricle. Whether this conclusion is a correct one or is caused by model limitations is not clear and requires further investigation.

We tested the model and predicted fiber orientation by comparing the model solution with the experimental data on left ventricle deformation and oxygen consumption. The largest difference between the model solution and experimental data that we have identified was the torsion of the apex (Fig. 10). Torsion of the apex was sensitive to the changes in the fiber orientation and the activation sequence of the ventricle. Both the fiber orientation and the activation sequence are influencing the apex rotation angle through the same mechanism-by changing the balance between epicardial and endocardial layers of the myocardium: the fiber orientation determines the direction of the developed force and the activation sequence determines when the force in particular direction is developed. In our simulations, we reproduced the apex rotation direction and obtained the apex rotation angle-pressure relationship in the form of the loop similar to the measurements (10). However, we were not able to reproduce the apex rotation angle-pressure loop quantitatively. This is most probably caused by the approximations used in the model: 1) linear approximation of _h and _t and/or 2) homogeneous mechanical activation.

In addition to the model limitations discussed above, several simplifications were made. First, we used an axisymmetric model of the left ventricle. Thus it is impossible to study the fiber orientation differences in the different parts of the left ventricle without modification of the model. Second, the boundary conditions used in the basal surface were quite simple, ignoring the constraints imposed by the basal skeleton. In addition, the constraints imposed by the right ventricle were ignored. Third, we have not considered the influence of the laminar structure of the myocardium to distributions of stress and strain in the left ventricular wall. However, because we studied the deformation of the left ventricle during the systole and the laminar sheets are of minor importance during this period of the cardiac cycle (33), the simplification used in our model was justified. Despite all of these simplifications, the found fiber orientation distribution was close the measured one, the deformation of the ventricle resembles the experimental measurements, and several important mechanoenergetic properties of the heart were reproduced (see below).

This is the first time we used the Huxley-type cross-bridge model to simulate the active properties of the cardiac muscle in the left ventricle model. Because it is possible to compute ATP consumption directly from the model equations, we were able to find ATP consumption of the ventricle, relate it to PVA (Fig. 12), and predict the distribution of ATP consumption in the left ventricle wall (Fig. 13).

One of the important properties of left ventricle is a linear relationship between the PVA and oxygen consumption (30). The similar property of myocardium has been identified on the tissue level-linear relationship between the stress-strain area (SSA) and oxygen consumption (14). Assuming that ATP consumption by excitation-contraction coupling and basal metabolism is almost constant regardless of PVA in a given contractile state, one can conclude from the linear PVA-oxygen consumption relationship (30) that the PVA-ATP consumption relationship should be linear, too. The computed PVA-ATP consumption and SSA-ATP consumption relationships are both linear and independent on the type of contractions (ejecting and isovolumetric contractions) reproducing the measured data quantitatively [see Fig. 12 and Vendelin et al. (36)]. Taking into account that we used the cross-bridge model, which reproduces the SSA-ATP consumption relationship to describe active stress development and ATP consumption of the left ventricle model, the PVA-ATP consumption relationship was predicted theoretically from the SSA-ATP consumption relationship in our simulations.

Transmural distribution of ATP consumption can be estimated from the measurements of high-energy phosphates in the cardiac wall. The PCr-to-ATP ratio measured by nuclear magnetic resonance is slightly higher in the epicardial layer than in the endocardial layer and with midwall layer value between these two (12, 39). The value of the PCr-to-ATP ratio in the endocardial layer is ~85% of the PCr-to-ATP ratio in the epicardial layer in control conditions. The level of inorganic phosphate at these conditions was too low to be detected, in correlation with the theoretical studies of cardiac intracellular energy transfer (37). Because a cardiac cell is metabolically stable (the PCr-to-ATP ratio changes slowly when the workload is changed) (26), the errors in PCr-to-ATP estimation lead to very large errors in estimation of the oxygen consumption. However, from the available PCr-to-ATP ratio measurements, one can conclude that oxygen consumption in epicardial layers is lower than in endocardial layers. According to our simulations, the region with the highest ATP consumption is in the middle of the left ventricular wall, slightly shifted toward subendocardial region (Fig. 13). The reason of the discrepancy between the computed and measured data is not clear and requires further investigation. The computed ATP consumption distribution is distorted in the apex and basal regions, which may be caused by the stress concentration close to the boundaries during the simulations. The question as to whether such a distortion of ATP consumption distribution takes place in vivo is still open.

In conclusion, we have shown in this study that there exists a local minimum of the sarcomere strain and stress variances in the region that corresponds to high ejection fraction and high efficiency of the left ventricle. If ATP consumption variance was used to find the fiber orientation angles, then the transverse fiber angle was underestimated. In addition, we have shown that the variances of sarcomere strain, developed stress, and ATP consumption are minimized by almost the same transmural course of the helix fiber close to the measured one.