# On continuum modeling of cell aggregation phenomena

Soheil Firooz<EMAIL_ADDRESS>Stefan Kaessmair  Vasily Zaburdaev  Ali
Javili  Paul Steinmann  Institute of Applied Mechanics, University of
Erlangen-Nuremberg, Egerland Str. 5, 91058 Erlangen, Germany  Siemens Industry
Software GmbH, Nordostpark 3, 90411 Nuremberg, Germany  Department of Biology,
University of Erlangen-Nuremberg, 91058 Erlangen, Germany  Max Planck Zentrum
für Physik und Medizin, 91058 Erlangen, Germany  Department of Mechanical
Engineering, Bilkent University, 06800 Ankara, Turkey  Glasgow Computational
Engineering Center, James Watt School of Engineering, University of Glasgow,
Glasgow G12 8QQ, United Kingdom

###### Abstract

Cellular aggregates play a significant role in the evolution of biological
systems such as tumor growth, tissue spreading, wound healing, and biofilm
formation. Analysis of such biological systems, in principle, includes
examining the interplay of cell-cell interactions together with the cell-
matrix interaction. These two interaction types mainly drive the dynamics of
cellular aggregates which is intrinsically out of equilibrium. Here we propose
a non-linear continuum mechanics formulation and the corresponding finite
element simulation framework to model the physics of cellular aggregate
formation. As an example, we focus in particular on the process of bacterial
colony formation as recently studied by Kuan et al. [1]. Thereby we describe
the aggregation process as an active phase separation phenomenon. We develop a
Lagrangian continuum description of the problem which yields a substantial
simplification to the formulations of the governing equations. Due to the
presence of spatial Hessian and Laplacian operators, a gradient-enhanced
approach is required to incorporate $\mathcal{C}^{1}$ continuity. In addition,
a robust and efficient finite element formulation of the problem is provided.
Taylor–Hood finite elements are utilized for the implementation to avoid
instabilities related to the LBB condition. Finally, through a set of
numerical examples, the influence of various parameters on the dynamics of the
cellular aggregate formation is investigated. Our proposed methodology
furnishes a general framework for the investigation of the rheology and non-
equilibrium dynamics of cellular aggregates.

###### keywords:

Cellular aggregates, Active phase separation, Continuum model, Eulerian
approach, Lagrangian approach

##  1 Introduction

Within the human body, most cells interact with their neighboring cells and
with their extracellular matrix to establish a unique organization. These
cell-cell and cell-matrix interactions form a complex network of mechanical,
biological and chemical signals which play a significant role in cell
physiology [2, 3, 4, 5]. The study of cellular interactions provides a
significant insight towards a better understanding of many biological
processes such as tumor growth [6, 7], tissue spreading [8, 9, 10, 11],
biofilm formation [12, 13, 14] and wound healing [15, 16]. Although there
exist numerous contributions on experimental methods to examine cellular
interactions and their final outcomes [17, 18, 19, 20, 21], there are still
certain details and information that remain hard to assess. Examples of such
details include the exact relation between the interactions and individual
cell’s adhesion properties, or the speed at which sorting between the cells
occurs. Clearly, determining the relationships between all the problem
variables via conducting experiments is an arduous and time consuming task.
Theoretical approaches, however, provide a valuable alternative as they do not
suffer from such limitations and facilitate parametric studies. Mathematical
modeling allows us to quantify parameters such as local speeds, fluxes and
their rates of change in a natural way. They also enable us to investigate the
effects of interactions between cells of different types or to allow the cell
properties to vary independently of one another. Given these benefits,
mathematical modeling proves to be an inclusive and efficient alternative for
analyses of cellular interactions [22].

There exist two major approaches for mathematical description of cellular
interactions: _agent-based_ and _continuum models_. In agent-based approaches,
the cells are modeled and treated individually using a set of biophysical
rules [23, 24, 25, 26]. Each cell is approximated, for example, as a
homogeneous, isotropic, elastic, spherical object parameterized by measurable
biophysical and biological quantities. The agent-based approach is
particularly useful when one wants to study the interaction of individual
cells with each other and with their environment [27, 28, 29]. Since this
method is based on a series of rules for each cell, translating biological
processes into a model is straightforward [30]. For small-scale studies or
cases in which the properties of the cells vary over distances comparable to
the size of a cell, a higher degree of spatial resolution is obtained via
agent-based models in comparison to continuum models. Despite all the
precision that the agent-based method offers, this approach is difficult to
study analytically and its computational cost greatly increases as the number
of cells increases. For instance, simulation of a tumor growth process
requires systems which evolve from a single progenitor cell to $10^{6}$ cells
_in vitro_ and $10^{11}$ cells _in vivo_. Carrying out a computational
analysis on such large cell population sizes is a cumbersome task, if possible
at all [31]. Additionally, it is often neither desirable nor necessary to
track each individual cell within a very large population. For larger-scale
applications, a continuum modeling proves to be a more viable alternative.
This approach is well suited to describe large scale phenomena where the cell
properties vary smoothly over a length scale of several cell diameters and
therefore the cell properties can be approximated by a local average.
Continuum models frequently involve ordinary and partial differential
equations which are usually in the reaction-diffusion form [32, 33, 34, 35,
36]. Many aspects of tumor and tissue growth have been studied using continuum
models [37, 38, 31]. A continuum description of cell motility due to cell-cell
and cell-matrix interaction was presented in [22, 39, 40]. They introduced a
non-local interaction term to account for adhesion between the cells and
between the cells and matrix. Coarse grained continuum approaches such as
hydrodynamic theories, have also provided a powerful tool to capture large
scale emergent behaviors in active cellular systems [41, 42, 43, 44, 45].
Further studies on continuum modeling of cellular aggregates are available in
[46, 47, 48]. We refer to [49] for a thorough comparison between agent-based
and continuum modeling of cellular aggregates.

A prototypical biological example of cellular aggregation is the formation of
bacterial microcolonies and biofilms. One of the first steps in the process of
bacterial colonization of biotic and abiotic surfaces is the formation of
aggregates or colonies consisting of several thousands of cells. Usually,
these microcolonies later evolve into much more complex bacterial communities,
known as biofilms [50, 51]. Bacterial infections involving biofilms are far
more resistant to anti-microbial treatments in many cases [52]. Thus,
investigation of the mechanism of bacterial microcolony formation is of
immense significance in the fields of medicine and engineering. A few well-
known examples of bacterial microcolonies causing dangerous microbial
infections are _Pseudomonas aeruginosa_ [53], _Neisseria meningitidis_ [54],
_Vibrio cholerae_ [55] and _Neisseria gonorrhoeae_ [56].

In this manuscript, we focus on microcolonies of _Neisseria gonorrhoeae (NG)_
bacteria. These bacterial microcolonies are the infectious units which form on
human epithelial tissue and cause gonorrhoeae, the second most common sexually
transmitted disease [57]. Multi-scale computational simulations have been
conducted recently to study biophysical aspects of NG microcolonies [58, 56,
59]. The NG bacteria, as well as many other bacteria species, use multiple
long and thin retractable filaments, called type IV pili, in order to interact
with the environment and with each other [60]. A series of studies have been
carried out to investigate the twitching motility of bacteria mediated by type
IV pili [61, 62, 63, 64, 65]. Pili can extend from the cell body, attach to
the substrate and retract. Pilus retraction generates forces which are then
translated into movement of cells. The magnitude of the forces generated by
the pilus retraction are in the range of $100-180\,\text{pN}$ which is
considered as one of the strongest active molecular forces known in nature
[66, 56]. Additionally, pili of one cell could also extend and attach to pili
of other cells. Retraction of the attached pili network attracts the cells
towards each other and leads to formation of an aggregate. These cycles of
growth, attachment, detachment and retraction drive the cell motility on
substrates and the aggregate formation process [67]. Pili mediated cell-cell
and cell-matrix interactions are crucial for the formation and maintenance of
microcolonies [68, 19, 21, 69].

The main objective of this contribution is to formulate and simulate the
process of cell aggregation phenomena within a nonlinear continuum mechanics
framework. We develop our framework in a Lagrangian setting which yields
considerable simplification of the equations and enables implicit time
integration which considerably increases the computational robustness. In
doing so, we take a prototypical example of NG bacteria and we describe the
process of colony formation as an active phase separation phenomenon. Our work
is mainly based on the coarse grained approach previously developed by Kuan et
al. [1, 70]. While we focus on intercellular interactions, our aim is to
provide a robust and efficient computational setting for generic cell-matrix
interaction problems, and to develop its fully nonlinear finite element
implementation. Our proposed framework provides a versatile and reliable
simulation technique that allows studying the processes of aggregate formation
under high forces and strong phase-separated regimes nearing much closer to
the physiologically relevant conditions.

Table 1: Summary of key definitions and notations. $\\{\bullet\\}$ |  an arbitrary quantity |  $\dot{\\{\bullet\\}}$ |  material time derivative of $\\{\bullet\\}$
---|---|---|---
$\nabla_{\boldsymbol{X}}\\{\bullet\\}$ |  material gradient of $\\{\bullet\\}$ |  $\nabla_{\boldsymbol{x}}\\{\bullet\\}$ |  spatial gradient of $\\{\bullet\\}$
$\nabla_{\boldsymbol{X}}\cdot\\{\bullet\\}$ |  material divergence of $\\{\bullet\\}$ |  $\nabla_{\boldsymbol{x}}\cdot\\{\bullet\\}$ |  spatial divergence of $\\{\bullet\\}$
$\Delta_{\boldsymbol{x}}\\{\bullet\\}$ |  spatial Laplacian of $\\{\bullet\\}$ |  $\nabla^{2}_{\boldsymbol{x}}\\{\bullet\\}$ |  spatial Hessian of $\\{\bullet\\}$
$\mathcal{L}_{t}\\{\bullet\\}$ |  Lie time derivative of $\\{\bullet\\}$ |  $\delta\\{\bullet\\}$ |  variation of $\\{\bullet\\}$
$R$ |  cell radius |  $E$ |  cell elastic modulus
$f^{\text{p}}$ |  pili-pili attractive force |  $f^{\text{s}}$ |  steric repulsive force
$n^{\text{p}}$ |  number of bound pili pairs |  $\xi$ |  cell-substrate friction coefficient
$k_{\text{on}}$ |  pili binding rate |  $k_{\text{off}}$ |  pili unbinding rate
$l$ |  pili length |  $\ell_{0}$ |  pili average length
$\boldsymbol{l}_{ij}$ |  distance vector between cells $i$ and $j$ |  $\boldsymbol{l}$ |  spatial velocity gradient
$\boldsymbol{X}$ |  material position vector |  $\boldsymbol{x}$ |  spatial position vector
$c_{0}$ |  material cell number density |  $c_{t}$ |  spatial cell number density
$p_{0}$ |  material bound pili number density |  $p_{t}$ |  spatial bound pili number density
$\boldsymbol{g}$ |  cell number density gradient |  $\boldsymbol{v}$ |  cell velocity
$\boldsymbol{y}$ |  non-linear deformation map |  $\boldsymbol{F}$ |  deformation gradient
$J$ |  Jacobian of the deformation gradient |  $\boldsymbol{K}$ |  cofactor of the deformation gradient
$\delta$ |  Kronecker delta |  $\Psi_{\text{tot}}$ |  total internal energy
$\boldsymbol{I}$ |  material second-order identity tensor |  $\boldsymbol{i}$ |  spatial second-order identity tensor
$\mathcal{B}_{0}$ |  material configuration |  $\partial\mathcal{B}_{0}$ |  boundary of the material configuration
$\boldsymbol{t}^{\text{a}}$ |  active traction on material configuration |  $\boldsymbol{t}^{\text{p}}$ |  passive traction on material configuration
$\delta\boldsymbol{y}$ |  linear momentum balance test function |  $\delta c$ |  cell number density conservation test function
$\delta\boldsymbol{g}$ |  cell density gradient continuity test function |  $\delta p$ |  bound pili number density evolution test function
$\boldsymbol{E}$ |  Green–Lagrange strain tensor |  $\boldsymbol{B}$ |  Piola deformation tensor
$\boldsymbol{S}^{\text{a}}$ |  active Piola–Kirchhoff stress |  $\boldsymbol{\tau}^{\text{a}}$ |  active Kirchhoff stress
$\boldsymbol{P}^{\text{a}}$ |  active Piola stress |  $\boldsymbol{\sigma}^{\text{a}}$ |  active Cauchy stress
$\boldsymbol{P}^{\text{p}}$ |  passive Piola stress |  $\boldsymbol{\sigma}^{\text{p}}$ |  passive Cauchy stress
$\boldsymbol{S}^{\text{f}}$ |  pili formation Piola–Kirchhoff stress |  $\boldsymbol{\sigma}^{\text{f}}$ |  pili formation Cauchy stress
$\overline{\boldsymbol{P}}$ |  cell density gradient continuity Piola stresses |  $\boldsymbol{N}$ |  material unit normal to the boundary
R |  assembled residual vector |  $\boldsymbol{U}$ |  global vector of unknowns
K |  assembled tangent stiffness |  $\\#e$ |  number of elements
$N$ |  shape function of the finite elements |  $\text{I}^{\text{sym}}$ |  symmetric fourth-order identity tensor
$\lambda$ |  penalty parameter |  $\Delta t$ |  time step

The remainder of this manuscript is organized as follows. Table 1 gathers the
key definitions and notations of the paper. Section 2 introduces the problem
definition and presents the governing equations. Finite element implementation
of the problem is elaborated in Section 3. Our proposed theory is illustrated
through a set of numerical examples in Section 4. Finally, Section 5 concludes
the work and provides further outlooks.

##  2 Governing equations

This section elaborates on the governing equations. First the problem of cell
aggregation is defined and all the parameters and their roles are introduced.
Afterwards, the continuum approach within the Lagrangian settings is detailed.
Note that for simplicity and readability, our formulations here are developed
for a two-dimensional case, but there are no conceptual limitations to
generalize it to three dimensions.

###  2.1 Problem definition

Figure 1:  Left: An image of N. gonorrhoeae bacterium and its pili obtained by
transmission electron microscopy in [1]. Right: Simplification of the cell
geometry to a circular shape for our analysis (right).  Figure 2:  A sketch of
N. gonorrhoeae bacteria interacting by their pili. The forces acting on each
bacterium are illustrated. Each bacterium $i$ is distinguished by its position
vector $\boldsymbol{r}_{i}$. The inter-bacterium forces are steric repulsion
forces between two attached bacteria which are shown in red, and the pili-pili
mediated attractive forces between two bacteria that have formed a pili
network are shown in blue. These two forces act in the same direction of the
line connecting the centers of two adjacent bacteria. Bound pili are depicted
by blue dotted lines whereas free pili are depicted by solid green lines. Note
that there could be multiple pili pairs pulling two neighboring cells
together, as counted by $n_{is}^{\text{p}}$.

Figure 1 (left) shows a transmission electron microscopy image of a single NG
bacterium together with its pili. In this contribution, we approximate the
cell geometry with a circular disk with the radius $R$ for the sake of
simplicity, as shown in Fig. 1 (right). Each NG bacterium is surrounded by
approximately $10-20$ pili which are isotropically distributed around the cell
[67, 64]. The length of each individual pilus $l$ was shown to be
exponentially distributed with the average value of $\langle
l\rangle=\ell_{0}$. Measurements report the average length $\ell_{0}$ to be
around $1-2\,\mu\text{m}$ with cell radius $R$ being around $1\,\mu\text{m}$
[62].

