Coupled agent‐based and hyperelastic modelling of the left ventricle post‐myocardial infarction

Abstract Understanding the healing and remodelling processes induced by myocardial infarction (MI) of the heart is important, and the mechanical properties of the myocardium post‐MI can be indicative for effective treatments aimed at avoiding eventual heart failure. MI remodelling is a multiscale feedback process between the mechanical loading and cellular adaptation. In this paper, we use an agent‐based model to describe collagen remodelling by fibroblasts regulated by chemical and mechanical cues after acute MI, and upscale into a finite element 3D left ventricular model. We model the dispersed collagen fibre structure using the angular integration method and have incorporated a collagen fibre tension‐compression switch in the left ventricle (LV) model. This enables us to study the scar healing (collagen deposition, degradation, and reorientation) of a rat heart post‐MI. Our results, in terms of collagen accumulation and alignment, compare well with published experimental data. In addition, we show that different shapes of the MI region can affect the collagen remodelling, and in particular, the mechanical cue plays an important role in the healing process.


INTRODUCTION
Myocardial infarction (MI), commonly known as a heart attack, occurs when blood flow decreases or stops to a part of the heart, causing damage to the heart muscle. The healing process in the heart poses a complex multiscale soft tissue problem, involving cardiac growth and remodelling (G&R). Cellular growth is often considered to be the cause for residual stress. 1 The myocardial stress and constitutive properties are considered to be the two key factors in myocardial G&R. 2 There are typically two types of G&R modelling approaches at the continuum level. One is the volumetric growth (the density remains unchanged) following the G&R, and the other is the density growth (volume remains unchanged). In the volumetric approach, G&R can be modelled using the growth tensor first introduced by Rödriguez et al. 3 For example, Göketepe et al 4 developed a multiscale framework to study myocardial growth via alignment of myocytes and parallel additional of sarcomere units. Their model was extended to study the mechanical effects on remodelling patterns of the left ventricle (LV) after certain surgical procedures. 5,6 Kerckhoffs et al 7 developed a strain-driven growth law to explain how the grown cardiac myocytes increased the fibre and cross-fibre strains. 8 Density growth, on the other hand, focuses on changing the constitutive properties of myocaridum in healthy remote and infarcted cardiac tissues, the latter being characterised with stiffer constitutive properties, simulating the replacement of fibrosis by collagen fibres. [9][10][11] For example, in a study, 12 the stiffness in the infarct zone is increased significantly compared with functional myocardium in a patient-specific infarcted LV model.
There are very few combined volumetric and density growth models for G&R, particularly for the heart. Eriksson et al 13 developed a framework to study the G&R process of fibre-reinforced arteries, while the evolution of constitutive properties was simulated via a density-and-volumetric growth approach for a mixture of constituents. However, the remodelling of the fibre structure was not included. Indeed, it is not clear how the collagen fibre structure is modified during the post-MI healing process, and how the macromechanical response of the infarcted LV changes post-MI.
Recent developments in soft tissue modelling in parallel with advances in experimental techniques have enabled us to predict mechanical states at, and obtain images of, the small scales relevant to biology. For example, the structure-based fibre-reinforced Holzapfel-Ogden (HO) model 14 for myocardium was developed based on simple shear experiments and has catalysed much of the recent activities in soft tissue modelling.
At the cellular level, collagen deposition and remodelling are regulated by fibroblast cell alignment. Environmental cues, such as mechanical and chemical cues, have been shown to influence cell migration and regulate the microfibre structures. 15 Recently, agent-based models that account for these effects have been developed and used to study a two-dimensional (2D) slab model of the myocardium infarction. 16,17 However, the extension of this approach to a three-dimensional (3D) LV model has not been reported.
In this paper, we combine an agent-based approach and a structure-based fibre-reinforced constitutive law to study MI healing in a 3D LV model for the first time. To account for the necessarily dispersed fibre structure post-MI, we modify the original HO constitutive law 14 by employing a distributed fibre model. 18 The specific fibre distribution is determined using an agent-based model similar to that of Fomovsky & Holmes. 16 We then incorporate the agent-based model within a finite element (FE) LV model. This new mathematical model is used to simulate the myocardium remodelling in terms of the collagen fibre structure and density. To account for the microfibre structure in the tissue constitutive laws, a commonly used up-scaling method is based on volumetric averaging. 19 However, since these microfibres in the soft tissue do not support compression, a switch is required to exclude the compressed fibres. For a dispersed fibre structure, however, such a "switch" is often applied inappropriately due to the volume averaging. 18 In our model, we describe the fibre structure using angular integration with an analytical expression for the fibre on-off switch, which is then implemented in the constitutive model.

