Residual Stress Estimates from Multi-cut Opening Angles of the Left Ventricle

Purpose Residual stress tensor has an essential influence on the mechanical behaviour of soft tissues and can be particularly useful in evaluating growth and remodelling of the heart and arteries. It is currently unclear if one single radial cut using the opening angle method can accurately estimate the residual stress. In many previous models, it has been assumed that a single radial cut can release the residual stress in a ring of the artery or left ventricle. However, experiments by Omens et al. (Biomech Model Mechanobiol 1:267–277, 2003) on mouse hearts, have shown that this is not the case. The aim of this paper is to answer this question using a multiple-cut mathematical model. Methods In this work, we have developed models of multiple cuts to estimate the residual stress in the left ventricle and compared with the one-cut model. Both two and four-cut models are considered. Given that the collagen fibres are normally coiled in the absence of loading, we use the isotropic part of the Holzapfel-Ogden strain energy function to model the unloaded myocardium. Results The estimated residual hoop stress from our multiple-cut model is around 8 to 9 times greater than that of a single-cut model. Although in principle infinite cuts are required to release the residual stress, we find four cuts seem to be sufficient as the model agrees well with experimental measurements of the myocardial thickness. Indeed, even the two-cut model already gives a reasonable estimate of the maximum residual hoop stress. We show that the results are not significantly different using homogeneous or heterogeneous material models. Finally, we explain that the multiple cuts approach also applies to arteries. Conclusion We conclude that both radial and circumferential cuts are required to release the residual stress in the left ventricle; using multiple radial cuts alone is not sufficient. A multiple-cut model gives a marked increase of residual stress in a left ventricle ring compared to that of the commonly used single-cut model.


INTRODUCTION
Living tissues in the heart continuously interact with their bio-environment, reshape and rearrange their constituents under chemical, mechanical or genetic stimuli during their life cycles. In the mature period, the tissues of a healthy heart remain in a homeostatic state. However, heart diseases disrupt this balance, and cause the tissues to grow and remodel. Physiologically, exercise may also induce healthy and reversible growth and remodelling. An important ingredient in evaluating the mechanics involved in the cardiovascular system is knowledge of the solid mechanical properties of the soft tissues involved, including the components of the heart, such as the left ventricle, henceforth abbreviated as LV. A particular aspect is that the tissues of the heart are residually stressed, so that when the external loading is removed, residual stresses remain in the material. However, residual stresses, which are generally assumed to result from growth and remodeling, are imprecisely characterized (experimentally) at present, and how best to include the important effects of residual stresses in cardiovascular applications therefore presents a modelling challenge.
Over the last century, 13 various hypotheses on the growth and remodelling response to mechanical loading have been put forward, with particular success in arteries. In conventional elasticity theory the existence of a stress-free reference configuration that coincides with the unloaded configuration is normally assumed. 19 In the present context, however, the unloaded configuration is not stress free, but is residually stressed. The residual stress can be estimated using the so-called opening angle method, 2,14 in which an opening angle indicative of the extent of the residual stress can be measured after a single radial cut of an unloaded arterial ring. Using the opened configuration as the reference configuration, the residual stress in a cylindrical artery model can then be estimated. 2,26 This methodology has been extended to multiple cuts by Taber and Humphrey, 27 and used in a two-layered arterial model by Holzapfel et al. 10 Inclusion of residual stress is important in modelling the mechanics of soft tissues for a number of reasons: (a) in nonlinear elasticity theory the stress state in the reference configuration can have a substantial effect on the subsequent response to loads, and omission of the residual stress can lead to significantly different stress predictions under load 9,16,24 ; (b) in biological tissues, the growth and remodelling process will significantly affect the (residual) stress statement in living tissue 1,15 ; (c) while the detailed process of local growth is difficult to measure, residual stress, on the other hand, may be estimated from experiments, as demonstrated by the opening angle measurement. Thus, estimates of the residual stress at particular time instants could provide useful information about the growth history of living tissues.
However, most work that includes residual stress in complex organs, such as the arteries and heart, assumed that a simple radial cut can release all the residual stresses. 23,22,7 However, this assumption is not supported by all experiments. For example, Omens et al. 21 showed that residual stress in a primary mouse heart could be further released by a circumferential cut following the initial radial cut, as illustrated in Fig. 1, They also showed that the opening angles are locationdependent, with greater values at heart apex. This implies that the single-cut opening angle configurations does not correspond to the stress-free configuration, and, as is well known, can only be considered as approximately stress-free. Holzapfel and Ogden 12 used a three-layer (adventitia, media and intimal) model to study residual stress in the artery, where they found each layer has a different opening angle when cutting open separately. In other words, it is impossible to release all the residual stress from a single radial cut across the three layers. By developing a model to count for the three-layer structure of the artery, the estimated residual stress is much greater than treating the artery as a single-layer model.
Myocardium does not have such a distinct layer structure; however, it is clear from the experiment by Omens et al. 21 that multiple cuts need to be considered when studying residual stress in the heart. Inspired by the finding in Ref. 21 and the three-layer modelling work by Holzapfel and Ogden 12 on arteries, in this paper, we estimate the residual stress distribution across the wall of an intact mature heart using multiple cuts in a simplified heart model. Unlike the approach in Ref. 12 where the three distinct artery layers are modelled with different material properties, we know that any material property change across the myocardium must be smooth and that the transmural stress distributions are continuous. These considerations are taken into account in the current work.

