A model to determine the effect of collagen fiber alignment on heart function post myocardial infarction

Abstract

Background

Adverse remodeling of the left ventricle (LV) following myocardial infarction (MI) leads to heart failure. Recent studies have shown that scar anisotropy is a determinant of cardiac function post-MI, however it remains unclear how changes in extracellular matrix (ECM) organization and structure contribute to changes in LV function. The objective of this study is to develop a model to identify potential mechanisms by which collagen structure and organization affect LV function post-MI.

Methods

A four-region, multi-scale, cylindrical model of the post-MI LV was developed. The mechanical properties of the infarct region are governed by a constitutive equation based on the uncrimping of collagen fibers. The parameters of this constitutive equation include collagen orientation, angular dispersion, fiber stiffness, crimp angle, and density. Parametric variation of these parameters was used to elucidate the relationship between collagen properties and LV function.

Results

The mathematical model of the LV revealed several factors that influenced cardiac function post-MI. LV function was maximized when collagen fibers were aligned longitudinally. Increased collagen density was also found to improve stroke volume for longitudinal alignments while increased fiber stiffness decreased stroke volume for circumferential alignments.

Conclusions

The results suggest that cardiac function post-MI is best preserved through increased circumferential compliance. Further, this study identifies several collagen fiber-level mechanisms that could potentially regulate both infarct level and organ level mechanics. Improved understanding of the multi-scale relationships between the ECM and LV function will be beneficial in the design of new diagnostic and therapeutic technologies.

Background

Coronary heart disease is the leading cause of death, accounting for over 400,000 lives in the United States every year [1]. Blockage of the diseased coronary arteries leads to myocardial infarction (MI). Post-MI, the left ventricle (LV) undergoes a complex remodeling process that results in the formation of a scar or infarct. The increased stiffness and diminished contractility of the scar reduce LV function and can lead to heart failure [2]. The decrease in LV pump function has been correlated to the size of the infarct, but infarct size alone cannot explain patient outcome [3]. One key factor affecting LV function is the mechanical properties of the scar tissue. Increased infarct stiffness reduces inflation under systolic pressures, but also impairs filling of the LV during diastole [4]. Therefore, it is imperative to elucidate the relationship between scar composition and mechanics in order to identify the properties that best preserve LV function.

Cardiac scar tissue is composed of collagen fibers, the primary determinant of the mechanical properties of the scar, as well as cells and other extracellular matrix proteins which are considered to be a ground substance [2]. Thus, cardiac scar tissue demonstrates anisotropic, non-linear material behavior [5]. Naturally, the density and alignment of collagen fibers affect the mechanical properties of scar tissue and thus the wall stress, with scar tissue demonstrating a high stiffness along the collagen fiber direction and a lower stiffness in the cross-fiber directions [6]. Furthermore, it has been shown that collagen deposition during remodeling is highly controlled by the local stresses or strains in the tissue [7–9]. In arteries it has been shown that the angle of mean collagen fiber alignment adapts to optimize a tissue’s load-bearing ability which leads to fibers that are aligned in a direction in between the directions of the two principal stresses [10]. In fact, recent studies by Fomovsky et al. have shown that mechanical cues are the primary determinant of collagen alignment [11]. Therapeutic strategies that aim to limit infarct expansion by reducing local stress are currently under development [12, 13], including efforts to anisotropically reinforce the infarct region [14]. Thus, it is of clinical interest to determine the relationship between collagen alignment, scar tissue mechanics, and LV function, in order to find a possible alignment that would better preserve LV function.

The objective of this study was to illustrate how collagen alignment affects LV function post-MI using a simple, thin-walled, cylindrical LV model with a scar region governed by a collagen fiber based constitutive equation. The cylindrical model is based on a concept, previously presented by Han et al. [15] and Oshinski et al. [16], to predict ejection fraction from MRI images, and the collagen fiber-based constitutive equation was based on the work of Grytz and Meschke [17]. This simple multi-scale model allows us to focus on identifying the properties of the collagen matrix that primarily determine LV function post-MI.

Main text

Methods

Constitutive model

The LV was modeled as a cylindrical membrane formed by folding a planar sheet consisting of three normal, healthy myocardial regions and a collagenous scar region into a cylinder (Figure 1). The mid-wall circumference of the cylinder was taken to be the length of the planar sheet in the circumferential direction. The healthy LV regions were modeled as sheets with nine layers having myocyte fiber angles ranging from −50 to +50 degrees relative to the circumferential direction. While reports of myocyte rotation in the literature are somewhat variable between species, location in the ventricle, and measurement technique [18–20], the distribution parameters chosen provide a reasonable physiological approximation. The healthy region was governed by a Fung-type exponential strain energy density equation [21, 22]:

$$W=\frac{1}{2}c\left(e^Q-1\right),$$

(1)

with

$$Q=b_{\mathit{ff}}{E_{\mathit{ff}}}^2+b_{\mathit{xx}}\left({E_{\mathit{cc}}}^2+{E_{\mathit{rr}}}^2\right)+b_{\mathit{fx}}\left({E_{\mathit{fc}}}^2+{E_{\mathit{cf}}}^2\right).$$

(2)

Figure 1

Schematic of the four-region cylindrical model. The heart wall was modeled as blocks of membrane of normal contraction (R: remote) and scar (S) with roman numeral supercripts used to denote the three different remote regions. Labels are given to denote the circumferential length (L), longitudinal height (H), and wall thickness (T) of the regions.

Where c is a material constant, b _ff , b _xx , and b _fx are non-dimensional material constants and E _ff , E _cc , E_rr,E _fc , and E _cf are the components of the Green strain tensor. The subscript f denotes the fiber direction, the subscript c denotes the in-plane cross fiber direction, the subscript r denotes the radial (out-of-plane) direction and the subscript x denotes either the cross-fiber or radial direction. The Fung-type constitutive equation has been widely used in soft tissue biomechanics and is still commonly used for modeling myocardial tissue [23, 24]. Since the collagen density of the healthy tissue is much lower in the healthy regions of the myocardium compared to the scar region, no explicit model of collagen is included for these regions. Cauchy stresses were determined from the strain energy density functions following [25],

