# Colloquium: Mechanical formalisms for tissue dynamics

## Abstract

The published version of this work,
Tlili et al., Eur. Phys. J. E 38, 33-63 (2015),
is freely accessible (PDF)
from the “Highlight and Colloquia” page of the journal website at
http://epje.epj.org/epje-news-highlights-colloquia

Changes made in the present version
as compared to the published version
are indicated in red color.
They all concern Appendix D.2
devoted to the dissipation function formalism at large deformations.
In Eq. (235) (formerly D50),
the notation
was corrected to .
In Eq. (236) (formerly D51),
the small deformation notation
was replaced with .
The paragraph concerning consitutive equations for large deformations
(D.2.2) was reformulated for clarity
and one error concerning Eqs. (227,228)
(formerly D37 and D39) was corrected.
Correspondingly, some different equations were cited
in Sections D.2.3 and D.2.4.

Abstract. The understanding of morphogenesis in living organisms
has been renewed by tremendous progress in experimental techniques
that provide access to cell-scale, quantitative information
both on the shapes
of cells within tissues
and on the genes being expressed.
This information suggests that our understanding of the respective
contributions of gene expression and mechanics, and of their crucial
entanglement, will soon leap forward.
Biomechanics increasingly benefits from models, which assist the design
and interpretation of experiments, point out the main ingredients and
assumptions, and ultimately lead to predictions.
The newly accessible local information thus calls
for a reflection on how to select suitable classes of mechanical models.
We review both mechanical ingredients
suggested by the current knowledge of tissue behaviour,
and modelling methods that can help generate
a rheological diagram or a constitutive equation.
We distinguish cell scale (“intra-cell”) and tissue scale (“inter-cell”) contributions.
We recall the mathematical framework developped for continuum materials
and explain how to transform a constitutive equation
into a set of partial differential equations
amenable to numerical resolution.
We show that when plastic behaviour is relevant, the dissipation function formalism appears
appropriate to generate constitutive equations; its variational nature
facilitates numerical implementation, and we discuss
adaptations needed in the case of large deformations.
The present article gathers theoretical methods that can readily
enhance the significance of the data to be
extracted from recent or future high throughput biomechanical experiments.

Contact: cyprien.gay@univ-paris-diderot.fr, francois.graner@univ-paris-diderot.fr

###### pacs:

87.19.R- Mechanical and electrical properties of tissues and organs 87.19.lx Development and growth 83.10.Gr Constitutive relations 83.60.La Viscoplasticity; yield stress## I Introduction

### i.1 Motivations

While biologists use the word “model” for an organism studied as an archetype, such as Drosophila or Arabidopsis, physicists rather use this word for models based on either analytical equations or numerical simulations. Biomechanical models have a century-old tradition and play several roles [1]. For instance, they assist experiments to integrate and manipulate quantitative data, and extract measurements of relevant parameters (either directly, or through fits of models to data). They also lead to predictions, help to propose and design new experiments, test the effect of parameters, simulate several realisations of a stochastic phenomenon, or simulate experiments which cannot be implemented in practice. They enable to illustrate an experiment, favor its interpretation and understanding. They point out the main ingredients and assumptions, test the sensitivity of an experiment to a parameter or to errors, and determine which assumptions are sufficient to describe an experimental result. Models can help to determine whether two facts which appear similar have a superficial or deep analogy, and whether two facts which are correlated are causally related or not.

Two themes have dominated the recent literature: modelling the mechanics of some specific adult tissues like bones or muscles, for which deformations and stresses are obviously part of the biological function [2]; [3]; and unraveling the role of forces in the generation of forms during embryonic development [4]; [5]. During the last decades, both the physics and the biology sides of the latter question have been completely transformed, especially by progress in imaging.

On the physics side,
new so-called “complex”
materials with an internal structure, such as foams, emulsions or gels,
have been thoroughly studied,
especially in the last twenty years,
with a strong emphasis on the difficult problem of the feedback between
the microscopic structure and the mechanical response [6]; [7].
The development
of new tools to image the changes in the microstructure
arrangement under
well-controlled global stresses or deformations has provided a wealth of data.
Modelling has played a crucial part
via the determination of so-called “constitutive equations”.
A constitutive equation characterises the local properties of a
material within the framework of continuum mechanics. It
relates dynamical quantities, such as the stress carried by the material, with
kinematical quantities, *e.g.*
the deformation (also called “strain”) or the deformation rate.

On the biology side, questions now arise regarding the interplay of cell scale behaviour and tissue scale mechanical properties, among which the following two examples.

A first question is: how does a collective behaviour, which is not obviously apparent at the cell scale, emerge at the tissue scale? Analyses of images and movies suggested that epithelia or whole embryos behave like viscous liquids on long time scales [8]. The physical origin and the value of the (effective) viscosity should be traced back to the cell dynamics: it can in principle incorporate contributions from ingredients such as cell divisions and apoptoses [9] or cell contour fluctuations [10]; [11], but also from orientational order, cell contractility, cell motility or cell rheological properties. All these local and sometimes changing ingredients become progressively accessible to experimental measurements. Biomechanical models can investigate the bottom-up relationship between local cell-scale structure and tissue-scale mechanical behaviour, unraveling the signature of the cellular structure in the continuum mechanics descriptions [10]; [11].

A second question is: how can the mechanical state of the tissue have an influence on the cell division rate [12]; [13], or on the orientation of the cells undergoing division [14]? In addition, the mechanical state of the tissue can generate cell polarity and hence an anisotropy of the local cell packing, which may affect the mutual influence between the local mechanics (forces and deformations) and the cell behaviour. Biomechanical models contribute to disentangle these complex feedback loops and address such top-down relationships.

To address these questions, a natural strategy is first to reconstruct the mechanics from the structural description, then to investigate the feedback between well-identified mechanical variables and the expression of specific genes. In particular, this interplay between genes and mechanics is expected to be the key to the spontaneous construction of the adult form in a developing tissue without an organising center. Such problem in its full complexity will probably require a “systems biology” approach based on large scale mapping of expression for at least tens of genes, coupled to a correct mechanical modelling on an extended range of scales in time and space, which in turn supposes experimental setups able to produce the relevant genetic and mechanical data.