METHODOLOGY
There is experimental evidence that the collagen fibres are coiled and wavy in their unloaded state in arteries 3,5 and heart. 22 Similarly to previous studies, 12,28 we assume that collagen fibres are not in tension in the unloaded configuration of the myocardium. Therefore, for modelling the residual stress, we use the isotropic part of the invariant-based constitutive law for the myocardium developed by Holzapfel and Ogden 11 : FIGURE 1. A typical short-axis apical segment of a mouse heart before and after cuts. 21 The initial intact segment, shown in a, was about 2 mm thick. The same segment after a single radial cut and a further circumferential cut is shown in b and c, respectively. In particular, the endocardial segment has reversed its curvature, in c. Note that the definition of the opening angle in Ref. 21 follows that in Chuong and Fung, 2 which is different from that used in the present paper. Reproduced from 21 with permission.
where I 1 ¼ tr C ¼ tr ðF T FÞ, and a, b are material constants. The Cauchy stress tensor is then 11 where B ¼ FF T is the left Cauchy-Green deformation tensor.
For simplicity, we model the LV as an incompressible single-layered cylindrical tube. We consider different scenarios based on different numbers of cuts and assume that the stress in the tube in the absence of loads on its curved surfaces can be released by either a single (radial) cut or multiple cuts (a radial cut followed by one or three circumferential cuts). We also assume that all the cut segments retain their cylindrical configurations, each with its own opening angle, i.e. each is a circular cylindrical sector.

One Cut: Radial
We take basis vectors fe r ; e h ; e z g to correspond to the local radial, circumferential and longitudinal directions, respectively, in the intact circular cylindrical ring (B 3 ). For a single (radial) cut, the opening angle approach has been well described 2,10 for arteries, but is summarised briefly here for completeness. Let the geometry of the right-hand panel in Fig. 2 represent a stress-free configuration B 2 , which is assumed to be a circular cylindrical sector described by cylindrical polar coordinates fR; H; Zg as where R ðiÞ , R ðoÞ , and L denote the inner and outer radii, and the tube length, respectively, and a is the opening angle. Let fE R ; E H ; E Z g be the associated cylindrical polar basis vectors in B 2 . The (isochoric) deformation from B 2 to the intact configuration B 3 is then expressed as where, by incompressibility, r ðiÞ being the inner radius in the configuration B 3 , k ð3Þ z ¼ l=L is the constant axial stretch from B 2 to B 3 , l is the cylinder length in B 3 , and k ¼ 2p=ð2p À aÞ is a measure of the opening angle in B 2 . The outer radius in B 3 is The corresponding deformation gradient (from B 2 to B 3 ), denoted F ð3Þ , is given by It follows that the invariant I 1 (with the superscript ð3Þ omitted temporarily for simplicity) is given by The components of the Cauchy stress tensor in B 3 are then where p is the Lagrangian multiplier. In the absence of body forces the stress components r rr and r hh in B 3 satisfy the equilibrium equation div r ¼ 0, which, for the considered deformation, reduces to and the associated zero-traction boundary conditions are r rr ¼ 0 for r ¼ r ðiÞ ; r ðoÞ . On integration and use of the latter boundary conditions Eq. (13) gives Z r ðoÞ r ðiÞ which, on substitution from (10) and (11,12), can be used to obtain r ðiÞ in B 3 (and r ðoÞ from (5)) when the initial radii R ðiÞ and R ðoÞ and k and c are known. Hence, given W, all the Cauchy stress components can be obtained explicitly in B 3 . It should be emphasized that r is not strictly a residual stress since the presence of the components r zz requires appropriate non-zero boundary conditions, whereas true residual stress is associated with zerotraction boundary conditions. However, the axial force N required to maintain the circular cylindrical configuration B 3 is non-zero but small. Following, Ref. 12 we adjust the axial stress so that r zz À N p½ðr ðoÞ Þ 2 Àðr ðiÞ Þ 2 is the approximate measure of the residual axial stress. The adjusted residual axial stress has its mean removed so that the resulting axial load vanishes.