Figure 2 depicts a group of cells interacting via their pili. In order to move
on a substrate, cells use their pili. Pili can grow, attach to the substrate
and retract, which produces the force required for the displacement. Pili-
substrate interactions play a significant role in determining the cell
motility which is essential for understanding the kinetics of aggregation.
Here, since we are mainly interested in examining the behavior of aggregates,
only a substrate friction is considered to represent the pili-substrate
interactions. Further details on pili-substrate interactions are available in
[71]. The growth and retraction of the pili occurs through the process of
polymerization and depolymerization which is powered by specific motor protein
complexes [72, 73, 74]. Apart from attaching to the substrate, pili can also
attach to pili of other cells. Retraction of an attached pili pair generates
an attractive force $f^{\text{p}}$ which pulls the cells towards each other.
In Fig. 2, bound pili are depicted by blue lines with a dot on them, whereas
free pili are shown by solid green lines. The number of the pili pairs that
have formed between two cells is denoted as $n^{\text{p}}$. Accordingly, as
illustrated for the cells $i$ and $s$, the overall attractive force between
two cells which have formed a network of bound pili is
$f^{\text{p}}n^{\text{p}}_{is}$. Assuming no bound pairs at the initial time,
the number of bound pili between the cells $i$ and $j$ (that do not move
relative to each other) can be obtained via the relation

$\displaystyle\frac{\mbox{d}n_{ij}^{\text{p}}}{\mbox{d}t}=\displaystyle\frac{k_{\text{on}}\,e^{-l_{ij}\ell_{0}}}{2\pi\ell_{0}^{2}}-k_{\text{off}}n_{ij}^{\text{p}}\qquad\Longrightarrow\qquad
n_{ij}^{\text{p}}=\displaystyle\frac{k_{\text{on}}\,e^{-l_{ij}\ell_{0}}}{2\pi\ell_{0}^{2}k_{\text{off}}}\left[1-e^{-k_{\text{off}}\,t}\right]\,,$
(1)

with $k_{\text{on}}$ being the pili binding rate and $k_{\text{off}}$ being
the pili unbinding/detachment rate. The distance between the cells $i$ and $j$
is $l_{ij}$ which is the magnitude of the vector pointing from cell $j$ to
cell $i$ as
$l_{ij}=||\boldsymbol{l}_{ij}||=||\boldsymbol{r}_{i}-\boldsymbol{r}_{j}||$.
This relation was obtained by Kuan et al. [70, 1] which is a mean field
approximation that ignores the discreteness of the pili number. The factor
$\displaystyle\frac{k_{\text{on}}\,e^{-l_{ij}/\ell_{0}}}{2\pi\ell_{0}^{2}}$
stems from the assumption of an exponential distribution of the pili length
and integration over all possible binding points along the line connecting two
bacteria. Figure 3 renders the variation of bound pili pairs between the cells
$i$ and $j$ with respect to their distance $l_{ij}$ and time for three
different values of binding rate $k_{\text{on}}$ (two different views are
given for better illustration). The average pili length is set to
$\ell_{0}=1.0\,\mu\text{m}$. It is confirmed that as time evolves, the number
of bound pili pairs reaches a steady state. Larger distance between the cells
yields a decreased number of bound pili pairs. This is justifiable since the
further the cells are from each other, the less chance their pili have to
attach to each other due to their limited length. Evidently, larger binding
rates also result in more bound pili pairs.

Figure 3:  Variation of the bound pili pairs $n_{ij}^{\text{p}}$ between the
cells $i$ and $j$ with respect the their distance $l_{ji}$ and time $t$
according to relation (1). Three different values $0.1$, $0.2$ and $0.3$ for
the pili binding rate are considered. Two different views are provided for
better illustration. Pili average length is set to $\ell_{0}=1.0\,\mu\text{m}$
and the unbinding rate is set to $k_{\text{off}}=0.005\,\text{s}^{-1}$.

Another force that plays a role in the cell aggregate formation dynamics is
the steric repulsive force $f^{\text{s}}_{ij}$ which occurs between the cells
$i$ and $j$ that are in direct contact. In the model, having this repulsive
force is necessary to prevent interpenetration of the cells. In summary, three
major forces determine the dynamics of the cells in the network; cell-
substrate friction, pili-pili mediated attractive force and steric repulsive
force. The forces arising from the contraction of the constantly remodeling
pili network attract the cells together and in the presence of the excluded
volume interactions, a densely packed colony is formed which behaves as an
active visco-elastic material [1].

Since the pili-pili interactions are driven by active retractions and growth
of pili, the system is inherently out of equilibrium. The force balance
equation for an individual cell in the aggregate can be obtained as

$f^{\text{p}}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}n_{ij}^{\text{p}}\hat{\boldsymbol{l}}_{ji}+\sum_{\begin{subarray}{c}j=1\\\
j\neq
i\end{subarray}}f^{\text{s}}_{ij}\hat{\boldsymbol{l}}_{ij}-\xi\boldsymbol{v}_{i}=0\,,$
(2)

with $\boldsymbol{v}_{i}$ being the velocity of cell $i$ and $\xi$ being the
friction coefficient between the cell and the substrate. Note the first
summation is over the neighboring cells with which the cell has formed pili
pairs and the second summation is over the neighboring cells with which the
cell is in direct contact. The unit vector
$\hat{\boldsymbol{l}}_{ij}=\boldsymbol{l}_{ij}/l_{ij}$ shows the direction of
the pili-pili mediated force acting on cell $i$ by cell $j$ which is in the
same direction of the steric repulsive force between these two cells. The
three terms in Eq. (2) represent pili-pili mediated forces, excluded volume
interactions and cell-substrate friction, respectively. The summations in Eq.
(2) imply that the equation is solved for each cell individually. While this
approach offers more accuracy together with a higher degree of spatial
resolution, it suffers from computational costs as the number of cells
increases. Moreover, if one is interested in examining the collective/overall
behavior of the aggregate, tracking each individual cell is neither desirable
nor necessary. In this section, we propose a continuum model to describe the
dynamics of dense bacterial colonies. In our methodology, the system is
characterized by two key features which distinguish it from previous models of
active systems. First, at time scales smaller than the pili detachment time,
the bound pili network endows aggregates with elastic-like material properties
for which a Lagrangian approach might be more suitable. At larger time scales,
pili can rearrange which allows for stress relaxation resulting in a fluid-
like behavior for which an Eulerian approach is the common choice. Second, the
attractive pili-pili mediated force dipoles are balanced by steric repulsion
forces which allow for considering dense cellular aggregates. These two key
features are captured by the proposed continuum equations.

###  2.2 Eulerian approach

In the Eulerian approach, the cell aggregate is treated as a compressible
fluid-like material. Therefore, our domain would be a fixed control window and
we observe the changes that the aggregate goes through over this window. This
section briefly details on the governing equations of the problem within an
Eulerian approach that were previously derived by coarse-graining microscopic
equations in [1].

####  2.2.1 Balance equations

The first balance equation is the overall cell number conservation equation
which reads

$\displaystyle\frac{\partial c_{t}}{\partial
t}+\nabla_{\\!\boldsymbol{x}}\\!\cdot(c_{t}\boldsymbol{v})=0\,,$ (3)

where $c_{t}$ is the spatial cell _number_ density and $\boldsymbol{v}$ is the
cell (ether) velocity. Note, the subscript $\boldsymbol{x}$ indicates that
spatial derivatives are taken with respect to the coordinates in the current
configuration. The two variables $c_{t}$ and $\boldsymbol{v}$ are the unknowns
of our problem. The second balance equation is the linear momentum balance
which arises from the microscopic balance of forces in Eq. (2) and reads

$\nabla_{\\!\boldsymbol{x}}\\!\cdot\boldsymbol{\sigma}^{\text{a}}+\nabla_{\\!\boldsymbol{x}}\\!\cdot\boldsymbol{\sigma}^{\text{p}}-\xi
c_{t}\boldsymbol{v}=\boldsymbol{0}\,,$ (4)

with $\boldsymbol{\sigma}^{\text{a}}$ and $\boldsymbol{\sigma}^{\text{p}}$
being the active and passive Cauchy stresses, respectively. The third balance
equation is the bound pili number density evolution equation which reads

$\displaystyle\frac{\partial p_{t}}{\partial
t}+\nabla_{\\!\boldsymbol{x}}\\!\cdot\left(p_{t}\boldsymbol{v}\right)-k_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\,\left[\Delta_{\boldsymbol{x}}c_{t}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\nabla_{\\!\boldsymbol{x}}c_{t}\big{|}^{2}\right]+k_{\text{off}}p_{t}=0\,,$
(5)

with $p_{t}$ the spatial bound pili number density. This equation incorporates
the pili turnover dynamics from Eq. (1). Note, in our continuum formulation we
denote the spatial pili number density as $p_{t}$ which is a coarse-grained
version of $n_{ji}^{\text{p}}$. The spatial Laplacian operator
$\Delta_{\boldsymbol{x}}$ is defined by the inner product
$\Delta_{\boldsymbol{x}}\\{\bullet\\}\\!:=\nabla_{\boldsymbol{x}}\cdot\left(\nabla_{\boldsymbol{x}}\\{\bullet\\}\right)$.
Equation (5) clearly states that the overall bound pili number is not
conserved since the pili network constantly changes due to the binding and
unbinding pairs.

####  2.2.2 Stress definitions

The next step is to state the definitions for the active and passive stresses.
For the passive stress, we have a straightforward explicit expression which
reads

$\boldsymbol{\sigma}^{\text{p}}=-\displaystyle\frac{E\pi R^{2}c_{t}}{1-\pi
R^{2}c_{t}}\boldsymbol{i}\,,$ (6)

with $\boldsymbol{i}$ being the spatial second-order identity tensor, $E$
being the bulk modulus of the cell and $R$ being the cell radius. The
expression for the passive Cauchy stress is similar to a _hydrostatic
pressure-like stress_. This expression is one of the simplest forms to
describe the excluded volume interactions. It resembles the van der Waals gas
law for pressure in the absence of attractive interactions. As a result, in
the absence of the active stress, the only forces acting on cells would be
steric repulsive forces and cell-substrate friction, thus while moving on the
substrate, the cells tend to repel each other until they reach a state of
equilibrium. As it will be elucidated in the numerical examples, this gives
rise to formation of a uniform distribution of cells throughout the domain.
Figure 4 depicts the behavior of the passive stress with respect to the cell
number density. The $x$-axis represents the dimensionless parameter
$c_{t}R^{2}$ with the left $y$-axis representing $1-\pi R^{2}c_{t}$ and the
right $y$-axis representing the pressure term in the passive stress. At
$c_{t}R^{2}=1/\pi$, the stress renders an asymptotic behavior since the
denominator tends to zero. As a result, this value represents the upper bound
for the cell number density.

Figure 4:  Illustration of the behavior of passive stress with respect to cell
number density. The denominator of the stress becomes zero at the value
$\text{max}(c_{t}R^{2})=1/\pi$ hence, the asymptotic behavior. This value
represents the upper bound for the cell number density.

The definition of the active Cauchy stress is given as coarse-grained average
of the microscopic expression (for details see [1]). For us, its time
evolution equation is of interest which reads

$\displaystyle\frac{\partial\boldsymbol{\sigma}^{\text{a}}}{\partial
t}+\nabla_{\\!\boldsymbol{x}}\\!\cdot(\boldsymbol{\sigma}^{\text{a}}\otimes\boldsymbol{v})-2\left[\boldsymbol{l}\cdot\boldsymbol{\sigma}^{\text{a}}\right]^{\text{sym}}=-\frac{1}{\ell_{0}p_{t}f^{\text{p}}}\left[\boldsymbol{\sigma}^{\text{a}}\otimes\boldsymbol{\sigma}^{\text{a}}\right]:\boldsymbol{l}^{\text{sym}}+\boldsymbol{\sigma}^{\text{f}}-k_{\text{off}}\,\boldsymbol{\sigma}^{\text{a}}\,,$
(7)

with $\boldsymbol{l}$ being the spatial velocity gradient
$\boldsymbol{l}=\mbox{grad}\boldsymbol{v}$ and
$\boldsymbol{\sigma}^{\text{f}}$ being the pili formation Cauchy stress tensor
which will be defined below. The active Cauchy stress results from the pili-
pili mediated attractive forces between the cells which pull the cells
together and lead to the formation of an aggregate. For the pili formation
Cauchy stress an explicit formulation is achievable. Due to the intrinsic
nature of pili-pili interactions as force dipoles, this stress tensor has a
nematic symmetry and reads

$\boldsymbol{\sigma}^{\text{f}}=\displaystyle\frac{f^{\text{p}}k_{\text{on}}\ell_{0}}{2}\left[c_{t}^{2}\boldsymbol{i}-\displaystyle\frac{3\ell_{0}^{2}}{4}\left[2\left[\nabla_{\boldsymbol{x}}c_{t}\otimes\nabla_{\boldsymbol{x}}c_{t}\right]+\big{|}\nabla_{\boldsymbol{x}}c_{t}\big{|}^{2}\boldsymbol{i}\right]+\displaystyle\frac{3\ell_{0}^{2}}{4}\left[c_{t}\left[\Delta_{\boldsymbol{x}}c_{t}\right]\boldsymbol{i}+2c_{t}\nabla^{2}_{\boldsymbol{x}}c_{t}\right]\right]\,,$
(8)

where $\nabla^{2}_{\boldsymbol{x}}$ denotes the spatial Hessian operator
defined by
$\nabla^{2}_{\boldsymbol{x}}\\{\bullet\\}:=\nabla_{\boldsymbol{x}}\left(\nabla_{\boldsymbol{x}}\\{\bullet\\}\right)$.

###  2.3 Lagrangian approach

This section elaborates the pertinent equations within a Lagrangian framework,
hence, all the equations are written in the material configuration. In the
Lagrangian approach, we follow the trajectory of the cells as time elapses.
Key geometric notions from nonlinear continuum kinematics are thereby the
integration of the cell velocity into the cell nonlinear deformation map
$\boldsymbol{x}=\boldsymbol{y}(\boldsymbol{X},t)$ with $\boldsymbol{x}$ and
$\boldsymbol{X}$ denoting the Eulerian and the Lagrangian coordinates. Based
on the notion of the cell deformation map $\boldsymbol{y}$, we introduce the
corresponding material deformation gradient
$\boldsymbol{F}:=\nabla_{\boldsymbol{X}}\boldsymbol{y}$, its Jacobian
(determinant) $J:=\det\boldsymbol{F}$ and its cofactor
$\boldsymbol{K}:=J\boldsymbol{F}^{-T}$. As will be elucidated, the main
advantage of the Lagrangian approach is that it yields a significant
simplification for the equations which facilitates finite element
implementation of the problem. In addition, employing the Lagrangian approach
enables utilization of implicit time integration schemes which is
computationally more robust compared to explicit time integration and less
prone to instability issues. A significant step towards the Lagrangian
formulation of the problem is parametrization of all fields in Lagrangian
coordinates $\boldsymbol{X}$. Especially the spatial cell number density is
parameterized as $c_{t}=c_{t}(\boldsymbol{X},t)$. This critical step allows us
to adopt $\mathcal{C}^{0}$ continuous or discontinuous schemes and to avoid
complications regarding the implementation of $\mathcal{C}^{1}$ continuous
elements.

####  2.3.1 Balance equations

Similar to the previous section, we start with the overall cell number
conservation equation. Using the material time derivative
$D\\{\bullet\\}/Dt=\dot{\\{\bullet\\}}=\partial\\{\bullet\\}/\partial
t+\nabla_{\boldsymbol{x}}\\{\bullet\\}\cdot\boldsymbol{v}$ and Eq. (3), the
Lagrangian form of the overall cell number conservation equation can be
written as

$\displaystyle\frac{D\left(Jc_{t}\right)}{Dt}=\dot{\overline{Jc_{t}}}=\dot{J}c_{t}+\dot{c_{t}}J=0\,.$
(9)

The second balance equation is the linear momentum balance. To define the
Lagrangian version of the linear momentum balance equation, the velocity field
is replaced by the time derivative of the deformation map as
$\boldsymbol{v}=\dot{\boldsymbol{y}}$. Using the relation
$\nabla_{\boldsymbol{X}}\cdot\boldsymbol{P}=J\nabla_{\boldsymbol{x}}\cdot\boldsymbol{\sigma}$
the material linear momentum balance reads

$\nabla_{\boldsymbol{X}}\cdot\boldsymbol{P}^{\text{a}}+\nabla_{\boldsymbol{X}}\cdot\boldsymbol{P}^{\text{p}}-\xi
Jc_{t}\dot{\boldsymbol{y}}=\boldsymbol{0}\,,$ (10)

