A Chemomechanobiological Model of the Long-Term Healing Response of Arterial Tissue to a Clamping Injury

Vascular clamping often causes injury to arterial tissue, leading to a cascade of cellular and extracellular events. A reliable in silico prediction of these processes following vascular injury could help us to increase our understanding thereof, and eventually optimize surgical techniques or drug delivery to minimize the amount of long-term damage. However, the complexity and interdependency of these events make translation into constitutive laws and their numerical implementation particularly challenging. We introduce a finite element simulation of arterial clamping taking into account acute endothelial denudation, damage to extracellular matrix, and smooth muscle cell loss. The model captures how this causes tissue inflammation and deviation from mechanical homeostasis, both triggering vascular remodeling. A number of cellular processes are modeled, aiming at restoring this homeostasis, i.e., smooth muscle cell phenotype switching, proliferation, migration, and the production of extracellular matrix. We calibrated these damage and remodeling laws by comparing our numerical results to in vivo experimental data of clamping and healing experiments. In these same experiments, the functional integrity of the tissue was assessed through myograph tests, which were also reproduced in the present study through a novel model for vasodilator and -constrictor dependent smooth muscle contraction. The simulation results show a good agreement with the in vivo experiments. The computational model was then also used to simulate healing beyond the duration of the experiments in order to exploit the benefits of computational model predictions. These results showed a significant sensitivity to model parameters related to smooth muscle cell phenotypes, highlighting the pressing need to further elucidate the biological processes of smooth muscle cell phenotypic switching in the future.

