A homeostatic-driven turnover remodelling constitutive model for healing in soft tissues

1. Background

Living matter has the ability to adapt and evolve in response to disease, external loads and environmental stimuli or aggressions. In particular, tissues are capable of healing in order to arrest the extent of the damage caused by injury or disease and, ultimately, restore the tissue's original function. The repair of soft tissues is known to be driven by a complex sequence of events involving cellular processes as well as biochemical and biomechanical factors [1,2]. The exact role of many of these factors is not completely understood yet. Nonetheless, from a physiological point of view, the healing process is classified into four distinct but overlapping phases: haemostasis, inflammation, proliferation and remodelling [3,4].

Following an injury, the loss of structural integrity in the tissue immediately activates the coagulation cascade that results in a platelet-rich fibrin clot (haemostasis phase). Within 1–3 days after injury, neutrophils and macrophages are attracted to the wound site to phagocytose bacteria and debris, preventing wound infection (inflammation phase). Cytokines are also released, which stimulate angiogenesis and enhance the production of fibroblasts. The ensuing fibroblast proliferation and migration results in the synthesis and deposition of collagen in the extracellular matrix (ECM), that leads to the formation of granulation tissue (proliferation phase). This phase can last up to several weeks. In skin wounds, epithelialization and wound retraction is also observed [5]. The final stage (remodelling phase) lasts from weeks to years and consists of a continuous synthesis and degradation of collagen as the ECM is remodelled, and the granulation tissue becomes the scar tissue. As this matrix turnover takes place, its composition shifts and reorganizes: the newly formed blood vessels regress, the ‘flaws' (such as fat cells and inflammatory pockets) are removed, and the collagen fibres become increasingly organized. Over time, and under adequate biochemical and biomechanical conditions, the remodelled tissue approaches the characteristics of the original undamaged tissue. However, the completely healed scar tissue often does not fully recover the characteristics of the uninjured tissue it replaces [1]. In addition, homeostasis may never be reached, such that remodelling may also appear as a continuous pathological transformation of diseased tissues during aneurysm expansion for example [6].

The mathematical modelling of wound healing has been widely addressed since the development of the first models in the 1990s [7–9]. These models focus on the underlying cellular and biochemical mechanisms to define and simulate dermal wound contraction [10–12] and angiogenesis [13–15] from a continuum-based approach. The inflammation and proliferation phases have also been modelled using a discrete or a hybrid discrete/continuum approach [16–18] and, more recently, a systems-biology multi-scale and multi-field approach has been proposed [19]. The reader is referred to the works by Tepole & Kuhl [20] and Valero et al. [21] for a comprehensive review of mathematical and computational dermal wound healing models.

On the other hand, the remodelling phenomena, understood as a change in the properties of the tissue, have also been extensively addressed. In general, it is treated together with growth and not necessarily in the biological context of the above-described tissue healing process [22–25]. Many of these models aim at characterizing collagen fibre reorientation through evolving structural tensors [26–31]. A different approach characterizes growth and remodelling as a continual turnover of tissue constituents by means of a constrained mixture theory [32–35], or, more recently, a mechanistic microstructural theory [36].

Over the past few years, advancements in the field resulted in sophisticated models with cellular [37,38] and molecular [39] processes being the driving forces of remodelling. In this sense, much effort is directed towards representing remodelling in vascular tissue [40–43], with a particular focus on the pathological remodelling observed in aortic aneurysm tissue [44–46]. The mathematical modelling of the inflammation, proliferation and remodelling phases in ligament tissue has also been addressed [47,48].

Numerous studies, both in animal models and in patients, have shown that mechanical loading has a significant impact on the speed and efficiency of healing [49–52]. However, the optimal loading regime remains unclear, and the detailed mechanobiogical mechanisms involved are not fully understood. Computational approaches have been widely used in bone healing mechanobiological modelling to enable predictions of bone healing and improve the understanding of both mechanical and biological mechanisms at play [53,54]. In order to apply this approach to soft tissue healing, a continuum constitutive model is required that can represent both the changing soft tissue mechanics during healing and the proposed biophysical stimuli for the cells involved.

