A finite element model of the infarcted left ventricle

A finite element model of the infarcted left ventricle

0021 -Y?W c 83 oloo4c I4 $03 00 0 1983 Perpamon Press Lid A FINITE ELEMENT MODEL OF THE INFARCTED LEFT VENTRICLE ALAN NEEDLEMAN Division of Engine...

1MB Sizes 18 Downloads 83 Views

0021 -Y?W c

83 oloo4c

I4 $03 00 0

1983 Perpamon Press Lid

A FINITE ELEMENT MODEL OF THE INFARCTED LEFT VENTRICLE ALAN NEEDLEMAN Division of Engineering, Brown

University,

Providence.

RI 02912. U.S.A.

and STUART A. RABINOWITZ,

Division

of Applied

DANIEL

K. BOGEN and THOMASA. MCMAHON

Sciences. Harvard

University,

Cambridge.

MA 02138. U.S.A.

Abstract-A membrane theory model of the infarcted left ventricle is described which represents the pressure-volume relation of the ventricle in diastole and end-systole, states in which the contraction velocity is zero. This membrane theory model is shown to give essentially the same diastolic pressure-volume relation as a more elaborate model which fully accounts for variations through the wall thickness. The end-systolic pressure-volume relation of an infarcted ventricle is obtained by means of a specially tailored finite element method. which is described in detail. From the end-systolic and diastolic pressure-volume relations the effect of infarction on the pumping performance of the ventricle can be determined. Representative numerical results are presented.

INTRODlJCTlOh’

Fig. 1. Starting at the end-diastolic point A, the pressure rises in isovolumic systole to the onset of

We develop a ventricular model which focuses on a particular mechanical consequence of myocardial infarction: the effect of an infarct on the pressure-volume relation of the ventricle. For the present purposes, the cardiac cycle is modelled as shown in

ejection at point B. We employ the idealization that ejection from the ventricle takes place at a constant pressure (between points Band C in Fig. 1). Isovolumic relaxation proceeds between C and D, and diastolic (passive) filling from D to A. The point A is determined by the diastolic filling curve and the filling pressure (10 mm Hg) in Fig. 1, where, as throughout this paper, all pressures are transmural pressures. The point C is determined by the end-systolic pressure-volume curve

Rrcriredfor

publicamn

13 .4uyusr

1982

--END-SYSTOLIC

1

i

_._-;s

DlA;LlC

20 NORMhLlZED

22

]

2.4

VOLUME

Ftg. 1. Diastolic and end-systolic pressure-volume curves for a normal ventricle. calculated from equations (17) and (19) of the text with Lp = 16, pp = 2. to/R, = 0.20 and (Zp~r,jR,) = 246. The stroke volume, SV,, is shown for the cardiac cycle ABCD with an end-diastolic pressure of 10 mmHg and an end-systolic pressure of 100 mmHg. Stroke volume. SV,. corresponds to the cycle AR’CD. with an end-diastolic pressure of 20 mmHg.

45

46

ALAN NEEDLEMAN, STUART A. RAEIINOWITZ,DANIEL K. B~GEN and THOMAS A. MCMAHON

and the outflow pressure (100 mm Hg in Fig. 1). Thus, given the diastolic and the end-systolic pressure-volume curves for a ventricle, various measures of cardiac performance (e.g. stroke volume) can be computed as functions of the filling pressure and outflow pressure as in Bogen et al. (1980), assuming the peak isovolumic and end-ejection pressure-volume relations coincide. Experimental evidence (e.g. Monroe and French, 1961; Suga and co-workers, 1971, 1973, 1974, 1977) indicates that in the normal canine left ventricle, for a given contractile state of the myocardium, there is essentially a unique end-systolic pressure-volume relation independent of the loading conditions on the ventricle. A recent thorough review of the experimental evidence relevant to this issue is given by Sagawa (1978). Thus, for a normal ventricle, myocardial activation and rates of contraction are sufficiently rapid to allow full ejection within the time limits imposed by the finite activation time. We will, in addition, assume that there is a unique end-systolic pressure-volume relation, independent of the loading conditions, for the infarcted ventricle. Violation of this assumption for the infarcted ventricle, implying less than full activation of the residual myocardium, would tend to degrade mechanical performance beyond what is predicted by the present model. We model the ventricle as an axisymmetric membrane, with the diastolic reference configuration taken to be a spherical membrane of uniform thickness. The infarcted myocardium is represented as a region of altered mechanical properties. The diastolic pressure-volume curve is determined numerically by passively inflating the membrane model of the ventricle. At various stages of inflation, the ventricle is allowed to contract isovolumically by changing the constitutive relation of the non-infarcted myocardium to its systolic form, while constraining the volume enclosed by the ventricle to remain constant. In this manner, points on the end-systolic pressure-volume curve can be determined. The appropriateness of this simple membrane model. for the present purposes. is justified by comparing the diastolic pressure-volume relation obtained from membrane theory with the one obtained by Janz and Waldron (1978), employing a model of realistic geometry and fully accounting for variations through the wall thickness. A description is given of the finite element numerical method employed to solve the governing equations, with particular emphasis on the approach used in the systolic phase of the analysis. Finally. some representative numerical results are presented.