The tonicity of vascular smooth muscle cells (SMC) after clamping was also studied extensively (Barone et al., 1989;Famaey et al., 2010;Geenens et al., 2016b). Experimental data on the effects of arterial clamping were collected in mice (Geenens et al., 2016a,b). In this study, descending thoracic aortas were clamped at different levels of loading. Then, the aorta was either excised immediately or excised after a fixed duration of healing. After excision, rings were cut and tested with a myograph to measure the vascular tone under vasoconstriction and vasodilatation stimulations, followed by histological analyses. An acute decline of endothelium-dependent vasodilatation was observed just after clamping, but the functional response was restored after 1 month (Geenens et al., 2016b). Arterial clamping was also followed by an inflammatory response leading to some degree of fibrosis.
The role of mechanobiology in the response to arterial clamping is not clearly understood. It is known that in many conditions, vascular remodeling is mediated by the mechanical stimuli sensed by vascular SMCs, permitting to maintain wall stresses at homeostasis (Humphrey, 2002). SMCs modulate their phenotype in response to changing local environmental cues (Epstein et al., 1994), possibly performing biosynthetic, proliferative, and contractile roles in the vessel wall (Thyberg et al., 1995). Contractile SMCs react to environmental changes on the short term by contracting and relaxing to restore a homeostatic state. On the longer term, biosynthetic vascular SMC produce, and degrade the extracellular matrix, thus enabling growth and remodeling (Owens et al., 2004).
In order to decipher the role of mechanobiology in the response to arterial clamping, in silico predictive models can be helpful. A number of computational models for damage through overloading of soft tissues have been developed and tested by Balzani et al. (2006Balzani et al. ( , 2012, Rodríguez et al. (2006), Gasser (2011), Peña (2011, Sáez et al. (2012), Famaey et al. (2013), Forsell et al. (2013), Schmidt and Balzani (2016), and Li and Holzapfel (2019). Most of these models are based on continuum damage mechanics (Kachanov, 1986;Simo and Ju, 1987), where the amount of damaged tissue is determined by a damage parameter. These models were successfully applied to predict acute damage after arterial clamping. However, most of them focused on short term fiber damage and modeled neither the active behavior of vascular SMCs nor the healing process occurring on the longer term.
Modeling vascular healing is also rather recent. Comellas et al. (2016) presented a computational model of tissue healing after mechanical overload, in which temporal evolutions of damage are homeostasis-driven. However, no discrimination was made between the different tissue constituents (elastin, collagen, cells) in terms of damage and mechanical behavior. This can be addressed by microstructurally-motivated growth and remodeling models based on the constrained mixture model introduced by Humphrey and Rajagopal (2002). In the constrained mixture theory, the different constituents of the tissue are constrained to move together in a mixture but all have different biologically relevant stress-free states. Tissue remodeling is governed by laws of production and degradation for each constituent based on stress states. This type of model has been used to predict different tissue adaptations such as aneurysm growth for instance by Baek et al. (2006), Alberto Figueroa et al. (2009), Watton and Hill (2009), Zeinali-Davarani and Baek (2012, Valentín et al. (2013), Cyron et al. (2016), Braeu et al. (2017), Famaey et al. (2018), Latorre and Humphrey (2018b), and Mousavi et al. (2019) or wound healing by Zuo et al. (2020). However, to the best of our knowledge, the constrained mixture theory has never been used to model healing after arterial clamping.
In the present work, we aim to computationally capture the mechanobiological effects of arterial clamping. Therefore, we introduce a chemomechanical model in a constrained mixture framework, considering inflammation, collagen deposition, SMC proliferation, SMC active response as well as SMC switch from contractile to synthetic phenotype, all depending on the mechanical and chemical environment. After introducing the details of the model, we simulate the response to arterial clamping after 1 and 2 months of healing and compare the results to experimental data.

Mouse Experiments
As reported by Geenens et al. (2016a), 108 wildtype mice were subjected to a surgical procedure, in which the descending thoracic aorta was isolated and clamped in vivo with a nonserrated, 2 mm wide clamp at either a loading level of 0.0 N (control group), 0.6 N or 1.27 N. The clamped tissue was then either immediately excised, or in vivo healing was allowed for 6 h, 2 weeks, or 1 month. After these four time points, histological analyses were carried out to assess the structural integrity of the tissue through CD105, CD45, Verhoeff 's-Van Gieson, and osteopontin-α-SMA stainings. After the immediate excision or after 1 month, myograph tests were carried out to assess the functional integrity of the tissue. The aorta segment was mounted onto two rods in an organ bath and, upon stretching of the tissue, a stable pre-load of 20 mN was reached. Afterwards, the vasoactive substances PE, ACh, and SNP were subsequentially added to the solution to assess endothelium dependent and independent vasodilation. In total, all mice that underwent surgery were divided into eighteen groups corresponding to a particular condition, depending on the clamping force and the healing time, and on the type of assessment, i.e., histology or myograph. More details on these animal experiments are given in Geenens et al. (2016a).

Passive Material Behavior
The anisotropic and nonlinear passive mechanical behavior of arterial tissue is often represented by a Gasser-Ogden-Holzapfel (Gasser et al., 2006) hyperelastic formulation. The deviatoric strain energy function is decomposed in an isotropic Neo-Hookean part, representing the elastin fibers in the tissue, and an exponential, anisotropic part, representing two collagen fiber families running in two symmetric directions. Assuming a fully incompressible material and ignoring the volumetric contribution, the strain energy function of the elastin and collagen contribution is respectively written aŝ (1) where C 10 and k 1 represent the stiffness of elastin and collagen. k 2 determines the exponential collagen behavior and κ quantifies the fiber dispersion.Ī elas 1 andĪ coll 1 are the first invariants or traces of the deviatoric right Cauchy-Green stretch tensors tr(C elas ) = tr(J −2/3 F elas T F elas ) and tr(C coll ) = tr(J −2/3 F coll T F coll ), where F elas and F coll are the deformation gradients of elastin and collagen respectively and J is the Jacobian of the deformation gradient F. More information on these different deformation gradients follows in section 2.2.4.Ī coll 4 andĪ coll 6 are the fourth and sixth invariants ofC coll and M i , representing the stretch along the preferred fiber direction, written as with M i the undeformed fiber vector defined by the fiber angle α i with respect to the circumferential direction. Therefore, M i = [0 cos α i sin α i ] T , assuming that the radial direction is the first direction, the circumferential direction the second and the axial the third.

Active Material Behavior
Contractile SMCs in the media actively generate vascular tone. An active component to the strain energy function, as described by Murtada et al. (2010) and used by Famaey et al. (2013) takes the formˆ where µ smc is a stiffness-like material parameter, n 3 and n 4 together are the fractions of the smooth muscle filaments in the force-producing states. u rs represents the normalized sliding between the filaments arising from the difference between the stress in the surrounding matrix P mat and the driving stresses of the cross-bridges of the filaments P smc . Murtada et al. (2010) give an in-depth explanation of these variables. This is also further elaborated in section 2.5.Ī smc 4 is the fourth invariant of M smc andC smc = J −2/3 F smcT F smc , the deviatoric right Cauchy-Green stretch tensor of the smooth muscle fibers associated with the smooth muscle deformation gradient F smc .Ī smc 4 can be written similarly to Equation (2), where M smc represents the orientation of the cells. Assuming that the cells are aligned along the circumferential direction, we write M smc = [0 1 0] T .

Strain Energy Function
Similarly to Famaey et al. (2018), the overall strain energy density stored in the material is calculated with a mass-averaged rule as where ρ elas k represents the elastin density at the current time step k. ρ csmc k is the density of SMCs in their contractile phenotype and ρ coll k, τ is the density of collagen cohort τ . These considered densities and constituent specific strain energy densities relate to the reference configuration. The deposition of collagen is discretized, such that different collagen cohorts can be identified, depending on the time of production. We consider all the initially present collagen as one cohort deposited at k = 0. On top of that, at every discrete time step, one cohort per collagen family is produced. At every time step, all existing cohorts, for example previously deposited at time step τ , are degraded through a slowly decaying survival fraction. In the present study, we consider two symmetric fiber families, each divided into k + 1 cohorts at every time step k.

Deformation Gradient
The strain energy defined in Equation (4) depends on the deformation gradients of the considered constituents. According to the constrained mixture growth and remodeling theory, the total deformation gradient of elastin F elas is written as where G elas is a deformation gradient containing the deposition stretches of elastin in the in vivo homeostatic reference state of the artery and F represents any deformation of the mixture as a whole with respect to this reference. The total deformation gradient of a certain collagen cohort τ is G coll represents the deformation of the collagen cohort at deposition. F coll dep (τ ) is the deformation of the mixture at the time of deposition with respect to the homeostatic reference state and F is the current deformation. In a steady state regime, the deformation at deposition of all collagen still present is equal to the current deformation, such that F coll dep = F and that F coll is simply equal to G coll , the collagen deposition stretch tensor. Collagen is assumed to be deposited at a constant stretch g coll (Bellini et al., 2014) along the main fiber direction. Therefore, for a particular fiber direction M (Cyron et al., 2016), Contractile SMCs are assumed to only feel the deformation with respect to the state at which they were deposited. Therefore, their deformation gradient is where F csmc dep (τ ) is the deformation gradient of the mixture at the time of deposition τ of the considered cohort.
All deformations are considered fully incompressible. Moreover, volumetric changes due to mass addition or loss are neglected, such that no deformation is observed as a result of growth. Figure 1 gives an overview of the considered damage effects. A short-term damage model for contractile SMCs and collagen inspired by Famaey et al. (2013) and Balzani et al. (2006) is considered. The fraction of damaged cells is modeled as a damage parameter d csmc , calculated as

Damage Model
where m csmc is a damage constant and where λ θ θ is the local circumferential stretch with respect to the in vivo reference stretch, assuming that the deformation gradient is known in a predefined local coordinate system whose axes are aligned with the local radial, circumferential, and axial directions. The fraction of damaged collagen d coll becomes where again m coll is a constant. ζ is calculated as the difference between the current and homeostatic fiber stresses (see also Equation 13). We assume that endothelium can be damaged as a result of overloading of the inner elastin of the media, since the endothelium itself bears almost no load. As stated by Jufri et al. (2015), endothelial cells react differently to physiological and pathological ranges of mechanical stretch, where the latter may induce apoptosis (Kou et al., 2009). We therefore assume that the local endothelium dies if a certain threshold m ec of the local β is exceeded.

Remodeling Model
A remodeling algorithm is defined considering six main components. In the following sections, remodeling pathways are introduced for the two main passive load-bearing constituents elastin and collagen. Two SMC phenotypes are considered, active load-bearing contractile cells and non-load-bearing synthetic cells that produce extracellular matrix. We also consider the healing of the endothelium and the infiltration of inflammatory agents. The scheme on Figure 2 is an overview of the remodeling pathways with the corresponding equations introduced in the following sections.

Elastin
We assume that the production or degradation of elastin is negligible over the considered time frame and that new elastin cannot be produced. The elastin density at each time point is therefore equal to the homeostatic density ρ elas 0 .
FIGURE 1 | Schematic representation of the damage effects presented in section 2.3.
Frontiers in Bioengineering and Biotechnology | www.frontiersin.org
m coll (τ ) represents the amount of collagen of the specified cohort at the time of deposition τ and q coll k, τ is the fraction of this cohort that survives until time k.
The degradation of collagen depends on the current fiber stress. Upon discretization of Equation (8) in Famaey et al. (2018) or (53) in Valentín et al. (2013), the survival fraction of a certain collagen cohort is (13) where K coll qh is the homeostatic decay constant and ζ k,τ represents the difference between the current and homeostatic fiber stresses as defined in Valentín et al. (2013) and Famaey et al. (2018).
The production of new collagen cohorts is proportional to the current density of synthetic cells. The rate at which they produce collagen depends on the presence of contractile cells and on the mechanical stimulus λ felt by these latter cells. The production rate at time τ is written as where ρ csmc and ρ csmc 0 are the current and homeostatic densities of contractile cells and ρ ssmc and ρ ssmc 0 are the corresponding densities of the synthetic cells. K coll m is a remodeling parameter and the mechanical criterion for remodeling is written as Note that λ is very similar to β csmc . It can however be positive or negative for circumferential stretch or compression, respectively.

Contractile Smooth Muscle Cells
Contractile SMCs dedifferentiate into synthetic cells upon mechanical triggering, for example as observed by Wang et al. (2018) or when losing grip to the surrounding extracellular matrix. We assume that these cells react to stretch in the circumferential direction as Through numerical integration, this becomes where K smc dd is a rate parameter for cell differentiation and t the considered time step. Regardless of this equation, the maximum relative amount of contractile cells ρ csmc /ρ csmc 0 is bounded by the relative amount of elastin ρ elas /ρ elas 0 since we assume that cells cannot be contractile if they are unable to grip the extracellular environment. The increase of contractile SMCs is also bounded by the available amount of synthetic cells to differentiate from.

Synthetic Smooth Muscle Cells
Whereas, contractile cells are quiescent in their normal state, synthetic SMCs are more proliferative (Hao et al., 2003). We therefore assume that their density can increase through dedifferentiation or proliferation based on the mechanical environment (Mantella et al., 2015). Moreover, their proliferation also increases as a reaction to inflammation (Yang et al., 2018). We write the evolution law of these cells as dρ ssmc dt = λK smc pl + K smc ic φ ic ρ ssmc + λK smc dd ρ csmc . (19) Therefore, in discretized form, the current synthetic cell density is where K smc pl and K smc ic are rate parameters related to the mechanical and inflammatory stimulus, respectively and φ ic is the current fraction of inflammation.

Endothelial Cells
After degradation, the endothelium heals following the logistic growth law in which K ec is a rate parameter and φ ec is the total fraction of endothelium present. φ ec = 1 means that the endothelium is fully recovered. If φ ec = 0 and no endothelial cells are present at all, no recovery is possible.

Inflammation
We model inflammation provoked when platelets and leukocytes adhere to the de-endothelialized artery and send inflammatory agents in the tissue. We therefore assume that the inflammation is directly related to the fraction of intact endothelium: where φ ic represents a relative level of inflammation with a maximum of 1.

Contractility Model
Equation (3) describes the energy generated by SMCs. As stated before, this energy depends on the muscle filament sliding and on the fraction of the filaments in their force-producing states. Similarly to Murtada et al. (2010) and Famaey et al. (2013), the driving equation for the evolution of the relative sliding of the myofilaments u rs iṡ with and for P mat < κ c n 3 P mat , for κ c (n 3 + n 4 ) ≥ P mat ≥ κ c n 3 κ c (n 3 + n 4 ) , for P mat > κ c (n 3 + n 4 ) In a steady-state or homeostatic condition, u rs is evolved to a situation where P smc = P mat . Mathematically, three situations can be discerned. Potentially, the u rs value from a previous state already allows κ c (n 3 + n 4 ) ≥ P mat ≥ κ c n 3 to be true in the new steady-state, such that P smc = P mat already holds, as can be seen from Equation (26). u rs then does not evolve further in the new steady-state. The previous u rs may also cause P mat to be smaller than κ c n 3 . In that case, u rs evolves until P mat = κ c n 3 . At the final steady-state, u rs is then written as u rs = κ c n 3 µ smc (n 3 + n 4 ) Alternatively, if the previous u rs causes P mat to be greater than κ c (n 3 + n 4 ), u rs evolves to The steady-state configuration of the muscle filaments can therefore be calculated at any level of deformation, quantified byĪ 4 . As defined by Hai and Murphy (1988) and Murtada et al. (2010), the myofilaments switch between their states n 1 , n 2 , n 3 , or n 4 by the set of differential equations where n 1 and n 2 represent the fractions of myofilaments in their detached state, while n 4 and n 3 represent the fractions of attached filaments, dephosphorylated, and phosphorylated, respectively. As explained by Murtada et al. (2010), the rate parameters k 1 and k 6 are dependent on the calcium concentration [Ca 2+ ] using Michaelis-Menten kinetics as The second equation represents the formation of Calcium-Calmodulin complex, where α Ca is a positive constant. K CaCaM in the first equation, inspired by Yang et al. (2003) is a CaCaMdependent phosphorylation rate parameter.
We assume that [Ca 2+ ] represents the intracellular calcium contraction. This concentration can be influenced by vasoactive agents, such as the vasodilator NO and the vasoconstrictor PE. The response to these agents is normalized with respect to the maximal possible response by the Hill equation: and where [Ca 2+ ] hom is the homeostatic intracellular calcium concentration and α PE is the maximum extra calcium concentration in response to PE. α NO is the maximal calcium concentration that is removed in response to NO. Plugging the resulting [Ca 2+ ] back into Equation (30), a NO and PE dependency of the rate parameters k 1 and k 6 is observed. k 2 and k 5 , that define the rate of dephosphorylation, are also directly affected by the NO concentration since NO activates myosin light chain phosphatase (Carvajal et al., 2000). We write where k 2,hom is the homeostatic rate of dephosphorylation and α 2 is the maximal increase in response to NO. NO is produced by a healthy endothelium in response to, for example, the vasodilating agent acetylcholine (ACh) or wall shear stress (WSS). Cohen et al. (1997)  In summary, the whole dependency of the contractile state of the contractile cells on the vasoactive agent PE, NO, and ACh is schematically represented in Figure 3.

Finite Element Model
The material, damage and remodeling models explained above are used for an in silico reproduction of the experiments carried out by Geenens et al. (2016a). A finite element model is set up in Abaqus/Standard 2017 to represent a mouse aorta. The diastolic geometry of the aorta is represented as a half cylinder with inner diameter 0.65 mm and thickness 0.04 mm (Bersi et al., 2016). Due to symmetry, only a length of 0.04 mm of the cylinder is modeled. The geometry consists of 12,852 full integration, hexahedral, hybrid elements (C3D8H). The simulation goes through a number of steps explained below, according to the steps followed in the actual experiments. An overview is also shown in Figure 4.

Homeostasis
In order to model the in vivo, mechanobiologically homeostatic condition of the mouse aorta, the prestressing algorithm explained in Mousavi and Avril (2017), Famaey et al. (2018), and Maes et al. (2019) is used. This algorithm looks for a suitable deposition stretch deformation gradient for elastin G elas in order to balance the diastolic in vivo reference geometry with the intraluminal diastolic pressure p = 10 kPa, while the top and bottom of the arterial section are fixed in axial direction. The collagen deposition stretch g coll and the axial elastin deposition stretch g elas ax are fixed as prior knowledge. This simulation is shown in Figure 4 as step 1.

Clamping
As shown in step 2 of Figure 4, the clamping of the aorta is simulated using two undeformable parallel plates, similarly to Famaey et al. (2013) and as shown Figure 4. Self-contact of the inner surface of the artery is defined, as well as contact between the plates and the outer surface of the artery. The plates first move toward each other until reaching the desired clamping force, while no damage to the material is allowed. Two different clamp forces are discerned: 0.6 or 1.27 N per 2 mm length of the clamp, respectively named load 1 and load 2. In a few next steps, the clamps are held at a constant distance and the damage model explained in section 2.3 is activated. When the damage to the FIGURE 3 | Schematic representation of the pathways presented from Equations (29) to (35). A green arrow with a plus sign represents a positive influence, a red arrow with a minus sign represents a negative influence. endothelium and contractile SMC has stabilized, the damage is held constant again, while the clamp plates are removed. During the whole process the intraluminal pressure is fixed at mouse aortic level p = 10 kPa.

Remodeling
After the releasing of the clamp, the remodeling algorithm explained in section 2.4 is activated, while keeping the pressure constant. This is shown as step 3 in Figure 4. The local collagen and SMC densities are initialized based on the previously calculated local damage. The initial value of φ ec is defined as the percentage of intact endothelial layer in the entire considered segment. Therefore, this is not a locally defined variable, and the same value is attributed to every integration point.
Due to the initial loss of contractile SMC and collagen, a dilatation of the artery is observed. This non-homeostatic mechanical state drives remodeling. Every remodeling step corresponds to 1 day of remodeling. During this process, all nodes are constrained to only move in the radial direction in order to avoid excessive shearing between the layers and failure of the simulation.

Myograph Test
After 31 days of remodeling (cases R1 and R2 for loads 1 and 2), immediately after clamping (cases A1 and A2 for loads 1 and 2) or immediately after obtaining the homeostatic configuration (cases A0 and R0), a myograph test is simulated, as in the mouse experiments explained in section 2.1. For further reference, an overview of these six cases is given in Table 1.
During the simulation, first, the axial boundary condition is released, as well as the intraluminal pressure to simulated excision, as shown by step 4 in Figure 4. Step 5 depicts the simulation of a myograph experiment as explained in section 2.1 similarly to Famaey et al. (2013). An undeformable rod with a radius of 0.15 mm in Abaqus is pushed into the arterty until the approximate required preload of 0.0133 N per mm length of the arterial segment is reached, assuming that the preload was set at 0.02 N for a length of approximately 1.5 mm in the actual experiments carried out by Geenens et al. (2016a). The rod is then fixed at this position, while 10 −6 M PE, 10 −5 M ACh, or 10 −6 M NO is virtually added through field variables. A contact definition is prescribed between the rod and the inner surface of the artery. Pre-constriction due to PE and relaxation due to ACh or NO can then be observed according to the smooth muscle contraction model described in section 2.5. FIGURE 4 | Overview of all simulation steps.
Step 1: obtaining homeostatic configuration, step 2: clamping the artery while allowing damage to the constituents, step 3: remodeling of the tissue after damage, step 4: releasing intraluminal pressure and axial tension, step 5: simulation of a myograph experiment. Steps 4 and 5 are done after step 1 (cases A0 and R0), after step 2 (cases A1 and A2) and after step 3 (cases R1 and R2).

Remodeling Beyond 31 Days
In order to further examine the behavior of the remodeling model, the finite element analysis of section 2.6.3 was extended to a remodeling period of 91 days. The effect of slight adaptations to the model was investigated as well. The first adaptation is the assumption that synthetic SMCs do not redifferentiate into their contractile phenotype, such that Equations (17) and (19) only hold when λ is greater than or equal to zero. In the opposite case, only the inflammation level influences the SMC densities: The second adaptation is the assumption that collagen production is solely dependent on the amount of synthetic SMCs, while their production rate is not directly affected by the mechanical environment. Equation (15) then simplifies to (37)

Model Parameters
An overview of all used parameter values is given in Tables 2,  3. A code number from 1 to 5 is attributed to every parameter, explaining the way its value was determined. We either used an exact value from the specified reference (1), or used a representative value from the reference, when for example a range was given based on the results of tests on multiple samples (2). Some parameter values are estimated (5), or the parameter is manually fitted to the experimental myograph data (see section 2.1), either with an idea of the order of magnitude from literature (4), or without (3). Figure 5 shows the distribution of β (see Equation 10), defining the local loss of contractile SMC due to clamping. The highest damage is concentrated at the inner side of the wall at the edge of the clamp. Table 4 gives an overview of the relative collagen, contractile SMC, and endothelium content acutely after clamping at the three different clamp loads (cases A0, A1, and A2). From this table it can be concluded that the difference between clamping at 0.6 and 1.27 N is small in terms of acute damage. There is approximately 70% loss of endothelium, 9% collagen loss, and 28% contractile SMC loss in both cases. The small difference in damage is due to a minimal required clamp displacement to increase the reaction force from 0.6 to 1.27 N, yielding only small stretch differences. Table 4 also shows the situation after the simulated in vivo healing period of 31 days (cases R0, R1, and R2) using the presented remodeling model, taking into account  (1) The exact value from the reference is used.

Remodeling
(2) A representative value from the reference is used.
(3) The parameter is manually fitted. (4) The parameter is manually fitted in the same order of magnitude as the reference. (5) The parameter is estimated.
(1) The exact value from the reference is used.
(2) A representative value from the reference is used.
(3) The parameter is manually fitted.  Note that the level of inflammation is zero in the normal artery wall state, and has a maximal value of 1. Cases A0, A1, and A2 refer to the acute situation after clamping at three load levels (0.0, 0.6, and 1.27 N). R0, R1, and R2 correspond to the respective cases after 31 healing days.
cell differentiation, ECM production by synthetic cells and inflammation after clamping injury. Figure 6 shows the evolution of the total content of each constituent relative to its normal amount over a remodeling time of 31 days after damage due to clamping at 1.27 N. In other words, it shows the evolution from case A2 to case R2. Due to the similar level of damage at cases A1 and A2, as is clear from Table 4, the evolution from cases A1 to R1 resembles the one depicted in this figure.
There is an initial dedifferentiation of the cells from their contractile to synthetic phenotype due to an initial overstretching of the wall. Along, with the high initial inflammation level, this also causes the synthetic cells to proliferate, such that the collagen content increases. At about 14 days, the initial stiffness loss is compensated, and the collagen, synthetic, and contractile cell contents slowly return to their normal levels. Figure 7 shows the normalized reaction force in the simulated rod while it moves toward the pre-load position before the addition of vasoactive substances. It indicates the overall stiffness of each material for each case. There was no discernible difference between cases A1 and A2 on the one hand and cases R1 and R2 on the other hand, as can be observed in Figure 7. Furthermore, the simulations do not show differences between cases A0 and R0. Figure 8 gives an overview of the results of the simulated myograph experiments upon the addition of vasoactive substances, compared to the results obtained on mouse arteries, as explained in section 2.1. The figure shows how the isometric force changes upon addition of PE, NO and ACh. PE drives an increased phosphorylation rate k 1 = k 6 of the myofilaments through an increased intracellular calcium level, inducing a vasoconstrictive effect. NO has the reverse effect on calcium and it also increases the dephosphorylation rate k 2 = k 5 . ACh does not act directly on the contractile SMC, but triggers the endothelium to produce NO. Therefore, the vasodilating effect of ACh is smaller than that of NO.

Remodeling Beyond 31 Days
The evolution of relative collagen, synthetic cells, and contractile cells density over a remodeling period of 91 days is shown in Figure 9, for the original remodeling model (A) and two adapted models (B and C) as explained in section 2.6.5. Beyond 1 month, unnatural periodic behavior emerges when using the original model, caused by the initial extra loss of contractile SMC upon overstretching, causing an extra stiffness loss, and a delay in the increased collagen production through the proliferation of synthetic cells.
When synthetic cells do not redifferentiate into the contractile phenotype, a more rapid stabilization of the remodeling is observed ( Figure 9B). However, the loss of contractile cells due to clamping overload will never be compensated in this case.
The results of the last variant of the remodeling model are shown in Figure 9C and show an increased oscillation of the synthetic cell density. Collagen production is not dependent on a mechanical stimulus anymore, such that a bigger increase in synthetic cell count is required to restore the collagen density, and along with it, restore the homeostatic mechanical environment.

DISCUSSION
The aim of this study is to introduce a computational model predicting healing in arterial tissues subjected to mechanical overloading and damage, for instance after clamping. Three models are introduced for the in silico simulation of the experiments carried out by Geenens et al. (2016a): a damage model for clamping, a remodeling model to predict healing and a contractility model to simulate myograph experiments.
The contractility model is original as it is the first to take the vasoactive substances PE, NO, and ACh into account. Their respective influence on the rate of phosphorylation and dephosphorylation of myosin light chain leads to a reliable response in the simulation of a myograph experiment, as shown in Figure 8. The model is based on signaling pathways on the cellular level, dependent as well as independent on intracellullar calcium, as shown in Figure 3. The approach is different from the recent model presented by Murtada et al. (2016), in which the smooth muscle tone prediction was based on a structurally motivated model of the contractile unit. In their implementation, the response to an external factor, such as a change in loading or in the concentration of a vasoactive agent, is modeled as an evolving scaling factor for the myosin filament length. Before us, the continuum mechanics-based model of Murtada et al. (2017) was the only one that accounted for the dependency of the phophorylation rates on the diffusion of the vasoconstrictor potassium chloride (KCl) from the adventitia, although diffusion itself is neglected in the present study.
The remodeling model includes novel aspects of cell differentiation upon mechanical stimulus and the production of extracellular matrix by synthetic SMCs. This production is also dependent on a certain level of tissue inflammation, as for example done by Latorre and Humphrey (2018b). In the present approach however, the inflammation level is directly related to the damage and healing of the endothelium. Inflammation increases the synthetic cell proliferation, thus indirectly enhancing collagen production (Davis et al., 2003), as summarized in Figure 2. Hence, our remodeling model includes all the relevant biological processes and pathways, in contrast to more phenomenological models, where collagen turnover is directly related to a mechanical stimulus, such as in Baek et al.   This more detailed description of SMC behavior in vascular healing and remodeling comes at an increased computational cost. Moreover, Figure 9 shows stability issues of the model in the form of unnatural temporal oscillations of the densities at longer time scales. A solution could be to neglect the transient effects and only consider the steady state, such as done by Latorre and Humphrey (2018a). Alternatively, we can include damping in the model to obtain a critical or overdamped dynamic system in order to avoid unnatural periodic behavior. From a mathematical point of view, the main limitation is the high number of parameters, as summarized in Tables 2, 3. Some parameters are determined based on previous works or based on their physical meaning, others were set in order to match experimental findings, mainly based on the tissue properties at 0 and 31 days of healing, which were however not sufficient to uniquely determine the parameter values.
Unfortunately, the currently available experimental data is not sufficient to proof its pilot application. It is likely that other parameter combinations would amount to the same results as shown in Figure 8. Nevertheless, the phenomenological nature of this new model is strongly reduced as compared to state-ofthe-art models. A high number of parameters can be qualified as physics-based, such that their values can be obtained through the design of dedicated biochemical experimental set-ups. This will allow these parameter values to be measured with more certainty, or with smaller confidence intervals, capturing the individual differences and differences between tissue types, allowing a better focus of the parameter fitting process.
Constrained mixture models are generally computationally expensive due to their high memory use, inversely related to the length of the time step. To ensure the feasibility, we chose to use a time step of 1 day, where a time convergence study showed errors of <5% with respect to the situation with a time step of half a day. Also in an attempt to limit the computational cost, only a very short segment of artery is modeled and the defined boundary conditions cause a plane strain situation. Considering a longer segment, possibly along with a more realistic patient-specific geometry, would improve the reliability of the model, mainly near the edges of the clamp and near the edges of the excised sample during the myograph simulation. Using this very short artery segment allows to use a non-localized variable φ ec that represents the overall intactness of the endothelium in the segment. Localizing the endothelial damage would greatly affect the complexity of the model, since diffusion of inflammatory agents and NO would need to be integrated. Similarly, taking into account the migration of SMC as an important mechanism in vascular remodeling, would increase the complexity as the remodeling in a certain integration point would be affected, not only by all variables defined in that specific location, but also by its surroundings. In a similar way, one could also consider reendothelialization as a non-localized process of proliferation and migration of nearby endothelial cells. Including all these processes would increase the biofidelity of the model, although it is unclear to what extent, given the already many unknowns in the present version.
Furthermore, to further improve the remodeling model, an improved understanding of biological and biochemical phenomena is required. To this day, some unknowns, uncertainties and controversies remain. For example, it is unclear to what kind of mechanical stimulus cells react. There are indications that SMCs and fibroblasts have a preferred structural stiffness of the extracellular matrix and react based on deviations from this ideal value (Humphrey, 2008). On the other hand, certain signaling pathways are thought to be triggered by so-called baroreceptors, sensitive to mechanical stretch (Lacolley et al., 2017). Multiple studies have investigated the effects of cyclic straining of arterial tissue, as reviewed by Mantella et al. (2015). Some apparently contradictory results emerge. For example, Chang et al. (2003) observed an increased SMC proliferation under in vitro cyclic strain, while Morrow et al. (2005) and Guha et al. (2011) observed a decreased proliferation, potentially due to a different experimental design that mimics in vivo loading conditions better (Mantella et al., 2015). The widely accepted theory that precursor cells differentiate into synthetic cells and subsequently become fully differentiated contractile cells has been challenged recently as well, given that both phenotypes can be present in healthy tissues while maintaining vascular tone and tissue architecture (Rensen et al., 2007).
In summary, the presented models provide a detailed description of vascular SMC behavior under conditions of damage as well as at different concentrations of vasoactive agents. This allows us to study tissue healing and the effects of, for example, vasoactive or anti-proliferative drugs. However, there are still many unknowns regarding these phenomena, which is why more detailed and carefully designed experiments are needed in order to fully capture SMC behavior in all its aspects.
To conclude, in this study, a damage model, as well as a remodeling and cell contractility model were introduced, taking into account endothelial damage and healing, tissue inflammation, mechanosensing, extracellular matrix production and phenotype switching of SMCs. Using these models, in vivo clamping tests on mice aortas and subsequent healing and myograph tests, were simulated through finite element modeling. The results of the simulated myograph tests showed great resemblance to the results of the actual experiments. This detailed mechanobiological description of vascular SMC behavior can be clinically relevant to enable in silico investigations of drug effects. However, the results show that there is still a need for an improved biological and biochemical fundamental understanding to reliably capture vascular SMC mechanobiology at all the relevant spatio-temporal scales.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.