# Entropy-driven cell-decision making predicts fluid-to-solid transition in
multicellular systems

Arnab Barua  Department of Systems Immunology and Braunschweig Integrated
Centre of Systems Biology, Helmholtz Centre for Infection Research, Rebenring
56, 38106 Braunschweig, Germany  Simon Syga  Technische Univesität Dresden,
Center for Information Services and High Performance Computing, Nöthnitzer
Straße 46, 01062, Dresden, Germany  Pietro Mascheroni  Department of Systems
Immunology and Braunschweig Integrated Centre of Systems Biology, Helmholtz
Centre for Infection Research, Rebenring 56, 38106 Braunschweig, Germany
Nikos Kavallaris  Department of Mathematical and Physical Sciences, University
of Chester, Chester  Michael Meyer-Hermann  Department of Systems Immunology
and Braunschweig Integrated Centre of Systems Biology, Helmholtz Centre for
Infection Research, Rebenring 56, 38106 Braunschweig, Germany  Institute for
Biochemistry, Biotechnology and Bioinformatics,Technische Universität
Braunschweig, Braunschweig, Germany  Andreas Deutsch  Technische Univesität
Dresden, Center for Information Services and High Performance Computing,
Nöthnitzer Straße 46, 01062, Dresden, Germany  Haralampos Hatzikirou
Department of Systems Immunology and Braunschweig Integrated Centre of Systems
Biology, Helmholtz Centre for Infection Research, Rebenring 56, 38106
Braunschweig, Germany  Technische Univesität Dresden, Center for Information
Services and High Performance Computing, Nöthnitzer Straße 46, 01062, Dresden,
Germany

###### Abstract

Cellular decision making allows cells to assume functionally different
phenotypes in response to microenvironmental cues, without genetic change. It
is an open question, how individual cell decisions influence the dynamics at
the tissue level. Here, we study spatio-temporal pattern formation in a
population of cells exhibiting phenotypic plasticity, which is a paradigm of
cell decision making. We focus on the migration/resting and the
migration/proliferation plasticity which underly the epithelial-mesenchymal
transition (EMT) and the go or grow dichotomy. We assume that cells change
their phenotype in order to minimize their microenvironmental entropy (LEUP:
Least microEnvironmental Uncertainty Principle) and study the impact of the
LEUP-driven migration/resting and migration/proliferation plasticity on the
corresponding multicellular spatio-temporal dynamics with a stochastic cell-
based mathematical model for the spatio-temporal dynamics of the cell
phenotypes. In the case of the go or rest plasticity, a corresponding mean-
field approximation allows to identify a bistable switching mechanism between
a diffusive (fluid) and an epithelial (solid) tissue phase which depends on
the sensitivity of the phenotypes to the environment. For the go or grow
plasticity, we show the possibility of Turing pattern formation for the
"solid" tissue phase and its relation with the parameters of the LEUP-driven
cell decisions.

###### keywords:

Cell-decision making; Phenotypic plasticity; Langevin equations; Mean-field
theory; Least microEnvironmental Uncertainty Principle (LEUP)

; Fluid-to-solid transition; Spatio-temporal patterns

##  1 Introduction

Cellular decision making allows cells to respond to microenvironmental cues
[1]. For instance, if an external ligand from the microenvironment binds to
the cell membrane, it may trigger a chain of biochemical reactions inside the
cell [2]. After the microenvironmental information has been processed, the
cell may respond by changing its properties implying cell differentiation,
proliferation, migration, apoptosis etc. Thus, cellular decision making
depends on intrinsic signal transduction pathways, the genetic cell network
[3], extrinsic stimuli, and molecular noise [4]. Cell-decision making plays a
key role in cell fate determination and the maintenance of phenotypic
plasticity. In the former, cells irreversibly acquire new fates by following a
hierarchical lineage, where pluripotent stem cells in a proper
microenvironment (stem cell niche) differentiate into e. g. bone, muscle,
epithelial and further specialised cells. On the other hand, phenotypic
plasticity allows a reversible adaptation to different microenvironmental
stimuli [5]. While cellular decision making is well-studied for single cells,
it is an open question how individual cell decisions influence the dynamics at
the tissue level.

In this paper, we analyse the implications of phenotypic plasticity on spatio-
temporal pattern formation in the cell population and focus on the
migration/resting and migration/proliferation plasticity. The
migration/resting type of plasticity is observed during the epithelial-
mesenchymal transition (EMT/MET)[6, 7, 8] while the migration/proliferation
plasticity characterizes the go-or-grow (GoG) dichotomy found in different
types of cancers, such as glioblastoma or breast cancer, but also during
normal development [9, 10, 11, 12]. Epithelial cells are immotile, connected
by strong adhesion bonds, exhibit regular shapes and show strong apico-basal
cell polarity. Epithelial and mesenchymal phenotypes can be interpreted as
representing two migration modes, resting and motile, respectively. On the
other hand, GoG refers to the mutual exclusion of cell migration and
proliferation, leading to the existence of proliferative/immotile and
nonproliferative/mobile cell phenotypes. We assume that a Least
microEnvironmental Uncertainty Princple (LEUP) dictates the phenotypic
switching in both types of phenotypic plastic behaviours. The LEUP principle
has been introduced in [13] and assumes that cells change their phenotypes in
order to minimize their local microenvironmental entropy. Formulated in the
language of statistical physics and Bayesian statistics, cells try to optimize
their phenotypic prior to reduce the joint entropy of their own phenotype and
microenvironmental variables. Similar ideas using information theoretic
measures have been proposed also in seminal works by W. Bialek [14].

The manuscript is structured as follows: In the next section (2), we define a
LEUP-driven cell decision model based on the migration/resting plasticity. In
Section 3, we develop a stochastic individual-based model (IBM) for moving and
resting cells where cells can change their phenotype according to the LEUP-
driven cell decision model and show the formation of aggregates of resting
cells in corresponding simulations. In section 4, we derive a macroscopic
mean-field approximation of our microscopic model which allows to explain the
formation of aggregates in the microscopic model. In section 5, we analyze the
stability of the steady states and the pattern formation potential of the
macroscopic go-or grow model. Finally, in section 6, we discuss the biological
implications of our results in terms of multicellular growth and pattern
control.

##  2 A minimal model of LEUP-driven cell migration plasticity

Here, we define a mathematical model for LEUP-based switching between moving
and resting phenotypes. We assume that cells can sense their microenvironment
and can change their own phenotype $X_{i}$ on a domain
$\mathcal{L}\subset\mathbb{R}^{2}$ accordingly, where $X_{i}=0,1$ corresponds
to the resting and migrating state, respectively. We define the
microenvironment of cell $i$ as the number $N_{i}^{0}$ of cells having
phenotype $X_{i}=0$ and by the number $N_{i}^{1}$ of cells having phenotype
$X_{i}=1$. By assuming a maximum cell capacity $N$ of the microenvironment ,
we can define

$N=N_{i}^{1}+N_{i}^{0}+N_{i}^{\phi},$ (1)

where $N_{i}^{\phi}$ are the free spaces/slots. The total number of cells
$N_{T}$ is defined as the sum of cells having phenotype $\left(X_{i}=1\right)$
and cells having phenotype $\left(X_{i}=0\right)$. So,

$N_{T}=N_{i}^{1}+N_{i}^{0}.$ (2)

Cells will select their phenotype in a Bayesian fashion, i.e.

$P\left(X_{i}\mid Y_{i}\right)=\frac{P\left(Y_{i}\mid
X_{i}\right)P\left(X_{i}\right)}{P\left(Y_{i}\right)},$ (3)

where $P\left(Y_{i}\mid X_{i}\right)$ can be interpreted as the probability
that the cell perceives all other cells in its surroundings, and
$P\left(X_{i}\right)$ is the prior probability distribution of the cell’s
phenotypes. The former models the microenvironmental sensing the intrinsic
programming of cellular phenotype. However, sensing other cells and evaluating
$P\left(Y_{i}\mid X_{i}\right)$ entails an energy cost. It is reasonable to
assume that the cell will try to optimize its prior $P\left(X_{i}\right)$ for
the sake of energetic frugality.

Here we define $Y_{i}$ as the extrinsic random variable of the $i$-th cell,
where $Y_{i}=N_{i}^{X_{i}}\mid N_{T},X_{i}=0,1$ represents the number of cells
of phenotype $X_{i}$ given the total local number of cells $N_{T}$. The
conditional probability of having $N_{i}^{1}$ number of cells present in the
microenvironment follows a binomial distribution $\left(\boldsymbol{B}\right)$
(for the derivation details see S.I.):

$P\left(N_{i}^{1}\mid N_{T}\right)=P\left(Y_{i}=N_{i}^{1}\mid
N_{T}\right)=\boldsymbol{B}\left(N_{T},p_{1}\right),$ (4)

where $p_{1}$ is the probability of $N_{i}^{1}$ number of cells having
phenotype $X_{i}$ = 1 out of $N_{T}$ cells:

$p_{1}=\frac{N_{i}^{1}}{N_{T}}.$ (5)

Using statistical mechanical arguments [13], the problem of finding the
probability of the cell’s phenotype is equivalent to finding the prior that
minimizes the entropy of a cell’s phenotype and its surroundings. Let
$S\left(X_{i},Y_{i}\right)$ be the entropy of the cell-surroundings,
$S\left(X_{i}\right)$ the internal entropy of the cell, and $S\left(Y_{i}\mid
X_{i}\right)$ the entropy of the information sensed by the cell. Then, the
entropies are connected by the relation
$S\left(X_{i},Y_{i}\right)=S\left(X_{i}\right)+S\left(Y_{i}\mid X_{i}\right)$.
The optimization problem is finding $P\left(X_{i}\right)$ that minimizes
$S\left(X_{i},Y_{i}\right)$, while making sure that $P\left(X_{i}\right)$ is
normalized; in other words we impose

$\begin{split}\frac{\delta}{\delta
P\left(X_{i}\right)}&\left\\{S\left(X_{i}\right)+\beta\left[\sum
P\left(X_{i}\right)S\left(Y_{i}\mid X_{i}\right)-\bar{S}\left(Y_{i}\mid
X_{i}\right)\right]\right.\\\ &\left.-\lambda\left[\sum
P\left(X_{i}\right)-1\right]\right\\}=0,\end{split}$ (6)

where $\frac{\delta}{\delta P\left(X_{i}\right)}$ is the functional
derivative, $\bar{S}\left(Y_{i}\mid X_{i}\right)$ is the expected ensemble
statistics [13], and $\lambda$ and $\beta$ are Lagrange multipliers. Taking
into account the relations among entropies, eqn. (6) yields

$P\left(X_{i}\right)=\frac{e^{\beta S\left(Y_{i}\mid
X_{i}\right)}}{Z}=\frac{1}{1+e^{\beta\Delta S}},$ (7)

where $Z$ is the normalization factor

$Z=\sum_{\alpha=0,1}e^{\beta S\left(Y_{i}\mid X_{i}=\alpha\right)}=e^{\beta
S(Y_{i}|X_{i}=0)}+e^{\beta S(Y_{i}|X_{i}=1)}.$ (8)

Please note that the parameter $\beta$ has a biological interpretation, since
it quantifies the intensity with which a cell senses and complies to the
microenvironment. To estimate the entropy we use the binomial distribution of
the resting and migrating cells. Given $N_{T}$ cells in the environment and
assuming the frequency of cells at a moving state
$p_{1}=\frac{N^{1}_{i}}{N_{T}}$, the probability for $N_{i}^{1}$ moving cells
is

