Abstract
In the present paper, the study of the ventricular motion during systole was addressed by means of a computational model of ventricular ejection. In particular, the implications of ventricular motion on blood acceleration and velocity measurements at the valvular plane (VP) were evaluated. An algorithm was developed to assess the force exchange between the ventricle and the surrounding tissue, i.e., the inflow and outflow vessels of the heart. The algorithm, based on the momentum equation for a transitory flowing system, was used in a fluidstructure model of the ventricle that includes the contractile behavior of the fibers and the viscous and inertial forces of the intraventricular fluid. The model calculates the ventricular center of mass motion, the VP motion, and intraventricular pressure gradients. Results indicate that the motion of the ventricle affects the noninvasive estimation of the transvalvular pressure gradient using Doppler ultrasound. The VP motion can lead to an underestimation equal to 12.4 ± 6.6%.
 fluidstructure interaction
 Doppler ultrasound
 pressure gradients
 blood acceleration
during cardiac ejection, biochemical energy is converted into mechanical energy by muscle fibers and transferred to blood. The blood is then delivered to the systemic circulation. The mass and momentum released with the outflowing blood (stroke volume) are responsible for ventricular motion in the opposite direction, similar to the recoil of a gun. This motion is restrained by the inflow and outflow vessels, which constrain the ventricle. To the best of our knowledge, few studies have been performed concerning this phenomenon. The center of mass (CoM) motion has been investigated with reference to the total heart (9) and the singleventricle heart (7). These studies demonstrated that, during systole, CoM motion is less than 3 mm and is directed toward the ventricle base. Cao and colleagues (3a) have calculated a maximal excursion of 10–20 mm for the valvular plane (VP) during the cardiac cycle.
To determine the influence of ventricular motion on ventricular mechanics, an algorithm based on the momentum equation for a transitory flowing system was developed. The algorithm was integrated into an axisymmetric fluidstructure model of the ventricle (4,19) to calculate the CoM motion during ejection under normal conditions.
Different elastic and viscous constraints were applied to the VP. The elastic constraint was varied to modify the extent of the maximum displacement. The viscous constraint was introduced to damp the oscillations due to a purely elastic constraint. The model was used to evaluate the effect of ventricular motion on the intraventricular pressure gradients, which were demonstrated in a previous work to be more sensitive to changes in ventricular contraction than the maximum ventricular pressure and the first time derivative of ventricular pressure (dP/dt) (20); furthermore, it was used to assess the noninvasive Doppler estimation of ventricular performance indexes, such as the transvalvular pressure drop and dP/dt.
MATERIALS AND METHODS
CoM motion algorithm.
The momentum equation for the ventricular CoM is written in scalar terms according to the vector assignment of Fig.1. The relative velocity of bloodw
_{VP} at the outlet is
The ventricle is defined by a control volume limited by the outer surface of the ventricular wall and by the outflow area at VP, and its mass varies with time. During systole, mass outflows through the VP;M is the mass of the ventricle, including blood and wall, and its time derivative divided by the density ρ is the volumetric flow entering the ventricle (M˙), which is the opposite of the volumetric flow leaving the ventricle
The overall momentum variation for a constant mass system equals the external forces. In this case, where a varying mass system is taken into account, the momentum inside the control volume changes also according to the momentum acquired (lost) throughout the mass entering (leaving) at the VP; the additional term to be included in the momentum equation is the unit volume momentum of the blood crossing the VP, ρv
_{VP}, multiplied for the volumetric transfer rate across the VP taken positive for entering flows (−Aw
_{VP}). The term to be added is hence
The first right term is responsible for the backward motion due to blood ejection.
Equation 4
can be modified by applying the chain rule to the time derivative of the lefthand term. When the mass derivative terms in the righthand side are grouped and the elastic and viscous forces are explicitly stated, it becomes
Finally, by substituting Eq. 1
in the first right term, the following equation is obtained
By combining Eqs. 2
and
6
, the following equation is obtained
Fluidstructure interaction model.
Equation 6 was integrated in an axisymmetric model of the left ventricle described in detail in previous publications (19, 20) and depicted in Fig.2. The model accounts for fluidstructure interaction between the intraventricular fluid and the myocardial wall.
In the reference state, the geometry of the left ventricle is assumed to be a prolate ellipsoid with a minor semiaxis equal to 25 mm, a major semiaxis equal to 53 mm, and an inner volume of 113 cm^{3}according to normal human values (5, 23). The ellipsoid is truncated at a distance of 79.9 mm from the ventricular apex. The ventricle wall mass is 300 g with a uniform wall thickness. The heart rate is assumed to be 72 beats/min.
The intraventricular blood behavior is simulated by means of the FIDAP generalpurpose fluid dynamics code (Fluent, Lebanon, NH). The main assumptions are that the fluid is incompressible, homogeneous, and Newtonian (ρ = 1.06 × 10^{3}kg/m^{3}, μ = 3 × 10^{−3} kg · m^{−1} · s), the walls are impermeable, and no slip occurs at the wall. Gravitational effects are ignored because a supine position is assumed. The fluiddynamics mesh consists of 1,251 nodes; ninenode (quadrilateral) isoparametric annular elements are used to discretize the fluid domain. The full Navier Stokes equation (including the transient term) and the continuity equation in their axisymmetric formulation are solved by using Galerkin's weighted residual approach in conjunction with finite element approximation. To account for myocardial inner surface movement, an additional free surface degree of freedom is taken into account in accordance with Saito and Scriven's method (21). This method requires the assignment of displacement and velocity for each moving boundary node through time functions. A quasi Newton solver was used for the solution method. The timeintegration technique adopted was the implicit backward Euler with a fixed time step of 1 ms.
The wall stress and motion were calculated by means of a structural model, constructed on the basis of the modified version of Wong's sarcomere model (27) proposed by Montevecchi and Pietrabissa (12) and Pietrabissa et al. (16). A detailed description of the sarcomere model is given in the . Sarcomeres are arranged to constitute a thinshell model composed of two sets of counterrotating fibers.
When gravitational effects and local wall inertia and viscosity are ignored, the equation of equilibrium is
The fiber thin shell is discretized into five threenode axisymmetric elements, with a total number of nodes equal to 11. Nodal forces and displacements are calculated by integrating Eq. 7 throughout each element, approximating the strain and stress fields quadratically. The boundary pressure field is also approximated quadratically. A solver based on a multiplane quasiNewton method proposed by Buzzi Ferraris and Tronconi (3) was used for the solution of the wall motion.
To calculate the aortic impedance, a lumped parameter model of the systemic circulation was used and integrated with the structural model. It supplies the aortic pressure at the outlet section (aortic valve plane) on the basis of the ventricular outflow. The lumped parameter values adopted are reported in a paper by Avanzolini et al. (1). A fourthorder RungeKutta method was used to solve the model with a fixed time step equal to 1 ms; the stability of the RungeKutta algorithm was evaluated by comparing the results obtained with different time steps (range 0.1–10 ms). Four cycles were performed to allow the lumped parameter model to reach the proper initial conditions.
Interface and solution procedures.
The calculation scheme is based on an uncoupled approach (4). On the basis of a tentative intraventricular pressure distribution on the moving boundary, the fiber shortening is calculated at the 11 nodes of the thin shell model to satisfy equilibrium: the fiber shortens so that, at these nodes, the force generated by the wall balances the force generated by the intraventricular fluid forces. The fluid dynamic stress tensor (ψ_{ij}), calculated by FIDAP at each node of the boundary, has the following form
The wall motion and u
_{VP} are the boundary conditions (Fig. 2) for the fluid dynamic model, which provides a new intraventricular pressure distribution along the boundary (P_{FIDAPi}; i = 1, … , number of nodes of the structural model). The latter is compared with the tentative one (P_{i}). If the difference is within the tolerance limit (0.1 mmHg), the time is incremented and a new time step is performed. Otherwise the tentative values for stress distribution are recomputed as
Computations were performed with an IBM RISC/6000 model 230 workstation. Computational time for each iteration was 105 s, for a total number of iterations per time step of ∼10 iterations (800 ÷ 1,200 s per step). The adopted time step was 1 ms with a total number of ∼140 time steps to simulate the ejection phase.
RESULTS
Two simulation sets were performed to investigate the role of the elastic and viscous constraints, respectively. In the first set, two simulations were performed with two different values for the elastic constraint k (200 and 800 kg/s^{2}) and without the viscous constraint η (equal to 0 kg/s). As control cases, two simulations have also been performed by leaving the VP free and maintaining VP fixed, respectively. In the second set, the influence of the viscous constraint was studied by comparing the results of two simulations performed with two different values of η (5 and 10 kg/s) and leaving k fixed equal to 200 kg/s^{2}. The simulation of the previous set with k equal to 200 kg/s^{2} and η equal to 0 kg/s was adopted as the control case for this set.
The values for k were chosen to obtain maximum VP displacements of ∼8 mm in the case of k = 200 kg/s^{2} and 2 mm in the case of k = 800 kg/s^{2} (3a). The values for η were chosen iteratively to halve the displacements of VP and CoM.
The overall features of the ventricular performance did not vary significantly throughout the simulations: the ventricular stroke was ∼58.4 ± 0.3 cm^{3} (mean ± SD) with an ejection fraction equal to 51.3 ± 0.2%. The pressure at the VP (P_{0}) varied in the range of 77 to 133.6 ± 0.6 mmHg. In turn, the intraventricular pressure gradients and VP and CoM motion varied appreciably as discussed below.
In Fig. 3, intraventricular pressure difference (IPD) time courses are given. They are calculated in four points (16, 32, and 48 mm from the VP and in the ventricular apex) as the difference between the local pressure (P_{i}) and P_{0}. IPD is shown because it was demonstrated in a previous paper that it is the most sensitive index of ventricular contraction (20). Results are reported with reference to the two elastic constraints (200 and 800 kg/s^{2}) and with free and fixed VP. IPDs show an oscillating pattern that fades during ejection, with the highest peak occurring in the earlier phase when blood acceleration is maximum. As previously reported (19,20), the computed IPD peak values at the ventricular apex are similar to IPD values observed in vivo (6,15, 18). The oscillating pattern is more pronounced than in in vivo observations (14,15) because of the viscous characteristic of the wall, which has not been accounted for in the present study (19). The IPDs in the four panels are similar and scaled moving from the ventricular apex toward the base. The IPD patterns are driven by ventricular contraction and VP and CoM motion. Their interpretation is not trivial; nevertheless, some considerations may be argued about the effects of k on IPD when the acceleration phase (the first 50 ms of the ejection phase) is considered. In this phase, two peaks occur. The effect of k is evident on the second peak with reference to the IPD time courses of the ventricles with free VP, k equal to 200 kg/s^{2}, andk equal to 800 kg/s^{2}; a lower kpromotes the motion of the ventricle in the direction opposite to the flow, thus facilitating the ventricular ejection and decreasing the value of the IPD peak. During the first peak, this occurrence is not detectable because the VP motion is still small and the forces due to the constraints are negligible. In turn, when the VP motion is impeded, a higher IPD (first) peak occurs along with a delay of the IPD oscillation. The difficulties encountered by the ventricle with fixed VP in the initial phase are detectable also by another occurrence: the comparison of an IPD peak at 30 ms close to the VP (Fig. 3,C and D). This occurrence is due to the fact that the wall contraction occurs mainly in the region close to the ventricular base where blood inertia is lower (19).
Figure 4 shows the displacements and the velocities of CoM and VP, respectively. At free VP, both CoM and VP displace increasingly throughout the ejection phase. In the two cases with elastic constraint, the motion of VP is initially negative. The motion of CoM is positive and goes toward VP because of the decrease of the ventricular volume. In middle/late systole, the behavior is different, depending on the value of k. The higher k is, the earlier the reversion of motion.
In the second simulation set, the effects of a viscous constraint were investigated. Results obtained with two η values of 5 and 10 kg/s, respectively, were compared with those obtained without viscous constraint. The value of k was fixed equal to 200 kg/s^{2}. Also in this simulation set, the overall features of the ventricular performance did not vary significantly throughout the simulations and from the previous set. The effect of the viscous constraint on IPDs is reported in Fig. 5. As expected, a value of η different from 0 decreases the oscillations of IPD according to a damping behavior. This occurrence is discernible with reference to the two negative peaks in the accelerative phase and the late systole. Also CoM and VP displacements and velocities in the direction opposite to blood outflow decrease when the η value is increased according to a damping behavior as reported in Fig.6.
DISCUSSION
In the present study, a model, described in detail previously (4, 19), was used to test the motion of the ventricular CoM. For this purpose, an algorithm was implemented to take into account the balance of momentum with external k and η constraints. Different values of k and η were used to simulate the VP motion within the physiological range.
IPD and VP and CoM motions in early systole.
During the first 20 ms of ejection, no changes can be noted when the mechanical features of the constraints are varied. This can be attributed to the fact that the elastic and viscoelastic constraints are not significantly loaded. A different behavior occurs when VP is fixed; in this case higher values of the first peak were calculated at the apex. In the other simulations, VP shows a little movement toward the apex during early ejection; this movement facilitates the ejection and reduces the initial pressure peak significantly. A fixed VP requires the generation of a higher fiber force and evokes a different dynamic behavior throughout ejection by increasing both the first IPD peak and the period of IPD oscillation.
IPD and VP and CoM motions in late systole.
In the subsequent part of the ejection phase, the time courses of VP and CoM motions and of IPD differ significantly; this allows evaluation of the effects of the constraints. For instance, a high value fork decreases VP and CoM motion in the direction opposite to blood outflow and increases the oscillating frequency ofs _{VP} (Fig. 4). Furthermore, the second positive peak of IPD increases markedly (Fig. 3). This phenomenon, however, has not been observed in in vivo recordings (14,15), thus suggesting an inhibition either by a viscoelastic constraint or by a nonlinear behavior of the elastic constraint. Indeed, the values of η used in the present study decrease the magnitude of IPD oscillations (Fig. 5). In general, the IPD course and CoM and VP motions are damped by increasing the value of η. Furthermore they become positive sooner (Figs. 4 and 5). However, this is not the case for the first IPD peak in the early ejection phase where the velocities are still too low to generate a substantial viscous force.
Doppler velocimetry implications.
There are at least two clinical implications of the present results. They concern the Doppler estimation of the acceleration and velocity of the outflowing blood used for the noninvasive estimation of ventricular performance. The first implication is the estimation of the maximum dP/dt (dP/dt _{max}) (2,24) through the measurement of the maximum aortic acceleration by using the formula dP/dt _{max}= ρc dw/dt _{max}, where ρ is the density and c is the velocity of the pulse wave; second is the estimation of the transvalvular pressure gradients by means of the simplified Bernoulli formula (10,13, 28) IPD_{trans} = 4w _{max} ^{2} where IPD_{trans} (mmHg) is the estimated transvalvular gradient and w _{max}(cm/s) is the maximum outflow blood velocity. Doppler velocimetry measures the absolute velocity and acceleration in a selected frame, usually located close to VP, thus accounting for both blood and VP velocities and accelerations. When the two velocities (accelerations) have opposite directions, Doppler ultrasound will underestimate the blood outflow velocity (acceleration); it will overestimate when both are in the same direction. This inaccuracy is not known a priori because VP motion depends on the elastic features of the fibers and the constraints. In the case of blood acceleration, the peak occurs at the very beginning of the ejection phase concomitant with the VP acceleration peak that has an opposite direction; however, simulation results show that VP acceleration is ∼1.2 ± 0.4% (mean ± SD) of the blood outflow acceleration; hence it does not significantly affect dP/dt estimation. In turn, as far as the transvalvular pressure gradient estimation is concerned,w _{max} occurs ∼40 ms after the valve opening. Also in this case, the outflow and VP velocity have opposite directions because the ventricular CoM moves backward (as does the VP); this would imply that Doppler ultrasound estimates a value for the velocity peak that is lower than the blood outflow value. Simulations show a difference between the blood outflow velocityw _{max} and the absolute velocityv _{max} equal to +6.2 ± 3.3% with respect to blood outflow velocity, depending on the constraints. This leads to an underestimation of transvalvular pressure gradients equal to 12.4 ± 6.6%, compared with recordings in which solidstate pressure sensors are used.
Conclusions.
This paper describes the physics of a ventriclelike model. The model accounts for bloodventricular wall interaction by assuming axisymmetry and the absence of fiberdamping behavior. As discussed previously (19, 20), the model is not able to illustrate local fiber impairment and flow field asymmetry and overestimates the oscillatory pattern of IPD in late ejection. In the present study, the model was modified to simulate the ventricleexternal structure interaction, with the assumption that the vessels connected to the ventricle have a linear viscoelastic behavior. Results show that 1) ventricular motion aids early ejection blood outflow; 2) the first intraventricular pressure peak is not affected by VP motion except when VP motion is neglected;3) the intraventricular pressure gradient time course is, in turn, sensitive to ventricular motion; 4) the ventricular motion does not affect ejection in terms of overall performance; and5) VP velocity is comparable to blood outflow velocity and can give inaccurate Doppler estimation of transvalvular pressure gradients, whereas VP acceleration is negligible with respect to blood acceleration.
Footnotes