METHODOLOGY
Here, we describe the model in detail. For convenience, the list of the variables used in the model is defined in Table 1.

The geometry of the LV model
An idealised half ellipsoid geometry is used to construct a FE model for a rat LV, with 2100 eight-node hexahedral elements and 2375 nodes ( Figure 1). Global Cartesian coordinates (X, Y, Z) are used to describe material points in the undeformed reference configuration, with the corresponding basis vectors denoted Myocardium is considered to be a fibre-reinforced material mostly composed of collagen fibres and myocytes. 14 A local coordinate system with the circumferential, longitudinal, and transmural basis vectors (c 0 , l 0 , n 0 ), is introduced to describe the layered fibre structure within the ventricular wall, as shown in Figure 1. Note that The myofibre architecture is generally described by a "fibre-sheet-normal" system (f, s, sn) 21 in the current configuration. Here, we assume that the fibre direction f always lies in the c − l plane, the sheet direction is transmural, and the sheet-normal sn = f × s, where c and l are the current circumferential and longitudinal directions: c = Fc 0 , l = Fl 0 , and F is the deformation gradient.

Agent-based model 2.2.1 Chemokine concentration
Fibroblast cells adjust the microstructure of heart tissue, including the density and orientation of the collagen fibre bundles, which determine the material properties of infarcted tissues. The collagen fibres are organized in a highly layered architecture in healthy tissues. Following an acute MI, sudden changes of the chemokine concentration and mechanical environment in and around the infarct area activate the fibroblasts, which will remodel the micro-collagen structure of heart tissues by changing the collagen orientations and its volume fraction through deposition and degradation. Subsequently, the material properties and mechanical behaviour are modified at the tissue level. In this study, we extend the 2D agent-based model developed in Fomovsky and Holmes 16 to 3D, to describe the collagen remodelling post-MI. The modified 3D agent-based model is explained in the coming section. We assume a spherical MI with its centre at X c . The static chemokine concentration at point X, representing the milieu of cytokine and chemokine surroundings, is described by a chemical diffusion equation as where Ω b is the volume of the infarct region, D c is the diffusion coefficient, C is the chemokine concentration, and k c,gen , k c,deg are the chemokine generation and degradation rates, respectively. The boundary condition for the chemokine concentration equation is Experimental data suggest that the infarct-induced chemical concentration reduces rapidly from the infarct area, and drops to nearly zero away from the infarct centre. Hence, in practice, we choose C(X) = 0 when ||X − X c || is 10 times of the infarct zone in order to approximate diffusion into an infinite space, following Rouillard and Holmes. 17 Continuity of the chemokine field requires that For a spherical and 3D infarct zone of radius r 0 , (2) has an analytical solution, where a 0 = k c,gen ∕k c,deg , The fibroblast activation parameters, including cell migration speed, collagen degradation rate, collagen deposition rate, and collagen reorientation rate, are all modulated by the local chemokine concentration C(r). To compute the activation parameters in the agent-based model, the local fibroblast activation parameters are assumed to vary linearly with C as where P i is the ith rate parameter, P i,max and P i,min are the maximum and minimum rates, and C max , C min are the maximum and minimum chemokine concentrations.

Fibroblast migrations regulated by environmental cues
We assume fibroblasts are regulated by local environmental cues. 16 These include the chemical, mechanical, persistence, and structural cues. [15][16][17] We further assume the fibroblasts are rigid spheres with radius R cell , and there is no interaction between fibroblasts. We also use Θ to denote the angles of the fibroblast cells at time t, and to indicate the collagen fibre angles.