with $\boldsymbol{P}^{\text{a}}$ and $\boldsymbol{P}^{\text{p}}$ being the
active and passive Piola stresses, respectively. The next step is to define
the Lagrangian version of the bound pili evolution equation. Using the
relations $p_{0}=Jp_{t}$ and
$\dot{\\{\bullet\\}}=\partial\\{\bullet\\}/\partial
t+\nabla_{\boldsymbol{x}}\\{\bullet\\}\cdot\boldsymbol{v}$, the bound pili
evolution equation can be written in terms of the material pili density as

$\dot{p_{0}}-Jk_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\Delta_{\boldsymbol{x}}c_{t}-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\nabla_{\boldsymbol{x}}c_{t}\big{|}^{2}\right]+k_{\text{off}}p_{0}=0\,.$
(11)

Note, in Eq. (11), the derivatives of the cell number density are still with
respect to the spatial coordinates. In Section 2.3.4, these derivatives will
be transformed completely to the material configuration hence, our fully
Lagrangian formulation.

####  2.3.2 Stress definitions

Similar to the passive Cauchy stress, derivation of the passive Piola stress
is straightforward. Using the relation
$\boldsymbol{P}=J\boldsymbol{\sigma}\cdot\boldsymbol{F}^{-T}$, the passive
Piola stress reads

$\boldsymbol{P}^{\text{p}}=-\displaystyle\frac{E\pi R^{2}c_{t}}{1-\pi
R^{2}c_{t}}\boldsymbol{K}\,.$ (12)

Derivation of the active Piola stress from the relation (7) is an intricate
task and requires further attention. The active Piola stress tensor, a two-
point tensor, must be converted to a fully material stress tensor. In doing
so, firstly the active Cauchy stress $\boldsymbol{\sigma}^{\text{a}}$ is
transformed to the active Kirchhoff stress $\boldsymbol{\tau}^{\text{a}}$ as

$\boldsymbol{\tau}^{\text{a}}=J\boldsymbol{\sigma}^{\text{a}}\,,$ (13)

with its material time derivative

$\dot{\boldsymbol{\tau}^{\text{a}}}=J\left[\displaystyle\frac{\partial\boldsymbol{\sigma}^{\text{a}}}{\partial
t}+\nabla_{\boldsymbol{x}}\\!\cdot(\boldsymbol{\sigma}^{\text{a}}\otimes\boldsymbol{v})\right]\,.$
(14)

Via multiplying Eq. (7) by the Jacobian $J$, we arrive at a corresponding
relation for Eq. (7) in terms of the Kirchhoff stress

$\dot{\boldsymbol{\tau}^{\text{a}}}-2\left[\boldsymbol{l}\cdot\boldsymbol{\tau}^{\text{a}}\right]^{\text{sym}}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{\tau}^{\text{a}}\otimes\boldsymbol{\tau}^{\text{a}}\right]:\boldsymbol{l}^{\text{sym}}+\boldsymbol{\tau}^{\text{f}}-k_{\text{off}}\,\boldsymbol{\tau}^{\text{a}}\,.$
(15)

To obtain a fully material form of this equation, the Lie time derivative of
the Kirchhoff stress needs to be introduced. To express the Lie time
derivative of the Kirchhoff stress, we need to calculate the pull-back of the
Kirchhoff stress, calculate its material time derivative and then push it
forward to the spatial configuration again which renders

$\displaystyle\mathcal{L}_{t}\boldsymbol{\tau}^{\text{a}}=\dot{\boldsymbol{\tau}^{\text{a}}}-2\left[\boldsymbol{l}\cdot\boldsymbol{\tau}^{\text{a}}\right]^{\text{sym}}\,.$
(16)

Note, the symmetry of the Kirchhoff stress is utilized in the above
derivation. Further details regarding the formulation of the Lie time
derivative are available in B. The pull-back of the active Kirchhoff stress
yields the active Piola–Kirchhoff stress as

$\boldsymbol{S}^{\text{a}}=\boldsymbol{F}^{-1}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{F}^{-T}\,.$
(17)

Using Eqs. (15) and (16), one can write

$\mathcal{L}_{t}\boldsymbol{\tau}^{\text{a}}=\boldsymbol{F}\cdot\dot{\boldsymbol{S}^{\text{a}}}\cdot\boldsymbol{F}^{T}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{\tau}^{\text{a}}\otimes\boldsymbol{\tau}^{\text{a}}\right]:\boldsymbol{l}^{\text{sym}}+\boldsymbol{\tau}^{\text{f}}-k_{\text{off}}\,\boldsymbol{\tau}^{\text{a}}\,.$
(18)

Accordingly, the material time derivative of the active Piola–Kirchhoff stress
in terms of the active Kirchhoff stress is obtained as

$\dot{\boldsymbol{S}^{\text{a}}}=\boldsymbol{F}^{-1}\cdot\left[\mathcal{L}_{t}\boldsymbol{\tau}^{\text{a}}\right]\cdot\boldsymbol{F}^{-T}=\boldsymbol{F}^{-1}\cdot\left[-\displaystyle\frac{1}{\ell_{0}\rho_{t}^{\text{p}}f^{\text{p}}}\left[\boldsymbol{\tau}^{\text{a}}:\boldsymbol{l}^{{}^{\text{sym}}}\right]\boldsymbol{\tau}^{\text{a}}+\boldsymbol{\tau}^{\text{f}}-k_{\text{off}}\,\boldsymbol{\tau}^{\text{a}}\right]\cdot\boldsymbol{F}^{-T}\,.$
(19)

Finally, via replacing $\boldsymbol{\tau}^{\text{a}}$ with
$\boldsymbol{F}\cdot\boldsymbol{S}^{\text{a}}\cdot\boldsymbol{F}^{T}$, we
arrive at the final form of the material time derivative of the active
Piola–Kirchhoff stress

$\dot{\boldsymbol{S}^{\text{a}}}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\dot{\boldsymbol{E}}\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}\,,$
(20)

with $\boldsymbol{E}$ being the Green–Lagrange strain tensor
$\boldsymbol{E}=\tfrac{1}{2}[\boldsymbol{F}^{T}\cdot\boldsymbol{F}-\boldsymbol{I}]$.
For further details regarding the derivation of Eq. (20), see B. The most
significant advantage of using the Lagrangian approach is that it yields
considerable simplification to the governing equations and the active stress
time evolution equation. The active stress can be numerically determined via
time-discretizing Eq. (20), see C for further details. Subsequently, after
solving for the active Piola–Kirchhoff stress, the active Piola stress can be
immediately obtained via
$\boldsymbol{P}^{\text{a}}=\boldsymbol{F}\cdot\boldsymbol{S}^{\text{a}}$.
Stating the material time derivative of the active Piola–Kirchhoff stress
(20), severely alleviates incrementally objective time integration of the
active stress evolution. The next step is to define the pili formation second
Piola–Kirchhoff stress $\boldsymbol{S}^{\text{f}}$. Inserting the relation
$\boldsymbol{S}=J\boldsymbol{F}^{-1}\cdot\boldsymbol{\sigma}\cdot\boldsymbol{F}^{-T}$
into Eq. (8), yields the pili formation Piola–Kirchhoff stress as

$\displaystyle\boldsymbol{S}^{\text{f}}=\frac{1}{2}f^{\text{p}}k_{\text{on}}\ell_{0}J\left[c_{t}^{2}\boldsymbol{B}+\frac{3\ell_{0}^{2}}{4}\left[-2\boldsymbol{F}^{-1}\\!\\!\cdot\left[\nabla_{\boldsymbol{x}}c_{t}\otimes\nabla_{\boldsymbol{x}}c_{t}\right]\cdot\boldsymbol{F}^{-T}\\!\\!\\!-\big{|}\nabla_{\boldsymbol{x}}c_{t}\big{|}^{2}\boldsymbol{B}+c_{t}\left[\Delta_{\boldsymbol{x}}c_{t}\right]\boldsymbol{B}+c_{t}\left[\boldsymbol{F}^{-1}\\!\\!\cdot\nabla^{2}_{\boldsymbol{x}}c_{t}\cdot\boldsymbol{B}+\boldsymbol{B}\cdot\nabla^{2}_{\boldsymbol{x}}c_{t}\cdot\boldsymbol{F}^{-T}\right]\right]\right]\,,$
(21)

with $\boldsymbol{B}=\boldsymbol{F}^{-T}\cdot\boldsymbol{F}^{-1}$ the Piola
deformation tensor. Note, in Eq. (21), the gradients and divergences are with
respect to the spatial coordinates which will be transformed completely to the
material configuration in Section 2.3.4.

####  2.3.3 Weak form

To obtain the weak form of the overall cell number conservation equation, its
strong form is multiplied by a scalar-valued test function $\delta c_{t}$ and
then integrated over the referential domain as follows

$\int_{\mathcal{B}_{0}}\left[\dot{J}c_{t}+\dot{c_{t}}J\right]\delta
c_{t}\,\mbox{d}V=0\,\qquad\forall\delta c_{t}\,.$ (22)

To obtain the weak form of the linear momentum balance, its strong form too is
multiplied by a vector-valued test function $\delta\boldsymbol{y}$ and then
integrated over the referential domain which yields

$\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{a}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V=\int_{\partial\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{t}^{\text{a}}\,\mbox{d}A+\int_{\partial\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{t}^{\text{p}}\,\mbox{d}A\,\qquad\forall\delta\boldsymbol{y},$
(23)

where $\boldsymbol{t}^{\text{a}}$ and $\boldsymbol{t}^{\text{p}}$ are active
and passive tractions acting on the boundary of the domain, respectively.
Similarly, the weak form of the bound pili evolution equation is obtained via
multiplying Eq. (11) by a scalar-valued test function $\delta p_{0}$ and then
integrating over the referential domain as

$\int_{\mathcal{B}_{0}}\left[\dot{p_{0}}-Jk_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\Delta_{\boldsymbol{x}}c_{t}-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\nabla_{\boldsymbol{x}}c_{t}\big{|}^{2}\right]+k_{\text{off}}p_{0}\right]\delta
p_{0}\,\mbox{d}V=0\,\qquad\forall\delta p_{0}.$ (24)

####  2.3.4 Gradient enhanced framework

So far, the unknowns of our problem in Lagrangian approach have been the cell
number density $c_{t}$, the cell deformation map $\boldsymbol{y}$ and bound
pili number density $p_{0}$. A further challenge that arises in our problem is
the presence of the Hessian of the spatial cell number density
$\nabla^{2}_{\boldsymbol{x}}c_{t}$ and the Laplacian of the spatial cell
number density $\Delta_{\boldsymbol{x}}c_{t}$ in the bound pili number density
evolution equation (11) and the pili density formation Piola stress equation
(21). Thus, second derivatives of the spatial cell number density $c_{t}$ are
required to be calculated. In doing so, different strategies have been
employed in the literature, among which the well-established ones are
$\mathcal{C}^{1}$ continuous elements [75, 76], isogeometric analysis [77,
78], micromorphic continuum approach and the gradient enhanced framework [79,
80, 81]. In order to stay within the realm of the classical finite element
method (FEM) associated with a $\mathcal{C}^{0}$ continuous interpolation
approach, we adopt the gradient enhanced framework. In doing so, we introduce
an additional independent spatial vector field $\boldsymbol{g}$ to represent
the spatial gradient of the cell number density $\nabla_{\boldsymbol{x}}c_{t}$
and _weakly enforce_
$\boldsymbol{g}\overset{!}{=}\nabla_{\boldsymbol{x}}c_{t}$. Note, similar to
the spatial cell number density, its gradient is also parametrized as a
function of the Lagrangian coordinates
$\boldsymbol{g}=\boldsymbol{g}(\boldsymbol{X},t)$. The first step is to
transform the Laplacian and Hessian to substitute the terms including spatial
derivatives of $c_{t}$ with $\boldsymbol{g}$ as

$\displaystyle\Delta_{\boldsymbol{x}}c_{t}=\nabla_{\boldsymbol{x}}\cdot(\nabla_{\boldsymbol{x}}c_{t})=\nabla_{\boldsymbol{x}}\cdot\boldsymbol{g}\,,$
(25)
$\displaystyle\nabla^{2}_{\boldsymbol{x}}c_{t}=\nabla_{\boldsymbol{x}}(\nabla_{\boldsymbol{x}}c_{t})=\nabla_{\boldsymbol{x}}^{\text{sym}}\boldsymbol{g}=\frac{1}{2}\left[\nabla_{\boldsymbol{x}}\boldsymbol{g}+\nabla_{\boldsymbol{x}}^{T}\boldsymbol{g}\right]\,,$

where a symmetric gradient of $\boldsymbol{g}$ follows from the symmetry of
the Hessian of the cell number density. Afterwards the resulting terms must be
pulled back to the material configuration so as to unify the Lagrangian
formalism. For an arbitrary vector field $\boldsymbol{a}$, the following
relations hold between the spatial and material divergence and gradient
$\nabla_{\boldsymbol{x}}\cdot\boldsymbol{a}=\nabla_{\boldsymbol{X}}\boldsymbol{a}:\boldsymbol{F}^{-T}$
and
$\nabla_{\boldsymbol{x}}\boldsymbol{a}=\nabla_{\boldsymbol{X}}\boldsymbol{a}\cdot\boldsymbol{F}^{-1}$.
Accordingly, the Laplacian and Hessian of the cell number density can be
written in terms of $\boldsymbol{g}$ in the material configuration as

$\displaystyle\Delta_{\boldsymbol{x}}c_{t}=\nabla_{\boldsymbol{x}}\cdot\boldsymbol{g}=\nabla_{\boldsymbol{X}}\boldsymbol{g}:\boldsymbol{F}^{-T}\,,$
(26)
$\displaystyle\nabla^{2}_{\boldsymbol{x}}c_{t}=\nabla_{\boldsymbol{x}}^{\text{sym}}\boldsymbol{g}=\frac{1}{2}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}\cdot\boldsymbol{F}^{-1}+\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}^{T}\boldsymbol{g}\right]\,.$

Using Eq. (26), the strong (11) and weak (24) forms of the bound pili number
density evolution equation can be stated as fully Lagrangian equations as

$\dot{p_{0}}-Jk_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}:\boldsymbol{F}^{-T}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\boldsymbol{g}\big{|}^{2}\right]+k_{\text{off}}p_{0}=0\,,$
(27)

and

$\int_{\mathcal{B}_{0}}\left[\dot{p_{0}}-Jk_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}:\boldsymbol{F}^{-T}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\boldsymbol{g}\big{|}^{2}\right]+k_{\text{off}}p_{0}\right]\delta
p_{0}\,\mbox{d}V=0\,,$ (28)

respectively. Accordingly the pili density formation Piola–Kirchhoff stress
equation (21) can be rewritten as

$\displaystyle\boldsymbol{S}^{\text{f}}=\frac{1}{2}f^{\text{p}}k_{\text{on}}\ell_{0}J\left[c_{t}^{2}\boldsymbol{B}+\frac{3\ell_{0}^{2}}{4}\left[-2\boldsymbol{F}^{-1}\\!\\!\cdot\left[\boldsymbol{g}\otimes\boldsymbol{g}\right]\cdot\boldsymbol{F}^{-T}\\!\\!\\!-\left[\boldsymbol{g}\right]^{2}\boldsymbol{B}+c_{t}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}:\boldsymbol{F}^{-T}\right]\boldsymbol{B}+c_{t}\left[\boldsymbol{F}^{-1}\\!\\!\cdot\nabla_{\boldsymbol{X}}\boldsymbol{g}\cdot\boldsymbol{B}+\boldsymbol{B}\cdot\nabla_{\boldsymbol{X}}^{T}\boldsymbol{g}\cdot\boldsymbol{F}^{-T}\right]\right]\right]\,.$
(29)

####  2.3.5 Fictitious variational approach

