Mechanical models of pattern and form in biological tissues: the role of stress-strain constitutive equations
Chiara Villa, Mark A. J. Chaplain, Alf Gerisch, Tommaso Lorenzi
MMechanical models of pattern and form in biological tissues: the roleof stress-strain constitutive equations
Chiara Villa ∗ Mark A. J. Chaplain † Tommaso Lorenzi ‡ Abstract
Mechanical and mechanochemical models of pattern formation in biological tissues have been used tostudy a variety of biomedical systems, particularly in developmental biology, and describe the physicalinteractions between cells and their local surroundings. These models in their original form consist of abalance equation for the cell density, a balance equation for the density of the extracellular matrix (ECM),and a force-balance equation describing the mechanical equilibrium of the cell-ECM system. Under theassumption that the cell-ECM system can be regarded as an isotropic linear viscoelastic material, the force-balance equation is often defined using the Kelvin-Voigt model of linear viscoelasticity to represent the stress-strain relation of the ECM. However, due to the multifaceted bio-physical nature of the ECM constituents,there are rheological aspects that cannot be effectively captured by this model and, therefore, depending onthe pattern formation process and the type of biological tissue considered, other constitutive models of linearviscoelasticity may be better suited. In this paper, we systematically assess the pattern formation potential ofdifferent stress-strain constitutive equations for the ECM within a mechanical model of pattern formation inbiological tissues. The results obtained through linear stability analysis and the dispersion relations derivedtherefrom support the idea that constitutive equations capturing viscous flow and permanent set, such asthe Maxwell model and the 3-parameter viscous model, have a pattern formation potential much higher thanthe others, such as the Kelvin-Voigt model and the standard linear solid model. This is confirmed by theresults of numerical simulations, which demonstrate that, all else being equal, spatial patterns emerge in thecase where the Maxwell model is used to represent the stress-strain relation of the ECM, while no patternsare observed when the Kelvin-Voigt model is employed. Our findings suggest that further empirical workis required to acquire detailed quantitative information on the mechanical properties of components of theECM in different biological tissues in order to furnish mechanical and mechanochemical models of patternformation with stress-strain constitutive equations for the ECM that provide a more faithful representationof the underlying tissue rheology.
Pattern formation resulting from spatial organisation of cells is at the basis of a broad spectrum of physiologicaland pathological processes in living tissues [23]. The development of mathematical models for this biologicalphenomenon started halfway through the twentieth century to elucidate the mechanisms that underly morpho-genesis and embryogenesis [29]. Since then, a number of mathematical models for the formation of cellularpatterns have been developed [58]. Amongst these, particular attention has been given to reaction-diffusionmodels and mechanochemical models of pattern formation [40]. ∗ University of St Andrews, UK ([email protected]) † University of St Andrews, UK ([email protected]) ‡ Politecnico di Torino, IT ([email protected]) a r X i v : . [ q - b i o . T O ] O c t eaction-diffusion models of pattern formation, first proposed by Turing in his seminal 1952 paper [57]and then further developed by Gierer and Meinhardt [16, 36], apply to scenarios in which the heterogeneousspatial distribution of some chemicals ( i.e. morphogens) acts as a template ( i.e. a pre-pattern) according towhich cells organise and arrange themselves in different sorts of spatial patterns. These models are formulatedas coupled systems of reaction-diffusion equations for the space-time dynamics of the concentrations of twomorphogens, with different reaction kinetics depending on the biological problem at stake. Such systems exhibitdiffusion-driven instability whereby homogenous steady states are driven unstable by diffusion, resulting in theformation of pre-patterns, provided that the diffusion rate of one of the morphogens is sufficiently higher thanthe other [28, 32, 33, 39].On the other hand, mechanochemical models of pattern formation, first proposed by Murray, Oster andcoauthors in the 1980s [44, 45, 46, 50], describe spatial organisation of cells driven by the mechanochemicalinteraction between cells and the extracellular matrix (ECM) – i.e. the substratum composed of collagen fibersand various macromolecules, partly produced by the cells themselves, in which cells are embedded [19, 20].These models in their original form consist of systems of partial differential equations (PDEs) comprising abalance equation for the cell density, a balance equation for the ECM density, and a force-balance equationdescribing the mechanical equilibrium of the cell-ECM system [42, 43]. When chemical processes are neglected,these models reduce to mechanical models of pattern formation [9, 42, 43].Over the years, mechanochemical and mechanical models of pattern formation in biological tissues have beenused to study a variety of biomedical problems, including morphogenesis and embryogenesis [8, 11, 30, 41, 43, 44,45, 46, 50, 52], angiogenesis and vasculogenesis [34, 53, 55], cytoskeleton reorganisation [1, 26], wound healingand contraction [22, 31, 49, 56], and stretch marks [17]. These models have also been used to estimate the valuesof cell mechanical parameters, with a particular focus on cell traction forces [2, 3, 4, 13, 37, 52]. The roles thatdifferent biological processes play in the formation of cellular patterns can be disentangled via linear stabilityanalysis (LSA) of the homogenous steady states of the model equations – i.e. investigating what parameters ofthe model, and thus what biological processes, can drive homogenous steady states unstable and promote theemergence of cell spatial organisation. Further insight into certain aspects of pattern formation in biologicaltissues can also be provided by nonlinear stability analysis of the homogenous steady states when long-rangeeffects are incorporated into the model [11, 26, 30].These models usually rely on the assumption that the cell-ECM system can be regarded as an isotropic linearviscoelastic material. This is clearly a simplification due to the non-linear viscoelasticity and anisotropy of softtissues [7, 21, 27, 48, 54, 59, 61], a simplification that various rheological tests conducted on biological tissueshave nonetheless shown to be justified in the regime of small strains [5, 27, 48, 59], which is the one usuallyof interest in the applications of such models. Under this assumption, the force-balance equation for the cell-ECM system is often defined using the Kelvin-Voigt model of linear viscoelasticity to represent the stress-strainrelation of the ECM [9, 43, 50]. However, due to the multifaceted bio-physical nature of the ECM constituents,there are rheological aspects that cannot be effectively captured by the Kelvin-Voigt model and, therefore,depending on the pattern formation process and the type of biological tissue considered, other constitutivemodels of linear viscoelasticity may be better suited [3]. In this regard, Byrne & Chaplain [9] demonstratedthat, ceteris paribus , using the Maxwell model of linear viscoelasticity to describe the stress-strain relation ofthe ECM in place of the Kelvin-Voigt model can lead to different dispersion relations with a higher patternformation potential. This suggests that a more thorough investigation of the capability of different stress-strainconstitutive equations of producing spatial patterns is required.With this aim, here we complement and further develop the results presented in [9] by systematically assessingthe pattern formation potential of different stress-strain constitutive equations for the ECM within a mechanicalmodel of pattern formation in biological tissues [9, 43, 50]. Compared to the work of [9] here we consider awider range of constitutive models, we allow cell traction forces to be reduced by cell-cell contact inhibition,and undertake numerical solutions of the model equations showing the formation of cellular patterns both inone and in two spatial dimensions. A related study has been conducted by Alonso et al. [1], who considered a2athematical model of pattern formation in the cell cytoplasm.The paper is structured as follows. In Section 2, we recall the essentials of viscoelastic materials and providea brief summary of the one-dimensional stress-strain constitutive equations that we examine. In Section 3,we describe the one-dimensional mechanical model of pattern formation in biological tissues that is used inour study, which follows closely the one considered in [9, 43, 50]. In Section 4, we carry out a linear stabilityanalysis (LSA) of a biologically relevant homogeneous steady state of the model equations, derive dispersionrelations when different stress-strain constitutive equations for the ECM are used, and investigate how the modelparameters affect the dispersion relations obtained. In Section 5, we verify key results of LSA via numericalsimulations of the model equations. In Section 6, we complement these findings with the results of numericalsimulations of a two-dimensional version of the mechanical model of pattern formation considered in the previoussections. Section 7 concludes the paper and provides an overview of possible research perspectives. In this section, we first recall the main properties of viscoelastic materials (see Section 2.1). Then, we brieflypresent the one-dimensional stress-strain constitutive equations that are considered in our study and summarisethe main rheological properties of linear viscoelastic materials that they capture (see Section 2.2). Most of thecontents of this section can be found in standard textbooks, such as [14, chapters 1 and 5] and [35], and arereported here for the sake of completeness. Specific considerations of and applications to living tissues can befound in [15].
As the name suggests, viscoelastic materials exhibit both viscous and elastic characteristics, and the interplaybetween them may result in a wide range of rheological properties that can be examined through creep andstress relaxation tests. During a creep test, a constant stress is first applied to a specimen of material and thenremoved, and the time dynamic of the correspondent strain is tracked. During a stress relaxation test, a constantstrain is imposed on a specimen of material and the evolution in time of the induced stress is observed [14].Here we list the main properties of viscoelastic materials that may be observed during the first phase of acreep test (see properties 1a-1c), during the recovery phase, that is, when the constant stress is removed fromthe specimen (see properties 2a-2c), and during a stress relaxation test (see property 3).1a
Instantaneous elasticity.
As soon as a stress is applied, an instantaneous corresponding strain is observed.1b
Delayed elasticity.
While the instantaneous elastic response to a stress is a purely elastic behaviour, dueto the viscous nature of the material a delayed elastic response may also be observed. In this case, underconstant stress the strain slowly and continuously increases at decreasing rate.1c
Viscous flow.
In some viscoelastic materials, under a constant stress, the strain continues to grow withinthe viscoelastic regime ( i.e. before plastic deformation). In particular, viscous flow occurs when the strainincreases linearly with time and stops growing at removal of the stress only.2a
Instantaneous recovery.
When the stress is removed, an instantaneous recovery ( i.e. an instantaneousstrain decrease) is observed because of the elastic nature of the material.2b
Delayed recovery.
Upon removal of the stress, a delayed recovery ( i.e. a continuous decrease of the strainat decreasing rate) occurs.2c
Permanent set.
While elastic strain is reversible, in viscoelastic materials a non-zero strain may persisteven when the stress is removed. 3
Stress relaxation.
Under constant strain, gradual relaxation of the induced stress occurs. In some cases,this may even culminate in total stress relaxation ( i.e. the stress decays to zero).
In this section, we briefly describe the different constitutive equations that are used in our study to represent thestress-strain relation of the ECM. In general, these equations can be used to predict how a viscoelastic materialwill react to different loading conditions, in one spatial dimension, and rely on the assumption that viscous andelastic characteristics of the material can be modelled, respectively, via linear combinations of dashpots andsprings. Different stress-strain constitutive equations correspond to different arrangements of these elementsand capture different subsets of the rheological properties summarised in the previous section (see Table 2). Inthe remainder of this section, we will denote the stress and the strain at position x and time t by σ ( t, x ) and ε ( t, x ), respectively. Linear elastic model.
When viscous characteristics are neglected, a linear viscoelastic material can bemodelled as a purely elastic spring with elastic modulus ( i.e.
Young’s modulus)
E >
0. In this case, thestress-strain constitutive equation is given by Hooke’s spring law for continuous media, that is, σ = Eε . (1)
Linear viscous model.
When elastic characteristics are neglected, a linear viscoelastic material can bemodelled as a purely viscous damper of viscosity η >
0. In this case, the stress-strain constitutive equation isgiven by Newton’s law of viscosity, that is, σ = η ∂ t ε . (2) Kelvin-Voigt model.
The Kelvin-Voigt model relies on the assumption that viscous and elastic character-istics of a linear viscoelastic material can simultaneously be captured by considering a purely elastic springwith elastic modulus E and a purely viscous damper of viscosity η in parallel. The corresponding stress-strainconstitutive equation is σ = Eε + η ∂ t ε . (3) Maxwell model.
The Maxwell model relies on the assumption that viscous and elastic characteristics of alinear viscoelastic material can be captured by considering a purely elastic spring with elastic modulus E anda purely viscous damper of viscosity η in series. The corresponding stress-strain constitutive equation is1 E ∂ t σ + ση = ∂ t ε . (4) Standard linear solid (SLS) model.
The SLS model relies on the assumption that viscous and elasticcharacteristics of a linear viscoelastic material can be captured by considering a Kelvin arm of elastic modulus E and viscosity η in series with a purely elastic spring of elastic modulus E , or equivalently a Maxwell armin parallel with a purely elastic spring. The corresponding stress-strain constitutive equation is [35]1 E ∂ t σ + 1 η (cid:18) E E (cid:19) σ = ∂ t ε + E η ε . (5)4able 1: Relations between the generic 4-parameter model (7) and the stress-strain constitutive equations (1)-(6) Generic 4-parameters model a a a b b b Linear elastic model 0 0 1 0 0 E Linear viscous model 0 0 1 0 η η E Marwell model 0 E η E η (cid:16) E E (cid:17) E η η η E η η E The 3-parameter viscous model, also known as the Jeffrey model, relies onthe assumption that viscous and elastic characteristics of a linear viscoelastic material can be captured byconsidering a Kelvin arm of elastic modulus E and viscosity η in series with a purely viscous damper ofviscosity η , or equivalently a Maxwell arm in parallel with a purely viscous damper. The correspondingstress-strain constitutive equation is (cid:18) η η (cid:19) ∂ t σ + E η σ = η ∂ tt ε + E ∂ t ε . (6) Generic 4-parameter model.
The following stress-strain constitutive equation encompasses all constitutivemodels of linear viscoelasticity whereby a combination of purely elastic springs and purely viscous dampers, upto a total of four elements, is considered a ∂ tt σ + a ∂ t σ + a σ = b ∂ tt ε + b ∂ t ε + b ε . (7)Here the non-negative, real parameters a , a , a , b , b , b depend on the elastic moduli and the viscosities ofthe underlying combinations of springs and dampers. When these parameters are defined as in Table 1, thegeneric 4-parameter constitutive model (7) reduces to the specific stress-strain constitutive equations (1)-(6).For convenience of notation, we define the differential operators L a := a ∂ tt + a ∂ t + a and L b := b ∂ tt + b ∂ t + b (8)so that the stress-strain constitutive equation (7) can be rewritten in the following compact form L a [ σ ] = L b [ ε ] . (9)A summary of the rheological properties of linear viscoelastic materials listed in Section 2.1 that are capturedby the one-dimensional stress-strain constitutive equations (1)-(6) is provided in Table 2. These propertiescan be examined through mathematical procedures that mimic creep and stress relaxation tests [14]. Noticethat, for all these constitutive models, instantaneous elasticity correlates with instantaneous recovery, delayedelasticity correlates with delayed recovery, and viscous flow correlates with permanent set.5able 2: Properties of linear viscoelastic materials captured by the stress-strain constitutive equations (1)-(6)Instantaneouselasticity Delayedelasticity Viscousflow Instantaneousrecovery Delayedrecovery Permanentset StressrelaxationLinear elastic model (cid:88) (cid:88) Linear viscous model (cid:88) (cid:88)
N. A.Kelvin-Voigt model (cid:88) (cid:88)
Maxwell model (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
SLS model (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88) (cid:88)
We consider a one-dimensional region of tissue and represent the normalised densities of cells and ECM at time t ∈ [0 , T ] and position x ∈ [ (cid:96), L ] by means of the non-negative functions n ( t, x ) and ρ ( t, x ), respectively. Welet u ( t, x ) model the displacement of a material point of the cell-ECM system originally at position x , whichis induced by mechanical interactions between cells and the ECM – i.e. cells pull on the ECM in which theyare embedded, thus inducing ECM compression and densification which in turn cause a passive form of cellrepositioning [60]. Following [43, 50], we consider a scenario where cells change their position according to a combination of: (i)undirected, random movement, which we describe through Fick’s first law of diffusion with diffusivity ( i.e. cellmotility)
D >
0; (ii) haptotaxis ( i.e. cell movement up the density gradient of the ECM) with haptotacticsensitivity α >
0; (iii) passive repositioning caused by mechanical interactions between cells and the ECM,which is modelled as an advection with velocity field ∂ t u . Moreover, we model variation of the normalised celldensity caused by cell proliferation and death via logistic growth with intrinsic growth rate r > n ( t, x ) ∂ t n = ∂ x [ n ( D ∂ x n − α ∂ x ρ − ∂ t u )] + r n (1 − n ) (10)subject to suitable initial and boundary conditions. As was done for the cell dynamics, in a similar manner we model compression and densification of the ECMinduced by cell-ECM interactions as an advection with velocity field ∂ t u . Furthermore, as in [43, 50], weneglect secretion of ECM components by the cells since this process occurs on a slower time scale compared tomechanical interactions between cells and the ECM. Under these assumptions, we describe the cell dynamicsthrough the following transport equation for ρ ( t, x ) ∂ t ρ = − ∂ x ( n ∂ t u ) (11)6ubject to suitable initial and boundary conditions. Following [43, 50], we represent the cell-ECM system as a linear viscoelastic material with low Reynolds number( i.e. inertial terms are negligible compared to viscous terms) and we assume the cell-ECM system to be inmechanical equilibrium ( i.e. traction forces generated by the cells are in mechanical equilibrium with viscoelasticrestoring forces developed in the ECM and any other external forces). Under these assumptions, the force-balance equation for the cell-ECM system is of the form ∂ x ( σ c + σ m ) + ρ F = 0 , (12)where σ m ( t, x ) is the contribution to the stress of the cell-ECM system coming from the ECM, σ c ( t, x ) is thecontribution to the stress of the cell-ECM system coming from the cells, and F ( t, x ) is the external force perunit matrix density, which comes from the surrounding tissue that constitutes the underlying substratum towhich the ECM is attached.The stress σ c is related to cellular traction forces acting on the ECM and is defined as σ c := τ f ( n ) nρ with f ( n ) := 11 + λ n . (13)Definition (13) relies on the assumption that the stress generated by cell traction on the ECM is proportionalto the cell density n and the ECM density ρ . The factor of proportionality is given by a positive parameter, τ ,which measures the average traction force generated by a cell, multiplied by a non-negative and monotonicallydecreasing function of the cell density, f ( n ), which models the fact that the average traction force generatedby a cell is reduced by cell-cell contact inhibition [40]. The parameter λ ≥ λ = 0 corresponds to neglecting the reduction in the cell traction forces causedby cellular crowding.The stress σ m is given by the stress-strain constitutive equation that is used for the ECM, which we chooseto be the general constitutive model (9) with the strain ε ( t, x ) being given by the gradient of the displacement u ( t, x ), that is, ε = ∂ x u . Therefore, we define the stress-strain relation of the ECM via the following equation L a [ σ m ] = L b [ ∂ x u ] , (14)where the differential operators L a and L b are defined according to (8).Assuming the surrounding tissue to which the ECM is attached to be a linear elastic material [40], the externalbody force F can be modelled as a restoring force proportional to the cell-ECM displacement, that is, F := − s u . (15)Here the parameter s > u ( t, x ), we apply the differential operator L a [ · ] tothe force-balance equation (12) and then substitute (13)-(15) into the resulting equation. In so doing, we find L a [ ∂ x ( σ m + σ c ) ] = − L a [ ρ F ] ⇒ L a [ ∂ x σ m ] + L a [ ∂ x σ c ] = L a [ sρu ] , ⇒ ∂ x L a [ σ m ] = L a [ sρu ] − L a [ ∂ x σ c ] , ⇒ ∂ x L b [ ∂ x u ] = L a [ sρu − ∂ x σ c ] , ⇒ L b [ ∂ xx u ] = L a [ sρu − ∂ x σ c ] , L b [ ∂ xx u ] = L a (cid:20) sρu − ∂ x (cid:18) τ nρ λn (cid:19) (cid:21) . (16)Finally, to close the system, equation (16) needs to be supplied with suitable initial and boundary conditions. We close our mechanical model of pattern formation defined by the system of PDEs (10), (11) and (16) withthe following boundary conditions n ( t, (cid:96) ) = n ( t, L ) , ∂ x n ( t, (cid:96) ) = ∂ x n ( t, L ) ,ρ ( t, (cid:96) ) = ρ ( t, L ) , ∂ x ρ ( t, (cid:96) ) = ∂ x ρ ( t, L ) u ( t, (cid:96) ) = u ( t, L ) , ∂ x u ( t, (cid:96) ) = ∂ x u ( t, L ) , for all t ∈ (0 , T ] . (17)The periodic boundary conditions (17) reproduce a biological scenario in which the spatial region considered ispart of a larger area of tissue whereby similar dynamics of the cells and the ECM occur. In this section, we carry out LSA of a biologically relevant homogeneous steady state of the system of PDEs (10),(11) and (16) (see Section 4.1) and we compare the dispersion relations obtained when the constitutive mod-els (1)-(6) are alternatively used to represent the contribution to the overall stress coming from the ECM, inorder to explore the pattern formation potential of these stress-strain constitutive equations (see Section 4.2).
Biologically relevant homogeneous steady state.
All non-trivial homogeneous steady states (¯ n, ¯ ρ, ¯ u ) (cid:124) of the system of PDEs (10), (11) and (16) subject to boundary conditions (17) have components ¯ n ≡ u ≡
0, and we choose ¯ ρ ≡ ρ . Hence, we consider the biologically relevant homogeneoussteady state ¯v = (1 , , (cid:124) . Linear stability analysis to spatially homogeneous perturbations.
In order to undertake linear sta-bility analysis of the steady state ¯v = (1 , , (cid:124) to spatially homogeneous perturbations, we make the ansatz v ( t, x ) ≡ ¯v + ˜v ( t ), where the vector ˜v ( t ) = (˜ n ( t ) , ˜ ρ ( t ) , ˜ u ( t )) (cid:124) models small spatially homogeneous perturba-tions, and linearise the system of PDEs (10), (11) and (16) about the steady state ¯v . Assuming ˜ n ( t ), ˜ ρ ( t ) and˜ u ( t ) to be proportional to exp ( ψt ), with ψ (cid:54) = 0, one can easily verify that ψ satisfies the algebraic equation ψ ( ψ + r )( ψ a + ψa + a ) = 0. Since r is positive and the parameters a , a and a are all non-negative, thesolution ψ of such an algebraic equation is necessarily negative and, therefore, the small perturbations ˜ n ( t ), ˜ ρ ( t )and ˜ u ( t ) will decay to zero as t → ∞ . This implies that the steady state ¯v will be stable to spatially homo-geneous perturbations for any choice of the parameter a , a , a , b , b and b in the stress-strain constitutiveequation (14) ( i.e. for all constitutive models (1)-(6)). 8 inear stability analysis to spatially inhomogeneous perturbations. In order to undertake linearstability analysis of the steady state ¯v = (1 , , (cid:124) to spatially inhomogeneous perturbations, we make the ansatz v ( t, x ) = ¯v + ˜v ( t, x ), where the vector ˜v ( t, x ) = (˜ n ( t, x ) , ˜ ρ ( t, x ) , ˜ u ( t, x )) (cid:124) models small spatially inhomogeneousperturbations, and linearise the system of PDEs (10), (11) and (16) about the steady state ¯v . Assuming ˜ n ( t, x ),˜ ρ ( t, x ) and ˜ u ( t, x ) to be proportional to exp ( ψt + ikx ), with ψ (cid:54) = 0 and k (cid:54) = 0, we find that ψ satisfies thefollowing equation c ( k ) ψ + c ( k ) ψ + c ( k ) ψ + c ( k ) = 0 , (18)with c ( k ) := ( b − τ ( λ + λ ) a ) k + a s (19) c ( k ) := D ( b − τ λ a ) k + (cid:0) b + a Ds + b r − τ ( λ + λ ) a − τ ( λ r + λ α ) a (cid:1) k + s ( a + a r ) (20) c ( k ) := D ( b − τ λ a ) k + (cid:0) b + a Ds + b r − τ ( λ + λ ) a − τ ( λ r + λ α ) a (cid:1) k + s ( a + a r ) (21)and c ( k ) := D ( b − τ λ a ) k + (cid:0) a Ds + b r − τ ( λ r + λ α ) a (cid:1) k + sa r , (22)where λ := 11 + λ and λ := (1 − λ )(1 + λ ) . Solving equation (18) for ψ gives the dispersion relation Re( ψ ( k )), where Re( · ) denotes the real part of( · ). For cell patterns to emerge, we need the non-trivial homogeneous steady state ¯v to be driven unstable byspatially inhomogeneous perturbations, that is, we need Re( ψ ( k )) > k >
0. Notice that a necessarycondition for this to happen is that at least one amongst c ( k ), c ( k ), c ( k ) and c ( k ) is negative for some k >
0. Hence, the fact that if τ = 0 then c ( k ), c ( k ), c ( k ) and c ( k ) are all non-negative for any valueof k allows us to conclude that having τ > c ( k ) = 0 and c ( k ) = 0, solving equation (18) for ψ gives the following dispersion relation ψ ( k ) = − c ( k ) c ( k ) (23)and for the condition Re( ψ ( k )) > c ( k ) > c ( k ) < c ( k ) < c ( k ) > . On the other hand, when the model parameters are such that only c ( k ) = 0, from equation (18) we obtainthe following dispersion relation ψ ( k ) = − c ( k ) ± (cid:113)(cid:0) c ( k ) (cid:1) − c ( k ) c ( k )2 c ( k ) , (24)and for the condition Re( ψ ( k )) > c ( k ) > c ( k ) < c ( k ) > , c ( k ) < c ( k ) > c ( k ) < c ( k ) > c ( k ) < , c ( k ) > c ( k ) < . c ( k ) (cid:54) = 0 as well, from equation (18)we obtain the following dispersion relation ψ ( k ) = (cid:26) q + (cid:104) q + (cid:0) m − p (cid:1) (cid:105) / (cid:27) / + (cid:26) q − (cid:104) q + (cid:0) m − p (cid:1) (cid:105) / (cid:27) / + p , (25)where p ≡ p ( k ), q ≡ q ( k ) and m ≡ m ( k ) are defined as p := − c c , q := p + c c − c c c , m := c c . In this case, identifying sufficient conditions to ensure that the real part of ψ ( k ) is positive for some k > Substituting the definitions of a , a , a , b , b and b corresponding to the stress-strain constitutive equa-tions (1)-(6), which are reported in Table 1, into definitions (19)-(22) for c ( k ), c ( k ), c ( k ) and c ( k ), andthen using the dispersion relation (23), (24) or (25) depending on the values of c ( k ) and c ( k ) so obtained,we derive the dispersion relation for each of the constitutive models (1)-(6). Base-case dispersion relations.
Figure 1 displays the dispersion relations obtained for the stress-strainconstitutive equations (1)-(6) under the following base-case parameter values E = E = E = 0 . , η = η = η = 0 . , D = 0 . , (26) α = 0 . , r = 0 . , s = 0 . , λ = 0 . , τ = 0 . . (27)The parameter values given by (26) and (27) are chosen for illustrative purposes only in order to highlightthe different qualitative behaviour of the dispersion relations obtained using different models. A comparisonbetween the plots in Figure 1 reveals that constitutive models capturing viscous flow, that is, the linear viscousmodel (2), the Maxwell model (4) and the 3-parameter viscous model (6) ( cf. Table 2), have a much higherpattern formation potential than the others, since they have a much broader range of values of k for whichRe( ψ ( k )) > i.e. a larger set of unstable modes).We now undertake a sensitivity analysis with respect to the different model parameters and discuss keychanges that occur in the base-case dispersion relations displayed in Figure 1. ECM elasticity.
The plots in Figure 2 illustrate how the base-case dispersion relations displayed in Figure 1change when different values of the parameters E , E and E ( i.e. the parameters modelling ECM elasticity)are considered. These plots show that larger values of these parameters: • correlate with a narrower range of values of k for which Re( ψ ( k )) > • leave the unstable modes unchanged but reduce the values of Re( ψ ( k )) for the 3-parameter viscousmodel (6), which corresponds to slowing down the formation of spatial patterns.Sufficiently small values of the parameters E , E and E appear to lead to singular dispersion relations ( cf. theplots for the linear elastic model (1), the Maxwell model (4) and the SLS model (5) in Figure 2), which suggeststhat linear stability theory may fail in the regime of low ECM elasticity.10igure 1: Base-case dispersion relations.
Dispersion relations corresponding to the stress-strain constitutiveequations (1)-(6) for the base-case set of parameter values given by (26) and (27).
ECM viscosity.
The plots in Figure 3 illustrate how the base-case dispersion relations displayed in Figure 1change when different values of the parameters η , η and η ( i.e. the parameters modelling ECM viscosity)are considered. These plots show that larger values of these parameters leave the range of values of k forwhich Re( ψ ( k )) > ψ ( k )). This supports the idea that a higher ECMviscosity may not change the pattern formation potential of the different constitutive models but may slowdown the corresponding pattern formation processes. Cell motility.
The plots in Figure 4 illustrate how the base-case dispersion relations displayed in Figure 1change when different values of the parameter D ( i.e. the parameter modelling cell motility) are considered.These plots show that larger values of this parameter may significantly shrink the range of values of k for whichRe( ψ ( k )) >
0. In particular, the linear viscous model (2), the Kelvin-Voigt model (3), and the SLS model (5)display – as D increases – infinitely many, a finite number of and no unstable modes. This is to be expecteddue to the stabilising effect of undirected, random cell movement and indicates that higher cell motility maycorrespond to lower pattern formation potential. Intrinsic growth rate of the cell density.
The plots in Figure 5 illustrate how the base-case dispersionrelations displayed in Figure 1 change when different values of the parameter r ( i.e. the intrinsic growth rateof the cell density) are considered. These plots show that considering larger values of this parameter reducesthe values of Re( ψ ( k )) for the linear viscous model (2), the Maxwell model (4) and the 3-parameter viscousmodel (6), whereas it shrinks the range of unstable modes for the other constitutive models. This supportsthe idea that higher growth rates of the cell density ( i.e. faster dynamics of cell proliferation and death) mayslow down pattern formation processes for constitutive models that capture viscous flow and reduce the pattern11ormation potential of the other constitutive models. Elasticity of the surrounding tissue and level of contact inhibition of the cell traction forces.
The plots in Figures 6 and 7 illustrate how the base-case dispersion relations displayed in Figure 1 changewhen different values of the parameter s ( i.e. the elasticity of the surrounding tissue) and the parameter λ ( i.e. the level of cell-cell contact inhibition of the cell traction forces) are considered. Considerations similar tothose previously made about the dispersion relations obtained for increasing values of the parameter r applyto the case where increasing values of the parameter s and the parameter λ are considered. In addition to suchconsiderations, the plots in Figure 6 indicate that higher values of s may also reduce the pattern formationpotential of the different constitutive models by making more likely that Re( ψ ( k )) < k ( i.e. low-frequency perturbation modes will be more likely to vanish). Moreover, the plots in Figure 7 suggestthat singularities in the dispersion relations may develop for λ = 0 ( cf. the plots for the linear elastic model (1),the Maxwell model (4) and the SLS model (5)), which suggests that linear stability theory may fail in theabsence of contact inhibition of cell traction forces due to cell crowding. Cell haptotactic sensitivity and cell traction forces.
The plots in Figures 8 and 9 illustrate how thebase-case dispersion relations displayed in Figure 1 change when different values of the parameter α ( i.e. thecell haptotactic sensitivity) and the parameter τ ( i.e. the cell traction force) are, respectively, considered. Asexpected [40], larger values of these parameters broad the range of values of k for which Re( ψ ( k )) > Effects of varying the ECM elasticity.
Dispersion relations corresponding to the stress-strainconstitutive equations (1)-(6) for increasing values of the ECM elasticity, that is, E = E = E = 0 .
01 (solidline), E = E = E = 0 . E = E = E = 1 (dash-dot lines). The values of the otherparameters are given by (26) and (27). 12igure 3: Effects of varying the ECM viscosity.
Dispersion relations corresponding to the stress-strainconstitutive equations (1)-(6) for increasing values of the ECM viscosity, that is, η = η = η = 0 .
01 (solid line), η = η = η = 0 . η = η = η = 10 (dash-dot lines). The values of the other parametersare given by (26) and (27).Figure 4: Effects of varying the cell motility.
Dispersion relations corresponding to the stress-strainconstitutive equations (1)-(6) for increasing values of the cell motility, that is, D = 0 .
001 (solid line), D = 0 . D = 1 (dash-dot lines). The values of the other parameters are given by (26) and (27).13igure 5: Effects of varying the intrinsic growth rate of the cell density.
Dispersion relations corre-sponding to the stress-strain constitutive equations (1)-(6) for increasing values of the intrinsic growth rate ofthe cell density, that is, r = 0 (solid line), r = 0 . r = 1 (dash-dot lines). The values of theother parameters are given by (26) and (27).Figure 6: Effects of varying the elasticity of the surrounding tissue.
Dispersion relations correspondingto the stress-strain constitutive equations (1)-(6) for increasing values of the elasticity of the surrounding tissue,that is, s = 0 (solid line), s = 0 . s = 2 (dash-dot lines). The values of the other parametersare given by (26) and (27). 14igure 7: Effects of varying the level of cell-cell contact inhibition of the cell traction forces.
Dispersion relations corresponding to the stress-strain constitutive equations (1)-(6) for increasing levels of cell-cell contact inhibition of the cell traction forces, that is, λ = 0 (solid line), λ = 0 . λ = 1(dash-dot lines). The values of the other parameters are given by (26) and (27).Figure 8: Effects of varying the cell haptotactic sensitivity.
Dispersion relations corresponding to thestress-strain constitutive equations (1)-(6) for increasing values of the cell haptotactic sensitivity, that is, α =0 .
001 (solid line), α = 0 . α = 2 (dash-dot lines). The values of the other parameters aregiven by (26) and (27). 15igure 9: Effects of varying the cell traction forces.
Dispersion relations corresponding to the stress-strainconstitutive equations (1)-(6) for increasing cell traction forces, that is, τ = 10 − (solid line), τ = 0 .
05 (dashedlines) and τ = 0 . In this section, we verify key results of LSA presented in Section 4 by solving numerically the system ofPDEs (10), (11) and (16) subject to boundary conditions (17). In particular, we report on numerical solutionsobtained in the case where equation (16) is complemented with the Kelvin-Voigt model (3) and the Maxwellmodel (4). A detailed description of the numerical schemes employed is provided in Appendix A.
Set-up of numerical simulations.
We carry out numerical simulations using the parameter values givenby (26) and (27). We choose the endpoints of the spatial domain to be (cid:96) = − L = 2, and the final time T is chosen sufficiently large so that distinct spatial patterns can be observed at the end of simulations. Weconsider first the initial conditions n (0 , x ) = 1 + 0 .
01 sin(10 πx ) , ρ (0 , x ) ≡ , u (0 , x ) ≡ n (0 , x ) = 1 + 0 . (cid:15) ( x ) , ρ (0 , x ) ≡ , u (0 , x ) ≡ , (29)where (cid:15) ( x ) is a normally distributed random variable with mean 0 and variance 1 for each x ∈ [ − , n = 1, ρ = 1and u = 0. This is the steady state considered in the LSA undertaken in Section 4.1. Numerical computationsare performed in MATLAB. 16 ain results. The results obtained are summarised by the plots in Figures 10-13. These results demonstratethat, in agreement with the dispersion relations displayed in Figure 1, for the parameter values given by (26)and (27), small perturbations present in the initial cell density, being them periodic (see Figures 10 and 11) orrandomly distributed (see Figures 12 and 13):- vanish in the case of the Kelvin-Voigt model, thus leading the cell density to relax to the homogeneoussteady state n = 1 and attain numerical equilibrium at t = 10 while leaving the ECM density unchanged;- grow in the case of the Maxwell model, leading to the formation of spatial patterns both in the cell density n and in the ECM density ρ .Notice that the formation of spatial patterns correlates with the growth of the cell-ECM displacement u . Infact, the displacement remains close to zero ( i.e. ∼ O (10 − )) for the Kelvin-Voigt model, whereas it grows withtime for the Maxwell model. Due to the presence of infinite unstable modes ( cf. the dispersion relation obtainedfor the Maxwell model in Figure 1), the spatial patterns observed reflect the nature of the small perturbationsthat are superimposed to the homogeneous steady state. In fact, periodic spatial patterns are formed in the caseof periodic small perturbations (see Figure 11), while random small perturbations bring about the formation ofdisordered spatial patterns (see Figure 13). Finally, we remark that the spatial patterns displayed in Figures 11and 13 are not stable as the displacement u ( t, x ) continues to grow and blows up for t → ∞ .Figure 10: Simulation results for the Kelvin-Voigt model (3) under initial conditions (28) . Celldensity n ( t, x ) (left), ECM density ρ ( t, x ) (centre) and cell-ECM displacement u ( t, x ) (right) at t = 0 (top)and t = 10 (bottom) obtained solving numerically the system of PDEs (10), (11) and (16) complemented withthe Kelvin-Voigt model (3), subject to boundary conditions (17) and initial conditions (28), for the parametervalues given by (26) and (27). 17igure 11: Simulation results for the Maxwell model (4) under initial conditions (28) . Cell density n ( t, x ) (left), ECM density ρ ( t, x ) (centre) and cell-ECM displacement u ( t, x ) (right) at t = 0 (top) and t = 40(bottom) obtained solving numerically the system of PDEs (10), (11) and (16) complemented with the Maxwellmodel (4), subject to boundary conditions (17) and initial conditions (28), for the parameter values given by (26)and (27).Figure 12: Simulation results for the Kelvin-Voigt model (3) under initial conditions (29) . Celldensity n ( t, x ) (left), ECM density ρ ( t, x ) (centre) and cell-ECM displacement u ( t, x ) (right) at t = 0 (top)and t = 10 (bottom) obtained solving numerically the system of PDEs (10), (11) and (16) complemented withthe Kelvin-Voigt model (3), subject to boundary conditions (17) and initial conditions (29), for the parametervalues given by (26) and (27). 18igure 13: Simulation results for the Maxwell model (4) under initial conditions (29) . Cell density n ( t, x ) (left), ECM density ρ ( t, x ) (centre) and cell-ECM displacement u ( t, x ) (right) at t = 0 (top) and t = 38(bottom) obtained solving numerically the system of PDEs (10), (11) and (16) complemented with the Maxwellmodel (4), subject to boundary conditions (17) and initial conditions (29), for the parameter values given by (26)and (27). In this section, we complement the results presented in the previous sections with the results of numericalsimulations of a two-dimensional mechanical model of pattern formation in biological tissues. In particular,we report on numerical solutions obtained in the case where the two-dimensional analogue of the system ofPDEs (10), (11) and (16) is furnished with a two-dimensional version of the one-dimensional Maxwell model (4).
A two-dimensional mechanical model of pattern formation.
The mechanical model of pattern forma-tion defined by the system of PDEs (10), (11) and (16) posed on a two-dimensional spatial domain representedby a bounded set Ω ⊂ R with smooth boundary ∂ Ω reads as ∂ t n = ∇ x · [ n ( D ∇ x n − α ∇ x ρ − ∂ t u )] + r n (1 − n ) ,∂ t ρ = −∇ x · ( ρ ∂ t u ) , ∇ x · ( σ m + σ c ) + ρ F = 0 , (30)with t ∈ (0 , T ] and x = ( x , x ) (cid:124) ∈ Ω. We close the system of PDEs (30) imposing the two-dimensional versionof the periodic boundary conditions (17) on ∂ Ω. Furthermore, we use the following two-dimensional analoguesof definitions (13) and (15) σ c := τ nρ λn I and F := − s u , (31)where I is the identity tensor. Moreover, in analogy with the one-dimensional case, we define the stress tensor σ m via the two-dimensional constitutive model that is used to represent the stress-strain relation of the ECM.19n particular, we consider the following two-dimensional version of the one-dimensional Maxwell model (4),which is derived in Appendix B, 1 η σ + 1 E (cid:48) ∂ t σ = ∂ t ε + ν (cid:48) ∂ t θ I . (32)Here, the strain ε ( t, x ) and the dilation θ ( t, x ) are defined in terms of the displacement u ( t, x ) as ε = 12 (cid:0) ∇ u + ∇ u (cid:124) (cid:1) and θ = ∇ · u . (33)Notice that both ε and θ reduce to ε = ∂ x u in the one-dimensional case. In the stress-strain constitutiveequation (32), η is the shear viscosity, E (cid:48) := E ν and ν (cid:48) := ν − ν , (34)where ν is Poisson’s ratio and E is Young’s modulus. As clarified in Appendix B, the two-dimensional Maxwellmodel (32) holds under the simplifying assumption that the quotient between the shear viscosity and the bulkviscosity of the ECM is equal to ν (cid:48) . Set-up of numerical simulations.
We solve numerically the system of PDEs (30) subject to the two-dimensional version of the periodic boundary conditions (17) and complemented with (31)-(34). Numericalsimulations are carried out using the following parameter values E = 0 . , η = 0 . , D = 0 . , ν = 0 . , (35) α = 0 . , r = 0 . , s = 0 . , λ = 0 . , τ = 0 . , (36)which are chosen for illustrative purposes only. We choose Ω = [0 , × [0 ,
1] and the final time T is chosensufficiently large so that distinct spatial patterns can be observed at the end of simulations. We consider firstthe following two-dimensional analogue of initial conditions (28) n (0 , x , x ) = 1 + 0 .
01 sin(20 πx ) sin(20 πx ) , ρ (0 , x , x ) ≡ , u (0 , x , x ) ≡ (37)and then the following two-dimensional analogue of initial conditions (29) n (0 , x , x ) = 1 + 0 . (cid:15) ( x , x ) , ρ (0 , x , x ) ≡ , u (0 , x , x ) ≡ , (38)where (cid:15) ( x , x ) is a normally distributed random variable with mean 0 and variance 1 for each ( x , x ) ∈ [0 , × [0 , Main results.
The results obtained are summarised by the plots in Figures 14 and 15, which demonstrate that,for the parameter values given by (35) and (36), both periodic and randomly distributed small perturbations ofthe initial cell density grow with time resulting in the formation of cellular patterns. Analogous spatial patternsbut with a smaller gradient are observed in the ECM density (results not shown). The black arrows in Figures 14and 15 highlight the cell-ECM displacement. Similar to the case of the one-dimensional patterns displayed inFigures 11 and 13, these patterns reflect the nature of the small perturbations that are superimposed on thehomogeneous steady state and are not stable as the displacement continues to grow and ultimately blows up.20igure 14:
Simulation results for the two-dimensional Maxwell model (32) under initial condi-tions (37) . Cell density n ( t, x , x ) at t ≈
17 obtained solving numerically the system of PDEs (30) subject tothe two-dimensional version of the periodic boundary conditions (17) and initial conditions (37), complementedwith (31)-(34), for the parameter values given by (35) and (36). The black arrows highlight the correspondingcell-ECM displacement u ( t, x , x ) at t ≈ Conclusions.
We have investigated the pattern formation potential of different stress-strain constitutiveequations for the ECM within a one-dimensional mechanical model of pattern formation in biological tissuesformulated as the system of PDEs (10), (11) and (16). The results of linear stability analysis of a biologicallyrelevant homogeneous steady state of the model equations presented in Section 4 and the dispersion relationsderived therein support the idea that stress-strain constitutive equations capturing viscous flow and permanentset ( i.e. the linear viscous model (2), the Maxwell model (4) and the 3-parameter viscous model (6)) have apattern formation potential much higher than the others ( i.e. the linear elastic model (1), the Kelvin-Voigtmodel (3) and the SLS model (5)). This is confirmed by the results of numerical simulations presented inSection 5, which demonstrate that, all else being equal, spatial patterns emerge in the case where the Maxwellmodel (4) is used to represent the stress-strain relation of the ECM, while no patterns are observed when theKelvin-Voigt model (3) is employed. In Section 6, as an illustrative example, we have also reported on theresults of numerical simulations of a two-dimensional version of the model, which is given by the system ofPDEs (30) complemented with the two-dimensional Maxwell model (32). These results confirm the capabilityof such a model of linear viscoelasticity of producing cellular patterns.Our findings corroborate the conclusions of [9] suggesting that prior studies on mechanochemical models of21igure 15:
Simulation results for the two-dimensional Maxwell model (32) under initial condi-tions (38) . Cell density n ( t, x , x ) at t ≈
21 obtained solving numerically the system of PDEs (30) subject tothe two-dimensional version of the periodic boundary conditions (17) and initial conditions (38), complementedwith (31)-(34), for the parameter values given by (35) and (36). The black arrows highlight the correspondingcell-ECM displacement u ( t, x , x ) at t ≈ cf. [15]. Research perspectives.
The dispersion relations given in Section 4 indicate that, when pattern formationoccurs, constitutive models of linear viscoelasticity that capture viscous flow and permanent set display infinitelymany unstable modes. Hence, it is natural to expect the spatial patterns that are formed to reflect the choiceof the small perturbations which are superimposed to the homogeneous steady state. This is confirmed bythe numerical solutions of the Maxwell model (4). Moreover, as has been remarked in Section 5, the spatialpatterns we have reported on in this paper are not stable, which is due to the fact that the cell-ECM displacementcontinues to grow and ultimately blows up. In this regard, it would be interesting to consider extended versionsof the mechanical model of pattern formation given by the system of PDEs (10), (11) and (16) that overcomethese limitations. This would enable one to explore possible differences in the spatial patterns obtained whendifferent stress-strain constitutive equations for the ECM are used – such as amplitude of patterns, perturbationmode selection and geometric structure in two spatial dimensions. For instance, it is known that including long-22ange effects can promote the formation of stable spatial patterns [38, 50], which could be explored throughnonlinear stability analysis, as previously done for the case in which the stress-strain relation of the ECM isrepresented by the Kelvin-Voigt model [11, 26, 30].It would also be interesting to construct numerical solutions for the mechanical model defined by the systemof PDEs (10), (11) and (16) furnished with the 3-parameter viscous model (6). For this to be done, suitableextensions of the numerical schemes presented in Appendix A need to be developed.It would also be relevant to systematically assess the pattern formation potential of different constitutive mod-els of viscoelasticity in two spatial dimensions. This would require to relax the simplifying assumption (A.9)on the shear and bulk viscosities of the ECM, which we have used to derive the two-dimensional Maxwellmodel (32), and, more in general, to find analytically and computationally tractable stress-strain-dilation rela-tions, which still remains an open problem [6, 18]. In order to solve this problem, new methods of derivationand parameterisation for constitutive models of viscoelasticity might need to be developed [59].As previously mentioned, the values of the model parameters used in this paper have been chosen for illustra-tive purposes only. Hence, it would be useful to re-compute the dispersion relations and the numerical solutionspresented here for a calibrated version of the model based on real biological data. On a related note, an inter-esting application that could be explored by varying parameter values in the generic constitutive equation (16)is cancer-associated fibrosis. This is a disease characterised by an excessive production of collagen, elastin andproteoglycans, which directly affects the structure of the ECM resulting in alterations of viscoelastic tissueproperties [12]. Such alterations in the ECM may facilitate tumour invasion and angiogenesis. Considering acalibrated mechanical model of pattern formation in biological tissues, whereby the values of the parameters inthe stress-strain constitutive equation for the ECM change during fibrosis progression, may shed new light onthe existing connections between structural changes in the ECM components and higher levels of malignancyin cancer [10, 51].
Conflict of interest
The authors declare that they have no conflict of interest.
A Numerical schemes for the system of PDEs (10) , (11) and (16) Numerical solutions for the system of PDEs (10), (11) and (16) are constructed using a uniform discretisationof the spatial domain consisting of N points and a discretisation of the time domain with a variable time step∆ t k determined by a CFL condition. The normalised cell density n ( t, x ), the normalised ECM density ρ ( t, x )and the displacement of a material point of the cell-ECM system u ( t, x ) are approximated as n ( t k , x i ) ≈ n ki , ρ ( t k , x i ) ≈ ρ ki , u ( t k , x i ) ≈ u ki for i = 0 , . . . , N − k = 1 , , . . . . Numerical scheme for the balance equation (10) . Using the notation v = ∂ t u , we rewrite the balanceequation (10) as ∂ t n = D∂ xx n − ∂ x ( φ n ) + rn (1 − n ) with φ = α ∂ x ρ + v , which we discretise as n k +1 i = n ki + ∆ t k +1 (cid:34) D ∆ x (cid:16) n ki +1 − n ki + n ki − (cid:17) − x (cid:16) F ki + − F ki − (cid:17) + r n ki (cid:16) − n ki (cid:17)(cid:35) (A.1)with φ ki = α ∆ x (cid:16) ρ ki +1 − ρ ki (cid:17) + 1∆ t k +1 (cid:16) u k +1 i − u ki (cid:17) . (A.2)23he value of u k +1 i in (A.2) is obtained using a predictor-corrector method as similarly done in [34]. Thenumerical fluxes F ki + and F ki − in (A.1) are defined as F ki + ≡ F (cid:16) φ ki + , n ki + (cid:17) := max (cid:16) , φ ki + (cid:17) n kLi + + min (cid:16) , φ ki + (cid:17) n kRi + ,F ki − ≡ F (cid:16) φ ki − , n ki − (cid:17) := max (cid:16) , φ ki − (cid:17) n kLi − + min (cid:16) , φ ki − (cid:17) n kRi − , with the extrapolated cell-edge variables being defined as φ ki + := 14 (cid:16) φ kRi + + φ kLi + (cid:17) , φ ki − := 14 (cid:16) φ kRi − + φ kLi − (cid:17) , with φ kRi + := φ ki +1 − ϕ ( f ki +1 ) (cid:0) φ ki +2 − φ ki +1 (cid:1) , φ kRi − := φ ki − ϕ ( f ki ) (cid:0) φ ki +1 − φ ki (cid:1) ,φ kLi + := φ ki + 12 ϕ ( f ki +1 ) (cid:0) φ ki +1 − φ ki (cid:1) , φ kLi − := φ ki − + 12 ϕ ( f ki − ) (cid:0) φ ki − φ ki − (cid:1) , and n kRi + := n ki +1 − ϕ ( g ki +1 ) (cid:0) n ki +2 − n ki +1 (cid:1) , n kRi − := n ki − ϕ ( g ki ) (cid:0) n ki +1 − n ki (cid:1) ,n kLi + := n ki + 12 ϕ ( g ki +1 ) (cid:0) n ki +1 − n ki (cid:1) , n kLi − := n ki − + 12 ϕ ( g ki − ) (cid:0) n ki − n ki − (cid:1) , where f ki := φ ki − φ ki − φ ki +1 − φ ki , g ki := n ki − n ki − n ki +1 − n ki and ϕ ( · ) is the OSPRE flux limiter function [25]. Numerical scheme for the transport equation (11) . Using the notation v = ∂ t u , we rewrite the transportequation (11) as ∂ t ρ = − ∂ x ( v ρ )which we discretise as ρ k +1 i = ρ ki − ∆ t k +1 ∆ x (cid:16) G ki + − G ki − (cid:17) (A.3)with v ki = 1∆ t k +1 (cid:16) u k +1 i − u ki (cid:17) . (A.4)The numerical fluxes G ki + and G ki − in (A.1) are defined as G ki + ≡ G (cid:16) v ki + , ρ ki + (cid:17) := max (cid:16) , v ki + (cid:17) ρ kLi + + min (cid:16) , v ki + (cid:17) ρ kRi + ,G ki − ≡ G (cid:16) v ki − , ρ ki − (cid:17) := max (cid:16) , v ki − (cid:17) ρ kLi − + min (cid:16) , v ki − (cid:17) ρ kRi − , v ki + := 14 (cid:16) v kRi + + v kLi + (cid:17) , v ki − := 14 (cid:16) v kRi − + v kLi − (cid:17) , with v kRi + := v ki +1 − ϕ ( h ki +1 ) (cid:0) v ki +2 − v ki +1 (cid:1) , v kRi − := v ki − ϕ ( h ki ) (cid:0) v ki +1 − v ki (cid:1) ,v kLi + := v ki + 12 ϕ ( h ki +1 ) (cid:0) v ki +1 − v ki (cid:1) , v kLi − := v ki − + 12 ϕ ( h ki − ) (cid:0) v ki − v ki − (cid:1) , and ρ kRi + := ρ ki +1 − ϕ ( m ki +1 ) (cid:0) ρ ki +2 − ρ ki +1 (cid:1) , n kRi − := ρ ki − ϕ ( m ki ) (cid:0) ρ ki +1 − ρ ki (cid:1) ,ρ kLi + := ρ ki + 12 ϕ ( m ki +1 ) (cid:0) ρ ki +1 − ρ ki (cid:1) , n kLi − := ρ ki − + 12 ϕ ( m ki − ) (cid:0) ρ ki − ρ ki − (cid:1) , where h ki := v ki − v ki − v ki +1 − v ki and m ki := ρ ki − ρ ki − ρ ki +1 − ρ ki . As mentioned earlier, the value of u k +1 i in (A.4) is obtained usinga predictor-corrector method as similarly done in [34], and ϕ ( · ) is the OSPRE flux limiter function [25]. Numerical scheme for the force-balance equation (16) . The most popular approach adopted in theliterature to solve numerically equation (16) is to use NAG library routines, such as the NAG solver D03PHF[47] or the NAG Mark 15 routine D03EBF [34]. Here, we employ the numerical scheme described hereafter,which can be used to construct numerical solutions for (16) in the case where b = a = 0 and b (cid:54) = 0 – i.e. whenthe linear viscous model (2), Kelvin-Voigt model (3), Maxwell model (4) and the SLS model (5) are considered.This scheme can be easily adapted to the case where b = a = b = a = 0 ( i.e. when equation (16) iscomplemented with the linear elastic model (1)), while solving numerically (16) in the case where b (cid:54) = 0 ( i.e. when the 3-parameter viscous model (6) is used) would require a different numerical scheme.When b = a = 0 and b (cid:54) = 0, the force-balance equation (16) reads as b ∂ t ∂ xx u + b ∂ xx u = a ∂ t ( sρu ) + a sρu − (cid:34) a ∂ t ∂ x (cid:16) τ nρ λn (cid:17) + a ∂ x (cid:16) τ nρ λn (cid:17)(cid:35) , which can then be rewritten as the following system of PDEs ∂ t w = a (cid:32) sρu − ∂ x (cid:16) τ nρ λn (cid:17)(cid:33) − b ∂ xx u − a ∂ t ∂ x (cid:16) τ nρ λn (cid:17) ,b ∂ xx u − a sρu = w . (A.5)We solve the parabolic PDE (A.5) for w ( t, x ) using the following numerical scheme w k +1 i = w ki + ∆ t k +1 (cid:34) a (cid:16) sρ ki u ki − H ki (cid:17) − b ∆ x (cid:0) u ki +1 − u ki + u ki − (cid:1) − a ∆ t k +1 (cid:0) H k +1 i − H ki +1 (cid:1)(cid:35) H ki is the discretisation of the term ∂ x (cid:16) τ nρ λn (cid:17) , that is, H ki := (cid:34) n ki x (cid:0) ρ ki +1 − ρ ki (cid:1) + (cid:0) − λ ( n ki ) (cid:1)(cid:0) λ ( n ki ) (cid:1) ρ ki x (cid:0) n ki +1 − n ki (cid:1)(cid:35)(cid:0) λ ( n ki ) (cid:1) . Moreover, we solve the elliptic PDE (A.5) for u ( t, x ) by iteration using the Jacobi method whereby given u atiteration j we compute u at iteration j + 1 as u j +1 i = 12 (cid:34) u ji +1 + u ji − − a ∆ x sb ρ k +1 i u ji − ∆ x b w k +1 i (cid:35) , (A.6)and we iterate until convergence. This scheme has proven stable and convergent for the parameter sets consideredhere. Notice that in the case where b = a = b = a = 0 ( i.e. when the linear elastic model (1) is considered),the force-balance equation (16) reduces to an elliptic equation for u ( t, x ), which can be directly solved withoutintroducing the dependent variable w ( t, x ). On the other hand, in the case where b (cid:54) = 0 ( i.e. when the 3-parameter viscous model (6) is considered) the above numerical scheme cannot be used due to the presence ofa second order derivative in t . B Derivation of the two-dimensional Maxwell model (32)
Landau & Lifshitz derived from first principles the stress-strain relations that give the two-dimensional versionsof the linear elastic model (1) and of the linear viscous model (2) in isotropic materials [24], which read,respectively, as σ e = E ν (cid:16) ε e + ν − ν θ e I (cid:17) and σ v = η ∂ t ε v + µ ∂ t θ v I . (A.7)Here, E is Young’s modulus, ν is Poisson’s ratio, I is the identity tensor, η is the shear viscosity and µ is thebulk viscosity. Moreover, ε e and θ e are the strain and dilation under a purely elastic deformation u e while ε v and θ v are the strain and dilation under a purely viscous deformation u v , which are all defined via (33).In the case of a linearly viscoelastic material satisfying Kelvin-Voigt model, the two dimensional analogueof (3) is simply given by σ = σ e + σ v . This is the stress-strain constitutive equation that is typically used to describe the contribution to the stress ofthe cell-ECM system coming from the ECM in two-dimensional mechanochemical models of pattern formation[11, 13, 22, 30, 34, 40, 43, 44, 45, 46, 49, 50, 52].On the other hand, deriving the two-dimensional analogues of Maxwell model (4), of the SLS model (5) andof the 3-parameter viscous model (6) is more complicated due to the presence of elements connected in series.In the case of Maxwell model, using the fact that the overall strain and dilation will be distributed over thedifferent components ( i.e. ε = ε e + ε v and θ = θ e + θ v ) along with the fact that the stress on each componentwill be the same as the overall stress ( i.e. σ = σ e = σ v ), one finds1 η σ + 1 E (cid:48) ∂ t σ = ∂ t ε + ν (cid:48) ∂ t θ I + (cid:18) ν (cid:48) − µη (cid:19) ∂ t θ v I , (A.8)with E (cid:48) and ν (cid:48) being defined via (34). Under the simplifying assumption that µη = ν (cid:48) (A.9)the stress-strain constitutive equation (A.8) reduces to the two-dimensional Maxwell model (32).26 eferences [1] S. Alonso, M. Radszuweit, H. Engel, and M. B¨ar , Mechanochemical pattern formation in simplemodels of active viscoelastic fluids and solids , Journal of Physics D: Applied Physics, 50 (2017), p. 434004.[2]
V. H. Barocas, A. G. Moon, and R. T. Tranquillo , The fibroblast-populated collagen microsphereassay of cell traction force—part 2: Measurement of the cell traction parameter , Journal of BiomechanicalEngineering, 117 (1995), pp. 161–170.[3]
V. H. Barocas and R. T. Tranquillo , Biphasic theory and in vitro assays of cell-fibril mechanicalinteractions in tissue-equivalent gels , in Cell Mechanics and Cellular Engineering, Springer, 1994, pp. 185–209.[4]
D. E. Bentil and J. D. Murray , Pattern selection in biological pattern formation mechanisms , AppliedMathematics Letters, 4 (1991), pp. 1–5.[5]
L. E. Bilston, Z. Liu, and N. Phan-Thien , Linear viscoelastic properties of bovine brain tissue inshear , Biorheology, 34 (1997), pp. 377–385.[6]
V. B. Birman, W. K. Binienda, and G. Townsend ,
2d maxwell model , Journal of MacromolecularScience, Part B, 41 (2002), pp. 341–356.[7]
J. E. Bischoff, E. M. Arruda, and K. Grosh , A rheological network model for the continuumanisotropic and viscoelastic behavior of soft tissue , Biomechanics and Modeling in Mechanobiology, 3 (2004),pp. 56–65.[8]
F. Brinkmann, M. Mercker, T. Richter, and A. Marciniak-Czochra , Post-turing tissue patternformation: Advent of mechanochemistry , PLoS Computational Biology, 14 (2018), p. e1006259.[9]
H. M. Byrne and M. A. J. Chaplain , The importance of constitutive equations in mechanochemicalmodels of pattern formation , Applied Mathematics Letters, 9 (1996), pp. 85–90.[10]
C. Chandler, T. Liu, R. Buckanovich, and L. Coffman , The double edge sword of fibrosis in cancer ,Translational Research, (2019).[11]
G. C. Cruywagen and J. D. Murray , On a tissue interaction model for skin pattern formation , Journalof Nonlinear Science, 2 (1992), pp. 217–240.[12]
T. Ebihara, N. Venkatesan, R. Tanaka, and M. S. Ludwig , Changes in extracellular matrix andtissue viscoelasticity in bleomycin–induced lung fibrosis: Temporal aspects , American journal of respiratoryand critical care medicine, 162 (2000), pp. 1569–1576.[13]
I. Ferrenq, L. Tranqui, B. Vailhe, P. Y. Gumery, and P. Tracqui , Modelling biological gelcontraction by cells: mechanocellular formulation and cell traction force quantification , Acta Biotheoretica,45 (1997), pp. 267–293.[14]
W. N. Findley, J. S. Lai, and K. Onaran , Creep and relaxation of nonlinear viscoelastic materials -with an introduction to linear viscoelasticity , Dover Publications, 1976.[15]
Y. C. Fung , Biomechanics Mechanical Properties of Living Tissues (2nd Edition) , Springer-Verlag NewYork, 1993.[16]
A. Gierer and H. Meinhardt , A theory of biological pattern formation , Kybernetik, 12 (1972), pp. 30–39. 2717]
S. J. Gilmore, B. L. Vaughan Jr, A. Madzvamuse, and P. K. Maini , A mechanochemical model ofstriae distensae , Mathematical Biosciences, 240 (2012), pp. 141–147.[18]
M. Haghighi-Yazdi and P. Lee-Sullivan , Modeling linear viscoelasticity in glassy polymers using stan-dard rheological models , in Proceedings of the COMSOL Conference, Boston, MA, USA, 2011.[19]
A. K. Harris , Tissue culture cells on deformable substrata: Biomechanical implications , Journal of Biome-chanical Engineering, 106 (1984), pp. 19–24.[20]
A. K. Harris, D. Stopak, and P. Wild , Fibroblast traction as a mechanism for collagen morphogenesis ,Nature, 290 (1981), p. 249.[21]
Y.-P. Huang, Y.-P. Zheng, and S.-F. Leung , Quasi-linear viscoelastic properties of fibrotic neck tissuesobtained from ultrasound indentation tests in vivo , Clinical biomechanics, 20 (2005), pp. 145–154.[22]
E. Javierre, P. Moreo, M. Doblar´e, and J. M. Garc´ıa-Aznar , Numerical modeling of a mechano-chemical theory for wound contraction analysis , International Journal of Solids and Structures, 46 (2009),pp. 3597–3606.[23]
J. Jernvall, S. A. Newman, et al. , Mechanisms of pattern formation in development and evolution ,Development, 130 (2003), pp. 2027–2037.[24]
L. D. Landau and E. M. Lifshitz , Theory of Elasticity , Pergamon Press, 1970.[25]
R. J. LeVeque , Finite difference methods for ordinary and partial differential equations: steady-state andtime-dependent problems , SIAM, 2007.[26]
M. A. Lewis and J. D. Murray , Analysis of stable two-dimensional patterns in contractile cytogel ,Journal of Nonlinear Science, 1 (1991), pp. 289–311.[27]
Z. Liu and L. Bilston , On the viscoelastic character of liver tissue: experiments and modelling of thelinear behaviour , Biorheology, 37 (2000), pp. 191–201.[28]
P. Maini, K. Painter, and H. P. Chau , Spatial pattern formation in chemical and biological systems ,Journal of the Chemical Society, Faraday Transactions, 93 (1997), pp. 3601–3610.[29]
P. K. Maini , Morphogenesis, biological , in Morphogenesis, Biological, Encyclopedia of Nonlinear Science,A. Scott, ed., Routledge, 2005, pp. 587–589.[30]
P. K. Maini and J. D. Murray , A nonlinear analysis of a mechanical model for biological patternformation , SIAM Journal on Applied Mathematics, 48 (1988), pp. 1064–1072.[31]
P. K. Maini, L. Olsen, and J. A. Sherratt , Mathematical models for cell-matrix interactions duringdermal wound healing , International Journal of Bifurcation and Chaos, 12 (2002), pp. 2021–2029.[32]
P. K. Maini and T. E. Woolley , The Turing model for biological pattern formation , in The Dynamicsof Biological Systems, Springer, 2019, pp. 189–204.[33]
P. K. Maini, T. E. Woolley, R. E. Baker, E. A. Gaffney, and S. S. Lee , Turing’s model forbiological pattern formation and the robustness problem , Interface Focus, 2 (2012), pp. 487–496.[34]
D. Manoussaki , A mechanochemical model of angiogenesis and vasculogenesis , ESAIM: MathematicalModelling and Numerical Analysis, 37 (2003), pp. 581–599.2835]
G. E. Mase , Continuum mechanics , McGraw-Hill New York, 1970.[36]
H. Meinhardt , Models of Biological Pattern Formation , Academic Press, London, 1982.[37]
A. G. Moon and R. T. Tranquillo , Fibroblast-populated collagen microsphere assay of cell tractionforce: Part 1. continuum model , AIChE Journal, 39 (1993), pp. 163–177.[38]
P. Moreo, E. A. Gaffney, J. M. Garcia-Aznar, and M. Doblar´e , On the modelling of biologicalpatterns with mechanochemical models: insights from analysis and computation , Bulletin of MathematicalBiology, 72 (2010), pp. 400–431.[39]
J. D. Murray , A pre-pattern formation mechanism for animal coat markings , Journal of TheoreticalBiology, 88 (1981), pp. 161–199.[40]
J. D. Murray , Mathematical biology. II Spatial models and biomedical applications { InterdisciplinaryApplied Mathematics V. 18 } , Springer-Verlag New York Incorporated New York, 2001.[41] J. D. Murray and P. K. Maini , A new approach to the generation of pattern and form in embryology ,Scientific Programme, Oxford, (1986).[42] ,
Pattern formation mechanisms–a comparison of reaction-diffusion and mechanochemical models , inCell to Cell Signalling, Elsevier, 1989, pp. 159–170.[43]
J. D. Murray, P. K. Maini, and R. T. Tranquillo , Mechanochemical models for generating biologicalpattern and form in development , Physics Reports, 171 (1988), pp. 59–84.[44]
J. D. Murray and G. F. Oster , Cell traction models for generating pattern and form in morphogenesis ,Journal of Mathematical Biology, 19 (1984), pp. 265–279.[45] ,
Generation of biological pattern and form , Mathematical Medicine and Biology: A Journal of theIMA, 1 (1984), pp. 51–75.[46]
J. D. Murray, G. F. Oster, and A. K. Harris , A mechanical model for mesenchymal morphogenesis ,Journal of Mathematical Biology, 17 (1983), pp. 125–129.[47]
M. R. Myerscough, P. K. Maini, and K. J. Painter , Pattern formation in a generalized chemotacticmodel , Bulletin of Mathematical Biology, 60 (1998), pp. 1–26.[48]
S. Nasseri, L. E. Bilston, and N. Phan-Thien , Viscoelastic properties of pig kidney in shear, experi-mental results and modelling , Rheologica acta, 41 (2002), pp. 180–192.[49]
L. Olsen, J. A. Sherratt, and P. K. Maini , A mechanochemical model for adult dermal woundcontraction and the permanence of the contracted tissue displacement profile , Journal of Theoretical Biology,177 (1995), pp. 113–128.[50]
G. F. Oster, J. D. Murray, and A. K. Harris , Mechanical aspects of mesenchymal morphogenesis ,Development, 78 (1983), pp. 83–125.[51]
J. Park, D. S. Kim, T. S. Shim, C. M. Lim, Y. Koh, S. D. Lee, W. S. Kim, W. D. Kim, J. S.Lee, and K. S. Song , Lung cancer in patients with idiopathic pulmonary fibrosis , European RespiratoryJournal, 17 (2001), pp. 1216–1219.[52]
A. S. Perelson, P. K. Maini, J. D. Murray, J. M. Hyman, and G. F. Oster , Nonlinear patternselection in a mechanical model for morphogenesis , Journal of Mathematical Biology, 24 (1986), pp. 525–541. 2953]
M. Scianna, C. G. Bell, and L. Preziosi , A review of mathematical models for the formation ofvascular networks , Journal of Theoretical Biology, 333 (2013), pp. 174–209.[54]
J. G. Snedeker, P. Niederer, F. Schmidlin, M. Farshad, C. Demetropoulos, J. Lee, andK. Yang , Strain-rate dependent material properties of the porcine and human kidney capsule , Journal ofBiomechanics, 38 (2005), pp. 1011–1021.[55]
L. Tranqui and P. Tracqui , Mechanical signalling and angiogenesis. the integration of cell–extracellularmatrix couplings , Comptes Rendus de l’Acad´emie des Sciences-Series III-Sciences de la Vie, 323 (2000),pp. 31–47.[56]
R. T. Tranquillo and J. D. Murray , Continuum model of fibroblast-driven wound contraction:inflammation-mediation , Journal of Theoretical Biology, 158 (1992), pp. 135–172.[57]
A. M. Turing , The chemical basis of morphogenesis , Philosophical Transactions of the Royal Society ofLondon B, 237 (1952), pp. 37–72.[58]
S. Urdy , On the evolution of morphogenetic models: mechano-chemical interactions and an integratedview of cell differentiation, growth, pattern formation and morphogenesis , Biological Reviews, 87 (2012),pp. 786–803.[59]
D. Valtorta and E. Mazza , Dynamic measurement of soft tissue viscoelastic properties with a torsionalresonator device , Medical Image Analysis, 9 (2005), pp. 481–490.[60]
S. Van Helvert, C. Storm, and P. Friedl , Mechanoreciprocity in cell migration , Nature cell biology,20 (2018), p. 8.[61]