## Abstract

The aim of this study was to investigate the influence of fiber orientation in the left ventricular (LV) wall on the ejection fraction, efficiency, and heterogeneity of the distributions of developed fiber stress, strain and ATP consumption. A finite element model of LV mechanics was used with active properties of the cardiac muscle described by the Huxley-type cross-bridge model. The computed variances of sarcomere length (SL_{var}), developed stress (DS_{var}), and ATP consumption (ATP_{var}) have several minima at different transmural courses of helix fiber angle. We identified only one region in the used design space with high ejection fraction, high efficiency of the LV and relatively small SL_{var}, DS_{var}, and ATP_{var}. This region corresponds to the physiological distribution of the helix fiber angle in the LV wall. Transmural fiber angle can be predicted by minimizing SL_{var} and DS_{var}, but not ATP_{var}. If ATP_{var}was minimized, then the transverse fiber angle was considerably underestimated. The results suggest that ATP consumption distribution is not regulating the fiber orientation in the heart.

- heart
- oxygen consumption
- pressure-volume area
- mathematical modeling

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

#### 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.

In describing the left ventricular geometry, ellipsoid coordinates (ξ, θ, φ) are used, which are related to Cartesian coordinates (*x*, *y*, *z*) according to
Equation 1
Equation 2
Equation 3Surfaces of constant radial ellipsoid coordinate ξ correspond to ellipsoids. Within the wall, ξ ranges from 0.37 to 0.68 at endocardial and epicardial surfaces, respectively. The longitudinal ellipsoid coordinate θ ranges from [3π/10] at the base to π at the apex. In addition, local normalized transmural
and longitudinal
coordinates are defined. These coordinates vary linearly with the distances in the wall:
= −1 at the endocardial surface,
= 1 at the epicardial surface,
= −1 at the apex,
= 0 at the equatorial plane, and
= 0.5 at the base.

The fiber orientation is quantified by two angles: the helix fiber angle (α_{h}) and the transverse angle (α_{t}) (Fig. 1). Here, α_{h} is defined as the angle between the φ-direction and the projection of fiber path on the (θ,φ)-plane. α_{t} is then an angle between the φ-direction and the projection of fiber path on the (ξ, φ)-plane. To keep the number of parameters as small as possible, α_{h} is approximated as a linear function of
Equation 4Because fibers do not end at the endocardial and epicardial surfaces, α_{t} is 0 at these surfaces
Equation 5The parameters *P*
_{1},*P*
_{2}, and *P*
_{3} were varied in this study.

#### 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 (−p**I**) of the fluid trapped in the solid
Equation 6where **I** is the unity tensor.

The myocardial tissue is simulated as a hardly compressible material with hydrostatic pressure given by
Equation 7where **F** is the deformation gradient tensor,*K*
_{b} is the bulk modulus, and det is determinant. The value of *K*
_{b} has been set at 50 kPa to keep volume changes within the ±5% range.

The constitutive properties of the passive tissue are modeled transversely isotropic, with stress increasing exponentially with strain (38). The ς_{p} is determined by a strain-energy function [W(**E**)] that relates the second Piola-Kirchoff stress tensor (**S**) to the Green-Lagrange strain tensor (**E**)
Equation 8The strain-energy function is adapted from Bovendeerd et al. (2)
Equation 9where *a _{i}
* (

*i*= 0, 1, 2, 3) are material parameters and Equation 10 Equation 11The ς

_{p}is found from

**S**according to Equation 12where the index

*c*denotes a tensor conjugate. To overcome convergence problems in apical and basal regions, the tissue has been stiffened by multiplication of

*a*in

_{i}*Eq. 9*with factor β shown in Fig.2 for normal and stiffened models.

The ς_{a} is computed by a Huxley-type cross-bridge model (36). The cross-bridge model is able to reproduce the following experiments performed on the cardiac muscle:*1*) active stress dependency on time and sarcomere length in isometric contraction (15), *2*) sarcomere shortening velocity as a function of afterload (35),*3*) end-systolic relationship in isometric and isotonic experiments (14), and *4*) the linear dependency of oxygen consumption on the stress-strain area (14). The active stress and ATP consumption are computed as functions of time if sarcomere length and sarcomere shortening dynamics are given. In the model, the sarcomere length and its shortening velocity are found with the use of tensor **F** and **e**
_{f}. The complete description of equations and parameters is given in Ref.36.