$P\left(N^{1}_{i}\mid
N_{T}\right)=\binom{N_{T}}{N^{1}_{i}}p_{1}^{N_{i}^{1}}(1-p_{1})^{N_{T}-N^{1}_{i}},$
(9)

and the entropy of the binomial distribution in the last equation is

$S(N^{1}_{i},p_{1})=-\sum_{N^{1}_{i}=0}^{N_{T}}P(N^{1}_{i},p_{1})\ln
P(N^{1}_{i},p_{1}).$ (10)

When there are no cells in the microenvironment the entropy
$S(N^{1}_{i},p_{1})$ is 0. To facilitate the evaluation of the entropy
distribution of the microenvironment, we have used a Gaussian distribution,
which approximates the binomial distribution (for $N_{T}>5$), i.e.

$P(N^{1}_{i},p_{1})\rightarrow{\boldsymbol{N}(N_{T}p_{1},N_{T}p_{1}(1-p_{1}))},$
(11)

for the microenvironment [15]. Now we can evaluate the equilibrium
distribution $P(X_{i}=0)$ as

$\begin{split}P\left(X_{i}=0\right)&=\frac{e^{\beta S\left(Y_{i}\mid
X_{i}=0\right)}}{Z}=\frac{1}{1+e^{\beta\Delta S}}=\\\
&=\frac{1}{1+\left[\frac{N_{i}^{0}\left(N_{i}^{1}-1\right)}{N_{i}^{1}\left(N_{i}^{0}-1\right)}\right]^{\frac{\beta}{2}}}=\frac{1}{1+\left[\frac{\rho_{i}^{0}\left(\rho_{i}^{1}-\frac{1}{V}\right)}{\rho_{i}^{1}\left(\rho_{i}^{0}-\frac{1}{V}\right)}\right]^{\frac{\beta}{2}}},\end{split}$
(12)

with the cell densities $\rho_{i}^{0,1}:=\frac{N_{i}^{0,1}}{V}$ (in the SI
Fig. (12), we show the dependency of $P\left(X_{i}=0\right)$ on resting
density $\rho_{0}$ and $\beta$). Accordingly, the probability of a $i$-th
cell’s phenotype $P(X_{i}=1)$ reads

$\begin{split}P\left(X_{i}=1\right)&=\frac{e^{\beta S\left(Y_{i}\mid
X_{i}=1\right)}}{Z}=1-\frac{1}{1+e^{\beta\Delta S}}=\\\
&=1-\frac{1}{1+\left[\frac{N_{i}^{0}\left(N_{i}^{1}-1\right)}{N_{i}^{1}\left(N_{i}^{0}-1\right)}\right]^{\frac{\beta}{2}}}=\frac{\left[\frac{\rho_{i}^{0}\left(\rho_{i}^{1}-\frac{1}{V}\right)}{\rho_{i}^{1}\left(\rho_{i}^{0}-\frac{1}{V}\right)}\right]^{\frac{\beta}{2}}}{1+\left[\frac{\rho_{i}^{0}\left(\rho_{i}^{1}-\frac{1}{V}\right)}{\rho_{i}^{1}\left(\rho_{i}^{0}-\frac{1}{V}\right)}\right]^{\frac{\beta}{2}}}.\end{split}$
(13)

Please note that, the difference between the microenvironmental entropy of the
$i$-th cell is

$\Delta
S=S(Y_{i}|X_{i}=1)-S(Y_{i}|X_{i}=0)=\frac{1}{2}\ln{\left[\frac{N_{i}^{0}\left(N_{i}^{1}-1\right)}{N_{i}^{1}\left(N_{i}^{0}-1\right)}\right]},$
(14)

$N_{i}^{0}$ and $N_{i}^{1}$ are the number of cell of phenotype ($X_{i}$ = 0)
and phenotype ($X_{i}$ = 1) in the microenvironment of $i$-th cell,
respectively (See S.I. for details).

(a)

(b)

Figure 1: Schematic diagram of (a) the microenvironment of a motile cell,
where $\ell$ is the sensing radius. The sensed number of cells is proportional
to the volume of the sensed microenvironment, i.e. $N\propto\ell^{d}$, where d
is the dimension. (b) Transitions in the “Go or Grow” model [9]. The switching
between motile $\left(X_{i}=1\right)$ and resting $\left(X_{i}=0\right)$
phenotype is shown, where the transtion probabilities are defined as $W_{10}$
and $W_{01}$. The proliferation rate is defined by $r$.

Fig. (1(a)), shows a sketch of our model. In short, we make the following
assumptions
$\left(A1\right)\hskip 14.22636pt$ The microenvironment is defined by the
numbers of two phenotypes only.
$\left(A2\right)\hskip 14.22636pt$ We have assumed a binomial distribution for
the occurrence of the two phenotypes.
$\left(A3\right)\hskip 14.22636pt$ Cells are making decisions at a fast time
scale. This justifies to assume an equilibrium distribution for the different
phenotypes.

In this section, we have defined a minimal LEUP-driven cell decision making
model and derived the corresponding phenotypic steady states. Based on this,
we subsequently develop a cell-based model to understand the resulting
multicelluar spatiotemporal dynamics.

##  3 An individual-based model (IBM) for LEUP-driven cell migration

We define a discrete stochastic, spatio-temporal IBM that incorporates the
phenotypic switch dynamics according to LEUP. To this end, we model the
movement of single cells with Langevin equations. Langevin equations are well-
suited to model cell-cell interactions and cell migration [16, 17].

Our Langevin’s equation is defined on a domain
$\mathcal{L}\subset\mathbb{R}^{2}$ with periodic boundary conditions. We
define an interaction radius $\ell\in\mathbb{R}$ around the $i$-th cell at
position $\mathbf{x_{i}}\in\mathcal{L}$. The expected interaction volume is
$V\propto\ell^{2}$. The time evolution of our model is defined by the
following rules:

  1. (R1)

Cells change their phenotype by sensing their microenvironment within the
interaction radius $\ell$ according to LEUP.

  2. (R2)

Moving cells change their orientations randomly (random walk).

  3. (R3)

Once cells become migratory they move with a constant speed $\bar{v}$.

The above are translated to the following Langevin’s equations

$\begin{split}&\frac{\mathrm{d}\mathbf{x_{i}}}{\mathrm{d}t}=\bar{v}\mathbf{v_{i}}\quantity(\theta_{i}),\\\
&\frac{\mathrm{d}\theta_{i}}{\mathrm{d}t}=\frac{1}{\bar{v}}\xi_{i}^{\theta}(t),\\\
&\frac{\mathrm{d}p_{i}}{\mathrm{d}t}=-\frac{1}{\tau}(p_{i}-p_{i}^{eq}),\\\
&v_{i}=\begin{cases}\bar{v},&\mbox{\quad if\quad}X_{i}=1,\\\ 0,&\mbox{\quad
if\quad}X_{i}=0,\end{cases}\end{split}$ (15)

where $\mathbf{v_{i}}=\quantity(\cos\theta_{i},\sin\theta_{i})^{T}$ is the
direction of movement of cell $i$. For the temporal evolution of the
probability $p_{i}$ of the motile state, we use the BGK
(Bhatnagar–Gross–Krook) operator technique [18]. We assume that the
probability $p_{i}$ evolves weakly out of its equilibrium probability
$p_{i}^{eq}$, which is the LEUP steady state probability $P(X_{i}=1)$ (see eq.
(13)). In turn, the parameter $\tau$ is the relaxation time towards the
corresponding probability distribution $p_{i}^{eq}$. Here noise is assumed to
have a zero-mean, white noise term, which has the statistical properties
$\left\langle\xi^{\theta}_{i}(t)\right\rangle=0$ and
$\left\langle\xi^{\theta}_{i}\left(t_{1}\right)\xi^{\theta}_{j}\left(t_{2}\right)\right\rangle=2D_{\theta}\delta\left(t_{1}-t_{2}\right)\delta_{ij}$,
where $t_{1}$ and $t_{2}$ are two time points, $D_{\theta}\in\mathbb{R}_{+}$
is the angular diffusion coefficient, $\delta(t)$ is the Dirac delta, and
$\delta_{ij}$ is the Kronecker delta. Parameter $\bar{v}$ is the constant
speed of every motile cell.

We simulate the Langevin model, on a two-dimensional domain, with a varying
interaction radius $\ell$, mean cell density and sensitivity $\beta$. We
assume an initial state that is approximately homogeneous in space. For
sufficiently large densities and sensitivities, we observe aggregation
patterns, as shown in Figure (2). To quantify cell clustering, we calculate
the radial distribution function for different sensitivities and different
interaction radii as shown in Figures (3(a)) and (3(b)). The key observations
from the simulation study are:

Figure 2: Clustering of resting cells in the IBM. Snapshot of an example
realization at $t=100$, where the total number of cells is fixed at 1000 and
$\beta=40$. Resting cells are shown in blue and migrating cells are marked by
red arrows indicating their direction of movement.

  1. 1.

There is a critical threshold for the parameter $\beta>0$, where clustering of
resting cells occurs (Fig. 3(a)). This is rather expected since a sufficient
sensing of the microenvironment is required.

  2. 2.

For an intermediate interaction radius $\ell$, we observe cell clusters as
shown in Fig. 3(b). Very high $\ell$ corresponds to a large number of sensed
cells. This leads to also equal steady state probabilities, as the entropy
difference in the microenvironment, associated with the single cell states
becomes negligible. On the other hand, the lower bound of $\ell$ is expected,
since enough sampling size of cellular microenvironment is required to induce
aggregation patterns.

  3. 3.

In Fig. 3(d), we show the phase diagram of the system, when parameter $\beta$
and the average density are varied. In particular, we observe that there is a
parametric regime where clustering behavior is emerging. Interestingly, high
densities do not always imply patterning and require higher values of $\beta$
to support cluster emergence (see Fig. 3(c)). Finally, as expected, low
densities also reduce the area of the patterning regime for low $\beta$.

(a)

(b)

(c)

(d)

Figure 3: Quantification of cell clustering. (a) Radial distribution function
for different values of $\beta$, for a fixed interaction radius $\ell=6$. The
number of cells was fixed at 1000, corresponding to a mean cell density of
0.4. (b) Radial distribution function for varying interaction radius $\ell$ at
fixed sensitivity $\beta=100$. There is an optimal sensing radius for a given
mean density, so that interaction radii that are too large or too low do not
lead to aggregation, indicated by an almost flat radial distribution, see
purple and blue lines. (c) Radial distribution function for different mean
densities, at $\beta=20$ and $\ell=4$. When the mean density becomes too
large, no clustering can be observed. (d) Phase space diagram for cluster
pattern formation. Here we define clustering when $\max_{r}g(r)>1.09$ (see
text for explaination). The interaction radius is fixed at $\ell=4$. All data
points in (a - d) correspond to the mean of 20 independent simulations.

##  4 Mean-field approximation

In this section, we derive a continuous approximation of the aforementioned
discrete model using a mean-field approach. The goal is to shed light of the
pattern formation mechanisms, i.e. cell clustering, as observed in the
microscopic simulations. The main idea of the mean-field approximation is to
replace the description of many-particle interactions by a single particle
description based on an average or effective interaction. Thereby, any multi-
particle problem can be replaced by an effective description, that can be
stated in the form of an ordinary (ODE) or partial differential equation
(PDE). In order to proceed, we will first treat the switch dynamics and the
migration process separately.

