A Cell-Regulatory Mechanism Involving Feedback between Contraction and Tissue Formation Guides Wound Healing Progression

Wound healing is a process driven by cells. The ability of cells to sense mechanical stimuli from the extracellular matrix that surrounds them is used to regulate the forces that cells exert on the tissue. Stresses exerted by cells play a central role in wound contraction and have been broadly modelled. Traditionally, these stresses are assumed to be dependent on variables such as the extracellular matrix and cell or collagen densities. However, we postulate that cells are able to regulate the healing process through a mechanosensing mechanism regulated by the contraction that they exert. We propose that cells adjust the contraction level to determine the tissue functions regulating all main activities, such as proliferation, differentiation and matrix production. Hence, a closed-regulatory feedback loop is proposed between contraction and tissue formation. The model consists of a system of partial differential equations that simulates the evolution of fibroblasts, myofibroblasts, collagen and a generic growth factor, as well as the deformation of the extracellular matrix. This model is able to predict the wound healing outcome without requiring the addition of phenomenological laws to describe the time-dependent contraction evolution. We have reproduced two in vivo experiments to evaluate the predictive capacity of the model, and we conclude that there is feedback between the level of cell contraction and the tissue regenerated in the wound.


Mechano-chemical model
The cellular species (fibroblasts and myofibroblasts) and chemical substances (growth factor and collagen) densities are obtained from a conservation law where Q denotes the cellular/chemical species, J Q denotes its net flux over the domain of interest (that may include terms due to random dispersal -migration or diffusion-, directed migration -chemotaxis-, but contains a passive convection term due to the extracellular matrix (ECM) deformation), and f Q its net production. The matrix deformation is obtained from the conservation of linear momentum where σ ecm denotes the passive resistant ECM stress, σ cell the ECM stress due to the cells-ECM adhesions and f subs the ECM-substrate anchoring forces resisting ECM deformation.

Species concentrations
Fibroblasts are connective tissue cells found in skin and the main cellular species involved in wound contraction. The main functions of fibroblasts are the synthesis of connective tissue in response to injury and the remodeling of the collagen extracellular matrix (ECM) by exerting traction forces on it [3]. Fibroblasts are motile cells whose migration is due to random dispersal, chemotaxis and passive convection caused by the ECM displacements. Hence, its net flux term can be written as where all the parameters are listed in Table 1. Fibroblasts production is due to proliferation, differentiation into myofibroblasts and differentiation back from myofibroblasts. We propose a similar expression to the one presented by [1] but with a different differentiation term. It is well known that fibroblasts sense the strains in the ECM [4], and the differentiation to myofibroblasts is regulated by mechanical loads [5][6][7]. Hence, we consider that the differentiation process is driven by the deformation of the tissue where the cells are allocated instead of depending on the mechanical stress of the matrix [8]. Moreover, the differentiation process is also enhanced by different growth factors (PDGF,TGF-β) [7,9]. Hence, the fibroblasts net production is given by the expression Fibroblasts differentiate into myofibroblasts under the influence of TGF-β and this resulting phenotype is capable of exerting and maintaining higher contraction forces in the tissue [6]. We propose a new differentiation mechanism guided by the cells volumetric strain. We consider that fibroblasts are able to differentiate into myofibroblasts only when their volumetric strain is negative, that is, when fibroblasts are under compression. In this situation cells are able to exert forces in the tissue. Myofibroblasts are smooth-muscle like cells [7], which means that they are not motile and its flux is uniquely due to passive convection and can be written as, Myofibroblasts evolution is mainly due to proliferation, differentiation from fibroblasts and inverse differentiation to fibroblasts and apoptosis [10]. Hence, myofibroblasts net production can be written as follows Cells in the skin are embedded in the ECM, whose main components are the collagen fibers, which are produced by fibroblasts and myofibroblasts. Hence, we model the ECM density through the collagen density. Collagen fibers are non motile, and hence their net flux is expressed in terms of the passive convection of the skin Following [2] we consider the role of fibroblasts and myofibroblasts on collagen synthesis [7,11]. Furthermore, collagen production is enhanced by the presence of growth factors like TGF-β [12], The wound healing process is regulated by several growth factors. Collagen-matrix contraction is regulated by PDGF [5] among others and fibroblasts differentiation is driven by TGFβ [7]. In this work, we consider a unique generic growth factor that regulates these processes for simplicity. The net flux of growth factor is due to passive convection and diffusion through the tissue. This can be written as The growth factor production is regulated by fibroblasts and myofibroblasts. Considering a simplification of the expression proposed by [2] and following [1] it can be written as where all the parameter values are included in Table 2.

Mechanosensing and mechanotransduction mechanism
This work proposes a new expression for calculating the stresses exerted by cells during wound healing.
The new expression presented in the main manuscript is where the traction force per cell, p cell , is a function of the ECM volumetric deformation, θ. In this expression we consider the stiffening effect of the ECM through p cell , that depends on the volumetric strain of the tissue [13]. Traction force per cell was introduced by Moreo et al. [13] and is defined as follows, where all the parameters and their values are included in Table 3.
In this expression the two mechanisms that a cell uses to generate the stress are considered. The actin mechanism only actuates between some deformation limits, while the passive mechanism is always activated.

Initial distributions of cells and substances
Initial species distributions have been set according to experimental observations. The wound site is occupied by hematoma tissue with low mechanical stiffness and null cell population at the beginning of the analysis. However, we consider that the healthy tissue is saturated with fibroblasts, collagen and growth factor. Nevertheless, the initial density of myofibroblasts in the healthy tissue is zero, as they appear only where there has been damage.

Constitutive relation for the skin material
The skin was assumed as a viscoelastic material [1,14]. For this kind of materials the passive resistant ECM stress σ ecm is defined where the skin's Young's modulus (E) evolves following the expression E = E 0 ρ ρ0 , that considers a stiffer ECM as the collagen density increases. All the model parameters are listed in Table 3.

Model implementation
The problem has been solved using the finite element analysis and implemented using an Abaqus User Subroutine R ⃝. Usual shape functions are used to interpolate the values of the primary unknowns (that is, n, m, ρ, c and u) from their nodal values. Time derivatives are approximated using a generalized trapezoidal method.

Weak formulation
Equations (7) and (11) can be written equivalently using Gauss' theorem to obtain the weak formulation of the problem ∫ where σ denotes the sum of σ ecm and σ cell . υ Q and υ represent the weighting functions and n ⊥ denotes the normal vector pointing outwards Ω and t := σn ⊥ denotes the tension vector. it is possible to obtain the weak formulation.

System of equations nonlinearly coupled
The model primary unknowns can be written as a function of the shape functions and their nodal values where the superscript h denotes the finite element solution. The fully discrete and nonlinear algebraic system of equations is found when substituting the above expressions in the weak form (14) and (15). We choose the weighting functions equal to shape functions. Denoting the time-dependent values of the primary variables as Z the system of equations can be expressed as a balance of internal and external forces as follows: where the subscript n + 1 denotes the time step on which the solution is being computed. Moreover, forces F int and F ext can be written in a vectorial form where The solution of the nonlinear system of equations is obtained using a standard Newton-Raphson method