### i.2 State of the art

Recent developments in *in vivo* microscopy
yield access
to the same richness of structural information for living tissues as
has already been the case for complex fluids.
The biology of cells and tissues
is now investigated in detail in terms of protein distribution and gene
expression, especially during development [15].
It is possible to image
the full geometry of a developing
embryonic tissue at cellular resolution
[16]; [17]; [18]; [19]; [20],
while
visualising the expression of various
genes of interest [21]; [22]; [23].
Mechanical fields such as
the deformation, deformation rate or plastic deformation rate are increasingly
accessible to direct measurement.
Several fields can be measured quantitatively at least up to an unknown prefactor.
This is the case for:
distributions of proteins (involved in cytoskeleton, adhesion or force production), via quantitative fluorescence [24];
elastic forces and stresses, either
by laser ablation of cell junctions [25] and
tissue pieces
[26], or through image-based force inference
methods
[27]; [28]; [29]; [30];
even viscous stress fields, indirectly estimated
[23].
Other methods
include absolute measurements of forces based
on micro-manipulation [31], in situ incorporation of deformable force sensors [32]
or fluorescence resonance energy transfer (FRET) [33].

*In vitro* assemblies of cohesive cells are useful experimental materials.
Within a reconstructed cell assembly,
each individual cell
retains its normal physiological behaviour:
it can grow, divide, die, migrate.
In the absence of any regulatory physiological context,
cells display small or negligible variation of gene expression.
Reconstructed assemblies
thus allow to separate
the mechanical behaviour of a tissue from its feedback
to and from genetics. Furthermore,
in absence of any coordinated variation of
the genetic identity of constituent cells,
spatial homogeneity may be achieved.
Simple, well-controlled
boundary conditions can be implemented by a careful choice of the geometry,
either in two or three dimensions.

In two dimensions, confluent monolayers are usually grown on a substrate used both as a source of external friction and as a mechanical sensor to measure local forces [34]; [35]; [36]. 2D monolayers facilitate experiments, simulations, theory and their mutual comparisons [37]; [38]; [39]; [40]; [41]. 2D images are easier to obtain and can be analysed in detail; data are more easily manipulated, both formally and computationally.

In three dimensions, multi-cellular spheroids in a well-controlled, *in vitro* setting
[10]; [12]; [42]; [43]
are a good material
to mimic
the mechanical properties of tumors, and of
homogeneous parts of whole organs, either adult or during development.
They are also useful for rheological studies [44]; [45],
especially since they are free from contact with a solid substrate.
Although the full reconstruction of the geometry of multi-cellular
spheroids at cellular resolution remains challenging,
it has progressed in recent years
[46].

### i.3 Outline of the paper

A tissue can be seen as a cellular material, active in the sense that it is out of equilibrium due to its reservoir of chemical energy, which is converted into mechanical work at or below the scale of the material’s constituents, the cells. A consequence of this activity is that force dipoles, as well as motion, are generated autonomously.

Our global strategy is as follows. We construct rheological diagrams based on insights concerning the mechanics of the biological tissue of interest. One of the main insights is a distinction between intra-cell mechanisms: elasticity, internal relaxation, growth, contractility, and inter-cell mechanisms: cell-cell rearrangement, division and apoptosis. The transcription of the rheological diagrams within the dissipation function formalism provides local rheological equations. We show how this local rheology should be inserted into the balance equations of continuum mechanics to generate a complete spatial model expressed as a set of partial differential equations. This procedure is conducted not only in the usually treated case of small elastic deformations, but also in the relevant, less discussed, case of large elastic deformations.

Our hope is to provide a functional and versatile toolbox for tissue modelling.
We would like
to guide the choice of approaches and of models according to
the tissue under consideration, the experimental set-up, or the scientific
question raised.
We propose
a framework for a
tensorial treatment of spatially heterogeneous tissues.
It is suitable to incorporate the data arising from
the analysis of
experimental data,
which are increasingly often live microscopy movies.
Although the simplest applications concern *in vitro* experiments,
often performed with epithelial cells, the same approach applies
to a wide spectrum of living tissues including animal tissues during development, wound healing, or tumorigenesis.

This article is organised as follows. Section II makes explicit our assumptions and arguments, then details the use of the dissipation function formalism. Section III reviews mechanical ingredients suitable for the theoretical description of a wide range of living tissues, both in vitro and in vivo, illustrated with worked out examples chosen for their simplicity. Section IV groups these ingredients to form models of more realistic applications. Section V summarizes and opens perspectives.

Appendix A examines the link between the scale of discrete cells and the scale of the continuous tissue. Appendix B provides more details on how to write and use equations within the dissipation function formalism. Appendix C provides further examples of coupling with non-mechanical fields. Appendix D examines the requirements to treat tissue mechanics at large deformation.

## Ii Choices and methods

In this Section, we explain our choices and our assumptions. Section II.1 compares discrete and continuum approaches. Section II.2 compares rheological diagrams, hydrodynamics, and the dissipation function formalisms. Section II.3 discusses specific requirements to model cellular materials. Section II.4 suggests how to incorporate space dependence in constitutive equations to write partial differential equations.

### ii.1 A continuum rather than a discrete description

Models that describe tissue mechanical properties may be broadly split into two main categories: bottom-up “cell-based” simulations and continuum mechanics models.

Direct cellular simulations are built upon the (supposedly known) geometry and rheology of each individual cell and membrane. They generate a global tissue behaviour through the computation of the large-scale dynamics of assemblies of idealised cells [47]; [48]; [49]; [50]; [51]; [52]; [53]; [54]; [55]; [56]; [57]. Simulations enable to directly test the collective effect of each cell-scale ingredient, and of their mutual feedbacks. Also, they work well over a large range of cell numbers, up to tens of thousands, and down even to a small number of cells, where the length scale of a single cell and that of the cell assembly are comparable.

A continuum approach requires the existence of an intermediate length scale, larger than a typical cell size, yet smaller than the tissue spatial extension, and beyond which the relevant fields vary smoothly. The continuum rheology is captured through a constitutive equation relating the (tensorial) stresses and deformations [58]; [59]; [60]; [61]; [62]; [63]; [64]; [65]; [66]. This rheological model is incorporated into the usual framework of continuum mechanics using fundamental principles such as material and momentum conservation. Note that continuum models have also been applied to vegetal tissues, as in plant growth [58]; [67].