###  4.1 Mean-field description of the phenotypic switch dynamics

Let’s consider cell $i$ at position $\mathbf{x_{i}}$ with states $X_{i}=0,1$.
The cell changes its state in dependence on the local microenvironment
according to the LEUP with rates $W_{01},W_{10}$ and the probabilities to be
in either states are denoted by $p_{0,1}$. We assume that the system is always
close to the steady state, so the master equation reads

$\derivative{t}p_{0}=W_{10}p_{1}-W_{01}p_{0}\approx 0,$ (16)

and $p_{0,1}$ are given by eq. 7. By rearranging the terms, we obtain for the
switching rates

$\frac{W_{10}}{W_{01}}=\frac{p_{0}}{p_{1}}.$ (17)

For simplicity, we set $W_{01}=p_{1}$ where the transition probability towards
motile phenotypes equals to the moving steady state probability $p_{1}$. This
coincides with the “detailed balance condition”. The steady state
probabilities $p_{0,1}$ depend on the number of cells $N_{i}^{0,1}$ in the
respective phenotypes in the microenvironment. Cell $i$ at position
$\mathbf{x_{i}}$ senses

$N_{i}^{X}(\mathbf{x_{1},\dots,x_{i},\dots,x_{N}},t):=\sum_{j=1}^{N}\xi\quantity(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t),t)\delta_{X_{j},X}\,,\quad
X=0,1.$ (18)

Here, we sum over all cells $j$, and
$\xi\quantity(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t),t)$ is a Boolean stochastic
variable that serves as the sensing function of a cell at $\mathbf{x_{i}}$. It
depends on the distance between cells and time, and $\xi=1$ indicates that the
cell $j$ at $\mathbf{x_{j}}$ is sensed by cell $i$ at position
$\mathbf{x_{i}}$. To match the IBM, we assume that

$\xi(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t),t)=\Theta\quantity(\ell-\absolutevalue{\mathbf{x}_{j}(t)-\mathbf{x_{i}}(t)}),\quad\ell>0,$
(19)

with the Heaviside step function $\Theta(x)$, so that all cells in a ball of
radius $\ell$ around $\mathbf{x_{i}}$ are sensed.

To proceed, we apply a mean-field approximation to calculate the expected
switching rate $\expectationvalue{W_{10}(N_{i}^{0},N_{i}^{1})}\approx
W_{01}\quantity(\expectationvalue{N_{i}^{0}},\expectationvalue{N_{i}^{1}})$.
Note that we dropped the dependence on space and time for better readability.
Let $\phi_{i}(X,\mathbf{x},t)$ denote the probability of finding cell $i$ in a
small volume $\differential V$ around $\mathbf{x}$ at time $t$ with the
phenotype $X=0,1$. Formally we have

$\phi_{i}(X,\mathbf{x},t)=\expectationvalue{\delta\quantity(\mathbf{x}-\mathbf{x_{i}}(t))\delta_{X,X_{i}}},$
(20)

where the average is the ensemble average, and the total densities of
resting/migrating cells are

$\displaystyle\Phi(\mathbf{x},t)$
$\displaystyle:=\sum_{i}\phi_{i}(0,\mathbf{x},t),$ (21)
$\displaystyle\Psi(\mathbf{x},t)$
$\displaystyle:=\sum_{i}\phi_{i}(1,\mathbf{x},t).$ (22)

For the expected density of sensed cells in the microenvironment around
$\mathbf{x}$ we obtain

$\displaystyle\rho_{0}(\mathbf{x},t):=\frac{\expectationvalue{N_{i}^{0}(\mathbf{x},t)}}{V}$
$\displaystyle=\frac{1}{V}\int_{\mathcal{B}_{\ell}(\mathbf{x})}\Phi(\mathbf{y},t)\,\differential
V,$ (23)
$\displaystyle\rho_{1}(\mathbf{x},t):=\frac{\expectationvalue{N_{i}^{1}(\mathbf{x},t)}}{V}$
$\displaystyle=\frac{1}{V}\int_{\mathcal{B}_{\ell}(\mathbf{x})}\Psi(\mathbf{y},t)\,\differential
V,$ (24)

where $V\propto\ell^{d}$, and where we have already used eq. eq. 19.
Consequently, we obtain the approximate switching rate

$\bar{W}_{10}=p_{0}(\rho_{0},\rho_{1}),$ (25)

and the switch dynamics is (dropping dependencies on space and time for
simplicity)

$\partialderivative{t}\phi_{i}(0)=p_{0}\phi_{i}(1)-p_{1}\phi_{i}(0).$ (26)

Summing over all cells, we obtain for the total density of resting cells

$\partialderivative{t}\Phi=p_{0}\Psi-p_{1}\Phi.$ (27)

Note that this is a non-local, non-linear set of PDEs, which are difficult to
treat analytically. However, we can further simplify our analysis by making
another approximation, assuming that the sensing radius $\ell$ is much smaller
than the total domain size. Then, we can replace the non-local sensing
function $\xi\quantity(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t),t)$ by a local
delta distribution

$\xi\quantity(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t),t)\to\delta\quantity(\mathbf{x_{j}}(t)-\mathbf{x_{i}}(t)).$
(28)

In this case, the expected number of sensed cells in the microenvironment
around $\mathbf{x}$ simply becomes the local density of the respective
phenotype

$\displaystyle\rho_{0}(\mathbf{x},t):=\expectationvalue{N_{i}^{0}(\mathbf{x},t)}$
$\displaystyle=\frac{\Phi(\mathbf{x},t)}{\ell^{d}},$ (29)
$\displaystyle\rho_{1}(\mathbf{x},t):=\expectationvalue{N_{i}^{1}(\mathbf{x},t)}$
$\displaystyle=\frac{\Psi(\mathbf{x},t)}{\ell^{d}}.$ (30)

Finally, the rate $\bar{W}_{10}$ reduces to

$\bar{W}_{10}=p_{0}(\rho_{0},\rho_{1})=\frac{1}{1+e^{\beta\Delta
S}}=\frac{1}{1+\left[\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right]^{\frac{\beta}{2}}},$
(31)

and the reversed transition probability $\bar{W}_{01}$ is

$\bar{W}_{01}=p_{1}(\rho_{0},\rho_{1})=1-\frac{1}{1+e^{\beta\Delta
S}}=1-\frac{1}{1+\left[\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right]^{\frac{\beta}{2}}},$
(32)

where the volume is defined as $V=\ell^{d}$.

###  4.2 Mean-field description of the cell migration process

In this section, we derive the macroscopic equation in two spatial dimensions
for the motile cell population. As before we can write the corresponding
stochastic Langevin’s equations as

$\begin{split}&\derivative{t}\mathbf{x_{i}}=\bar{v}\mathbf{v_{i}}\quantity(\theta_{i}),\\\
&\derivative{t}\theta_{i}=\frac{1}{\bar{v}}\xi_{i}^{\theta}(t).\\\
\end{split}$ (33)

This process can be considered as a special kind of active Brownian motion. In
this kind of Langevin’s equation, the stochastic force creates variations of
orientation. According to [19] we can derive the corresponding Fokker-Planck
equation for migrating cells using adiabatic elimination and averaging it over
different noise realizations obtaining the following diffusion equation

$\frac{\partial\rho_{1}(x,t)}{\partial
t}=\frac{\bar{v}^{4}}{2D_{\theta}}\nabla^{2}\rho_{1}(x,t).$ (34)

Here we define $D_{1}$ as diffusion coefficient which is
$\frac{\bar{v}^{4}}{2D_{\theta}}$.

###  4.3 Coupling migration and switching dynamics

Combining our results in the previous sections, we can easily formulate a
system of PDEs

$\displaystyle\partial_{t}\rho_{0}(x,t)$ $\displaystyle=\nu
E(\rho_{0},\rho_{1})$ (35) $\displaystyle\partial_{t}\rho_{1}(x,t)$
$\displaystyle=D_{1}\laplacian{\rho_{1}}-\nu E(\rho_{0},\rho_{1})$ (36)
$\displaystyle E(\rho_{0},\rho_{1})$
$\displaystyle=\bar{W}_{10}(\rho_{0},\rho_{1})\rho_{1}-\bar{W}_{01}(\rho_{0},\rho_{1})\rho_{0}$
(37)

where $E(\rho_{0},\rho_{1})$ is the phenotypic exchange term and $\nu$ is the
corresponding timescale ratio of the switching and diffusion processes, i.e.
$\nu=\frac{\tau_{S}}{\tau_{D}}$. To ensure the numerical consistency of the
above system, we assume that resting cells diffuse in a very slow manner i.e.
$D_{0}\ll D_{1}$, which results in the following reaction-diffusion system of
equations

$\displaystyle\partial_{t}\rho_{0}(\mathbf{x},t)$
$\displaystyle=D_{0}\laplacian{\rho_{0}}+\nu E(\rho_{0},\rho_{1}),$ (38)
$\displaystyle\partial_{t}\rho_{1}(\mathbf{x},t)$
$\displaystyle=D_{1}\laplacian{\rho_{1}}-\nu E(\rho_{0},\rho_{1}).$ (39)

##  5 Spatio-temporal dynamics of the LEUP-driven migration/proliferation
plasticity

In this section, we study the mean-field approximation of the aforementioned
stochastic plasticity dynamics for different regimes of the sensitivity
$\beta$ and the interaction radius $\ell$. In particular, we initially study
the system dynamics, when cells have a very large interaction radius ($\ell\gg
1)$ and then a finite one. Finally, we study also the special case $\beta=0$,
i.e. cells decide independently of their microenvironment.

###  5.1 Large interaction radius case: existence of a critical sensitivity

Here, we focus on the very large interaction radius system ($\ell\gg 1\implies
V\gg 1$). Although, this parameter regime is not biologically relevant, it is
very instructive since it allows us to derive analytical estimates for our
systems dynamics. Here, we generalize the macroscopic system by adding
proliferation dynamics (logistic growth). Our phenotypic switch dynamics are
recapitulated by setting the proliferation rate to zero, i.e. $r=0$. The full
system reads

$\begin{split}&\frac{\partial\rho_{0}}{\partial
t}=D_{0}\nabla^{2}{\rho_{0}}+\nu
E_{0}\left(\rho_{0},\rho_{1}\right)+r\rho_{0}\left(1-\rho_{1}-\rho_{0}\right),\\\
&\frac{\partial\rho_{1}}{\partial t}=D_{1}\nabla^{2}{\rho_{1}}-\nu
E_{0}\left(\rho_{0},\rho_{1}\right).\end{split}$ (40)

In turn, we conduct a non-dimensionalization of eqns. (40) to identify the
variables. Moreover, this helps us to gain a knowledge about the relationships
between the different model parameters. By assuming that the system size is
fixed at $L$ the non-dimensional quantities read

$\begin{split}x^{*}&=\frac{x}{L}\implies\frac{\partial^{2}}{\partial
x^{*2}}=L^{2}\frac{\partial^{2}}{\partial x^{2}},\\\
t^{*}&=\frac{D_{0}t}{L^{2}}\implies\frac{\partial}{\partial
t^{*}}=\frac{L^{2}}{D_{0}}\frac{\partial}{\partial t},\\\
\gamma&=\frac{L^{2}\nu}{D_{0}},\\\ D&=\frac{D_{1}}{D_{0}},\\\
r^{{}^{\prime}}&=\frac{r}{\nu}.\end{split}$ (41)