This work introduces a novel constitutive model that captures the continuous turnover of tissue observed in the remodelling phase that ultimately leads to the recovery of injured tissue. This homeostatic-driven turnover remodelling (HTR) model is capable of modelling the last stage in the above-described healing process, but can also capture the pathological remodelling of tissue observed in certain diseases [55], such as the abdominal aortic aneurysm (AAA) [56].

The proposed HTR model is consistent with an open-system thermodynamics framework [57]. It is formulated in accordance with the principles of continuum damage mechanics (CDM), following ideas proposed in bone remodelling by Doblaré & García [58]. Consequently, damage is assimilated to an apparent density that evolves in response to mechanical loading. In the present formulation, however, damage (and its repair) can be physically associated with the injury (and healing) observed in live tissues. Thus, the granulation tissue formed in the proliferation phase is akin to a damaged material in CDM, characterized by microfissures and voids that result in a loss of stiffness and strength. Then, through remodelling, the ‘flaws' in the granulation tissue are removed (healed) and it becomes increasingly organized, approaching the original characteristics of the ECM material. This process is comparable to a reversal of damage in a CDM framework, in which the load-carrying capacity of the tissue is gradually recovered.

The HTR model describes the overall change in material behaviour at tissue level of healing/remodelling tissues. Healing in this model is not only driven by mechanical loading, but also by biological stimuli. In particular, the underlying metabolism in healing tissues is represented by phenomenological parameters.

The thermodynamic basis and derivation of the HTR constitutive model is developed in §2. The damage evolution and healing rate equations are also outlined in this section as are the details of the numerical implementation in an in-house finite-element (FE) software [59]. Additional information for this numerical implementation is detailed in appendix A. Section 3 includes several examples with the aim of illustrating the basic constitutive characteristics of the model and validate it with experimental data from the literature. The characteristics of the model, in addition to its advantages and shortcomings, are discussed in §4. Finally, conclusions are addressed in §5.

2. A constitutive model for homeostatic-driven turnover remodelling

Soft tissue is known to be highly deformable, yet experience negligible volume changes, and exhibit a characteristic J-shaped nonlinear response [60]. For many applications, its behaviour can be described by means of decoupled quasi-incompressible hyperelastic models [61–63], with damage affecting solely the deviatoric term [64–68]. Albeit their limitations, discussed in §4, these simplifying hypotheses are also assumed in this work.

The HTR model captures the homeostatic-driven turnover remodelling of soft tissues, i.e. the last phase of tissue healing. The model is developed within the framework of CDM and is based on the thermodynamics of irreversible processes with internal state variables [69–72].

Unlike inert materials, tissues have an underlying metabolism, which is essential to the growth, healing and remodelling processes characteristic of living organisms. From a continuum mechanics standpoint, this metabolism introduces energy into the system, allowing for the ‘recovery’ of the energy dissipated during damage and, thus, permitting a ‘reversal’ of the damage produced in the material. Consequently, the total specific Helmholtz free energy density introduced into the system is the sum of the initial strain energy Ψⁱⁿⁱ contained in the tissue and the strain energy introduced by the metabolism such that 2.1

Here, the energies are given with respect to the density in the reference configuration, and the tilde indicates the deviatoric or volume-preserving part of the free energy. The subscript ‘vol’ refers to the volumetric part.

The recovery energy reverses the damage in the tissue such that the internal damage variable is no longer accumulative in nature, i.e. as in classic CDM models. Specifically, we postulate the deviatoric part of the specific Helmholtz free energy density to be of the form 2.2where is the original (undamaged) hyperelastic specific Helmholtz free energy density given in terms of the deviatoric part of the right Cauchy–Green strain tensor, . The effective damage D_eff is the internal (recoverable) damage variable, which only affects the deviatoric part of the tissue's strain energy. Its rate is given by 2.3where is the rate of , an explicit Kachanov-like mechanical damage variable, and is the rate of R, the repair or healing term. From a CDM point of view, D may be associated with the microvoids and small fissures that appear and extend as damage initiates and evolves. D = 0 corresponds to a compact material with no voids or fissures, whereas D = 1 is a completely damaged material whose amount of voids and fissures is such that it can no longer carry any load. The healing term R represents the reversal or ‘filling’ of these microvoids and small fissures such that the original load-carrying capacity of the material is recovered. Then, R = 1 corresponds to a mass deposition equivalent to the original undamaged material stiffness and, thus, coincides with a recovery energy , i.e. the initial pre-injury strain energy. Hence, the healing term can be defined in terms of the specific strain energy densities as .