Both categories are complementary and have respective advantages.
For the purpose of the present paper, we favor the continuum approach,
which incorporates more easily precise details of cellular rheology.
When it succeeds, a continuum approach yields a synthetic grasp
of the relevant mechanical
variables on an intermediate scale (*i.e.* averaged over many individual
cells), and helps dealing with large tissues.
It often involves a smaller number of independent parameters
than a discrete approach, and this
helps comparing with experimental observations.

In order to test and calibrate a continuum model, it is generally necessary to extract continuum information from other sources such as discrete simulations or experiments on tissues with cell-scale resolution. Analysis tools have been developed in recent years to process segmented experimental movies in order to extract tensorial quantities from cell contours, which might include the elastic deformation and the plastic deformation rate [68]; for completeness, these tools are recalled in Appendix A.1.

The continuum models describing amorphous cellular materials can be tensorial and can incorporate viscoelastoplastic behaviour [58]; [64]; [69]; [70]; [71]. In addition to the viscous and elastic behaviour expected for a complex fluid that can store elastic deformation in its microstructure, one incorporates the ingredient of plasticity. It captures irreversible structural changes, more specifically here (i) local rearrangements of the individual cells (see Fig. 2), (ii) cell division or (iii) plasticity within the cells. A viscoelastoplastic model assembling these ingredients could capture the different “short-time” phenomena described in Section III and at the same time display a viscous liquid-like behaviour on longer time scales.

### ii.2 Choice of the dissipation function formalism

In this Section, we discuss three complementary general formalisms used to construct constitutive equations, indifferently in two or three dimensions: rheological diagrams (Section II.2.1), hydrodynamics in its general sense (Section II.2.2) and dissipation function (Section II.2.3). For simplicity we consider here stress and other mechanical tensors only with symmetrical components, while in general these three formalisms can include antisymmetric contributions when required.

#### Rheological diagram formalism

Fig. 1 represents a classical example of rheological diagram: the Oldroyd viscoelastic fluid model [73]. It consists in a dashpot with viscosity and deformation carrying the stress in parallel with a Maxwell element carrying the stress . This Maxwell element is itself made of a spring (stiffness , deformation ) in series with a dashpot (viscosity , deformation ). The elementary rheological equations read:

(1) |

Eliminating , , and between Eqs. (1) yields:

(2) |

Eq. (2) is a constitutive equation for the ensemble. Such a straightforward method is useful when physical knowledge or intuition of the material and its mechanical properties is sufficient to determine the topology (nodes and links) of the diagram.

Note that the relationship between a rheological diagram
and a constitutive equation
is *not* one-to-one.
For instance, a Maxwell element (linear viscoelastic liquid)
in parallel with a dashpot is a diagram distinct from a Voigt element (linear viscoelastic solid)
in series with a dashpot, but both are associated to the same constitutive equation.
Section III.3 presents another example.

#### Hydrodynamic formalism

When non-mechanical variables are present, the rheological diagram formalism (Section II.2.1) is not sufficient to establish the constitutive equation. Another formalism is necessary to include couplings between mechanical and non-mechanical variables.

A possible formalism is linear out-of-equilibrium thermodynamics, also called “hydrodynamics” [74] although its range of application is much larger than the mechanics of simple fluids. This approach has been highly successful, leading for instance to the derivation of the hydrodynamics of nematic liquid crystals (with the nematic director field as an additional, non-conserved hydrodynamic variable) [75]; [76]; the collective movements of self-propelled particles [77]; [78] (which have been suggested to be analogous with tissues dynamics [56]); or, more recently, of soft active matter (where a chemical field typically couples to an orientational order parameter to be modeled separately, e.g. cytoskeletal mechanics) [79]; [80].

Broadly speaking, hydrodynamics may be defined as the description of condensed states of matter on slow time scales and at large length scales. Macroscopic behaviour is characterized by the dynamics of a small number of slow fields (so-called “hydrodynamic fields”), related to conservation laws and broken continuous symmetries [74]. On time scales long compared to the fast relaxation times of microscopic variables, the assumption of local thermodynamic equilibrium leads to the definition of a thermodynamic potential as a function of all relevant (long-lived) thermodynamic variables and their conjugate quantities. Standard manipulations lead to the expression of the entropy creation rate as a bilinear functional of generalized fluxes and forces. In the vicinity of thermodynamical equilibrium, generalized fluxes are expressed as linear combinations of generalized forces. Due to microreversibility, the Onsager symmetry theorem implies that cross-coefficients must be set equal (respectively opposite) when fluxes and forces have equal (respectively opposite) sign under time reversal [81]; [82].

The hydrodynamic formalism is physically intuitive. It is flexible and can accommodate a broad spectrum of physical quantities, as long as deviations from equilibrium can be linearised. This approach is quite general since it relies on thermodynamic principles and on the invariance properties of the problem under consideration. Constitutive equations may thus be written within the domain of linear response as linear relationships between generalized fluxes and forces.

#### Dissipation function formalism

However, not all materials are described by a linear force-flux relationship. In tissues, plastic events such as cell rearrangements (Fig. 2) have thresholds which break down the linearity, and hydrodynamics (Section II.2.2) becomes inadequate. Deriving constitutive equations requires a more general formalism.

In the dissipation function formalism, the state of the material is described by the total deformation and by additional, independent, internal variables , with . These variables may be scalar, vectorial or tensorial. So-called “generalized standard materials” are defined by the existence of the energy function and the dissipation function , which are continuous (not necessarily differentiable) and convex functions of their respective arguments [83]; [84]; [85]:

(3) | |||||

(4) |

Here denotes the total deformation rate, and is the (Lagrangian) time derivative of . Although it is rarely explicitely stated, these energy and dissipation functions should increase when the norms of their arguments tend to infinity, so that they admit one (and only one) minimum, reached for a finite value of their arguments. Constitutive and evolution equations are obtained through the following rules, where denotes the stress:

(5) | |||||

(6) |

Appendix B provides more details on how to use the dissipation function formalism. In particular, Appendix B.1 explains how to manipulate the corresponding equations. Appendix B.2 explicitly treats the tensorial case, and shows that the tensorial variables in Eqs. (3,4) should be decomposed into their trace and deviatoric parts:

(7) | |||||

(8) |

considered as independent variables. Appendix B.3 explicits the incompressible case.

The formalism of Eqs. (3-6) is a convenient tool for building complex models and obtaining in a systematic way the full set of partial differential equations from simple and comprehensive graphical schemes. The coupling coefficients arise as cross partial derivatives, with the advantage that they derive from a smaller number of free parameters than in the hydrodynamics formalism.

For purely mechanical diagrams made of springs, dashpots and sliders, the dissipation function formalism yields the same equations as when directly writing the dynamical equations from the rheological model, as shown in Appendix B.4. The same formalism also applies to systems with non-mechanical variables, see Section III.4 and Appendix C.

A direct link between
the hydrodynamic and dissipation function formalisms can be established
when the dissipation function is a quadratic function of its arguments.
For a given variable , quadratic terms
in the energy function or in the dissipation function are harmonic,
*i.e.* they yield a linear term
or in the derived dynamical equations,
exactly like in the hydrodynamics formalism (Section II.2.2).
In this case, is proportional
to the rate of entropy production [82]
(this is also true for viscoplastic flows [86]).

The dissipation function formalism is also suitable for non-linear terms. Terms of the form or , with ( integer or real) yield non-linear terms or in the dynamical equations.

Interestingly, the dissipation function can even include the particular case which corresponds to terms like or . This yields terms of the form or in the derived equations: these are non-linear terms which dominate over the linear ones. This lowest-order case is useful to include plasticity (see Section III.1 for an example), whose treatment thus becomes straightforward [84]; [85]; [86]: the dissipation function formalism has been successfully applied to viscoelastoplastic flows [87]; [88]. As discussed in Section III.1, regularizing such terms would suppress yield stress effects.

Since the cross coupling coefficients arise as cross partial derivatives, they are by construction always equal by pairs. The Onsager symmetry theorem [81] is thus immediately obeyed when fluxes and forces behave similarly under time reversal, but not if they behave differently. Note that this is compatible with the constitutive equations of living tissues: since microreversibility may not apply at the (cell) microscale, the Onsager symmetry theorem does not need to apply.

Active ingredients which impose a force, a deformation rate, or a combination thereof, can be included in the dissipation function formalism, as shown for instance in Section III.3. The functions and remain convex and still reach a minimum for a finite value of their arguments. As expected, the entropy creation rate is no longer always positive.

From a mathematical point of view, the large set of nonlinear differential equations is known to be well posed in the Eulerian and small deformation setting [89]. Since the free energy function and the dissipation function are both convex, in the case of small deformations the existence and uniqueness of solution is guaranteed, while the second law of thermodynamics is automatically satisfied [83]; [85]. This is a major advantage of the dissipation function framework.

In addition, the formalism is also effective from a computational point of view. For problems which involve multi-dimensional and complex geometries together with large deformations of the tissue, there is no hope to obtain an explicit expression of the solution: its computation should be obtained by an approximation procedure. The resolution of the large set of nonlinear differential equations and its convergence at high accuracy require both a dedicated algorithm and a large computing time with the present computers. The convexity of and functions enables to use robust optimization algorithms to solve efficiently the set of dynamical equations thanks to variational formulations [90]; [91]. This second major advantage of the dissipation function framework has been widely used in small deformation, for applications in solid mechanics and plasticity, and has allowed the development of robust rocks and soils finite element modeling softwares (see e.g. [86]; [89] and references therein).

In summary, since the dissipation function formalism allows to treat plasticity and is convenient for numerical resolution, we recommend to adopt it for living tissues.

### ii.3 Specificity of cellular material modelling

While continuum mechanics is standard, cellular material modelling requires care on specific points. They include: the separation of the deformation between its contribution arising from inside each cell and from the mutual cell arrangement (Section II.3.1); the choice of Eulerian rather than Lagrangian description for a viscous, elastic, plastic material (Section II.3.2); and the treatment of large elastic deformations (Section II.3.3).

#### Intra-cell and inter-cell deformation

Different deformation rates can be measured simultaneously and independently (see Appendices A.1 and D.1).

The total deformation rate can be measured by tracking the movements of markers, moving with the tissue as if they were pins attached to the tissue matter. This total deformation rate originates from the following two contributions at the cellular level.

The intra-cellular deformation rate is the average of the deformation rate as perceived by individual cells, where each cell is only aware of the relative positions of its neighbours. The intra-cellular deformation can be measured by observing the anisotropy of a group of tracers attached to a reference cell and its neighbours, followed by an average over reference cells. By contrast with , the tracers are not attached permanently to the tissue itself: when neighbours rearrange and lose contact with the reference cell, the corresponding tracers are switched immediately to the new neighbours. The intra-cell deformation rate is then obtained as the rate of change of this intra-cell deformation measure .

The inter-cellular deformation rate reflects the cell rearrangements and relative movements. It can be measured by tracking the rearrangements themselves.

In tissues made of cells which tile the space, the stress at the tissue scale is the stress carried by the cells themselves (like in foams, but as opposed e.g. to the case of plants, where the rigid walls are as important as pressure for stress transmission). We thus advocate a decomposition in series, where intra- and inter-cellular stresses are equal, while intra- and inter-cellular deformation rates add up (Fig. 3):

(9) |

Choosing this decomposition into intra-cell and inter-cell contributions in series has consequences on the arguments of the dissipation function. When defining the variables, see Eqs. (3,4), it is relevant to choose one of them as equal to . Appendix D.1.2 discusses the case of large deformations.

#### An Eulerian rather than a Lagrangian approach

We now compare the Lagrangian and Eulerian points of view, and explain why we choose the latter.

For materials that retain information about their initial state, it is natural (and common) to use a deformation variable, often denoted , that compares the current local material state to the initial state of the same material region. This is called the Lagrangian description and is usually preferred for elastic solids [92]. The deformation rate of the Lagrangian description is the time derivative of the deformation , where the dot denotes the material derivative used for small deformations, :