THE MEMBRANE

MODEL OF THE VENTRICLE

The model of the infarcted ventricle employed by Bogen et a/. (1980) is depicted in Fig. 2, with Fig. 2(a) illustrating the spherical unstressed reference configur-

ation and Fig. 2(b) illustrating a representative deformed configuration. In each figure, the unshaded region is normal myocardium while the upper, shaded. cap-shaped region represents the infarct. The formulation of the membrane theory equations employed here is the specialization of the general nonlinear membrane equations given by Budiansky (1968) to initially spherical membranes. Material points on the surface of the sphere are labelled by the coordinates 0 and 4, where 6,is the polar angle and 4 is the longitudinal angle (see Fig. 3). Since attention is restricted to axisymmetric deformations of the membrane, all field quantities are independent of 4 and the only nonvanishing displacement components are w, normal to the undeformed sphere,and u, tangent to the undeformed sphere in a meridian plane, as illustrated in Fig. 3. Principal axis notation is employed with the following conventions adopted; subscripts range from 1 to 2, subscript 1 denoting principal values associated with the e-direction, subscript 2 denoting corresponding values in the &direction and ( )rgdenoting differentiation with respect to 6. The deformation of an initially spherical membrane is specified by the following three quantities, e,

=+(w+u,,)r,=$-(w+ocotfI) 0

$

I##=

0

().t’,@- u)

(1)

0

where e, and e2 are the principal values of the linear strain tensor and + is a measure of the rotation of the normal to the membrane. The principal stretches, i.i, and the principal values of the Lagrangian strain tensor. rli, are given in terms of e,, e2 and $ by, i., = [(l +e,)2+$2]1’2

i., = 1 +e,

(2)

p/r = er ++e: + +$’

t!2 = e2 ++e:.

(3)

The materials considered here are taken to be incompressible and the incompressibility condition is expressed by the relation 1 A’=t,=i.,12

1 (4)

where t is the current thickness. The principle of virtual work. on which the numerical methods to be described subsequently are based, takes the form, I

27tR; I,