$$\left[\mathit{\sigma }\right]=-p\left[\mathit{I}\right]+2\left[\mathit{F}\right]\cdot \frac{\partial W}{\partial \left[\mathit{C}\right]}\cdot {\left[\mathit{F}\right]}^T,$$

(3)

where p is a Lagrange multiplier used to enforce incompressibility, [I] is the identity matrix, [F] is the deformation gradient, and [C] is the right Cauchy-Green deformation tensor.

Contraction of the healthy tissue was modeled using the elastance model of active cardiac contraction proposed by Guccione et al. in the form of [26, 27],

$${\sigma }_{\mathit{active}}={{\sigma }^0}_{\left(l_s\right)}\cdot \left(1-\mathrm{co}{\mathrm{s\omega }}_{\left(l_s,t\right)}\right)$$

(4)

where l _s and t are the the sarcomere length and time during a cardiac cycle, respectively. ω is a function of the sarcomere length and time. σ⁰ is the peak contraction stress that depends on cellular calcium concentration and sarcomere length [27]. The total stress in the healthy regions is calculated as the sum of the active stress which is non-zero solely along the direction of myocyte orientation, and the passive Cauchy stress tensor. The individual layers of the healthy tissue were assumed to have the same planar deformation (i.e. the transmural shear deformation was ignored).

The scar tissue was modeled as a thin sheet composed of a distribution of collagen fibers embedded within a ground matrix. The mechanical properties of the scar region were governed by the microstructure-based nonlinear model proposed by Grytz and Meschke [17, 28]. Briefly, the stretch of collagen fibers was modeled as the uncrimping of a coiled collagen filament. This is similar to the manner in which a spring uncoils as it is stretched. The strain energy density function for the collagen fiber component of the tissue, W _col , was calculated by integrating the first Piola-Kirchoff stress, P _col which is a function of the stiffness of the collagen filament, E _col , the ratio of the fiber diameter, D _col , to the filament diameter, d _col , and the initial crimp angle, θ _o , (a measure of how tightly coiled the fiber is, with 0° representing a fully straightened filament and 90° representing a filament that has been compacted into a disc, see Figure 2) over the deformation of the fiber, λ.

$$W_{\mathit{col}}\left(\mathrm{\lambda }\right)={\displaystyle {\int }_1^{\mathrm{\lambda }}}{P}_{\mathit{col}}\left(\mathrm{\lambda },E_{\mathit{col},}\frac{D_{\mathit{col}}}{d_{\mathit{col}}},{\theta }_o\right)d\mathrm{\lambda }$$

(5)

where λ depends both on the probability density function of the collagen fiber alignment, ρ, and the right Cauchy-Green deformation tensor [C] [17, 28]. The ground matrix was modeled as an isotropic, neo-Hookean material with a strain energy density

$$W_g=c_g\left(I_1-3\right)$$

(6)

where c _g is a material constant and I₁ is the first invariant of the right Cauchy-Green deformation tensor. The total strain energy density, W, was calculated from the volume weighted summation of the strain energy densities of the collagen fibers, W _col , and the ground substance W _g ,

$$W=\left(1-{\omega }_{\mathit{col}}\right)W_g+{\omega }_{\mathit{col}}{W}_{\mathit{col}},$$

(7)

where ω _col is the collagen volume fraction. For a detailed description of the calculation of W _col please refer to the original work by Grytz and Meshke [17, 28]. Accordingly, the Cauchy stresses were determined from this strain energy density functions using Eq. (3).

Figure 2

Schematic of crimped collagen fiber. The collagen filament coils around a central axis that is aligned in the fiber direction. The diameter of the filament, the diameter of the fiber helix, and the crimp angle are illustrated.

Boundary conditions and equilibrium equations

For clarity, regions will be referred to as R^I, R^II, R^III, representing the healthy, remote regions, and S representing the scar region (Figure 1). In order to solve for the deformed volume of the LV it is necessary to establish boundary conditions and equilibrium equations to relate deformation to the applied loads. The base of the LV is fixed in the longitudinal direction, but not in the radial or circumferential directions. Further, it is assumed that the LV will remain cylindrical throughout deformation. Given this boundary condition and the assumptions of thin walls and incompressibility, the deformation gradient, [F], for each region has only 3 unique components:

