An updated Lagrangian constrained mixture model of pathological cardiac growth and remodelling

Progressive left ventricular (LV) growth and remodelling (G&R) is often induced by volume and pressure overload, characterized by structural and functional adaptation through myocyte hypertrophy and extracellular matrix remodelling, which are dynamically regulated by biomechanical factors, inﬂammation, neurohormonal pathways, etc. When prolonged, it can eventually lead to irreversible heart failure. In this study, we have developed a new framework for modelling pathological cardiac G&R based on constrained mixture theory using an updated reference conﬁguration, which is triggered by altered biomechanical factors to restore biomechanical homeostasis. Eccentric and concentric growth, and their combination have been explored in a patient-speciﬁc human LV model under volume and pressure overload. Eccentric growth is triggered by overstretching of myoﬁbres due to volume overload, i.e. mitral regurgitation, whilst concentric growth is driven by excessive contractile stress due to pressure overload, i.e. aortic stenosis. Different biological constituent’s adaptations under pathological conditions are integrated together, which are the ground matrix, myoﬁbres and collagen network. We have shown that this constrained mixture-motivated G&R model can capture different phenotypes of maladaptive LV G&R, such as chamber dilation and wall thinning under volume overload, wall thickening under pressure overload, and more complex patterns under both pressure and volume overload. We have further demonstrated how collagen G&R would affect LV structural and functional adaption by providing mechanistic insight on anti-ﬁbrotic interventions. This updated Lagrangian constrained mixture based myocardial G&R model has the potential to understand the turnover processes of myocytes and collagen due to altered local mechanical stimuli in heart diseases, and in providing mechanistic links between biomechanical factors and biological adaption at both the organ and cellular levels. Once calibrated with patient data, it can be used for assessing heart failure risk and designing optimal treatment therapies.