Ideally, the original properties of the uninjured tissue should be recovered at the end of the healing process such that the healed tissue is indistinguishable from the pre-injured tissue. In practice, some healed tissues are softer than their corresponding healthy uninjured tissue [1,73,74], whereas others become stiffer, often loosing functionality. The latter is the case of fibrotic scar tissue [75–77], which has been associated with pathological conditions caused by an aberrant ECM production that results from perturbed homeostasis in the tissue [55]. In this work, the former case will be addressed and, thus, is assumed, i.e. at most, the original properties can be recovered. The full recovery (R = 1) corresponds to a successful restoration of the tissue's homeostatic state.

The evolution of both D and R will be defined in more detail in §§2.2 and 2.3, respectively. However, because R will be seen to implicitly depend on the tissue damage, it is anticipated that .

2.3. Healing rate

The evolution of the repair or healing variable R is defined in accordance with the biochemical and biomechanical observations of healing soft tissue. It is inferred from the description of the phases of the healing process that damage is a trigger of this process, but healing only occurs when the metabolism allows for it (figure 1). In addition, in many cases, the mechanical properties of the completely healed tissue remain inferior to uninjured tissue [1,73,74]. Based on this experimental evidence, the healing rate, 2.9is proposed. Here, represents the Macaulay brackets [80], is a function that regulates how fast healing occurs (introduces a time scale) and ξ defines the percentage of stiffness that is not recovered at the end of the healing process. Note the implicit character of the healing rate, because D_eff is a function of R. In addition, because D_eff is also a function of D, the healing rate is also implicitly dependent on the mechanical loading of the tissue.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Interpretation of the healing process in a CDM framework and contribution of the HTR model in this context. (Online version in colour.)

The irreversible stiffness loss parameter is a given value that dictates the amount of stiffness lost, with respect to the uninjured tissue's stiffness, at the end of the healing process. In other words, ξ establishes the remnant effective damage that is not recovered in the completely healed tissue. For example, ξ = 0.2 indicates that, after complete healing, the tissue will have recovered 80% of its original stiffness, namely there will remain a D_eff = 0.2.

The function regulates the healing speed, which is directly related to the system's metabolism or biological availability. Here, the biological availability is understood as the complete set of internal biochemical elements (proteins, enzymes, growth factors, etc.) necessary for healing to take place [81]. Owing to lack of experimental data and for the sake of simplicity, a constant healing rate has been defined, . The healing rate parameter k is a given value that determines the healing time scale and is measured in [time]⁻¹.

Thus, the healing rate function proposed here complies with the basic biomechanical conditions that under absence of injury (D_eff = 0) or in the case of no biological availability (k = 0 days⁻¹) healing will not occur.

3. Validation examples

The main characteristics of the HTR model are illustrated by means of a simple uniaxial tensile test example. Then, data on ligament healing taken from the literature are used to validate the model. Finally, an AAA is numerically reproduced under different healing conditions to demonstrate the applicability of the model in reproducing experimental set-ups and the capability of the formulation to analyse geometrically complex models.

3.1. Homogeneous uniaxial tension

An eight-noded cubic element with 1 cm length sides is subjected to a displacement-driven pure tensile load applied in steps of 0.1 mm, as shown in figure 2. Each load step corresponds to a time increment of 0.05 days. Two sets of hyperelastic and damage material properties have been considered (listed in table 2), one reproduces a neo-Hookean-like behaviour and the other, an Ogden-like one. A penalizer value 10⁹ times the maximum value of the shear moduli has been considered for the bulk modulus κ in all cases.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Prescribed displacements applied on an eight-noded hexahedral linear element with a single pressure integration point (Q1P0) used in the homogeneous uniaxial tensile test example.