Two Cuts: Radial and Circumferential
For the two-cut model we consider the separation of the sector in B 2 into two separate circular cylindrical sectors (inner and outer) by means of a circumferential cut around the mid-wall in B 2 at radius R ¼ ðR ðiÞ þ R ðoÞ Þ=2. The two new sectors form the configuration B 1 , which is now taken as the stress-free reference configuration. Thus, B 2 is no longer stress free but requires torsional and axial loads to maintain its circular cylindrical shape. It is, however, residually stressed in the sense that there is no traction on its curved surfaces. The transition from B 2 to B 3 is now different from that in the one-cut model. The stress in B 3 is calculated from the constitutive laws based on the reference configuration B 1 with the appropriate deformation gradients from the two sectors in B 1 to B 3 . The transition from B 1 to B 2 to B 3 is depicted in Fig. 3.
We note, in particular, that the inner sector in B 1 has a negative curvature. The geometry in B 1 is described in terms of cylindrical polar coordinates fR j ; H j ; Z j g, with subscripts j = 1 and 2 corresponding to the inner and outer sectors, respectively. Thus, translate to R ðiÞ and R ðoÞ , respectively, the opening angle is a and the axial length L.
For each sector, the isochoric deformation from B 1 to B 2 can be written as FIGURE 3. The two-cut model a cylindrical model of the LV as the intact ring in B 3 , b after a radial cut to B 2 , and c followed by a circumferential cut to B 1 . Notice that the inner sector in B 1 has a negative curvature, as in Ref. 21. The red curve, at the mid-wall radius R ¼ ðR ðiÞ þ R ðoÞ Þ=2 in (b), separates the inner and outer sectors which become the separate inner and outer sectrors in B 1 after the circumferential cut. Appropriate axial and torsional loads are required to maintain the shapes in B 2 and B 3 : for the inner and outer sectors, respectively, where 1 Þ; j ¼ 1; 2, and k 2j is the axial stretch of sector j = 1 and 2 in B 1 respectively. Note that the negative curvature of the inner sector in B 1 depicted in Fig. 3 is captured by the expression for R in (19), and that structural compatibility in B 2 is ensured since the two expressions for R match at R. As indicated in Fig. 3 the deformation gradient from B 1 to B 2 is denoted F ð2Þ , which is shorthand notation for the two separate deformation gradients from the two sectors in B 1 to B 2 . These are denoted F ð2jÞ ; j ¼ 1; 2, and given by Similarly to the one-cut model, the equilibrium equation in B 2 yields Equation (23) 1 can be rearranged as from which it follows, on use of the zero-traction boundary conditions r RR ¼ 0 on R ðiÞ and R ðoÞ , that i.e. the mean value of r HH through the thickness is zero. It is assumed that is there is no bending moment on the faces H ¼ a=2 and Substitution of r HH from (24) into this equation followed by integration by parts and a further application of the zero-traction boundary conditions leads to Z R ðoÞ R ðiÞ Eqs. (23) 2 , (26) and (27) are solved with Eqs. (19) and (20) to obtain the radii R ðiÞ and R ðoÞ , and the angle a of the sector in B 2 . The deformation gradient associated with the transition F 1!3 from B 1 to B 3 has the form F ð3Þ F ð2Þ , where F ð3Þ is given by Eq. (7) from the one-cut approach, and F ð2Þ is either F ð21Þ or F ð22Þ , as given in Eq. (21). Note that F ð3Þ F ð21Þ and F ð3Þ F ð22Þ are the values obtained for the inner and outer sectors, respectively, and these must match at the interface R, i.e. the deformation gradient must be continuous in B 2 , which implies that k 1 =R 2 at the interface. Enforcing of this requirement will ensure that when traction continuity is applied p is continuous and hence that all the stress components are continuous, in particular that r HH is continuous.
Hence, for the two-cut model, we require geometric information, e.g. R in B 1 , from experiments (Fig. 1c). Then the geometric information of B 2 is obtained from the two-cut model. The residual stress of the intact-ring configuration, B 3 , can now be solved using the same equilibrium equations as (23), (26) and (27), except the deformation gradient is now F ð3Þ F ð2Þ .
An expression for the required axial load N for the intact ring is obtained from a formula similar to that in (15), and the corresponding residual axial stress is adjusted to make the resulting axial load vanish.