Introduction
Heart disease is known as the leading cause of human mortality worldwide [1] . The incidence of heart failure, such as after myocardial infarction, has remained persistently high due to maladaptive growth and remodelling (G&R) of heart in the activation of compensatory mechanisms, which are aimed at maintaining normal pump function. Cardiac G&R involves physiological, neurohormonal, molecular and cellular adaptations in response to internal or external environmental changes, such as mechanical loading conditions, i.e. volume or pressure overload [2] . These adaptive processes under disease-induced stimuli are usually referred to as pathological cardiac G&R.
Of the various pathological conditions associated with cardiac G&R, one such example is an overload of left ventricular (LV) pressure in systole. This can be caused by aortic valve stenosis or systemic hypertension. In order to accommodate this resultant high after-load, an increase in myocardial contractile force is necessary, with the associated adaptation of myocytes realised through sarcomerogenesis due to the very low rate of myocyte turnover. More specifically, myocytes will react to the increased systolic wall stress by adding sarcomeres in parallel, leading to wall thickening at the organ level, so-called concentric growth [3] . A second common condition is an overload of LV cavity volume in diastole, e.g. due to mitral valve regurgitation, which results in a high myocardial stretch at end-diastole. When this occurs, myocytes react by adding sarcomeres in serial, which leads to LV chamber dilation and wall thinning, so-called eccentric growth. Not only structural and functional adaption in myocytes, the extracellular matrix, mainly the collagen network, also experiences extensive G&R with altered cardiac fibroblast function, which is further dynamically regulated by biomechanical factors [4] . Consequently, systolic stress and diastolic stretch have commonly been considered to be the triggers of pathological cardiac G&R in previously published studies [5][6][7] .
In this study, we will only briefly discuss two main frameworks of modelling cardiac G&R: the kinematic growth theory [8] and the constrained mixture theory [9] . Details on both methods are given in Niestrawska et al. [10] , Sharifi et al. [11] , Yoshida and Holmes [12] , Humphrey [13] , while a summary of recent G&R studies on the cardiovascular system is included in Table 1 . The kinematic growth theory, first proposed by Rodriguez et al. [8] , involves a multiplicative decomposition of the deformation gradient tensor F , such that F = F e F g , following the concept of plasticity. The inelastic growth tensor F g describes changes in shape and size of a stress-free tissue, and the elastic deformation gradient tensor F e assembles the added and original mass into a new compatible grown state, in which residual stress may be generated due to incompatible growth. The evolution of F g , commonly referred to as the "growth law", is necessary to define the relationship between growth extent and stimuli from changed local environments (stretch, stress, inflammation, etc.), while we note that the tissue mechanical stress is singularly derived from F e (i.e. no contribution from F g ). Since its first introduction, kinematic growth theory has been widely used to model cardiac G&R under physiological and pathological conditions. For example, the generic eccentric and concentric growth [5,6,[14][15][16] , reversal G&R [17,18] , explaining experimental observations [19,20] , as well as post myocardial infarction [21] .
Kinematic growth theory is a simple yet elegant approach to describe and predict soft tissue G&R, with additional benefits of simple concept and fast computation. A downside of this simplicity is that the phenomenological growth law usually uses one growth tensor to describe overall G&R in biological tissue, but not individual F g for different constituents. Thus, the fundamental mechanobiological mechanisms of biological constituents are usually missing. To address that, constrained mixture theory was proposed by Humphrey and Rajagopal [9] to take into account the deposition and degradation of individual constituents. The main hypothesis of this fundamentally different framework is that the added mass of a constituent i can be deposited at an intermediate time τ with a certain pre-stretch F i pre (τ ) along with degradation of the existing mass. Known as the classical constrained mixture model, this has primarily been used to study vascular G&R [22][23][24][25][26][27] . However, tracking the reference configuration of each individual mass comes at a very high computational cost. To overcome this, Cyron et al. [28] introduced a temporally homogenized constrained mixture model which does not require tracking the reference configuration of each added mass. Instead, for each constituent i , the total G&R tensor F i gr is assumed to be a product of an inelastic remodelling tensor F i r and an inelastic growth tensor F i g . As such the tensor F can be decomposed into F = F i e F i r F i g where, as before, F i e is the elastic deformation tensor. Different configurations of added mass at different times can thus be homogenised into F i r meaning that the initial reference configuration can instead be used to model G&R, saving considerable computational costs. This method was then adopted by Gebauer et al. [7] who applied it to an idealised LV to study both mechanobiological stability and reversal G&R.
In the temporally homogenised constrained mixture model, the evolution of the inelastic remodelling tensor, F i r , must still be established, which can in itself prove challenging. Furthermore, a single reference geometry is irrelevant to newly added mass [14] and could also be complicated if a topological change is introduced in the middle of G&R process. Hence, in this study, we aim to develop a G&R model using updated reference based on the constrained mixture model by taking into account different mechanical properties and adaptation of individual myocardial constituents, which are the ground matrix, myofibres (bundles of myocytes) and collagen network. Instead of tracking historical configurations of each added mass, we will introduce an updated Lagrangian approach by updating the reference configuration of the homogenized mixture at each G&R cycle. This method of updating the reference configuration has previously been applied to kinematic growth models of cardiac G&R [14,18] , but it has not been studied in constrained mixture models of cardiac G&R. We then test this new framework using a subject-specific LV geometry under pressure and volume overload. Finally, qualitative comparisons with experimental studies and clinical observations will be analyzed, including wall thickening/thinning and residual stress.

Method
The Method section is structured as follow: first, Section 2.1 briefly describes the updated reference configuration approach; secondly, Section 2.2 introduces the constrained mixture theory within one G&R cycle; thirdly, Section 2.3 describes the new G&R framework we propose, including updated reference configuration with residual stress, growth laws and strain energy function evolution; fourthly, Section 2.4 lists various G&R models; and finally, Section 2.5 describes a human left ventricular model for cardiac G&R modelling. A summary of nomenclature is included in Table 2 for all notations used throughout this paper.

Updated reference configuration
As different constituents in soft tissue experience a continuous turnover (deposition and degradation), using a fixed reference configuration is irrelevant for newly added mass and unnecessary for degraded mass. For that reason, an evolving reference configuration is needed in soft tissue G&R, which is naturally incorporated in the constrained mixture theory. Using kinematic growth theory, Kroon et al. [14] was the first group to investigate the difference between fixed and updated reference configurations in an LV G&R model, and they suggested that updated reference configurations would be more physiological instead of fixing it, and further concluded that the reference configuration plays a crucial role when Table 1 Summary of recent G&R studies on the cardiovascular system. I : the identity matrix, ϑ: the growth amount, n : the normal direction. Refer to [10][11][12] for excellent recent reviews.  1 the orthogonal fibre-sheet-normal system at the beginning of one G&R cycle f u − s u − n u the updated orthogonal fibre-sheet-normal system at the end of one G&R cycle f 0 , s 0 , n 0 the fibre, sheet, sheet-normal directions in a fictitious configuration B 1 , 0 ( Fig. 3 ) c s , c n the cross-fibre directions with respect to f 0 in a fictitious configuration B 1 , 0 ( Fig. 3 1. Schematic illustration of the update reference configuration framework in the first two G&R cycles. The configuration B 0 is stress-free for the first G&R cycle, whilst B 1 and B 2 with non-zero residual stress due to incompatible G&R, G 1 and G 2 are pure growth tensors, F 1 e and F 2 e are assembling tensors after the pure growth step that are responsible for non-zero residual stress, F 1 eg and F 2 eg are the mapping deformation tensors between two grown configurations without loading, and F 1 r is the residual stretch tensor. This updated reference configuration is similar to Kroon et al. [14] . modelling G&R in soft tissue. Subsequently, Lee et al. [18] employed the updated reference scheme when simulating an integrated electromechanical-growth heart model. The classical constrained mixture theory (CMT) requires tracking the evolution of reference configurations of individual constituents [9] with time, which has limited its wide applications [10] . In this study, to overcome the difficulty of tracking all historical reference configurations of grown materials, an updated reference framework similar to Kroon et al. [14] is adopted as shown in Fig. 1 . In this new framework, the initial configuration B 0 is assumed to be stress-free ( σ r = 0 ) before growing to a new compatible configuration B 1 with nonzero residual stress ( σ r = 0 ). Following the kinematic growth framework, this process of G&R from B 0 to B 1 can be split into two sub processes. In the first subprocess, B 0 grows purely with a growth tensor G 1 towards a fictitious grown configuration B 0,gr , an incompatible configuration with zero residual stress by assuming pure growth will not generate stress. In the second sub-process, an elastic deformation ( F 1 e ) is followed to assemble incompatible grown materials in B 0,gr into a compatible grown configuration B 1 , in which a non-zero residual stress field could be induced by F 1 e . The tensor F 1 eg describes the total deformation in one growth step that contains both the pure growth and the elastic deformation. In the next growth step, B 1 is used as the new reference configuration, in other words, B 1 is the updated reference configuration for the next growth cycle. By pulling back B 1 with its residual stretch tensor ( F 1 r ), a corresponding fictitious incompatible configuration B 1 , 0 with zero residual Fig. 2. Schematic illustration of G&R and loading process for three main constituents of the myocardium (the ground matrix, myofibres and collagen fibres) in one G&R cycle, consisting of the long-time G&R step and the short-time loading step. The configuration B 1 is the updated reference, which will grow into an incompatible configuration B 1,g under the pure growth tensor G i , followed by a remodelling process towards B 1,gr with an inelastic remodelling tensor P i r . Since B 1,gr is still incompatible, finally an elastic tensor F i r assembles all grown constituents to a compatible state B 2 . That is the G&R step. The loading step begins at B 2 that deforms into B ζ under external loads, the corresponding deformation tensor is F i ζ .
stress can be found. We then assume all material growth occurs in this stress-free incompatible configuration with a pure growth tensor G 2 , leading to, as before, an incompatible grown configuration B 1,gr , then followed by an elastic deformation F 2 e to reach the final compatible grown configuration B 2 with none-zero residual stress. Note that F r is the corresponding residual deformation related to σ r , and that the two configurations B 0,gr and B 1,0 are virtually identical. All subsequent growth steps follow the same way from B 1 to B 2 by using the previous grown configuration as the updated reference configuration. Furthermore, we note that the assumption of zero σ r in the very first configuration B 0 can be relaxed if the initial residual stress is known.

Constrained mixture theory
Following the classical constrained mixture theory, we consider the myocardium to be composed of three main constituents, which are the ground matrix, myofibres (myocyte aggregates) and collagen network, denoted by { g,m,c } , respectively. In constrained mixture theory, although the deformation of each component is quantified in its own reference configuration, all components share the same geometrical configuration before and after a G&R process, i.e. sharing the same total deformation mapping tensor F [9,28] . The constrained G&R process proposed in this study is shown in Fig. 2 and consists of two time-scale steps: Fig. 3. G&R in the updated reference configuration with residual stress. The configuration B 1 is now with non-zero residual stress due to the residual stretch F i r . Fictitious incompatible stress free configurations ˜ B 1 , 0 , B 1 , 0 , B 1,g and B 1,gr are included to depict the G&R process. In particular, ˜ B 1 , 0 is pulled back from B 1 by the residual stretch F i r . After polar decomposition of F i r , an equivalent stress-free B 1 , 0 is obtained from ˜ B 1 , 0 through a rotation process ( R i r ). Note by applying the corresponding left stretch tensor V i r , we can also arrive B 1 . In this study, we assume the pure growth starts from B 1 , 0 to B 1,g and followed by an inelastic remodelling to B 1,gr . Finally the incompatible fictitious state B 1,gr is assembled into B 2 by an elastic tensor F i r , or deforms into B ζ under external loadings with the total elastic deformation F i e . The index i represents the ground matrix, myofibres and collagen fibres.
• the G&R step ( B 1 → B 2 ), the long-time scale occurring within days and weeks; • the loading step ( B 2 → B ζ ), the short-time scale occurring within seconds, i.e. every heartbeat.
In detail, a G&R step starts from a compatible configuration B 1 , each individual constituent then grows into a fictitious incompatible grown configuration B 1,g with individual growth tensors ( G i ). A further remodelling process, denoted by P i r , follows which involves e.g. rearranging collagen fibres and results in an incompatible configuration B 1,gr . In this study, the growth tensor is denoted by G i and the remodelling tensor by P i r and we note that neither of these will introduce stress, meaning that no stress exists in B 1,g and B 1,gr . By assuming there exist elastic tensors ( F i r ) for each constituent based on constrained mixture theory, a compatible configuration B 2 can be obtained where the residual stress may exist depending on G i and P i r . Note B 2 is the unloaded state that could be measured experimentally. The above is the long-time G&R step. For the short-time loading step, the unloaded configuration B 2 is constantly loaded with external forces under significant deformation B ζ , i.e. LV blood pressure. Thus B 2 → B ζ represents a purely mechanical process without G&R due to the very short time. Finally, the total deformation from B 1 → B ζ can be written as where the total elastic deformation is F i e = F i ζ F i r . We define the inelastic part F i gr = P i r G i that doesn't directly lead to stress generation, and F eg = F i r P i r G i for the mapping from B 1 to B 2 .

Update reference configuration with residual stress
As illustrated in Fig. 1 , except for the first G&R cycle, subsequent G&R commences from an updated reference configuration with non-zero residual stress in general. We now consider the reference configuration ( B 1 ) in Fig. 2 with non-zero residual stress. In order to include residual stress in B 1 , we rely on the residual elastic deformation tensor F r . Figure 3 summarises how a G&R cycle is modelled using an updated reference configuration ( B 1 ). Firstly, B 1 is pulled back to a fictitious incompatible state ˜ B 1 , 0 using the residual elastic deformation tensor F i r of each constituent. Secondly, from the polar decomposition, we have F i r = V i r R i r , which means we can apply a rotation tensor R i r to ˜ B 1 , 0 to obtain a zerostress incompatible configuration B 1 , 0 . Finally, the left stretch tensor V i r provides a direct mapping from B 1 , 0 to B 1 . Thus the two fictitious stress-free configurations ˜ B 1 , 0 and B 1 , 0 can be obtained by using the residual deformation tensor.
Existing G&R modelling research usually assumes that pure growth occurs in a stress-free state [5] . We also follow this by considering pure growth will only occur in B 1 , 0 but not in B 1 . In each G&R cycle, the long-time growth step is B 1 , 0 → B 1 → B 2 , and the short-time loading step is B 2 → B ζ . Under this assumption of restricting pure growth to B 1 , 0 , an equivalent fictitious pathway for both the short-time and long-time steps can be established as shown in Fig. 3 , specifically • the long-time G&R step: • the short-time loading step: B 1,gr → B ζ .
From B 1 , 0 to B ζ , the mapping for each individual constituent can be written as in which G i is the pure growth tensor for i th constituent, P i r is the inelastic remodelling tensor, and F i e is the total elastic tensor consisting of the residual deformation ( F i r ) and the loading ( F i ζ ). Therefore, which is required for computing Cauchy stress at B ζ . Similarly, the residual deformation from B 1,gr to B 2 is given by The algorithm for obtaining F i r can be found in Appendix E ( Algorithm 2 ). Here we note that the geometries in B 1 , B 2 and B ζ can be observable.
In the next G&R cycle, B 2 becomes the updated reference configuration as highlighted by the red rectangle in Fig. 3 . By pulling it back using F i r , we can obtain ˜ B 2 , 0 where next pure growth occurs.

2.3.2.
The fibre-sheet-normal system after G&R As can be seen from Fig. 3 , B 1 is the reference configuration for the current G&R cycle, and B 2 is the updated reference configuration for the next G&R cycle. In this study, we assume myofibres always form a layered architecture with three orthogonal directions in the (updated) reference configuration, the so-called fibre-sheetnormal system. The assumption of a layered fibre-sheet-normal system in (updated) reference configurations is based on various experimental studies [29][30][31] and has been widely used in the cardiac modelling community [32,33] . Further discussion can be found in Section 4.6 . In B 1 , the unit fibre direction is f 1 , the unit sheet direction is s 1 , and the unit sheet-normal is n 1 . After a G&R process with F eg , we have It is worth mentioning that f 1 , s 1 and n 1 are not orthogonal in general. When B 2 is used as the updated reference configuration for the next G&R cycle, we implicitly assume a further remodelling process to ensure an orthogonal fibre-sheet-normal system ( f u − s u − n u ), such that In an updated reference configuration, we further assume the ground matrix is isotropic, and the collagen fibres follow the myocytes with the same fibre-sheet-normal system for the local material coordinate system.

Growth tensor G i
We now consider the directional growth that is mainly determined by myocytes since they make up 70% of the myocardial mass [34] . The orthogonal fibre-sheet-normal system f 1 − s 1 − n 1 at B 1 for myocytes can be pulled back to B 1 , 0 using the residual stretch tensor V i r , Instead of defining G i using f i 0 − s i 0 − n i 0 because it is not an orthogonal system in general, we assume the primary directions of growth are along the myofibre direction f i 0 and the cross-myofibre directions which are defined as in which c i s is along the radial (transmural) direction in general. Finally, the growth tensor related to B 1 , 0 is defined as where ϑ f ,i and ϑ s ,i are the incremental growth ratios along the myofibre and cross-myofibre directions, respectively. Following [5] , the evolution equation for stretch-driven fibre growth ( θ f ) with respect to the initial non-grown state (i.e. B 0 in in which τ f is the parameter related to the duration of one G&R cycle with a unit of time , θ max,f is the maximum growth ratio along the fibre direction, λ crit is the critical threshold of fibre stretch, F e is the total elastic stretch tensor, and H(·) is the Heaviside function. The total growth ratio from the first growth cycle (the initial non-grown state) to the current (n + 1) th growth cycle is θ f n +1 , and the incremental growth ratio from the n th growth cycle to the (n + 1) th growth cycle is In a similar way, the stress-driven cross-fibre growth ratio related to the initial non-grown state [5] in which τ s is also the parameter related to one G&R cycle with a unit of time , θ max,s is the maximum growth ratio along the sheet direction, p crit is the critical threshold value, and σ can be the peak stress in systole. Accordingly, the incremental growth ratio is Equations (10) and (12) are discretized with finite difference for the first-order material time derivative. At each G&R cycle, we allow θ f and θ s to achieve the most possible growth. The exact unit for τ f and τ s needs to be determined with experimental data, which could be day or week , etc. Details of how these two equations are solved numerically can be found in Appendix E ( Algorithm 1 ).

Strain energy functions
For the stress-free configuration B 1 , 0 in Fig. 3 , the strain energy function of the myocardium is where φ g n , φ m n and φ c n are volume fractions of the ground matrix, myofibres and collagen fibres in the n th growth cycle, respectively. Here, we also assume each constituent is incompressible with the same constant density. The strain energy functions for the ground matrix ( g ) and myofibre ( m ) are defined as where a g , b g , a m , b m are material parameters, e , and f m 0 is the mean fibre direction in the reference con- only the stretched fibres can bear load passively.
Following the strain energy function proposed by Holzapfel and Ogden (HO) [35] , we consider the collagen fibre network is responsible for the orthogonal material response in the myocardium, thus its strain energy function is where a cf , b cf , a cs , b cs , a cn , b cn are material parameters, and where C c = (F c e ) T F c e , f c 0 , s c 0 and n c 0 are the mean fibre, sheet and sheet-normal directions of the collagen network in the reference configuration B 1 , 0 . We further assume the collagen network has transverse isotropic mechanical property, i.e., a cs = a cn and b cs = b cn .
After each growth cycle, the total myocardial volume and individual volumetric fraction need to be updated. The incremental kinematic growth ratio at the end of the n th growth cycle is where V n −1 is the total volume, V i n −1 is the individual volume of constituent i in the (n − 1) th growth cycle, and J i n = det (G i ) . The updated volume fraction φ i n of each constituent at the end of the n th growth cycle is given by As introduced in Section 2.2 , each growth cycle consists of a long-term G&R step and a short-term loading step. In the shortterm loading step, the growth tensor G i and the remodelling tensor P i r are assumed to be constant, then the passive stress response of the myocardium is only a function of the elastic deformation tensor F i e as suggested in Luo and Zhuan [6] , that is where . Following the additive stress approach, the total Cauchy stress at the n th growth cycle is with the active stress tensor in which ˆ f m is the unit myofibre direction at the current configuration, and T a is the active tension generated by myocytes. The active tension is determined by a well-established length-tension model [36,37] , in which T max is the maximum isometric active tension, l 0 is the minimum sarcomere length to produce active stress, l r is the stress-free sarcomere length, ω is a time function after onset of contraction, and B is a constant. Details of this active contraction model can be found in Guccione and McCulloch [36] , Sack et al. [37] .

Eccentric, concentric and hybrid G&R
In this study, we will mainly focus on generic pure eccentric growth driven by stretch, pure concentric driven by stress, and the combination of the two types of growth (hybrid). The actual myocardial G&R will be much more complex than what we have considered here.

Myofibres
In eccentric growth, the myofibre growth tensor G m is defined as where ϑ f,m is the incremental growth ratio along myofibres, i.e. sarcomerogenesis in series. In contrast, for concentric growth, the myofibre growth tensor is where ϑ s,m is the incremental growth ratio for the cross-myofibre direction, i.e. sarcomerogenesis in parallel. Accordingly, for hybrid growth, we define Note inelastic remodelling is not considered here, i.e. P m r = I , thus, the total elastic tensor for myofibre from B 1,gr to B ζ shown in Fig. 3

Ground matrix
We assume that the ground matrix does not grow ( G g = I ) but only experiences an inelastic remodelling process with conserved volume to be compatible with myofibre G&R. This suggests that the same directional stretch for the ground matrix as myofibres grow, while the directions not experiencing myofibre G&R will be compressed due to the incompressible constraint. Note myofibre G&R refers to sarcomerogenesis either in parallel or in series. Based on such an assumption, the inelastic remodelling tensor for the ground matrix is The total elastic tensor for the ground matrix from B 1,gr to B ζ shown in Fig. 3

Collagen fibres
The collagen network, the main component of the extra-cellular matrix (ECM) of the myocardium, is tightly regulated by the deposition of new collagen from myofibroblasts and degradation of old collagen by matrix metalloproteinases (MMPs) [38] . Under pathological G&R, myocardial ECM is mediated by the activation of myofibroblasts and MMP expression accompanied by myocyte G&R through sarcomerogenesis, together they will contribute to the overall ventricular G&R. Studies have shown that increased MMP activity will lead to partial degradation of the interstitial ECM that precedes a chronic volume overload or rapid pacing [39,40] . It also has been demonstrated that increased MMP activity could lead to progressive LV G&R in diseased hearts through degradation or damage of the collagen network, resulting in myocyte slippage and myocardial wall thinning [41] . Activation of myofibroblasts is another response under pathological conditions. Ultimately, it will lead to excessive collagen deposition and ECM expansion [42] . We point the reader to Frangogiannis [42] for a recent review on ECM G&R in heart. In this study, modelling of specific biologically-informed collagen G&R patterns is out of scope, thus we only consider three simplified scenarios for collagen G&R aiming to predict general cardiac G&R patterns. In specific, 1. CF-R : Corresponds to an extreme condition when MMP only degrades collagen crosslink to allow myofibre slippage, but the total mass of collagen network is preserved. Thus the collagen network only experiences inelastic remodelling to be compatible with myofibre G&R. We have G c = I , and the inelastic remodelling tensor for the collagen network Then the total elastic tensor for the collagen network from B 1,gr to B ζ shown in Fig. 3 is F c e = FV c r (P c r ) −1 .
2. CF-E : We assume the collagen network does not experience any growth and remodelling, which would correspond to an extreme condition when both MMP and myofibroblasts are not activated. Then G c = I and P c r = I , and the total elastic tensor for the collagen network from B 1,gr to B ζ shown in Fig. 3 is F c e = FV c r . 3. CF-G : The collagen network is assumed to grow exactly the same as the myofibres with P c r = I , corresponding to an ideal scenario when the net G&R of the collagen network is the same as the myofibre G&R under well controlled myofibroblast activation and MMP expression. Then and the total elastic tensor of the collagen network from B 1,gr to B ζ shown in Fig. 3 is F c e = FV c r (G c ) −1 .

A human LV G&R model
A healthy human LV geometry reconstructed at early-diastole adapted from [33] is used here. Figure 4 (a) shows this LV geometry with a rule-based fibre structure. Myofibres linearly rotate from epicardium ( −60 • ) to endocardium ( 60 • ) [43] . A simplified systemic circulation system [44] is attached to the LV model to provide physiologically accurate boundary conditions (see Fig. 4 (b)). We consider the LV geometry at early diastole to be the unloaded reference configuration though in reality it is not because of constant pressure loading inside LV. See further discussion about un- loaded configuration in Section 4.4 . Details of this LV model can also be found in Gao et al. [33] , Guan et al. [44,45] .
The hypothesis that there exists a homeostatic state, an internal steady state maintained by living tissue, has become the cornerstone in biological soft tissue G&R. Biomechanically, it refers to a preferred mechanical state targeted by various mechanobiological activities [4] , i.e. myocyte hypertrophy and fibrosis. It has been demonstrated that mechanical stimuli such as stress/strain play a significant role in soft tissue G&R though many open questions are yet to be answered. As discussed in Eichinger et al. [4] , homeostasis shall be in a range near the preferred state but not a fixed value, and very adaptive under different loading and pathological conditions. In this study, average values of end-diastolic stretch and peak active tension have been used to represent the homeostatic states, similar to other modelling studies [18,46] . Specifically, we have chosen the average tr (σ a ) at peak systole as the only stimulus to trigger concentric growth, neither the Mandel stress [5] nor the trace of total Cauchy stress [6] . It is worthwhile to mention that future work shall consider the total stress tensor over cardiac cycles.
The LV model is firstly simulated without G&R in order to determine λ crit in Eq. (10) and p crit in Eq. (12) , denoted as the baseline case. In detail, the LV is firstly preloaded with a linearly increased pressure up to normal end-diastolic pressure (EDP), i.e. 8 mmHg, and λ crit is calculated as the spatial mean myofibre stretch at enddiastole. Then, the iso-volumetric contraction begins, followed by the systolic ejection when the LV pressure exceeds a prescribed aortic pressure (80 mmHg), and finally, the ejection ends when the LV pressure is lower than the aortic pressure. Values of λ crit and p crit are obtained from the baseline case as where N is the total number of elements in the LV model, λ max i is the maximum myofibre stretch in diastole for the i th element, and p max i is the maximum of the trace of active stress tr (σ a ) in systole. Parameters of the baseline model can be found in Table 3 .
The details for determining the passive material parameters can be found in Appendix A by ensuring that the baseline model behaves physiologically [47] . Young et al. [48] reported that the LV EDP was twice or triple higher in hearts with eccentric G&R under mitral regurgitation than the value from normal hearts. For that reason, increased EDP is usually applied when modelling eccentric cardiac growth [17,46] . To trigger eccentric growth, the LV EDP is increased to 16 mmHg with no change in contractility. Aortic stenosis has been considered the paradigm for pressure overload in LV, triggering a hypertrophic response with increased wall thickness [49] . To trigger concentric growth, only the aortic resistance C AV V is increased from 1.0 to 10 MPa mm 2 s/tonne while the contractility ( T max in Eq. (23) ) remains constant again, see Fig. 4

(b). Ten-folder increase in C AV
V is chosen such that the end-systolic pressure at the final growth cycle is within the reported range in ( 197 ± 5 mmHg) [50] , see Fig. 9 (ac). In hybrid growth, both EDP and the aortic resistance are increased to the values of eccentric and concentric growth. Table 4 lists the parameter values for the eccentric growth and concentric growth laws, respectively.
The nonlinear finite-element software ABAQUS EXPLICIT (Dassault Systemes, Johnston RI, USA) is used to simulate LV G&R, and the strain energy function defined in Eq. (14) is implemented using a Fortran user-subroutine. In-house developed MATLAB codes are used to update the reference configuration of the grown LV geometry. All simulations are run on a Windows workstation with 42 CPU cores (Intel(R) Xeon(R) CPU E5-2680 v3 @ 2.50 GH, 64 GB), and one growth cycle takes around 3.5 h.
Please note in our simulations, cycle 0 is the LV model in the baseline state, and the pathological boundary conditions are applied from cycle 1, thus there is no G&R in cycle 0, and the initial growth ratio for cycle 1 is one. Subsequently, the growth tensor for each constituent is determined from the previous cycle. We found most simulations finish in 10 cycles when the average incremental growth to the previous cycle is within 1%. We further compared the mean myofibre stretch at end-diastole and peak-systolic active stress from the 10th cycle to the thresholds λ crit , p crit , and both values are close to the predefined thresholds, see Fig. D.1 in Appendix D for a detailed comparison. We note here that using the average incremental growth as the criteria may not be the optimal one, and some cases after 10 growth cycles may not meet the chosen criteria. Considering the computational cost of running all simulations and to have a like-for-like comparison, we ran all simulations with 10 growth cycles in the following analysis. Additional growth cycles can be simulated, see Section 4.5 for further discussion. Figure 5 (a) shows the final grown LV geometries under three G&R scenarios of collagen. In general, the endocardium experiences much higher θ f than the epicardium (see Fig. 5 (a)), and the most growth occurs at the base and mid-ventricle but least at the apex.

Eccentric growth
Large spatial variations of θ f can be seen in the three cases, which are resulted from spatially varied mechanical environments. For example, the apex usually experiences the least stretch and stress, thus with the least growth. The final grown LV geometries are similar for cases CF-R and CF-G, whose cavity volumes are much larger than the baseline case. The least growth is found in case CF-E, which could be potentially explained by the elastic stretch in collagen fibres induced by myofibre G&R, and this in turn hinders further myofibre G&R. The average growth ratio θ f is shown in Fig. 5 (b). It can be found that θ f gradually increases until it reaches a plateau when no growth can occur either because of approaching the growth limit θ max,f or below the critical threshold of fibre stretch λ crit . The final θ f at the end of the 10th cycle is 1.28 for case CF-R, 1.26 for case CF-G, and 1.12 for case CF-E. The average incremental growth ratio θ f decreases with time as shown in Fig. 5 (c), and θ f approaches 1 much faster in case CF-E compared to other two cases. Figure 5 (c) further shows the standard devia- Table 3 Parameter values for the lumped circulation model, and material parameters for passive and active responses for the baseline case.
Parameters of the lumped circulation model [45] C AV   tion of ϑ f , from which we can see that its value is getting much smaller in late growth cycles. The much reduced standard deviation of ϑ f also suggests the LV reaches a stable grown state at the 10th cycle since θ f does not change much temporally with further growth cycles. Figure 6 (a-c) shows the pressure-volume (P-V) loops of the three cases within 10 growth cycles. In general, the P-V loops gradually shift to the right with reduced systolic pressure but much increased end-diastolic volume (EDV) and end-systolic volume (ESV). Case CF-E has the least LV cavity compared to other scenarios. Figure 6 (d) shows ejection fraction (EF), and Fig. 6 (e) shows the corresponding stroke volume (SV). Both EF and SV decrease significantly after the first few G&R cycles, except for case CF-R at the second cycle. That is presumably because inelastic remodelled collagen fibres in case CF-R allow the LV to be filled more, leading to increased stroke volume according to the Frank-Starling law. As the wall gets thinner and thinner, from the third cycle, the overall contractile force from the thin LV wall would be reduced, thus leading to decreased SV. At the final G&R cycle, EF is 31.4% for case CF-R, 29.7% for case CF-G, and much higher in case CF-E (40.7%). Although case CF-E has the least SV (60 mL), but still within the lower bound of reported values in healthy humans (95 ±14 mL) [51] . The average wall thickness at the unloaded state is shown in Fig. 6 (f). LV wall thinning occurs in cases CF-R and CF-G, but not in case CF-E that surprisingly has a slightly thicker wall compared to the baseline case with 2% increase. Detailed comparison of LV pump function can be found in Table 5 .
Because of no growth in the ground matrix, its volume fraction is gradually reduced as shown in Fig. 7 (a). The volume fraction of myocytes increases for all three scenarios as shown in Fig. 7 (b), while its increase is much reduced when the collagen does not experience G&R. The highest myofibre growth occurs when the collagen only experiences inelastic remodelling. Figure 7 (c) shows the change in the collagen volume fraction. In case CF-G, a slightly increased collagen volume fraction is observed since it grows exactly the same as myocytes, but for the other two scenarios, its volume fraction decreases, in particular for case CF-R because of no growth in collagen. Figure 7 (e) shows the uni-axial stretch-stress response along the myofibre direction at the end of the 10th growth cycle. Compared to the baseline case, myofibre stiffness is slightly higher for case CF-G because of growth in both myofibres and collagen, but slightly lower for the other two cases because of the reduced volume fraction of collagen fibres. Detailed comparison of volume fraction can be found in Table 6 .    In summary, eccentric growth would lead to increased enddiastolic and end-systolic cavity volume with decreased ejection fraction and stroke volume. Collagen G&R could play an important role in overall LV wall G&R. Specifically, by eliminating collagen G&R, the grown LV has the least increase in its cavity volume without wall thinning and nearly the same passive stiffness along the myofibre direction.

Concentric growth
Similar to the cases of eccentric growth, myocardial G&R mainly occurs at the endocardium and the base as shown in Fig. 8 (a), and nearly no growth in the apex. LV cavities in three scenarios are all reduced with a thicker ventricular wall compared to the baseline case. Interestingly, case CF-E grows longitudinally with much higher θ s near the base compared to the other two cases. The average θ s at each G&R cycle is plotted in Fig. 8 (b) where case CF-E has a maximum value of 1.44 with reduced values in cases CF-G (max: 1.29) and CF-R (max: 1.27). Figure 8 (c) shows the average incremental growth ratios ( θ s ). It can be seen that cases CF-G and CF-R have much slower incremental growth compared to case CF-E. Again, the standard deviation of θ s is getting smaller and smaller at late growth cycles.
The P-V loops for the three cases are very similar as shown in Fig. 9 (a-c), all with decreased ESV and EDV, and significantly increased end-systolic pressure (ESP). For example, ESP is 182.1 mmHg in case CF-E at the 10th cycle, 155.3 mmHg for CF-R and 161.9 mmHg for CF-G. EDVs are reduced to 73.6 ml for CF-E, 58.6 ml for CF-R, and 62.0 ml for case CF-G. Similarly, ESVs are 30.4 ml for CF-E, and around 23 ml for the other two cases. Figure 9 (d) shows EFs at each growth cycle. EF at the first G&R cycle (43.4%) is much lower compared to the baseline case (58.3%). Because of concentric growth, the thickened wall would generate higher active tension to overcome the increased resistance due to aortic stenosis, EF gradually increases to the normal range ( > 50%) after the first 2 growth cycles. Interestingly, EF peaks at the 5th growth cycle and then decreases for all three cases, which may indicate that at the 5th growth cycle, the heart is at its optimal compensatory condition, further G&R might lead to a decompensated state, i.e. decreased EF and SV. At the final growth cycle, EFs for all three scenarios are all above 55%, while case CF-E has the least EF. Similar trends can be found for SVs as shown in Fig. 9 (e). Specifically, SV peaks after the first several growth cycles and then decreases significantly till the 10th cycle due to the restricted LV cavity, which is much smaller compared to the baseline case (59.6 mL). Only case CF-E reaches a higher SV than the baseline case from the 3rd to the 5th growth cycles but with the highest ESP. Figure 9 (f) shows the average wall thickness, which is the least in CF-E, followed by CF-G, and the highest in CF-R. Detailed comparison of LV pump function can also be found in Table 5 .
The volume fractions of the ground matrix decrease for all three scenarios ( Fig. 10 (a)), and the myofibre volume fraction increases as shown in Fig. 10 (b), while the volume fraction of collagen only increases in case CF-G ( Fig. 10 (c)) as expected. Figure 10 (d) further shows the uni-axial stretch-stress responses of the myocardium at the end of the 10th cycle. Again, case CF-G is slightly stiffer compared to other cases because both myofibres and collagen grow simultaneously. Detailed comparison of volume fraction can be found in Table 6 .
In summary, concentric growth leads to decreased end-diastolic and end-systolic cavity volume with increased ejection fraction and end-systolic pressure. The stroke volume initially increases but then decreases significantly at later G&R cycles, which may suggest an intervention could be necessary to prevent adverse G&R when LV achieves the highest EF and SV, for example at the 4th G&R cycle ( Fig. 9 (e)). Compared to cases CF-R and CF-G, case CF-E promotes the most growth in myocytes with the least collagen fibre volume fraction and slightly softer passive property. This leads to a larger stroke volume with the least wall thickening, although ESP is the highest. Again our results might suggest that collagen G&R could also play an important role in concentric G&R. Specifically, by eliminating collagen G&R, the grown LV has the least decrease in its cavity volume with moderate wall thickening and nearly the same passive stiffness along the myofibre direction compared to the baseline case.

Hybrid growth
As shown in Fig. 11 (a), case CF-R leads to a much shorter longitudinal length compared to the baseline case, but not in case CF-E, which experiences both longitudinal lengthening and wall thickening. The much shortened LV geometry in case CF-R may be explained by compaction in the longitudinal direction caused by elongation along the fibre and cross-fibre directions in the collagen network since both θ f and θ s are greater than 1, shown in Fig. 11 (b, c). On the contrary, the additional elastic stretch in the collagen network along the fibre and cross-fibre directions in case CF-E leads to elongated LV because the precompressed collagen network in the longitudinal direction facilitates longitudinal lengthening to accommodate growth along f m 0 and c m s . In general, myocardial growth mainly occurs at the base and the mid-ventricle with the least growth in the apex. As shown in Figs. 11 (b, c), case CF-E experiences the highest growth both along myofibre and cross-fibre directions, with the least crossfibre growth in CF-R that is different when compared to eccentric growth (the highest for CF-R and the least for CF-E). Further simulations were carried out with reduced growth in the myofibre direction, details can be found in Appendix B ( Fig. B.1 ). We found by reducing myofibre directional growth, case CF-E is more concentric growth dominated with limited chamber dilation and much increased systolic pressure, similar to the other two cases. Therefore, it seems stretch-driven growth competes with stress-driven growth in hybrid growth, and such competing effect would lead to different phenotypes of LV G&R compared to pure eccentric or concentric growth. Figure 12 (a-c) show P-V loops under hybrid growth. It can be found that EDV initially increases due to increased EDP, and further ESP also increases to a higher level with a larger ESV than the baseline because of aortic stenosis. Gradually, EDV decreases in both cases CF-R and CF-G, but not in case CF-E. ESP keeps increasing in case CF-G but decreases in case CF-R. Interestingly, EDV and ESV keep increasing in case CF-E while ESP increases very slowly compared to the other two cases. At the 10th cycle, case CF-E has EDV: 193.4 mL and ESV: 104.7 mL and ESP: 151.2 mmHg. EFs of Fig. 9. Evolution of the cardiac function with concentric growth. P-V loops for cases CF-R (a), CF-E (b) and CF-G (c). The black solid line is the baseline case with the aortic valve resistance C V = 1 . 0 , the blue dash line is from the first loading step after increasing aortic stenosis to C V = 10 . 0 , the red dot line is from the 10th G&R cycle, and the green lines are from in-between G&R cycles. Corresponding ejection fraction in (d) with stroke volume in (e) and average wall thickness in (f). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) three cases can be found in Fig. 12 (d), EF for case CF-R peaks at the 4th cycle (58.8%) with an SV of 70.6 mL ( Fig. 12 (e)), though case CF-E has always increased SV, due to increased EDV, its EF is below 45%. Case CF-G has increased EF towards 58.1% at the 10th cycle, its SV increases initially until the 4th cycle, then decreases slightly. Wall thickness keeps increasing for all three cases as shown in Fig. 12 (f). Detailed comparison of LV pump function can be found in Table 5 . Our results seem to suggest that case CF-R experiences the most concentric growth, followed by case CF-G, while case CF-E experiences the most eccentric growth ( Fig. 12 (b)).
Volume fractions of the ground matrix, myofibres and collagen fibres are similar to those of concentric growth as shown in Fig. 13 (a-c). In particular, case CF-E has the highest myofibre volume with the least collagen and ground matrix, which may explain its softer stiffness along myofibre direction due to its least collagen content, see Fig. 13 (d). In summary, the competition between stretch-driven and stress-driven, and the interactions between myofibres and collagen G&R would play an important role in this hybrid G&R, which may lead to either eccentric or concentric-dominated growth. Eliminating collagen G&R may not be the best option as found in pure eccentric or concentric growth. Instead, when collage grows simultaneously with myocytes, the final grown LV geometry and its pump function could be closer to the baseline case as shown in Fig. 12 (c). Therefore, careful control of G&R in myofibres and collagen networks would be critical to prevent adverse remodelling in hybrid growth.

Residual stress and myofibre evolution
At the end of each growth cycle, the LV model is unloaded to a compatible state where residual stress usually exists due to incompatible growth. Figure 14 shows the average circumferential residual stress ( τ cc ) from endocardium (0%) to epicardium (100%) at the 10th cycle. For eccentric growth ( Fig. 14 (a)), τ cc of cases CF-R and CF-G is positive at the endocardium (0.75 kPa and 0.5 kPa), then becomes negative ( −0.37 kPa and −0.24 kPa) at the middle wall, and finally increase to 0.1 kPa at the epicardium. While for case CF-E, τ cc remains very low from the endocardium to epicardium. For concentric growth, all three scenarios experience negative τ cc at the endocardium and positive τ cc at the epicardium, see Fig. 14 (b). Similar trends of τ cc can be found in hybrid growth ( Fig. 14 (c)) Table 7 Myofibre angles at endocardium and epicardium for cases CF-R, CF-E and CF-G at the beginning of the first cycle and the end of the 10th cycle.  [6,[52][53][54] . Negative τ cc at the middle wall and positive at the endocardium and epicardium in eccentric growth ( Fig. 14 (a)) could be explained by higher growth in the middle-wall (stretch-driven), see Fig. 5 . Similar results have been reported in Genet et al. [54] , Choi et al. [55] . We further summarize myofibre angles at the endocardium and epicardium at the end of the 10th cycle in Table 7 . The largest variation can be found in cases CF-R with more circumferentially aligned myofibres, i.e. smallest angle values at the endocardium and epicardium, in particular in hybrid growth, which may also partially explain the most squashed LV shape when unloaded for case CF-R, see Fig. 11 (a).

Comparisons with existing studies
In this study, we have studied generic pure concentric growth, pure eccentric growth and the combination of the two types of growth using this updated Lagrangian constrained mixture theory based G&R model. Parameters in growth laws are not estimated from experimental data but adapted from existing studies [5] . Thus, the eccentric, concentric and hybrid growth shall not be considered as the actual myocardial G&R under pathological conditions. Nevertheless, our new G&R model with updated reference is able to capture concentric growth with sacrcomegenesis in parallel and eccentric growth with sacrcomegenesis in series. Specifically, for pure eccentric growth, the LV at the final G&R state has a much larger LV cavity and lower EF compared to the baseline state, see Fig. 6 . There are a few experimental-based studies on cardiac eccentric growth. For example, Tsamis et al. [56] measured infarct-induced cardiac dilation in vivo in eleven adult sheep up to 8 weeks, and they reported that infarct could lead to longitudinal growth of more than 10% and wall thinning of more than 25%. Costabal et al. [57] studied myocardial G&R under volume overload due to moderate to severe mitral regurgitation in six pigs. They measured myocardial G&R both at the cellular and organ levels within eight-week time. They found that the average myocyte length increased by 26%, corresponding to θ f = 1 . 26 , while the average myocyte width remained almost constant. Both the enddiastolic and end-systolic volumes increased by 47% and 70%, respectively, but with decreased ejection fraction from 55 . 6 ± 4 . 9 % to 47.7 ± 9.5%. More recently, Fan et al. [20] inversely estimated myofibre and transverse growth in 4 three-month old swine who underwent treadmill exercise training for 12 weeks. Using the kinematic growth theory, they inferred myofibre growth θ f = 1 . 58 ± 0 . 196 with a negligible cross-fibre growth θ s = 1 . 04 ± 0 . 08 at week 6, and θ f = 1 . 793 ± 0 . 156 with θ s = 1 . 078 ± 0 . 094 at week 12. The higher θ f and θ s in Fan's study [20] could be potentially explained by the fact that young swine also grow in size physiologically. The average θ f from our pure eccentric growth is 1.28 for case CF-R, 1.26 for CF-G and 1.12 for CF-E, which are in the range reported in those experimental studies mentioned above.
In the clinic, few studies have tried to measure myocardial G&R in vivo but mainly focused on global changes in LV structure and function, such as LV cavity volume and wall thickness. With the development of cardiac magnetic resonance imaging, quantifying local myocardial change is made possible by tracking G&R post-MI. For example, O'Regan et al. [58] studied LV G&R in MI patients within one-year follow-up post-MI in 47 patients. Recently, Li et al. [59] estimated apparent growth tensor using cardiac magnetic resonance images from 16 patients within seven months post-MI. Growth multipliers were summarized in terms of dilation, no-change and shrinkage groups based on LV end-diastolic volume [60] . They reported that myofibre lengthening can be found  In a series of studies, we have shown myofibre and collagen microstructures can play an important role in both passive filling and active contraction [61,62] . Studies have shown that fibre structures are also under continuous change both physiologically and patho-logically. For example, Fomovsky et al. [63] found collagen fibres were more circumstantially oriented in post-MI heart, and similar results were reported in Mojsejenko et al. [64] . Avazmohammadi et al. [16] observed the reorientation of myofibres towards longitudinal direction under pulmonary hypertension in rat hearts, and they further modelled myofibre reorientation by assuming myofibres would adapt their directions according to the largest stretch. Although myofibre orientation in this study is not remodelled according to its stretch state, we still observe a significant myofibre reorientation for all cases, especially for case CF-R with hy-brid growth as shown in Table 7 . The transmural rotation angle of case CF-R at the 10th cycle is 76 . 7 • , a reduction of 43 . 3 • from the baseline case ( 120 • ). Future cardiac G&R studies shall include fibre structure re-orientation regulated by certain evolution mechanisms under various stimuli [65] .
To overcome the high computational cost of the classical constrained mixture model [9] , a temporally homogenized constrained mixture model has been proposed to take advantage of both the kinematic growth theory and the constrained mixture theory with a single time-independent reference configuration [28] . This model has recently been applied to an idealized LV under hypertension [7] . Our proposed model also makes use of both kinematic growth theory and the constrained mixture theory, while there exist several differences to the temporally homogenised constrained mixture model cited above. These include: (1) In our framework, the reference configuration is updated at every growth cycle, obviating the need for tracking the temporal evolution of inelastic deformations of added mass with respect to a constant geometric configuration; (2) It is unnecessary for us to adhere to the prior assumption that the incompressible remodelling tensor, through mass turnover, aligns with the same direction in the reference configuration. Rather, we define both the growth and remodelling tensors in the updated reference configuration, which may be formulated independently in accordance with the underlying mechanobiology; (3) Residual stress/stretch is explicitly included in updated reference configurations; (4) Our proposed framework has the potential to allow for geometrical change. One example of this would be that it may be possible to predict how the heart will respond to ventricular reconstructive surgery in patients with ventricular aneurysms [66] . In summary, both this updated Lagrangian constrained mixture model and the temporally constrained mixture model [7] can model generic cardiac G&R well with the specific choice depending on the specifics of the problem being examined, and the available data, etc.
Kroon et al. [14] and Lee et al. [18] both have modelled cardiac G&R with updated references. However, neither of these models included the turnover of each individual constituent due to the phenomenological nature of the kinematic growth theory. Moreover, residual stress/stretch was discarded at the beginning of each growth cycle. In other words, stress-free was assumed in all updated reference configurations in Kroon et al. [14] , Lee et al. [18] . We further investigated the impact of neglecting residual stretch and revealed that it gives rise to a more concentrically dominant final grown LV geometry, featuring larger growth ratios, EDV, and ESV compared to the outcomes obtained when accounting for growth-induced residual stress. Details can be found in Appendix C. Therefore, our results may indicate that including G&R induced residual stretch/stress is important in capturing a more physiologically realistic G&R.

Growth laws targeting mechanical homeostasis
Mechanical homeostasis in the myocardium is largely undefined, in particular in vivo . In this study, we define the homeostatic state using the average end-diastolic stretch and the average peak systolic active tension from the baseline case, which is similar to [18,46] . Other options could also be possible, such as locationdependent stretch/stress. In this regard, it is urgently essential to experimentally quantify the mechanobiological mechanism of different heart cells under physiological and pathological conditions, in order to bring cardiac G&R models into the clinic.
Currently, there is no universal consensus on how to formulate a growth law to characterize soft tissue G&R. Various growth laws exist in the literature, such as stretch/stress-driven [5] or hypertrophic hormone-driven laws [67] . Even for stretch/stressdriven laws, different components have been used, including fi-bre stretch [5] , cross-fibre stretch [46] , the trace of Mandel stress [5,21] , and the trace of Cauchy stress [6] . By testing eight different phenomenological myocardial G&R laws in a kinematic growth model, Witzenburg et al. [68] found that only three laws with multiple inputs could capture experimentally observed G&R under both pressure and volume overload, but none of them could lead to a completely homeostatic state. They further recommended that poorly coupled mechanical stimuli (end-diastolic and end-systolic stretch/stress) from multiple directions (fibre, cross-fibre, etc.) will be needed to improve the predictability of cardiac G&R models, especially under both pressure and volume overload, which is usually found in patients with heart disease. By preliminarily combining eccentric and concentric growth triggered by both the enddiastolic stretch and peak-systolic active tension, we have studied the hybrid growth using the same LV model. Because the grown LV geometry will affect subsequent LV biomechanics, the hybrid growth is not a simple addition of the pure eccentric and concentric growth. Our results have shown that the competition between eccentric and concentric growth could affect how LV adapts its function and structure. For example, cases CF-R and CF-G are more concentrically dominated while case CF-E is more eccentric, even though they start from the same LV geometry with the same passive and active properties. This further suggests that not only stretch-driven and stress-driven compete with each other, but also the coupling between myofibre and collagen G&R could lead to different LV G&R phenotypes.
It is still sparse to quantitatively describe the temporal course of LV G&R. For example, Costabal et al. [57] carried out an eightweek-long volume overload study of six pigs, and they observed a steep increase in myocyte length from week 6 to week 8 with a myofibre growth ratio of 1.26. Similar growth amounts can be found in case CF-E with hybrid growth, and cases CF-R and CF-G with eccentric growth. This might suggest one G&R cycle in our model roughly corresponds to a week based on the experimental measurements from [57] by considering 10 growth cycles corresponding to 8 weeks. Figure D.1 in Appendix D shows that both the mean diastolic fibre stretch and active tension are very close to the set-point values with some local variations. This further indicates the LV model is roughly in a mechanically balanced state after 10 growth cycles, or close to a mechanical homeostatic state. We are not aware of human cardiac G&R longitudinal data except [59] , in which apparent growth tensors were estimated in 10 days and 7 months post-MI.
The growth laws in this study can also be formulated following the constrained mixture approach according to intra-and extracellular turnover processes in the myocardium. For example, by assuming a constitutive relation of the net mass production of myocardium including both deposition driven by stress and degradation following a simple Poisson process, the growth can be enforced by penalizing the myocardial strain energy function on the spatial density [7] . Interested readers can refer to [69] for further details. It is also possible to include hormonal drivers in growth laws as proposed in Estrada et al. [67] , where the radial component of the growth tensor was determined by a cellular-level network model of hypertrophic signalling pathways that took into account both mechanical and hormonal stimuli. Furthermore, the myocardium also experiences reversal remodelling either physiologically (adaptive G&R i.e. during pregnancy or exercise) or pathologically, e.g. fibrosis post-MI. Studies have found that medications and medical devices can stop or even reverse cardiac G&R e.g. left ventricle assistance device [70] . Thus it is essential that a clinically-useful cardiac G&R model shall capture reverse remodelling. A few studies have shown that both the kinematic growth theory [17] and the constrained mixture theory [7] can be used to predict reverse G&R. Our framework can also be extended to model reverse remodelling either phenomenologically following the kinematic growth theory or the constrained mixture theory, which will be our next step of work.

Extracellular matrix G&R
ECM is a highly adaptive structure that is dynamically regulated by the local micro-environment [2] . Clinically, pathological myocardial G&R in heart wall has often been divided into three stages, they are (1) the acute stage with rapid transitory myocyte adaptation and rapid ECM degradation followed by new matrix formation, (2) the compensated stage with continued ECM turnover with potential ECM net degradation leading to further LV chamber dilation and wall thinning; (3) the decompensated failure with continued increase of ECM turnover leading to ECM net deposition, thus increased myocardial stiffness and pump failure eventually. Volume overload-induced heart failure is usually associated with wall thinning and chamber dilation, and excessive ECM deposition [2] . In pressure overload-induced heart failure, excessive deposition of collagen fibres and increased cross-linking are usually present [71] . In this study, we have only considered three simple scenarios of ECM G&R with no new deposition and degradation in both CF-E and CF-R cases, and a controlled net deposition for case CF-G with the same growth tensor as myocytes. In this aspect, all three scenarios are within the compensated stages since ECM net deposition is limited and myocardial passive stiffness does not increase significantly at the final growth cycle.
Although with simplified assumptions of collagen G&R, our results are in general consistent with experimental studies and clinical observations with some interesting results regarding collagen G&R [2,72] . For instance, CF-E with eccentric growth can limit the increase of LV chamber, wall thinning and EF reduction although its stroke volume is slightly smaller compared to the other two scenarios. In concentric growth, CF-E leads to the least reduction in LV chamber and wall thickening with the highest SV although the peak systolic pressure is the highest. This may suggest by reducing collagen G&R, the LV pump function could be closer to its healthy state than in other scenarios under either volume overload or pressure overload. While for hybrid growth, because of the simultaneous eccentric and concentric G&R, case CF-E leads to a very different phenotype of LV G&R with increased EDV and ESV, the largest SV and the least wall thickening and ESP ( < 150 mmHg). Since fibrosis is often associated with end-stage heart failure, limiting collagen G&R (deposition or degradation) would be beneficial to prolong myocardium adaption in the compensatory state. Results of case CF-E might further provide evidence and in-silico tools for anti-fibrotic interventions in patients with cardiac fibrosis [71] , a major focus in the past two decades but with limited success [73] . We expect that once the proposed G&R model is calibrated with different patterns of collagen G&R, it would have the potential to improve our understanding of the dynamic role of cardiac fibrosis in myocardial mechanics, and further aid in the development of effective therapeutic strategies by promoting the adaptive function of collagen network but not compromising myocyte contraction.

G&R-induced residual stress
It is widely hypothesized that incompatible G&R will lead to residual stress [74] . Genet et al. [54] have further demonstrated that growth-induced residual stresses were comparable with experimental studies both qualitatively and quantitatively. For that reason, we have adopted a similar approach to determine the residual stresses after G&R. The exact origin of residual stress is still an open question. For example, Cyron et al. [28] suggested that the residual stress is more likely to result from deposited pre-stresses rather than differential growth. In this study, the growth-induced residual deformation is used to take into account residual stress but not the residual stress tensor directly [53,75] . Residual deformation has been widely used in other studies [54,76,77] because it can be directly included in the multiplicative decomposition of the observed mapping F , making using of the continuum theory of fictitious configurations.
In both concentric and hybrid growth, the average circumferential residual stresses are negative at the endocardium and positive at the epicardium with transmural variations, which is consistent with experimental findings [52,78] and modelling studies [6,17,54] . Genet et al. [54] reported the fibre residual stress varied from -0.5 kPa to 0.3 kPa that was computed from an isotropic growth law in porcine hearts. Wang et al. [53] estimated residual stress at the mid-anterior free wall of a human LV using measured residual strain data from Costa et al. [52] , which varied from −0 . 2 kPa at endocardium to 1 kPa at epicardium. In a recent study, Zhuan and Luo [6] reported residual stress in a kinematic growth model of a small animal LV, and the circumferential residual stress was in a range of −2 kPa to 3 kPa. In general, quantitative agreement in circumferential residual stress can be found between both concentric and hybrid growths of our models and published results [6,54] . Non-physiological residual stress has been reported in Luo and Zhuan [6] by using the classic kinematic growth theory with a fixed-reference configuration, which was addressed by using an incremental growth law defined in the current configuration. In this aspect, the incremental G&R defined in the updated reference configuration, see Eqs. (10) and (11) , is similar to the current-growth frame proposed in Luo and Zhuan [6] .
As shown in Fig. 14 , the residual stress distribution from eccentric growth is very different from published results [52,54] , i.e. compressed in the middle wall while stretched near endocardium and epicardium in cases CF-R and CF-G. As discussed in Genet et al. [54] , Choi et al. [55] , the diastolic myofibre stretch is larger in the middle wall but not in the sub-endocardial region, thus stretch-induced G&R (pure eccentric growth) will be higher in the mid-wall compared to sub-endocardial/epicardial wall, which would lead to compression (negative τ cc ) in the mid-wall but stretch (positive τ cc ) in the sub-endocardial/epicardial wall, see Fig. 14 . The associated residual stress field under pure stretchdriven eccentric growth is considered not physiological [78] . This in turn might indicate the distributions of residual stress constraint plausible function forms of growth law.
The evolution of residual stress in the heart under pathological G&R is still unclear. Omens et al. [79] found the opening angle of a rat LV decreased during ageing, but it was unclear how the decreased opening angle would impact the residual stress because ageing is usually associated with increased tissue stiffness. Lee et al. [17] showed that residual stress increased during ventricular G&R while decreased during reverse remodelling using a computational model based on the kinematic growth theory. Zhuan and Luo [6] recently reported that the residual stress increased initially and then decreased in either eccentric or concentric growth when formulating the growth tensor in the current configuration. Nevertheless, future experimental studies are needed to quantify myocardial residual stress evolution under both physiological and pathological conditions along with the changes in myocardial stiffness, to test and validate various G&R models, and to even identify appropriate growth laws [80] .

Computational cost
We further compared the computational time of one G&R cycle between our method and the kinematic growth theory by using the same growth tensor for the ground matrix, myofibres and collagen network. We found extra 20 minutes were needed to compute F r in our current implementation compared to the kinematic growth model. We further note that higher computational cost would be needed in our framework because each cycle is a finite time G&R, whilst the kinematic growth theory could reach the final grown state by simply using the total growth tensor [20] . Thus, it is expected that the computational time of our approach will be longer than the kinematic growth theory. In a recent study, Gebauer et al. [7] reported their computational time of LV G&R using a temporally homogenized constrain mixture model. It was within 2 hours for the first 250 days of G&R on one Intel Xeon E5-2630 "Haswell" (16 cores, 2.5 GHz, 64 GB memory) using an idealized LV geometry with 3970 quadratic tetrahedral elements. Note the LV model in this study was discretized with 133,042 linear tetrahedral elements, 30 times denser compared to the LV mesh in Gebauer et al. [7] . Translating cardiac G&R models into the clinical setting will require real-time prediction, existing 3-dimensional G&R models rarely can be solved in minutes, especially considering the time needed for model construction and calibration. Therefore, alternative approaches are needed for real-time predictions, such as reduced-order models [81] , or statistical emulators [19,82,83] .
In this study, all simulations finish in 10 growth cycles. We are aware of three simulations that do not meet the prescribed convergence criteria (1% difference of the average incremental growth ratio between consecutive growth cycles), they are cases CF-E with concentric and hybrid growth, and case CF-G with hybrid growth. The average incremental growth ratios to the previous cycle for the three cases are 1.013, 1.017 and 1.012, respectively, close to the 1% difference criteria. We further ran the case CF-E with hybrid growth for 20 cycles because it has the highest average incremental growth ratio. We found it can meet the criteria at the 20th cycle, while the differences in LV pump function between the 10th cycle and the 20th cycle were not that significant, i.e. 4% in ESP, 3% in EDV, and 7.6% in ESV. Hence, we do not expect further simulations with extra growth cycles will change the results.

Further limitations
In addition to what has been discussed previously, we note further limitations here. The LV geometry was constructed from a healthy volunteer at early diastole but not a patient LV under pressure or volume overload. Also because of non-zero pressure loading at early diastole, the LV geometry is not a truly unloaded geometry with zero residual stress. Various methods have been developed to determine unloaded geometries from a loaded in vivo configuration, such as the so-called inverse elastic problem [84,85] , or the inverse displacement approach [53,54] . Hadjicharalambous et al. [86] studied how different choices of the unloaded reference configurations influence personalized cardiac models, including image-derived geometry at end-systole, inversely estimated geometry from end-diastole, and unloaded in vitro geometry. Their results showed that all three configurations could lead to consistent pressure-volume loop and myofibre stress-strain responses despite variations in absolute values, with little qualitative effect on myocardial active contraction. They further suggested that imagederived reference configuration (i.e. at end-systole) could be used in personalized cardiac modelling for reliable comparative analysis provided that all models employ the same reference configuration. In fact, the early-diastole LV geometry is very close to end-systolic geometry because early-diastole is the time when the mitral valve first opens just after LV iso-volumetric relaxation. Future studies shall investigate how different reference configurations influence myocardial G&R.
At each grown state, we have assumed the orthogonal fibresheet-normal ( f − s − n ) layered architecture is preserved. This layered architecture has been predominately reported from various experimental studies [30,31] since its early description by LeGrice et al. [29] . However, the orthogonal f − s − n may become dis-rupted in the diseased state, i.e. fibrosis after myocardial infarction. In this regard, longitudinal studies of myocyte and collagen organizations under normal ageing and pathological conditions would be very necessary to provide insights into cardiac function and dysfunction, and also be great of interest to the cardiac modelling community.
We have mainly focused on myocyte G&R by assuming no growth in the ground matrix and simplified G&R scenarios in collagen fibres. While the actual pathological cardiac G&R is a multifactor driven process, stress/strain factors together with hormonal and inflammatory signals for different constituents need to be incorporated in G&R laws [72] . Other parts of the heart are also not considered in this study, like the right ventricle, the pericardial wall, the atria and valves, etc. Hemodynamic loads also do not adapt over time. Furthermore, myocardial G&R not just occurs biomechanically but also involves G&R in electrophysiology and electro-mechanics coupling, etc. However, such a model will be extremely challenging both mathematically and numerically.

Conclusion
In this study, we have developed an updated Lagrangian constrained mixture theory framework to model pathological cardiac G&R with updated reference configuration at each growth cycle. This new framework can integrate critical mechanisms of individual biological constituent's adaptations to altered biomechanical environments, which are the ground matrix, myofibres and collagen network in this study. Eccentric and concentric growth, and their combination have been explored in a patient-specific human left ventricular geometry under volume and pressure overload. Eccentric growth is triggered by the overstretching of myofibres with increased end-diastolic pressure, whilst concentric growth is driven by excessive contractile stress caused by aortic stenosis. We further investigate the interactions between collagen and myofibre G&R. Results have shown that this constrained mixture theory motivated cardiac G&R model can capture different phenotypes of maladaptive LV G&R under pathological conditions, consistent with experimental studies and clinical observations, i.e. chamber dilation and wall thinning/thickening. This cardiac G&R model has the potential to serve as a basis for developing more advanced predictive models informed by patient data to assess heart failure risk, predict disease progression, select the optimal treatments, and eventually towards a truly precision cardiology using in-silico models.

Data availability
The datasets supporting this article have been uploaded to GitHub as part of the electronic supplementary material, https: //github.com/HaoGao/UpdatedLagrangianCMT.git .

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. To reduce growth in myofibre direction, τ f was increased from 0.8 to 3.2, and the other parameters in Table 4 remain the same. The updated parameters for hybrid growth are further listed in Table B1 . Because of reduced myofibre growth, all three cases ap- pear more concentrically dominated growth with much increased systolic pressure as shown in Fig. B.1 (a-c). In particular, case CF-E has much less chamber dilation compared to the result in Comparing to Fig. 12 in Section 3.3 , all cases have smaller θ f with nearly the same values of θ s . It may suggest stretch-driven and stress-driven compete in hybrid growth, and such competing effects could lead to different phenotypes of LV G&R compared to pure eccentric or concentric growth.

Appendix C. The hybrid growth case CF-R with zero residual stress
We set pre-stretch tensor F i r = I for each component at the updated reference configuration of case CF-R with hybrid growth, thus the residual stress is zero at the beginning of each G&R cycle. Zero residual stress in the updated reference configuration has been applied in Kroon et al. [14] , Lee et al. [18] . Figure C.1 (a) shows the final grown LV geometry for case CF-R with hybrid growth. The  14) ) against the biaxial and simple shear data measured by Sommer et al. [32] along (a) the mean fibre direction (MFD) and (b) the cross-fibre direction (CFD), and (c) the simple shear stress in the n-f plane along the myofibre direction, i.e. (nf) case.

Table A1
Passive parameter estimated by fitting to the experimental data of one shear response (nf) and one biaxial stretch (ratio of 1 MFD : 1 CFD).
Estimated material parameters by fitting to experimental data a g (kPa) b g a m (kPa) b m a cf (kPa) b cf a cs = a cn (kPa) b cs = b cn 0.7054   average growth ratios along myofibre θ f and cross-fibre θ s directions are shown in Fig. C.1 (b) and (c), respectively. Both growth ratios are higher than the same model with residual stress at the reference configuration, 1.26 (without residual stretch) vs. 1.2 (with residual stretch) for θ f , and 1.5 (without residual stretch) vs. 1.31 (with residual stretch) for θ s . P-V loops within the 10 growth cycles are shown in Fig. C.1 (d), which gradually shift to the right with increased ESP, EDV and ESV, but opposite when the residual stretch is included in the updated reference configuration as shown in Fig. 12 (a).

Appendix D. Diastolic myofibre stretch and peak active tension in the last growth cycle
Mechanobiological stability or biomechanical homeostasis is assessed by comparing the maximum myofibre diastolic stretch and peak active stress (the trace of active stress) in the 10th growth cycle with corresponding threshold values estimated from the baseline case. As can be seen in Fig. D.1 , the average diastolic λ max i and p max i are very close to the chosen threshold values. Thus, our LV models are roughly in a biomechanical homeostasis state though locally some regions may still experience G&R.

Appendix E. Algorithms
E1. Numerical scheme for Eqs. (10) and (12) Equations (10) and (12) are discretized with finite difference for the first-order material time derivative. Following the implicit Euler backward-type time stepping schemes, Eq. (10) is reformulated in terms of the discrete residual R f with respect to the unknown fibre grown θ f at ( n + 1 ) growth cycle, Its linearization gives the tangent for local Newton iteration K = d R f / d θ f n +1 . Details of numerical implementation of Eq. (10) can be found in Algorithm 1. Here we set t = 1 , and Newton iteration stops when R / K is less than 10 −9 . In fact t is for the convenience of computation rather than the real time, its choice can be arbitrary as long as the solution converges. In other words, we assume for each G&R cycle, maximum potential G&R will be achieved. Equation (12) is solved in the same way. See Algorithm 1 for details.