Chemical cues
We define the chemokine vector as the product of the chemokine concentration and the outward normal vector at the fibroblast boundary. The chemical cue v c is obtained from where n c (Θ) is the unit outward normal vector of the cell surface in the current configuration, r is the distance from the boundary of a fibroblast cell to the infarction centre X c .

Mechanical cues
The mechanical cue v m is defined as where n = n c · Dn c is the normal strain, is the Almansi strain tensor, and I is the identity matrix.

Structural cues
The structural cue v s is defined as where N( ) is the number of collagen fibres at angle , f( ) is the unit vector in the current configuration defined by f = Ff 0 ∕||Ff 0 ||, where f 0 = cos c 0 + sin l 0 is the fibre vector in the reference configuration.

Persistence cue
The persistence cue v p is defined as the fibroblast cell migration velocity where U(Θ) is the unit vector along Θ in the current configuration, given by The fibroblast position at time t is tracked from

Resultant cue
The fibroblast resultant cue ( Figure 2) is a weighted and normalised sum of all the cues (7) to (10), ie, where v i is the ith cue vector, c i is the ith normalised and weighted cue vector, M i and W i are ith scaling and weighting factors for the ith cue as defined in Fomovsky and Holmes, 16 and is the persistence tuning factor. The orientation of provides the mean fibre angleΘ. We assume that the fibroblasts obey a von Mises distribution where I 0 is the Bessel function of the first kind of order zero

Remodelling of the collagen fibre structure
The local collagen fibre distribution in the current configuration is remodelled by fibroblasts in the following way.
1. Deposition and degradation: The total collagen fibre number at time t is The rate of change of the number N( ) of collagen fibres in the direction is determined from where is the maximum collagen fibre number per unit area, R cell is the radius of the fibroblast cell, is the delta function to ensure that new collagen fibres are aligned with the fibroblast migration direction, and k cf,gen and k cf,deg are the rates of generation and degradation of collagen fibres, respectively. 2. Rotation: The angle of a collagen fibre changes according to where k cf,rot is the rate of rotation estimated from (6).
We assume that the material properties are the same for f and −f, so we only need to consider values of ∈ ( − ∕2, ∕2). To achieve this, the following transformation is used: The deformed and remodelled fibre structure needs to be pushed back to the reference configuration, 22,23 ie,

Upscaling the fibre structure to the tissue level
The updated volume fraction Φ cf,t is where Φ 0 cf and N 0 tot are the volume fraction and numbers of the fibres in the reference configuration. The fibre orientation density function, when pushed back to the reference configuration, is estimated from In the initial reference configuration, we assume that the fibres have -periodic von Mises distribution 24,25 : where 0 is the initial concentration parameter, fitted from experimental observation, 17 and̄0 is the initial mean angle of the collagen fibre structure.

Constitutive laws for myocardium 2.3.1 Modified HO model with fibre orientation density function
For the remote healthy myocardium (r(X) > r 0 ), we use the standard HO model, 14 ie, where in which a, b, a cf , b cf are material parameters, and I 1 , I 4 are the invariants of the right Cauchy-Green tensor C = F T F, Here, ":" denotes a double contraction.
To describe the mechanical behaviour of the infarcted tissue (r(X) ≤ r 0 ), we modify the HO model as where Φ m,t ( = 1 − Φ cf,y ) is the volume fraction of the ground matrix, Ψ m is the same as in (22) but Ψ cf is changed to Notice that the integration in (25) is over a domain Σ in which the fibres are in tension, not ( − ∕2, ∕2). This is because collagen fibres bear load only when they are stretched. We will address this point in Section 2.3.2.
A rule-based myocardium fibre generation algorithm 21,26 is adopted to describe the local mean fibre angle. We first calculate the normalised distance parameter d for an arbitrary material point X as where d endo , d epi are the distances from X to the endocardial and epicardial surfaces. Then, the local mean fibre anglē0 at X is defined as̄0 where max = ∕3, because it has been observed that the mean angle of fibres in a healthy left ventricle rotates transmurally across the heart wall from /3 at the endocardial surface to − /3 at the epicardial surface. 27

Exclusion of the compressed fibres
To exclude the compressed fibres, we choose Σ in (25) so that where C ij are the Cartesian components of C. The inequality (27) gives the range of Σ for the stretched collagen fibres in the following two cases, as exemplified in Figure 3.
where B = FF T is the left Cauchy-Green deformation tensor. All the constitutive parameters are given in Table 2.

