remarkRemark hypothesisHypothesis claimClaim Coarse-graining multi-agent
stochastic systemsD. Stepanova, H. M. Byrne, P. K. Maini and T. Alarcón
[supp-]SupplementaryMaterial

# A method to coarse-grain multi-agent stochastic systems with regions of
multistability††thanks: . This work is supported by a grant of the Obra Social
La Caixa Foundation on Collaborative Mathematics awarded to the Centre de
Recerca Matemàtica through a scholarship awarded to D.S. D.S. and T.A. have
been partially funded by the CERCA Programme of the Generalitat de Catalunya.
They also acknowledge MINECO (https://www.ciencia.gob.es/) for funding under
grants MTM2015-71509-C2-1-R and RTI2018-098322-B-I00. D.S. and T.A.
participate in project 2017SGR01735 which was awarded by AGAUR
(https://agaur.gencat.cat/en/inici/index.html) but with no actual funding. The
funders had no role in study design, data collection and analysis, decision to
publish, or preparation of the manuscript. H.M.B. and P.K.M. received no
specific funding for this work.

Daria Stepanova222Centre de Recerca Matemàtica, Bellaterra (Barcelona) 08193,
Spain 333Departament de Matemàtiques, Universitat Autònoma de Barcelona,
Bellaterra (Barcelona) 08193, Spain 666 Helen M. Byrne444Wolfson Centre for
Mathematical Biology, Mathematical Institute, University of Oxford, Oxford OX2
6GG, UK Philip K. Maini444Wolfson Centre for Mathematical Biology,
Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK Tomás
Alarcón555Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona
08010, Spain 222Centre de Recerca Matemàtica, Bellaterra (Barcelona) 08193,
Spain 333Departament de Matemàtiques, Universitat Autònoma de Barcelona,
Bellaterra (Barcelona) 08193, Spain<EMAIL_ADDRESS>
###### Abstract

Hybrid multiscale modelling has emerged as a useful framework for modelling
complex biological phenomena. However, when accounting for stochasticity in
the internal dynamics of agents, these models frequently become
computationally expensive. Traditional techniques to reduce the computational
intensity of such models can lead to a reduction in the richness of the
dynamics observed, compared to the original system. Here we use large
deviation theory to decrease the computational cost of a spatially-extended
multi-agent stochastic system with a region of multi-stability by coarse-
graining it to a continuous time Markov chain on the state space of stable
steady states of the original system. Our technique preserves the original
description of the stable steady states of the system and accounts for noise-
induced transitions between them. We apply the method to a bistable system
modelling phenotype specification of cells driven by a lateral inhibition
mechanism. For this system, we demonstrate how the method may be used to
explore different pattern configurations and unveil robust patterns emerging
on longer timescales. We then compare the full stochastic, coarse-grained and
mean-field descriptions via pattern quantification metrics and in terms of the
numerical cost of each method. Our results show that the coarse-grained system
exhibits the lowest computational cost while preserving the rich dynamics of
the stochastic system. The method has the potential to reduce the
computational complexity of hybrid multiscale models, making them more
tractable for analysis, simulation and hypothesis testing.

###### keywords:

Large deviation theory, coarse-graining, phenotype pattern formation,
multiscale modelling, hybrid modelling

60F10 (Large deviations), 92C15 (Developmental biology, pattern formation),
92C42 (Systems biology, networks), 92B05 (General biology and biomathematics),
92-08 (Computational methods for problems pertaining to biology)

##  1 Introduction

When modelling a biological process, one has to make choices on how detailed
the model should be in order to capture the characteristic features of the
system. At the same time, the model should be as simple as possible in order
to facilitate its analysis and numerical simulations. The evolution of systems
with large numbers of agents (e.g. molecules, cells, species) can be described
by the average behaviour of their agents, or their mean-field limits using
(ordinary or partial) differential equations ([6, 9, 34]). Dynamical systems
theory provides methods and techniques for the analysis and numerical
simulations of such systems. This description might become insufficient when
the system comprises agents with internal variables that change in time, thus
altering the agents’ behaviour, or when the system is not ‘large enough’ to be
described accurately by the mean-field equations. For these systems,
stochastic descriptions are employed [38] (for example, continuous time Markov
chains, CTMCs, or stochastic differential equations, SDEs). In biological
systems, the number of agents is finite and some level of noise is always
present which can affect the system dynamics [38]. While exhibiting richer
dynamics than deterministic systems, stochastic models are more
computationally intensive.

Furthermore, in order to formulate a theoretical model of a biological
phenomenon, it is often necessary to account for dynamics that act on
different temporal and/or spatial scales [2, 17]. This has led to the
development of hybrid multiscale models, in which different modelling
techniques may be applied at each scale and then efficient coupling algorithms
are used to integrate these models (see, e.g., [8, 16, 36] and references
therein). In many of these models, individual entities (cells, species, etc.)
are considered as discrete agents which are, themselves, equipped with models
for their internal states determining the behaviour (e.g. subcellular
signalling, cell cycle, response to extracellular stimuli). Such models have
great potential for generating insights into the behaviour of a system (e.g.,
endothelial cell rearrangements [3], cell differentiation and tissue
organization in intestinal crypts [8], and multiscale cancer modelling [10]).
However, they frequently become numerically intractable because of their
complexity (e.g. the internal dynamics of agents) [2]. This limits possible
applications of these models.

(a)

(b)

(c)

Figure 1: Cell phenotype specification. (a) Phenotype (Delta-high and Delta-
low cells) patterning of cells induced by a mechanism of lateral inhibition in
two different domains: a cell monolayer and a branching network. (b) Dynamic
time evolution of phenotype adaptation of an individual cell. Using a
phenotype proxy, e.g. level of Delta, allows for identification of a
continuous cell phenotype. (c) Phenotype switches, as in (b) (dashed vertical
lines), occur due to either a change in a cell’s microenvironment or naturally
present noise in intracellular signalling.

In this work, we explain how to reduce the computational complexity of a
hybrid model by coarse-graining the internal dynamics of its agents when these
are described by a stochastic system with multiple steady states. The method
involves applying large deviation theory (LDT) to reduce the dynamics of the
stochastic system to a continuous time Markov chain (CTMC) on the state space
of its stable steady states. LDT provides a theoretical framework with which
to quantify how small time-dependent fluctuations can lead to significant
deviations from the mean-field behaviour (rare events) such as transitions
between stable steady states which cannot occur in deterministic systems [23].
This approach has previously been used to study rare, noise-induced events in
individual stochastic systems [14, 15, 19, 38, 39, 41], but to our knowledge,
this is its first application to a multi-agent model.

In previous work, we developed a multiscale model of angiogenesis [44], the
process of growth of new blood vessels from pre-existing ones [28], which
accounts for gene expression patterns (phenotypes) of endothelial cells (ECs)
at the subcellular scale. For prescribed levels of extracellular stimuli, the
system is either monostable (i.e. only one cell phenotype exists) or bistable
(i.e. two stable steady states, cell phenotypes, coexist). Cell phenotype is
specified via contact-dependent cross-talk with neighbouring ECs via the VEGF-
Delta-Notch signalling pathway [5, 24]. VEGF, or vascular endothelial growth
factor, is the activating external stimulus; Delta and Notch are transmembrane
ligands and receptors, respectively, which can trans-bind, (i.e. a ligand on
one cell can bind to a receptor on another cell, thus allowing the two cells
to ‘communicate’). Cells adjust their gene expression in order to maintain a
pattern of two distinct phenotypes, Delta-high and Delta-low cells (see
LABEL:PhenotypeSwitch_Motivation_Config and
LABEL:PhenotypeSwitch_Motivation_Trajectory). We use the internal level of
Delta as a proxy to distinguish between the phenotypes. In angiogenesis, the
Delta-high (Delta-low) cells are referred to as tip (stalk) cells [5]. The
number of transmembrane proteins in this signalling pathway is on the order of
thousands for each cell [6]. Therefore, in order to formulate a mathematical
model, it is tempting to use deterministic mean-field equations to describe
the kinetic reactions of this signalling pathway. However, deterministic
descriptions cannot account for noise-induced transitions between stable
steady states or, in the case of this signalling pathway, phenotypic switches,
which can occur in regions of bistability (see
LABEL:PhenotypeSwitch_Motivation_Trajectory and
LABEL:PhenotypeSwitch_Motivation_Switch). Since branching patterns of vascular
networks are affected by the distribution of cells with different phenotypes,
such phenotype transitions are potentially significant. Therefore, we modelled
the subcellular signalling pathways stochastically, which increased the
computational cost of the model. This example illustrates a general problem
associated with computational and, in particular, hybrid models: in order to
preserve emergent features of the system, such as continuous cell phenotypes
and noise-induced phenotype switches, the model becomes computationally
intractable for large lattice simulations.

We illustrate the coarse-graining method by reference to the subcellular model
of the VEGF-Delta-Notch signalling pathway that defines cell phenotype. The
core Delta-Notch signalling pathway plays a key role in phenotype adaptation
in cell types which can form cell monolayers, such as epithelial sheets [35,
43], bristle patterning in Drosophila [11, 30, 13], and neural precursor cells
[22]. In all of these biological processes, the lateral inhibition mechanism
of Delta-Notch signalling generates spatial patterns of cells with alternating
fates (phenotypes). In the VEGF-Delta-Notch model, the stationary distribution
of VEGF serves as an activating extracellular stimulus for the particular case
of endothelial cells. In other cell types, which use the lateral inhibition
mechanism to communicate, the external stimulus may differ from VEGF. In this
paper we perform our simulations for two spatial geometries: a cell monolayer
and a branching network (LABEL:PhenotypeSwitch_Motivation_Config). For our
model of multicellular VEGF-Delta-Notch signalling, we show typical simulation
results of the coarse-grained system which allows us to explore different
configurations of spatial patterns in a single realisation of the model (due
to phenotypic switches). We then demonstrate how this dynamic exploration of
possible patterns may be used to uncover robust patterns emerging at long
timescales. We finally compare the spatio-temporal dynamics and computational
cost of the full stochastic CTMC, the coarse-grained and the deterministic
mean-field descriptions. Our results show that the coarse-grained model, while
preserving the continuous description of cell phenotype and rare events of
phenotype switching, is more computationally efficient than the other two
systems. Thus, it significantly reduces the computational complexity of the
model without sacrificing the rich dynamics of the original stochastic system.

The remainder of the paper is organised as follows. In Section 2, we review
the hybrid (multiscale) modelling approach (Section 2.1) and summarise large
deviation theory (Section 2.2). This provides us with the information needed
to formulate the coarse-grained model in Section 3. In Section 3.1, we start
by coarse-graining the individual agent system and checking the accuracy of
the method. We then extend the technique to a multi-agent system in Section
3.2 where we outline a general algorithm for formulating and simulating the
coarse-grained model. In Section 4, we present typical simulation results for
the model of the VEGF-Delta-Notch signalling pathway (Section 4.2) and compare
the full stochastic, coarse-grained and mean-field systems via metrics which
quantify the spatial patterns formed by the two cell phenotypes and we also
compare computational cost of simulations (Section 4.3). The paper concludes
in Section 5 with a summary of our findings and suggestions for future
research directions.

##  2 Theoretical background

###  2.1 Hybrid models

Biological systems are often highly complex, involving processes that may
interact across multiple spatial and temporal scales (see Figure 2). From a
general perspective, the subcellular scale is characterised by intracellular
chemistry (e.g. gene expression, signal transduction and receptor/ligand
dynamics). Subcellular processes determine behaviour at the cellular scale and
may generate emergent properties at the tissue scale. In addition to this
upward coupling across spatial scales, there is downward coupling whereby
extracellular chemicals and biomechanical cues influence the subcellular
chemistry/mechanics within a cell. In this way, dynamic interactions,
encompassing all the scales, can occur (Figure 2).

(a)  Figure 2: A schematic diagram illustrating characteristic spatial and
temporal scales of a typical biological process and coupling between them. The
VEGF-Delta-Notch signalling pathway, which serves as an illustrative example
for application of the CG method, acts at the subcellular scale (highlighted
in blue) on a timescale shorter than other processes (e.g. cell migration,
cell-extracellular matrix interaction at the tissue scale) involved in the
multiscale model of angiogenesis [44]. As a result, we may use LDT theory to
coarse-grain its dynamics.

From the theoretical perspective, models which consider only processes at a
single spatial/temporal scale do not allow for investigation of emergent
features which manifest at other scales (for example, collective migration or
phenotype patterning which arise from individual cell dynamics and govern
tissue scale organisation). Equally, difficulties associated with the physical
interpretation of parameters in phenomenological models, i.e. large scale
models which capture the overall evolution of a biological process, make it
challenging to fit the model to biological data. In particular, this abstract
parameter construct hinders model calibration/validation and limits potential
applications of the models. Multiscale models, which couple processes at
different spatial and/or temporal scales, have the potential to address these
issues [4].

A challenge in formulating a multiscale model relates to the number of
entities (protein, cells, extracellular components, etc.) that should be
included at each scale of interest. Using the same mathematical formalism to
model processes involving entities which vary in number by several orders of
magnitude may lead to the omission of essential features or make the model
computationally intractable. Hybrid approaches are increasingly being
recognised as suitable tools for trying to overcome problems of this type and
have become a key part of multiscale modelling [16, 17]. The central idea is
to employ the modelling framework most suitable to each subprocess and then to
couple them. For example, the extracellular environment and signalling cues
are usually modelled deterministically due to the large number of proteins
involved. On the other hand, cells may be treated as individual entities,
equipped with a subcellular model which determines their behaviour (e.g.
proliferation, cell polarity and migration). This framework has been used to
develop multiscale models of cancer (see reviews [16, 40] and references
therein), angiogenesis [28], collective cell migration [17], among other
examples [2].

Hybrid modelling allows for efficient parameter estimation and model
visualisation, forging interdisciplinary collaboration between researchers in
theoretical modelling and experimental biology [2, 36]. There is also the
potential of using high-throughput experimental data to develop more detailed
multiscale models. As an example, one of the aspects of biological systems
that has received little attention in theoretical modelling is the effect of
stochasticity in the response of individual entities to external stimuli [17].
Hybrid modelling allows investigation of this effect on the collective,
emergent behaviour. However, increasing computational complexity makes these
models intractable for large-scale simulations [16].

This challenge motivated us to develop a technique which reduces the
computational complexity of a model while preserving its stochasticity. The
method is applicable to systems characterised by stochastic processes which
exhibit multistability and which evolve on timescales shorter than those
associated with other system processes. The example that we study in this
paper is of this type: the subcellular dynamics of cell fate determination via
lateral inhibition (a bistable, stochastic system) act on a shorter timescale
than those associated with, for example, cell migration, and tissue scale
processes such as the dynamics of extracellular soluble factors (e.g.
diffusion, secretion by cells, degradation) [28] (Figure 2). This observation
motivates us to use large deviation theory to coarse-grain the dynamics
associated with intracellular signalling to produce a jump process (i.e. a
Markov chain) on the stable state space of the steady states of the original
system which describes the VEGF-Delta-Notch pathway.

###  2.2 Large deviation theory (LDT)

In the presence of noise, small fluctuations can drive significant deviations
from mean-field behaviour such as, for example, transitions from one stable
steady state to another. These transitions are usually referred to as rare
events since their likelihood is small. LDT is predicated on the assumption
that when rare events occur, the system follows the least unlikely paths.
Deviations from these paths occur with very small probability (i.e. smaller
than the probability of a rare event). Specifically, Freidlin-Wentzell’s
theory of large deviations predicts that the deviations are exponentially
suppressed [23], making such transitions ‘predictable’. LDT provides the means
to analyse the frequency of rare events and to identify the maximum likelihood
path (minimum action path, MAP) along which these transitions can occur.

A stochastic differential equation (SDE) of a diffusion process,
$x^{\epsilon}\in\mathbb{R}^{n}$, has the following form

(1)
$\mathrm{d}x^{\epsilon}(t)=b(x^{\epsilon})\mathrm{d}t+\sqrt[]{\epsilon}\sigma(x^{\epsilon})\mathrm{d}W,$

where $b:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ is a drift vector,
$a(x^{\epsilon})=(\sigma\sigma^{T})(x^{\epsilon})$ is a diffusion tensor
($\sigma:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}\times\mathbb{R}^{m}$, $m$
corresponds to the number of kinetic reactions in the system), $W$ is a Wiener
process in $\mathbb{R}^{m}$ and $\epsilon=\Omega^{-1}$ is noise amplitude.

The mean-field limit of Equation 1, $x(t)\in\mathbb{R}^{n}$, solves the
following differential equation:

(2) $\derivative{x}{t}=b(x).$

Assume that Equation 2 has two stable steady states,
$x_{1},~{}x_{2}\in\mathbb{R}^{n}$, whose basins of attraction form a complete
partition of $\mathbb{R}^{n}$. We are interested in transitions from
$x_{1}\rightarrow x_{2}$ (and $x_{2}\rightarrow x_{1}$) which cannot be
accounted for unless noise is present in the system.

A key player in LDT is the action functional

$S_{T}(\psi)=\begin{cases}\displaystyle\int_{0}^{T}L(\psi,\dot{\psi})\mathop{dt},&\text{if
$\psi\in C(0,T)$ is absolutely continuous and}\\\\[-7.0pt] &\text{the integral
converges,}\\\\[5.0pt] +\infty,&\text{otherwise,}\end{cases}$

which is computed for a transition path $\psi:[0,T]\rightarrow\mathbb{R}^{n}$
from $x_{1}$ to $x_{2}$ ($\psi(0)=x_{1}$ and $\psi(T)=x_{2}$, $T$ is the
transition time). Here
$L(x,y)=\displaystyle\sup_{\theta\in\mathbb{R}^{n}}\left(\langle
y,\theta\rangle-H(x,\theta)\right)$ is the large deviation Lagrangian, with
$\langle\cdot,\cdot\rangle$ being the Euclidean scalar product in
$\mathbb{R}^{n}$ and $H(x,\theta)$ being the Hamiltonian associated with
$L(x,y)$. The particular form of the Hamiltonian depends on the dynamical
system under consideration (in LABEL:supp-appendix:gMAM, we explain how to
define the Hamiltonian for an SDE such as Equation 1 and a general birth-death
CTMC).

The action functional is used to estimate the probability that a trajectory
$x^{\epsilon}(t)$ lies in a narrow neighbourhood, of width $\delta>0$, of a
given path $\psi\in C(0,T)$ (see Figure 3 for an illustration):

(3) $\operatorname{P}\left\\{\sup_{0\leq t\leq T}\mid
x^{\epsilon}(t)-\psi(t)\mid<\delta~{}\middle|~{}x^{\epsilon}(0)=x_{1}\right\\}\approx\exp{-\epsilon^{-1}S_{T}(\psi)}.$

(a)  Figure 3: An illustration of a transition path between two stable steady
states of an arbitrary bistable system. The two stable steady states, $x_{1}$
and $x_{2}$, are marked by filled red circles; an unstable saddle point is
marked by an unfilled red circle. The transition path, $\psi(t)$, from $x_{1}$
to $x_{2}$ is shown by a thick green line, whereas a single stochastic
trajectory, $x^{\epsilon}(t)$, is indicated by a thin black path. The shaded
blue region indicates a $\delta$-neighbourhood around $\psi(t)$ ($\delta$ as
defined in Equation 3).

Since the probability function in Equation 3 decreases as the action
functional, $S_{T}(\psi)$, increases, the maximum likelihood path, $\psi^{*}$,
is the minimiser of $S_{T}(\cdot)$. This leads naturally to the idea of the
quasipotential:

(4)
$V(x_{1},x_{2})=\displaystyle\inf_{T>0}~{}\inf_{\psi\in\overline{C}_{x_{1}}^{x_{2}}(0,T)}S_{T}(\psi).$

Here $\overline{C}_{x_{1}}^{x_{2}}(0,T)$ is the space of absolutely continuous
functions $f:[0,T]\rightarrow\mathbb{R}^{n}$ such that $f(0)=x_{1}$ and
$f(T)=x_{2}$. Roughly speaking, the quasipotential gives an estimate of how
‘difficult’ it is to move from $x_{1}$ to $x_{2}$.

On timescales which are much longer than those associated with relaxation to a
stable steady state, the dynamics of Equation 1 can be reduced, or coarse-
grained, to that of a CTMC on the state space of the two stable steady states,
$\left\\{x_{1},x_{2}\right\\}$, with transition rates

(5) $k_{x_{1}\rightarrow
x_{2}}\asymp\exp\left(-\epsilon^{-1}V(x_{1},x_{2})\right),\qquad
k_{x_{2}\rightarrow
x_{1}}\asymp\exp\left(-\epsilon^{-1}V(x_{2},x_{1})\right).$

Here $\asymp$ denotes log-asymptotic equivalence so that $f(\epsilon)\asymp
g(\epsilon)$ if and only if
${\lim_{\epsilon\rightarrow 0}\frac{\log f(\epsilon)}{\log g(\epsilon)}=1}$.

In practice, most double minimisation problems, such as Equation 4, do not
have a solution for finite $T>0$. Furthermore, closed-form Lagrangians exist
for SDEs of the type defined by Equation 1 but not for general birth-death
CTMCs. Equation 4 can be reformulated in terms of a Hamiltonian system of the
form

$\derivative{\phi}{t}=\partialderivative{H(\phi,\theta)}{\theta},\qquad\derivative{\theta}{t}=-\partialderivative{H(\phi,\theta)}{\phi}.$

This problem must be solved as a boundary-value problem, i.e. $\phi(0)=x_{1}$
and $\phi(T)=x_{2}$, on an infinite time interval, $T\rightarrow\infty$, [26]
which makes it a non-trivial numerical problem. Thus the traditional LDT
methods are inapplicable in most cases.

One way to resolve these problems is to reformulate the minimisation problem
defined by Equation 4 on the space of curves (i.e. transition paths from one
stable steady state to another). In [29], Heymann and Vanden-Eijnden proved
that the minimisation problem defined by Equation 4, is equivalent to

(6)
$V(x_{1},x_{2})=\displaystyle\inf_{\phi}\widehat{S}(\phi),~{}~{}\text{with}~{}~{}\widehat{S}(\phi)=\displaystyle\sup_{\begin{subarray}{c}\hat{\theta}:[0,1]\rightarrow\mathbb{R}^{n}\\\
H(\phi,\hat{\theta})=0\end{subarray}}\displaystyle\int_{0}^{1}\langle\phi^{\prime},\hat{\theta}\rangle\mathop{d\alpha},$

where $\phi:[0,1]\rightarrow\mathbb{R}^{n}$ is a curve from $x_{1}$ to $x_{2}$
parametrised by standard arc length.

The geometric reformulation, Equation 6, resolves analytically the issue of
the infinite time, $T$, in the original minimisation problem. Furthermore,
only the Hamiltonian is needed. In this respect, the method is more general as
it can be applied to SDEs, CTMCs and other systems for which the Hamiltonian
is known (see LABEL:supp-appendix:gMAM in Supplementary Material).

In [29], an algorithm was developed to efficiently compute $V(x_{1},x_{2})$
and the corresponding minimiser, $\phi^{*}$, from the geometric reformulation.
The algorithm is known as the geometric minimum action method (gMAM) and the
minimiser, $\phi^{*}$, of the action functional is referred to as the minimum
action path (MAP) (for more details see LABEL:supp-appendix:gMAM).

Once the quasipotential has been computed, the coarse-grained system is given
by a CTMC, with rates defined by Equation 5.

##  3 Coarse-graining (CG)

We now illustrate how the theory described in the previous section can be used
to coarse-grain a specific hybrid multiscale model, one for which the internal
dynamics of the agents are described by multistable stochastic systems. This
property is characteristic of, for example, systems driving cell fate
(phenotype) determination. We begin by using LDT to formulate a CG model for a
system comprising a single agent (here a cell). The subcellular signalling
pathway, which we use to illustrate the method, is the VEGF-Delta-Notch
pathway (see LABEL:supp-appendix:VEGFDeltaNotch in Supplementary Material and
[44] for details). This pathway regulates phenotypic adaptation via lateral
inhibition [12, 35]. This system meets the requirements for application of the
CG technique: (a) it is bistable; its stable steady states are associated with
cellular phenotypes (Delta-high and Delta-low cells); (b) we are interested in
its evolution on timescales longer than the typical time for relaxation to an
equilibrium since other processes (e.g. cell migration and dynamics of
extracellular matrix) act on longer timescales (see Figure 2).

We then extend the method to the general case of multi-agent systems. Here the
dynamics of each entity is coarse-grained to a CTMC on the state space of its
stable states, and coupling between the internal dynamics of individual agents
is achieved via the external variables whose dynamics depend on the states of
neighbouring agents and/or the time evolution of these variables. We outline
below how we apply this method to a monolayer of cells (motivated by phenotype
patterning via the core Delta-Notch pathway in cell monolayers [35]) and a
branching network (angiogenesis-motivated application [44]) that interact via
VEGF-Delta-Notch signalling.

###  3.1 Individual agent system

(a)  Figure 4: A flowchart of the procedure used to coarse-grain a multistable
stochastic system for an individual entity. The steady state solutions,
quasipotential and prefactor depend on the model parameters and external
variables, $v\in\mathrm{R}^{V}$ ($V$ indicates the dimension of the vector of
external variables). Here the transition rates, $k_{x_{s}\rightarrow x_{l}}$,
are defined by Equation 7, the prefactor, $C_{x_{s}\rightarrow x_{l}}$, is
determined from Equation 8b, and $\overline{\Omega}$ is given by Equation 9.

Our algorithm for coarse-graining a stochastic system with a region of
multistability involving a single entity is illustrated in Figure 4. For the
particular case of VEGF-Delta-Notch signalling, a cell’s internal state
(phenotype) depends on two model parameters (inputs) corresponding to the
extracellular levels of Delta and Notch,
$v=\left(d_{ext},n_{ext}\right)\in\mathrm{R}^{2}$ (see LABEL:supp-
appendix:VEGFDeltaNotch). We fix the values of the model parameters and the
external variables, $v$ (see LABEL:supp-Params). We then use the mean-field
system defined by LABEL:supp-eq:DN_single_nondimensional to compute the steady
state solutions. For this example, the values of the external variables, $v$,
are chosen so that the system is bistable; the two stable steady states
correspond to Delta-high and Delta-low cell phenotypes,
$\left\\{x_{1},x_{2}\right\\}=\left\\{\text{Delta-high},~{}\text{Delta-
low}\right\\}$, and the unstable steady state is an unstable saddle. Our goal
is to compute the transition rates of the CG system which we approximate as
follows:

(7) $k_{x_{s}\rightarrow x_{l}}\approx C_{x_{s}\rightarrow
x_{l}}\exp\left(-\Omega V(x_{s},x_{l})\right),\quad
s,~{}l\in\left\\{1,2\right\\},~{}s\neq l.$

We note that the prefactor, $C_{x_{s}\rightarrow x_{l}}$, arises from the
asymptotic equivalence relation defined by Equation 5. The system size is
given by $\Omega=\epsilon^{-1}$, where $\epsilon$ is the noise level.

We use the gMAM to compute the quasipotential values and corresponding paths
(MAPs) for transitions between the Delta-high and Delta-low phenotypes (for
more details, see LABEL:supp-appendix:MAP in Supplementary Material). An
illustrative example is shown in Figure 5, where we compare the MAPs and
sample paths of the full stochastic CTMC for an individual cell (see also
LABEL:supp-FullStochasticSystem in LABEL:supp-appendix:VEGFDeltaNotch).
Several characteristic features of the phenotype transitions are noteworthy.
First, the dynamics of the MAP can be split into two parts: the transition
from the steady state of origin to the saddle point (for example, from the
Delta-low phenotype to the saddle point, indicated by the blue circle in
LABEL:MAP_Path_s2t) which is possible due to the presence of noise. The main
contribution to the quasipotential comes from this transition. The MAP from
the unstable saddle point to the stable steady state of destination (from the
saddle point indicated by the blue circle to the Delta-high phenotype in
LABEL:MAP_Path_s2t) follows the fastest route given by the deterministic
heteroclinic orbit connecting the steady states (i.e. the unstable saddle and
the stable Delta-high cell state). The second noteworthy feature of the
phenotype transitions is that, as the level of noise, $\epsilon$, decreases,
the stochastic sample path follows the MAP more closely (compare
LABEL:MAP_Path_s2t and LABEL:MAP_Path_t2s for which $\Omega=\epsilon^{-1}=70$
and $\Omega=\epsilon^{-1}=450$, respectively).

(a)

(b)

Figure 5: An illustration of the minimum action paths (MAPs) and stochastic
sample paths for transitions between the Delta-high and Delta-low cell
phenotypes. We computed the MAPs (indicated by the dotted magenta lines) for
the subcellular VEGF-Delta-Notch system in an individual cell using the gMAM
for transitions from (a) Delta-low to Delta-high cell and (b) Delta-high to
Delta-low cell. The stochastic sample paths obtained by simulating the full
stochastic CTMC model (LABEL:supp-FullStochasticSystem) with the system sizes
(a) $\Omega=70$, (b) $\Omega=450$, are plotted in black. The thin grey lines
indicate streamlines of the corresponding mean-field system (LABEL:supp-
eq:DN_single_nondimensional). The Delta-high (Delta-low) cell stable steady
state is indicated by a green (red) filled circle; the unstable saddle by a
blue unfilled circle. The plots represent three-dimensional projections of the
full five-dimensional system as defined by LABEL:supp-
eq:DN_single_nondimensional. Parameter values are fixed as indicated in
LABEL:supp-Params.

To fully determine the CG transition rates, the prefactor value,
$C_{x_{s}\rightarrow x_{l}}$, must be estimated. From Equation 7, for
$s,~{}l\in\left\\{1,2\right\\},~{}s\neq l$, we have

(8a) $\displaystyle\log\langle T^{\Omega}_{x_{s}\rightarrow x_{l}}\rangle$
$\displaystyle\approx\Omega V(x_{s},x_{l})-\log C_{x_{s}\rightarrow
x_{l}}~{},$ (8b) $\displaystyle\log C_{x_{s}\rightarrow x_{l}}$
$\displaystyle\approx\Omega V(x_{s},x_{l})-\log\langle
T^{\Omega}_{x_{s}\rightarrow x_{l}}\rangle~{},$

where $\langle T^{\Omega}_{x_{s}\rightarrow
x_{l}}\rangle=1/k_{x_{s}\rightarrow x_{l}}$ is the mean passage time between
the stable steady states, $x_{s}$ and $x_{l}$ (Delta-high and Delta-low
phenotypes), for a fixed value of the system size, $\Omega$. $\langle
T^{\Omega}_{x_{s}\rightarrow x_{l}}\rangle$ can be determined from direct
simulation of the full stochastic model using the reaction kinetics given in
LABEL:supp-FullStochasticSystem.

(a)

(b)

Figure 6: Convergence of the quasipotential, $V(x_{s},x_{l})$, as the system
size, $\Omega$, increases. We ran 1000 realisations of the stochastic VEGF-
Delta-Notch model for an individual cell (see LABEL:supp-FullStochasticSystem)
for fixed values of $d_{ext}=0.2$, $n_{ext}=0.5$ and increasing system size,
$\Omega$. We plotted the convergence to the quasipotential value (a)
$V(\text{Delta-low},\text{Delta-high})$ and (b) $V(\text{Delta-
high},\text{Delta-low})$ as a function of $\Omega$ (black circle markers). For
these parameter values, transitions from the Delta-low to Delta-high phenotype
are less likely to occur (higher noise levels, $\epsilon=\Omega^{-1}$, and/or
longer transition times are needed) than transitions from the Delta-high to
Delta-low phenotype (see Equation 8a). Therefore, the perturbations of this
random event are smaller and convergence is reached for higher values of
noise. This is why lower values of $\Omega$ in (a) suffice to accurately
determine the prefactor value from Equation 8. The blue dashed lines indicate
the value of the corresponding quasipotential computed via the gMAM; the red
dotted lines indicate $\overline{\Omega}$ from Equation 9. All other parameter
values are fixed as indicated in LABEL:supp-Params.

An accurate estimate of the quasipotential (as obtained via the gMAM) allows
us to obtain the prefactor given the mean passage time, $\langle
T^{\Omega}_{x_{s}~{}\rightarrow~{}x_{l}}\rangle$, for a single value of the
system size, $\Omega$. However, the approximate relation in Equation 8 is
valid in the limit $\Omega\rightarrow\infty$ (see Figure 6). Thus, $\Omega$
should be chosen sufficiently large to achieve convergence in Equation 8 and,
at the same time, not too large in order to ensure that transitions between
the phenotypes occur in a computationally feasible time, since the waiting
times for transitions between stable steady states increase exponentially as
$\Omega$ grows. Specifically, we fix a maximum simulation time, $T_{max}$, and
an average prefactor value, $\bar{C}$ (determined computationally from
simulations), and approximate the corresponding system size,
$\overline{\Omega}$, as:

(9) $\overline{\Omega}\approx\frac{\log T_{max}+\log\bar{C}}{V(x_{s},x_{l})}.$

Then the prefactor, $C_{x_{s}\rightarrow x_{l}}$, can be approximated using
Equation 8b with $\Omega=\overline{\Omega}$.

From Equation 8a, we know that $\log\langle
T^{\Omega}_{x_{s}~{}\rightarrow~{}x_{l}}\rangle$ is a linear function of
$\Omega$ whose slope and intercept are given by the quasipotential,
$V(x_{s},x_{l})$, and $\left(-\log C_{x_{s}~{}\rightarrow~{}x_{l}}\right)$,
respectively. Thus, in order to check the accuracy of our estimate for the
system size, $\overline{\Omega}$ (Equation 9), we compared linear fitting of
data obtained from the full stochastic CTMC model for increasing $\Omega$,
with the estimate obtained from the gMAM quasipotential and the prefactor
extracted from simulations with system size, $\overline{\Omega}$. The results
presented in Figure 7 show that the estimates converge as $\Omega$ increases,
confirming the accuracy of the two methods.

(a)

(b)

Figure 7: Prefactor estimation. Comparison of prefactor estimates obtained
from simulations of the full stochastic CTMC model (black circles) and
estimates obtained using the gMAM-quasipotential and mean passage times for a
single value of the system size, $\overline{\Omega}$ (blue line), see Equation
8a. The linear fit of the full stochastic data (red line) was performed for
values of $\Omega$ such that the corresponding sample
$\left\\{T^{\Omega}_{x_{s}~{}\rightarrow~{}x_{l}}\right\\}$ is exponentially
distributed (high levels of noise might affect the distribution of these
transitions). Panel (a) corresponds to the transition from Delta-low to Delta-
high phenotype; panel (b) corresponds to the transition from Delta-high to
Delta-low phenotype. The red dotted lines indicate $\overline{\Omega}$ from
Equation 9. All other parameter values are fixed as indicated in LABEL:supp-
Params.

To summarise, we coarse-grain the stochastic VEGF-Delta-Notch dynamics as
follows (see Figure 4):

  1. I

Fix the model parameter values and the vector of external variables, $v$,
which, for this system, is given by the extracellular levels of Delta and
Notch, $v=\left(d_{ext},n_{ext}\right)$.

  2. II

Compute the steady states of the corresponding mean-field system (LABEL:supp-
eq:DN_single_nondimensional).

  3. III

Formulate the CG model:

    1. i

If, for the given $v=\left(d_{ext},n_{ext}\right)$, the system is monostable
(either Delta-high or Delta-low cell steady state exists), then the
quasipotential value to arrive at this state is 0. The value of the other
quasipotential can be assumed infinite (since the system is monostable, this
transition is impossible). For example, if the only stable steady state is the
Delta-high cell, then $V(\text{Delta-low},\text{Delta-high})=0$ and
$V(\text{Delta-high},\text{Delta-low})=\infty$. The CG model is defined by its
unique stable steady state.

    2. ii

If the system is within the bistable regime (both Delta-high and Delta-low
steady states are stable), then the CG model is defined as a CTMC on the state
space of $\\{x_{s},x_{l}\\}=\\{\text{Delta-high},$ $\text{Delta-low}\\}$. The
transition rates are given by Equation 7. The quasipotential,
$V(x_{s},x_{l})$, is approximated using the gMAM; the prefactor value,
$C_{x_{s}~{}\rightarrow~{}x_{l}}$, is obtained via Equation 8b from stochastic
simulations of the full VEGF-Delta-Notch model for a fixed value of the system
size, $\overline{\Omega}$, defined by Equation 9.

  4. IV

The CG model can be simulated using any variant of the SSA, such as, for
example, the classical Gillespie algorithm [25].

The above method generalises naturally for systems with an arbitrary number of
stable steady states (see Figure 4). In this case, the quasipotential and the
corresponding prefactor must be approximated for each pair of stable steady
states. The method can also be applied to systems which possess other
attractors, e.g. limit cycles [15, 23].

###  3.2 Multi-agent system

(a)  Figure 8: A flowchart of the procedure to coarse-grain a multi-agent
stochastic system with a region of multistability. A pseudocode of the
simulation algorithm for the multi-agent CG model is presented in LABEL:supp-
appendix:SimulationAlgorithm. The simulation part of the diagram illustrates
an iteration of the Gillespie algorithm for simulation of multi-agent CG
systems. Here $T_{final}$ stands for the final simulation time;
$\mathrm{Exp}(\lambda)$ is an exponential distribution of intensity,
$\lambda$.

In this section we show how the CG method can be applied to multi-agent
systems with a region of multistability. In this case, the dynamics of each
agent is coarse-grained to that of a CTMC between its stable steady states for
given values of the external variables, $v$, which establish the coupling
between the internal dynamics of individual agents ($v$ depends on the state
of agents in the local environment of the focal agent and/or time, and defines
its internal state, e.g. phenotype). If the dynamics of an individual agent
are independent of its neighbours and time (i.e. the values of the external
variables are constant) then we use the CG method described in Section 3.1
(see also Figure 4). A suitable range of values for the external variables,
$v\in\mathcal{V}$, where $\mathcal{V}\subset\mathrm{R}^{V}$, can be determined
by simulating the original multiscale model. Here $V$ indicates the dimension
of the vector of external variables, $v$. In order to reduce the computational
cost in the multi-agent CG system, it is convenient to calculate a priori
look-up tables for the steady states, quasipotential and prefactor values for
a discretisation, $\left\\{v_{j}\right\\}_{j\in\mathcal{J}}\subset\mathcal{V}$
(here, $j$ indexes entries in the generated discretisation; $\mathcal{J}$ is
the size of the discretisation). Interpolation routines can then be used to
establish an input-output relationship between an arbitrary $v\in\mathcal{V}$
and the values of the corresponding steady states and the transition rates
between them. Therefore, we split the general CG method for multi-agent
systems into two steps (see Figure 8):

  1. (i)

Pre-simulation: calculate look-up tables for the system steady states, quasi-
potential and prefactor values for each entry in a discretisation,
$\left\\{v_{j}\right\\}_{j\in\mathcal{J}}$, for a range of values of the
external variables, $\mathcal{V}\subset\mathrm{R}^{V}$.

  2. (ii)

Simulation: the CG model is simulated (via, e.g., the Gillespie algorithm) as
a CTMC on a state space defined by the steady states of all of its entities,
with the coupling maintained via the external variables, $v$, updated at each
simulation step according to entities’ local environments and/or time.

We now provide more details on the pre-simulation and simulation steps.

####  3.2.1 Pre-simulation: look-up tables

Pre-computed look-up tables of system steady states, quasipotential and
prefactor values are used to interpolate the values of the system steady
states and the CG transition rates between them for an arbitrary set of values
of the external variables, $v$, without calculating them explicitly at each
step during simulations of the CG model. In a general setting, the dimension
of each table is equal to $V$, the dimension of the vector of external
variables.

The steady states must be computed numerically for each entry $v_{j}$ in the
discretisation, $\left\\{v_{j}\right\\}_{j\in\mathcal{J}}$, using the mean-
field limit for an individual entity (as described in Section 3.1). For values
of $v_{j}$ that fall within the multistability region, the quasipotential is
computed via the gMAM in a pair-wise manner, for each pair of stable steady
states, $\left\\{x_{s}\right\\}_{s=1}^{\mathcal{S}}$. The last look-up table
corresponds to the prefactor, $C_{x_{s}\rightarrow x_{l}}$, $x_{s},~{}x_{l}\in
1\ldots\mathcal{S}$, which must be approximated for each $v_{j}$ within the
multistability region. The prefactor values are obtained from Equation 8b as
before, using the mean passage times, $\langle
T^{\Omega}_{x_{s}~{}\rightarrow~{}x_{l}}\rangle$, which are determined by
simulating the full stochastic model with the system size,
$\Omega=\overline{\Omega}$, defined by Equation 9.

####  3.2.2 Simulation algorithm

Once all the look-up tables have been computed, the multi-agent CG system can
be simulated as a standard Gillespie algorithm (or one of its variants, e.g.,
Next Subvolume method [20]) in which the total propensity, $P$, at each time
step is computed as a sum of transitions, $k^{e}_{x_{s}\rightarrow x_{l}}$,
for each entity, $e$, to switch its (stable) state (see Figure 8). The steady
states corresponding to each entity (and the transition rates between them)
for the exact value of the external variables, $v^{e}\in\mathcal{V}$, ($v^{e}$
has to be computed for each entity, $e$, according to its microenvironment)
are interpolated via appropriate numerical routines. We present pseudocode for
the simulation procedure in LABEL:supp-appendix:SimulationAlgorithm.

Note that our CG method does not account for the initial, relatively short
(compared to the LDT timescale), relaxation time during which the system
relaxes onto the timescale on which the CG approximation is valid. Thus, it is
necessary to obtain an initial stable steady state configuration, i.e to pre-
pattern the system, using either the full stochastic CTMC or the mean-field
model (see Figure 8 and line 5 in LABEL:supp-SimulationAlgorithm). The final
simulation time for the pre-patterning should be large enough to ensure that
the system relaxes to an equilibrium. Since this procedure is performed only
once, it does not affect the computational complexity of the CG simulations.
We have chosen to use the mean-field system to pre-pattern our simulations
since it is less time-consuming and the stochasticity (i.e. transitions
between phenotypes) is preserved later in the CG simulation loop.

##  4 Results

For illustrative purposes, we consider the specific example of spatial
phenotype patterning via the Delta-Notch lateral inhibition mechanism in
response to an external signalling cue (VEGF). First, we provide more details
about our implementation of the CG model and present typical simulation
results and the robust patterns that emerge at long times. We then discuss the
relative merits of the CG method, using a variety of metrics to compare its
performance with the original stochastic and mean-field systems. We used the
Next Subvolume method [20] for simulations of the full stochastic CTMC and the
Euler-Lagrange method (explicit scheme) for the numerical integration of the
mean-field equations.

###  4.1 CG model of spatial cell phenotype patterning

The multicellular VEGF-Delta-Notch (i.e. the Delta-Notch signalling pathway
coupled with external VEGF stimulation) model is bistable (see LABEL:supp-
appendix:VEGFDeltaNotch in Supplementary Material). When simulated in a two-
dimensional geometry, it produces ‘salt-and-pepper’ patterns in which the
phenotypes of neighbouring cells alternate between Delta-high and Delta-low
states [44]. For this model, cross-talk between individual cells is achieved
via external variables, $d_{ext}$ and $n_{ext}$, which represent the levels of
Delta and Notch, respectively, summed over cells in a circular neighbourhood
with a fixed interaction radius, $R_{s}$ (see Figure 9 and LABEL:supp-
appendix:VEGFDeltaNotch). Hence, for this system,
$v=\left(d_{ext},n_{ext}\right)$ defines a cell’s internal state (phenotype)
and the dimension of the pre-computed look-up tables is 2 (see Section 3.2.1).
We determined a suitable range,
$\mathcal{V}=\left[0,d_{ext}^{max}\right]\times\left[0,n_{ext}^{max}\right]\subset\mathrm{R}^{2}$,
for these variables by running 100 realisations of the multiscale model of
angiogenesis (the number of realisations depends on the model of interest).

(a)  Figure 9: A schematic diagram showing the non-local interactions in the
multicellular VEGF-Delta-Notch model. Cell-to-cell interactions may be non-
local (i.e. beyond immediate neighbours on a given lattice) provided they lie
within an interaction radius, $R_{s}$. The diagram illustrates the weights of
interactions between the focal cell (highlighted in blue) and cells in its
neighbourhood, for a regular hexagonal lattice (the weights are defined as a
normalised area of the overlap between a neighbouring voxel and the circular
neighbourhood of the focal cell, see LABEL:supp-ExtDN in LABEL:supp-
appendix:VEGFDeltaNotch).

We then generated a regular discretisation of $\mathcal{V}$,
$\left\\{v_{j}\right\\}_{j\in\mathcal{J}}$, with a grid $100\times 100$. For
each $v_{j}$ in this grid, we computed the steady states for the mean-field
limit defined by LABEL:supp-eq:DN_single_nondimensional using non-linear
solvers from the C++ GNU Scientific Library (GSL). We note that, once the
steady states of the full system have been computed, the subcellular variables
$\iota$, $r_{2}$ and $r_{2}^{*}$, corresponding to the Notch intracellular
domain, VEGF receptor 2 (VEGFR2) and VEGF-VEGFR2 complexes, respectively, (see
definitions in LABEL:supp-appendix:VEGFDeltaNotch in Supplementary Material)
are redundant; it is not necessary to track these variables because the input-
output relationship between $v=\left(d_{ext},n_{ext}\right)$, and the steady
states completely defines the configuration of the system.

(a)

(b)

(c)

(d)

Figure 10: An illustration of the quasipotential surfaces. Upper panels: a
noise-induced transition from Delta-high (in magenta) to Delta-low (in black)
phenotype of a single cell during a simulation of the angiogenesis model [44]
plotted as a function of the focal cell’s (a) Delta and (b) Notch levels. The
external Delta, $d_{ext}$, (Notch, $n_{ext}$) for the focal cell is computed
as a weighted sum of the Delta (Notch) levels of its neighbours as defined by
LABEL:supp-MulticellularMeanField. Lower panels: 2D projections of the
quasipotential surfaces (c) $V(\text{Delta-low},\text{Delta-high})$ and (d)
$V(\text{Delta-high},\text{Delta-low})$ as functions of $d_{ext}$ and
$n_{ext}$. The monostability region in which the unique stable steady state
corresponds to a Delta-high (Delta-low) cell is coloured green (red). The
colour bar indicates the value of the corresponding quasipotential. The
trajectory (as in panels (a) and (b)) plotted on the quasipotential surfaces
(in (c) and (d)), illustrates that phenotype switches are more likely to occur
for lower values of the quasipotential. Parameter values are fixed as
indicated in LABEL:supp-Params.

For values of $v_{j}$ that fall within the bistability region, we computed the
quasipotential values of the transitions between phenotypes (see Figure 10),
using the gMAM (see LABEL:supp-appendix:gMAM in Supplementary Material). We
also used the full stochastic system to check those values of the
quasipotential for which a phenotype switch is more likey to occur. As
expected, most phenotype transitions occur close to the boundary of the
bistability region, where values of the quasipotential are lower. For example,
LABEL:QasipotentialIllustration_trajectory_D and
LABEL:QasipotentialIllustration_trajectory_N show a sample path of the full
stochastic system for an individual cell during a simulation of the multi-
agent model [44]. The cell undergoes a noise-induced switch from a Delta-high
to a Delta-low phenotype. LABEL:QasipotentialIllustration_s2t and
LABEL:QasipotentialIllustration_t2s show the same sample path projected onto
the quasipotential surfaces. These plots show that phenotypic switches are
more likely to occur when the values of external Delta and Notch,
$\left(d_{ext},n_{ext}\right)$, are such that the quasipotential,
$V(x_{1},x_{2})=V(\text{Delta-high},\text{Delta-low})$, is small.

We constructed a look-up table of prefactor values, $C_{x_{s}\rightarrow
x_{l}}$, $x_{s},~{}x_{l}\in\\{\text{Delta-high},$ $\text{Delta-low}\\}$, by
approximating the mean passage times, $\langle
T^{\overline{\Omega}}_{x_{s}\rightarrow x_{l}}\rangle$, (sample size of 1000
realisations) for an individual cell to switch its phenotype from simulations
of the full stochastic CTMC (LABEL:supp-FullStochasticSystem) with the system
size, $\overline{\Omega}$, given by Equation 9.

We then implemented the CG model in C++ using LABEL:supp-SimulationAlgorithm.
In order to establish an input-output relationship between an arbitrary
$v=\left(d_{ext},n_{ext}\right)$ and the corresponding cell phenotypes and
transition rates, we used bilinear interpolation routines from the C++ GNU
Scientific Library (GSL) (gsl_interp2d routines). The model was then simulated
using the standard Gillespie algorithm. We used no-flux boundary conditions to
compute for each cell the extracellular levels of Delta and Notch in all our
simulations.

###  4.2 Spatial patterning in the CG model

In order to illustrate the CG model, we first ran numerical simulations on a
small cell monolayer ($10\times 12$ voxels). The results presented in
LABEL:PatternConfigurations_S1, LABEL:PatternConfigurations_S2,
LABEL:PatternConfigurations_S3, and LABEL:PatternConfigurations_S4 show how
the distribution of Delta-high and Delta-low cells changes over time during a
typical CG realisation (see also Movie S1). Starting from an initial pre-
pattern (LABEL:PatternConfigurations_S1), noise-induced phenotype transitions
enable the system to explore different pattern configurations for the given
geometry, while the proportion of Delta-high cells remains on average constant
(see LABEL:PatternConfigurations_TipProportion).

(a)

(b)

(c)

(d)

(e)

Figure 11: Different pattern configurations explored by the CG model. (a)-(d)
Series of plots showing how the distribution of cell phenotypes changes over
time during a single simulation of the CG model. The colour bar indicates the
level of Delta. (a) $t=0$; (b) $t=40$; (c) $t=260$; (d) $t=410$ minutes. (e)
Time evolution of the Delta-high cell proportion (defined as a ratio of cells
with the Delta-high phenotype to the total cell number) for a single
simulation of the CG model (blue line) and averaged over $1000$ realisations
(red line). For these simulations, the interaction radius and system size were
fixed at $R_{s}=15\mu m$ and $\Omega=100$, respectively; the values of the
remaining parameters were fixed as indicated in LABEL:supp-Params.

The mean proportion of Delta-high cells (and, thus, the spatial pattern)
during simulations of the CG system depends on the interaction radius,
$R_{s}$. For values of $R_{s}$ corresponding to nearest-neighbours interaction
($R_{s}\leq 1.5h$, where $h$ is the voxel width), we observe classical
patterns of alternating Delta-high and Delta-low cells (i.e. the so-called
salt-and-pepper pattern [12]; see LABEL:supp-PatternVaryingR10). As $R_{s}$
increases, the number of Delta-low cells that may be inhibited by a focal
Delta-high cell increases, causing the proportion of Delta-high cells in the
spatial patterns to decrease [44]. Thus, for larger values of $R_{s}$
($R_{s}>1.5h$), Delta-high cells are separated by larger distances (see
LABEL:supp-PatternVaryingR20, LABEL:supp-PatternVaryingR30, and LABEL:supp-
PatternVaryingR40). These results for CG simulations are consistent with those
obtained for the full multicellular stochastic model of the VEGF-Delta-Notch
signalling pathway [44]. The ability of the CG system to explore different
spatial patterns increases as the size of the interaction radius, $R_{s}$,
grows, and the corresponding emerging patterns are more diverse (see
LABEL:PatternConfigurations_S1, LABEL:PatternConfigurations_S2,
LABEL:PatternConfigurations_S3, LABEL:PatternConfigurations_S4, LABEL:supp-
PatternVaryingR20, LABEL:supp-PatternVaryingR30, and LABEL:supp-
PatternVaryingR40).

It is noteworthy that spatial patterns explored in simulations of the CG model
differ in their robustness to noise. In particular, the mean passage time for
a phenotype switch, and, thus, a change in the pattern, to occur, which is
equal to the inverse of the total propensity, $P$, depends on the values of
the quasipotential, $V(x_{s},x_{l})$, for all entities in the system. Here,
the total propensity, $P$, for a phenotype switch event is defined as a sum of
transition rates, $k^{e}_{x_{s}\rightarrow x_{l}}$, for each cell with index,
$e$, to change its state from $x_{s}$ to $x_{l}$, see Figure 8. When, via
random exploration, the system finds a configuration for which the values of
$V(x_{s},x_{l})$ are larger, the waiting time for a phenotype switch increases
and the configuration is more resilient to further changes.

This feature of the CG method facilitates exploration of new robust spatial
patterns which cannot practically be achieved using other numerical
frameworks: (i) simulations of the full stochastic model are too
computationally intensive, which makes the exploration of these patterns
infeasible because of the longer timescales needed; (ii) the deterministic
framework does not allow for transitions between stable steady states, which
makes this exploration impossible; (iii) the complexity of analytic methods
needed to verify the stability of a pattern of a system with non-local
interactions does not permit exploration of complex pattern configurations
[37].

(a)

(b)

(c)

(d)

Figure 12: Emergence of robust pattern configurations in simulations of the CG
model. At long times, via exploration of different pattern configurations, the
dynamics of the CG system evolve to a robust pattern in which any further
phenotype switches are unlikely. (a) A typical emergent pattern for a single
realisation of the CG model (the colour bar indicates the level of Delta, $d$,
for each cell). (b) The time evolution of the total propensity, $P$, for a
phenotype switch to occur. Cells in the border rim (three-cell width) are
excluded from $P$ since, due to the model geometry, they do not possess a
‘robust’ configuration of neighbours. As $P$ decreases to $0$, the waiting
time for a phenotype switch to occur approaches infinity, and the pattern
becomes more robust to change. (c)-(d) The dynamics of an individual cell
(outlined in cyan in (a)) during this simulation. (c) Temporal evolution of
the internal level of Delta, $d$, (defining cell phenotype: high (low) values
of $d$ correspond to Delta-high (Delta-low) phenotype) and that in its
microenvironment, $d_{ext}$. (d) Temporal evolution of transition rates for a
phenotype switch for this cell. We note that the large difference in the order
of values for transition rates for the total propensity, $P$, of the lattice
($O(10^{6})$), plot (b), and for an individual cell
($O(10^{-17})-O(10^{-3})$), plot (d), comes from the contribution to $P$ of
transition rates for cells which are, for the given values of the external
variables, on the border of the bistability region (see Figure 10). For these
simulations, the interaction radius and system size were fixed at $R_{s}=15\mu
m$ and $\Omega=1000$, respectively; the values of all remaining parameters
were fixed as indicated in LABEL:supp-Params.

We now present simulation results which illustrate the ability of the CG
method to uncover new spatial patterns for the VEGF-Delta-Notch system at long
times. We fixed the interaction radius at $R_{s}=3.0h=15\mu m$ ($h=5\mu m$ is
the voxel width), so that interactions occur between cells that are first and
second order neighbours in the lattice; the noise amplitude was fixed at
$\epsilon=\Omega^{-1}=0.001$. We ran a CG simulation on a medium size
monolayer of cells (see LABEL:OptimalPattern_CellHighlighted and Movie S2).
Starting from the initial pre-pattern, the CG model explores various patterns
until it eventually settles on a more robust configuration (shown in
LABEL:OptimalPattern_CellHighlighted). In order to confirm our prediction
regarding pattern robustness, we plotted the temporal evolution of the total
propensity of the lattice, $P$, in LABEL:OptimalPattern_TotalPropensity. As
its value decreases, $P\rightarrow 0$, the mean waiting time for a change in
the spatial pattern becomes infinite, which accounts for the robustness of the
emerging pattern. We also considered the dynamics of an individual cell (its
position in the monolayer is highlighted by a cyan line in
LABEL:OptimalPattern_CellHighlighted). LABEL:OptimalPattern_DeltaLevel shows
how the phenotype of this cell changes over time: at early times, the cell
switches between Delta-high and Delta-low phenotypes (low (high) values of
subcellular Delta, $d$, correspond to Delta-low (Delta-high) phenotype). As
the spatial pattern settles to a robust configuration, the cell’s environment,
i.e. the levels of Delta of its neighbours, $d_{ext}$, stop changing and the
cell acquires a Delta-high phenotype that remains unchanged for the rest of
the simulation. The transition rates for phenotype switches for this cell
(LABEL:OptimalPattern_TransitionRate) exhibit similar dynamics to the total
propensity, $P$, of the whole lattice (LABEL:OptimalPattern_TotalPropensity).

Our CG simulation results show that this robust pattern configuration is not
unique. However, we note that the spatial patterns tend to have a regular
structure; for example, Delta-high cells may be organised in similar clusters
comprising two or three cells as in the pattern shown in
LABEL:OptimalPattern_CellHighlighted. These configurations have lower values
of the total propensity, $P$. Cells on the border of the lattice undergo
phenotype switches (see Movie S2), since they cannot attain this ‘more robust’
combination of neighbours for the given geometry (since we use no-flux
boundary conditions in our simulations).

###  4.3 Comparison of the full stochastic, coarse-grained and mean-field
frameworks

We compared the dynamics of the multicellular VEGF-Delta-Notch model using
three frameworks: (i) full stochastic CTMC, (ii) CG, and (iii) mean-field
descriptions. Simulated (using any of these frameworks) on a 2D domain, the
model produces a characteristic pattern of ECs with two cell phenotypes (see,
for example, Figures 11 and LABEL:supp-PatternVaryingR). Since the CG
approximation describes the long-term behaviour of the system, when its
evolution is dominated by the timescale associated with phenotypic switches,
it does not account for the initial relaxation onto a quasi-steady state
pattern. Thus, the three frameworks cannot be compared with respect to their
behaviour at early evolution times. Instead, we quantified the final pattern
and the computational cost of simulations. The final simulation time,
$t=T_{final}$, was chosen sufficiently large to ensure that a steady state
pattern had been established for the mean-field simulations (since stochastic
systems do not have a steady state pattern in a classical sense). In order to
systematically compare the three frameworks, we used the same final simulation
time, $t=T_{final}$, for the other two systems.

We used the following set of metrics to compare the dynamics of the three
mathematical descriptions (Supplementary Material):

  * •

Delta-high cell proportion, which is defined as the ratio of the number of
cells with Delta-high phenotype to the total number of cells in the system;

  * •

distribution of Delta-high cell clusters, which provides a breakdown of sizes
of Delta-high cell clusters (adjacent cells with Delta-high phenotype, e.g., a
single Delta-high cell, two adjacent Delta-high cells, etc.) in a steady
pattern configuration;

  * •

computational cost, which is defined as the average CPU time (in seconds) to
perform a single realisation of model simulation.

Since the pre-calculated look-up tables for the CG simulations (Section 3.2.1)
were computed for a fixed set of model parameters (see LABEL:supp-Params), we
held them fixed for all simulations. However, the cell-to-cell interaction
radius, $R_{s}$, which is used in the multicellular simulations to determine
for each cell, $e$, the vector of extracellular variables,
$v^{e}=\left(d^{e}_{ext},n^{e}_{ext}\right)$, may vary. In our simulations, we
used $R_{s}\in\left\\{5,~{}7.5,~{}10,~{}12.5,~{}15\right\\}~{}\mu m$ which
correspond to experimental observations of the distance over which cell-to-
cell interaction can occur in endothelial cells [18] (which corresponds to up
to three cells in the interaction circle). Nonetheless, from a theoretical
point of view, this quantity can take any value greater than the half-width of
a voxel, $R_{s}>0.5h$, where $h$ is the voxel width (we fix $h=5\mu m$ in our
simulations). In addition, for the full stochastic CTMC and CG descriptions,
we vary the noise amplitude, $\epsilon=1/\Omega$, by changing the system size
parameter, $\Omega$. We used
$\Omega\in\left\\{50,~{}100,~{}200,~{}500,~{}1000\right\\}$. The larger the
value of $\Omega$, the closer will be the dynamics of a stochastic system to
its mean-field description. For each numerical setup ($R_{s}$ and $\Omega$),
we ran 100 realisations.

We considered two simulation geometries: a 2D cell monolayer and a branching
network.

##### Setup 1: a cell monolayer

We first ran numerical simulations on a cell monolayer (see LABEL:supp-
InitialSetup_M). This spatial geometry was motivated by the biological process
of cell fate specification induced by lateral inhibition via Delta-Notch
signalling in flat domains. Examples of such cell fate specification include
bristle patterning in Drosophila notum [11, 30, 13], and differentiation of
neural precursors in neurogenesis [22] (see [7, 35] and references therein for
other examples). The fixed stationary distribution of the VEGF serves as an
external stimulus which enhances lateral inhibition via Delta-Notch
signalling. We chose VEGF as an illustrative example, although, depending on
the specific system, other extracellular signals will provide cell stimulus.

(a)

(b)

Figure 13: Comparison of the dynamics of the multicellular VEGF-Delta-Notch
model simulated on a cell monolayer using the full stochastic (CTMC), CG, and
mean-field descriptions. (a) The Delta-high cell proportion as a function of
the cell-to-cell interaction radius, $R_{s}$, for varying noise amplitude,
$\epsilon=1/\Omega$ (the value of $\Omega$ is indicated in the title of each
plot), for the full stochastic CTMC (black), CG (blue) and mean-field (red)
descriptions. To explore different possible patterns in the deterministic
mean-field system, we created a small initial perturbation to the initial
configuration (LABEL:supp-InitialSetup_M). (b) A series of barplots showing
how the long-time distribution of Delta-high cell clusters changes as the
interacton radius, $R_{s}$, varies for the full stochastic CTMC (left panel),
CG (middle panel), and mean-field (right panel) systems. The number of single
Delta-high cells in the final pattern (i.e. at a fixed final simulation time)
is shown in blue; the number of clusters with 2, 3, and 4 adjacent Delta-high
cells is shown in yellow, green, and red, respectively. For these simulations,
we fixed $\Omega=1000$ ($\epsilon=0.001$). The results are averaged over 100
realisations. The remaining parameter values were fixed as indicated in
LABEL:supp-Params.

We began by considering the dynamics of the Delta-high cell proportion for
this spatial geometry (see LABEL:ResultsMonolayer_TipProportion). Consistent
with the previous results [44], for all simulation frameworks (i.e. the full
stochastic (CTMC), CG, and mean-field descriptions), the Delta-high cell
proportion decreases as the cell interaction radius, $R_{s}$, increases.
LABEL:ResultsMonolayer_TipProportion confirms that, as expected, differences
in this metric between the three systems decrease as the level of noise is
reduced (i.e. as $\Omega$ increases). In particular, for high noise levels
(i.e. lower values of $\Omega$), the patterns generated by the stochastic
systems (full CTMC and CG frameworks) are more diverse, and the Delta-high
cell proportions differ from those for the associated mean-field description.
We note that the dynamics of the Delta-high cell proportion for the mean-field
system (red lines) are identical in all subplots in
LABEL:ResultsMonolayer_TipProportion since noise is absent in deterministic
systems (i.e. the system size parameter, $\Omega$, is irrelevant).

We also quantified the size distribution of the Delta-high cell clusters
associated with the final patterns established on the cell monolayers. Since
the dynamics of the three systems converge for larger values of the system
size, $\Omega$ (as shown in LABEL:ResultsMonolayer_TipProportion),
LABEL:ResultsMonolayer_Clusters shows results for this metric computed for
simulations with $\Omega=1000$. The distributions are in good quantitative
agreement for the three systems. The discrepancy for simulations with larger
cell interaction radius (e.g. $R_{s}=15\mu m$) arises because (for this value
of $\Omega$) the CG system is more likely to explore long timescale patterns
which have a more ‘regular’ structure and are more robust to noise (cells with
Delta-high phenotype organised in similar clusters, see Section 4.2).

##### Setup 2: a branching network

We next considered a more complex spatial geometry of a small branching
network (see LABEL:supp-InitialSetup_N) extracted from a simulation of a
hybrid model of angiogenesis [44]. LABEL:supp-SnapshotsCGVascularNetwork shows
a series of patterns explored by the CG system at different time points during
a typical simulation for this configuration (for the full simulation, see
Movie S3).

For this spatial configuration, we compared the three simulation frameworks
using the same metrics as for the cell monolayer. The results for the Delta-
high cell proportion are presented in LABEL:supp-
ResultsVascularNetwork_TipProportion. We find that the number of possible
patterns generated by lateral inhibition is lower for the branching network
geometry than for the cell monolayer (see LABEL:supp-
SnapshotsCGVascularNetwork). Consequently, the Delta-high cell proportions
converge for smaller values of $\Omega$ (compare
LABEL:ResultsMonolayer_TipProportion and LABEL:supp-
ResultsVascularNetwork_TipProportion). We also note that, since, in the
network configuration, cells have fewer neighbours, the values of this metric
are higher than those computed for a cell monolayer

LABEL:supp-ResultsVascularNetwork_Clusters shows the size distribution of
Delta-high cell clusters for simulations on the branching network. We note
that, for this configuration, isolated Delta-high cells (i.e. cells not
adjacent to another Delta-high cell) are predominant in the final spatial
patterns and the patterns generated by the three frameworks are comparable.

Regarding the computational cost (see technical specifications of computers
used in File S1), the CG method showed a great reduction in the average CPU
time compared to the original stochastic system when performing a single
realisation (see LABEL:supp-ResultsCPU). Whereas the numerical cost of
simulations of the full stochastic system (LABEL:supp-ResultsCPU, left panels)
increases exponentially as the system size, $\Omega$, grows, simulations of
the CG system decrease in average computational time as $\Omega$ increases
(LABEL:supp-ResultsCPU, middle panels). This is because, as the noise level
decreases (i.e. $\Omega$ increases), fewer transitions occur in a CG
simulation for a fixed final simulation time. Interestingly, the CG
simulations are also faster than the mean-field system (LABEL:supp-ResultsCPU,
right panels). The numerical integration of the mean-field system (we used the
explicit scheme for the Euler-Lagrange method, although other schemes for
numerical integration may show better performance) required evaluation of the
non-linear right-hand-side of the equations of the mean-field description (see
LABEL:supp-MulticellularMeanField in Supplementary Material) at each time step
for every voxel in the lattice, whereas for the CG simulations only one voxel
undergoes a change (i.e. a phenotype switch) at each iteration.

To summarise, the CG method, while preserving stochasticity of transitions
between cell phenotypes and producing spatial patterns comparable to those
generated using the original stochastic and mean-field descriptions,
significantly reduces computational time of simulations.

##  5 Discussion and conclusions

Hybrid (multiscale) models of complex biological phenomena are often
computationally inefficient, which hinders their potential utility. To address
this issue, we have developed a coarse-graining (CG) method that reduces the
numerical cost of simulations of multi-agent stochastic systems with multiple
stable steady states. The CG technique is based on large deviation theory that
allows the dynamics of a stochastic system to be reduced to a jump process
(i.e. a continuous time Markov chain) on a discrete state space which
comprises the stable steady states of all agents in the system. The CG system
operates on a timescale on which transitions between these steady states take
place. This allows the method to be applied to models whose dynamics act on
timescales longer than the typical timescale for relaxation to an equilibrium
(e.g., molecular or subcellular processes act on longer timescales when
compared to higher spatial scales such as cell migration, dynamics of
extracellular cues, etc.). Our results show good qualitative and quantitative
agreement between CG simulations and other simulation methods (Figures 13 and
LABEL:supp-ResultsVascularNetwork). Furthermore, the CG algorithm is
numerically more efficient in terms of CPU time even when compared with the
corresponding mean-field simulations (see LABEL:supp-ResultsCPU). Likewise,
the CG framework allows exploration of new emergent properties of the system,
such as long timescale patterns in multicellular systems (Figure 12).

The implementation of the CG method requires pre-calculation of several look-
up tables (for stable steady state solutions of the system that is being
coarse-grained, quasipotential values for transitions between them and the
corresponding prefactor of these transitions) which are used later in
simulations. To do this, the values of model parameters must be fixed (except
for the external variables). However, in order to perform sensitivity analysis
with respect to any specific parameter, this parameter may be added to the set
of external variables (thus, adding a new dimension to the look-up tables).
Since the procedure of pre-calculating the look-up tables is done once, prior
to model simulation, it does not increase the numerical cost of the algorithm.
Likewise, the computational cost of computing the quasipotential via the
geometric minimum action method (gMAM) is independent of the system size,
$\Omega$, and an estimate for the required prefactor can be obtained from
simulations of the full stochastic model for a single value of the system size
parameter, $\Omega$, for which we provided an accurate estimate (see Equation
9 and Figure 7). Then the CG model can be efficiently simulated using the
standard Gillespie algorithm for any value of $\Omega$ (or, equivalently,
noise level, $\epsilon=1/\Omega$).

After introducing the CG method (Section 3), we applied it to a multi-agent
model of phenotypic specification of cells via the VEGF-Delta-Notch signalling
pathway. For this system, we demonstrated how the spatial patterning of cells
with different phenotypes changes as CG transitions between these steady
states (phenotypes) occur (Figure 11). We then confirmed that the patterns
generated by the CG system are quantitatively similar to steady state
configurations of the original stochastic system and the associated mean-field
limit for this model (see Figures 13 and LABEL:supp-ResultsVascularNetwork).
We conclude that the CG method preserves the continuous cell phenotypes and
stochasticity of the original system, while reducing the computational cost of
simulations by several orders of magnitude (as compared to the numerical cost
of simulations of the full stochastic system, see LABEL:supp-ResultsCPU).

In this paper, we used the VEGF-Delta-Notch model to illustrate the benefits
of the CG method. We note, however, that the CG method can be applied to a
wider class of multi-agent models in which the behaviour of the agents is
regulated by stochastic models with multiple stable attractors (e.g. steady
states, limit cycles) and whose dynamics are controlled by external cues (e.g.
morphogens, growth factors, levels of specific ligands/receptors in
neighbouring cells, etc.). Examples of systems with subcellular dynamics which
satisfy the requirements for application of the CG method include fate
specification of cells in intestinal crypts [32, 8], epithelial to mesenchymal
phenotypic transition (and its reverse) in cancer invasion [31] and
development [42], cell differentiation in neurogenesis [22], and a general
class of models describing cell decision switches [27]. These models are
multistable and the timescale of simulations is longer than the timescale of
the relevant subcellular signalling pathway. Nonetheless, the spectrum of
models which are suitable for coarse-graining via the CG algorithm is not
restricted to intracellular signalling pathways in animal cells; other
examples include vegetation patterning in arid ecosystems [33] or plant
morphogenesis mediated via the auxin hormone [1, 21]. The exact implementation
of the CG system for the aforementioned models is beyond the scope of this
paper.

To conclude, the CG method developed in this paper paves the way for a
systematic reduction of the dynamics of a wide class of multistable stochastic
models. It allows for investigation of their behaviour on longer timescales
than is possible with other frameworks (e.g. full stochastic simulations or
deterministic equations). To our knowledge, this is the first example in which
large deviation theory has been used to coarse-grain the dynamics of a multi-
agent system. In future work we intend to further investigate the performance
of the CG method by incorporating the CG system for the VEGF-Delta-Notch
signalling into a multiscale model of angiogenesis [44].

## Data management

All of the computational data output is included in the manuscript and/or in
the supplementary material. The code of the numerical procedures used in this
work is available upon request.

## Supplementary materials

Text

##### Supplementary Material

The file contains a more detailed description of the VEGF-Delta-Notch model,
implementation of the CG method, and additional figures and tables.

##### File S1

Technical specifications of the computers used to perform simulations in this
work.

##### Movie S1

A simulation movie showing different pattern configurations explored by the CG
system in a small 2D cell monolayer. The colour bar indicates the levels of
Delta. For this simulation, the interaction radius and system size were fixed
at $R_{s}=15\mu m$ and $\Omega=100$, respectively; the values of the remaining
parameters were fixed as indicated in LABEL:supp-Params.

##### Movie S2

A simulation movie showing emergence of robust pattern configurations in
simulations of the CG system. The colour bar indicates the levels of Delta. A
single cell, whose dynamics is shown in LABEL:OptimalPattern_DeltaLevel and
LABEL:OptimalPattern_TransitionRate, is outlined in cyan. For this simulation,
the interaction radius and system size were fixed at $R_{s}=15\mu m$ and
$\Omega=1000$, respectively; the values of all remaining parameters were fixed
as indicated in LABEL:supp-Params.

##### Movie S3

A simulation movie showing different pattern configurations explored by the CG
system in a branching network. The colour bar indicates the levels of Delta.
For this simulation, the interaction radius and system size were fixed at
$R_{s}=15\mu m$ and $\Omega=100$, respectively; the values of the remaining
parameters were fixed as indicated in LABEL:supp-Params.

## References

  * [1] V. Baldazzi, N. Bertin, H. de Jong, and M. Génard, Towards multiscale plant models: integrating cellular networks, Trends in plant science, 17 (2012), pp. 728–736.
  * [2] R. Bardini, G. Politano, A. Benso, and S. Di Carlo, Multi-level and hybrid modelling approaches for systems biology, Computational and Structural Biotechnology Journal, 15 (2017), pp. 396–402.
  * [3] K. Bentley, C. A. Franco, A. Philippides, R. Blanco, M. Dierkes, V. Gebala, F. Stanchi, M. Jones, I. M. Aspalter, G. Cagna, et al., The role of differential ve-cadherin dynamics in cell rearrangement during angiogenesis, Nature cell biology, 16 (2014), p. 309.
  * [4] S. Bernard, How to build a multiscale model in biology, Acta biotheoretica, 61 (2013), pp. 291–303.
  * [5] R. Blanco and H. Gerhardt, Vegf and notch in tip and stalk cell selection, Cold Spring Harbor Perspectives in Medicine, 3 (2013), p. a006569.
  * [6] M. Boareto, M. K. Jolly, M. Lu, J. N. Onuchic, C. Clementi, and E. Ben-Jacob, Jagged–delta asymmetry in notch signaling can give rise to a sender/receiver hybrid phenotype, Proceedings of the National Academy of Sciences, 112 (2015), pp. E402–E409.
  * [7] F. Bocci, J. N. Onuchic, and M. K. Jolly, Understanding the principles of pattern formation driven by notch signaling by integrating experiments and theoretical models, Frontiers in Physiology, 11 (2020).
  * [8] P. Buske, J. Galle, N. Barker, G. Aust, H. Clevers, and M. Loeffler, A comprehensive model of the spatio-temporal stem cell and tissue organisation in the intestinal crypt, PLoS Comput Biol, 7 (2011), p. e1001045.
  * [9] H. Byrne and M. Chaplain, Mathematical models for tumour angiogenesis: numerical simulations and nonlinear wave solutions, Bulletin of mathematical biology, 57 (1995), pp. 461–486.
  * [10] M. A. Chaplain, Multiscale modelling of cancer: Micro-, meso-and macro-scales of growth and spread, in Approaching Complex Diseases, Springer, 2020, pp. 149–168.
  * [11] M. Cohen, M. Georgiou, N. L. Stevenson, M. Miodownik, and B. Baum, Dynamic filopodia transmit intermittent delta-notch signaling to drive pattern refinement during lateral inhibition, Developmental cell, 19 (2010), pp. 78–89.
  * [12] J. R. Collier, N. A. Monk, P. K. Maini, and J. H. Lewis, Pattern formation by lateral inhibition with feedback: a mathematical model of delta-notch intercellular signalling, Journal of Theoretical Biology, 183 (1996), pp. 429–446.
  * [13] F. Corson, L. Couturier, H. Rouault, K. Mazouni, and F. Schweisguth, Self-organized notch dynamics generate stereotyped sensory organ patterns in drosophila, Science, 356 (2017).
  * [14] R. de la Cruz, P. Guerrero, J. Calvo, and T. Alarcón, Coarse-graining and hybrid methods for efficient simulation of stochastic multi-scale models of tumour growth, Journal of computational physics, 350 (2017), pp. 974–991.
  * [15] R. De La Cruz, R. Perez-Carrasco, P. Guerrero, T. Alarcon, and K. M. Page, Minimum action path theory reveals the details of stochastic transitions out of oscillatory states, Physical review letters, 120 (2018), p. 128102.
  * [16] T. S. Deisboeck, Z. Wang, P. Macklin, and V. Cristini, Multiscale cancer modeling, Annual review of biomedical engineering, 13 (2011), pp. 127–155.
  * [17] A. Deutsch, P. Friedl, L. Preziosi, and G. Theraulaz, Multi-scale analysis and modelling of collective migration in biological systems, 2020.
  * [18] Y. Du, S. C. Herath, Q.-g. Wang, D.-a. Wang, H. H. Asada, and P. C. Chen, Three-dimensional characterization of mechanical interactions between endothelial cells and extracellular matrix during angiogenic sprouting, Scientific Reports, 6 (2016), p. 21362.
  * [19] M. I. Dykman, E. Mori, J. Ross, and P. Hunt, Large fluctuations and optimal paths in chemical kinetics, The Journal of chemical physics, 100 (1994), pp. 5735–5750.
  * [20] J. Elf and M. Ehrenberg, Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases, Systems Biology, 1 (2004), pp. 230–236.
  * [21] E. Farcot, C. Lavedrine, and T. Vernoux, A modular analysis of the auxin signalling network, PLoS One, 10 (2015), p. e0122231.
  * [22] P. Formosa-Jordan, M. Ibañes, S. Ares, and J.-M. Frade, Lateral inhibition and neurogenesis: novel aspects in motion, International Journal of Developmental Biology, 57 (2013), pp. 341–350.
  * [23] M. I. Freidlin and A. D. Wentzell, Random perturbations, in Random perturbations of dynamical systems, Springer, 1998, pp. 15–43.
  * [24] H. Gerhardt, M. Golding, M. Fruttiger, C. Ruhrberg, A. Lundkvist, A. Abramsson, M. Jeltsch, C. Mitchell, K. Alitalo, D. Shima, et al., Vegf guides angiogenic sprouting utilizing endothelial tip cell filopodia, The Journal of Cell Biology, 161 (2003), pp. 1163–1177.
  * [25] D. T. Gillespie, A general method for numerically simulating the stochastic time evolution of coupled chemical reactions, Journal of Computational Physics, 22 (1976), pp. 403–434.
  * [26] T. Grafke, T. Schäfer, and E. Vanden-Eijnden, Long term effects of small random perturbations on dynamical systems: Theoretical and computational tools, in Recent Progress and Modern Challenges in Applied Mathematics, Modeling and Computational Science, Springer, 2017, pp. 17–55.
  * [27] R. Guantes and J. F. Poyatos, Multistable decision switches for flexible control of epigenetic differentiation, PLoS Comput Biol, 4 (2008), p. e1000235.
  * [28] T. Heck, M.-M. Vaeyens, and H. Van Oosterwyck, Computational models of sprouting angiogenesis and cell migration: towards multiscale mechanochemical models of angiogenesis, Mathematical Modelling of Natural Phenomena, 10 (2015), pp. 108–141.
  * [29] M. Heymann and E. Vanden-Eijnden, The geometric minimum action method: A least action principle on the space of curves, Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Institute of Mathematical Sciences, 61 (2008), pp. 1052–1117.
  * [30] G. L. Hunter, Z. Hadjivasiliou, H. Bonin, L. He, N. Perrimon, G. Charras, and B. Baum, Coordinated control of notch/delta signalling and cell cycle progression drives lateral inhibition-mediated tissue patterning, Development, 143 (2016), pp. 2305–2310.
  * [31] M. K. Jolly, M. Boareto, B. Huang, D. Jia, M. Lu, E. Ben-Jacob, J. N. Onuchic, and H. Levine, Implications of the hybrid epithelial/mesenchymal phenotype in metastasis, Frontiers in oncology, 5 (2015), p. 155.
  * [32] S. K. Kay, H. A. Harrington, S. Shepherd, K. Brennan, T. Dale, J. M. Osborne, D. J. Gavaghan, and H. M. Byrne, The role of the hes1 crosstalk hub in notch-wnt interactions of the intestinal crypt, PLoS computational biology, 13 (2017), p. e1005400.
  * [33] S. Kéfi, M. B. Eppinga, P. C. de Ruiter, and M. Rietkerk, Bistability and regular spatial patterns in arid ecosystems, Theoretical Ecology, 3 (2010), pp. 257–269.
  * [34] P. Macklin, S. McDougall, A. R. Anderson, M. A. Chaplain, V. Cristini, and J. Lowengrub, Multiscale modelling and nonlinear simulation of vascular tumour growth, Journal of mathematical biology, 58 (2009), pp. 765–798.
  * [35] N. A. Monk, J. A. Sherratt, and M. R. Owen, Spatiotemporal patterning in models of juxtacrine intercellular signalling with feedback, Mathematical Models for Biological Pattern Formation, (2001), pp. 165–192.
  * [36] J. M. Osborne, A. Walter, S. Kershaw, G. Mirams, A. Fletcher, P. Pathmanathan, D. Gavaghan, O. Jensen, P. Maini, and H. Byrne, A hybrid approach to multi-scale modelling of cancer, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 368 (2010), pp. 5013–5028.
  * [37] R. D. O’Dea and J. R. King, Continuum limits of pattern formation in hexagonal-cell monolayers, Journal of mathematical biology, 64 (2012), pp. 579–610.
  * [38] R. Perez-Carrasco, P. Guerrero, J. Briscoe, and K. M. Page, Intrinsic noise profoundly alters the dynamics and steady state of morphogen-controlled bistable genetic switches, PLoS computational biology, 12 (2016), p. e1005154.
  * [39] G. L. Poppe Jr, Physical applications of the geometric minimum action method, (2018).
  * [40] K. A. Rejniak and A. R. Anderson, State of the art in computational modelling of cancer, Mathematical medicine and biology: a journal of the IMA, 29 (2012), pp. 1–2.
  * [41] D. M. Roma, R. A. O’Flanagan, A. E. Ruckenstein, A. M. Sengupta, and R. Mukhopadhyay, Optimal path to epigenetic switching, Physical Review E, 71 (2005), p. 011902.
  * [42] Y. Sha, D. Haensel, G. Gutierrez, H. Du, X. Dai, and Q. Nie, Intermediate cell states in epithelial-to-mesenchymal transition, Physical biology, 16 (2019), p. 021001.
  * [43] D. Sprinzak, A. Lakhanpal, L. LeBon, L. A. Santat, M. E. Fontes, G. A. Anderson, J. Garcia-Ojalvo, and M. B. Elowitz, Cis-interactions between notch and delta generate mutually exclusive signalling states, Nature, 465 (2010), p. 86.
  * [44] D. Stepanova, H. M. Byrne, P. K. Maini, and T. Alarcón, A multiscale model of complex endothelial cell dynamics in early angiogenesis, PLOS Computational Biology, 17 (2021), p. e1008055.