The final step to complete our formulation is to update the governing
equations given the new continuity requirement. The advantage of using a
(fictitious) variational approach is that it provides the admissible forms of
tractions and external forces without any prior knowledge or assumptions. In
addition, with no extra effort or derivations, the variational approach
furnishes the weak form which is central to the finite element implementation
of the problem. In order to obtain the updated governing equations, a total
energy functional is minimized via setting its first variation to zero. To
begin, we assume that there exist a fictitious energy function $\Psi_{1}$
whose variation with respect to the degrees of freedom yields the already
existing weak forms (22), (23) and (28). The term “fictitious” for this energy
function implies that the original weak forms were directly derived from the
strong forms and not from any specific energy function and, as mentioned
before, the strong forms themselves are obtained using the coarse-grained
continuum approach developed in [1, 70]. To enforce the continuity of the
spatial cell number density gradient, we introduce a new energy function
$\Psi_{2}$ whose definition will be provided shortly. The total energy
$\Psi_{\text{tot}}$ consists of the initial fictitious energy and the newly
defined energy

$\Psi_{\text{tot}}=\Psi_{1}+\Psi_{2}\,.$ (30)

The energies are the integrals of their corresponding internal energy
densities over their associated domains as

$\Psi_{\text{tot}}=\int_{\mathcal{B}_{0}}\psi_{\text{tot}}\,\mbox{d}V\,,\qquad\Psi_{1}=\int_{\mathcal{B}_{0}}\psi_{1}\,\mbox{d}V\,,\qquad\Psi_{2}=\int_{\mathcal{B}_{0}}\psi_{2}\,\mbox{d}V\,.$
(31)

To minimize $\Psi_{\text{tot}}$, its first variation is set to zero. That is

$\delta\Psi_{\text{tot}}\stackrel{{\scriptstyle\boldsymbol{.}}}{{=}}0\qquad\Longrightarrow\qquad\delta\Psi_{1}+\delta\Psi_{2}\stackrel{{\scriptstyle\boldsymbol{.}}}{{=}}0\,.$
(32)

The field variables in our problem are the cell number density $c_{t}$, cell
deformation map $\boldsymbol{y}$, bound pili number density $p_{0}$ and cell
number density gradient $\boldsymbol{g}$. The energy density function to
impose the continuity of the spatial cell density gradient reads

$\psi_{2}=\psi_{2}\left(c_{t},\boldsymbol{y},\boldsymbol{g}\right)=\displaystyle\frac{1}{2}J\lambda\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{g}\right]^{2}\,,$
(33)

with $\lambda$ being the penalty parameter. The penalty parameter determines
how strongly the condition
$\boldsymbol{g}\overset{!}{=}\nabla_{\boldsymbol{x}}c_{t}$ is satisfied. Note,
the dependence of $\psi_{2}$ on $\boldsymbol{y}$ is through the deformation
gradient $\boldsymbol{F}$ and its determinant $J$. Accordingly, one could
write

$\int_{\mathcal{B}_{0}}\\!\\!\delta\psi_{\text{tot}}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\delta\psi_{1}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\delta\left(\displaystyle\frac{1}{2}J\lambda\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{g}\right]^{2}\right)\,\mbox{d}V\,.$
(34)

Calculating the variation of the overall energy density with respect to all
degrees of freedom and their gradients, yields the weak forms as

$\int_{\mathcal{B}_{0}}\delta\psi_{\text{tot}}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\\!\\!\\!\delta\psi_{1}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\delta\psi_{2}\,\mbox{d}V\,,$
(35)

with

$\displaystyle\int_{\mathcal{B}_{0}}\delta\psi_{1}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{c_{t}}\psi_{1}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{\boldsymbol{y}}\psi_{1}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{p_{0}}\psi_{1}\,\mbox{d}V\,,$
(36)
$\displaystyle\int_{\mathcal{B}_{0}}\delta\psi_{2}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{c_{t}}\psi_{2}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{\boldsymbol{y}}\psi_{2}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\delta_{\boldsymbol{g}}\psi_{2}\,\mbox{d}V\,.$

The variations of the newly defined energy function read

$\displaystyle\int_{\mathcal{B}_{0}}\\!\delta_{c_{t}}\psi_{2}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\lambda
J\left[\boldsymbol{B}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{F}^{-1}\cdot\boldsymbol{g}\right]\cdot\nabla_{\boldsymbol{X}}\delta
c_{t}\,\mbox{d}V\,,$ (37)
$\displaystyle\int_{\mathcal{B}_{0}}\\!\delta_{\boldsymbol{y}}\psi_{2}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\\!\underbrace{\left[\displaystyle\frac{1}{2}\lambda
J\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{g}\right]^{2}\boldsymbol{F}^{-T}-\lambda
J\left[\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}\right]\otimes\left[\boldsymbol{B}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{F}^{-1}\cdot\boldsymbol{g}\right]\right]\right]}_{\overline{\boldsymbol{P}}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V\,,$
$\displaystyle\int_{\mathcal{B}_{0}}\\!\delta_{\boldsymbol{g}}\psi_{2}\,\mbox{d}V=-\int_{\mathcal{B}_{0}}\lambda
J\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{g}\right]\cdot\delta\boldsymbol{g}\,\mbox{d}V\,,$

with $\overline{\boldsymbol{P}}$ being the cell-density-gradient-continuity-
induced Piola stress which arises due to the dependence of $\psi_{2}$ on the
deformation. Finally, we arrive at the final weak form of the governing
equations which form our residual system as

overall cell number conservation:
$\displaystyle\int_{\mathcal{B}_{0}}\\!\\!\\!\left[\dot{J}c_{t}+\dot{c_{t}}J\right]\delta
c_{t}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\\!\\!\\!\lambda
J\left[\boldsymbol{B}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{F}^{-1}\cdot\boldsymbol{g}\right]\cdot\nabla_{\boldsymbol{X}}\delta
c_{t}\,\mbox{d}V=0\,,$ (38) linear momentum balance:
$\displaystyle\int_{\mathcal{B}_{0}}\\!\\!\\!\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V+\\!\\!\\!\int_{\mathcal{B}_{0}}\\!\\!\\!\left[\boldsymbol{P}^{\text{a}}+\boldsymbol{P}^{\text{p}}+\overline{\boldsymbol{P}}\right]:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V-\int_{\partial\mathcal{B}_{0}}\\!\\!\\!\\!\delta\boldsymbol{y}\cdot\left[\boldsymbol{t}^{\text{a}}+\boldsymbol{t}^{\text{p}}\right]\,\mbox{d}A=\boldsymbol{0}\,,$
bound pili number density evolution:
$\displaystyle\int_{\mathcal{B}_{0}}\left[\dot{p_{0}}-Jk_{\text{on}}\left[c_{t}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}:\boldsymbol{F}^{-T}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\boldsymbol{g}\big{|}^{2}\right]+k_{\text{off}}p_{0}\right]\delta
p_{0}\,\mbox{d}V=0\,,$ cell number density gradient continuity:
$\displaystyle\int_{\mathcal{B}_{0}}\\!\\!\\!-\lambda
J\left[\boldsymbol{F}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t}-\boldsymbol{g}\right]\cdot\delta\boldsymbol{g}\,\mbox{d}V=\boldsymbol{0}\,.$

##  3 Finite element implementation

In this section, we present a general finite element formulation for the
implementation of our developed methodology. The first step towards finite
element implementation of our problem is the derivation of the discretized
weak forms which has been concluded in Section 2.3.5. For time integration,
the time interval $\mathds{T}$ is subdivided into a set of intervals $\Delta
t$ with

$\mathds{T}=\bigcup_{n=0}^{\\#\text{ts}-1}\left[t_{n},\,t_{n+1}\right]\,,$
(39)

where $\\#\text{ts}$ denotes the number of time steps and the time increment
is defined by $\Delta t=t_{n+1}-t_{n}$. As mentioned earlier, an implicit time
integration associated with the backward Euler method is employed here in
which an unknown $y_{n+1}$ is determined via the relation
$y_{n+1}-y_{n}=f(y_{n+1})\Delta t$ based on the continuous evolution equation
$\dot{y}=f(y)$.

Here, the spatial discretization of the problem domain is carried out using
the Bubnov–Galerkin finite element method. The geometry together with the
variables are approximated using isoparametric coordinates
$\boldsymbol{\xi}\in[-1,1]^{2}$. Using standard interpolations together with
the isoparametric concept, the geometry is discretized as

$\displaystyle\boldsymbol{X}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\boldsymbol{X}(\boldsymbol{\xi})=\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})\boldsymbol{X}^{i}\,,$
(40)

with $\\#\text{e}$ denoting the number of elements and $\boldsymbol{\xi}$
denoting the natural space coordinates. Accordingly the discretized fields
read

$\displaystyle c_{t}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx
c_{t}(\boldsymbol{\xi})=\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})c_{t}^{i}\,,$
$\displaystyle\delta c_{t}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\delta
c_{t}(\boldsymbol{\xi})=\\!\\!\\!\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})\delta
c_{t}^{i}\,,$ (41)
$\displaystyle\boldsymbol{y}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\boldsymbol{y}(\boldsymbol{\xi})=\sum_{i=1}^{\\#\text{e}}M^{i}(\boldsymbol{\xi})\boldsymbol{y}^{i}\,,$
$\displaystyle\delta\boldsymbol{y}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\delta\boldsymbol{y}(\boldsymbol{\xi})=\\!\\!\\!\sum_{i=1}^{\\#\text{e}}M^{i}(\boldsymbol{\xi})\delta\boldsymbol{y}^{i}\,,$
$\displaystyle p_{0}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx
p_{0}(\boldsymbol{\xi})=\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})p_{0}^{i}\,,$
$\displaystyle\delta p_{0}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\delta
p_{0}(\boldsymbol{\xi})=\\!\\!\\!\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})\delta
p_{0}^{i}\,,$
$\displaystyle\boldsymbol{g}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\boldsymbol{g}(\boldsymbol{\xi})=\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})\boldsymbol{g}^{i}\,,$
$\displaystyle\delta\boldsymbol{g}\Big{|}_{\mathcal{B}_{0}}\\!\\!\\!\approx\delta\boldsymbol{g}(\boldsymbol{\xi})=\\!\\!\\!\sum_{i=1}^{\\#\text{e}}N^{i}(\boldsymbol{\xi})\delta\boldsymbol{g}^{i}\,,$

with $N^{i}$ and $M^{i}$ being different shape functions, which, however, can
also be chosen identical in particular cases. For the cases where a mixed
finite element formulation is required, these shape functions are polynomials
of different orders.

The discretized weak form of the overall cell number conservation equation
$\eqref{eq:weak_Lagrangian}_{1}$ reads

$\int_{\mathcal{B}_{0}}\left[c_{t_{n}}\displaystyle\frac{J_{n+1}-J_{n}}{\Delta
t}+J_{n+1}\displaystyle\frac{c_{t_{n+1}}-c_{t_{n}}}{\Delta
t}\right]N^{i}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\lambda
J_{n+1}\left[\boldsymbol{B}_{n+1}\cdot\nabla_{\boldsymbol{X}}c_{t_{n+1}}-\boldsymbol{F}_{n+1}^{-1}\cdot\boldsymbol{g}_{n+1}\right]\cdot\nabla_{\boldsymbol{X}}N^{i}\,\mbox{d}V=0\,.$
(42)

The discretized weak form of the linear momentum balance equation
$\eqref{eq:weak_Lagrangian}_{2}$ reads

$\int_{\mathcal{B}_{0}}\xi
J_{n+1}c_{t_{n+1}}\displaystyle\frac{\boldsymbol{y}_{n+1}-\boldsymbol{y}_{n}}{\Delta
t}M^{i}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\left[\boldsymbol{P}^{\text{a}}_{n+1}+\boldsymbol{P}^{\text{p}}_{n+1}+\overline{\boldsymbol{P}}_{n+1}\right]\cdot\nabla_{\boldsymbol{X}}M^{i}\,\mbox{d}V=\int_{\partial\mathcal{B}_{0}}\\!\\!\\!\boldsymbol{t}^{\text{a}}M^{i}\,\mbox{d}A+\\!\\!\\!\int_{\partial\mathcal{B}_{0}}\\!\\!\\!\boldsymbol{t}^{\text{p}}M^{i}\,\mbox{d}A\,.$
(43)

The discretized weak form of the bound pili number density evolution equation
$\eqref{eq:weak_Lagrangian}_{3}$ reads

$\int_{\mathcal{B}_{0}}\left[\displaystyle\frac{p_{0_{n+1}}-p_{0_{n}}}{\Delta
t}-J_{n+1}k_{\text{on}}\left[c_{t_{n+1}}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t_{n+1}}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}_{n+1}:\boldsymbol{F}_{n+1}^{-T}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\boldsymbol{g}_{n+1}\big{|}^{2}\right]+k_{\text{off}}p_{0_{n+1}}\right]N^{i}\,\mbox{d}V=0\,,$
(44)

and finally, the discretized weak form of the cell density gradient continuity
equation $\eqref{eq:weak_Lagrangian}_{4}$ reads

$\int_{\mathcal{B}_{0}}-\lambda
J_{n+1}\left[\boldsymbol{F}_{n+1}^{-T}\cdot\nabla_{\boldsymbol{X}}c_{t_{n+1}}-\boldsymbol{g}_{n+1}\right]N^{i}\,\mbox{d}V=\boldsymbol{0}\,.$
(45)

The fully-discrete form of balance equations can be obtained by applying the
spatial approximation which form a residual system associated with the global
node $I$ as

$\displaystyle\text{R}^{I}_{c}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\left[c_{t_{n}}\frac{J_{n+1}-J_{n}}{\Delta
t}+J_{n+1}\frac{c_{t_{n+1}}-c_{t_{n}}}{\Delta
t}\right]N^{i}\,\mbox{d}V+\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\lambda
J_{n+1}\left[\boldsymbol{B}_{n+1}\cdot\nabla_{\boldsymbol{X}}c_{t_{n+1}}-\boldsymbol{F}_{n+1}^{-1}\cdot\boldsymbol{g}_{n+1}\right]\cdot\nabla_{\boldsymbol{X}}N^{i}\,\mbox{d}V=0\,,$
(46)
$\displaystyle\text{R}^{I}_{\boldsymbol{y}}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\xi
J_{n+1}c_{t_{n+1}}\frac{\boldsymbol{y}_{n+1}-\boldsymbol{y}_{n}}{\Delta
t}M^{i}\,\mbox{d}V+\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\left[\boldsymbol{P}^{\text{a}}_{n+1}+\boldsymbol{P}^{\text{p}}_{n+1}+\overline{\boldsymbol{P}}_{n+1}\right]\cdot\nabla_{\boldsymbol{X}}M^{i}\,\mbox{d}V-\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#se}\int_{\partial\mathcal{B}_{0}}\\!\\!\\!\left[\boldsymbol{t}^{\text{a}}+\boldsymbol{t}^{\text{p}}\right]M^{i}\,\mbox{d}A=\boldsymbol{0}\,,$
$\displaystyle\text{R}^{I}_{p}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\left[\frac{p_{0_{n+1}}-p_{0_{n}}}{\Delta
t}-J_{n+1}k_{\text{on}}\left[c_{t_{n+1}}^{2}+\displaystyle\frac{3\ell_{0}^{2}}{4}c_{t_{n+1}}\left[\nabla_{\boldsymbol{X}}\boldsymbol{g}_{n+1}:\boldsymbol{F}_{n+1}^{-T}\right]-\displaystyle\frac{3\ell_{0}^{2}}{4}\big{|}\boldsymbol{g}_{n+1}\big{|}^{2}\right]+k_{\text{off}}p_{0_{n+1}}\right]N^{i}\,\mbox{d}V=0\,,$
$\displaystyle\text{R}^{I}_{\boldsymbol{g}}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}-\lambda
J_{n+1}\left[\boldsymbol{F}_{n+1}^{-T}\nabla_{\boldsymbol{X}}c_{t_{n+1}}-\boldsymbol{g}_{n+1}\right]N^{i}\,\mbox{d}V=\boldsymbol{0}\,,$