FIGURE 3
The shaded areas show the range of Σ of stretched collagen fibres within − 2 < < 2 for selected scenarios (dash line denotes

Change of basis from Cartesian to local coordinates
It is convenient to transform the description of fibre orientation between the global Cartesian coordinates and the local polar basis vectors, so that right Cauchy-Green tensor are given as where Q is the rotation tensor from Cartesian basic vectors to the local polar basis vectors in the reference configuration.

The coupled agent-based and FE model
The LV diastolic dynamics is simulated using a coupled agent-based FE LV model. We assume that mechanical cues are based on the end-diastolic state and the remodelling only occurs in the MI region. The FE model is solved by using the FE software FEAP, 29 in which the incompressibility is enforced with a penalty method. 30 Specifically, we consider the myocardium as a quasi-incompressible material. Using the multiplicative decomposition of the deformation gradient F, where J = det(F), we write the strain energy function Ψ as whereC =F TF , andΨ m ,Ψ cf are isochoric contributions of the matrix and fibres, respectively, expressed in terms of the modified invariants,Ī 1 (=C ∶ I) andĪ 4 (= f 0 ·Cf 0 ), The volumetric penalty term is chosen to be The Cauchy stress is then whereb =FF T , and the operator dev(•) is defined as The Hu-Washizu mixed variational approach is used to solve the FE equations to avoid locking, 29 as briefly listed in Appendix A. The computational flowchart is shown in Figure 4.

MI healing case studies
To validate the model, we follow the experimental study of the infarcted LV of rats by Fomovsky et al, 31 in which cryoinfarctions were created by sewing steel cylinders filled with liquid nitrogen into tissues under the epicardial surfaces of rat ventricles. After weeks of healing, the hearts were harvested to study the collagen structures at the surfaces of the infarcted zones. We select two cases from Fomovsky et al 31 to simulate. In the first case, a transmural circular MI is  Figure 5A, and in the second case, a transmural elliptical MI is as shown in Figure 5B, each in the mid-ventricular wall. We modelled the healing process for 5 weeks while the mature scar was formed and the ventricular filling pressure increased. 31 A linearly increasing end-diastolic pressure profile is assumed, with values of 12 mmHg at time 0 to 18 mmHg at 5 weeks. The LV base was constrained in the longitudinal direction, but free movements in the radial and circumferential directions were allowed. The model parameters are listed in Table 2, taken from Fomovsky et al. 31

Evolution of the fibre structure post-MI
The simulated collagen accumulations for both the circular and elliptical MI shapes in the infarct centre at 0, 1, 2, and 5 weeks post-MI are shown in Figure 6. The volume fraction increases from 3% to around 28%. This agrees well with the measured data. 16 We note that the difference of the volume fraction increase is very little between the different MI shapes (28.2% for circular MI vs 28.0% for elliptical MI after 5 weeks). The distributions of the collagen fibres at different weeks post-MI for the circular case are illustrated in Figure 7. At 5 weeks post-MI, the simulated result is also compared with the corresponding measurement. 16 Given the large error bars in the measurements, there is a reasonable agreement in terms of the overall distribution and the maximum collagen fibre fraction predicted. Figure 8A shows that for both MI shapes, the mean angle of collagen fibre families at 20% of the distance from the epicardium to the endocardium surface decreases significantly at the fifth week (from the initial 45 • to 22.5 • for the circular MI, and to 10 • for the elliptical MI).
As the anisotropy in collagen fibre structure is a critical determinant for the pump function of LV, 17,31 we use the kurtosis of the collagen fibre structure to describe the anisotropy level, which is defined as where E(•) is the expectation of •, and̄is the mean fibre angle. A smaller value of kurtosis indicates a lower level of anisotropy. Figure 8B shows that the anisotropy level at the infarcted region decreases initially but eventually reaches a plateau. It is interesting to see that although the fibre volume fraction is similar for both the circular and elliptical MI shapes, the distribution of the fibres is quite different. The evolution of the fibre structure for the elliptical MI is faster than that for the circular MI. The anisotropy level of the fibre structure, in terms of the kurtosis, decreases from the initial value of 3.38 to 3.29 for the circular MI, and to 3.21 for the elliptical MI. In general, the elliptical MI causes greater  changes in terms of the mean fibre angle and kurtosis over the healing time. However, in both cases, the volume fraction increased to a stable level of around 28% ( Figure 6) at the fifth week. The decreased anisotropy of the fibre distribution over time is better shown in Figure 9, where the local fibre structure near the epicardial surface at different weeks post-MI is illustrated. Note that to make the fibre structure visible, we have used the averaged fibre vectors of small local clusters, each consisting of hundreds of the simulated fibre vectors.

Evolution of the stress and strain post-MI
The first principal stress distribution for the circular MI at different times post-MI is illustrated in Figure 10, where the infarcted region is indicated by an arrow. It is evident that the stress distribution is greatly influenced by the evolution of the local collagen fibre structure. The stress increases at the endocardial surface and gradually smears out towards the epicardial surface. The local stress level around the MI region is caused by collagen accumulation and reorientation, but the overall stress also increases due to the increased end-diastolic pressure. The evolution of the stress level with time is shown in Figure 11, and the corresponding mean stress values at the infarct centre are also listed in Table 3. During the healing process, the circumferential stiffness is reinforced by the newly deposited collagen fibres. Hence, the maximum circumferential stress increases dramatically from 4.3 to 15.1 kPa. The local longitudinal stress has not changed much (from 0.78 to 0.76 kPa at 5 weeks). This is because most of the collagen fibre remodelling happens in the circumferential direction. The elevated circumferential stress level around the MI region can lead to further adverse remodelling in the heart, such as myocyte apoptosis, because the local abnormal stress might cause cell degradation. 32 The stress distribution for the elliptical MI is similar to that of the circular MI, but the absolute stress levels are slightly different, as shown in Table 3.
Fomovsky et al 31 also measured the evolution of strain during the healing process in the epicardial surface of the infarcted region. They found that the circumferential strain in the infarct region continuously decreases with an average decrease of 65%, and an average decrease of 13% for the longitudinal strain at 6 weeks. Our simulation results give similar trends, as shown in Table 4 and in Figure 12 for both MI shapes. For example, in the case of the circular MI, the strain shows an average decrease of 49% in the circumferential strain (0.03 ± 0.01 at 0 week, 0.04 ± 0.024 at 1 week, 0.042 ± 0.01 at 3 weeks, 0.03 ± 0.02 at 4 weeks, and 0.022 ± 0.02 at 5 weeks), and an average decrease of 35% in the longitudinal strain at 5 weeks.
The evolution of the strains at the endocardial surface of the elliptical MI is similar to that of the circular MI ( Figure 12). Both circumferential and longitudinal strains decrease to the homeostatic state after the healing process, with 31% decrease in the circumferential strain and 43% decrease in the longitudinal strain.

FIGURE 12
Strain evolution over time: preferential accumulation in the circumferential direction reduces the anisotropy. The predicted strain trends agree reasonably with the experimentally measured strains 5 weeks after cryoinfarction 16

Influence of different cues
We now isolate the effect of different cues by switching off individual cues in turn. The distributions of the collagen fibres are shown in Figure 13. The mechanical cue was suggested to be the most important factor in regulating the collagen alignment. 17 This is confirmed in our study. When the mechanical cue is switched off, the collagen fraction in the infarcted region decreases more quickly, and the mean angle of the fibre structure at the endocardial surface remains 60 • , which is similar to that in healthy heart tissue. 14 In other words, with the mechanical cue, the mean angle of the fibre structure at the endocardial surface reduces from 60 • to 25 • , which is closer to the experimentally measured fibre value. 16 The effects of other cues are less prominent. For instance, the mean angles of the fibre structure at the endocardial surface change from 37 • to 25 • if the persistence cue is turned off; from 25 • to 21 • if the structure cue is turned off. The fact that the

DISCUSSION
Our model is the first to couple an agent-based model and a 3D idealised LV model with local MI and is a significant extension of the previous study based on a 2D FE slab. 34 In particular, a modified HO constitutive law is employed to simulate the mechanical behaviour of the myocardium post-MI, where fibre dispersion due to wound healing needs to be accounted for. To exclude the effects of the compressed fibres, an analytical fibre switch is implemented in this FE LV model. To use a model that includes the evolving collagen fibre structure post-MI requires one to estimate many material parameters at different length scales. In this study, the model parameters are taken from published cryoinfarction measurements, 16 and the predictions of our model showed good quantitative agreement. The only parameter not based on the measurements is the persistence tuning factor (Table 2). To see how much effect this has, we also ran the model with the value of increased. Our results (not shown) suggest that an increased leads to a greater standard deviation , reducing the randomness of the new fibroblast angles through (13). The fibre structure reached the steady state faster, and matured at around 4 weeks when was doubled, which is one week earlier than by using the default value of . The final fibre structures, however, are very similar.
Environmental cues, such as the chemical, mechanical and persistence cues, play important roles in regulating collagen fibre structures during MI healing. Our results show the poorest agreement with measurements if the mechanical cue is switched off. According to Zimmerman's experimental observation on the effect of initial collagen fibre structure, 35 the orientations of existing collagen fibres act as the template for alignments of new fibres during infarction healing process. This agrees with our model prediction when the structural cue is switched off. In this case, the collagen structure is more isotropic after scar formation; however, the remodelled mean angle of the fibre structure does not change much, with or without the structural cue. If the persistence cue is switched off, the distribution of the fibre structure is similar, but the shift of the mean angle is smaller, and the evolution is much slower than that in the experimental measurement. 16 McDougall et al 36 suggested that the distribution of the chemical gradient established by the chemokines generated around dermal wounds, determines patterns of local collagen alignment. In our simulations, different chemokine distributions are realised through the MI shapes; the chemokine field is mostly flat inside the infarction, but with a sudden drop near the infarct border. The influences of chemical gradient are greater in the region near the border. Therefore, the mean angle and kurtosis of the fibre distributions for different infarct shapes are different. Our simulations are also supported by the observations in Rouillard and Holmes. 17 Finally, we mention the limitations of this study. Firstly, the transmural fibroblast migration is not considered. Secondly, our model geometry is based on an idealised half ellipsoid and without active contraction. Therefore, the mechanical cue is simplified. Thirdly, we have not included volumetric growth in our model, and therefore, effects such as wall thinning are not included here. Furthermore, no transient chemokine concentration is included. Finally, we understand that MI in the human heart may involve different biological processes and parameters compared with cryoinfarction in the rat heart. Estimation of these parameters for patient-specific human hearts post-MI, therefore, remains a daunting task. These limitations may affect the quantitative findings presented in the paper because physiologically correct anatomy and microstructure affect the stress state in the heart and therefore the remodelling phase. However, we believe our agent-based approach for the time-dependent density growth and fibre distribution is a step forward towards cellular-based modelling of MI.

CONCLUSION
We have developed an agent-based MI model that takes into consideration collagen remodelling with a 3D LV model in diastole. In this model, the evolution of the microstructure of collagen fibres is simulated under the framework of an agent-based model based on the statistical behaviour of cellular movements. Then, the microstructure change is upscaled into the tissue level to describe the material properties of infarcted heart tissues. To describe the remodelled material properties of myocardium in the infarcted region, a collagen fibre tension-compression "switch" is incorporated in the FE LV model for the first time. The time-dependent model also captures the interaction and information exchange processes between the mechanical behaviour and collagen tissue remodelling guided by various external cues. The model results show similar trends in collagen accumulation and collagen alignment to those found in experimental measurements of infarct healing in the rat heart over 5 weeks. The shapes and the locations of the infarctions could affect the local collagen accumulation and the LV dynamics. For example, the mean angle of the fibre structure decreases from about 45 o to 22.5 o near the endocardial surface for the circular infarction but decreases less in the elliptical infarction. Stresses are affected in both the circular and the elliptical infarctions. The reduction of the circumferential strains at 5 weeks post-MI agrees well with the experimental observations. With further development, we expect that this agent-based homogenized approach could provide useful insights for the clinical practice of MI patient management.
where is derived from the derivative of the second Piola-Krirchoff stress P, For the given strain energy function (33), we have