View this table:

• View inline
• View popup

Table 2.

Hyperelastic and damage material parameters used in the homogeneous uniaxial tensile test example. The fracture energy per unit area is computed as: , where L₀ is the localization or characteristic length in the reference configuration [85,86].

In the first set of examples (figure 3), an irreversible stiffness loss parameter ξ = 0 has been used, such that the initial stiffness properties will be completely recovered by healing. The healing rate parameter k changed between 0 and 1000 days⁻¹. A high healing rate (k = 1000 days⁻¹ in figure 3) is undistinguishable from the hyperelastic model, because healing immediately compensates for the damage produced. This can be understood as a representation of the continuous turnover known to occur in living tissues. In addition, a null healing rate (k = 0 days⁻¹ in figure 3) results in a passive damage response, i.e. accumulation of damage in an inert material.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Second Piola–Kirchhoff stress versus stretch (a) and effective damage versus stretch (b) responses of the homogeneous uniaxial tensile test example using the neo-Hookean material parameters (table 2) with linear damage for an irreversible stiffness loss parameter ξ = 0 and varying values of the healing rate parameter k (values given in days⁻¹).

The next set of examples (figure 4) show the effect of varying the parameter ξ, which dictates the final effective damage in the completely healed tissue. As expected, for a value ξ = 1.0, a behaviour analogous to the passive damage model is obtained, because no stiffness can be recovered and damage continuously accumulates.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Second Piola–Kirchhoff stress versus stretch (a) and effective damage versus stretch (b) responses of the homogeneous uniaxial tensile test example using the Ogden material parameters (table 2) with exponential damage for a healing rate parameter k = 0.25 days⁻¹ and varying values of the irreversible stiffness loss parameter ξ.

Finally, a loading–unloading–reloading case is reproduced for different values of the healing variables (figure 5) to illustrate how healing may continue while unloading takes place, such that damage progression and recovery (healing) may or may not occur simultaneously.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Second Piola–Kirchhoff stress versus time (a) and effective damage versus time (b) responses of the homogeneous uniaxial tensile loading–unloading–reloading test example using the neo-Hookean material parameters (table 2) with linear damage for varying values of the healing rate parameter k (values given in days⁻¹) and irreversible stiffness loss parameter ξ.

3.2. Healing ligament

Quantitative experimental data on healing are difficult to find in the literature and, when available, are not always in a form that can be readily used and reproduced to validate numerical models. As one of the rare examples, the experimental work by Abramowitch et al. [87] on healing medial collateral ligaments (MCL) in goat knees provides excellent data to validate the HTR model. In their experiments, the MCL is surgically sectioned, and the free ends of the ligament are realigned but not sutured, leaving a gap of about 0.5 cm between the free ends [87,88]. The wound is then closed, and the animals are allowed to recover for six weeks, after which they are humanely euthanized and their knees are prepared for testing. Typical tensile stress–strain curves are provided for the healed ligament and a healthy (uninjured) ligament used as control (see grey lines in figure 6). Because there are no specific geometry and boundary conditions associated with these curves, the data have been used to calibrate material properties with the cubic element of the previous set of examples (figure 2). However, the length of the element sides has been reduced to 0.5 cm to match the experimental data provided.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Cauchy stress versus engineering strain responses to uniaxial loading of a healthy and a six week healed MCL tissue following a surgical sectioning. FE results (solid black lines) were obtained using the material properties given in table 3. The grey curves illustrate the response from the experimental data in Abramowitch et al. [87].

View this table:

• View inline
• View popup

Table 3.

Material parameters, estimated from experimental data [87], used in the MCL healing example (figure 6). The fracture energy per unit area is computed as: , where L₀ is the localization or characteristic length in the reference configuration [85,86].