(10) |

On the other hand, when plastic or viscous flows erase most or all memory of past configurations, it is common practice to use only the current velocity field as the main variable, with no reference to any initial state. This is called the Eulerian description, used for instance when writing the Navier-Stokes equations.

Both descriptions are tightly connected: the deformation rate of the Lagrangian description is equal to the (symmetric part of the) gradient of the velocity field of the Eulerian description:

(11) |

See Eqs. (189,191) for the complete expression at large deformations.

The separation of intra- and inter-cellular dynamics, discussed in Section II.3.1, can be viewed as mixing the Lagrangian and Eulerian points of view. Since cells retain their integrity, the quantity is similar to a deformation variable in a Lagrangian approach. By contrast, the relative motions of cells (described by ) is similar to the relative motions of material points in usual fluids, described from an Eulerian point of view.

A globally Eulerian description has been implemented for liquid foams in a direct manner, based on the argument that the rearrangement deformation rate progressively erases from the material the memory of the initial configuration and thus progressively wipes, like in common fluids, the relevance of the material deformation for predicting the future evolution of the material [69].

Here, the same argument should apply. In the examples provided in Section III the variables and naturally disappear from the final constitutive equation, and only their corresponding deformation rates and contribute. This reflects the absence of any structure holding cells together beyond the first neighbours. It confirms the relevance of an Eulerian description for a material such as a tissue.

#### Large elastic deformations

For pedagogical reasons, in the present article, all equations are written within the limit of small elastic deformations. Yet, in living tissues, large elastic deformations are encountered.

Appendix D explains in details how to model large deformations in the specific context of tissue mechanics, and reformulate accordingly the dissipation function formalism. In particular, it discusses the volume evolution, the elasticity and its transport, and the intra-cell deformation.

An example of changes due to large deformations is the distinction between two quantities which are equal in the limit of small deformations (Eq. (11)): the deformation rate and the symmetrised velocity gradient (Eqs. (189,191)). The transport of large elastic deformations involves objective derivatives. Several such derivatives exist, for instance lower- and upper-convected derivatives, as well as Gordon-Schowalter derivatives which interpolate between them. In rheological studies of complex fluids, the selection of the derivative is often motivated by formal reasons, or is empirical. In Appendix D.1.1, for physical reasons, we describe the deformation using tensors attached to the cellular structure; we show that this choice selects univocally the upper-convected derivative.

### ii.4 Set of partial differential equations

A tissue may be spatially heterogeneous: its material properties, its history, its interaction with its environment are under genetic control and may depend on the position . For instance the tissue may comprise different cell types, or it may be placed on a spatially modulated substrate.

The parameters and variables which describe the tissue are fields that may vary spatially. Here, we consider only tissues amenable to a continuum mechanics description, namely tissues whose relevant fields are smooth and slowly variable over the scale of a group of cell (the “representative volume element” of the continuum mechanics description). In what follows, we assume that the fields are continuous and differentiable.

The evolution of the tissue is then expressed as a set of partial differential equations (PDE), consisting in conservations laws and constitutive equations. To make this article self-contained, we show in the present Section how constitutive relations, such as derived using the dissipation function formalism, can be embedded in the rigorous framework of continuum mechanics in order to obtain a closed set of evolution equations.

In continuum mechanics, one usually starts from the conservation equations of mass, momentum and angular momentum. The mass conservation equation reads:

(12) |

where is the mass density (or mass per unit area in 2D); and represents material sources or sinks which, in the context of a tissue, can be linked with cell growth and apoptosis, respectively, see Section III.2.2.

In general, the conservation of momentum reads:

(13) |

which relates the acceleration to the internal stress tensor and the external forces . For instance, the external force may contain a friction component [26]; [38]; [41]. Note that, in a tissue, the inertial term is generally negligible when compared to the stress term . The validity of this approximation has to be checked in specific examples by estimating the value of the relevant dimensionless number, e.g. the Reynolds number for a purely viscous material, or the elastic number for a purely elastic solid. In this case, the conservation of momentum (Eq. (13)) reads:

(14) |

Finally, the conservation of angular momentum implies that the stress tensor is symmetric [74]: .

We obtain a set of evolution equations (Eqs. (5,6,11, 12,14)). There are unknown fields: , , , and . For any value of the space dimension, the number of coordinates of the unknown fields always equals the number of equations. This set of partial differential equations is closed by suitable initial conditions on the same variables and boundary conditions in terms of velocity, deformation and/or stress components. Its solution can be estimated by numerical resolution: see e.g. [70]; [90]; [91]; [93]; [94] for such numerical methods in the context of liquid foam flows.

## Iii Ingredients included in tissue modelling

A (non-exhaustive) list of ingredients for tissue modelling includes viscosity, elasticity, plasticity, growth, contractility, chemical concentration fields, cell polarity, and their feedbacks. Note that other tissue-specific ingredients such as (possibly active) boundary conditions [41] do not contribute to constitutive equations themselves: they are used to solve the set of partial differential equations (established in Section II.4).

In Sections III.1 to III.4, we present four worked out examples showing how such ingredients are taken into account within the dissipation function formalism. Each choice in this Section is motivated by the simplicity (rather than by the formalism, as in Section II, or by the realism, as in Section IV). Section III.5 combines individual ingredients into a composite rheological model by classifying them in terms of shape or volume contributions at the intra-cell or inter-cell level, and derives the corresponding set of equations.

### iii.1 Plasticity

#### Rearrangements and plastic deformation rate

Recent experiments performed on cell aggregates and cell monolayers have shown that these tissues can have a yield stress [39]; [45] and display a plastic behavior [10]; [95]. The origin of this plasticity includes cell rearrangements (Fig. 2) [96]; [97] which also play an important role during development, as in e.g. convergence-extension [15].

At the cell scale, and independently of its biological origin and regulation, a cell rearrangement is mathematically speaking a discontinuous topological process in a group of neighbouring cells. The associated mechanical description is decomposed into several steps. Before the rearrangement, the cells deform viscoelastically. During the rearrangement, two cells in contact get separated. Both other cells establish a new contact. They all eventually relax towards a new configuration with the consequence that two cells have got closer while the others have moved apart.