In the limit $V\gg 1$, the eqn. (40) can also be written as

$\begin{split}&\frac{\partial\rho_{0}}{\partial
t^{*}}=\nabla^{*2}{\rho_{0}}+\gamma\left(\left(\rho_{1}-\rho_{0}\right)\left(\frac{1}{2}-\tilde{\beta}\frac{\left(\rho_{1}+\rho_{0}\right)}{\rho_{1}\rho_{0}}\right)+r^{{}^{\prime}}\rho_{0}\left(1-\rho_{1}-\rho_{0}\right)\right),\\\
&\frac{\partial\rho_{1}}{\partial
t^{*}}=D\nabla^{*2}{\rho_{1}}-\gamma\left(\left(\rho_{1}-\rho_{0}\right)\left(\frac{1}{2}-\tilde{\beta}\frac{\left(\rho_{1}+\rho_{0}\right)}{\rho_{1}\rho_{0}}\right)\right).\end{split}$
(42)

To understand the behaviour of the system at long times, we conduct a fixed
point analysis. Initially we assume a well-stirred system, i.e. no spatial
interactions. Then, eqs. (42) can be written as coupled non-linear ODEs which
have three fixed points

$\begin{split}&\left(\rho_{0}^{A*},\rho_{1}^{A*}\right)=\left(\frac{1}{2},\frac{1}{2}\right),\\\
&\left(\rho_{0}^{B*},\rho_{1}^{B*}\right)=\left(\frac{1}{2}\left(1+\sqrt{1-8\tilde{\beta}}\right),\frac{1}{2}\left(1-\sqrt{1-8\tilde{\beta}}\right)\right),\\\
&\left(\rho_{0}^{C*},\rho_{1}^{C*}\right)=\left(\frac{1}{2}\left(1-\sqrt{1-8\tilde{\beta}}\right),\frac{1}{2}\left(1+\sqrt{1-8\tilde{\beta}}\right)\right).\end{split}$
(43)

The above imply a pitchfork bifurcation for the
$\tilde{\beta}=\frac{\beta}{8V}$ parameter, i.e. there exists a critical value
$\tilde{\beta}_{c}$ that introduces the bistable state. From eqs.(43), we can
easily deduce that the critical sensitivity value

$1-8\tilde{\beta_{c}}=0\Longleftrightarrow\beta_{c}=V.$ (44)

This is an acceptable approximation even for the finite interaction radius
system (see next section). In the following, we analyze if the system is able
to produce spatial patterns. Applying linear stability analysis for the
spatially resolved system eqs. (42), we can deduce that no pattern formation
is possible (for details check next section). Any perturbations to the
homogeneous state lead always to a spatially homogeneous steady state. When
$r=0$ this result is consistent with our findings in discrete IBM simulations
where very large interaction radii do not confer any clustering, as shown in
Fig. 3(b).

###  5.2 The finite interaction radius case: emergence of a bistable switch
between "fluid" and "solid" tissue phases and pattern formation

Now we turn to the full system for intermediate interaction radius also
(assuming proliferation). This implies a finite interaction volume $V$ and for
analytical feasibility we are interested in the Gaussian approximation of the
switching probabilities and their corresponding mean-field terms. The full
system of PDEs assuming also proliferation reads

$\begin{split}\frac{\partial\rho_{0}}{\partial
t^{*}}&=\nabla^{*2}{\rho_{0}}+\gamma\left(\frac{\rho_{1}-\rho_{0}\left(\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right)^{\frac{\beta}{2}}}{1+\left(\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right)^{\frac{\beta}{2}}}+r^{{}^{\prime}}\rho_{0}\left(1-\rho_{1}-\rho_{0}\right)\right),\\\
\frac{\partial\rho_{1}}{\partial
t^{*}}&=D\nabla^{*2}{\rho_{1}}-\gamma\frac{\rho_{1}-\rho_{0}\left(\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right)^{\frac{\beta}{2}}}{1+\left(\frac{\rho_{0}\left(\rho_{1}-\frac{1}{V}\right)}{\rho_{1}\left(\rho_{0}-\frac{1}{V}\right)}\right)^{\frac{\beta}{2}}}.\end{split}$
(45)

Since closed expression of steady states are not analytically feasible, we
obtain the bifurcation diagram numerically, see Fig (4(a)).

(a)

(b)

Figure 4: (a) Bifurcation diagram for moving cells in the finite interaction
radius limit with respect to $\beta$, where proliferation rate $r$ is fixed to
0 and $V^{-1}$ at 0.02. The bifurcation diagram is symmetric for resting cells
i.e. $\rho_{0}^{*}$. (b) The critical point value of the bifurcation diagram
(where the solution splits in two branches) is plotted against the inverse
volume $\frac{1}{V}$.

We observe the existence of a supercritical pitchfork bifurcation and the
existence of a critical $\beta_{c}$. For $\beta\geq\beta_{c}$, the systems
depart from balanced state ($\rho_{0}=\rho_{1}=\frac{1}{2})$ to the
coexistence of a “fluid” (most cells migrate) and a “solid” (most cells are
resting) phase. The switch is controlled by the perturbation on the ratio of
migratory and resting cells. Interestingly, we can compare the analytic
estimate of $\beta_{c}=V$, from the $V\gg 1$ case, with the one calculated for
the finite $V$ system. Fig. (4(b)) shows that the infinite system
approximation provides an upper bound for $\beta_{c}$, which is not too far
from the real value of finite systems.

The existence of a critical $\beta_{c}$ with respect to the two phases is
evident in the IBM simulations as well. In particular if we quantify the ratio
of stationary over motile cells, we observe a similar behavior of the critical
sensitivity for increasing $V$, since it decreases (see Fig.13 in SI ).
However, we cannot conduct a strict quantitative comparison since our
bifurcation analysis does not involve any diffusion as opposed to the IBM
simulation.

Figure 5: Turing space diagram of sensitivity and proliferation rate, i.e.
$\frac{\beta}{2}$ vs. $r^{{}^{\prime}}$, where different values of
$\log{(d_{c})}$ and $k_{c}$ have been plotted. Here we denote $a=\frac{1}{V}$.

In turn, we apply linear stability analysis to identify parameter regimes that
promote diffusion driven pattern formation or Turing instability. To analyze
the Turing instability [20] we have to find the system’s steady state, (i.e.
when diffusion is not present in the systems of eqn.(40)). It has been shown
that $d=\frac{D_{1}}{D_{0}}\gg 1$ is a necessary condition for the emergence
of spatially heterogeneous solution i.e. patterns. Now, we can write the
system of PDEs (i.e eqn.(40)) in a generalized matrix form

$\begin{split}\frac{\partial\boldsymbol{\rho}}{\partial
t}&=\boldsymbol{D}\nabla^{2}\boldsymbol{\rho}+\gamma\boldsymbol{R\rho},\hskip
7.11317pt\boldsymbol{D}=\begin{pmatrix}1&0\\\ 0&d\end{pmatrix},\hskip
7.11317pt\\\ \boldsymbol{R}&=\begin{pmatrix}\frac{\partial
E_{0}\left(\rho_{0},\rho_{1}\right)}{\partial\rho_{0}}+r^{{}^{\prime}}\left(1-\rho_{1}-\rho_{0}\right)-r^{{}^{\prime}}\rho_{0}&\frac{\partial
E_{0}\left(\rho_{0},\rho_{1}\right)}{\partial\rho_{1}}-r^{{}^{\prime}}\rho_{0}\\\
-\frac{\partial
E_{0}\left(\rho_{0},\rho_{1}\right)}{\partial\rho_{0}}&-\frac{\partial
E_{0}\left(\rho_{0},\rho_{1}\right)}{\partial\rho_{1}}\end{pmatrix}_{\left(\rho_{0}^{*},\rho_{1}^{*}\right)},\\\
&=\begin{pmatrix}f_{u}&f_{v}\\\
g_{u}&g_{v}\end{pmatrix}_{\left(\rho_{0}^{*},\rho_{1}^{*}\right)},\end{split}$
(46)

where $\boldsymbol{\rho}$ is defined as
$\left(\rho_{0}-\rho_{0}^{*},\rho_{1}-\rho_{1}^{*}\right)^{T}$ and
$\boldsymbol{R}$ is the Jacobian at $\rho^{*}=(\rho_{0}^{*},\rho_{1}^{*})$.
Using the Turing conditions of instability [20], we found patterns when $N$ is
finite in zero-flux boundary conditions. We have checked all the Turing
instability conditions:

$\begin{split}&f_{u}+g_{v}<0,\hskip 22.76228ptf_{u}g_{v}-f_{v}g_{u}>0,\\\
&Df_{u}+g_{v}>0,\hskip
22.76228pt(Df_{u}+g_{v})^{2}-4D(f_{u}g_{v}-f_{v}g_{u})>0.\end{split}$ (47)

Interestingly, only when the density of resting cells is larger than that of
the moving cells, patterns are formed under Turing instability conditions. In
order to investigate the system’s potential to exhibit pattern formation, we
check the range of validity of the Turing instability criteria (47).
Diffusion-driven instability conditions are satisfied for a large portion of
the parameter space. In turn for these parameters, we calculate the critical
diffusion coefficient $d_{c}$. For values $d>d_{c}$ we are able to observe
patterns, The existence of $d_{c}$ is associated with the existence of a
critical wavenumber $k_{c}$ [20].

$\begin{split}&d_{c}=\frac{-2(2f_{v}g_{u}-f_{u}g_{v})\pm\sqrt{\left(2(2f_{v}g_{u}-f_{u}g_{v})\right)^{2}-4f_{u}^{2}g_{v}^{2}}}{2f_{u}^{2}},\\\
&k_{c}=\gamma\left[\frac{Det(A)}{d_{c}}\right]^{\frac{1}{2}}.\end{split}$ (48)

In Fig. (5) we identify the parameter regime that allows us to observe
patterns and then corresponding $d_{c}$ and $k_{c}$.

By fixing the initial conditions to a cosine wave, we observe in one dimension
the existence of regular spikes as shown in Fig (6). Please note that the
exact pattern is sensitive to the initial conditions. In turn, we simulate our
system in 2D and for the same parameters and we observe patterning in the form
of dots. Finally, we investigate if the type of patterns changes for
variations in parameters $\beta$ and $r^{{}^{\prime}}$ (see Fig. (7)). If
$d>d_{c}$ then we need a large domain or size of the system to observe the
patterns. In the 2D case, we have observed the discoidal patterns of resting
phenotype which resemble the 1D case. So, the radius of the circles of the
patterns are increased if we fix the domain size. Moreover, we observe that in
both dimensions (i.e. 1D and 2D) the critical spatial frequency increases with
decreasing $r^{{}^{\prime}}$.

Figure 6: Densities of proliferating (a) and migrating (b) cells of the 1D
mean-field simulation after 50000 time steps. Figure 7: One-dimensional
simulation of the mean-field approximation for finite interaction radius
system. In (b-c) densities of proliferating cells are plotted over time
corresponding to two points (TS1 and TS2) inside the Turing space shown in
(a). Here $\beta^{{}^{\prime}}=\frac{\beta}{2}$ and $a$ is defined as the
inverse volume $\frac{1}{V}$. Figure 8: Two-dimensional simulation of the
mean-field approximation for finite interaction radius system. The patterns
are obtained for parameters that correspond to the first point (TS1) in the
Turing space, as displayed in Fig.(7)