#### 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 Equation 13The 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).

The hemodynamic coupling of the left ventricle to the aorta is described by a three-element aortic input impedance (2). Axial displacement of the nodes in the basal surface and circumferential displacement of subepicardial basal ring are suppressed. A uniform left ventricular pressure is applied to the entire endocardial surface. Epicardial pressure is assumed to be zero during a cardiac cycle.

#### 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
Equation 14where **
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 Equation 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*

_{3}= 5° over a large (

*P*

_{1},

*P*

_{2}) design space:

*P*

_{1}and

*P*

_{2}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*/2), var(ς

_{s}_{a}), and var(

*V*) to the changes in the mechanical activation sequence.

Near the minima corresponding to the highest ejection fraction of the left ventricle, a more detailed optimization was performed in the left ventricle with normal stiffness and simultaneous activation. To find the optimal *P*
_{1} and *P*
_{2}values as a function of *P*
_{3}, the following procedure was used. First, the variance of stress, strain, or ATP consumption was computed at (*P*
_{1},*P*
_{2}) values, which corresponded to the nodal coordinates of the mesh (7 by 6 elements) in (*P*
_{1}, *P*
_{2}) space with the following corners: (15°, −72°), (27°, −72°), (27°, −62°), and (15°, −62°). The smallest value of the variance was then identified and the mesh was refined around the corresponding node. The procedure was repeated three times. The value of*P*
_{1} and *P*
_{2} with the smallest variance was then recorded together with the corresponding value of the variance. The distance between the neighbor nodes in the final mesh was recorded as *P*
_{1} and*P*
_{2} estimation error. It turned out that only one local minimum of the variances was found in our simulations.

The model was tested by comparing the deformation of the left ventricle to the experimental data. To test the computed*V*
distribution in the left ventricle, we compared the computed distribution with the phosphocreatine (PCr)-to-ATP ratio measured by nuclear magnetic resonance (11,39). In addition, total ATP consumption by the left ventricle was computed and related to the pressure-volume area (PVA).

## RESULTS

#### 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*
_{1} and *P*
_{2} represent the α_{h} in the midwall (*P*
_{1}) and the slope (*P*
_{2}) of transmural α_{h} changes. There are several local minima of the variances in the (*P*
_{1},*P*
_{2}) plane. In this simulation, the smallest variance of sarcomere length was found at *P*
_{1} = 20° and*P*
_{2} = −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).

In the following simulations, we focused on the behavior of local minima that correspond to the high ejection fraction of the left ventricle. We traced the position of the minima in the (*P*
_{1},*P*
_{2}) plane at different *P*
_{3} angles (Fig.5). The values of the corresponding variances are shown as functions of *P*
_{3} angles in Fig. 6. The smallest variances of sarcomere length, ς_{a}, and ATP consumption were at*P*
_{3} angles of 5°, 10°, and 0°, respectively. The absolute values of the smallest variances were as follows: var(*l*
_{s}/2) = 0.015 μm, var(ς_{a}) = 3.7 kPa, and var(*V*
) = 0.22 ATP molecules per myosin head per beat. It is important to note that the position of the local minima of the variances in the (*P*
_{1},*P*
_{2}) plane were close to the position determined by stiffened model (Fig. 3), indicating relatively small sensitivity of the optimal (*P*
_{1},*P*
_{2}) to the increase of the stiffness. However, the decrease in passive stiffness resulted in an increase of the ejection fraction from 32% to 48%. Ventricular ejection fraction and efficiency on α_{t} are both maximized at a *P*
_{3} value of 16° (Fig.7).

The comparison of the found optimal fiber orientations with the available experimental data is presented in Figs.8 and 9. Taking into account the linear approximation of the transmural distribution of α_{h} used in our model, the measured and predicted α_{h} are relatively close. Distribution of transverse angle α_{t} predicted by the model is in the error range at the left ventricle region between apex and equator if strain or active stress distributions are optimized. At the base and apex region, α_{t} is underestimated. The angle α_{t} predicted by optimization of ATP consumption distribution did not correspond with the measured data (see Fig. 8).