Upon coarse-graining spatially at a scale of several cells, and temporally over a time scale much larger than the relaxation time, the discontinuities at the cell scale are wiped out. The net result is an irreversible change in the stress-free configuration of the tissue, with convergence along one axis and extension along the perpendicular one. It is thus best described as a tensor with positive and negative eigenvalues [68]; [69], which tends towards the plastic deformation rate in the continuum limit. This tensor is the difference between the total deformation rate and the elastic deformation rate , so that the cumulated effect of elastic and plastic contributions add up:

(15) |

#### Plasticity and dissipation function

In the dissipation function formalism, elastoplastic and viscoelastoplastic materials are classically described by adding in the dissipation energy a yield stress term, for instance proportional to the norm of the deformation rate (Section II.2.3). Such non-analytic term is fully compatible with the formalism. The convexity of the energy function is preserved, so that equations are readily written and can be solved using known numerical approaches. For details, see Ref. [98].

When writing equations, physicists often favor smooth (analytic) expressions rather than discontinuous (singular) ones. Formally, it would of course be possible to regularize the plasticity equations and obtain a differentiable dissipation function by replacing each non-analytical term with analytic, strong non-linearities. However, this would lead to a completely different category of models, from which yield stress effects are absent and the solid behaviour vanishes in the long time limit.

#### A viscoelastoplastic example

We treat explicitly an example obtained by adding an elastic element of modulus in series with a diagram representing a Bingham fluid, namely the combination in parallel of a dashpot of viscosity and a slider of yield stress (Fig. 4). A more realistic (and therefore more complex) rheological model of a tissue that includes plasticity is treated in Section IV.1.

According to Fig. 4, we have:

(16) |

Choosing and as independent variables, the energy function reads:

(17) |

and the dissipation function:

(18) |

(19) | |||||

(20) |

which together yields the constitutive equation:

(21) |

When , takes a value in the interval : for a rigorous mathematical analysis, see Ref. [88], and in particular its Eqs. (9,10).

### iii.2 Growth

#### Conservation equations

The growth of a tissue has two aspects: mass growth and cell concentration growth. The mass growth affects the mass conservation (Eq. (12)) through its source term , where is the rate of variation of mass density, which has the dimension of an inverse time:

(22) |

Similarly, the cell concentration (number of cells per unit volume) evolves according to:

(23) |

where is the rate of variation of the cell concentration.

If we consider a tissue without gaps between cells, it is reasonable to write that the mass density is a constant, which simplifies the description. Equation (22) reduces to:

(24) |

Eq. (115) indicates how to locally measure . For a monolayer of varying thickness but homogeneous and constant density , Eq. (22) can be rewritten using the two-dimensional divergence applied in the plane of the monolayer, :

(25) |

#### Cell growth modes

Several cell processes alter the tissue volume, the cell concentration, or both (Fig. 5). They are compatible with the dissipation function formalism, see for instance Section IV.2.1. All rates are noted and have the dimension of the inverse of a time (for instance typically a few hours or a day for the division rate of epithelial cells in vitro [99]). For each process, the corresponding is the proportion of cells that undergo the process per unit time. Each process becomes relevant as soon as the duration of an experiment is of the order of, or larger than, the inverse of its rate .

Cells may grow, i.e. swell, with rate , which creates volume and decreases . They may undergo cytokinesis, i.e. split into daughter cells, which increases with rate ( doubles in a time ) without altering the volume. In the case where the swelling and cytokinesis rates are equal, and their common value is the cell-cycle rate , the long-time average of cell size is constant. For simplicity, the description of cytokinesis proposed in this section is scalar. More generally, one may need to define a tensor , with (see Section III.5.1).

Cells may undergo necrosis (rate ), which does not alter the tissue volume but decreases the concentration of living cells. The description of apoptosis (rate ), namely cell death under genetic control, requires some care. If the content of the apoptosed cell material is eliminated (for instance diffuses away, or is cleaned by macrophages) without being taken up by the neighbouring cells; and if we further assume that the tissue remains connective (i.e., neighbouring cells move to span the emptied region); then apoptosis causes the tissue volume to decrease while remains unaltered.

The mass and cell concentration growth rates (defined in Section III.2.1) are:

(26) | |||||

(27) |

A first special case is the situation where cells undergo cytokinesis without growing, hence leading to a decrease of the average cell size, as encountered for example during the first rounds of cytokineses in a developing embryo. In that case, only the cell concentration has a non-zero source, and Eqs. (26,27) reduce to:

(28) | |||||

(29) |

In the situation that combines cell swelling and cytokinesis at equal rates, in equal amounts, and in the presence of apoptosis, Eqs. (26,27) reduce to:

(30) | |||||

(31) |

#### A one-dimensional example

We illustrate with a one-dimensional example the case where cell swelling and cytokinesis rates are equal, resulting in a constant long-time average cell size. Here, we assume that after a full cell cycle the rest length of each daughter cell, defined as the length where their elastic energy is lowest, is eventually identical to that of the mother cell. Appendix A.3 describes a tissue made of elastic cells which divide and/or die; it shows that the stress evolution equation reads:

(32) |

Eq. (32) expresses that as expected, the stress (counted positive when tensile), increases when the tissue is subjected to elongation and decreases when the growth rate increases the tissue rest length [100]; [101].

Eq. (32) can be derived from the rheological diagram shown on Fig. 6, as follows. Consider a motor working at constant deformation rate, here the growth rate , in series with a spring of stiffness . It creates a difference between the total deformation rate and the deformation rate of the spring:

(33) |

Combined with the elasticity:

(34) |

this yields (32).

### iii.3 Contractility

In tissues, cell-scale contractility is often determined by the distribution of molecular motors such as myosin II. Upon coarse-graining, this distribution translates into tissue-scale contractility. In one dimension, such contractility may be modelled by a constant stress, as done in classical models of muscle mechanics [102]. As an example, we study the rheological diagram of Fig. 7a, where within a Maxwell model (spring and dashpot in series) a contractile element is placed in parallel with the dashpot.