A uniaxial tensile loading is reproduced in order to estimate the Ogden and damage material properties that fit best the healthy stress–strain curve. These material properties are then used in a simulation with a forced initial damage D_eff = 1, in which no load is applied but healing is allowed to progress for six weeks. An irreversible stiffness loss parameter ξ = 0.65 has been considered, because MCL scar tissue is known to regain at most 30–40% of its normal stiffness [73]. The healing rate parameter k is adjusted such that, after a six week healing period, the stress–strain curve obtained for uniaxial tensile loading fits the experimental data. Table 3 summarizes the material parameters used in this numerical example, and figure 6 compares the numerical results with the experimental data. The set of parameters used was achieved by a manual trial and error approach and is not unique nor satisfies the minimum of an objective function. A penalizer value κ = 10⁶ Pa has been considered as the bulk modulus.

3.3. Remodelling abdominal aortic aneurysm

An AAA is a permanent localized dilatation of the abdominal aorta which, if left untreated, progresses over time and can eventually rupture, leading to death. AAA rupture is a multifactorial process that involves interacting biomechanical, biochemical, cellular and proteolytic aspects. An irreversible remodelling is known to occur in the connective tissue of an aneurysm's aortic wall, characterized by a progressive imbalance between the synthesis and degradation of collagen and elastin in the ECM.

The degradation of elastin is linked to a decreased load-bearing capacity of the wall tissue, which leads to the initial arterial dilatation. A compensatory increase in collagen synthesis, associated with the overall hardening of the aortic tissue, is observed in the latter stages of AAA evolution. Beyond a certain threshold, however, the aneurysm becomes at high risk of rupture. It is believed that this final progression to rupture involves the proteolytic degradation of the tissue's collagen fibres. The reader is referred to [6,45,89,90] and references therein for further details on the many factors involved in the progression and rupture of AAAs, some of which are not completely understood yet.

The proteolytic degradation of elastin described above may be regarded, from a macroscopic point of view, as a degradation of the tissue's properties. Thus, the HTR model has the potential to characterize this particular factor in the complex evolution of AAAs, linking the pathological arterial dilatation observed in the initial stages of AAA formation to the ‘healing’ capacity of the tissue.

A three-dimensional reconstruction of an AAA was obtained through segmentation of computer-tomography images (A4research, VASCOPS GmbH [91]) and meshed using 4707 hexahedral Q1P0 elements. A single element was included across the wall thickness with an approximately constant value of 1.5 mm throughout the aneurysm. Therefore, bending effects are neglected in the simulation. The model is fully fixed at the top slice and allowed vertical displacements at the bottom one. A blood pressure of 100 mmHg (13.33 kPa) is applied in 200 load increments on the inner surface of the wall by means of a deformation-dependent follower pressure load on the face of each element. Material properties were estimated from the experimental tensile test data available in Gasser [92] using a single element (figure 2). A penalizer value κ = 10¹² Pa has been considered as the bulk modulus. The set of parameters used (listed in table 4) was achieved by a manual trial and error approach and is not unique nor satisfies the minimum of an objective function. The corresponding constitutive response is plotted in figure 7. The distal and proximal extents of the aneurysm are excluded from damage evolution, i.e. assigned the purely hyperelastic response shown in figure 7.

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

First Piola–Kirchhoff stress versus stretch responses to uniaxial loading of an AAA tissue. The FE results (black lines) were obtained using the material properties given in table 4. The grey dots illustrate the response from the experimental data provided in Gasser [92].

View this table:

• View inline
• View popup

Table 4.

Ogden and damage material parameters, estimated from experimental data [92], used in the remodelling AAA example. The fracture energy per unit area is computed as: , where L₀ is the localization or characteristic length in the reference configuration [85,86].

The example was studied with two different values of the healing rate parameter k, and an irreversible stiffness loss parameter ξ = 0 was assumed in both cases. Under non-pathological conditions, the aortic wall is continuously remodelling and, thus, for a high healing rate its behaviour should be that of a healthy tissue. Figure 8 shows the deformed shape of the same AAA at identical loading and boundary conditions but considering two different healing rate parameters: k = 0.01 (high rate) and k = 0.002 years⁻¹ (low rate). The high healing rate resulted in deformations comparable to a sole hyperelastic simulation, because damage is healed quasi-simultaneously. However, for the low healing rate, the simulation failed at a blood pressure of 71.5 mmHg (9.53 kPa). At this loading value, a high damage concentration localizes in a narrow band of elements (figure 9), which leads to structural instability and numerical failure in the following load step.