In comparison with the discrete IBM simulations, we have to state that the
simulation clusters can be identified by the discoidal mean-field patterns.
Under this statement, we find that:

  1. 1.

There exists a critical $\beta$ that allows for the emergence of patterns as
in the IBM simulations.

  2. 2.

For very low and very high interaction radius $\ell$, we observe no patterns,
which is consistent with our IBM results (see Fig.3(b)).

###  5.3 Microenvironment independent phenotypic switching leads to
uncontrolled growth dynamics

Finally we investigate the case $\beta$ = 0 where cells do not sense their
microenvironment and so we end up with the following system,

$\begin{split}&\frac{\partial\rho_{0}}{\partial
t^{*}}=\nabla^{*2}{\rho_{0}}+\gamma\left(\frac{\left(\rho_{1}-\rho_{0}\right)}{2}+r^{{}^{\prime}}\rho_{0}\left(1-\rho_{1}-\rho_{0}\right)\right),\\\
&\frac{\partial\rho_{1}}{\partial
t^{*}}=D\nabla^{*2}{\rho_{1}}-\gamma\frac{\left(\rho_{1}-\rho_{0}\right)}{2}.\end{split}$
(49)

We can do a fixed point analysis similar to equation (42). Assuming no spatial
dynamics, we can find two fixed points i.e.

$\begin{split}\left(\rho_{0}^{P*},\rho_{1}^{P*}\right)&=\left(\frac{1}{2},\frac{1}{2}\right),\\\
\left(\rho_{0}^{Q*},\rho_{1}^{Q*}\right)&=(0,0).\end{split}$ (50)

For further details see Fig. 9 in S.I. We can clearly show that the fixed
point (0,0) is unstable and $\left(\frac{1}{2},\frac{1}{2}\right)$ is a saddle
point. We can also write the coupled PDE equations in a single PDE as well, if
we consider $\rho_{0}=\rho_{1}=\tilde{\rho}$. Then we obtain

$\frac{\partial\tilde{\rho}}{\partial
t^{*}}=\frac{(1+D)}{2}\nabla^{*2}{\tilde{\rho}}+\frac{r}{2}\tilde{\rho}\left(1-2\tilde{\rho}\right).$
(51)

Now we can clearly see that the eqn.(51) is similar to the Fisher-Kolmogorov
equation [20]. It is known that the Fisher-Kolmogorov equation does not
exhibit any pattern formation instabilities. The latter observation is
consistent with our discrete IBM simulations, where no clustering is observed.

##  6 Discussion

Recently, the Least microEnvironmental Uncertainty Principle (LEUP) has been
proposed as an organization principle for cell decision-making in
multicellular systems. In this paper, we apply this principle to shed light on
the effects of phenotypic plasticity for tissue dynamics with a mathematical
model. We focus on two types of plasticity: a go or rest, and a go or grow
phenotypic dichotomy, which play key roles in important processes in
biological development and pathological situations as cancer invasion. We
assume that in any given spatial mesoscopic sample the presence of cell
phenotypes follows a binomial distribution. Using this assumption, we are able
to calculate the microenvironmental entropy and the LEUP-driven probability
distribution of each phenotype. On the basis of this distribution, we defined
an appropriate microscopic stochastic model for the spatio-temporal dynamics
of both phenotypes, and in turn derived the corresponding mean-field
description resulting in a system of coupled reaction-diffusion equations. The
main results of our study are: (i) in the case of go or rest plasticity, there
exists a supercritical pitchfork bifurcation that defines a switch between a
"fluid" and a "solid" phase, and (ii) in the case of go or grow plasticity,
for the "solid" phase, we can derive conditions for the emergence of Turing
patterns. Interestingly, the existence of Turing patterns requires a critical
LEUP sensitivity to the microenvironment and a minimum interaction range.

Our model assumes binary transitions between two discrete phenotypes, which is
an extreme simplification of biological reality. Certainly, one could include
a continuous spectrum of motile/proliferative phenotypes that is expected to
imply even richer spatio-temporal dynamics. Indeed, assuming a continuous
state (velocity) space would potentially lead to further interesting
bifurcations, such as the metastable EMT state as found in [21, 22], or
complex spatiotemporal patterns as indicated in [23]. But the very abstracted
case chosen here is still instructive to indicate how entropy-driven
phenotypic decisions can lead to particular types of spatio-temporal pattern
formation.

Interestingly, our LEUP-driven IBM is an extension of a Vicsek-type model
[24], formulated in the context of self-propelled particles [25]. Our model
exhibits a novel collective behavior when compared to the past published
results from Vicsek-type of models. In particular, to our knowledge it is the
first time to produce with such models Turing patterns, i.e. dynamics clusters
of non-motile cells of specific characteristic wavelength. Typically, in
Viscek-type models we may observe moving clusters of swirling cells (e.g. the
milling Viscek model) but never static ones.

At this point, let us focus on the biological assumptions and implications of
our study. The molecular regulatory mechanisms involved in EMT or GoG remain
largely unknown, where the latter can be viewed as an EMT with proliferation
constrained to the epithelial/resting phase. Here, we assume that the
phenotypic regulation of both mechanisms is based on the minimization of
microenvironmental entropy in physiological tissues, which allows us to
predict the multicellular spatiotemporal dynamics. This assumption is
supported by the fact that healthy physiological processes, biological
development or processes like wound healing, where EMT/MET or GoG are present,
typically lead to an ordered (low entropy) tissue from a disordered initial
condition. On the other hand, deregulation of EMT and GoG have been already
identified as pivotal elements in invading cancers [26, 9], where genetic and
phenotypic heterogeneity, characterized by high entropy, is a key
characteristic. Assuming a LEUP-driven migration/proliferation phenotypic
regulation allows us to understand how cells control multicellular dynamics in
terms of growth and patterning.

Implications in multicellular growth control: The central finding of our study
is associated with the bifurcation diagram of the LEUP sensitivity parameter
$\beta$ in Fig. 4(a). As stated before, the parameter $\beta$ quantifies how
prone cells are in sensing and responding to their microenvironmental stimuli.
When $\beta=0$ cells migrate and proliferate independent of their
microenvironmental cues, which corresponds to one of the cancer hallmarks
[27]. In this case, the systems grows uncontrollably, resembling a cancerous
tissue [27]. The resulting Fisher-Kolmogorov macroscopic behavior has been
prototypically used to model invading tumours [28]. On the other hand, by
adding a death process for any motility state in the GoG model, we can
recapitulate the Allee effect (bistability between extinction and growth) as
found in Boettger et al [29]. A cell population at the “fluid” state will go
extinct (motile cells do not proliferate but still have a probability to die,)
where as systems in the “solid” state will always grow until carrying
capacity. Therefore, in the bistability regime cell sensing properties lead to
multicellular growth control.

Implications in multicellular pattern control: Increased cell sensing $\beta$
represents the physiological tissue dynamics, since it allows the system to
control its behavior. By tuning $\beta$ and the ratio of motile/resting cells,
the system exhibits a bistable behavior between a "fluid"/mesenchymal-type and
a “solid”/epithelial-type tissue phase. This kind of tissue level switch is of
utmost importance in physiological processes such as wound healing or
embryogenesis [30]. After a tissue injury, the healing is characterized by a
"fluid" diffusive expansion of fibroblast cells, that adopt a migratory
phenotype via EMT [31]. After covering the wound, the "solid" phase emerges as
cells stop the migration program and proliferate to finalize tissue repair. In
the abscence of proliferation, the bistable switch from the "solid" to the
"fluid" state could potentially explain the jamming phase transition observed
in epithelial colonies, under EGF modulation and Rab5a knock-out [32].
Typically, the "solid" multicellular phase is prone to the emergence of
pattern formation, which frequently occurs in physiological epithelial tissues
[31, 33]. When EMT is combined with a Notch-Delta cell-cell communication then
epithelial/immotile cell clusters emerge [22], as observed in our IBM and
mean-field simulations. Adding proliferation in the GoG model, the type and
the size of such emerging patterns require a tight regulation of
microenvironmental sensing $\beta$ and proliferation rate $r$, as indicated by
Fig. 11(a), since the critical wavelength $k_{c}$ depends on the ratio
$\frac{\beta}{r^{{}^{\prime}}}$. Such Turing patterns are in agreement with
previous GoG studies [34], where Turing patterns where emerging.

In conclusion, our study shows how individual LEUP-driven cell decisions the
dynamics at the tissue level and how knowledge of collective cell decision-
making can be used to control of growth and pattern.

### Acknowledgment

AB thanks the International Graduate School of HZI, Braunschweig. HH, MMH and
PM gratefully acknowledge the funding support of the Helmholtz Association of
German Research Centers—Initiative and Networking Fund for the project on
Reduced Complexity Models (ZT-I-0010). HH and PM acknowledge the funding
support of MicMode-I2T (01ZX1710B) by the Federal Ministry of Education and
Research (BMBF). HH is supported by MulticellML (01ZX1707C) of the Federal
Ministry of Education and Research (BMBF) and the Volkswagenstiftung for the
its support within the "Life?" programm (96732). SS acknowledges financial
support co-financed by the European Social Fund (ESF) and tax funds in
accordance with the budget adopted by the members of the Saxon State
Parliament. Part of the current work was inspired and initiated when NK was
visiting the Helmholtz Centre for Infection Research and he would like to
express his gratitude for the warm hospitality of the institute. NK would also
like to acknowledge financial support from the Faculty of Science and
Engineering of University of Chester. The authors thank the Centre for
Information Services and High Performance Computing (ZIH) at TU Dresden for
providing an excellent infrastructure.