Choosing and for instance as independent variables (with ), the energy and dissipation functions can be written in the form of Eqs. (3,4), with :

(35) | |||||

(36) |

where denotes the active stress: it is positive in the case of a contracting tissue.

Eqs. (35,36) injected into Eqs. (5,6) yield:

(37) | |||||

(38) |

Differentiating Eq. (37) yields . Combining it with Eq. (38) yields the stress evolution equation:

(39) |

where is the viscoelastic relaxation time. Eq. (39) is the evolution equation of a classical Maxwell element modified by a constant shift in stress due to the active stress. This is reminiscent of the active force included in [38]. Note that the same equation also describes an active gel of polar filaments, as introduced in [103]. The active stress can of course be tensorial, for instance when the spatial distribution of motors is anisotropic. This can readily be taken into account by the formalism, here at tissue scale (analogous continuum descriptions have already been performed at the scale of the cytoskeleton [80]).

In the rheological diagram of Fig. 7b, represents a constant deformation rate (counted negative when the tissue contracts). Such a rheological diagram has been introduced at the sub-cellular length scale, in the context of the actin-myosin cortex [104]; the active deformation rate is then interpreted in terms of the myosin concentration , the step length of the molecular motors and the binding rate , yielding: .

Strikingly, the rheological diagram of Fig. 7b leads to the same stress evolution equation as the diagram of Fig. 7a. This can be directly checked by decomposing the total deformation rate as

(40) |

and by defining the energy and dissipation functions again in the form of Eqs. (3,4) with :

(41) | |||||

(42) |

Injecting Eqs. (41,42) into Eqs. (5,6) and differentiating with respect to time yields:

(43) |

Eq. (43) is the same as Eq. (39) provided that is replaced with . Both rheological diagrams of Fig. 7 are thus equivalent. This non-uniqueness also exists in rheological diagrams with only passive elements, see Section II.2.1.

### iii.4 Coupling non-mechanical fields to a rheological model

Suppose we need to describe an additional field which is non-mechanical and cannot be included in rheological diagrams. The dissipation function formalism, which allows to postulate forms of the energy and dissipation functions that respect the symmetries of the problem, provides a framework within which various couplings between fields may be introduced in a systematic manner.

A consistent continuum modelling of a cellular material sometimes requires to include tensors, for instance as variables of the energy and dissipation functions. Within a variational framework, Sonnet and coworkers [105]; [106] have performed a detailed study of the derivation of constitutive equations involving a tensorial order parameter. In epithelia, an example is given by noting that planar cell polarity proteins exhibit tissue-scale ordered domains that are often best described by a tensor field [22]; [23].

Inspired by this last example, we treat the case of a viscous liquid (Fig. 8), and choose to couple the deformation rate tensor to a second-order tensor in the dissipation function. We could have chosen a non-mechanical field which is scalar (Appendix C.1), or vectorial: polar, for a usual oriented vector (Appendix C.2), or axial, for a nematic-like vector. Of course, more complex couplings may be considered whenever needed, that also involve other ingredients such as tissue growth (Sections III.2 and IV.2) or cell contractility (Section III.3).

Formally, Eqs. (3,4) should be written with an additional internal variable , so that , and with tensorial coupling parameters: , . Since trace and deviators of tensors have complementary symmetries, it is convenient to treat them as separate variables with distinct, scalar coupling parameters (see Appendix B.2); Eqs. (3,4) thus read, with :

(44) | |||

(45) |

where the colon denotes the double contracted product between tensors :

The parameters , , , , , are non-negative and the inequalities , ensure the convexity of the dissipation function. From Eq. (5) we first compute the stress tensor:

(46) | |||||

(47) |

where a linear coupling to modifies the usual constitutive equation of a viscous liquid. From Eq. (6), we next obtain the evolution equations:

(48) | |||||

(49) | |||||

which yield the evolution equations for the tensor :

(50) | |||||

(51) |

Note that the relaxation times for the deviator and trace, respectively and , can in principle be different. Inserting Eqs. (50,51) into Eqs. (46,47) yields the stress tensor:

(52) | |||||

(53) |

In the long time limit, the tensor tends to:

(54) | |||||

(55) |

and the viscous stress tensor tends to:

(56) | |||||

(57) |

As an example, the spatial distribution of a myosin (called Dachs) in the dorsal thorax of fruitfly pupae was studied quantitatively in [23]. Fluorescence microscopy images reveal a tissue-scale organization of Dachs along lines, allowing to measure an orientation and an amplitude, from which a deviatoric tensor is defined. Dachs orientation correlates at the tissue scale with the direction of contraction quantified by the corresponding eigenvector of the velocity gradient tensor, as predicted by Eq. (54).

### iii.5 Combining ingredients

Tissue volume | Tissue shape | ||||

intra | inter | intra | inter | ||

cell | cell | cell | cell | ||

vol. | num. | sha. | pos. | ||

cell swelling | + | ||||

apoptosis | - | ||||

cytokinesis | - | + | + | ||

contractility | + | ||||

rearrangements | + |

The respective contributions of plasticity, growth and contractility to the rate of change in elastic deformation have been expressed by Eqs. (15,33,40).

To model a given experiment, the relevant ingredients, for instance those listed in the introduction of Section III or in Fig. 5, can be assembled at will. Section III.5.1 suggests how to classify them into intra- and inter-cell ingredients and Section III.5.2 presents an example of such a combination.

#### Classification of ingredients

The distinction between “intra-cell” and “inter-cell” contributions can be complemented by a distinction between contributions which alter the shape and/or the volume of the tissue. Such distinction helps understanding the biological meaning of equations. For tensorial ingredients, Eq. (9) becomes:

(58) | |||||

(59) |

This classification (Table 1 and Fig. 9) is merely indicative and should be adapted for any specific tissue under consideration according to the available knowledge or intuition. For instance, Table 1 assumes (see ) that contractility does not change the actual volume of each cell, whether in a 3D tissue or in an epithelium, but that it may change the apparent surface area of cells in an epithelium.

Let us review some ingredients expected to contribute to the four parts of the total deformation rate tensor (Eqs. (58), (59)).

