Constitutive model for the rheology of biological tissue

The rheology of biological tissue is key to processes such as embryo development, wound healing and cancer metastasis. Vertex models of confluent tissue monolayers have uncovered a spontaneous liquid-solid transition tuned by cell shape; and a shear-induced solidification transition of an initially liquid-like tissue. Alongside this jamming/unjamming behaviour, biological tissue also displays an inherent viscoelasticity, with a slow time and rate dependent mechanics. With this motivation, we combine simulations and continuum theory to examine the rheology of the vertex model in nonlinear shear across a full range of shear rates from quastistatic to fast, elucidating its nonlinear stress-strain curves after the inception of shear of finite rate, and its steady state flow curves of stress as a function of strain rate. We formulate a rheological constitutive model that couples cell shape to flow and captures both the tissue solid-liquid transition and its rich linear and nonlinear rheology.

The rheology of biological tissue is crucial to processes such as morphogenesis, wound healing and cancer metastasis.On short timescales, tissues withstand stress in a solid-like way.On longer timescales, they reshape via internally active processes such as cell shape change, rearrangement, division and death [1,2].Tissues are thus viscoelastic [3].Power law stress relaxation [4,5] and slow oscillatory cell displacements [6] after straining underline their rate dependent mechanics.Tissues furthermore undergo spontaneous solid-liquid transitions [7][8][9][10][11] driven by both active processes, such as fluctuations of cell-edge tensions, motility and alignment, and geometric constraints [12], with important implications for morphogenesis and cancer progression.Nonlinear rheological response to tensile stretching includes stiffening [13] or fluidization [14] of single cells, and stiffening then rupture of tissue monolayers [15].Internal activity can likewise induce nonlinear phenomena such as superelasticity [16] and fracture [17].
Understanding tissue rheology theoretically is thus of major importance.Well studied vertex and Voronoi models [9,18,19] of confluent tissue, with no gaps between cells, represent a 2D tissue monolayer as a tiling of polygonal cells.They capture a density-independent solidliquid transition tuned by a parameter characterising the target cell shape, which in turn embodies the competition between cortex contractility and cell-cell adhesion [7][8][9].Vertex models have also been used to study the linear mechanics of tissues [20][21][22], and their response to nonlinear stretch [23] and shear [24][25][26][27].Recently, vertex model simulations of a tissue that is fluid-like in zero shear demonstrated a shear-induced rigidity transition above a critical strain, applied quasistatically [27].
While vertex models and other mesoscopic models have played an important role in advancing our understanding of tissue mechanics, it is also helpful to develop coarse grained continuum rheological constitutive models.Early work formulated a continuum model that couples cell shape and cell motility, capturing some of the glassy dynamics of tissue [28].Inspired by early hydrodynamic theories of active fluids and gels [29,30], con-tinuum constitutive models have been developed to characterize the role of cell shape change, rearrangements, division and death in morphogenesis [2,25,[31][32][33][34][35].
Still lacking, however, is a continuum hydrodynamic constitutive model capable of describing both the spontaneous solid-liquid transition of confluent tissues and its rheological response to external deformation and flow.Inspired by mean-field theories of cell-shape driven transitions [22,27,28] and by fluidity models of the rheology of dense soft suspensions [36], we introduce such a model.
The key new insights of our approach are as follows.First, we distinguish the role of geometric frustration (encoded in the cell perimeter p), from that of T1 topological rearrangements (encoded in our fluidity variable a).The former is key to the zero-shear liquid-solid transition and (when coupled to our orientation tensor σ ij ) strain stiffening at small to modest imposed strains [27].The latter cause the plasticity associated with the stress overshoot at imposed strains O(1), and the ultimate steady flowing state.Second, in modeling the geometric frustration, we distinguish a tensor characterizing individual cell shape (of which p is the trace), and a tensor characterizing the average cell orientation at the tissue scale [28].
We furthermore submit this new continuum model to stringent comparison with simulations across a full range of shear rates from quasi-static to fast.We demonstrate our continuum model to capture both the zero-shear solid-liquid transition and strain stiffening transitions reported in Ref. [27], the full nonlinear stress vs. strain behavior after the inception of shear, and the steady state flow curves of stress vs. shear rate.
Vertex model simulations -The vertex model [18,19] represents the tightly packed confluent cells of a 2D tissue monolayer as c = 1 • • • N c polygons that tile the plane.Each cell is defined by the location of its n c = 1 • • • ν c vertices, with any two neighbouring vertices α and β connected by an edge of length ℓ αβ .The elastic energy of the tissue is controlled by the interplay of pressure within each cell and tension along the cell edges.Assuming the cell-edge tension per unit length is uniform across the tissue, the energy can be written as where each cell experiences an energy cost for deviation of its area A c and perimeter P c from target values A c0 and P c0 , with area and perimeter stiffness κ A and κ P .The first term on the RHS models 3D cell volume incompressibility via an effective 2D area elasticity [19,37].
The second describes the competition between cell cortical contractility and adhesion between neighbouring cells in controlling cell-edge tension and perimeter [7,19,37].We denote by ⃗ F n = − δE δ⃗ xn the total force on the n th vertex of the tiling at position ⃗ x n due to interactions with all other vertices.In an applied shear of rate γ, with flow direction x and shear gradient y, we assume over-damped dynamics with drag ζ, d⃗ xn dt = ζ −1 ⃗ F n + γy n ⃗ x, with Lees-Edwards periodic boundary conditions.The cells also undergo T1 topological neighbor exchanges that allow the tissue to plastically relax stresses [9,[38][39][40].
To focus on amorphous tissue structures, we simulate a 50 : 50 bidisperse tiling of N c = 4096 cells of target areas A 0 = 1, 1.4, which sets our length unit.We adjust P c0 for the two cell populations to maintain the target cell shape p 0 = P c0 / √ A c0 the same for all cells.We choose units in which κ A = 1 and ζ = 1 and set κ p = 1.0 throughout.We vary p 0 and the imposed shear rate γ.As an initial condition, we seed a planar Voronoi tiling then evolve the above dynamics to steady state in zero shear.At time t = 0, we switch on shear and measure the shear stress Σ ij (t) = 1 N N n=1 F ni x nj , where the sum is over all N vertices in the tiling, and the mean cell perimeter p(t) = 1 Nc Nc c=1 p c .Denoting by ⃗ t nc the unit vector along the edge of length l nc between the n c th and (n c + 1)th vertices of cell c, we define a single-cell shape tensor σ c ij = 1 νc νc n=1 l nc t nc i t nc j , where the sum is over the ν c vertices of the c-th cell, and the tissue-scale averaged orientation tensor We use the same notation Σ ij , σ ij , p for the counterpart coarse-grained quantities in our constitutive model below.
In the absence of external stress, the vertex model exhibits a liquid solid transition as a function of the target shape p 0 [7,41].For p 0 < p * 0 the energy barriers to T1 transitions are finite and the system is a solid with a finite zero-frequency linear shear modulus.At the critical value p * 0 , the mean energy barrier for T1 transitions vanishes, giving liquid response for p 0 > p * 0 .For our bidisperse tiling, p * 0 = 3.85.For monodisperse disordered polygons p * 0 ≃ 3.81, a value close to that of a regular pentagon [7].This value is renormalized by motility [9] and by cell alignment with local spontaneous shear [42].It was recently realized that this transition has a geometric origin associated with the underconstrained nature of the energy in Eq. 1 [20,22,43].For regular hexagons the transition occurs at the isoperimetric value p iso = 8 √ 3 ≃ 3.722.Below this value it is not possible to satisfy both target area and perimeter and the ground state has p = p * 0 and finite energy.This is the solid or incompatible state.For p 0 > p iso there is a family of zero energy area and perimeter preserving ground states, with p = p 0 .The system can accommodate an externally applied linear shear by adjusting its shape within this degenerate manifold [22].The compatible system is therefore a liquid with zero shear modulus, although it stiffens and acquires rigidity at finite strains [27].
Constitutive model -We now construct a continuum model that accounts for the mean-field liquid-solid transition, and also captures the key rheological features of the vertex model: (i) reversibility of linear response to small strains, (ii) strain stiffening at intermediate strains, (iii) plastic relaxation at larger strains, due to T1 cell rearrangements, and (iv) a yield stress in the steady state flow curve Σ( γ), as obtained in Ref. [27].Although our model below is cast in frame invariant form, capable of addressing any flow, we focus on response to simple shear, to compare with our vertex model simulations.
We assume dynamics of the cell perimeter governed by: with α and τ p constants and invariant strain rate γ = 2D ij D ij .In the absence of shear, p relaxes on a timescale τ p to a steady state that displays a transcritical bifurcation as a function of the target cell perimeter p 0 , with p = p * 0 in the solid phase p 0 < p * 0 and p = p 0 in the liquid phase p 0 > p * 0 , capturing the liquid-solid transition [7].The same transcritical structure emerges by writing exact equations for the relaxation of a single cell modeled as a regular n−sided polygon according to the vertex model dynamics prescribed above.
In shear, the perimeter is advected by flow and stretched by the shear rate γ.In addition, the coupling ασ ij σ ij captures a key intuition of our approach: that a shear-induced global cell orientation σ ij provides an effective mean field that distorts the individual cell's shape p away from its zero-shear value.As a result, in the solid phase p increases relative to its zero shear value p = p * 0 from the outset of straining.In the liquid phase, p increases relative to its zero shear value p = p 0 only after a critical strain amplitude γ c , capturing the strain-induced stiffening transition [27].The behavior introduced by the coupling of single-cell shape, as quantified by the mean perimeter p, to the tissue-scale cell shape σ ij is analogous to the influence of cell alignment due to internally generated stresses in Drosophila germband extension [42].Indeed, the form of coupling of p to σ ij in Eqn. 2 is justified both by experiment [42] and mean field theory [22,27].
The cell orientation tensor is taken to obey an evolution equation of the widely used Maxwellian form, where K ij = ∂ j v i is the strain rate tensor and D ij = 1 2 (K ij + K ji ).The last term in Eq. 3 describes plastic relaxation.It vanishes in linear response (small strains), where a = 0 (see below), allowing the orientation tensor σ ij to build linearly and reversibly with strain, as expected in the absence of plastic T1 events.Consistent with previous studies of the vertex model [19,22,38] we write the deviatoric stress tensor Here C is constant and p 0 the target cell perimeter.
In linear response (small strains), the effective modulus G 0 = C(p − p 0 ) is non-zero in the solid phase, where p > p 0 , and zero in the liquid phase, where p = p 0 .
Were the factor a on its RHS a constant inverse relaxation time, Eq. 3 would be the widely used Maxwell model, capturing viscoelasticity, but not the irreversible plasticity of T1 events.To model plasticity, we take a to be a fluidity-like variable [36] with dynamics: with f ( γ) = β γ/(1 + 1 2 τ 0 γ), in which β is constant and τ 0 a microscopic time.As suited to an athermal tissue, with no relaxation events induced by temperature or activity (no cell motility, division or death), this is a purely straindriven dynamics.In linear response, a = 0, giving a reversible dependence of σ ij on strain.In weak shear, a builds on a strain O(1) to model the plasticity of T1 events via the final term in Eq. 3. In steady weak shear a = f ( γ) ≈ β γ, giving a divergent relaxation time 1/a as γ → 0, and a yield stress in the steady state flow curve.
We explore different values of shear rate, γ, and the target perimeter p 0 relative to the transition p * 0 .(See Appendix for model parameters.)We prescribe as initial condition to shear a perimeter p(t = 0) equal to its steady state value in zero shear, an orientation tensor σ ij (t = 0) = 0, and fluidity a(t = 0) = 0. We then switch on a simple shear K ij = γδ iy δ jx at time t = 0 and track the evolution of p, σ xy and Σ xy as a function of time t or equivalently (to within a constant factor γ) accumulating strain γ = γt.Hereafter we drop the xy subscript, writing σ xy = σ and Σ xy = Σ.
Results -Our constitutive model captures the liquidsolid transition as a function of target cell shape in zero shear [7] and the shear-induced rigidity transition of the liquid-like tissue, above a critical shear strain, applied quasistatically γ → 0 [27].See Fig. 1a, which shows the shear stress Σ vs. strain γ in shear at rate γ = 10 −6 .At small strains, just after the inception of shear, the modulus G 0 = dΣ/dγ| γ=0 is finite (solid-like) for p 0 < p * 0 but zero (liquid-like) for p 0 > p * 0 (Fig. 1b, dashed line).In the liquid phase, the stress Σ and slope dΣ/dγ first become non-zero above a nonlinear critical shear strain γ c , heralding a strain-induced stiffening transition (solid line in Fig. 1b, defined as the strain at which the stress first exceeds 10 −5 at any p 0 .) Having explored quasistatic shear, we now consider nonlinear shear flow across a full range of shear rates from quasi-static to fast.The evolution of Σ, σ and p as a function of strain since the inception of shear is shown in Fig. 2, for a range of p 0 below and above p * 0 .The left column shows the results of vertex model simulations.The right shows the predictions of our constitutive model, which performs well in capturing all the qualitative features of the simulations.
At small strains, just after shearing starts, the effective modulus G 0 = dΣ/dγ| γ=0 is finite in the solid phase, p 0 < p * 0 , but small in the liquid phase, p 0 > p * 0 .Indeed, repeating the simulations for progressively lower strain rates γ → 0 in the solid phase, G 0 tends to a non-zero constant, G 0 (p 0 , γ → 0), consistent with the quasistatic results discussed above.In the liquid phase, G 0 → 0 as γ → 0, again consistent with the quasistatic results.
At higher strains, γ = O(1), strain stiffening is observed: the slope of Σ vs γ increases with increasing γ.This is particularly pronounced in the liquid phase, p 0 > p * 0 , where the effective modulus dΣ/dγ was very small at small strains (tending to zero as γ → 0, as just discussed), but becomes appreciable after a strain γ = O(1) (even in the limit γ → 0).After this regime of strain stiffening, the stress overshoots slightly before declining to a constant in the final state of steady flow.
This rich behaviour is readily understood within our simple constitutive model.The initial fluid-like behaviour for p 0 > p * 0 arises because p = p 0 before shearing commences, giving zero effective modulus in Eq. 4. As strain increases, tissue deformation is captured by the growth of σ, which in turn yields an increase of p relative to its equilibrium value due to the coupling term in α in Eq. 2. This is also responsible for the less pronounced strain stiffening in the solid phase, p 0 < p * 0 .The subsequent overshoot in stress Σ (and perimeter p) at larger strains is caused by the overshoot in the cell orientation σ seen in the middle panels of Fig. 2. The stress decline after overshoot arises in the vertex model from plastic relaxation via T1 events, an effect captured in the constitutive model via an increase of fluidity a with shear.The tissue shape tensor σ is essentially independent of p 0 in the vertex model (at low strain rates), consistent with the lack of any coupling of the evolution equation for σ ij to p in the constitutive model.
At long times, t → ∞, after many strain units γ = γt → ∞, a state of final plastic flow is reached in which each of Σ, σ and p attains a steady value.This is reported as a function of γ in Fig. 3, for the vertex model (left column), and constitutive model (right), with good semi-quantitative agreement.In rheological parlance, the steady state relationship Σ = Σ( γ) is termed the "flow curve".The vertex model flow curves show a dynamical yield stress: a non-zero limiting intercept lim γ→0 Σ( γ) = Σ Y ̸ = 0. Importantly, this is true both for p 0 < p * 0 and for p 0 > p * 0 : whereas liquid and solid states are distinct and separated by a transition at small strains, in steady nonlinear shear, however slow, the vertex model displays a non-zero yield stress up to a larger p 0 = p * * 0 > p * 0 [27], as also seen in Fig. 1.This is easily understood within our constitutive model.In steady shear, Eq. 5 predicts the fluidity a = f ( γ) = β γ/(1 + 1 2 τ 0 γ).Combining with Eq. 3 for the orientation gives σ = γ/a = 1 β (1 + 1 2 τ 0 γ).Were we to assume p − p 0 = 1, independent of strain rate, we would obtain a flow curve Σ( γ) = C β (1 + 1 2 τ 0 γ), with a yield stress σ Y = C β as γ → 0 and Newtonian behaviour Σ ∝ γ as γ → ∞.The actual flow curve is modified somewhat in comparison, due to the strain rate dependence of p − p 0 .Importantly, however, it retains a yield stress because p ̸ = p 0 in steady flow, even in the limit γ → 0: the perimeter is always strongly perturbed from its unsheared value, due to the coupling ασ ij σ ij in Eq. 2. Intuitively, the key effect of a steady shear, even when applied quasistatically, is to deform cells away from their target shape such that they carry a stress and the liquid phase seen at small strains is destroyed.
Conclusions -We have presented a continuum constitutive model for the rheology of confluent 2D biological tissue and demonstrated it to capture the rich rheophysics seen in simulations of the vertex model under applied shear.This includes strain-stiffening of the liquid above a critical strain, a stress overshoot at larger strains due to the plasticity of T1 rearrangements, and a finite yield stress in steady shear, even in the (zero-shear) liquid phase.Our model includes the effects of cell shape change and rearrangements on mechanical behaviour, and will provide a useful phenomenological framework for modeling the rheology of biological tissue.Elucidating its predictions in deformation protocols besides simple shear is left to future work, as are extensions to incorporate other active processes such as cell motility, division and death.
Model parameters are the modulus-like quantity C, the microscopic time τ 0 and the parameter β in the function f for the fluidity, the transition value of the target perimeter p * 0 , the coupling of perimeter to orientation α, and the perimeter relaxation time τ p .We choose units C = 1 and τ 0 = 1, and treat p * 0 , α, β and τ p as fitting parameters in comparing our constitutive model with the vertex model simulations.We have found p * 0 = 3.85, α = 0.36, β = 2.0 and τ p = 0.1 to give the best fit.Among these, p * 0 is the value of p 0 at the (zero-shear) liquid-solid transition.Accordingly, we set the value of p * 0 in our continuum model to that value found in our vertex model simulations.β sets the quasistatic limit of the shear component of cell orientation tensor, lim γ→0 σ xy = 1/β, with β = 2.0 in our vertex model simulations.α sets the effective modulus G(p − p 0 ) in the shear induced solid phase, with p − p 0 = p * 0 − p 0 + ασ ij σ ij as γ → 0, and accordingly sets the flow curve's yield stress, lim γ→0 Σ( γ).We choose α to give the best fit of the continuum model's yield stress to that of the vertex model simulations.Finally, τ p controls the steepness of the flow curve at high strain rates (where the vertex model is likely to become less reliable ) and the small finite value ∼ γτ p of the stress before the true quasistatic strain-stiffening transition.
The numerical timestep is Dt = Dtl min /F max with l min the minimum edge length, F max the maximum vertex force and Dt = 0.01.T1 events are triggered below a critical edge length l c = 0.01.

Vertex model simulations
The vertex model [18,19] represents the confluent cells of a tissue monolayer via c = 1 • • • N c polygons that tile the plane.Each cell is defined by the location of its n c = 1 • • • ν c vertices, with any two neighbouring vertices connected by an edge.Each vertex belongs to three neighbouring cells (or transiently four, during a T1 event, see below), with three (or four) edges stemming from it accordingly.Each edge belongs to two neighbouring cells.
Consider the n c th and (n c + 1)th vertices of cell c. Cell c contributes to these two vertices an equal and opposite force of magnitude, κ P (P c −P 0 ), acting as a tension along the edge that connects them.This models a competition between cell cortical contractility and adhesion between neighbouring cells [37], with κ P an elastic constant [19], P c the cell perimeter and P 0 its target value [7].Cell c furthermore contributes to the same two vertices a force of magnitude κ A (A c − A 0 )l nc , acting as a pressure in the direction of the edge normal outwards from cell c, with l nc the length of the edge connecting the vertices.Physically, this models 3D cell volume incompressibility via an effective 2D area elasticity of constant κ A , with A c the cell area and A 0 its target value [19,37].Fig. 4 shows a sketch of these forces.Each of the two vertices also belongs to two other neighbouring cells (or three, during T1 events), which contribute forces likewise.
For the nth vertex of all N in the tiling, we denote the sum of the forces from each of its associated cell edges by F n .In shear of rate γ, with flow direction x and shear gradient y, the vertex position ⃗ x n obeys over-damped dynamics with drag coefficient ζ as a function of time t: with Lees-Edwards periodic boundary conditions.
Consider the vertex at the junction between cells αβγ and a neighbouring vertex between cells αβδ.When the edge connecting these vertices shrinks below a small length l c , a T 1 transition removes these two vertices and replaces them with new ones at the junctions of cells βγδ and αγδ, conferring plastic cell rearrangement.

Componentwise constitutive equations in simple shear
In the main text, we presented our constitutive model in tensorial, frame-invariant form, capable of addressing any imposed deformation or flow protocol.Here we extract the components of those equations relevant to homogeneous imposed simple shear flow with K ij = γδ iy δ jx .
The cell perimeter evolves according to: where α and τ p are constants, and the trace of the cell orientation tensor σ ii = σ xx + σ yy .
The orientation tensor obeys an evolution equation of the widely used Maxwellian form, which componentwise is written as: The shear stress where C is constant and p 0 the target cell perimeter.
To model plasticity, we take the quantity a in the evolution of the orientation tensor to be a fluidity-like variable [36] with its own dynamics: with and f ( γ) = β γ/(1 + 1

FIG. 2 .
FIG. 2. Rheological behaviour of the vertex model (left) and constitutive model (right) in shear startup at a shear rate γ = 10 −3 for values of the target perimeter p0 = 3.50, 3.55, 3.60 • • • 4.00 (in black, red, green • • • orange curves downwards; curve for p * 0 = 3.85 in purple).Shown is the evolution of the shear stress (top), shear component of the orientation tensor (middle) and cell perimeter (bottom) as a function of accumulating strain γ = γt.

FIG. 3 .
FIG. 3. Steady state (t → ∞) dependence of the shear stress (top), shear component of the orientation tensor (middle) and cell perimeter (bottom) for the same values of the target perimeter p0 as in Fig. 2, with the same line colour coding.Results are shown for the vertex model in the left column and the constitutive model in the right column.