## References

  * [1] Clive G. Bowsher and Peter S. Swain  “Environmental sensing, information transfer, and cellular decision-making”  In _Curr. Opin. Biotechnol._ 28 Elsevier Ltd, 2014, pp. 149–155  DOI: 10.1016/j.copbio.2014.04.010
  * [2] B. Alberts  “Molecular Biology of the Cell”  New York: W. W, 2015
  * [3] Laura Prochazka, Yaakov Benenson and Peter W. Zandstra  “Synthetic gene circuits and cellular decision-making in human pluripotent stem cells”  In _Curr. Opin. Syst. Biol._ 5.September Elsevier Ltd, 2017, pp. 93–103  DOI: 10.1016/j.coisb.2017.09.003
  * [4] Gábor Balázsi, Alexander Oudenaarden and James J. Collins  “Cellular decision making and biological noise: from microbes to mammals”  In _Cell_ 144.6, 2011, pp. 910–925  DOI: 10.1016/j.cell.2011.01.030
  * [5] Trevor D. Price, Anna Qvarnström and Darren E. Irwin  “The role of phenotypic plasticity in driving genetic evolution”  In _Proc. R. Soc. London. Ser. B Biol. Sci._ 270.1523, 2003, pp. 1433–1440  DOI: 10.1098/rspb.2003.2372
  * [6] Raghu Kalluri and Robert A. Weinberg  “The basics of epithelial-mesenchymal transition”  In _The J. of Cli. Inves._ 120.5 The American Society for Clinical Investigation, 2010, pp. 1786–1786  DOI: 10.1172/JCI39104C1
  * [7] Derek C. Radisky  “Epithelial-mesenchymal transition”  In _J. Cell Sci._ 118.19 The Company of Biologists Ltd, 2005, pp. 4325–4326  DOI: 10.1242/jcs.02552
  * [8] Jean Paul Thiery  “Epithelial–mesenchymal transitions in development and pathologies”  In _C. Opi. Cell Bio._ 15.6, 2003, pp. 740–746  DOI: https://doi.org/10.1016/j.ceb.2003.10.006
  * [9] H. Hatzikirou et al.  “’Go or Grow’: the key to the emergence of invasion in tumour progression?”  In _Math. Med. Biol._ 29.1, 2010, pp. 49–65  DOI: 10.1093/imammb/dqq011
  * [10] Livnat Jerby et al.  “Metabolic associations of reduced proliferation and oxidative stress in advanced breast cancer”  In _Cancer Res._ 72.22, 2012, pp. 5712–5720  DOI: 10.1158/0008-5472.CAN-12-2215
  * [11] Keith S. Hoek et al.  “In vivo switching of human melanoma cells between proliferative and invasive states”  In _Cancer Res._ 68.3, 2008, pp. 650–656  DOI: 10.1158/0008-5472.CAN-07-2491
  * [12] Abraham Q. Kohrman and David Q. Matus  “Divide or Conquer: Cell Cycle Regulation of Invasive Behavior”  In _T. Cell Bio._ 27.1 Elsevier Ltd, 2017, pp. 12–25  DOI: 10.1016/j.tcb.2016.08.003
  * [13] Haralampos Hatzikirou  “Statistical mechanics of cell decision-making: the cell migration force distribution”  In _J. Mech. Behav. Mater._ 27.1-2, 2018, pp. 1–7  DOI: 10.1515/jmbm-2018-0001
  * [14] William Bialek  “Biophysics: Searching for principles”  In _Biophys. Search. Princ._ Princeton University Press, 2012
  * [15] George E.. Box, William G. Hunter and Joanna S. Hunter  “Statistics for Experimenters: An Introduction to Design, Data Analysis, and Model Building”, 1979
  * [16] Vincent Hakim and Pascal Silberzan  “Collective cell migration: a physics perspective”  In _Rep. on Prog. in Phys._ 80.7 IOP Publishing, 2017, pp. 076601  DOI: 10.1088/1361-6633/aa65ef
  * [17] M. Schienbein and H. Gruler  “Langevin equation, Fokker-Planck equation and cell migration”  In _Bull. Math. Bio._ 55.3, 1993, pp. 585–608  DOI: 10.1007/BF02460652
  * [18] P.. Bhatnagar, E.. Gross and M. Krook  “A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems”  In _Phys. Rev._ 94 American Physical Society, 1954, pp. 511–525  DOI: 10.1103/PhysRev.94.511
  * [19] S. Milster, J. Nötel, I.. Sokolov and L. Schimansky-Geier  “Eliminating inertia in a stochastic model of a micro-swimmer with constant speed”  In _The Euro. Phys. Jour. Sp. Topics_ 226.9, 2017, pp. 2039–2055  DOI: 10.1140/epjst/e2017-70052-8
  * [20] J.. Murray  “Mathematical Biology II. Spatial Models and Biomedical Applications” 18, Interdisciplinary Applied Mathematics  New York, NY: Springer New York, 2003, pp. 1067  DOI: 10.1007/b98869
  * [21] Mohit Kumar Jolly et al.  “Towards elucidating the connection between epithelial - mesenchymal transitions and stemness”  In _J. Roy. Soc. Int._ 11.101, 2014, pp. 20140962  DOI: 10.1098/rsif.2014.0962
  * [22] Marcelo Boareto et al.  “Jagged–Delta asymmetry in Notch signaling can give rise to a Sender/Receiver hybrid phenotype”  In _PNAS_ 112.5 National Academy of Sciences, 2015, pp. E402–E409  DOI: 10.1073/pnas.1416287112
  * [23] Arnab Barua, Josue M. Nava-Sedeño and Haralampos Hatzikirou  “A least microenvironmental uncertainty principle (LEUP) as a generative model of collective cell migration mechanisms”  In _bioRxiv_ Cold Spring Harbor Laboratory, 2019  DOI: 10.1101/404889
  * [24] Tamás Vicsek et al.  “Novel type of phase transition in a system of self-driven particles”  In _Phys. rev. lett._ 75.6 APS, 1995, pp. 1226
  * [25] Robert Grossmann, Lutz Schimansky-Geier and Pawel Romanczuk  “Self-propelled particles with selective attraction–repulsion interaction: from microscopic dynamics to coarse-grained theories”  In _New J. Phys._ 15.8 IOP Publishing, 2013, pp. 085014
  * [26] Jacqueline Banyard and Diane R. Bielenberg  “The role of EMT and MET in cancer dissemination” PMID: 26291767  In _Conn. Tis. Res._ 56.5 Taylor & Francis, 2015, pp. 403–413  DOI: 10.3109/03008207.2015.1060970
  * [27] Douglas Hanahan and Robert A. Weinberg  “The Hallmarks of Cancer”  In _Cell_ 100.1 Elsevier, 2000, pp. 57–70  DOI: 10.1016/S0092-8674(00)81683-9
  * [28] E Dahlman and Y Watanabe  “SU-F-T-109: A Shortcoming of the Fisher-Kolmogorov Reaction-Diffusion Equation for Modeling Tumor Growth”  In _Med. Phys._ 43.6Part14, 2016, pp. 3486–3487  DOI: 10.1118/1.4956245
  * [29] Katrin Böttger et al.  “An Emerging Allee Effect Is Critical for Tumor Initiation and Persistence”  In _PLOS Comp. Bio._ 11.9, 2015, pp. e1004366  DOI: 10.1371/journal.pcbi.1004366
  * [30] Guislaine Barriere et al.  “Epithelial Mesenchymal Transition: a double-edged sword”  In _Clinic.Trans. Med._ 4.1, 2015, pp. 14  DOI: 10.1186/s40169-015-0055-4
  * [31] Le Yan and Dapeng Bi  “Multicellular Rosettes Drive Fluid-solid Transition in Epithelial Tissues”  In _Phys. Rev. X_ 9 American Physical Society, 2019, pp. 011029  DOI: 10.1103/PhysRevX.9.011029
  * [32] Chiara Malinverno et al.  “Endocytic reawakening of motility in jammed epithelia”  In _Nature Mat._ 16.5, 2017, pp. 587–596  DOI: 10.1038/nmat4848
  * [33] Evan Heller and Elaine Fuchs  “Tissue patterning and cellular mechanics”  In _J. Cell Bio._ 211.2, 2015, pp. 219–231  DOI: 10.1083/jcb.201506106
  * [34] Kara Pham et al.  “Density-dependent quiescence in glioma invasion: instability in a simple reaction–diffusion model for the migration/proliferation dichotomy”  In _J. Bio. Dyn._ Taylor & Francis, 2011, pp. 1–18  DOI: 10.1080/17513758.2011.590610

## Supplementary Information

Symbol | Explanation
---|---
$X_{i}$ | Phenotype of $i$-th cell
$Y_{i}$ | Phenotype of $i$-th neighbourhood cells
$N_{i}^{0}$ | Number of cells in $i$-th cell’s microenvironment having phenotype $(X_{i}=0)$
$N_{i}^{1}$ | Number of cells in $i$-th cell’s microenvironment having phenotype $(X_{i}=1)$
$N_{i}^{\phi}$ | Number of free slots in $i$-th cell’s microenvironment
$V$ | Total capacity
$N_{T}$ | Total number of phenotypes
$\ell$ | Radius of the microenvironment
$\nu$ | Exchange rate
$r$ | Growth rate
$\beta$ | LEUP sensitivity
$\rho_{0}(x)$ | Mean density of the resting cells at position x
$\rho_{1}(x)$ | Mean density of the migratory cells at position x
$D_{0}$ | Diffusion coefficient of resting cells
$D_{1}$ | Diffusion coefficient of migratory cells

## S.1. Calculation of microenvironmental entropy

We assume that a maximum number of N cells is present inside the cell’s
microenvironment, where l is the radius of the microenvironment (Fig.1(a)).
Since there is the possibility of free slots, if there are less than N cells
in the cell neighbourhood, we have started by assuming a trinomial
distribution for the $i$-th cell is $P\left(N^{1}_{i},N^{\phi}_{i}\right)$
where the number $N_{i}^{0}$ of cells having phenotype $X_{i}=0$ and by the
number $N_{i}^{1}$ of cells having phenotype $X_{i}=1$. In addition, a number
$N_{i}^{\phi}$ of empty slots are included inside the microenvironment. The
joint probability is defined by

$P\left(N^{1}_{i},N^{\phi}_{i}\right)=\frac{N!}{N^{1}_{i}!N^{\phi}_{i}!\left(N-N^{\phi}_{i}-N^{1}_{i}\right)!}p^{N^{1}_{i}}\theta^{N^{\phi}_{i}}\left(1-p-\theta\right)^{\left(N-N^{\phi}_{i}-N^{1}_{i}\right)},$
(52)

where $p$ and $\theta$ are the probabilities of having a number $N_{i}^{1}$ of
cells with phenotype $(X_{i}=1)$ out of $N$ cells and having a number
$N_{i}^{\phi}$ of free slots. The conditional probability of having a number
$N_{i}^{1}$ of cells present in the microenvironment given a number
$N_{i}^{\phi}$ of free slots is

$\begin{split}&P\left(N_{i}^{1}\mid
N_{i}^{\phi}\right)=\frac{P\left(N_{i}^{1},N_{i}^{\phi}\right)}{P\left(N_{i}^{\phi}\right)}\\\
&=\frac{\frac{N!}{N^{1}_{i}!N^{\phi}_{i}!\left(N-N^{\phi}_{i}-N^{1}_{i}\right)!}p^{N^{1}_{i}}\theta^{N^{\phi}_{i}}\left(1-p-\theta\right)^{\left(N-N^{\phi}_{i}-N^{1}_{i}\right)}}{\frac{N!}{N^{\phi}_{i}!\left(N-N^{\phi}_{i}\right)!}\theta^{N^{\phi}_{i}}\left(1-\theta\right)^{\left(N-N^{\phi}_{i}\right)}}\\\
&=\binom{N-N^{\phi}_{i}}{N^{1}_{i}}\left(\frac{p}{1-\theta}\right)^{N^{1}_{i}}\left(1-\frac{p}{1-\theta}\right)^{N-N^{1}_{i}-N^{\phi}_{i}}\\\
&=\boldsymbol{B}\left(N-N_{i}^{\phi},\frac{p}{1-\theta}\right)\\\
&=\boldsymbol{B}\left(N_{T},\frac{p}{1-\theta}\right)\end{split}$ (53)

Now, we use the form of binomial distribution to calculate the entropy of the
Binomial distribution. Please note that $P(N_{i}^{1}\mid N_{i}^{\phi})$ can be
written as $P(N_{i}^{1}\mid N_{T})$ due to Binomial distribution.

$S\left(N^{1}_{i}\mid
N_{T}\right)=S\left(Y_{i}\right)=\frac{1}{2}\log_{2}\left(2\pi
e\frac{N_{T}p}{1-\theta}\left(1-\frac{p}{1-\theta}\right)\right),$ (54)