• Download figure
• Open in new tab
• Download powerpoint

Figure 8.

Deformed shape of an AAA considering the material properties given in table 4 and subjected to a blood pressure of 71.5 mmHg (9.53 kPa). Two different healing rate parameters k have been considered in addition to an irreversible stiffness loss parameter ξ = 0. The distal and proximal extents were excluded from damage and assigned a hyperelastic material behaviour. (Online version in colour.)

• Download figure
• Open in new tab
• Download powerpoint

Figure 9.

Effective damage distribution in an AAA considering the material properties given in table 4 and subjected to a blood pressure of 71.5 mmHg (9.53 kPa). Two different healing rate parameters k have been considered in addition to an irreversible stiffness loss parameter ξ = 0. The distal and proximal extents were excluded from damage and assigned a hyperelastic material behaviour. (Online version in colour.)

4. Discussion

The effective damage D_eff (2.3) drives healing in the HTR model. This variable is a direct representation of the tissue's state because it dictates the stiffness of the healing tissue, which can be measured through experimental tests. In contrast with previous remodelling models (see §1), the present description does not attempt to capture the realignment of collagen or the processes taking place at cellular or microscopic level. Instead, it is a phenomenological model that aims to describe the overall change in material behaviour (stiffness) at tissue level of a healing tissue.

The driving internal variable D_eff accounts for both mechanical and biological stimuli. Mechanical loading induces damage in the tissue as D is a function of the stress. The injury produces a biological response such that, if the metabolism allows for it, then healing occurs and the effective damage in the tissue is reversed (figure 1). The metabolism's action is quantified through the two healing parameters, k and ξ. Then, a healed tissue that has completely recovered the original properties is undistinguishable from the original tissue. The model is able to capture this, as seen in figure 5 for the uniaxial tensile loading–unloading–reloading case with k = 1.00 days⁻¹; ξ = 0.0, where the reloading curve is exactly the same as the first loading curve.

In contrast, when the healed tissue does not recover the original properties (ξ > 0), it is assumed to have a remnant damage such that it is permanently softer than the initial pre-injured tissue. This is observed in figure 5 for the cases with k = 0.25 days⁻¹; ξ = 0.2 and k = 0.25 days⁻¹; ξ = 0.5. In both cases, the reloading curves have a lower stiffness in their initial elastic portion than the corresponding portion in the loading curves. Note how a healed scar tissue that suffers additional injury will heal back to the first scar tissue properties, i.e. the stiffness loss is not accumulative over successive injuries in the same tissue.

The issue arises, then, whether this (new) healed material should maintain the updated damage threshold corresponding to the remnant damage value. An alternative would be to redefine the healed tissue as a completely new material by eliminating the remnant damage and affecting the hyperelastic (and, possibly, the damage) material properties instead. Then, the same effect would be achieved (lower stiffness), but the material would be considered simply as new and ‘undamaged’. In this case, if ξ > 0, an additional injury would result in a further reduction of stiffness in the scar tissue. This approach entails certain difficulties, namely the calculation of the new material properties owing to the nonlinear nature of the Ogden material definition. In addition, the idea of keeping a remnant damage seems to fit well with the concept of a healed tissue that has not recovered completely from injury. That is, the ECM in the healing tissue tends to reorganize and remodel towards the original configuration but does not quite achieve it. Hence, the denomination of the model as homeostatic-driven turnover remodelling.

The pathological remodelling that is known to result in fibrotic scars, which are stiffer than the original pre-injured tissue, could be accounted for by introducing negative values of the irreversible stiffness loss parameter ξ. In this way, the ‘healed’ tissue would have a final negative value of the effective damage variable, resulting in a final material behaviour stiffer than the original pre-injured one. The HTR model would no longer be homeostatic-driven and would require substantial modifications to ensure that the ‘extra’ stiffening only occurs as a result of injury. This could probably be addressed by defining a variable value of ξ in terms of the present damage in the tissue. Furthermore, the fulfilment of the Clausius–Duhem inequality (2.4) could no longer be possible, because the energy introduced into the system to remodel the tissue would be larger than the energy dissipated owing to damage.