The rate of change in the cell volume can be written in terms of the isotropic elastic deformation, the cell swelling rate and the cytokinesis rate:

(60) |

The number of cells increases due to cytokinesis and decreases due to apoptosis:

(61) |

The cell shape deformation rate can be expressed in terms of the deviatoric elastic deformation, and the active contractility rate:

(62) |

Finally, the arrangement of cell positions is affected by cytokinesis and by the cell rearrangement contribution to plasticity, which is purely deviatoric:

(63) |

Each ingredient listed above then provides a term either in the energy or in the dissipation function , except motor elements (Sections III.2.3 and III.3) which correspond to .

Within the dissipation function framework (Section II.2.3), Eqs. (58-63) play the role of the topological relations between deformation rate variables, see Eq. (16). Combined together, Eqs. (58-63) enable to split the evolution of the elastic deformation (Eq. (33)) into the following two equations:

(64) | |||||

(65) |

#### A complex example

To open the way towards more realistic, complex descriptions, we now present an example of a tissue rheology model that incorporates most ingredients listed in the introductions of Sections III and III.5. Figure 10 is decomposed into four blocks. The upper line corresponds to the deviatoric part of the deformation, and the lower to the volume-related rheology. Each part is further decomposed into intra-cell rheology (left block) and inter-cell processes (right block).

Fig. 10 indicates the topology of the diagram, and the numerous parameters involved: , , , , , , , , , , . It yields the energy and dissipation functions:

(66) | |||||

(67) | |||||

where the following independent variables have been chosen (see Section II.2.3 and Appendix B.1): and for the springs, for the dashpot , for the cell volume, and as usual and for the total deformation rate.

## Iv Applications to the mechanics of cell aggregates

In the present Section, we combine ingredients introduced in Section III to write and solve the dynamical equations in two more realistic examples. These examples are inspired by the rheology of cellular aggregates, first when deformed on a timescale short compared to the typical cell cycle time (Section IV.1), second when growing between fixed walls on a timescale long compared to (Sections IV.2 and IV.3). In both cases, we separate the contributions of intra- and inter-cellular processes to aggregate rheology. Both examples derive from a rheological diagram, but more complex situations, including for instance non-mechanical fields as in Section III.4, can be solved within the same formalism.

### iv.1 Without divisions: creep response

Although an actual cell aggregate is complex, we crudely model it by combining an intra-cellular viscoelasticity and an inter-cellular plasticity (Fig. 11): one aim of the present Section is to illustrate their interplay. For simplicity, we neglect aggregate volume changes, and use scalar variables; we defer a tensorial treatment to Section IV.2.

Each cell is modelled as a fixed amount of viscous liquid enclosed in an elastic membrane which prevents it from flowing indefinitely at long times. We thus model the cell as a viscoelastic solid, more precisely a Kelvin-Voigt element: a spring which reflects the effective cell shear elastic modulus (typically of order of the cell cortex tension divided by the cell size) is in parallel with a dashpot which reflects the cytosplasm viscosity .

Cells need to undergo a finite deformation before triggering a rearrangement. When the stress exerted on the aggregate exceeds the yield stress , cells rearrange; the aggregate flows like a liquid with a viscosity much larger than [11]. If the aggregate presents a high level of active cell contour fluctuations, whatever their origin, the cells undergo rearrangements more easily: in practice, this biological activity will lower the yield stress and the viscosity [10].

For simplicity, we turn to a one-dimensional description, and assume that . Their time-derivatives and correspond to the respective 1D projections of the plastic deformation rate and of the elastic deformation rate .

From Fig. 11, the energy and dissipation functions of the independent variables and read:

(68) | |||||

(69) | |||||

where vanishes when . From Eqs. (5,6) we obtain:

(70) | |||||

(71) | |||||

We can now turn to predictions. For instance, in a “creep experiment”, the total stress is zero until , then is set at a constant value during the time interval from to . Figure 12 represents the corresponding creep curves, i.e. the time evolution of the total deformation , obtained as the analytical solution of Eqs. (70,71).

If the magnitude of the applied stress is lower than the yield stress , the deformation reaches a plateau value equal to . When is brought back to zero, the deformation relaxes back to zero over a time : this viscoelastic time is a natural timescale of the material and reflects the individual cell rheology.

If is larger than the yield stress , after a typical time cell shapes reach their maximal deformation, so that the aggregate thereafter flows only as a result of cell rearrangements: the aggregate deformation increases steadily at a rate . When is brought back to zero, the cell shapes relax to equilibrium within a time . The total deformation correspondingly relaxes, yet not back to zero.

The creep curves shown in Fig. 12 are similar to quantitative measurements performed during the micropipette aspiration of cell aggregates which were assumed to be viscoelastic [63]. A similar yield behavior was observed when stretching a suspended cell monolayer: the monolayer deformation reached a plateau at low applied stress, while a creep behavior appeared at higher stress [39]. The authors observe no divisions or rearrangements in the course of deformations that reach circa 70%, which suggests that in their experiment, the creep behavior arises from an intra-cell contribution.

### iv.2 With divisions: inhomogeneous proliferation

When the duration of the experiment becomes much longer than , cell divisions must be taken into account and the aggregate flows even under a weak stress [9].

Figure 13 represents an experiment [107] where a cell aggregate is confined between parallel rigid plates. The aggregate length grows from initially to in days and in . We thus estimate the growth rate as , where the factor reflects the fact that growth occurs effectively in only two dimensions: cells divide mostly at the aggregate periphery. The main hypothesis proposed by [107] is that the stress induced by growth in the confined aggregate could mechanically inhibit mitosis.

#### Model with divisions

We propose to qualitatively model the aggregate growth coupled to its mechanical response (determined in Section IV.1). Assuming translational invariance along , we treat this problem in the plane (). We introduce separate rheological diagrams for the trace and the deviator.

According to the mass conservation (Eq. (24)), the trace of the deformation rate is equal to the growth rate of the aggregate. Within a linear approximation [12], we assume that decreases with the pressure , from its value at zero pressure, as:

(72) |

where we use the pressure actively generated by the aggregate at zero deformation rate, e.g. when confined between fixed walls. We deduce the rheological equation for the trace (Fig. 14, top part):