S[/h+g& 0

1

12

I

2 ii

1

sm0d0 = PS(AV).

(3

Here, err and tr2 are the principal values of the true stress, P is the internal pressure (transmural) acting on the membrane, AV is the difference between the volume enclosed by the membrane in the current configuration and that enclosed by the membrane in

A finite element model of the infarcted left ventricle

47

r

(a)

0 t noninfarcted

region

-1 t

(b)

I

I

-1

0

1

-1

0

1

I

I

I

I

l-

o-

-_( noninforcted

region

-1 1

0

-1

Fig. 2. Schematic drawing of the ventricular model when the infarct encompasses 25 y0 of the wall volume. (a) diastolic reference configuration. (b) the calculated deformed shape at a transmural pressure of 101 mmHg when the noncontracting region has k, = 16, pp = 2.

the spherical

reference

SV, = (1 +e,)6e,

configuration ++a+

and

6~7s = (1 +e#e,

convenient (6)

where hei, de, and S$ are related to the displacement variations 6u and 6w by equation (1). The expression for s(A V) appearing in the principle of virtual work can be written in the computationally

s(AI’)

form (see e.g. Needleman,

= 2nR;

In diastole,

s

[Sw(l

+e,

+e,+ete,)

-6u(l

+e#]

1977).

0

the myocardium

sined@.

is represented

(7) as an

48

ALAN NEEDLEMAN,STUARTA. RAB~NOWIR,DANIEL K.

REFERENCE SPHERE

B~GENand THOMAS A. MCMAHON

-

DISPLACEMENT VECTOR TO THE CURRENT CONFIGURATION

\

+

/X

I \

‘\\ I I

‘i \ i

Fig. 3. Coordinate system and displacement vector from the spherical referenceconfiguration to the current configuration.

isotropic, incompressible non-linear elastic material. Traditionally, the strain energy function for a nonlinear elastic material is written in terms of algebraic invariants of some deformation measure. However, as noted by Hill (1969), in many applications it is more convenient to write the strain energy function directly in terms of the principal stretches 1,) A2 and I,. The specific form adopted here,

where i., = 1ii, ilz, has been proposed in varying contexts and in varying degrees of generality by Hill (1969). Ogden (1972) and Blatz er a/. (1969). Throughout this paper we will confine our attention to a one-term version of equation (8) which appears adequate for characterizing the passive behavior of both normal and infarcted myocardium (Rabinowitz, 1978; Janz and Waldron, 1978). In formu!ating the systolic model. attention is confined to end-systole, when the contraction velocity is zero. Thus. a great many details of muscle mechanics that would affect. for example, the time course of pressure development in the ventricle. need not be incorporated into the mode) for the present purposes. In end-systole the myocardium is regarded as having undergone a change in reference configuration. Within the framework of the present model, the deformation from the end-systolic reference configuration to the diastolic reference configuration is taken to be one of equal biaxial stretching. Accordingly, for a normal ventricle, the end-systolic reference configuration corresponds to a spherical membrane of smaller radius than the diastolic reference sphere. We denote by

&(i = 1,2, 3) the principal stretches from the endsystolic reference configuration to the current configuration. These are related to the stretches from the diastolic reference configuration to the current one, Ai(i = 1,2,3), by 2, = ici2

x 1 = A&

it3 = 1;213

(9)

where i, is the stretch necessary to deform the endsystolic reference configuration into the diastolic reference configuration. The tension developed at end-systole is taken to be derivable from a strain energy function, b,, which depends on the xi. The specific form adopted here is a one-term version of equation (8), so that

(10) or, alternatively,

by employing equation (9),

The total strain energy is denoted by 4, where in diastole at end-systole.

(12)

In either case, the principal stresses are given by

ai=I.%+a o 1ali

(i = 1,2,3; no sum)

where u,, is the mean normal framework of membrane theory 0

-43

0-

stress. Within

’ a2,

(13) the

(14)

A finite element model of the infarcted

so that the membrane stresses u, and u2 can be written as ui=Ai~-~~$

3

I

(i= 1,2; no sum).

P =

2(;)+-);,

Here, t/R is the current thickness to radius ratio of the sphere and in this spherically symmetrical state u1 = u2 = u. The incompressibility condition (4) leads directly to the second equality in equation (16). We now define the normalized volume, u, to be the ratio of the current volume to the diastolic reference volume. For a spherical membrane, v is simply given by 13. Carrying out the differentiation in equation (15) and substituting u for A3 gives, for a one-term version of equation (8), the diastolic pressure-volume relation

The pressure developed at end-systole, from equations (1 I), (15) and (16) is

Thus, the end-systolic given by p

(15)

Before considering in detail an infarcted ventricle, which requires a numerical solution of the governing equations, we consider the predictions of the present model for a normal ventricle. The model ventricle responds to the increasing pressure in diastole by expanding in a spherically symmetrical manner. Hence, from the well known law of Laplace

obtained

49

left ventricle

_

2;;o

+jz

pressure-volume

relation

[dr,v(k,-31/3_~;2k,t.-(2k,-3)‘3

+,(k,-

3)/3 _

c - (2k,-

3)13)].

1.2

I .4

1.6 1.8 NORMALIZED

(19)

Here we have introduced the factor5 defined by f=

{

1

Cal

0

r
since we take the diastolic stress (and therefore the diastolic pressure) to be zero for normalized volumes smaller than the diastolic reference volume (point 0 in Fig. 1). Figure 1 shows the diastolic and end-systolic pressure-volume curves obtained from equations (17) and (19) when the data described in Bogen et nl. (1980) for the canine left ventricle is used for the elastic constants, k, = 16and pP = 2. In this case to/R0 = 0.20. for reasons to be discussed later. As indicated in Bogen et al. (1980), the value of 2ps to/R0 = 246;and, based on the observations of Spotnitz et 01.(1966) and Grimm et al. (1970), >.
I,’ CHRONI

1.0

is

I 2.0 VOLUME

I 2.2 V/V,

I 2.4

I 2.6

2.8

Fig. 4. Comparison of thediastolic pressure-volume curves computed from the membrane model (solid line) and by the Janz and Waldron (1978) model (broken line). The ventricle with a 417” stiff infarct has an infarct comprising 41.2 on of the wall volume in the Janz-Waldron model and 41.3 “/, in the membrane model. The insert shows Janz and Waldron’s geometry.

50

ALAN NEEDLEMAN,STUARTA. RABINOWW

spherical membrane model (17) is compared with the numerical results of Janz and Waldron (1978), who employed an axisymmetric finite element model of realistic geometry and variable wall thickness. The results for the spherical membrane model of the normal ventricle shown in Fig. 4 are given by equation (17) using the same values of k, and p’pemployed by Janz and Waldron (1978; k, = 14.69 and pp = 0.64), taking the material to be fully incompressible and taking to/R0 to be 0.20 in order to achieve close coincidence with the Janz and Waldron (1978) result. Note, from equation (17), that the ratio to/R, appears linearly, so that the fact that the shapes of the pressure-volume curves are similar cannot be attributed to the choice of this parameter. Part of the coincidence between the predictions of the two models of the normal ventricle can be attributed to the wellknown fact that the actual mammalian ventricle geometry (and that used by Janz and Waldron, 1978) is thicker where the curvature is smaller, thus making the stretch relatively uniform over the whole ventricular surface. Of course, for a sphere the stretch is precisely uniform on the surface. In fact, the pressure-volume curve obtained from a membrane model, consisting of a material having a strain energy function described by a one-term version of equation (8), can be shown to be of the same general shape as that obtained from a thick-walled sphere employing the same material model. Consider a thick-walled sphere with A and B denoting the initial internal and external radii, respectively, and with a and b denoting the corresponding radii in the deformed configuration. Then,

give the normalized volumes enclosed by the inner and outer surfaces of the spherical shell. The pressure versus normalized internal volume curve is given by the relation (Haughton and Ogden, 1978):

where, for a one-term version of equation (8), f&(v)

=

p

[2”kP’3

-

*-2k,‘3

_

3-J.

(23)

KP

Substituting

equation (23) into equation (22) gives

- (0,

-(2k,+31/3_I:

Incompressibility

b-(tk,+31/31].

(24)

implies (25)

Now, consider the limit v, B 1, for the large values of the second term in par-

k, typical of myocardium,

DANIEL K.

B~CENand THOMAS A. MCMAHON

entheses in equation (24) will be negligible compared to the first. Then, equation (25) can be employed to express ub in terms of u, and the expression for dP/dv, can be written as a power series in terms of u,- ‘. The leading term of this power series is dP - = 2p,u, ‘k.--6)J3[;_;(;)-k’+3]. dua

(26)

The analogous formula obtained directly from equation (17) is dP _ = 2clpu(k+‘13 [(+do

l)z].

(27)

Hence, identifying u with the enclosed normalized volume, the slopes of the pressure obtained from membrane theory and the full axisymmetric theory are proportional. On the other hand, for normalized volumes (v,) near unity, a straightforward expansion of equation (24) in powers of u, gives g

= 2p$.,[;-;(;)-3]

(28)

while the membrane theory result, equation (17), yields dP = 2p,k,

du

$

[

1 .

0

(29)

Both the membrane theory and the exact result for a sphere give a linear pressure-normalized volume curve at small values of the normalized volume and a power law pressure-normalized volume curve (with the same power) for large normalized volumes. This analysis suggests that the pressure-normalized volume curve for a relatively thick walled sphere can be fitted, approximately, by the membrane theory formula (17) over a finite range by choosing an appropriate initial thickness to radius ratio, to/R,, which depends on the value of k, as well as on the geometry of the sphere. This serves to rationalize the close coincidence of the diastolic pressure-volume curve obtained from the spherical membrane model with the one obtained by Janz and Waldron (1978). However, additionally, this coincidence also holds for an infarcted ventricle as illustrated in Fig. 4. The infarct comprises 41 y0 of the wall volume (41.2 y0 in the Janz-Waldron model and 41.3 % in the membrane model) and is characterized by k, = 136.0 and ,u~ = 118.0. In the membrane theory analysis, which was carried out employing the finite element method described subsequently, the initial thickness to radius ratio, to/R,, was taken to be uniform and equal to 0.20; the same value employed in fitting the pressure-normalized volume curve for the normal ventricle. Hence, no further adjustable parameter was employed to obtain the close fit for the infarcted ventricle exhibited in Fig. 4. This shows that the change in ventricular compliance due to an infarct can be well represented by the initially spherical membrane model.

51

A finite element model of the infarcted left ventricle FINITE ELEMENT METHOD

Here, certain features of the numerical methods employed to implement a solution to the spherical membrane ventricle model are described. For the diastolic analysis, the starting point is the principle of virtual work (5). Given the pressure, the deformations and the stresses corresponding to some state, not necessarily an equilibrium state, and expanding equation (5) about this state, gives (to the lowest order)

- 2nR; P = 27rR;P

The element stiffness matrices are constructed by integrating equation (30) over the length of each element. Standard finite element procedures, further details of which are given in Rabinowitz (1978). lead to a discretized version of equation (30) of the form A,%& = Pbhx - rhx.

son[re,+e,+

erP2+P,ez)6”-{~(1+e,)+J/P2~Su]sinBdH

O*[(I+e~+eI+

e,e,)6w-$(l

(33)

Here, x is the vector of nodal displacements, Ad is the diastolic finite element stiffness matrix, which is calculated from the left hand side of equation (30). The

(30)

+e,)6u]sinCldH

s

Here. and subsequently, (‘) is termed an increment and denotes the derivative of a quantity with respect to some conveniently chosen monotonically increasing parameter that characterizes the deformation history. An approximate solution is obtained to the linearized variational equation of equilibrium (30) by means of the finite element method. The axisymmetric spherical membrane mode1 of the left ventricle is subdivided into one-dimensional finite elements. Within each element the displacements are taken to be linear functions of the polar angle 0. The size of the elements is allowed to vary with position, so that a fine grid can be used near the boundary between the normal and infarcted muscle without excessively increasing the number of elements. In the calculations to be reported on here, 38 elements were employed. A linear displacement model is assumed in each element so that the displacements, u and w, are related to the nodal displacements u,, wi and u2, w2, by

vector b is obtained by discretizing the first integral on the right hand side of equation (30) and b6x represents the change in volume enclosed by the membrane due to a small displacement variation 6x. Finally, r represents the difference between the internal and external virtual work of the given state, so that if the given state is an’ equilibrium state, c vanishes. Since (33) holds for any displacement variation 6x, it implies the set of algebraic equations A,% = ljb-r.

(34)

In addition to equation (34) the displacement boundary conditions must be satisfied. Two of the displacement boundary conditions arise from the axisymmetric constraint, which requires that the tangential displacement. u, vanish at the poles 0 = 0, II. A third displacement boundary condition is needed to prevent a rigid body translation parallel to the axis of symmetry. Here, this was imposed by setting the tangential displacement, u, equal to zero at one additional point on the membrane. (31) Now, let the given state correspond to an equilibw = wt 1, (0) -t wtfZ (e) rium state at some pressure, P. For a prescribed where fi (8) = (0 - e,)/(e, - 0,) and fi (0) = (0 - e,)i increment of pressure, P, the system of equations (34) (e2 - e, 1. can be solved for the vector of nodal displacement The strain components are related to the nodal increments x. From t, an updated vector of nodal displacements via equation (1) as displacements x is obtained in a linear incremental fashion. The deformations and stresses in each element in the configuration corresponding to the updated nodal displacements can be calculated and, hence, the vector r. Then, equation (34) can be solved repeatedly, ez= 0 with P zero, until an equilibrium state is obtained. + (coWf,u, + (cot@fz~21 (32) This procedure, often termed the incremental Newton-Raphson method, provides an efficient procedure for determining non-linear deformation hisIL=~Lw~f;+w211-u~f,-Ulf~l 0 tories. However, in highly non-linear problems such as wheref; = l/(0, -0,) andf; = l/&-e,). those considered here, this procedure may diverge.

$[w,fl+wzfi

52

ALAN NEEDLEMAN, STUARTA. RABINOWITZ,DANIEL K. B~GEN and THOMAS A. MCMAHON

Two devices were employed in this study to assist convergence. One involved linearly extrapolating from the known states at (P-P) and P to obtain an estimate of the equilibrium state at (P+ P). This estimate was then employed as the initial state in the iterative procedure. The other device was based on the observation that, when the Newton-Raphson procedure diverged, it did so in an oscillatory manner. The oscillations were damped in each iteration by multiplying the solution, x, of equation (34) by some factor less than unity. In the present problems, a factor of one half or one quarter was generally sufficient. The systolic analysis was initiated from the zero stress spherical configuration or from any stage during diastole and was also carried out by means of a finite element based incremental Newton-Raphson method. During systole the rest extension of the non-infarcted myocardium, i,, defined by equation (9) was increased incrementally to its end-systolic value of 1.18, while the change in volume enclosed by the membrane, A V, was required to vanish. Thus, the principle of virtual work (5) takes the form (~i~l+ez~

ZnR;t,

1

1

sinede = 0. (35)

2

equation derived from equation (38) must be satisfied, namely bk = - (A Us.

w

The constraint (40) is imposed directly by employing equation (40) to eliminate one nodal variable in terms of all the others. Denoting this nodal variable by x:rand the vector of the remaining nodal variables by f,, (40) is rewritten as .xe = y+@x, where y = - (A V),/b,. Now, rewrite equation notation as

(39) in partitioned

(41) matrix

Then, employing equation (41) in equation yields the system of equations A:%, = -r;

(42) (43)

where A: = (A,),/ + (A&/J + 8r(&)c/ + (&),,BrB r: = (0, + Br(r,), + Y(A,),, + Y(&8r

(44) (45)

In

equation (35), in addition to the displacement boundary conditions, all displacement fields are subject to the constraint (AV), = 0

(36)

where (A V), denotes the change in volume enclosed by the membrane from that at the onset of systole. We consider some given state, satisfying equation (30), and expand equations (35) and (36) about this state. To the lowest order we obtain

The system of equations (43) can then be solved. Note, however, that the induced pressure increment, P remains to be determined. The method for determining the pressure increment rests on the observation that the governing equations admit an alternative formulation. The constraint (36) can be incorporated directly into the principle of virtual work by introducing a Lagrange multiplier, namely the pressure acting on the membrane. Then, the incremental form of the principle of virtual work is identical to equation (30), so that in discretized form, we obtain f, i 6x = PbGx - F,6x

+s

(e,&, + I&) 1

and (ApI =

J=

1

(46)

where fi,, b and i, are calculated from the integrals appearing in equation (30) as discussed in reference to equation (33). Now, by setting equal to zero the three elements of b that correspond to the three elements of x required to vanish by the displacement boundary conditions, b is an admissible displacement field. Substituting b for 6x in equation (46) and solving for P gives

sined

[(l +e, +e,+e,e,)3

0

-$(l

= - (A V),.

+e&i]

sined (38)

As in the diastolic analysis, the finite element method was employed to obtain a discretized version of the incremental principle of virtual work (40). In systole, the discretized equation is of the form A,tbw = -rr,6x.

(39)

Furthermore, in addition to the displacement boundary conditions which are identical to those employed in the diastolic analysis, the constraint

(47) The solution of (43) and the calculation of P proceed iteratively until convergence is obtained.

NUMERICAL RESULTS

The finite element program was tested for accuracy by comparing the numerical solution for a normal ventricle with the analytical solutions, equations (17) and (19) for the diastolic and end-systolic press-

A finite element model of the infarcted ure-volume relations. The numerical results agreed with the analytical formulas to within 0.5 “/,. An additional check on the diastolic phase of the finite element analysis is provided by the agreement between the results obtained using the present finite element program and the results of Janz and Waldron (1978), as depicted in Fig. 4. The results of the systolic finite element method were compared with those obtained by Bogen (1977) and Bogen and McMahon (1979) who employed a shooting point method to solve an equivalent, but differently formulated, set of membrane equations. An infarct comprising 257; of the wall volume, with the normal and infarcted wall material characterized by the same material parameters used by Bogen (1977), was analyzed. The end-systolic pressure obtained by means of the finite element analysis differed by less than 5 “/afrom that obtained from the shooting point method. In addition to these comparisons an internal consistency check can be made. As can be seen in Figs 7(b) and 8(b), at the border between the infarcted and normal wall material the stresses are discontinuous. However, equilibrium requires continuity of the traction vector at this border and this continuity was maintained in the finite element analysis in all cases checked. Figure 5 displays the diastolic and end-systolic pressure-volume curves for ventricles with infarcts occupying 41.3 7, of the wall volume, but covering a range of possible infarct stiffnesses. In each case,

53

left ventricle

k, = 16.0, while pP = 1.0,2.0, or 16.0. Since for normal diastolic myocardium, pP = 2.0 and k, = 16.0. these three cases represent infarcts with either half, the same. or eight times the stiffness of normal diastolic myocardium. The diastolic pressure-volume curves in Fig. 5 were calculated by employing numerous small pressure increments, typically of 0.2 mmHg; the endsystolic pressure-volume curves were calculated by employing larger volume increments, typically of 0.1 normalized volume units. Figure 5 indicates that the end-systolic pressure-volume curves are displaced further to the right than normal, with decreasing slope, as the infarct extensibility is increased. In addition, the diastolic pressure-volume curves are displaced upwards or downwards from normal with increasing or decreasing infarct stiffness. To answer the question as to whether extensible infarcts interfere with heart function more or less than stiff infarcts, it is necessary to observe that while systolic function suffers with increasing infarct extensibility, diastolic function suffers with decreasing infarct extensibility. The stiff infarct produces a noncompliant ventricle which cannot be filled to a normal end-diastolic volume without exceeding the range of physiological end-diastolic pressures. Thus to evaluate the balance between diastolic and systolic effects, it is useful to calculate the ventricular function curves for the infarcted ventricles, as shown in Fig. 6. The ventricular function curve plots enddiastolic pressure, the input variable to the heart,

0

0.6

0.8

1.0

1.2

14 Normallred

1.6

20

2.2

2.4

valume

Fig. 5. Calculated diastolic and end-systolic pressure-volume curves for ventricles with infarcts occupying 41% of the wall volume. The infarcted myocardium is characterized by k, = 16 and pP = 1.0,2.0 or 16.0;the residual viable myocardium is characterized by the same parameters employed in Fig. 1.

54

ALAN NEEDLEMAN. STUARTA. RABINOWIT~, DANIEL K. B~GEN and THOMASA. MCMAHON

End-diastolic pressure (mmlig)

Fig. 6. Ventricular function curves for 25 “/and 41 “/ infarcted ventricles. Again three infarct stiffnessesare examined; pP = 1.0 ( ),2.0 (---), 16.0(---); k, = 16.The relative disadvantage ofextensible over stiff infarcts at low end-diastolic pressures is reversed at high pressures.

against stroke volume, the output variable, which is simply the difference between end-diastolic and endsystolic volume. Fig. 6 indicates ventricular function curves for the three 419, infarcts described in Fig. 5, as well as three smaller, 25 “/, infarcts covering the same stiffness range. From this figure it appears that for the particular range of infarct stiffness studied here, infarct size is a much stronger determinant of cardiac performance than infarct stiffness. However, a more thorough investigation of infarct stiffness is described in Bogen rr al. (I 980). Nevertheless Fig. 6 indicates that the relative advantages of extensible and inextensible infarcts depend on end-diastolic pressure; if the heart with the extensible infarct is permitted to operate at a high enough filling pressure, it can partially compensate for the impairment in systolic function. Unfortunately, however, this compensation is limited by a physiological limit on end-diastolic pressure of approximately 24 mmHg, beyond which impairment in pulmonary function will develop. The degradation in ventricular performance due to a 41 ‘I,, infarct with k, = 16 and pP = 16.0, (the stiffest infarct considered here) can be attributed to two effects other than the loss of contractile units per se. In diastole the infarct may prevent the surrounding noninfarcted myocardium from stretching to its normal end-diastolic length; and in systole the bulging infarct interferes with the shortening of the adjacent noninfarcted myocardium. These effects are demonstrated in Figs 7 and 8.

The diastolic stretches at 12.0 mm Hg as a function of polar position on the heart are given in Fig. 7(a). Diastolic stretches in a normal heart at this pressure are 1.23. The circumferential and longitudinal stretches are significantly lower in the infarcted region compared to the normal ventricular stretches. In addition. the normal myocardium at the border zone shows a significant decrease in circumferential stretch (iZ = 1.12) from that seen in the residual myocardium far removed from the infarct (AZ= 1.23). Normal circumferential stretches are not achieved until approximately 0 = 120’. Therefore 40”/, of the normal myocardium experiences subnormal stretches as a result of the geometric constraint imposed by the infarct. The diastolic stress distribution in the infarcted ventricle at 12.0mm Hg is shown in Fig. 7(b). The diastolic stress in a normal heart at 12.0 mm Hg is 55.6 mm Hg. The longitudinal stresses show a slight decrease in the infarcted region, while higher than normal circumferential stresses are seen in this area. End-systolic stretches at 101 mm Hg are given in Fig. 8(a). A normal ventricle at this pressure would have a stretch of 0.89. The graph shows that under the force of the normally contracting muscle the infarcted region bulges out, and large extensions of 1.23 are seen in the infarcted region. The bulging infarct region also prevents the adjacent normal myocardium from contracting and the stretches in the border zone of the normal muscle reach as high as 1.05. Thus the infarct

A finite element

I------

(a)

model of the infarcted

55

left ventricle

RESIDUAL

INFARCT

MYOCARDIUM

1.25 1

. ..’ ..*. ....

1.20

.

.:.

. ...’

. . . . . . . ...***..

-1

,:’ /’ ..: /\

.-‘D l.lSr Lo 5 “x w

Q,. circumferential

*:’ ;

j

I

.:’

1 .I0 -

. ..*

_:. . ...

*i ..a*

i

I A,, longitudinal

I.00,

I

I

I

20

40

60 Position

I

I

I

I

100

120

140

160

membrane

(degrees)

80 III

undeformed

(b) l-------INFARCT

RESIDUAL

p,,

MVXARDIUM--

lonqitudincl

i . . ..

1.30

/

. . . . . ..I

..,.... ....*

clrcumferentiol



0’ 0

20

40 Angular

(

60 position

:’

in undeformed

membrane

[degrees)

Fig. 7. Distributions of diastolic stretches (a] and stresses (b) in a 41 ‘I,, infarcted ventricle with k, = 16, pP = 16at 12 mmHg. At this end-diastolic pressure in a normal ventride. diastolic stretch is 1.23and diastolic stress IS 55.6 mm&.

acts as a geometric constraint affecting normal myocardial movement in both diastole and systole. The end-systolic stress distribution in the infarcted ventricle at 101 mm Hg is shown in Fig. 8(b). End-

systolic stress in a normal heart at this pressure IS 196mm Hg. Both the circumferential and longitudinal stresses in the normal myocardium in the border zone are increased. and in fact. higher than normal stresses

56

ALANNEEDLEMAN, STUART

A.

RABINOWITZ, DANIELK.

BOGEN and THOMAS A.

MCMAHON

(a 1

**a....*.*.

.‘...... *...

.. **., .. *. -. .%. :.* : :

:

: **.* ‘..*

X2, circumferential

“*.< ....

0.9 -

0.9\

a

I

I

I

20

40

60

80

Position in undeformed

lb) IiNFARCT

1

I

1

I

109

120

110

169

membrane

(degrees)

-RESIDUAL

MYOCARDIUM-

: : : : : :

, circumferential

!.7

‘2. X* f.* ‘*.* *..*

I

I

40

90

90

Angular position in undeformed

I

I

I

,

loo

120

140

180

membrone (degrees)

Fig. 8. Distributions of end-systolic stretches (a) and stresses (b) in a 41 0,” infarcted ventricle with k, = 16, mmHg. At this end-systolic pressure in a normal ventricle, the end-systolic stretch is 0.89 and the end-systolic stress is 196 mmHg.

pP = 16al 101

occur in about 60 “/, of the residual heart muscle. The increase in circumferential stress at the border zone is 2.7 times the normal level, while the longitudinal stresses are increased 1.2 times in this area. The physiological implications of this finding are discussed in Bogen er al. (1980).

CONCLUDING

REMARKS

By including a region of altered mechanical properties within the simplest model of a ventricle in diastole the mechanical effect of an infarct on the diastolic pressure-volume relation of the ventricle can be

57

A finite element model of the infarcted left ventricle

modelled. Despite the simplifications embodied in such a model, there is good agreement between the diastolic pressure-volume relations given by this model and those given by the model of Janz and Waldron (1978), which fully account for variations through the thickness. In addition, the extension of this membrane model to represent the ventricle in end-systole has been described. A specially tailored finite element method has been developed to enable the end-systolic pressure-volume relation of an infarcted ventricle to be obtained. Thus our model enables us to address the question of how an infarct effects the pumping performance

of a ventricle.

The effect

of infarct

size and

curves has been illustrated. Furthermore. our results show that the presence of an infarct can significantly alter the extension and stress pattern in the surrounding viable myocardium. A detailed study of the impairment of cardiac performance due to myocardial infarction, based on results obtained using the finite element method described here. is given in Bogen et al. ( 1980). stiffness

on

ventricular

function

Acl\no~f~dgemunts-This work was supported in part by grant HL21425. Some of the computations reported on here were carried out on the Brown University, Division of Engineering. VAX-l l/780 computer. The acquisition of this computer &as made possible by grants ?rom the U.S. National Science Foundation (Grant ENG78-193781. the General Electric Foundation, and the Digital Equipment Corporation. REFERENCES

Blatz. P. J.. Chu, B. M. and Wayland, H. (1969) On the mechamcal behavior of elastic animal tissue. Trans. Sot. Rheol. 13, 83-102. Bogen. D. K. (1977)The mechanical disadvantage of myocardial Infarction: a model of the infarcted ventricle. Ph.D. thesis. Harvard University. Boeen. D. K. and McMahon. T. A. (1979) Do cardiac aneurysms blow out? Biophys. J. 27, 301-316. Bogen. D. K.. Rabinowitz, S. A.. Needleman. A.. McMahon. T. A. and Abelmann, W. (1980) An analysis of the mechanical disadvantage of myocardial infarction m the canine left ventricle. Circular~on Res.. ‘47, 728-741. Budiansky. B. f 1968) Notes on nonlinear shell theory. J. uppl. Mrc,h. 35, 393401. Grimm. A. F.. Katele. K. V.. Kubota. R.and Whitehorn. W. V. (1970) Relation of sarcomere length and muscle length rn resting myocardium. Am. J. PI1ysiol. 218. 1412~1416. Haughton. D. M. and Ogden. R. W. (1978) On the mcremental equations m non-linear elasticity--II. Bifurcation of pressurized spherical shells. J. Mech. Ph~s. Solid.t 26, Ill- 138. HIII. R. (1969) Some aspects of the incremental behavior of Isotropic elastic solids after finite strain. Problems in Mrc~hunics. .Amr~ cr.wry

Deformation oj Solid Bodies, Volume, pp. 459-466. Leningrad.

Now-_hiloc

Janz. R. F. and Waldron. R. 1. (1978) Predicted erect of chronic apical aneurysms on the passive stiffness of the human left ventricle. Circularion Res. 42, 255--263. Monroe. R. G. and French. G. N. (1961) Left ventricular pressure-volume relationships and myocardial oxygen consumption in the isolated heart. Circulorion Res. 9, 362-374.

Needleman. A. (1977) Inflation of spherlcal rubber balloon:. Inr. J. Solids Sfrucr. 13. 409-442. Ogden. R. W. (1972) Large deformation lsotroplc elastIcIt> on the correlation of theory and experiment for Incornpresstble rubber-like solids. Proc,. R SCV . L~md A326, 56% 584. Rabinowltz. S A. 1197X) Myocardlal mechanics. con\tltutnc propertles and pumping performance of the ~nti~rctcd ventricle. Ph.D. thesis. Harvard Unibersltb. Rabin0witz.S. A.. Radvan). P.. McMahon.T.and Abelmann. W. (19781 Passive elastic propertles ofnormaland Infarcted myocardium. Cwmicctiotr 58 Suppl. 11. I I- 1959 Sagawa. K. (1978) The rentrlcular pressure \olumc dlasranl re~lslted. C~n~trl~rrro,~Rc\ 43. 677 687 Spotnitz. H. M.. Sonnenbhck. E. H and Spiro. D (1Vhh1 Relation of ultrastructure IO function In the Intact hcarr sarcomere structure relative lo pressure \olumc cur\cs ~~1~ intact left ventricles of dogs (‘lr( I~/~IIOIIRo. 18. 4l) hh Suga. H. (1971) Theoretical anal~ss of’ a iel’r lcntrtcular pumping model based on the systohc time-\ar)lng pres>ure’volume ratio. /EEE Trtrm. B~mwd 18: 4: 55 Suga, H.. Sagawa. K and Shoukas. A. A. cl9731 Load Independence of the Instantaneous pressure \olumr ratlo of the canme left ventricle and effects of epmephrme and heart rate on the ratio. C~rc~rfr~f~or~Rr,r. 32. 314 32’ Suga. H and Sagawa, K. (1974) Instantaneoub prerhure-volume relationshlps and their ratio In the ewsed. supported

canmr

left

ventricle.

C’I~(.I~/~/IIO~

Rt,\

35.

117-126. Suga. H. and Yamakoshl. K. (1977) En‘ects of strokebolumc and velocity of ejection on end-systohc pressure ol canine left ventricle. C~rc~ulorro~~Rr.x. 40. 445 45X

NOMENCLATURE undeformed internal radius of a thick walled sphere deformed internal radius of a thick u’alled sphere assembled diastolic stiffness matrix assembled systolic stiffness matrix undeformed external radius of a thick walled sphere deformed external radius of a thick walled sphere finite element discretization of the Integral giving the change in enclosed volume of a spherical membrane principal values of the linear strain tensor finite element shape functions power law exponent of passive myocardium

power law exponent of active myocardium pressure, transmural (mmHg) radius of spherical membrane in the diastolic reference configuration Newton-Raphson residual vector thickness of spherical membrane In the diastolic reference configuration displacements of points on the spherical membrane (see Fig. 3) nodal values of displacements vector of nodal displacements enclosed volume

normalized enclosed volume (current enclosed volume/reference enclosed volume) normalized volume enclosed by inner radius of a thick-walled sphere normalized volume enclosed by outer radius of a thick-walled sphere vector employed to impose the isovolumet-

ALAN NEEDLEMAN. STUART A. RABINOWIT~ DANIEL K.

58

a(

)

AV

4 li(i = 1, 2, 3) li(i = 1, 2, 3)

ric constraint in the systolic finite element analysis variation change in volume measure of stretch from systolic reference configuration to diastolic reference configuration principal stretches from diastolic reference configuration to current configuration principal stretches from systolic reference

BOGENand THOMASA.

fi’r Pr ;(i = 1, 2, 3)

; Oi(i = 1, 2, 3) *cl

MCMAHON

configuration to current configuration passive elastic modulus (mmHg) active elastic modulus (mmHg) principal Green strains polar angle strain energy function measure of rotation of normal to spherical membrane principal values of the true stress mean normal stress