Four cuts: one radial and three circumferential
The effect of two further circumferential cuts, one in each of the two separated sectors, is now considered in order to assess if there is any significant change in the resulting calculated residual stress compared with that obtained with a single circumferential cut, although this is not a test that has been carried out experimentally. For definiteness we consider taking a circumferential cut along the mid-wall of each of the two sectors of the two-cut model, leading to the four separate sectors depicted in Fig. 4. In this model we assume that the resulting configuration B 0 is stress-free with no further reversal of the curvature, so that the two sectors in B 1 are no longer stress free.
Each of the four sectors in B 0 is described in terms of cylindrical polar coordinates fq; /; fg according to where q ðiÞ n , q ðoÞ n , a ðnÞ 0 , n ¼ 1; 2; 3; 4, and L 0 are the internal radii, the external radii, opening angles, and the lengths of the four sectors in B 0 , and the notations I and II are identified in Fig. 4.
In B 1 , the geometries of the two sectors are described in terms of polar coordinates fR j ; H j ; Z j g, with indices 1 and 2, as in Eqs. (16) and (17). Next, in B 2 , R ðiÞ , R ðoÞ , a, L denote the internal and external radii, the opening angle and the length of the single sector according to (19)- (20).
The deformations from the four sectors in B 0 to the two sectors in B 1 are described by H 2 ¼k 24 ð/ À pÞ þ p; Z 2 ¼ k   The deformation gradients from B 0 in Fig. 4 to B 1 in Fig. 3 are The radial traction r R j R j should be continuous across each interface R j ; j ¼ 1; 2, and hence, since the deformation gradient is required to be continuous, p is continuous and the other stress components are also continuous. For the considered deformation, with k Z j given, continuity of both the radial and circumferential stresses guarantees that both p and the deformation are continuous. Note that continuity of the circumferential stress at the interfaces in B 1 requires that in the inner (outer) sector : Also, similarly to the two-cut model, The corresponding deformation gradient from B 1 to  (23) 2 , (26) and (27). This model involves eight independent Eqs. (38)-(41), and eight unknown geometrical parameters in B 0 : q ðiÞ n ; a ðnÞ 0 , n ¼ 1; 2; 3; 4. The required external axial load is then calculated similar to (15), and made to vannish by adjuing the residual axial stress.
In summary, for the four-cut model, we require geometric information of both B 1 and B 2 from experiments. Once we obtain all the details, e.g, opening angles and radii of all the sectors in B 0 , we estimate the (residual) stress components r rr and r hh in the intact-ring configuration B 3 , with the total deformation gradient F 0!3 ¼ F ð3Þ F ð2Þ F ð1Þ .