with $\\#se$ denoting the number of surface elements and
$\mathop{\mbox{\Large{{{A}}}}}$ being the assembly operator. The global
residual vector consist of the above four residual vectors as

$\text{R}_{\text{tot}}=\left[\begin{matrix}\text{R}_{c_{t}}\\\
\text{R}_{\boldsymbol{y}}\\\ \text{R}_{p_{0}}\\\
\text{R}_{\boldsymbol{g}}\end{matrix}\right]\,,\qquad\text{with}\qquad\text{R}_{c_{t}}=\left[\begin{matrix}\text{R}^{1}_{c_{t}}\\\
\text{R}^{2}_{c_{t}}\\\ \vdots\\\ \text{R}^{\text{nn}}_{c_{t}}\\\
\end{matrix}\right]\,,\qquad\text{R}_{\boldsymbol{y}}=\left[\begin{matrix}\text{R}^{1}_{\boldsymbol{y}}\\\
\text{R}^{2}_{\boldsymbol{y}}\\\ \vdots\\\
\text{R}^{\text{nn}}_{\boldsymbol{y}}\\\
\end{matrix}\right]\,,\qquad\text{R}_{p_{0}}=\left[\begin{matrix}\text{R}^{1}_{p_{0}}\\\
\text{R}^{2}_{p_{0}}\\\ \vdots\\\ \text{R}^{\text{nn}}_{p_{0}}\\\
\end{matrix}\right]\,,\qquad\text{R}_{\boldsymbol{g}}=\left[\begin{matrix}\text{R}^{1}_{\boldsymbol{g}}\\\
\text{R}^{2}_{\boldsymbol{g}}\\\ \vdots\\\
\text{R}^{\text{nn}}_{\boldsymbol{g}}\\\ \end{matrix}\right]\,,$ (47)

where nn denotes the total number of nodes. Finally, the fully discrete
nonlinear system of governing equations becomes

$\text{R}_{\text{tot}}=\text{R}_{\text{tot}}(\boldsymbol{U})\stackrel{{\scriptstyle\boldsymbol{!}}}{{=}}\boldsymbol{0}\,,$
(48)

with $\boldsymbol{U}$ being the global vector of unknowns. To find the
solution of the system (48), the Newton–Raphson scheme is employed. The
consistent linearization of the resulting system of equations yields

$\text{R}_{\text{tot}}(\boldsymbol{U}_{n+1})=\text{R}_{\text{tot}}(\boldsymbol{U}_{n})+\text{K}_{\text{tot}}\cdot\Delta\boldsymbol{U}_{n}\stackrel{{\scriptstyle\boldsymbol{!}}}{{=}}\boldsymbol{0}\qquad\text{with}\qquad\text{K}_{\text{tot}}=\displaystyle\frac{\partial\text{R}_{\text{tot}}}{\partial\boldsymbol{U}}\big{\lvert}_{n}\qquad\text{and}\qquad\boldsymbol{U}_{n+1}=\boldsymbol{U}_{n}+\Delta\boldsymbol{U}_{n}\,,$
(49)

where the subscript $n$ indicates the step number.

##  4 Numerical results

This section aims to illustrate our proposed theory through a set of numerical
examples. In doing so, the evolution of the cell number density under various
conditions is investigated. Moreover, parametric studies are carried out in
order to highlight the influence of different parameters on the evolution of
the cell number density. Throughout all the examples, the domain is a
$80\times 80$ square subject to periodic boundary conditions. Implicit time
integration is adopted with no stability condition which leads to robust
solution of our problem and larger time steps could be used. All the numerical
results are obtained from our in-house finite element code. The solution
procedure is robust and for all examples, we obtain convergence with a
quadratic rate associated with the Newton–Raphson scheme.

###  4.1 Cell number density evolution

Figure 5:  Snapshots of the cell number density evolution in the presence of
the active stress. The first and second rows represent the cell number density
distribution in undeformed and deformed configurations, respectively. The
third and the fourth rows illustrate how an Eulerian observer would see the
cell number density evolution through a fixed window associated with the
Eulerian framework.  Figure 6:  Another illustration of the cell aggregate
evolution in Fig. 5 with variable color-bar range.

This section investigates the evolution of the cells under two different
scenarios. In the first scenario, the cellular aggregate behavior is
investigated in the presence of all the driving forces. That is, the pili-pili
mediated attractive forces, steric repulsion forces and the cell-substrate
friction. In the second scenario, the problem is simplified to the case where
the pili-pili mediated attractive forces are eliminated. In this case, we
investigate the cell number density evolution when only the steric repulsion
forces and the cell-substrate friction are the forces acting on the cells.

Figure 5 renders six different snapshots of the cellular aggregate evolution
at different times associated with the first scenario. The color represents
the magnitude of cell number density throughout the domain. The first and
second rows represent the cell number density distribution in the undeformed
and deformed configurations, respectively. The third and the fourth rows
illustrate how an Eulerian observer would see the cellular aggregate evolution
through a fixed window associated with the Eulerian framework. At the initial
time $t=0$, the problem starts with a uniform distribution of cells with
initial $c_{t}=0.079$ perturbed with $\pm 0.001$ relative random fluctuation.
The cell radius is set to $R=1\,\mu\text{m}$, the cell bulk modulus is set to
$E=1\,\text{N}/\text{m}^{2}$ and the friction coefficient is set to
$\xi=10\,\text{Ns}/\text{m}$. In addition, the pili pair binding rate is
$k_{\text{on}}=0.0178\,\text{s}^{-1}$, the pili pair unbinding rate is
$k_{\text{off}}=0.01\,\text{s}^{-1}$, the average pili length is
$\ell_{0}=2.0\,\mathrm{\mu}\text{m}$ and the pili-pili mediated attractive
force is $f^{\text{p}}=18.0\,\text{pN}$. Since the backward Euler time
integration is adopted, an adaptive time stepping is employed with initial
time step of $\Delta t=0.5\,\text{s}$. If convergence is obtained at very few
iterations, the time step is enlarged by a factor of $1.2$ and vice versa.
Through the aggregate evolution, first an initial homogeneous smooth cell
distribution is observed due to the action of the passive stress. This is
justifiable since the formation of bound pili pairs (and thus emergence of
active stress) requires time hence, the dominance of the passive stress at the
starting steps. As time elapses, more bound pili pairs form as depicted in
Fig. 3. Thus cells tend to attract each other which is reflected in the
dominance of the active stress which triggers the onset of the phase
separation. Therefore, pili-mediated attractive forces act as the major
driving force which lead to the formation of an aggregate. Figure 6 provides
another illustration of the cell aggregate evolution in Fig. 5 with variable
color range in each step. The purpose of the figure is to elucidate the
initial random distribution of cells and highlight the early smooth
distribution of cells due to the passive stress which is reflected in the
tighter color bar range.

Figure 7:  A parametric study on the cell density evolution. The $y$-axis
shows the difference between the maximum and minimum cell density throughout
the domain while the $x$-axis shows the elapsed time. The left figure studies
the effects of the pili average length $\ell_{0}$ where
$k_{\text{on}}=0.05\,\text{s}^{-1}$ and $f^{p}=12\,\text{pN}$. The middle
figure studies the effects of the pili pair binding rate $k_{\text{on}}$ where
$\ell_{0}=2.0\,\mu\text{m}$ and $f^{p}=12\,\text{pN}$. The right figure
studies the effects of the pili pair attractive force $f^{p}$ where
$k_{\text{on}}=0.05\,\text{s}^{-1}$ and $\ell_{0}=2.0\,\mu\text{m}$.

In Fig. 7, a parametric study is carried out to investigate the influence of
the average pili length $\ell_{0}$, pili-pili binding rate $k_{\text{on}}$ and
pili-pili mediated attractive force $f^{\text{p}}$ on the formation of an
aggregate. Each figure renders the difference between the maximum and the
minimum cell number density versus time. The left figure studies the effects
of the pili average length $\ell_{0}$ where
$k_{\text{on}}=0.05\,\text{s}^{-1}$ and $f^{p}=12\,\text{pN}$, the middle
figure studies the effects of the pili pair binding rate $k_{\text{on}}$ where
$\ell_{0}=2.0\,\mu\text{m}$ and $f^{p}=12\,\text{pN}$ and the right figure
studies the effects of the pili pair attractive force $f^{p}$ where
$k_{\text{on}}=0.05\,\text{s}^{-1}$ and $\ell_{0}=2.0\,\mu\text{m}$. For all
figures, an initial decrease in the cell number density difference is observed
which is associated with function of the passive stress. This behavior is
vividly observed in the second step in Fig. 6 as the color bar limits tend to
tighten. It is observed that, if the values of $f^{p}$, $\ell_{0}$ and
$k_{\text{on}}$ are small, the pili-pili mediated forces cannot overcome the
steric repulsive forces and the difference continues to decrease until a
uniform homogeneous distribution is obtained. For these cases the phase
separation does not occur and the dynamics of the cell network is mainly
driven by the steric repulsive forces and cell-substrate friction. The lower
the values of $f^{p}$, $\ell_{0}$ and $k_{\text{on}}$, the faster the
equilibrium state is reached. However, if the values of $f^{p}$, $\ell_{0}$
and $k_{\text{on}}$ are large enough, the pili-pili mediated forces tend to
dominate the other forces as time elapses which is reflected in the larger
cell number density difference. In these cases, the active stress plays a more
considerable role compared to the passive stress which yields the onset of the
phase separation. Increasing the values of $f^{p}$, $\ell_{0}$ and
$k_{\text{on}}$ results in higher rates of phase separation as indicated by
larger slopes of graphs in Fig. 7. This behavior is natural to expect as
larger values of $f^{p}$ means that the attractive force between the pili
pairs is stronger which causes the cells to form a colony at a higher rate.
Larger values of $\ell_{0}$ indicate longer pili implying that pili are more
capable to reach out and attach to pili of other cells which facilitates the
formation of a colony. And finally, larger values of $k_{\text{on}}$ signify
faster binding rate between the pili leading to a faster formation of bound
pairs which accelerates the colony formation. Note, the parameters in the
graphs in Fig 7, totally follow the condition obtained by Kuan et al. [1, 70]
using linear stability analysis which reads

$-\displaystyle\frac{1}{\xi}\displaystyle\frac{E\pi
R^{2}}{\left[1-c_{t}^{\text{initial}}\pi
R^{2}\right]^{2}}+\displaystyle\frac{c_{t}^{\text{initial}}\ell_{0}f^{\text{p}}k_{\text{on}}}{k_{\text{off}}\xi}>0\,.$

Satisfying the above condition indicates the onset of the phase separation.

Figure 8: Illustration of the Taylor–Hood element utilized for the simulation.
For this element type we used quadratic approximation for the deformation map
field and linear approximation for the cell number density in order to observe
the LBB condition.

Now we consider the second scenario where the cell number density evolution is
examined while its dynamics is driven only by the steric repulsive forces and
the pili-substrate friction. In the absence of the the pili-pili mediated
attractive forces, the complexity of the problem reduces considerably since
the bound pili number density evolution equation is eliminated from the system
of equations. In addition, there is no need to adopt a gradient enhanced
framework to deal with the Laplacian and Hessian of the cell number density.
As a result, the field variables simply become the original cell number
density $c_{t}$ and the deformation map $\boldsymbol{y}$ and, under periodic
boundary condition, the fully discrete residual system (46) simplifies to