$$\left[\mathit{F}\right]=\left[\begin{array}{ccc}\hfill {\lambda }_{\theta }\hfill & \hfill {\kappa }_{\mathit{\theta z}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_z\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{{\lambda }_{\theta }{\lambda }_z}\hfill \end{array}\right].$$

(8)

The components λ _θ and λ _z represent the stretch ratios in the circumferential and longitudinal directions respectively and κ _θz represents the shear deformation. Since uniform stress distributions over each region are assumed, internal forces can be directly related to the Cauchy stress through the definition [25]

$${\mathit{q}}_i={\mathit{n}}_i{a}_i\cdot \left[\mathit{\sigma }\right]\mathit{for}i=\theta ,z,$$

(9)

where n _i is the normal to the i^th surface and a _i is the surface area of the i^th surface. The internal forces must balance the forces generated by the internal pressure. Under the assumptions of the Law of Laplace for thin-walled closed-end, cylindrical pressure vessels the magnitude of the internal forces are

$${||\mathit{q}||}_{\theta }=\mathit{Prh},$$

(10.a)

$${||\mathit{q}||}_z=\mathit{\pi P}{r}^2,$$

(10.b)

where q _θ and q _z are the internal circumferential and longitudinal forces, P is the pressure, r is the deformed radius of the LV, and h is the deformed longitudinal height of the LV.

To model the coupling of the regions and the division of load amongst adjoining regions, the LV wall was modeled as a 2-D system of springs (Figure 3). Based on this arrangement, the internal loads described in Eq. (10) are shared as

$${||\mathit{q}||}_{\theta }=||{{\mathit{q}}_{\theta }}^I||+||{{\mathit{q}}_{\theta }}^{\mathit{III}}||=||{{\mathit{q}}_{\theta }}^{\mathit{II}}||+||{{\mathit{q}}_{\theta }}^S||=\mathit{Prh},$$

(11.a)

$$||{\mathit{q}}_z||=||{{\mathit{q}}_z}^I||+||{{\mathit{q}}_z}^{\mathit{II}}||=||{{\mathit{q}}_z}^{\mathit{III}}||+||{{\mathit{q}}_z}^S||=\mathit{\pi P}{r}^2.$$

(11.b)

Figure 3

Schematic of force division and balance. The LV is modeled as a system of springs with a longitudinal and circumferential spring for each of the four regions. Springs shown in parallel share displacements and springs shown in series share forces. The applied forces are calculated from the Law of Laplace.

The spring system model was also used to describe which regions share displacements. The deformations are described by the components of the deformation gradient tensor, [F], as seen in Eq. (8). In the circumferential direction, regions R^I and R^III are fixed in parallel so that they will always have the same circumferential deformation λ _θ ^R. The same goes for regions R^II and S which share circumferential deformation λ _θ ^S. In the longitudinal direction, regions R^I and R^II are fixed in parallel so that they will share longitudinal deformation λ _z ^R and shear deformation κ _θz ^R. The same is true for regions R^III and S that share λ _z ^S, and κ _θz ^S. From Eq. (11) a full set of equilibrium equations can be derived which in conjunction with the constitutive equations can be solved numerically to determine the unknown displacements. A more complete derivation appears in the appendix.

Numerical solution methods

The LV volume contained by each region was determined independently as:

$$V^i=\left(\frac{l^i\cdot r}{2}-\frac{l^i\cdot t^i}{2}+\frac{l^i\cdot {\left(t^i\right)}^2}{2r}\right)\cdot h^i,$$

(12)

where, r is the radius of the LV, and lⁱ, hⁱ, and tⁱ are the deformed length, height and thickness of each of the four regions R^I, R^II, R^III and S assuming that the deformed LV shape remains cylindrical. The total volume of the ventricle was determined by combining the component volumes. All numerical operations were carried out using Matlab (Mathworks, Natick, MA). An iterative solver based on Newton’s method with derivatives calculated by complex Taylor series expansion was used to simulate diastolic filling from 0 to an end diastolic pressure of 1.3 kPa (9.8 mm Hg) over 100 pressure increments. Convergence testing showed that end-diastolic volumes calculated using 50 pressure increments were only 0.1% less than those calculated using 100 increments. To simulate systole, 10 points on the end-systolic pressure volume relationship (ESPVR) curve were calculated using a Newton’s method solution technique to solve the force balance equations of Eq. (11) with the added time-varying myocyte contraction forces and an added isovolumic constraint so that pressure could be calculated as a function of time. Each systolic contraction was calculated using 100 time points. From this ESPVR the final systolic volume was taken to be the point on the curve at 13 kPa (98 mm Hg). Stroke volume (SV) and ejection fraction (EF) were then calculated from the end diastolic and systolic volumes.

Numerical examples

A set of typical material properties and geometric dimensions was chosen as an example problem for our simulation. LV dimensions were chosen so as to roughly mimic the dimensions of a canine left ventricle [29]. Total LV height was chosen to be half of the mid-line circumference. The thickness of the infarct was chosen to be half the thickness of the healthy region as well as half of the total circumference and half of the total height. In this configuration the scar represents 25% of the total myocardial area representing a moderately sized scar. The final dimensions were chosen so that the LV volume in the undeformed state matched that reported by Fomovsky et al. [24]. The material constants of the Fung equation were taken from Fomovsky et al.[24], with incompressibility assumed. The properties governing myocyte contraction were the same as those reported by Guccione, and an initial sarcomere length of 1.65 μm was chosen as this resulted in a baseline ESPVR that roughly matched the baseline case reported by Fomovsky et al. [24]. The collagen fibers of the scar region were aligned following a Von Mises distribution wrapped between −90° and +90°:

$$\rho \left(\phi \right)=\frac{\exp \left(\mathrm{\gamma }*cos\left(2\mathrm{\phi }\right)\right)}{\mathrm{\pi }{\mathrm{I}}_{\mathrm{o}}\left(\mathrm{\gamma }\right)}$$

(13)

where γ is a fiber alignment factor such that high values of γ correspond to highly aligned fibers, and I _o is the modified Bessel function of order zero. The properties of the scar region were determined by assuming a collagen fiber density of 40% and a circumferential mean fiber alignment and fitting the remaining parameters (crimp angle, young’s modulus, neo-Hookean constant, ratio of helix diameter to filament diameter, and fiber alignment parameter) to the stress strain relationship reported by Gupta et al. [5] for a 6 week old infarct in a porcine model. A circumferential alignment was chosen as this was the axis of the largest stress under equibiaxial stretch in the experimental data. The values for these parameters are listed in Table 1. Least-square fitting was performed using a gradient based method and the R² value was 0.996 for the circumferential stress and 0.992 for the longitudinal stress (Figure 4). Using the properties described in Table 1 as a reference point, parametric studies were run to determine the effects of varying the mean collagen alignment angle, collagen density, collagen fiber alignment, collagen fiber stiffness, collagen crimp angle, scar size, and scar thickness. Additionally, a baseline case with 4 healthy regions was run for comparison.

Table 1 Parameters of the example infarcted LV model

Figure 4

Fitting of material properties to experimental data. The properties of the scar tissue were fit to experimental data obtained by Gupta, et al. 1994 for a 6-week old porcine infarct. (A) The fitting to the circumferential stress resulted in an R² value of 0.996. (B) The fitting to the longitudinal stress resulted in an R² value of 0.994.

Results

The SV and EF of the baseline (pre-MI) case were 22.0 mL and 0.486 respectively. In order to investigate the effect of varying the mean fiber angle, an example infarcted LV was run 19 times with the mean fiber angle ranging from 0° (circumferential alignment) to 90° (longitudinal alignment). Simulation results showed that SV and EF increased with fiber angle and the maximum SV and EF were achieved when the fibers were aligned at 90° (Figure 5A-B). SV and EF were reduced compared to baseline through lowered diastolic filling and a right shift of the ESPVR (Figure 5C). Longitudinally aligned fibers improved diastolic filling when compared to circumferentially aligned fibers. This extra diastolic deformation also allowed myocytes in the healthy regions to generate more contractile force which slightly reduced end-systolic volume as seen in the ESPVR’s in Figure 5C. At angles intermediate to 0° and 90°, SV and EF rise steeply from around 20° to 60° before plateauing.

Figure 5

Cardiac function of an example infarcted LV. Changes in EF (A) and SV (B) with the mean alignment angle of the collagen in scar tissue. Dotted lines are the baseline values corresponding to the healthy LV. Both EF and SV were maximized when collagen fibers were aligned longitudinally. The pressure volume relationships (C) show that SV is reduced through both reduced diastolic filling and a right shift in the ESPVR. Longitudinal alignment of fibers reduces both of these effects.

The deformations of the scar and remote regions (S and R^I), as measured by the components of the deformation gradient (see Eq. (8)), were also examined as the mean collagen angle was varied (Figure 6). When collagen aligns towards the longitudinal direction (90°), the circumferential deformation of the scar region increased while the longitudinal deformation at end-diastole slightly decreased (Figure 6A). Longitudinal fiber alignment also reduced the longitudinal deformation of the R^I region while having only a slight effect on the circumferential deformation at end-diastole (Figure 6B). At end-systole, longitudinal fiber alignment greatly reduced longitudinal deformation in the scar region while only slightly increasing the circumferential deformation (Figure 6C). Longitudinal fiber alignment reduced longitudinal deformation in the remote R^I region but also slightly increased the circumferential deformation at end-systole (Figure 6D).

Figure 6

Deformation and stress. Diastolic deformation of the scar (A) and remote region R^I(B) and systolic deformations of the scar (C) and R^I regions (D) plotted as functions of mean collagen alignment angle. The normal components of the deformation gradient (circumferential and longitudinal) are shown on the left vertical axis and the shear component is shown on the right axis.

Several parametric variations in material properties and geometry were simulated and are presented as contour plots (Figure 7). Increasing the size of the scar (as measured by the percentage of the total LV area greatly reduced SV (Figure 7A). The effect of increased longitudinal alignment was greatest for scars that were about 20-30% of the total LV area. Reducing the thickness of the scar tended to improve SV by allowing for greater diastolic deformation (Figure 7B). However, decreasing scar thickness below 20% of the thickness of the healthy regions resulted in decreased SV. Increasing the collagen fiber density improved SV for longitudinally aligned fibers, due primarily to increased diastolic deformation in the cross-fiber direction (Figure 7C). For circumferential alignments, increased collagen density only slightly reduced stroke volume. Increasing the neo-Hookean constant governing the ground substance of the scar region reduced the SV (Figure 7D). This effect was larger for longitudinally aligned fibers. When the Young’s modulus of the filaments was low, SV was correspondingly low (Figure 7E). Optimal stroke volumes were found when fibers were aligned longitudinally and the Young’s modulus was about 1 MPa. Increasing the Young’s modulus to 10 MPa tended to reduce SV, however increasing the modulus further to about 100 MPa increased the SV. Varying the initial crimp angle revealed that optimal SV occurred at around 20 degrees (Figure 7F). Reducing the alignment factor γ, a change that represents a more disperse distribution of fibers, increased the SV when the mean angle was near 0° but reduced the SV when the angle was near 90° (Figure 7G and H).

Figure 7

Parametric study of SV. Contour plots of stroke volume are presented with collagen alignment angle on the x-axis and one of the following parameters of interest on the y-axis. (A) Infarct size, presented as percent area of total LV area. (B) Wall thickness, presented as the ratio of the thickness of the infarct region to the thickness of the remote region in percent. (C) Collagen density presented as the percent of scar volume. (D) The log of the neo-hookean constant for the ground material of the infarct in log(Pa). (E) The log of the Young’s modulus of the collagen filament in log(Pa). (F) The initial crimp angle of the collagen filament in degree. (G) The log of the fiber alignment factor γ. (H) To aid in the interpretation of the contour plots, a plot of stroke volume vs. mean collagen alignment angle is presented for three levels of the fiber alignment factor γ.

Discussion

We developed a cylindrical LV model that links the wall stress and pumping function of the dysfunctional LV to the density and alignment of collagen fibers in the scar tissue. Using this model, we examined the effect of collagen alignment and density on LV function post-MI. Our results demonstrate that the orientation of collagen fibers in the scar region of the myocardium has a large influence on the performance of the LV. It was found that longitudinal fiber alignment limits scar deformation during systole and maximizes LV performance.

The current model clearly demonstrates how the mechanics of the healthy contractile tissue are coupled to the properties of the collagen network in the scar region. The conclusion of longitudinal aligned fibers maximizing the LV function agrees with reports in the literature [24]. This study also ellucidates the mechanism behind this improvement. The results suggest that rather than increased longitudinal stiffness improving function it is actually increased circumferential compliance that has the largest benefit. Increased circumferential compliance allows for the greatest diastolic filling, which in turn allows for greater contractile ability through the Frank-Starling mechanism. It is worth noting that an infarct that is too compliant may eventually lead to reduced systolic function as it inflates under systolic pressure. Examining the ESPVR’s in Figure 5C, it is seen that the longitudinal and circumferential alignment curves cross once the pressure exceeds about 16 kPa indicating that the increased circumferential compliance is leading to increased infarct bulging under high pressures.

While the parametric studies demonstrated that smaller scars preserve LV function they also revealed that for certain scar sizes, thinner scars increase LV function, a far less intuitive conclusion. This is due to the fact that thinner scars allow for more diastolic filling. There is a limit to this though, as scars with an infarct thickness to remote thickness ratio of less than 20% had a reduced SV due to increased systolic bulging. We do not endorse infarct thinning as an advantageous phenomenon as it increases the risk of cardiac rupture and infarct expansion.

The finding that increased collagen density has a positive effect on SV for longitudinally aligned fibers but very little effect for circumferential alignments is explained by the fact that as more collagen fibers replace the ground substance, the tissue becomes more compliant in the cross-fiber direction but stiffer in the direction of the fibers, which leads to increased diastolic filling. The importance of increased compliance is further seen in the findings that increasing the neo-Hookean constant decreased SV. One important takeaway from the parametric studies is that there appears to be an ideal stiffness at which diastolic filling is increased but systolic contraction is not greatly impaired. This is seen in the results for the infract thickness, Young’s modulus, and initial crimp angle. It is also interesting that the fitted values we found for the Young’s modulus and the crimp angle of porcine scar tissue appear to be near these ideal stiffness regions. The finding that increased fiber dispersion is beneficial for circumferentially aligned fibers but detrimental for longitudinally aligned fibers further highlights the importance of scar anisotropy in preserving ventricular function.

This work was strongly motivated by the work from the Holmes group [24]. They studied the effect of scar anisotropy on LV function using a finite element model with mechanical properties of both the healthy and infarct regions described by a Fung-type constitutive equation. By comparing isotropic scars, longitudinally stiff scars, and circumferentially stiff scars, they found that SV would be maximized when the scar tissue was stiffer in the longitudinal direction than in the circumferential direction. Fomovsky’s work identified the direct role that infarct stiffness plays in controlling LV mechanics so as to inform the design of infarct restraint or stiffening therapies. Our current study developed a simple model to relate cardiac function to the features of the scar extracellular matrix rather than simply to stiffness through the use of a microstructure based constitutive model. The agreements that we see between our results and the results of the more complex finite element model presented by Fomovsky, in terms of increased longitudinal stiffness leading to increased SV give us confidence in our findings. A more recent study from the Holmes group has demonstrated that regional mechanics control collagen alignment [11]. In complement, we have demonstrated here that collagen alignment affects the scar mechanical properties and thus LV function. The model we have presented could be a simple, easy to use, analytical tool for studying how initial collagen alignment controls remodeling and infarct expansion through mechanical feedback.

While this model has offered much insight into the importance of collagen alignment in determining heart function, there are several limitations to the study. This model assumes a thin walled geometry which ignores the effects of transmural variation in stresses. While the choice of the thin-walled assumption may create additional error in the calculated stresses, it was deemed acceptable since the focus of the study was on identifying potential mechanisms rather than conducting detailed local stress analysis. Based on the fact that infarct scars are thinner than normal LV wall, we chose to model the infarct as a single layer with a disperse fiber alignment. In reality the infarct is composed of several layers of collagen fibers with a mean fiber orientation that varies through the wall thickness similar to the helical pattern of myocyte alignment in normal myocardium. This simplification in itself should have little effect on the results due to the thin wall assumption of the model. Further, the choice of a single fiber distribution simplified the fitting of the scar properties to experimental data. This model also neglects the effects of ventricular torsion, which plays an important role in systolic ejection [30].

In addition, the assumption that adjacent regions share deformations along the shared edges allowed for a simple method of enforcing material compatibility. On the downside, this means that the stress field is non-continuous. This assumption was made to capture the phenomenon of infarct coupling whereby increasing the stiffness of the infarct reduces the deformation of the neighboring regions. In reality this effect is much smaller in regions that are far from the infarct region than in the small border zones that surround the infarct with properties that transition from scar to normal, making it possible that the magnitude of the fiber orientation effects are slightly overstated by the model. Given that the properties in the normal regions should be independent of the infarct properties we believe that the mechanisms through which scar architecture influence SV identified by this model are correct. In whole we believe that the model is not oversimplified as it matches behavior reported both by Fomovsky, (improved SV with increased longitudinal stiffness) and by Bogen (reduced stiffness allowing for greater diastolic filling and contraction for small infarcts) [4, 24].

The ground-up modeling approach we have chosen will allow us to refine the model in future work. While the collagen network is the primary determinant of infarct mechanics in the fibrotic scar, there are other structural proteins such as elastin and ECM modifications such as cross-linking that also play a role in LV mechanics [2]. Refinements of the constitutive model to include the effects of these ECM features would be an interesting direction for future work. In addition, the ischemic LV is highly inhomogeneous and further effort to model the border zones as well as the transmural variations in infarct properties would add to the physiological relevance of the model. Another future direction could examine the mechanics of the LV during the ischemic or inflammatory phases occurring immediately following MI. During these phases mechanical properties of the infarct are controlled less by collagen and more by the turnover of existing ECM and intracellular structural proteins such as titin as well as edema. The development of the collagen network and in turn the mechanical properties of the scar following this initial period are highly controlled by the expression of growth factors such as TGF-β and proteases such as MMPs [3, 31]. Coupling the mechanical model to a dynamic model of protein regulation is one potential direction for future research. Another dynamic process that could be modeled in future work is the change in the orientation of the fibers that occur under load, which can dynamically change the mechanical properties of the tissue [32, 33].

The post-MI LV is a complex and dynamic physiological environment and much work is ongoing to create highly accurate models. The goal of our paper was not to improve upon the accuracy of the existing numerical models, but to instead clarify the role that collagen plays in controlling infarct mechanics. Compared to some of the complex finite element models found in the literature, the simplicity of this model is actually an advantage as it focuses attention to the most important variables. This makes it a quick tool for identifying the key structural determinants of LV function post-MI. For example, the current simplicity of the model may make it a good choice for investigating many of the temporal changes that occur post-MI.

In summary, our model results showed that EF and SV are both optimized when fibers are aligned longitudinally, which further supports findings in the literature. The finding of an ideal fiber angle may aid in the development of diagnostic indices and treatment strategies based on determining fiber angle (perhaps through angle-sensitive MRI [34]) and directing it towards an optimal orientation through controlled remodeling. Since the model is based on the mechanical properties of actual infarct tissue it will also serve as an excellent starting point for creating multi-scale models that explore the interactions of genome, proteome, cell, tissue, and organ level changes in the post-MI remodeling process.

Conclusions

A simple model of the post-MI LV was developed to demonstrate that a longitudinal alignment of collagen fibers in the scar region maximizes cardiac function. This improvement is the result of increased compliance in the circumferential direction which allows for greater diastolic filling and greater systolic contraction through the Frank-Starling mechanism.

Appendix

1. Derivation of equilibrium equations

Many of the basic continuum mechanics equations used below can be found in Humphrey’s Cardiovascular Solid Mechanics text [25]. As described in the paper, there are six unknown deformation components. In order to solve for these six deformation components it is necessary to derive six independent equilibrium equations to relate stress to deformation. The internal forces in each region that balance the intraventricular pressure can be calculated from the Cauchy stress given in Eq. (9). If a matrix [α] is defined with rows ${\displaystyle n_i^T}$a _i and a matrix [ξ] is defined with rows q _i then Eq. (9) can be written as

$$\left[\mathit{\xi }\right]=\left[\mathit{\alpha }\right]\cdot \left[\mathit{\sigma }\right].$$

(A1)

Similarly a matrix [A] can be defined with rows N _i A _i representing the areas in the undeformed configuration and a matrix [Ξ] with columns Q _i representing the forces in the undeformed configuration. Both [Ξ] and [A] can be related to [ξ] and [α] respectively through the deformation gradient.

$$\left[\mathit{\xi }\right]=\left[\mathit{\Xi }\right]\cdot {\left[\mathit{F}\right]}^T,$$

(A2)

$$\left[\mathit{\alpha }\right]=\left[\mathit{{\rm A}}\right]\cdot {\left[\mathit{F}\right]}^{-1}.$$

(A3)

In the undeformed configuration, it is assumed that the surface normals N _i are aligned with the circumferential, longitudinal, and radial axes and that the internal forces in the undeformed configuration Q _i are normal to the undeformed surface normals N _i since pressure acts normal to surfaces. Thus,

$$\left[\mathit{\Xi }\right]=\left[\begin{array}{ccc}\hfill Q_{\mathit{\theta \theta }}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill Q_{\mathit{zz}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right],$$

(A4)

and

$$\left[\mathit{{\rm A}}\right]=\left[\begin{array}{ccc}\hfill \mathit{HT}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathit{LT}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathit{LH}\hfill \end{array}\right],$$

(A5)

where, H is the undeformed, longitudinal height of the region, L is the undeformed, circumferential length, and T is the thickness of the region. Now solving Eq. (A2) in matrix form,

$$\left[\mathit{\xi }\right]=\left[\begin{array}{ccc}\hfill Q_{\mathit{\theta \theta }}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill Q_{\mathit{zz}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]\left[\begin{array}{ccc}\hfill {\lambda }_{\theta }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill {\kappa }_{\mathit{\theta z}}\hfill & \hfill {\lambda }_z\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{{\lambda }_{\theta }{\lambda }_z}\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill Q_{\mathit{\theta \theta }}{\lambda }_{\theta }\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill Q_{\mathit{zz}}{\kappa }_{\mathit{\theta z}}\hfill & \hfill Q_{\mathit{zz}}{\lambda }_z\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill q_{\mathit{\theta \theta }}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill q_{\mathit{z\theta }}\hfill & \hfill q_{\mathit{zz}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right],$$

(A6)

it is seen that there are only 3 non-zero internal force components of [ξ]. Eq. (A3) can also be solved in matrix form for [α]

$$\left[\mathit{\alpha }\right]=\left[\begin{array}{ccc}\hfill \mathit{HT}\hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \mathit{LT}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathit{LH}\hfill \end{array}\right]{\left[\begin{array}{ccc}\hfill {\lambda }_{\theta }\hfill & \hfill {\kappa }_{\mathit{\theta z}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\lambda }_z\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \frac{1}{{\lambda }_{\theta }{\lambda }_z}\hfill \end{array}\right]}^{-1}=\left[\begin{array}{ccc}\hfill \frac{\mathit{HT}}{{\lambda }_{\theta }}\hfill & \hfill \frac{-\mathit{HT}{\kappa }_{\mathit{\theta z}}}{{\lambda }_{\theta }{\lambda }_z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\mathit{LT}}{{\lambda }_z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathit{LH}{\lambda }_{\theta }{\lambda }_z\hfill \end{array}\right].$$

(A7)

Substituting Eq. (A7) into Eq. (A1) and assuming negligible radial stresses yields

$$\left[\mathit{\xi }\right]=\left[\begin{array}{ccc}\hfill \frac{\mathit{HT}}{{\lambda }_{\theta }}\hfill & \hfill \frac{-\mathit{HT}{\kappa }_{\mathit{\theta z}}}{{\lambda }_{\theta }{\lambda }_z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill \frac{\mathit{LT}}{{\lambda }_z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \mathit{LH}{\lambda }_{\theta }{\lambda }_z\hfill \end{array}\right]\left[\begin{array}{ccc}\hfill {\sigma }_{\mathit{\theta \theta }}\hfill & \hfill {\sigma }_{\mathit{\theta z}}\hfill & \hfill 0\hfill \\ \hfill {\sigma }_{\mathit{\theta z}}\hfill & \hfill {\sigma }_{\mathit{zz}}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right]=\left[\begin{array}{ccc}\hfill \left({\sigma }_{\mathit{\theta \theta }}-{\sigma }_{\mathit{\theta z}}\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}\right)\frac{\mathit{HT}}{{\lambda }_{\theta }}\hfill & \hfill \left({\sigma }_{\mathit{\theta z}}-{\sigma }_{\mathit{zz}}\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}\right)\frac{\mathit{HT}}{{\lambda }_{\theta }}\hfill & \hfill 0\hfill \\ \hfill {\sigma }_{\mathit{\theta z}}\frac{\mathit{LT}}{{\lambda }_z}\hfill & \hfill {\sigma }_{\mathit{zz}}\frac{\mathit{LT}}{{\lambda }_z}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill 0\hfill \end{array}\right].$$

(A8)

Relating the result of Eq. (A8) with the components of [ξ] as given by Eq. (A6) provides a system of 4 equations

$$q_{\mathit{\theta \theta }}=\left({\sigma }_{\mathit{\theta \theta }}-{\sigma }_{\mathit{z\theta }}\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}\right)\frac{\mathit{HT}}{{\lambda }_{\theta }},$$

(A9.a)

$$0=\left({\sigma }_{\mathit{\theta z}}-{\sigma }_{\mathit{zz}}\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}\right)\frac{\mathit{HT}}{{\lambda }_{\theta }},$$

(A9.b)

$$q_{\mathit{z\theta }}={\sigma }_{\mathit{\theta z}}\frac{\mathit{LT}}{{\lambda }_z},$$

(A9.c)

$$q_{\mathit{zz}}={\sigma }_{\mathit{zz}}\frac{\mathit{LT}}{{\lambda }_z}.$$

(A9.d)

From Eq. (A9.b), it is seen that

$${\sigma }_{\mathit{\theta z}}={\sigma }_{\mathit{zz}}\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}.$$

(A10)

Combining this fact with Eq. (A9.c) and Eq. (A9.d) yields the following relationship between force and deformation

$$\frac{q_{\mathit{z\theta }}}{q_{\mathit{zz}}}=\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}.$$

(A11)

Now that relationships between force, stress, and deformation have been established for individual regions it is necessary to examine the four region model with forces and deformations being shared as described in Figure 3. In order to determine the magnitudes of the internal force vectors described in Eq. (11) it is necessary to examine the components of these vectors again.

From Eq. (A6) it is seen that

$${\mathit{q}}_{\theta }=q_{\mathit{\theta \theta }}{\mathit{i}}_{\theta },$$

(A12.a)

$${\mathit{q}}_z=q_{\mathit{z\theta }}{\mathit{i}}_{\theta }+q_{\mathit{zz}}{\mathit{i}}_z,$$

(A12.b)

where i _θ and i _z are unit vectors. From Eq. (A12.a) it is seen that

$$||{\mathit{q}}_{\theta }||=q_{\mathit{\theta \theta }},$$

(A13)

which allows for Eq. (11.a) to be rewritten as

$$q_{\mathit{\theta \theta }}={q_{\mathit{\theta \theta }}}^I+{q_{\mathit{\theta \theta }}}^{\mathit{III}}={q_{\mathit{\theta \theta }}}^{\mathit{II}}+{q_{\mathit{\theta \theta }}}^S=\mathit{Prh}.$$

(A14)

Substituting Eq. (A9.a) for each region into Eq. (A14) yields 2 equilibrium equations

$$\left({{\sigma }_{\mathit{\theta \theta }}}^I-{{\sigma }_{\mathit{z\theta }}}^I\frac{{{\kappa }_{\mathit{\theta z}}}^I}{{{\lambda }_z}^I}\right)\frac{H^I{T}^I}{{{\lambda }_{\theta }}^I}+\left({{\sigma }_{\mathit{\theta \theta }}}^{\mathit{III}}-{{\sigma }_{\mathit{z\theta }}}^{\mathit{III}}\frac{{{\kappa }_{\mathit{\theta z}}}^{\mathit{III}}}{{{\lambda }_z}^{\mathit{III}}}\right)\frac{H^{\mathit{III}}{T}^{\mathit{III}}}{{{\lambda }_{\theta }}^{\mathit{III}}}=\mathit{Prh},$$

(A15.a)

$$\left({{\sigma }_{\mathit{\theta \theta }}}^{\mathit{II}}-{{\sigma }_{\mathit{z\theta }}}^{\mathit{II}}\frac{{{\kappa }_{\mathit{\theta z}}}^{\mathit{II}}}{{{\lambda }_z}^{\mathit{II}}}\right)\frac{H^{\mathit{II}}{T}^{\mathit{II}}}{{{\lambda }_{\theta }}^{\mathit{II}}}+\left({{\sigma }_{\mathit{\theta \theta }}}^S-{{\sigma }_{\mathit{z\theta }}}^S\frac{{{\kappa }_{\mathit{\theta z}}}^S}{{{\lambda }_z}^S}\right)\frac{H^S{T}^S}{{{\lambda }_{\theta }}^S}=\mathit{Prh}.$$

(A15.b)

The vector q _z can be written in polar notation as

$${\mathit{q}}_z=||{\mathit{q}}_z||\sin \left(\phi \right){\mathit{i}}_{\theta }+||{\mathit{q}}_z||\cos \left(\phi \right){\mathit{i}}_z.$$

(A16)

where

$$\phi =\mathrm{atan}\left(\frac{q_{\mathit{z\theta }}}{q_{\mathit{zz}}}\right).$$

(A17)

Substituting Eq. (A11) into Eq. (A17) yields

$$\phi =\mathrm{atan}\left(\frac{{\kappa }_{\mathit{\theta z}}}{{\lambda }_z}\right).$$

(A18)

Equating Eq. (A16) with Eq. (A12.b) gives

$$q_{\mathit{z\theta }}=||{\mathit{q}}_z||\sin \left(\phi \right),$$

(A19.a)

$$q_{\mathit{zz}}=||{\mathit{q}}_z||\mathit{cos}\left(\phi \right),$$

(A19.b)

Solving Eq. (11.b) for the magnitude and equating with Eq. (A19) yields

$$||{\mathit{q}}_z||=\frac{q_{\mathit{z\theta }}}{\sin \left(\phi \right)}=\frac{{q_{\mathit{z\theta }}}^I}{\sin \left({\phi }^I\right)}+\frac{{q_{\mathit{z\theta }}}^{\mathit{II}}}{\sin \left({\phi }^{\mathit{II}}\right)}=\frac{{q_{\mathit{z\theta }}}^{\mathit{III}}}{\sin \left({\phi }^{\mathit{III}}\right)}+\frac{{q_{\mathit{z\theta }}}^S}{\sin \left({\phi }^{\mathit{IV}}\right)}=\mathit{\pi P}{r}^2,$$

(A20.a)

$$||{\mathit{q}}_z||=\frac{q_{\mathit{zz}}}{\cos \left(\phi \right)}=\frac{{q_{\mathit{zz}}}^I}{\cos \left({\phi }^I\right)}+\frac{{q_{\mathit{zz}}}^{\mathit{II}}}{\cos \left({\phi }^{\mathit{II}}\right)}=\frac{{q_{\mathit{zz}}}^{\mathit{III}}}{\cos \left({\phi }^{\mathit{III}}\right)}+\frac{{q_z}^S}{\cos \left({\phi }^{\mathit{IV}}\right)}=\mathit{\pi P}{r}^2.$$

(A20.b)

Substitution of Eq. (A9.c) for each region into Eq. (A20.a) gives 2 more equilibrium equations

$${{\sigma }_{\mathit{\theta z}}}^I\frac{L^I{T}^I}{{{\lambda }_z}^I\sin \left({\phi }^I\right)}+{{\sigma }_{\mathit{\theta z}}}^{\mathit{II}}\frac{L^{\mathit{II}}{T}^{\mathit{II}}}{{{\lambda }_z}^{\mathit{II}}\sin \left({\phi }^{\mathit{II}}\right)}=\mathit{\pi P}{r}^2,$$

(A21.a)

$${{\sigma }_{\mathit{\theta z}}}^{\mathit{III}}\frac{L^{\mathit{III}}{T}^{\mathit{III}}}{{{\lambda }_z}^{\mathit{III}}\sin \left({\phi }^{\mathit{III}}\right)}+{{\sigma }_{\mathit{\theta z}}}^S\frac{L^S{T}^S}{{{\lambda }_z}^S\sin \left({\phi }^S\right)}=\mathit{\pi P}{r}^2.$$

(A21.b)

Substitution of Eq. (A9.d) for each region into Eq. (A20.b) gives the last 2 equilibrium equations

$${{\sigma }_{\mathit{zz}}}^I\frac{L^I{T}^I}{{{\lambda }_z}^I\cos \left({\phi }^I\right)}+{{\sigma }_{\mathit{zz}}}^{\mathit{II}}\frac{L^{\mathit{II}}{T}^{\mathit{II}}}{{{\lambda }_z}^{\mathit{II}}\mathit{cos}\left({\phi }^{\mathit{II}}\right)}=\mathit{\pi P}{r}^2,$$

(A22.a)

$${{\sigma }_{\mathit{zz}}}^{\mathit{III}}\frac{L^{\mathit{III}}{T}^{\mathit{III}}}{{{\lambda }_z}^{\mathit{III}}\cos \left({\phi }^{\mathit{III}}\right)}+{{\sigma }_{\mathit{zz}}}^S\frac{L^S{T}^S}{{{\lambda }_z}^S\cos \left({\phi }^S\right)}=\mathit{\pi P}{r}^2.$$

(A22.b)

Collecting Eq. (A15), Eq. (A21), and Eq. (A22) and noting the shared lengths, heights, and deformations as shown in Figure 1 produces the final set of equilibrium equations.

$$\left({{\sigma }_{\mathit{\theta \theta }}}^I-{{\sigma }_{\mathit{z\theta }}}^I\frac{{{\kappa }_{\mathit{\theta z}}}^R}{{{\lambda }_z}^R}\right)\frac{H^R{T}^I}{{{\lambda }_{\theta }}^R}+\left({{\sigma }_{\mathit{\theta \theta }}}^{\mathit{III}}-{{\sigma }_{\mathit{z\theta }}}^{\mathit{III}}\frac{{{\kappa }_{\mathit{\theta z}}}^S}{{{\lambda }_z}^S}\right)\frac{H^S{T}^{\mathit{III}}}{{{\lambda }_{\theta }}^R}=\mathit{Prh},$$

(A23.a)

$$\left({{\sigma }_{\mathit{\theta \theta }}}^{\mathit{II}}-{{\sigma }_{\mathit{z\theta }}}^{\mathit{II}}\frac{{{\kappa }_{\mathit{\theta z}}}^R}{{{\lambda }_z}^R}\right)\frac{H^R{T}^{\mathit{II}}}{{{\lambda }_{\theta }}^S}+\left({{\sigma }_{\mathit{\theta \theta }}}^S-{{\sigma }_{\mathit{z\theta }}}^S\frac{{{\kappa }_{\mathit{\theta z}}}^S}{{{\lambda }_z}^S}\right)\frac{H^S{T}^S}{{{\lambda }_{\theta }}^S}=\mathit{Prh},$$

(A23.b)

$${{\sigma }_{\mathit{\theta z}}}^I\frac{L^R{T}^I}{{{\lambda }_z}^R}+{{\sigma }_{\mathit{\theta z}}}^{\mathit{II}}\frac{L^S{T}^{\mathit{II}}}{{{\lambda }_z}^R}=\mathit{\pi P}{r}^2\sin \left({\phi }^R\right),$$

(A23.c)

$${{\sigma }_{\mathit{\theta z}}}^{\mathit{III}}\frac{L^R{T}^{\mathit{III}}}{{{\lambda }_z}^S}+{{\sigma }_{\mathit{\theta z}}}^{\mathit{IV}}\frac{L^S{T}^S}{{{\lambda }_z}^S}=\mathit{\pi P}{r}^2\sin \left({\phi }^S\right).$$

(A23.d)

$${{\sigma }_{\mathit{zz}}}^I\frac{L^R{T}^I}{{{\lambda }_z}^R}+{{\sigma }_{\mathit{zz}}}^{\mathit{II}}\frac{L^S{T}^{\mathit{II}}}{{{\lambda }_z}^R}=\mathit{\pi P}{r}^2\cos \left({\phi }^R\right),$$

(A23.e)

$${{\sigma }_{\mathit{zz}}}^{\mathit{III}}\frac{L^R{T}^{\mathit{III}}}{{{\lambda }_z}^S}+{{\sigma }_{\mathit{zz}}}^S\frac{L^S{T}^S}{{{\lambda }_z}^S}=\mathit{\pi P}{r}^2\cos \left({\phi }^S\right),$$

(A23.f)

where

$${\phi }^R=\mathrm{atan}\left(\frac{{{\kappa }_{\mathit{\theta z}}}^R}{{{\lambda }_z}^R}\right),$$

(A24.a)

$${\phi }^S=\mathrm{atan}\left(\frac{{{\kappa }_{\mathit{\theta z}}}^S}{{{\lambda }_z}^S}\right).$$

(A24.b)

This set of 6 independent equilibrium equations can be used in conjunction with the constitutive equations to solve for the 6 independent deformation components, ${\lambda }_{\theta }^R,{\lambda }_{\theta }^S,{\lambda }_z^R,{\lambda }_z^S,{\kappa }_{\mathit{\theta z}}^R\mathit{and}{\kappa }_{\mathit{\theta z}}^S$.