Modelling Parameters
The reference configuration is different for the different models, so we need to define the parameters according to each specific model considered. From B 2 in the one-cut model, we assume that the axial stretch has the constant value k ð3Þ z ¼ 1:14, and R ðiÞ ¼ 2:06, R ðoÞ ¼ 3:20 a ¼ 65 are estimated from the experiments. 21 For the two-cut model, in addition to the parameters used in the one-cut model, we use additional measurements in B 1 of the two-cut model in Ref. 21: We also assume that there is no axial stretch in the transformation from B 1 to B 2 , i.e. k ð2jÞ For the four-cut model, in addition to the parameters used in the two-cut model, we need more geometrical information in B 1 , which is again estimated from the measurements 21 as listed in Table 2. We also assume that k Initially, we consider a homogeneous myocardium model, for which the material parameters of the con-stitutive law (1) are fitted to the data of mice. 20 This gives a ¼ 2:21kPa, and b ¼ 1:8.
However, the dramatic difference in the maximum hoop stress between the single-cut and multiple-cut models raises a question about the rationale of considering homogeneous material properties. Novak et al. 18 showed that canine myocardium is heterogeneous with location dependent material properties. In particular, they showed that the mid myocardium is softer than epicardium or endocardium, although they didn't find differences in regional (anterior wall vs septum) stiffness. We therefore also consider an inhomogeneous myocardium model, and fit the data from 18 with the strain-energy function (1).
The parameters are spatially dependent, as shown in Fig. 5. Since our model is for mice, and there are no experimental data on the heterogeneous properties of mice myocardium, we take the spatial variation of the canine data, but keep the mean values of the mice data from our fitted parameters, and include these in our calculations.

Results for the Homogeneous Myocardium Model
The geometry of the intact ring predicted in the three different models is summarized in Table 3, and compared with the measured data in Ref. 21.
It is clear that the agreement of the estimated radius and thickness of the unloaded configuration gets better as the number of cuts increases. Although in principle the true zero-stress configuration requires infinite cuts, Table 3 suggests that two or four cuts provides a good approximation of the zero-stressed configuration, given that the measured geometry is not exactly circular, although it is assumed to be circular in each model.
The components of the residual stress distribution in the intact-ring configuration B 3 from the three different models are shown in Fig. 6. The residual axial stress for each model is adjusted to remove the impact  of the non-zero values of N (=5:67 Â 10 À3 , 1:71 Â 10 À3 , and 1:42 Â 10 À3 , respectively, in the onecut, two-cut, and four-cut models.) Although the overall distributions are similar in all these models, there is marked difference in the magnitudes of the hoop stresses. In particular, the maximum r hh is 1.75, 17.13, and 17.15 kPa, respectively, for the one-cut, two-cut, and four-cut models. Comparing to the onecut model, the ratio of the maximum hoop stresses over the single cut is about 9.78 times for the two-cut model, and 9.80 times for the four-cut model. We also notice that although the four-cut model gives much smoother stress distribution, the two-cut model leads to a similar magnitude of the maximum hoop stress. This suggests that the significant rise in the residual stress is due to the negative curvature at the first circumferential cut.

Results for the Heterogeneous Myocardium Model
The results of the intact ring from the heterogeneous myocardium model are computed. Again, the residual axial stress for each model is adjusted to remove the impact of the non-zero values of N (=4:59 Â 10 À3 ,  6. Distribution of the residual stress components in the intact ring from a single cut, b two-cut, and c four-cut models based on the homogenous material assumption. 1:35 Â 10 À3 , and 1:09 Â 10 À3 , respectively, in the onecut, two-cut, and four-cut models.) . All the stess compoenets are plotted in Fig. 7, which shows that the maximum r hh is 1.85, 16.12, and 16.63 kPa, respectively, for the one-cut, two-cut, and four-cut models. The ratio of the maximum hoop stresses over the single cut is slightly lower than that of the homogeneous material, at about 8.71 times for the two-cut model and 8.73 times for the four-cut model. Again, apart from smoother and somewhat different stress distributions, the four-cut model predicts very similar maximum hoop stress as the two-cut model.