$\displaystyle\text{R}^{I}_{c}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\left[c_{t_{n}}\frac{J_{n+1}-J_{n}}{\Delta
t}+J_{n+1}\frac{c_{t_{n+1}}-c_{t_{n}}}{\Delta t}\right]N^{i}\,\mbox{d}V=0\,,$
(50)
$\displaystyle\text{R}^{I}_{\boldsymbol{y}}=\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\xi
J_{n+1}c_{t_{n+1}}\frac{\boldsymbol{y}_{n+1}-\boldsymbol{y}_{n}}{\Delta
t}M^{i}\,\mbox{d}V+\mathop{\mbox{\Large{{{A}}}}}_{\alpha=1}^{\\#e}\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}_{n+1}\cdot\nabla_{\boldsymbol{X}}M^{i}\,\mbox{d}V=\boldsymbol{0}\,.$

As mentioned earlier, pili-pili attractive forces are absent and the repulsive
forces between the cells are reflected in the passive stress. To cope with the
mixed nature of the problem, i.e. to respect the LBB inf-sup condition, we
allow interpolation of the cell number density $c_{t}$ and the deformation map
$\boldsymbol{y}$ with different polynomial order. We exploit the Taylor–Hood
element for our finite element implementation which is linear in $c_{t}$ and
quadratic in $\boldsymbol{y}$. Subsequently, the corresponding shape functions
$N^{i}$ and $M^{i}$ are polynomials of different orders. Figure 8 provides an
illustrative schematic of the Taylor–Hood elements adopted in our simulations.

Figure 9:  Different snapshots of the cell network evolution in the absence of
the active stress. The first and second rows represent the cell number density
distribution in undeformed and deformed configurations, respectively. The
third and the fourth rows illustrate how an Eulerian observer would see the
cell number density evolution through a fixed window associated with the
Eulerian framework. The dynamics of cell is only driven by the steric
repulsion forces and the cell-substrate friction. Uniform color throughout the
domain indicates uniform homogeneous distribution of the cells.

Figure 9 renders six different snapshots of the cell density evolution at
different times. The boundary and the initial conditions are similar to the
previous case study. Similarly the deformed and undeformed bodies are provided
in the first and second rows where the third and fourth rows are associated
with the Eulerian window. Since there exists no pili-pili attractive force, as
time elapses, the cells tend to repel each other until they reach to the state
of equilibrium. This equilibrium state is obtained when the cells are
uniformly distributed throughout the domain. It is observed that throughout
time evolution, the difference between the cell densities tend to vanish which
signifies the uniformity of the cell distribution. Moreover, due to the
absence of the active stress, the domain does not undergo any considerable
deformation where the deformation gradient remains close to identity.

In Fig. 10, a parametric study is carried out to investigate the influence of
the friction coefficient $\xi$ and cell bulk modulus $E$ on the cell number
density evolution. Each figure renders the difference between the maximum and
the minimum cell number density versus time. The left figure studies the
effects of the friction coefficient with $E=1\,\text{N}/\text{m}^{2}$ and
$R=1\,\mu\text{m}$ and the right figure studies the effects of the cell bulk
modulus with $\xi=10\,\text{Ns}/\text{m}$ and $R=1\,\mu\text{m}$. Smaller
density difference in each figure indicates more uniform distribution. It is
observed that increasing the friction coefficient delays the uniform
distribution which is understandable since the friction coefficient impedes
the cells movement. Thus, for larger friction coefficients, the cells require
more time to disperse due to the repulsive forces. On the other hand, we
observe that a uniform distribution is reached at a faster rate for larger
cell bulk modulus. This is also justifiable since larger bulk modulus
indicates more rigidity of the cells which gives rise to their quicker
dispersion after colliding with other cells.

Figure 10:  A parametric study on the cell density evolution in the absence of
the active stress. The $y$-axis shows the difference between the maximum and
minimum cell density throughout the domain as a measure of homogeneity of the
cell distribution while the $x$-axis shows the elapsed time. The left figure
renders the effects of friction coefficient on the cell density evolution when
$E=1\,\text{N}/\text{m}^{2}$ and $R=1\,\mu\text{m}$. The right figure renders
the effects of cell bulk modulus on the cell density evolution when
$R=1\,\mu\text{m}$ and $\xi=10\,\text{Ns}/\text{m}$.

###  4.2 Overall cell number conservation

It is noteworthy to mention that our mixed four-field
$\left(c_{t},\boldsymbol{y},p_{0},\boldsymbol{g}\right)$ implementation
conserves the total cell number in a Lagrangian solution domain exactly. To
ensure the discrete conservation of the total cell number throughout the
domain, we rewrite the strong form of the continuity equation (9), integrate
it over the referential domain and expand it as

$\int_{\mathcal{B}_{0}}\dot{\overline{Jc_{t}}}\,\mbox{d}V=\int_{\mathcal{B}_{0}}\left[\dot{J}c_{t}+\dot{c_{t}}J\right]\,\mbox{d}V\approx\int_{\mathcal{B}_{0}}\left[\displaystyle\frac{J_{n+1}-J_{n}}{\Delta
t}c_{t_{n}}+\displaystyle\frac{c_{t_{n+1}}-c_{t_{n}}}{\Delta
t}J_{n+1}\right]\,\mbox{d}V=0\,.$ (51)

Note, to satisfy the overall cell number conservation equation in the last
integral, the cell density in the first term must be chosen from step $n$
whereas the Jacobian in the second term must be chosen from step $n+1$.
Accordingly, Eq. (51) reads

$\int_{\mathcal{B}_{0}}J_{n+1}c_{t_{n+1}}\mbox{d}V-\int_{\mathcal{B}_{0}}J_{n}c_{t_{n}}\,\mbox{d}V=0\,,$
(52)

which can in turn be written as the integrals over the spatial domain as

$\int_{\mathcal{B}_{t_{n+1}}}c_{t_{n+1}}\mbox{d}v_{n+1}=\int_{\mathcal{B}_{t_{n}}}c_{t_{n}}\mbox{d}v_{n}\,,$
(53)

which is equivalent to the conservation of the total cell number. Figure 11
illustrates the percentage of the overall cell number change versus time for
the two case studies associated with Figs. 5 and 9. It is of great
significance to point out that no cell number loss occurs within our systems,
thus the variation of the total cell number remains zero without any
fluctuation, see [82, 83] for further discussions on the issue regarding mass
loss/production in Lagrangian formulations.

Figure 11:  Illustration of the discrete conservation of the overall cell
number for both case studies associated with Figs. 5 and 9. The $x$ axis
represents the elapsed time whereas the $y$ axis represents the percentage of
change in the overall cell number. The left figure corresponds to the case
where all forces are present in the dynamics of the system whereas the left
figure corresponds to the case where pili-pili mediated forces are absent.

###  4.3 Colony coalescence

Figure 12:  Illustration of the colony coalescence process. Five snapshots of
the coalescence process in both undeformed and deformed configurations are
depicted. The parameter $h$ measures the bridge length between the two merging
colonies. The center plots render the effects of $k_{\text{on}}$ and
$f^{\text{p}}$ on the evolution of the bridge between the colonies.

In this section, two colonies are put next to each other and the process of
their coalescence is examined. In doing so, we define a new parameter $h$
which measures the bridge length between the two merging colonies. From a
computational point of view, we define the boundary of the bridge where the
the cell number density gradient in $y$ direction is $96\%$ of its maximum
value throughout the domain. Figure 12 investigates the influence of the pili-
pili binding rate $k_{\text{on}}$ and the pili-pili mediated attractive force
$f^{\text{p}}$ on the size of the bridge. Five snapshots of the coalescence
process in both undeformed and deformed configurations are shown on top and
bottom, respectively. The snapshots correspond to the green line in the left
plot with $k_{\text{on}}=0.021\,\text{s}^{-1}$ and
$f^{\text{p}}=18\,\text{pN}$. The left plot renders the bridge length versus
time for three different values of the pili-pili binding rate with
$f^{\text{p}}=18\,\text{pN}$. The blue line corresponds to
$k_{\text{on}}=0.019\,\text{s}^{-1}$, the green line corresponds to
$k_{\text{on}}=0.021\,\text{s}^{-1}$ and the red line corresponds to
$k_{\text{on}}=0.023\,\text{s}^{-1}$. It is observed that the bridge grows
with a faster rate as we increase the pili-pili binding rate. The right plot
renders the bridge length versus time for three different values of the pili-
pili mediated attractive forces with $k_{\text{on}}=0.023\,\text{s}^{-1}$. The
blue line corresponds to $f^{\text{p}}=18\,\text{pN}$, the green line
corresponds to $f^{\text{p}}=31\,\text{pN}$ and the red line corresponds to
$f^{\text{p}}=43\,\text{pN}$. Similarly, larger values of $f^{\text{p}}$
result in a faster growth of the bridge. These two observations are
justifiable since larger binding rates and larger pili-pili attractive forces
imply quicker bound pili formation and stronger attractive forces between the
cells, respectively, and thus they yield quicker aggregate formation. In both
cases, the growth starts with smaller rate due to the function of the passive
stress. Further evolution of time leads the pili-pili mediated forces to
become more dominant which yields higher rate growths.

###  4.4 Aggregate position and periodicity

Figure 13:  Formation of aggregates under different initial random
distributions. In each block, the top figures render the undeformed
configuration whereas the bottom figures render the deformed bodies. Although
the initial distributions are the different, aggregates with similar sizes and
same degree of phase separations are obtained.  Figure 14:  Illustration of
the periodicity of the domain. Although the aggregate are formed in different
positions in the domain, due to periodicity, a single centered aggregate can
always be extracted if enough samples are put next to each other.

This section shows the formation of aggregates under different initial random
distribution of cells. Figure 13 renders three different cases with various
initial random distribution of cells. Each block contains the undeformed and
deformed configurations with five different snapshots from the domain. It is
observed that, despite different initial conditions, similar aggregates with
the same size and degree of the phase separation are obtained at the end. It
is noteworthy that due to the different initial distribution, the bodies have
deformed differently in order to yield the same aggregate. Although the
aggregates are formed in different locations in the three cases shown in Fig.
13, due to periodicity, a single centered aggregate can always be extracted if
enough samples are put next to each other. Figure 14 sheds light on this issue
more vividly. The final snapshots of the three cases in Fig. 13 are shown on
the top row in Fig. 14 . In the middle row, these samples are put together in
order to form a periodic structure. In the bottom row, it is shown that
similar centered aggregates can always be extracted.

##  5 Summary and outlook

A continuum framework to model and simulate the behavior of biological
cellular aggregates has been established. The process of micro-colony
formation has been described as an active phase separation phenomenon. It
turns out that employing the Lagrangian approach yields considerable
simplification of the equations in particular for the active stress time
evolution as compared to previously introduced Eulerian approach [1]. In
addition to satisfying the conservation of the total cell number, our proposed
Lagrangian formulation enabled implicit time integration which considerably
increases the computational robustness. We demonstrate that three major forces
determine the dynamics of the cells in an aggregate network; the pili-pili
mediated attractive forces, the steric repulsion forces and the cell-substrate
friction. In the absence of pili-pili mediated forces, the repulsive forces
simply distribute the cells uniformly throughout the domain whereas in the
presence of the pili-pili mediated forces, we observe a phase separation
leading to the formation of micro-colonies. A parametric study has been
carried out to study the influence of various parameters on the cell density
behavior. Our proposed methodology furnishes a general framework for the
continuum modeling of the non-equilibrium dynamics of dense cellular
aggregates. We believe that this contribution provides significant insights
towards the dynamics of cell aggregates which in turn can be exploited to
better understand the behavior of infectious diseases. Further extension of
this work include analysis of different cell species, as well as examining
other biological systems with pronounced cell-matrix interactions.

## Acknowledgment

Soheil Firooz and Paul Steinmann gratefully acknowledge the support provided
by EAM cluster. Also Soheil Firooz would like to thank Hui–Shun Kuan for
fruitful discussions regarding the parametric study. Vasily Zaburdaev would
like to acknowledge the support by Volkswagen foundation “Life?” initiative.

##  Appendix A Detailed derivations of the weak forms

This section provides further details regarding the intermediate steps in the
derivation of the linear momentum balance equation (23) as follows

$\displaystyle\int_{\mathcal{B}_{0}}\nabla_{\boldsymbol{X}}\cdot$
$\displaystyle\boldsymbol{P}^{\text{a}}\cdot\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\nabla_{\boldsymbol{X}}\cdot\boldsymbol{P}^{\text{p}}\cdot\delta\boldsymbol{y}\,\mbox{d}V-\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V$ (54)
$\displaystyle=$
$\displaystyle\int_{\mathcal{B}_{0}}\nabla_{\boldsymbol{X}}\cdot\left(\boldsymbol{P}^{a^{T}}\cdot\delta\boldsymbol{y}\right)\,\mbox{d}V-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{a}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\nabla_{\boldsymbol{X}}\cdot\left(\boldsymbol{P}^{p^{T}}\cdot\delta\boldsymbol{y}\right)\,\mbox{d}V-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V-\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V$
$\displaystyle=$
$\displaystyle\int_{\mathcal{B}_{0}}\left[\boldsymbol{P}^{a^{T}}\cdot\delta\boldsymbol{y}\right]\cdot\boldsymbol{N}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{a}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\left[\boldsymbol{P}^{p^{T}}\cdot\delta\boldsymbol{y}\right]\cdot\boldsymbol{N}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V-\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V$
$\displaystyle=$
$\displaystyle\int_{\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{P}^{\text{a}}\cdot\boldsymbol{N}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{a}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{P}^{\text{p}}\cdot\boldsymbol{N}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V-\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V$
$\displaystyle=$
$\displaystyle\int_{\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{T}^{\text{a}}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{a}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V+\int_{\mathcal{B}_{0}}\delta\boldsymbol{y}\cdot\boldsymbol{T}^{\text{p}}\,\mbox{d}A-\int_{\mathcal{B}_{0}}\boldsymbol{P}^{\text{p}}:\nabla_{\boldsymbol{X}}\delta\boldsymbol{y}\,\mbox{d}V-\int_{\mathcal{B}_{0}}\xi
Jc_{t}\dot{\boldsymbol{y}}\cdot\delta\boldsymbol{y}\,\mbox{d}V=0\,.$

##  Appendix B Detailed derivations of the active stress

In this section, the derivation of the material time derivative of the active
Piola–Kirchhoff stress is provided in detail. To begin, the Lie time
derivative of the Kirchhoff stress is obtained as follows

$\displaystyle\mathcal{L}_{t}\boldsymbol{\tau}^{\text{a}}$
$\displaystyle=\boldsymbol{F}\cdot\left[\dot{\overline{\boldsymbol{F}^{-1}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{F}^{-T}}}\right]\cdot\boldsymbol{F}^{T}=\boldsymbol{F}\cdot\left[\dot{\boldsymbol{F}^{-1}}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{F}^{-T}+\boldsymbol{F}^{-1}\cdot\dot{\boldsymbol{\tau}^{\text{a}}}\cdot\boldsymbol{F}^{-T}+\boldsymbol{F}^{-1}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\dot{\boldsymbol{F}^{-T}}\right]\cdot\boldsymbol{F}^{T}$
(55)
$\displaystyle=\boldsymbol{F}\cdot\left[-\boldsymbol{F}^{-1}\cdot\boldsymbol{l}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{F}^{-T}+\boldsymbol{F}^{-1}\cdot\dot{\boldsymbol{\tau}^{\text{a}}}\cdot\boldsymbol{F}^{-T}-\boldsymbol{F}^{-1}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}^{-T}\right]\cdot\boldsymbol{F}^{T}$
$\displaystyle=\dot{\boldsymbol{\tau}^{\text{a}}}-\boldsymbol{l}\cdot\boldsymbol{\tau}^{\text{a}}-\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{l}^{T}=\dot{\boldsymbol{\tau}^{\text{a}}}-\boldsymbol{l}\cdot{\boldsymbol{\tau}^{\text{a}}}^{T}-\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{l}^{T}=\dot{\boldsymbol{\tau}^{\text{a}}}-2\left[\boldsymbol{l}\cdot\boldsymbol{\tau}^{\text{a}}\right]^{\text{sym}}\,.$

Afterwards, using the relation
$\boldsymbol{S}=\boldsymbol{F}^{-1}\cdot\boldsymbol{\tau}^{\text{a}}\cdot\boldsymbol{F}^{-T}$
on could write

$\displaystyle\dot{\boldsymbol{S}^{\text{a}}}=-\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\left[\boldsymbol{F}\cdot\boldsymbol{S}^{\text{a}}\cdot\boldsymbol{F}^{T}\right]:\boldsymbol{l}^{T}\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}=-\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}\right]\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}\,.$
(56)

Since the Piola–Kirchhoff stress is a symmetric tensor, for an arbitrary
second order tensor $\boldsymbol{A}$ we can write

$\boldsymbol{S}:\boldsymbol{A}=\boldsymbol{S}:\left[\displaystyle\frac{1}{2}\left[\boldsymbol{A}+\boldsymbol{A}^{T}\right]\right]\,.$
(57)

Thus, Eq. (56) can be rewritten as

$\displaystyle\dot{\boldsymbol{S}^{\text{a}}}$
$\displaystyle=-\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}\right]\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}$
(58)
$\displaystyle=-\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\frac{1}{2}\left[\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}\right]+\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}\right]^{T}\right]\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}$
$\displaystyle=-\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{\text{sym}}\cdot\boldsymbol{F}\right]\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}$

Finally, utilizing the relation

$\boldsymbol{F}^{T}\cdot[\boldsymbol{l}]^{\text{sym}}\cdot\boldsymbol{F}=\boldsymbol{F}^{T}\cdot\left[\displaystyle\frac{1}{2}\left[\boldsymbol{l}^{T}+\boldsymbol{l}\right]\right]\cdot\boldsymbol{F}=\displaystyle\frac{1}{2}\left[\boldsymbol{F}^{T}\cdot\boldsymbol{l}^{T}\cdot\boldsymbol{F}+\boldsymbol{F}^{T}\cdot\boldsymbol{l}\cdot\boldsymbol{F}\right]=\displaystyle\frac{1}{2}\left[\dot{\boldsymbol{F}^{T}}\cdot\boldsymbol{F}+\boldsymbol{F}^{T}\cdot\dot{\boldsymbol{F}}\right]=\dot{\boldsymbol{E}}\,,$
(59)

we can derive the fully Lagrangian form of the active Piola–Kirchhoff stress
as

$\dot{\boldsymbol{S}^{\text{a}}}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\dot{\boldsymbol{E}}\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}\,.$
(60)

##  Appendix C Time integration of the active stress

In this section we detail on the time integration technique in order to
calculate the active second Piola–Kirchhoff stress. The non-linear relation
for the active second Piola–Kirchhoff stress reads

$\dot{\boldsymbol{S}^{\text{a}}}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}}\left[\boldsymbol{S}^{\text{a}}:\dot{\boldsymbol{E}}\right]\boldsymbol{S}^{\text{a}}+\boldsymbol{S}^{\text{f}}-k_{\text{off}}\,\boldsymbol{S}^{\text{a}}\,,$
(61)

which could be written in the form

$\displaystyle\frac{\boldsymbol{S}^{\text{a}}_{n+1}-\boldsymbol{S}^{\text{a}}_{n}}{\Delta
t}=-\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}\Delta
t}\boldsymbol{S}^{\text{a}}_{n+1}\left[\boldsymbol{S}^{\text{a}}_{n+1}:\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]+\boldsymbol{S}_{n+1}^{\text{f}}-k_{\text{off}}\boldsymbol{S}^{\text{a}}_{n+1}\,.$
(62)

To linearize this equation, we put all the terms on one side and treat them as
a residuum $\boldsymbol{R}$ that must vanish

$\boldsymbol{R}=\displaystyle\frac{\boldsymbol{S}^{\text{a}}_{n+1}-\boldsymbol{S}^{\text{a}}_{n}}{\Delta
t}+\displaystyle\frac{1}{\ell_{0}p_{0}f^{\text{p}}\Delta
t}\boldsymbol{S}^{\text{a}}_{n+1}\left[\boldsymbol{S}^{\text{a}}_{n+1}:\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]-\boldsymbol{S}_{n+1}^{\text{f}}+k_{\text{off}}\boldsymbol{S}^{\text{a}}_{n+1}\,.$
(63)

The linearization of $\boldsymbol{R}$ reads

$\boldsymbol{R}_{n+1}\approx\text{Lin}\boldsymbol{R}_{n+1}=\boldsymbol{R}_{n}+\displaystyle\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{S}^{\text{a}}}\big{\lvert}_{n}\cdot\Delta\boldsymbol{S}^{\text{a}}_{n}\stackrel{{\scriptstyle!}}{{=}}\boldsymbol{0}\,.$
(64)