An interesting feature of the HTR model is that there exist two completely different scales for the generation of the mechanical stimulus that produces damage (load step) and the effectivization of the biochemical part of the healing process (time step). Hence, the evolution of the mechanical damage D is dictated solely by the loading pattern imposed in the numerical simulation. Yet, the healing variable R is driven by both the load increment and the time step considered for that load increment. Then, the healing rate can be sped up or slowed down to match experimental observations independently of the loading speed imposed. Owing to this characteristic, the mechanical damage loses its physical meaning. In particular, a high value of D may be computed for a given load but, if the healing rate is high enough, the effective damage D_eff could be, in fact, practically null. As a result, a tissue can be completely healed, even when the value of D is significant. This ties in well with the fact that, in the HTR model, D_eff is the variable that describes the actual state of the tissue, as stated at the beginning of this section.

In this regard, a value D_eff = 1 corresponds to a newly formed granulation tissue, whereas D_eff = 0 corresponds to a healed scar tissue that has recovered the properties of the original uninjured tissue. This is in accordance with the simplifying hypothesis introduced in §2, namely that tissue behaviour is reproduced with a quasi-incompressible hyperelastic model, with damage affecting only the deviatoric term. The quasi-incompressible behaviour in tissues is attributed to the high volume fraction of water present in most soft tissues [22]. For supra-physiological loadings, the injury produced could potentially introduce changes in the water content of the tissue, resulting in a compressible material. In some cases, the adequacy of the quasi-incompressibility hypothesis in soft tissues subjected to physiological loading has also been debated [93,94]. The possibility of cavitational damage arising in soft tissues has also been put forth [95]. Nonetheless, the HTR model accounts for remodelling from the granulation tissue obtained in the proliferative phase of the healing process to the scar tissue resulting at the end of the remodelling phase (figure 1). Hence, a complete damage D_eff = 0 does not correspond to vacuum or inexistent tissue, but to a newly formed granulation tissue. Albeit the granulation tissue has barely any resistance to loading, it may be considered to have a fixed content of water and, thus, certain quasi-incompressibility. This would correspond to the volumetric part of the specific strain energy density, not affected by damage.

Healing is influenced by many factors such as age, severity of injury and location of the injury, among others [1,4], and the healing parameters ξ and k should account for this. At present, they are constant throughout the healing process and manually adjusted. It would be interesting to automatically adjust their value at the moment of injury, although this would require a comprehensive database quantifying the influence of the above factors on the value of the parameters. Unfortunately, the type of data required to produce this type of study is not abundant in the literature. On the other hand, the healing rate function could be made variable through the healing process. This would allow accounting for the regression of the blood vessels, i.e. the reduction of biological availability, observed in the final phase of soft tissue healing. For example, the healing rate function could be coupled to a convection/diffusion system such that the biochemical contribution to the healing rate would change as healing occurs, allowing for an adaptive biological availability distribution.

Further improvements to the HTR model include coupling it to a continuum growth model [96] to account for the tissue growth seen in hypertrophic scars. This is relatively straightforward for a volumetric growth model based on the multiplicative decomposition of the deformation gradient tensor, . Here, corresponds to the elastic part and to the incompatible growth part. The evolution of the growth stretch may be defined in terms of the mechanical loading and a function that accounts for the biological availability [81,97]. Then, the coupling of the HTR model to such a growth model would simply require determining first the behaviour owing to growth and, then, computing the healing effect on the updated grown configuration.

Nonetheless, the present HTR model is capable of reproducing experimental data on healing (figure 6) and has potential to reproduce certain characteristics of pathological remodelling, as has been exemplified in an AAA dilatation case (figures 8 and 9). This example aims at demonstrating a possible application of the HTR model to a complex three-dimensional problem. In this case, the model is shown to capture the degradation of the aortic wall's structural properties owing to the pathological remodelling in the initial stages of AAA disease.

The low healing rate case in this example reproduces the abnormal dilatation of an AAA, where the elastin degradation results in larger deformations and reduced load-bearing capacity. On the other hand, the results for the high healing rate case are comparable to a tissue with a higher structural integrity, more akin to a healthy tissue.