DISCUSSION
The issue of multiple cuts has been studied before. In particular, Fung suggested that one radial cut might be sufficient to release all the residual stress. This was supported by his experiments which showed that after two or more radial cuts, no obvious deformation occurs from the one-radial-cut configuration of the arteries. 6 Using our models we can show, however, that multiple radial cuts do not indeed release more residual stress, since due to the symmetry of the considered geometry, once a radial cut is made, no further elastic deformation can occur after more radial cuts.
In other words, the deformation gradient in each of our different models is independent of the azimuthal angle and the solutions are also independent of this angle. However, this does not indicate that the one-cut configuration is a stress-free one. To further release residual stress circumferential cuts following a radial cut are necessary. We note that even with multiple cuts, we may not release all the residual stresses. In principle, only infinite cuts can release all the residual stresses. This explains why the stress distributions are not smooth and appear to have deflections around where the cuts are, since we have to assume the configuration after two or four cuts is stress-free. This process could be improved with more cuts, though more than four cuts cannot be modelled in this paper without further experimental data. However, the curves in the four-cut model are almost smooth, the required axial load to maintain the cylindrical shape of the intact ring reduces is smaller, and the predicted intact ring thickness agrees very well with the measurements. All these tentatively suggest that the fourcut model is already a good approximation for the residual stresses. Indeed, in terms of estimating the residual stress magnitude, even the two-cut model seems to be good enough.
The fact that our two-cut and four-cut models predict a much higher (about 8-9 times!) hoop stress is perhaps not unexpected given the large negative curvature revealed by the experiments. 21 The increased value of the residual hoop stress agrees with the rough estimation by Omens et al. 21 based on a simplified concentric cylindrical shell model, where they showed that the maximum loop stress was about 20kPa, which is close to our estimation of about 17kPa.
We argue that the significant stress underestimate of the one-cut model does not just occur in the heart models. In the residual stress modelling of arteries, by treating the artery wall as three separate layers (intima, media and adventitia), and measuring the opening angles for each of the three layers, Holzapfel and Ogden 12 have essentially developed a three-cut model (one radial cut followed by two circumferential cuts, in this case separating the layers with different properties). Their model is also heterogeneous, as different material parameters are used for different layers. We now compare the residual stress distribution across the artery wall in Fig. 8 using their three-cut approach and the one-cut model. The results from the one-cut (or one-layer) model are reproduced here using the same parameters as in Ref. 12. The maximum values of the hoop stress and their ratio to the one-cut model results in different layers of the three-cut model are listed in Table 4, which shows that the ratio of the hoop stress in the different layers ranges from 24 to 50 times. This difference is even greater than that of the mouse heart.
In this study, we assume that the deformation gradient is diagonal. Note this is in general not true for the heart under loading. 8 However, since residual stresses are estimated from the unloaded configuration, the fibres are in general coiled and do not bear the load. Therefore the myocardium behaves like an isotropic material, and the radial, circumferential and axial directions are the principal directions of the deformation, i.e. the deformation gradient is diagonal in these directions. This is different from artery modelling, when the two families of fibres are assumed symmetrical about the axial direction. Hence, even when loaded with pressure, tensioned fibres do not alter the principal directions of the deformation. Therefore, the diagonal deformation gradient (e.g. Ref. 12) is assumed in arteries for a different reason.
Finally, we would like to state the limitations of the study. We have assumed that cross-sections of the heart are cylindrical, and all the cut segments retain their cylindrical configurations, each with its own  opening angle, i.e. each is a circular cylindrical sector. Needless to say, this is not always true, particularly in opening angle studies of arteries, where a single cut of an arterial ring can have a non-uniform curvature. 25 In addition, with the cylindrical assumption, we cannot model the experimental observation that opening angles are location-dependent, with higher values at heart apex. 21 Indeed, the ''true'' residual stress can only be achieved through infinite cuts using the opening angle method, which is not practical. We believe the residual stress estimation, although an approximation, is a step forward to the physiological range compared with using the single cut opening angle method.

CONCLUSION
Based on experimental observations that a single radial cut does not release all residual stress, we have used multiple cuts to estimate the residual stress distributions in a mouse left ventricle model. Our results show that both radial and circumferential cuts are required to release the residual stresses in the middle wall of the left ventricle. Remarkably, using radial cuts alone leads to a significant underestimate of the residual stress, which will be around 8 to 9 times greater if estimated on the basis of combined radial and circumferential cuts. Similar findings apply to arteries based on the model in Ref. 12. We remark that the results are not significantly different using the homogeneous or heterogeneous material models. In addition, although the stress distributions are different and much smoother in the four-cut model, the two-cut model can already estimate the maximum hoop residual stress quite satisfactorily.