The tangent reads

$\displaystyle\text{K}=\frac{\partial\boldsymbol{R}}{\partial\boldsymbol{S}^{\text{a}}}$
$\displaystyle=\left[\frac{1}{\Delta
t}+k_{\text{off}}\right]\text{I}^{\text{sym}}+\frac{1}{\ell_{0}p_{0}f^{\text{p}}\Delta
t}\bigg{[}\left[\boldsymbol{S}^{\text{a}}_{n+1}:\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]\text{I}^{\text{sym}}+\boldsymbol{S}^{\text{a}}_{n+1}\otimes\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\bigg{]}$
(65) $\displaystyle=\left[\frac{1}{\Delta
t}+k_{\text{off}}+\frac{1}{\ell_{0}p_{0}f^{\text{p}}\Delta
t}\left[\boldsymbol{S}^{\text{a}}_{n+1}:\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]\right]\text{I}^{\text{sym}}+\frac{1}{\ell_{0}p_{0}f^{\text{p}}\Delta
t}\bigg{[}\boldsymbol{S}^{\text{a}}_{n+1}\otimes\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\bigg{]}\,,$

with $\text{I}^{\text{sym}}$ being the symmetric fourth-order identity which
reads
$\text{I}^{\text{sym}}=1/2\left[\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right]$.
Finally, an iterative Newton–Raphson scheme is employed to solve for
$\boldsymbol{R}_{n+1}$ and thus, for $\boldsymbol{S}^{\text{a}}_{n+1}$ as

$\left[\boldsymbol{S}^{\text{a}}\right]_{n+1}=\left[\boldsymbol{S}^{\text{a}}\right]_{n}+\left[d\boldsymbol{S}^{\text{a}}\right]\quad\text{with}\quad\left[\Delta\boldsymbol{S}^{\text{a}}\right]_{ij}=-\left[\text{K}^{-1}\right]_{ijkl}\left[\boldsymbol{R}\right]_{kl}\,.$
(66)

To calculate $\text{K}^{-1}$ we use the Sherman–Morrison formula which states
that for an arbitrary fourth-order tensor A that can be written in the form

$\text{A}=\beta\text{B}+\alpha\boldsymbol{C}\otimes\boldsymbol{D}\,,$ (67)

with $\alpha$ and $\beta$ being scalars, $\boldsymbol{C}$ and $\boldsymbol{D}$
being second-order tensors and B being a fourth-order tensor, the inverse of A
reads

$\text{A}^{-1}=\displaystyle\frac{1}{\beta}\text{B}^{-1}-\displaystyle\frac{\alpha}{\beta^{2}+\alpha\beta\boldsymbol{D}:\text{B}^{-1}:\boldsymbol{C}}\left[\text{B}^{-1}:\boldsymbol{C}\otimes\boldsymbol{D}:\text{B}^{-1}\right]\,.$
(68)

To proceed, we write our tangent in the form of Eq. (67) as

$\text{K}=\beta\text{I}^{\text{sym}}+\alpha\bigg{[}\boldsymbol{S}_{n+1}\otimes\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\bigg{]}\,,$
(69)

with

$\alpha=\displaystyle\frac{1}{\ell_{0}\rho^{\text{p}}_{0}f^{\text{p}}\Delta
t}\qquad\text{and}\qquad\beta=\left[\displaystyle\frac{1}{\Delta
t}+k_{\text{off}}+\displaystyle\frac{1}{\ell_{0}\rho^{\text{p}}_{0}f^{\text{p}}\Delta
t}\left[\boldsymbol{S}_{n+1}:\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]\right]\,.$
(70)

Therefore, using the identity
$\text{I}^{\text{sym}^{-1}}=\text{I}^{\text{sym}}$, $\text{K}^{-1}$ reads

$\text{K}^{-1}=\displaystyle\frac{1}{\beta}\text{I}^{\text{sym}}-\displaystyle\frac{\alpha}{\beta^{2}+\alpha\beta\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]:\text{I}^{\text{sym}}:\boldsymbol{S}_{n+1}}\left[\text{I}^{\text{sym}}:\boldsymbol{S}_{n+1}\otimes\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]:\text{I}^{\text{sym}}\right]\,,$
(71)

or in index notation

$\left[\text{K}^{-1}\right]_{ijkl}\\!\\!\\!=\displaystyle\frac{1}{\beta}\left[\text{I}^{\text{sym}}\right]_{ijkl}-\displaystyle\frac{\alpha}{\beta^{2}+\alpha\beta\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]_{mn}:\left[\text{I}^{\text{sym}}\right]_{mnrs}:\left[\boldsymbol{S}_{n+1}\right]_{rs}}\left[\left[\text{I}^{\text{sym}}\right]_{ijmn}:\left[\boldsymbol{S}_{n+1}\otimes\left[\boldsymbol{E}_{n+1}-\boldsymbol{E}_{n}\right]\right]_{mnrs}:\left[\text{I}^{\text{sym}}\right]_{rskl}\right]\,.$
(72)

##  Appendix D Some useful relations

This section provides some useful derivatives in index notation relation which
will be helpful regarding the derivation of the tangents for the finite
element implementation.

$\left[\displaystyle\frac{\partial
J}{\partial\boldsymbol{F}}\right]_{ij}=J\left[\boldsymbol{F}^{-T}\right]_{ij}\,,\qquad\left[\displaystyle\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{F}}\right]_{ijkl}=\delta_{ik}\delta_{jl}\,,\qquad\left[\displaystyle\frac{\partial\boldsymbol{F}^{T}}{\partial\boldsymbol{F}}\right]_{ijkl}=\delta_{il}\delta_{jk}\,,\qquad\left[\displaystyle\frac{\partial\boldsymbol{F}^{-1}}{\partial\boldsymbol{F}}\right]_{ijkl}=-\left[\boldsymbol{F}^{-1}\right]_{ik}\left[\boldsymbol{F}^{-T}\right]_{jl}\,,$
$\left[\displaystyle\frac{\partial\boldsymbol{F}^{-T}}{\partial\boldsymbol{F}}\right]_{ijkl}=-\left[\boldsymbol{F}^{-1}\right]_{il}\left[\boldsymbol{F}^{-T}\right]_{jk}\,,\qquad\left[\displaystyle\frac{\partial\boldsymbol{B}}{\partial\boldsymbol{F}}\right]_{ijkl}=\left[\displaystyle\frac{\partial\boldsymbol{F}^{-1}}{\partial\boldsymbol{F}}\right]_{imkl}\left[\boldsymbol{F}^{-T}\right]_{mj}+\left[\boldsymbol{F}^{-1}\right]_{im}\left[\displaystyle\frac{\partial\boldsymbol{F}^{-T}}{\partial\boldsymbol{F}}\right]_{mjkl}\,,$
$\left[\displaystyle\frac{\partial\boldsymbol{E}}{\partial\boldsymbol{F}}\right]_{ijkl}=\displaystyle\frac{1}{2}\left[\left[\displaystyle\frac{\partial\boldsymbol{F}^{T}}{\partial\boldsymbol{F}}\right]_{imkl}\left[\boldsymbol{F}\right]_{mj}+\left[\boldsymbol{F}^{T}\right]_{im}\left[\displaystyle\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{F}}\right]_{mjkl}\right]\,,\\\
$
$\left[\displaystyle\frac{\partial\boldsymbol{P}}{\partial\boldsymbol{F}}\right]_{ijkl}=\left[\displaystyle\frac{\partial\boldsymbol{F}}{\partial\boldsymbol{F}}\right]_{imkl}\left[\boldsymbol{S}\right]_{mj}+\left[\boldsymbol{F}\right]_{im}\left[\displaystyle\frac{\partial\boldsymbol{S}}{\partial\boldsymbol{F}}\right]_{mjkl}\,,\qquad\left[\displaystyle\frac{\partial[\boldsymbol{g}\otimes\boldsymbol{g}]}{\partial\boldsymbol{g}}\right]_{ijk}\\!\\!\\!\\!=\left[\left[\boldsymbol{I}\right]_{ik}\left[\boldsymbol{g}\right]_{j}+\left[\boldsymbol{g}\right]_{i}\left[\boldsymbol{I}\right]_{jk}\right]\,,$
$\left[\displaystyle\frac{\partial[\nabla_{\boldsymbol{X}}\,\boldsymbol{g}:\boldsymbol{F}^{-T}]}{\partial\nabla_{\boldsymbol{X}}\,\boldsymbol{g}}\right]_{ij}\\!\\!\\!\\!=\left[\boldsymbol{F}^{-T}\right]_{ij}\left[\displaystyle\frac{\partial[\nabla_{\boldsymbol{X}}\,\boldsymbol{g}:\boldsymbol{F}^{-T}]}{\partial\,\boldsymbol{F}}\right]_{ij}\\!\\!\\!\\!=\left[\displaystyle\frac{\partial\boldsymbol{F}^{-T}}{\partial\boldsymbol{F}}\right]_{klij}\left[\nabla_{\boldsymbol{X}}\,\boldsymbol{g}\right]_{kl}\,.$

##  Appendix E Calculation of the derivatives at the element level

In this section a brief example to calculate the derivatives of a scalar field
$\alpha$, a vector field $\boldsymbol{a}$ and a second-order tensor field
$\boldsymbol{A}$ with respect to the nodal values at the element level is
elaborated.

$\left[\displaystyle\frac{\partial\alpha}{\partial\alpha^{J}}\right]=\displaystyle\frac{\partial\left(\alpha^{S}N^{S}\right)}{\partial\alpha^{J}}=\displaystyle\frac{\partial\alpha^{S}}{\partial\alpha^{J}}N^{S}=\delta^{SJ}N^{S}=N^{J}\,,$
$\left[\displaystyle\frac{\partial\nabla_{\boldsymbol{x}}\alpha}{\partial\alpha^{J}}\right]_{i}=\displaystyle\frac{\partial\left[\nabla_{\boldsymbol{x}}\left(\alpha^{S}N^{S}\right)\right]_{i}}{\partial\alpha^{J}}=\displaystyle\frac{\partial\left[\alpha^{S}\nabla_{\boldsymbol{x}}N^{S}_{i}\right]}{\partial\alpha^{J}}=\displaystyle\frac{\partial\alpha^{S}}{\partial\alpha^{J}}\,\nabla_{\boldsymbol{x}}N^{S}_{i}=\delta^{SJ}\,\nabla_{\boldsymbol{x}}N^{S}_{i}=\left[\nabla_{\boldsymbol{x}}N^{J}\right]_{i}\,,$
$\left[\displaystyle\frac{\partial\boldsymbol{a}}{\partial\boldsymbol{a}^{J}}\right]_{ij}=\displaystyle\frac{\partial\left(\boldsymbol{a}^{S}_{i}N^{S}\right)}{\partial\boldsymbol{a}^{J}_{j}}=\displaystyle\frac{\partial\boldsymbol{a}^{S}_{i}}{\partial\boldsymbol{a}^{J}_{j}}N^{S}=\delta_{ij}\delta^{SJ}N^{S}=\delta_{ij}N^{J}\,,$
$\left[\displaystyle\frac{\partial\nabla_{\boldsymbol{x}}\boldsymbol{a}}{\partial\boldsymbol{a}^{J}}\right]_{ijk}=\displaystyle\frac{\partial\left[\nabla_{\boldsymbol{x}}\left(\boldsymbol{a}^{S}_{i}N^{S}\right)\right]_{j}}{\partial\boldsymbol{a}^{J}_{k}}=\displaystyle\frac{\partial\left[\boldsymbol{a}^{S}_{i}\otimes\nabla_{\boldsymbol{x}}N^{S}_{j}\right]}{\partial\boldsymbol{a}^{J}_{k}}=\displaystyle\frac{\partial\boldsymbol{a}^{S}_{i}}{\partial\boldsymbol{a}^{J}_{k}}\,\nabla_{\boldsymbol{x}}N^{S}_{j}=\delta_{ik}\delta^{SJ}\,\nabla_{\boldsymbol{x}}N^{S}_{j}=\delta_{ik}\,\left[\nabla_{\boldsymbol{x}}N^{J}\right]_{j}\,,$
$\left[\displaystyle\frac{\partial\boldsymbol{A}}{\partial\boldsymbol{A}^{J}}\right]_{ijkl}\\!\\!\\!\\!=\displaystyle\frac{\partial\left(\boldsymbol{A}^{S}_{ij}N^{S}\right)}{\partial\boldsymbol{A}^{J}_{kl}}=\displaystyle\frac{\partial\boldsymbol{A}^{S}_{ij}}{\partial\boldsymbol{A}^{J}_{kl}}N^{S}=\delta_{ik}\delta_{jl}\delta^{SJ}N^{S}=\delta_{ik}\delta_{jl}N^{J}\,,$
$\left[\displaystyle\frac{\partial\nabla_{\boldsymbol{x}}\boldsymbol{A}}{\partial\boldsymbol{A}^{J}}\right]_{ijklm}=\displaystyle\frac{\partial\left[\nabla_{\boldsymbol{x}}\left(\boldsymbol{A}^{S}_{ij}N^{S}\right)\right]_{k}}{\partial\boldsymbol{A}^{J}_{lm}}=\displaystyle\frac{\partial\left[\boldsymbol{A}^{S}_{ij}\otimes\nabla_{\boldsymbol{x}}N^{S}_{k}\right]}{\partial\boldsymbol{A}^{J}_{lm}}=\displaystyle\frac{\partial\boldsymbol{A}^{S}_{ij}}{\partial\boldsymbol{A}^{J}_{lm}}\,\nabla_{\boldsymbol{x}}N^{S}_{k}=\delta_{il}\delta_{jm}\delta^{SJ}\,\nabla_{\boldsymbol{x}}N^{S}_{k}=\delta_{il}\delta_{jm}\,\left[\nabla_{\boldsymbol{x}}N^{J}\right]_{k}\,.$

## References