#### Model testing.

In the test simulations, we used the fiber orientation angles predicted by minimizing sarcomere length variance: *P*
_{1}= 22.5°, *P*
_{2} = −69°, and*P*
_{3} = 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.

Torsion of the left ventricle is determined by the balance of epicardial and endocardial fibers. Thus it should be sensitive to the fiber orientation angles α_{h} and α_{t}. Torsion of the apex with respect to the base, computed for simultaneous and nonsimultaneous activation, is demonstrated on the apex rotation angle-pressure relationship (see Fig.10). The computed apex rotation-pressure loop proceeds in the same direction as found experimentally. Apex rotation was sensitive not only to the changes in fiber orientation, but to the mechanical activation timing too. In case of nonsimultaneous activation (we used the activation that was three times faster than the electrical activation), the computed apex rotation angle-pressure loop was wider than the measured one. In case of simultaneous activation, the computed span is less than the measured one.

It is possible to estimate the efficiency of the left ventricle from the amount of consumed ATP molecules per myosin head and mechanical work performed by the ventricle. Taking into account the myosin ATPase concentration (34) of 0.18 mM or 0.18 mol/m^{3}, free energy change during ATP hydrolysis of 60 kJ mol^{−1}(11, 30), and assuming that the efficiency of the oxidative phosphorylation to the free energy change of ATP hydrolysis is equal to 60–70% (30), the mechanical efficiency (external work/total V˙o
_{2}) of the left ventricle computed by the model was equal to 21–24%. According to the data reviewed by Suga (30), the mechanical efficiency is 10–30%, in accordance with our simulation.

Pressure-volume relationship computed by the model for various diastolic filling pressures is summarized in Fig.11. This relationship was used to find the PVA and to relate the PVA to total ATP consumption by the ventricle (Fig. 12). To compute the PVA, we used the same procedure as by Suga et al. (31). First, the volume-axis intercept of the end-systolic pressure-volume relationship was found by using end-systolic points of isovolumetric and by ejecting contractions (Fig. 11). A straight line was then drawn between the found intercept point and end-systolic point on a specific pressure-volume trajectory. The area between this straight line, end-diastolic pressure-volume relationship and the systolic segment of pressure-volume loop is the PVA.

The ATP consumption-PVA relationship obtained for isovolumetric contractions was fitted by line (Fig. 12)
By expressing ATP consumption and PVA in joules and by using amount of the myosin head concentration together with the free energy change during ATP hydrolysis, the relationship can be transformed to*V*
(J/beat) = 1.71 · PVA(J/beat) + 0.17. In our simulations, the contractile efficiency (not mechanical efficiency computed above), defined as reciprocal of the slope, was 58%. Taking into account the efficiency of the oxidative phosphorylation (see above), the contractile efficiency related to the oxygen consumption of the left ventricle was 35–41%, in good correlation with the data reviewed by Suga (30).

According to our simulations, ATP consumption within the wall was not homogeneous but slightly higher at the midwall region (Fig.13).

## DISCUSSION

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*
_{1},*P*
_{2}) 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*
_{1},*P*
_{2}) 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*
_{1},*P*
_{2}) 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*
_{3} is varied. Namely, the amount of consumed ATP is maximal at*P*
_{3} = 20° (not shown) and is ∼15% smaller at *P*
_{3} = −10° when only*P*
_{3} is varied at fixed *P*
_{1}and *P*
_{2} values close to the optimal (*P*
_{1},*P*
_{2}), 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*
_{3}= 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*
_{1},*P*
_{2}) when either *P*
_{1} or *P*
_{2} 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.

## Acknowledgments

This work was supported in part by Estonian Science Foundation Grant 4704.

## Footnotes

Address for reprint requests and other correspondence: M. Vendelin, Institute of Cybernetics, Akadeemia 21, 12618 Tallinn, Estonia (E-mail: markov{at}ioc.ee).

The costs of publication of this article were defrayed in part by the payment of page charges. The article must therefore be hereby marked “

*advertisement*” in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.May 9, 2002;10.1152/ajpheart.00874.2001

- Copyright © 2002 the American Physiological Society