where

$\begin{split}&\frac{p}{1-\theta}=\frac{N^{1}_{i}}{N}\frac{1}{1-\frac{N^{\phi}_{i}}{N}}=\frac{N^{1}_{i}}{N-N^{\phi}_{i}}=\frac{N_{i}^{1}}{N_{T}},\\\
&1-\frac{p}{1-\theta}=\frac{N^{0}_{i}}{N-N^{\phi}_{i}}=\frac{N_{i}^{0}}{N_{T}}.\end{split}$
(55)

According to LEUP we have to evaluate the microenvironmental entropy of
$S\left(Y_{i}\mid X_{i}=0\right)$ and $S\left(Y_{i}\mid X_{i}=1\right)$ to
calculate the probability of the internal states $\left(X_{i}=0\right)$ and
$\left(X_{i}=1\right)$

$\begin{split}&S\left(Y_{i}\mid X_{i}=1\right)=\frac{1}{2}\log_{2}\left(2\pi
e\left(\frac{\left(N_{T}-1\right)\left(N^{1}_{i}-1\right)}{\left(N_{T}-1\right)}\right)\left(\frac{N^{0}_{i}}{\left(N_{T}-1\right)}\right)\right),\\\
&S\left(Y_{i}\mid X_{i}=0\right)=\frac{1}{2}\log_{2}\left(2\pi
e\left(\frac{\left(N_{T}-1\right)\left(N^{1}_{i}\right)}{\left(N_{T}-1\right)}\right)\left(\frac{N^{0}_{i}-1}{\left(N_{T}-1\right)}\right)\right).\end{split}$
(56)

From the Gaussian approximation we can write the entropy difference as

$\begin{split}&\Delta S=S(Y_{i}|X_{i}=1)-S(Y_{i}|X_{i}=0),\\\
&=\frac{1}{2}\ln{\left[\frac{N_{i}^{0}\left(N_{i}^{1}-1\right)}{N_{i}^{1}\left(N_{i}^{0}-1\right)}\right]}.\\\
\end{split}$ (57)

## S.2. Calculation of cell proliferation rate

In our go-or-grow model only resting cells are allowed to proliferate but not
the moving cells. The growth rate depends in a specific way on the number of
moving and resting cells in the microenvironment. In particular, we assume
that the per- capita growth rate for resting cells is a linearly decreasing
function of $N_{i}^{0}$ and $N_{i}^{1}$, and is also decreasing with the
number of migratory cells and a constant per-capita death rate $d_{1}$.
Accordingly,

$\displaystyle\frac{d\langle N_{i}^{0}\rangle}{dt}=\left(\langle
W^{+}\rangle_{N_{i}^{0},N_{i}^{1}}-\langle
W^{-}\rangle_{N_{i}^{0},N_{i}^{1}}\right)$ (58)
$\displaystyle=h_{1}\rho_{0}-q\left(\rho_{0}+\rho_{1}\right)\rho_{0}-d_{1}\rho_{0}$
$\displaystyle=r\rho_{0}\left(1-\rho_{0}-\rho_{1}\right)$

where $\frac{q}{r}=1$ and $r=h_{1}-d_{1}$.

## S.4. Phase space diagram for null sensitivity case

Figure 9: Phase space diagram of $\rho_{0}^{*}$ and $\rho_{1}^{*}$ for $r=1$

where two fixed points are marked by red circles.

## S.5. Turing space for pattern formation

(a) $a$ is fixed at 0.05. (b) $a$ is fixed at 0.065. (c) $a$ is fixed at 0.08.
Figure 10: Phase space diagram of $\beta$ vs. $r^{{}^{\prime}}$ for different
values of $a$ (a-c), where the Turing space is marked in blue. $a$ is defined
by $\frac{1}{V}$.

## S.6. Critical wavelengths inside the Turing space for pattern formation

(a) $a$ is fixed at 0.05. (b) $a$ is fixed at 0.065. (c) $a$ is fixed at 0.08.
Figure 11: Phase space diagram of $\beta$ vs. $r^{{}^{\prime}}$ for different
values of $a$ (a-c), where critical wavelengths $(\omega_{c})$ have been
plotted inside the Turing space. $a$ is defined by $\frac{1}{V}$.

## S.7. Plot of resting probability

Figure 12: Plot of resting probability with respect to $\beta$ and $\rho_{0}$.
We have fixed the value of $\frac{1}{V}$ to 0.01.

## S.8. Critical sensitivity vs. inverse capacity graph from IBM

Figure 13: Critical sensitivity with respect to the inverse capacity of the
IBM over 5 simulations. Throughout the simulations, we kept the total density
at 0.2 and the total number of cells was 500.


# Aesthetics Driven Autonomous Time-Lapse Photography Generation by Virtual
and Real Robots

Xiaobo Gao  Beihang UniversityBeijingChina<EMAIL_ADDRESS>,  Qi Kuang
0000-0002-1675-1197 Beihang UniversityBeijingChina<EMAIL_ADDRESS>,  Xin
Jin  0000-0003-3873-1653 Beijing Electronic Science and Technology
InstituteBeijingChina<EMAIL_ADDRESS>,  Bin Zhou  0000-0003-2091-1090
Beihang UniversityBeijingChina<EMAIL_ADDRESS>,  Boyan Dong  Beijing
Electronic Science and Technology InstituteBeijingChina<EMAIL_ADDRESS>and
Xunyu Wang  Beijing Electronic Science and Technology InstituteBeijingChina
<EMAIL_ADDRESS>
###### Abstract.

Time-lapse photography is employed in movies and promotional films because it
can reflect the passage of time in a short time and strengthen the visual
attraction. However, since it takes a long time and requires the stable
shooting, it is a great challenge for the photographer. In this article, we
propose a time-lapse photography system with virtual and real robots. To help
users shoot time-lapse videos efficiently, we first parameterize the time-
lapse photography and propose a parameter optimization method. For different
parameters, different aesthetic models, including image and video aesthetic
quality assessment networks, are used to generate optimal parameters. Then we
propose a time-lapse photography interface to facilitate users to view and
adjust parameters and use virtual robots to conduct virtual photography in a
three-dimensional scene. The system can also export the parameters and provide
them to real robots so that the time-lapse videos can be filmed in the real
world. In addition, we propose a time-lapse photography aesthetic assessment
method that can automatically evaluate the aesthetic quality of time-lapse
video. The experimental results show that our method can efficiently obtain
the time-lapse videos. We also conduct a user study. The results show that our
system has the similar effect as professional photographers and is more
efficient.

time-lapse video, aesthetics

††journal: TOG††ccs: Human-centered computing Systems and tools for
interaction design††ccs: Computing methodologies Computational
photography††ccs: Computing methodologies Image and video acquisition

##  1\. Introduction

With the popularity of smartphones and cameras, more and more people are
demanding the high-quality photography. Therefore, many research works focus
on image and video shooting to provide useful guidance (Rawat and Kankanhalli,
2016; Rawat et al., 2019; Lou et al., 2021). Time-lapse photography is one of
the special photographic techniques because of its ability to present the low-
speed changes in a scene with the high playback speed (Bennett and McMillan,
2007). However, photographers have to spend hours or even days shooting to get
the final time-lapse video clips. And during the shooting, the camera needs to
be stationary or move slowly and steadily. The photographers need to adjust
the camera pose frequently, which is time-consuming and laborious for the
photographers. As a result, robots are more suitable for this kind of time-
consuming work that requires stable movement because of the stable pan-tilt
and precise control system.

However, it is difficult to make the robots shoot autonomously. Successful
photography usually requires the photographers to have rich shooting
experience and strong aesthetic ability. How to determine a series of
parameters and let the robots shoot automatically or assist the users in
shooting is a difficult point. Computational aesthetics is the use of machines
to achieve a computable human aesthetic process (Neumann et al., 2005), so
that computational methods can make applicable aesthetic decisions in a
similar way to humans. Therefore, the system with the combination of robots
and computational aesthetics makes automatic time-lapse photography possible.

This paper proposes an automatic time-lapse photography framework, which
conducts the virtual and real-world time-lapse photography with robots by
calculating the time-lapse photography parameters of the three-dimensional
scene. We propose a parameter optimization method to obtain the optimal
shooting parameters. The innovation is to parameterize the complex process of
video generation of time-lapse photography. We classify the parameters into
image aesthetics-related, video aesthetics-related, and time-lapse aesthetics-
related according to the characteristics of time-lapse photography, and
propose parameter optimization methods for these three types of parameters.
According to the characteristics of different parameters, we design and employ
different aesthetic models for parameter optimization.

Meanwhile, we propose a time-lapse photography interface for users to view or
adjust the parameters generated automatically by our method and use virtual
robots to take virtual photography in three-dimensional scenes. Users can
preview the time-lapse photography videos, and adjust the parameters according
to their own preferences. The system also supports the export of parameters
and provides them to real robots to enable robots to move and control the
camera in the real world, shooting time-lapse videos with specific time
parameters.In addition, we propose a time-lapse photography aesthetic
evaluation method that can automatically evaluate the aesthetic quality of
time-lapse video. We apply it to validate our parameter optimization method.

The main contributions of this work can be summarized as follows:

  * •

a time-lapse photography framework that can automatically and efficiently
shoot time-lapse videos.

  * •

a method for time-lapse photography parameter optimization using different
aesthetic quality assessment models.

  * •

a time-lapse video assessment method to evaluate the aesthetic quality of
time-lapse video.

  * •

a time-lapse photography interface for users to conduct the virtual or real-
world time-lapse photography.

##  2\. Related Work

Existing work on time-lapse photography has focused more on generating
optimised time-lapse videos from existing images or videos, i.e. the post-
processing of the data. There is no research on how to go about capturing
time-lapse video or the acquisition of the data.

Figure 1. Illustration of the proposed framework. The 3D scene is
reconstructed with the images. The shooting parameters are then determined by
analyzing the 3D scene in the simulation environment. The interactive
interface enables users to modify or view parameters to complete virtual
photography. The system can also provide parameters for real robots to
complete real-world time-lapse photography.

##  3\. Method

To obtain the satisfactory time-lapse videos, careful preparations are needed
for photographers in advance, including scenery selection, shooting plan, etc.
Moreover, time-lapse photography requires photographers to take a large number
of images over a long period of time, which requires a great deal of patience
and carefulness. Therefore, in order to generate the optimal shooting plan so
that the robot can automatically generate the time-lapse videos or assist the
user in photography, our goal is to design a time-lapse shooting system that
can find a series of optimal parameters resulting in the high quality of the
final video. So we define the overall quality function $Q$ as follows:

(1) $Q=Q_{I}(\mathcal{P_{I}})+Q_{V}(\mathcal{P_{V}})+Q_{T}(\mathcal{P_{T}}),$

where $Q_{I}$, $Q_{V}$, $Q_{T}$ are the quality functions for image, video and
time-lapse video respectively. $\mathcal{P_{I}}$, $\mathcal{P_{V}}$,
$\mathcal{P_{T}}$ are the shooting parameters.

Secondly, different users may have different preferences when shooting time-
lapse videos. The system should therefore support parameter personalization to
suit different users. Moreover, in case of automatic robot shooting, it also
involves a series of parameters of the robot.

In order to solve the above difficulties, we propose an automatic time-lapse
photography framework, as shown in Figure 1.