## References

  * [1] H. S. Kuan, W. Pönisch, F. Jülicher, V. Zaburdaev, Continuum Theory of Active Phase Separation in Cellular Aggregates, Phys. Rev. Lett. 126 (2021) 18102.
  * [2] Y. Futaki, I. Amimoto, M. Tanaka, T. Ito, Y. Hirano, Discovery of cell aggregate-inducing peptides, Processes 9 (2021) 538.
  * [3] L. G. Griffith, M. A. Swartz, Capturing complex 3D tissue physiology in vitro, Nat. Rev. Mol. Cell Biol. 7 (2006) 211–224.
  * [4] S. j. Kim, E. M. Kim, M. Yamamoto, H. Park, H. Shin, Engineering Multi-Cellular Spheroids for Tissue Engineering and Regenerative Medicine, Adv. Healthc. Mater. 9 (2020) 1–18.
  * [5] F. Pampaloni, E. G. Reynaud, E. H. K. Stelzer, The third dimension bridges the gap between cell culture and live tissue, Nat. Rev. Mol. Cell Biol. 8 (2007) 839–845.
  * [6] W. Mueller-Kleiser, Multicellular spheroids, J. Cancer Res. Clin. Oncol. 113 (1987) 101–122.
  * [7] T. Eguchi, C. Sogawa, Y. Okusha, K. Uchibe, R. Iinuma, K. Ono, K. Nakano, J. Murakami, M. Itoh, K. Arai, T. Fujiwara, Y. Namba, Y. Murata, K. Ohyama, M. Shimomura, H. Okamura, M. Takigawa, T. Nakatsura, K. i. Kozaki, K. Okamoto, S. K. Calderwood, Organoids with cancer stem cell-like properties secrete exosomes and HSP90 in a 3D nanoenvironment, PLoS One 13 (2018) 1–34.
  * [8] S. Douezan, J. Dumond, F. Brochard-Wyart, Wetting transitions of cellular aggregates induced by substrate rigidity, Soft Matter 8 (2012) 4578–4583.
  * [9] G. Beaune, T. V. Stirbat, N. Khalifat, O. Cochet-Escartin, S. Garcia, V. V. Gurchenkov, M. P. Murrell, S. Dufour, D. Cuvelier, F. Brochard-Wyart, How cells flow in the spreading of cellular aggregates, Proc. Natl. Acad. Sci. U. S. A. 111 (2014) 8055–8060.
  * [10] G. Beaune, G. Duclos, N. Khalifat, T. V. Stirbat, D. M. Vignjevic, F. Brochard-Wyart, Reentrant wetting transition in the spreading of cellular aggregates, Soft Matter 13 (2017) 8474–8482.
  * [11] S. Douezan, K. Guevorkian, R. Naouar, S. Dufour, D. Cuvelier, F. Brochard-Wyarta, Spreading dynamics and wetting transition of cellular aggregates, Proc. Natl. Acad. Sci. U. S. A. 108 (2011) 7315–7320.
  * [12] H. Clevers, Modeling Development and Disease with Organoids, Cell 165 (2016) 1586–1597.
  * [13] L. R. Johnson, Microcolony and biofilm formation as a survival strategy for bacteria, J. Theor. Biol. 251 (2008) 24–34.
  * [14] E. Ben-Jacob, I. Cohen, D. L. Gutnick, Cooperative organization of bacterial colonies: From genotype to morphotype, Annu. Rev. Microbiol. 52 (1998) 779–806.
  * [15] M. T. Armstrong, P. B. Armstrong, Mechanisms of epibolic tissue spreading analyzed in a model morphogenetic system: Roles for cell migration and tissue contractility, J. Cell Sci. 102 (1992) 373–385.
  * [16] D. Bi, X. Yang, M. C. Marchetti, M. L. Manning, Motility-driven glass and jamming transitions in biological tissues, Phys. Rev. X 6 (2016) 1–13.
  * [17] N. E. Freitag, H. S. Seifert, M. Koomey, Characterization of the pilF—pilD pilus‐assembly locus of Neisseria gonorrhoeae, Mol. Microbiol. 16 (1995) 575–586.
  * [18] L. Brossay, G. Paradis, R. Fox, M. Koomey, J. Hebert, Identification, localization, and distribution of the PilT protein in Neisseria gonorrhoeae, Infect. Immun. 62 (1994) 2302–2308.
  * [19] W. J. Todd, G. P. Wray, P. J. Hitchcock, Arrangement of pili in colonies of Neisseria gonorrhoeae, J. Bacteriol. 159 (1984) 312–320.
  * [20] M. Klausen, A. Heydorn, P. Ragas, L. Lambertsen, A. Aaes-Jorgensen, S. Molin, T. Tolker-Nielsen, Biofilm formation by Pseudomonas aeruginosa wild type, flagella and type IV pili mutants, Mol. Microbiol. 48 (2003) 1511–1524.
  * [21] A. F. Imhaus, G. Duménil, The number of Neisseria meningitidis type IV pili determines host cell interaction, EMBO J. 33 (2014) 1767–1783.
  * [22] N. J. Armstrong, K. J. Painter, J. A. Sherratt, A continuum approach to modelling cell-cell adhesion, J. Theor. Biol. 243 (2006) 98–113.
  * [23] D. Drasdo, S. Höhme, A single-cell-based model of tumor growth in vitro: Monolayers and spheroids, Phys. Biol. 2 (2005) 133–147.
  * [24] J. Griffie, R. Peters, D. M. Owen, An agent-based model of molecular aggregation at the cell membrane, PLoS One 15 (2020) 1–17.
  * [25] A. R. A. Anderson, M. A. J. Chaplain, E. L. Newman, R. J. C. Steele, A. M. Thompson, Mathematical Modelling of Tumour Invasion and Metastasis, J. Theor. Med. 2 (2000) 129–154.
  * [26] M. Block, E. Schöll, D. Drasdo, Classifying the expansion kinetics and critical surface dynamics of growing cell populations, Phys. Rev. Lett. 99 (2007) 3–6.
  * [27] J. MOREIRA, A. DEUTSCH, Cellular Automaton Models of Tumor Development: a Critical Review, Adv. Complex Syst. 05 (2002) 247–267.
  * [28] M. S. Alber, M. A. Kiskowski, J. A. Glazier, Y. Jiang, On cellular automaton approaches to modeling biological cells, in: Math. Syst. theory Biol. communi- cation, Financ., Springer, New York, 2002, pp. 1–40.
  * [29] T. Alarcón, H. M. Byrne, P. K. Maini, A mathematical model of the effects of hypoxia on the cell-cycle of normal and cancer cells, J. Theor. Biol. 229 (2004) 395–411.
  * [30] L. Geris, J. M. A. Ashbourn, T. Clarke, Continuum-level modelling of cellular adhesion and matrix production in aggregates, Comput. Methods Biomech. Biomed. Engin. 14 (2011) 403–410.
  * [31] P. Macklin, J. Lowengrub, Nonlinear simulation of the effect of microenvironment on tumor growth, J. Theor. Biol. 245 (2007) 677–704.
  * [32] R. P. Araujo, D. L. S. McElwain, A history of the study of solid tumour growth: The contribution of mathematical modelling, Bull. Math. Biol. 66 (2004) 1039–1091.
  * [33] T. Roose, S. J. Chapman, P. K. Maini, Mathematical models of avascular tumor growth, SIAM Rev. 49 (2007) 179–208.
  * [34] D. Horstmann, K. J. Painter, H. G. Othmer, Aggregation under local reinforcement: From lattice to continuum, Eur. J. Appl. Math. 15 (2004) 545–576.
  * [35] K. Anguige, Multi-phase Stefan problems for a non-linear one-dimensional model of cell-to-cell adhesion and diffusion, Eur. J. Appl. Math. 21 (2010) 109–136.
  * [36] K. Anguige, C. Schmeiser, A one-dimensional model of cell diffusion and aggregation, incorporating volume filling and cell-to-cell adhesion, J. Math. Biol. 58 (2009) 395–427.
  * [37] H. M. Byrne, M. A. J. Chaplain, Modelling the role of cell-cell adhesion in the growth and development of carcinomas, Math. Comput. Model. 24 (1996) 1–17.
  * [38] V. Cristini, J. Lowengrub, Q. Nie, Nonlinear simulation of tumor growth, J. Math. Biol. 46 (2003) 191–224.
  * [39] A. Gerisch, M. A. J. Chaplain, Mathematical modelling of cancer cell invasion of tissue: Local and non-local models and the effect of adhesion, J. Theor. Biol. 250 (2008) 684–704.
  * [40] I. Ramis-Conde, M. A. J. Chaplain, A. R. A. Anderson, Mathematical modelling of cancer cell invasion of tissue, Math. Comput. Model. 47 (2008) 533–545.
  * [41] A. C. Callan-Jones, F. Jülicher, Hydrodynamics of active permeating gels, New J. Phys. 13 (2011) 093027.
  * [42] J. Ranft, M. Basan, J. Elgeti, J. F. Joanny, J. Prost, F. Jülicher, Fluidization of tissues by cell division and apoptosis, Proc. Natl. Acad. Sci. U. S. A. 107 (2010) 20863–20868.
  * [43] F. Jülicher, K. Kruse, J. Prost, J. F. Joanny, Active behavior of the Cytoskeleton, Phys. Rep. 449 (2007) 3–28.
  * [44] M. Popović, A. Nandi, M. Merkel, R. Etournay, S. Eaton, F. Jülicher, G. Salbreux, Active dynamics of tissue shear flow, New J. Phys. 19 (2017) 033006.
  * [45] J. Prost, F. Jülicher, J. F. Joanny, Active gel physics, Nat. Phys. 11 (2015) 111–117.
  * [46] M. H. Köpf, L. M. Pismen, A continuum model of epithelial spreading, Soft Matter 9 (2013) 3727–3734.
  * [47] L. Bao, Z. Zhou, Lattice and continuum models analysis of the aggregation diffusion cell movement, arXiv Prepr. (2018).
  * [48] J. C. Arciero, Q. Mi, M. F. Branca, D. J. Hackam, D. Swigon, Continuum model of collective cell migration in wound healing and colony expansion, Biophys. J. 100 (2011) 535–543.
  * [49] H. Byrne, D. Drasdo, Individual-based and continuum models of growing cell populations: A comparison, J. Math. Biol. 58 (2009) 657–687.
  * [50] L. R. Johnson, Microcolony and biofilm formation as a survival strategy for bacteria, J. Theor. Biol. 251 (2008) 24–34.
  * [51] G. O’Toole, H. B. Kaplan, R. Kolter, Biofilm formation as microbial development, Annu. Rev. Microbiol. 54 (2000) 49–79.
  * [52] N. Hoiby, T. Bjarnsholt, M. Givskov, S. Molin, O. Ciofu, Antibiotic resistance of bacterial biofilms, Int. J. Antimicrob. Agents 35 (2010) 322–332.
  * [53] F. Jin, J. C. Conrad, M. L. Gibiansky, G. C. L. Wong, Bacteria use type-IV pili to slingshot on surfaces, Proc. Natl. Acad. Sci. U. S. A. 108 (2011) 12617–12622.
  * [54] A. Charles-Orszag, E. Lemichez, G. Tran Van Nhieu, G. Duménil, Microbial pathogenesis meets biomechanics, Curr. Opin. Cell Biol. 38 (2016) 31–37.
  * [55] C. Toma, H. Kuroki, N. Nakasone, M. Ehara, M. Iwanaga, Minor pilin subunits are conserved in Vibrio cholerae type IV pili, FEMS Immunol. Med. Microbiol. 33 (2002) 35–40.
  * [56] W. Pönisch, C. A. Weber, G. Juckeland, N. Biais, V. Zaburdaev, Multiscale modeling of bacterial colonies: How pili mediate the dynamics of single cells and cellular aggregates, New J. Phys. 19 (2017) 015003.
  * [57] A. B. Jonsson, D. Llver, P. Falk, J. Pepose, S. Normark, Sequence changes in the pilus subunit lead to variation of Neisseria gonorrhoeae to human tissue, Mol. Microbiol. 14 (1994) 1103–1103.
  * [58] D. Bonazzi, V. Lo Schiavo, S. Machata, I. Djafer-Cherif, P. Nivoit, V. Manriquez, H. Tanimoto, J. Husson, N. Henry, H. Chaté, R. Voituriez, G. Duménil, Intermittent Pili-Mediated Forces Fluidize Neisseria meningitidis Aggregates Promoting Vascular Colonization, Cell 174 (2018) 1–13.
  * [59] W. Pönisch, K. B. Eckenrode, K. Alzurqa, H. Nasrollahi, C. Weber, V. Zaburdaev, N. Biais, Pili mediated intercellular forces shape heterogeneous bacterial microcolonies prior to multicellular differentiation, Sci. Rep. 8 (2018) 16567.
  * [60] H. L. Howie, M. Glogauer, M. So, The N. gonorrhoeae type IV pilus stimulates mechanosensitive pathways and cytoprotection through a pilT-dependent mechanism, PLoS Biol. 3 (2005) 0627–0637.
  * [61] J. S. Mattick, Type IV pili and twitching motility, Annu. Rev. Microbiol. 56 (2002) 289–314.
  * [62] C. Holz, D. Opitz, L. Greune, R. Kurre, M. Koomey, M. A. Schmidt, B. Maier, Multiple pilus motors cooperate for persistent bacterial movement in two dimensions, Phys. Rev. Lett. 104 (2010) 1–4.
  * [63] N. Biais, D. L. Higashi, J. Brujić, M. So, M. P. Sheetz, Force-dependent polymorphism in type IV pili reveals hidden epitopes, Proc. Natl. Acad. Sci. U. S. A. 107 (2010) 11358–11363.
  * [64] J. Eriksson, O. S. Eriksson, L. Maudsdotter, O. Palm, J. Engman, T. Sarkissian, H. Aro, M. Wallin, A. B. Jonsson, Characterization of motility and piliation in pathogenic Neisseria Microbial biochemistry, physiology and metabolism, BMC Microbiol. 15 (2015) 1–13.
  * [65] A. J. Merz, M. So, M. P. Sheetz, Pilus retraction powers bacterial twitching motility, Nature 407 (2000) 98–102.
  * [66] B. Maier, L. Potter, M. So, H. S. Seifert, M. P. Sheetz, Single pilus motor forces exceed 100 pN, Proc. Natl. Acad. Sci. 99 (2003) 16012–16017.
  * [67] V. Zaburdaev, N. Biais, M. Schmiedeberg, J. Eriksson, A. B. Jonsson, M. P. Sheetz, D. A. Weitz, Uncovering the mechanism of trapping and cell orientation during Neisseria gonorrhoeae twitching motility, Biophys. J. 107 (2014) 1523–1531.
  * [68] C. A. Weber, Y. T. Lin, N. Biais, V. Zaburdaev, Formation and dissolution of bacterial colonies, Phys. Rev. E - Stat. Nonlinear, Soft Matter Phys. 92 (2015) 1–8.
  * [69] J. Taktikos, Y. T. Lin, H. Stark, N. Biais, V. Zaburdaev, Pili-induced clustering of N. Gonorrhoeae Bacteria, PLoS One 10 (2015) 1–16.
  * [70] H. S. Kuan, W. Pönisch, F. Jülicher, V. Zaburdaev, see supplemental material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.126.018102.
  * [71] W. Pönisch, C. A. Weber, V. Zaburdaev, How bacterial cells and colonies move on solid substrates, Phys. Rev. E 99 (2019) 1–21.
  * [72] R. Marathe, C. Meel, N. C. Schmidt, L. Dewenter, R. Kurre, L. Greune, M. Alexander Schmidt, M. J. I. Müller, R. Lipowsky, B. Maier, S. Klumpp, Bacterial twitching motility is coordinated by a two-dimensional tug-of-war with directional memory, Nat. Commun. 5 (2014) 1–10.
  * [73] L. Craig, K. T. Forest, B. Maier, Type IV pili: dynamics, biophysics and functional consequences, Nat. Rev. Microbiol. 17 (2019) 429–440.
  * [74] L. Craig, M. E. Pique, J. A. Tainer, Type IV pilus structure and bacterial pathogenicity, Nat. Rev. Microbiol. 2 (2004) 363–378.
  * [75] S. A. Papanicolopulos, A. Zervos, A method for creating a class of triangular C1 finite elements, International Journal for Numerical Methods in Engineering 89 (2012) 1437–1450.
  * [76] P. Fischer, J. Mergheim, P. Steinmann, On the C1 continuous discretization of non-linear gradient elasticity: A comparison of NEM and FEM based on Bernstein–B´ezier patches P., International Journal for Numerical Methods in Engineering 82 (2010) 1282–1307.
  * [77] T. J. R. Hughes, J. A. Cottrell, Y. Bazilevs, Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement, Computer Methods in Applied Mechanics and Engineering 194 (2005) 4135–4195.
  * [78] Y. Bazilevs, V. M. Calo, J. A. Cottrell, J. A. Evans, T. J. R. Hughes, S. Lipton, M. A. Scott, T. W. Sederberg, Isogeometric analysis using T-splines, Computer Methods in Applied Mechanics and Engineering 199 (2010) 229–263.
  * [79] N. Kirchner, P. Steinmann, A unifying treatise on variational principles for gradient and micromorphic continua, Philosophical Magazine 85 (2005) 3875–3895.
  * [80] P. Neff, I. D. Ghiba, A. Madeo, L. Placidi, G. Rosi, A unifying perspective: The relaxed linear micromorphic continuum, Continuum Mechanics and Thermodynamics 26 (2014) 639–681.
  * [81] S. Kaessmair, P. Steinmann, Computational Mechanics of Generalized Continua, in: Calculus of Variations, Springer-Verlag, 2020, pp. 343–356.
  * [82] E. Oñate, A. Franci, J. M. Carbonell, Lagrangian formulation for finite element analysis of quasi-incompressible fluids with reduced mass losses Eugenio, International Journal for Numerical Methods in Fluids 74 (2014) 699–731.
  * [83] M. Cremonesi, S. Meduri, U. Perego, A. Frangi, An explicit Lagrangian finite element method for free-surface weakly compressible flows, Computational Particle Mechanics 4 (2017) 357–369.