Address for reprint requests and other correspondence: A. Redaelli, Dipartimento di Bioingegneria, Politecnico di Milano, P.za L. da Vinci, 32, 20133 Milano, Italy (Email:redaelli{at}biomed.polimi.it).

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. §1734 solely to indicate this fact.
 Copyright © 2000 the American Physiological Society
Appendix
According to Wong's model (27), the behavior of the sarcomere is defined by the sum of two terms representing the behavior of two nonlinear elastic springs, the former in series (SE) with the contractile unit (CE) and the latter in parallel (PE) with the complex SECE. The first term accounts for the passive elastic characteristic of the myofiber (PE); the second term differs from zero during contraction and accounts for the active elastic response (SE)
The Δl _{SE} is defined as the sum of the sarcomere shortening (l_{s} −l_{s} _{0}) and the contractile unit shortening Δl _{CE} and is normalized with respect to the difference between the sarcomere length at which maximum tension is generated and its resting length (12).
CE contracts with time according to a modified version of Julian's (11) model of skeletal muscle activation. This curve is parametric with respect to the frequency and the maximum shortening of CE (Δl _{CE max}).
To calculate the value of Δl
_{CE max}, the following two steps are necessary. First, the value of F_{CE max} is computed: the relationship between the normalized maximum tension F_{CE max} and the preload (Δl
_{CE})_{t=0} is calculated according to the method of Pollack and Krueger (17). Then the relationship between F_{CE max} and Δl
_{CE max} is computed from Eq. EA1
in the hypothesis of isometric contraction. According to this hypothesis, the preload is
In the calculations, the effective normalized stress defined as ς_{eff} = φ · ς is used instead of ς. φ Is the fiber number that varies with respect to the section considered. In the reference condition, e _{f} and φ are calculated for the 11 mesh nodes under the assumption of an intraventricular pressure of 7 mmHg and a fiber preload length equal to the length that gives the maximum isometric force (27,29).