3D Reconstruction. When shooting a scene, experienced photographers will
usually investigate the scene in advance, which takes a lot of time to
determine the location and plan the shot. The transformation of light and
shadow is a very important factor in time-lapse photography, so the shooting
time is also important. In order to capture the desired light and shadow, some
photographers need to rehearse the plan in advance to ensure that they can
efficiently conduct the photography at a specific time during the actual
shooting. Therefore, adequate preparation is very important for photography,
especially time-lapse photography.

3D reconstruction allows us to quickly understand an unknown scene and react
accordingly. The analysis of a 3D model in a simulation environment is an
effective way to fleetly determine the parameters for time-lapse photography
based on the shooting scene.

First we reconstruct the shooting scene. For large-scale scenes, image-based
3D reconstruction is a low-cost and efficient reconstruction method. The
images can be collected on the Internet or captured by drone or robot. We then
use structure-from-motion, an image-based 3D reconstruction method, to obtain
a textured and realistic 3D scene. In order to obtain reachable shooting
points, we limit the shooting height and generate reachable areas.

The image-based reconstructed 3D scene is static, but dynamic objects are
important in time-lapse photography. So we add dynamic elements to the static
3D scene, mainly including people, traffic, weather and lighting. Accessible
regions on the ground are chosen to simulate the flow of people and vehicles
to ensure that the dynamic elements are similar to those in the real world.

Shooting Parameter. For time-lapse photography we first determine the relevant
parameters, including shooting location $L$, shooting angle $A$, shooting path
$P$, time range $[S_{t},E_{t}]$, frame extraction interval $\Delta t$. As
shown in Equation 1, the overall video quality is mainly related to the image
quality, ordinary video quality, and time-lapse video quality. Therefore, we
classify the relevant parameters into three categories as well.

The first category is the viewfinder parameters $\mathcal{P_{I}}=\\{L,A\\}$,
which are related to the image quality, including the shooting location $L$
and shooting angle $A$. The second category is the camera movement parameters
$\mathcal{P_{V}}=\\{P\\}$ that are related to the ordinary video quality since
the camera motion of time-lapse photography is very similar to the camera
motion of ordinary photography. The last category is the time-warp parameters
$\mathcal{P_{T}}=\\{[S_{t},E_{t}],\Delta t\\}$, including time range
$[S_{t},E_{t}]$ and time interval $\Delta t$, which are the characteristic
parameters of time-lapse video to reflect the time change.

After obtaining the 3D model we can analyze the geometry and texture of the 3D
model, as well as the image or video rendered under different lighting
conditions, in order to determine the optimal set of shooting parameters.

Video Post. After determining the relevant parameters, we can determine the
initial position of the camera according to the shooting position and
orientation, then determine the camera movement according to the shooting
path, and finally make the camera take a set of photos according to the time
range and interval to get the final time-lapse video.

The flicker might occur when the final time-lapse videos are created (Martin-
Brualla et al., 2015). Essentially instead of the camera’s settings between
shots remaining perfectly constant or purposefully changing in a slight and
gradual way, large unintended exposure jumps occur in a few of the frames
which create images that look out of place when compiled with other images.
One of the ways that reduce the likelihood of flicker to occur is shooting in
manual mode (Chylinski, 2012). However, since the light is changing in our
scenes and manual mode might quickly over or under expose the images as the
scene changes, we still use auto exposure mode, which makes it necessary to
post-process the captured images to reduce flicker. We smooth out the exposure
in the time-lapse video with histogram equalization (Pizer et al., 1987) to
create a good visual effect. Finally we generate the time-lapse video with the
processed image sequences.

It can be seen that the results using the parameters selected by the system
are better than those selected by the amateurs. In the first two scenarios,
the time-lapse videos generated by our method are of higher quality than those
captured manually. The video quality captured by the professional
photographers in the latter two scenes is the highest. The overall results
show that our approach is better than the amateurs and comparable to the
professional photographers. Therefore, it can be considered that our system
can help amateurs to get more professional time-lapse videos. And our system
can automatically generate high-quality time-lapse videos, which also saves
photographers a lot of time.

##  4\. Results

In this section we present two applications based on our system, namely
virtual and real-world shooting.

###  4.1. Virtual Shooting

The system supports the import of 3D scenes, which can be obtained either
automatically and quickly through 3D reconstruction or through manual modeling
(which is time-consuming but more detailed). In the imported 3D scene, users
can take virtual videos, which can be applied to film and television
production, scene roaming, etc. The virtual photography can be classified into
automatic photography and user-defined photography.

####  4.1.1. Automatic Shooting

Users can just import the appropriate 3D scenes and the system can
automatically obtain the optimal parameters and generate video results for
preview with our method.

####  4.1.2. Customization Parameters

In addition, the system supports user-defined parameters and modification of
parameters so that the photography process can be done manually or semi-
automatically. As shown in the Figure 2, users can click and drag to complete
the virtual scene selection, and then choose the appropriate mode of camera
movement from the four shooting modes. After that, users can set the time
parameter to preview the lighting change instantly without waiting for a long
time. The instant feedback from the interactive system makes it easier for
users to set parameters.

###  4.2. Real-World Shooting

The system also supports outputting parameters and providing them to the
robotics platform to shoot time-lapse videos in real environments.

As shown in the Figure 3, the parameters of the viewfinder are the camera
position and orientation, which are converted into GPS coordinates and the
pitch and yaw angles of the pan-tilt. The parameter of the camera motion is
the motion path, which is converted into the motion trajectory of the robot.
The time parameters are the time range and time interval, which we convert
into the start and end time of the robot and the shooting interval. After the
robot receives the relevant parameters, it will shoot on the spot according to
the shooting plan.

###  4.3. Evaluation

####  4.3.1. Self-Evaluation.

The parameters associated with time-lapse photography are classified into
three types of parameters: image, video and time-lapse aesthetics-related. In
order to verify the validity of each type of parameter, we conduct the self-
evaluation experiment.

We choose five simulated environments, each divided into four areas to be
filmed, so that a total of 20 time-lapse videos are taken. The resulting
videos are scored using our time-lapse video aesthetic quality assessment
method. We contrast our method of parameter selection with random parameter
selection. To get the overall results for random parameters, we randomly
select 10 sets of parameters and average them. We first test the results of
random selection of all parameters, and then use our method to select
parameters in the order of image, video, and time.

Figure 2. Virtual shooting. The first line is the visual interface, and the
second line is the interface for parameter adjustment. From left to right:
camera position, camera orientation, camera movement and time parameters.
Figure 3. Real-World shooting. The parameters generated by the system need to
be converted into the parameters of the robot.

Table 1 shows the experimental results of the self-evaluation. It can be seen
that our system progressively reaches the optimum in determining the
parameters. For every type of parameter, our optimization method outperforms
the random selection method.

####  4.3.2. Comparison with Photographers.

In order to verify the effectiveness of the overall system, we invited two
professional photographers with rich experience and two novices who lack
photography experience to manually set the parameters in the simulation scene
and take photography. We selected five different shooting areas. The user
needs to select the appropriate position and orientation in the area to be
shot, and determine the way of camera movement and the shooting time. Finally,
the system is used for virtual photography to generate a time-lapse video. At
the same time, the system also automatically takes virtual photography in the
scene with our optimized parameters. The time-lapse photography aesthetic
quality evaluation model is utilized to score the final time-lapse photography
video, and the results are shown in the Table 2.

It can be seen that the results using the parameters selected by the system
are better than those selected by the amateurs. The overall results show that
our approach is better than the amateurs and comparable to the professional
photographers. Therefore, it can be considered that our system can help
amateurs to get more professional time-lapse videos. And our system can
automatically generate high-quality time-lapse videos, which also saves
photographers a lot of time.

Table 1. Alation study. The check mark indicates that our method was used for parameter optimization. The unmarked term means random selection of parameters. | Image | Video | Time | Quality
---|---|---|---|---
1 |  |  |  | 0.508
2 | ✓ |  |  | 0.522
3 | ✓ | ✓ |  | 0.574
4 | ✓ | ✓ | ✓ | 0.635
Table 2. Comparison results with the shots of photographers. Scene | Amateur | Professional | Ours
---|---|---|---
1 | 0.420 | 0.577 | 0.582
2 | 0.485 | 0.567 | 0.596
3 | 0.486 | 0.580 | 0.537
4 | 0.479 | 0.547 | 0.505
Avg. | 0.468 | 0.568 | 0.555

##  5\. Summary

In this paper, we propose an aesthetics driven framework for automatic time-
lapse photography, which replaces photographers with virtual and real robots
to complete the time-consuming work of time-lapse photography. First, we
parameterized the complex time-lapse photography process and proposed a method
to optimize the parameters with the aesthetic models. Then we proposed a
specific time-lapse video aesthetic quality evaluation model to automatically
evaluate the aesthetic quality of time-lapse photography. In addition, we also
proposed an automatic time-lapse photography system that supports users to
manually, semi-automatically or automatically complete virtual photography.
The system can also import parameters to real robots to complete real
photography. Experimental results and user studies show that our method can
efficiently complete time-lapse photography automatically, and the effect is
equivalent to that of professional photographers, which proves that our system
can help inexperienced users to shoot high-quality time-lapse videos.

## References

  * (1)
  * Bennett and McMillan (2007) Eric P Bennett and Leonard McMillan. 2007.  Computational time-lapse video.  In _ACM SIGGRAPH 2007 papers_. 102–es.
  * Chylinski (2012) Ryan Chylinski. 2012\.  _Time-lapse photography: A Complete Introduction to Shooting, Processing, and Rendering Time-lapse Movies with a DSLR Camera_.  LearnTimelapse. com.
  * Lou et al. (2021) Hao Lou, Heng Huang, Chaoen Xiao, and Xin Jin. 2021\.  Aesthetic Evaluation and Guidance for Mobile Photography. In _Proceedings of the 29th ACM International Conference on Multimedia_. 2780–2782.
  * Martin-Brualla et al. (2015) Ricardo Martin-Brualla, David Gallup, and Steven M Seitz. 2015\.  Time-lapse mining from internet photos.  _ACM Transactions on Graphics (TOG)_ 34, 4 (2015), 1–8.
  * Neumann et al. (2005) L Neumann, M Sbert, B Gooch, W Purgathofer, et al. 2005\.  Defining computational aesthetics.  _Computational aesthetics in graphics, visualization and imaging_ (2005), 13–18.
  * Pizer et al. (1987) Stephen M Pizer, E Philip Amburn, John D Austin, Robert Cromartie, Ari Geselowitz, Trey Greer, Bart ter Haar Romeny, John B Zimmerman, and Karel Zuiderveld. 1987.  Adaptive histogram equalization and its variations.  _Computer vision, graphics, and image processing_ 39, 3 (1987), 355–368.
  * Rawat and Kankanhalli (2016) Yogesh Singh Rawat and Mohan S Kankanhalli. 2016.  ClickSmart: A context-aware viewpoint recommendation system for mobile photography.  _IEEE Transactions on Circuits and Systems for Video Technology_ 27, 1 (2016), 149–158.
  * Rawat et al. (2019) Yogesh Singh Rawat, Mubarak Shah, and Mohan S Kankanhalli. 2019\.  Photography and Exploration of Tourist Locations Based on Optimal Foraging Theory.  _IEEE Transactions on Circuits and Systems for Video Technology_ 30, 7 (2019), 2276–2287.