The present model only addresses the dilatation of the aortic wall, linked to the progressive degradation of elastin in the ECM. It does not include other known factors that contribute to the evolution and rupture of AAAs such as growth in the abdominal wall tissue, changes in phenotypes and chemomechanical responses of the cells composing said tissue, or the effect of thrombus development and maturation of the AAA, among others [45]. Even so, the inclusion of the HTR model in a general model for AAAs that accounts for the aforementioned factors could potentially contribute to better understand the complex processes involved in the evolution of AAAs.

In particular, the proteolytic degradation of the tissue's collagen fibres that has been linked to the final progression to rupture in AAA disease could be reproduced with the same HTR model. However, it would require a separate material characterization and time-scale application than that of the elastin degradation presented in the example. Introducing a modification similar to the one previously proposed for the case of fibrotic scars, the hardening effect observed in the latter stages of AAA evolution, previous to the final progression to rupture, might also be captured.

Finally, the structural instability encountered in the low healing rate case of the AAA example is attributable to the numerical limitations of the generalized damage model [79]. From a numerical point of view, in problems with negligible healing effects the HTR model is limited by stress-locking owing to the smeared approach of the damage formulation. This has been widely addressed in the literature [83,84], and a known solution to the problem is to use higher-order FE formulations. Otherwise, the HTR formulation is robust and, in any case, thermodynamically consistent.

5. Conclusion

A constitutive model for homeostatic-driven turnover remodelling in soft tissues has been presented and discussed. This model captures the stiffness recovery that occurs as a consequence of the ECM turnover observed in both the last phase of healing in tissues and the pathological remodelling of certain tissues. During remodelling, the tissue composition shifts and reorganizes, approaching the characteristics of the original undamaged material. Thus, healing is understood as a recovery or reversal of damage in the tissue, which is driven by both mechanical and biochemical stimuli. Set in a CDM framework, the driving internal variable of the HTR model is the effective damage, whose rate is the sum of a Kachanov-like mechanical damage rate and a healing rate. The former is purely driven by mechanical loading, as observed in the proposed damage softening laws. The latter is defined as an implicit healing rate, which depends on the effective damage and two healing material parameters that account for the biochemical aspects of the healing process. The model is formulated in accordance with open-system thermodynamics to account for the energy introduced into the system by the metabolism.

Numerical implementation of the HTR model is straightforward and may be particularized for any hyperelastic model. The formulation is flexible and versatile, because both the damage softening laws and the healing rate can be easily changed or modified to fit particular biological observations. However, the advantage of the phenomenological evolution laws proposed here is that they require few (damage and healing) material parameters, which is especially useful when fitting experimental data. In addition, a physical meaning can be attributed to these parameters, conferring a more functional character to the model.

Albeit the HTR model's simplicity, it has the potential to represent the active properties of complex tissues. Usage of this model in conjunction with mixture theory [98] would allow the inclusion of fibrous components and, in this way, introduce anisotropy in the overall behaviour of the tissue. Furthermore, coupling with formulations which model other biomechanical aspects such as tissue growth and necrosis [99] could result in a powerful numerical tool to represent live soft tissue behaviour.

Data accessibility

The PLCd source code containing the numerical implementation of the model is available at http://www.cimne.com/PLCd.

Authors' contributions

S.O. and F.J.B. conceived and designed the study. E.C. and T.C.G. developed the theoretical model formulation with the collaboration of S.O. and F.J.B. E.C. implemented the computational model, performed the numerical simulations and drafted the manuscript. F.J.B. collaborated in the numerical simulations and T.C.G. helped draft the manuscript. All authors revised and approved the manuscript for publication.

Competing interests

We declare we have no competing interests.

Funding

This work was partially supported by the European Research Council under the Advanced grant ERC-2012-AdG 320815 COMP-DES-MAT ‘Advanced tools for computational design of engineering materials'. E.C. gratefully acknowledges the support of the Universitat Politècnica de Catalunya (BarcelonaTECH) through an FPI-UPC scholarship.