#  Pushing coarse-grained models beyond the continuum limit using equation
learning

Daniel J. VandenHeuvel  School of Mathematical Sciences, Queensland University
of Technology, Brisbane, Australia.  Pascal R. Buenzli  School of Mathematical
Sciences, Queensland University of Technology, Brisbane, Australia.  Matthew
J. Simpson  School of Mathematical Sciences, Queensland University of
Technology, Brisbane, Australia.  Corresponding author:
<EMAIL_ADDRESS>
###### Abstract

Mathematical modelling of biological population dynamics often involves
proposing high fidelity discrete agent-based models that capture stochasticity
and individual-level processes. These models are often considered in
conjunction with an approximate coarse-grained differential equation that
captures population-level features only. These coarse-grained models are only
accurate in certain asymptotic parameter regimes, such as enforcing that the
time scale of individual motility far exceeds the time scale of birth/death
processes. When these coarse-grained models are accurate, the discrete model
still abides by conservation laws at the microscopic level, which implies that
there is some macroscopic conservation law that can describe the macroscopic
dynamics. In this work, we introduce an equation learning framework to find
accurate coarse-grained models when standard continuum limit approaches are
inaccurate. We demonstrate our approach using a discrete mechanical model of
epithelial tissues, considering a series of four case studies that consider
problems with and without free boundaries, and with and without proliferation,
illustrating how we can learn macroscopic equations describing mechanical
relaxation, cell proliferation, and the equation governing the dynamics of the
free boundary of the tissue. While our presentation focuses on this biological
application, our approach is more broadly applicable across a range of
scenarios where discrete models are approximated by approximate continuum-
limit descriptions. All code and data to reproduce this work are available at
https://github.com/DanielVandH/StepwiseEQL.jl.

##  1 Introduction

Mathematical models of population dynamics are often constructed by
considering both discrete and continuous descriptions, allowing for both
microscopic and macroscopic details to be considered [1]. This approach has
been applied to several kinds of discrete models, including cellular Potts
models [2, 3, 4, 5], exclusion processes [6, 7, 8, 9], mechanical models of
epithelial tissues [10, 11, 12, 13, 14, 15, 16, 17], hydrodynamics [18, 19],
and a variety of other types of individual-based models [1, 20, 21, 22, 23,
24, 25, 26, 27]. Continuum models are useful for describing collective
behaviour, especially because the computational requirement of discrete models
increases with the size of the population, and this can become computationally
prohibitive for large populations, which is particularly problematic for
parameter inference [28]. In contrast, the computational requirement to solve
a continuous model is independent of the population size, and generally
requires less computational overhead than working with a discrete approach
only [15]. Continuum models are typically obtained by coarse-graining the
discrete model, using Taylor series expansions to obtain continuous partial
differential equation (PDE) models that govern the population densities on a
continuum or macroscopic scale [10, 11, 29, 30].

One challenge with using coarse-grained continuum limit models is that while
the solution of these models can match averaged data from the corresponding
discrete model for certain choices of parameters [10, 17, 31], the solution of
the continuous model can be a very poor approximation for other parameter
choices [13, 15, 32, 33]. More generally, coarse-grained models are typically
only valid in certain asymptotic parameter regimes [34, 31, 33]. For example,
suppose we have a discrete space, discrete time, agent-based model that
incorporates random motion and random proliferation. Random motion involves
stepping a distance $\Delta$ with probability $P_{\mathrm{m}}\in[0,1]$ per
time step of duration $\tau$. The stochastic proliferation process involves
undergoing proliferation with probability $P_{\mathrm{p}}\in[0,1]$ per unit
time step. The continuum limit description of this kind of discrete process
can be written as [33]

$\partialderivative{q}{t}=\partialderivative{x}\left(D(q)\partialderivative{q}{x}\right)+R(q),$
(1)

where $q$ is the macroscopic density of individuals, $D(q)$ is the nonlinear
diffusivity that describes the effects of individual migration, and $R(q)$ is
a source term that describes the effects of the birth process in the discrete
model [33]. Standard approaches to derive (1) require
$D(q)=\mathcal{O}(P_{\mathrm{m}}\Delta^{2}/\tau)$ and
$R(q)=\mathcal{O}(P_{\mathrm{p}}/\tau)$ in the limit that $\Delta\to 0$ and
$\tau\to 0$. To obtain a well-defined continuum limit such that the diffusion
and source terms are both present in the macroscopic model, some restrictions
on the parameters in the discrete model are required [34, 33]. Typically, this
is achieved by taking the limit as $\Delta\to 0$ and $\tau\to 0$ jointly such
that the ratio $\Delta^{2}/\tau$ remains finite, implying that
$P_{\mathrm{p}}=\mathcal{O}(\tau)$ so that both the diffusion and source terms
in (1) are $\mathcal{O}(1)$. In practice, this means that the time scale of
individual migration events has to be much faster than the time scale of
individual proliferation events, otherwise the continuum limit description is
not well defined [34, 33]. If this restriction is not enforced, then the
solution of the continuum limit model does not always predict the averaged
behaviour of the discrete model [33], as the terms on the right-hand side of
(1) are no longer $\mathcal{O}(1)$ so that the continuum limit is not well
defined [34].

Regardless of whether choices of parameters in a discrete model obey the
asymptotic restrictions imposed by coarse-graining, the discrete model still
obeys a conservation principle, which implies that there is some alternative
macroscopic conservation description that will describe population-level
features of interest [35, 36]. Equation learning is a means of determining
appropriate continuum models outside of the usual continuum limit asymptotic
regimes. Equation learning has been used in several applications for model
discovery. In the context of PDEs, a typical approach is to write $\partial
q/\partial t=\mathcal{N}(q,\mathcal{D},\bm{\theta})$, where $q$ is the
population density, $\mathcal{N}$ is some nonlinear function parametrised by
$\bm{\theta}$, $\mathcal{D}$ is a collection of differential operators, and
$\bm{\theta}$ are parameters to be estimated [37]. This formulation was first
introduced by Rudy et al. [37], who extended previous work in learning
ordinary differential equations (ODEs) proposed by Brunton et al. [38].
Equation learning methods developed for the purpose of learning biological
models has also been a key interest [39, 40]. Lagergren et al. [41] introduce
a biologically-informed neural network framework that uses equation learning
that is guided by biological constraints, imposing a specific conservation PDE
rather than a general nonlinear function $\mathcal{N}$. Lagergren et al. [41]
use this framework to discover a model describing data from simple in vitro
experiments that describe the invasion of populations of motile and
proliferative prostate cancer cells. VandenHeuvel et al. [42] extend the work
of Lagergren et al. [41], incorporating uncertainty quantification into the
equation learning procedure through a bootstrapping approach. Nardini et al.
[32] use discrete data from agent-based models to learn associated continuum
ODE models, combining a user-provided library of functions together with
sparse regression methods to give simple ODE models describing population
densities. Regression methods have also been used as an alternative to
equation learning for this purpose [43].

These previous approaches to equation learning consider various methods to
estimate the parameters $\bm{\theta}$, such as sparse regression or nonlinear
optimisation [37, 38, 39, 41, 42, 32], representing $\mathcal{N}$ as a library
of functions [37, 38, 39], neural networks [41], or in the form of a
conservation law with individual components to be learned [41, 42]. In this
work, we introduce a stepwise equation learning framework, inspired from
stepwise regression [44], for estimating $\bm{\theta}$ from averaged discrete
data with a given $\mathcal{N}$ representing a proposed form for the continuum
model description. We incorporate or remove terms one at a time until a
parsimonious continuum model is obtained whose solution matches the data well
and no further improvements can be made to this match. Our approach is
advantageous for several reasons. Firstly, it is computationally efficient and
parallelisable, allowing for rapid exploration of results with different
discrete parameters and different forms of $\mathcal{N}$ for a given data set.
Secondly, the approach is modular, with different mechanistic features easily
incorporated. This approach enables extensive computational experimentation by
exploring the impact of including or excluding putative terms in the continuum
model without any great increase in computational overhead. Lastly, it is easy
to examine the results from our procedure, allowing for ease of diagnosing and
correcting reasons for obtaining poor fitting models, and explaining what
components of the continuum model are the most influential. We emphasise that
a key difference between our approach and other work, such as the methods
developed by Brunton et al. [38] and Rudy et al. [37], is that we constrain
our problem so that we can only learn conservation laws rather than allow a
general form through a library of functions, and that we iteratively eliminate
variables from $\bm{\theta}$ rather than use sparse regression. These
important features are what support the modularity and interpretability of our
approach.

To illustrate our procedure, we consider a discrete, individual-based one-
dimensional toy model inspired from epithelial tissues [17, 10]. Epithelial
tissues are biological tissue composed of cells, organised in a monolayer, and
are present in many parts of the body and interact with other cells [45],
lining surfaces such as the skin and the intestine [46]. They are important in
a variety of contexts, such as wound healing [47, 48] and cancer [49, 50].
Many models have been developed for studying their dynamics, considering both
discrete and continuum modelling [10, 11, 12, 13, 14, 15, 16], with most
models given in the form of a nonlinear reaction-diffusion equation with a
moving boundary, using a nonlinear diffusivity term to incorporate mechanical
relaxation and a source term to model cell proliferation [12, 13, 16]. These
continuum limit models too are only accurate in certain parameter regimes,
becoming inaccurate if the rate of mechanical relaxation is slow relative to
the rate of proliferation [13, 15, 33]. To apply our stepwise equation
learning procedure, we let the nonlinear function $\mathcal{N}$ be given in
the form of a conservation law together with equations describing the free
boundary. We demonstrate this approach using a series of four biologically-
motivated case studies, considering problems with and without a free boundary,
and with and without proliferation,with each case study building on those
before it. The first two case studies are used to demonstrate how our approach
can learn known continuum limits, while the latter two case studies show how
we can learn improved continuum limit models in parameter regimes where these
known continuum limits are no longer accurate. We implement our approach in
the Julia language [51], and all code is available on GitHub at
https://github.com/DanielVandH/StepwiseEQL.jl.

##  2 Mathematical model

Figure 1: Discrete model and schematics for each case study (CS). (a) A fixed
boundary problem with $x_{1}=0$ and $x_{n}=L$ fixed. (b) A free boundary
problem with $x_{1}=0$ and $x_{n}(t)=L(t)$, show in red, free. (c)
Proliferation schematic, showing a cell $(x_{i}(t),x_{i+1}(t))$ dividing into
$(x_{i}(t+\Delta t),x_{i+1}(t+\Delta t))$ and $(x_{i+1}(t+\Delta
t),x_{i+2}(t+\Delta t))$ following a proliferation event, where
$x_{i+1}(t+\Delta t)=(x_{i}(t)+x_{i+1}(t))/2$. (d)–(g) show schematics for the
four case studies considered in the paper, where the first row in each panel
is a representation of the initial configuration of cells at $t=0$ and the
second row a representation at a later time $t>0$.

Following Murray et al. and Baker et al. [16, 10], we suppose that we have a
set of nodes $x_{1},\ldots,x_{n}(t)$ describing $n$ cell boundaries at a time
$t$. The interval $(x_{i}(t),x_{i+1}(t))$ represents the $i$th cell for
$i=1,\ldots,n-1$, where we fix $x_{1}=0$ and $x_{1}<x_{2}(t)<\cdots<x_{n}(t)$.
The number of nodes, $n$, may increase over time due to cell proliferation. We
model the mechanical interaction between cells by treating them as springs, as
indicated in Figure 1, so that each node $i$ experiences forces $F_{i,i\pm 1}$
from nodes $i\pm 1$, respectively, except at the boundaries where there is
only one neighbouring force. We further assume that each of these springs has
the same mechanical properties, and that the viscous force from the
surrounding medium is given by $\eta\mathrm{d}x_{i}(t)/\mathrm{d}t$ with drag
coefficient $\eta$. Lastly, assuming we are in a viscous medium so that the
motion is overdamped, we can model the dynamics of each individual node
$x_{i}(t)$, fixing $x_{1}=0$, by [16]

$\displaystyle\eta\derivative{x_{i}(t)}{t}$
$\displaystyle=F_{i,i-1}+F_{i,i+1},\quad i=1,\ldots,n-1,$ (2)
$\displaystyle\eta\derivative{x_{n}(t)}{t}$ $\displaystyle=F_{n,n-1},$ (3)

where

$F_{i,i\pm 1}=F\left(|x_{i}(t)-x_{i\pm 1}(t)|\right)\frac{x_{i}(t)-x_{i\pm
1}(t)}{|x_{i}(t)-x_{i\pm 1}(t)|}$ (4)

is the interaction force that the $i$th node experiences from nodes $i\pm 1$
(Figure 1). In Case Studies 1 and 3 (see Section 3, below), we hold
$x_{n}(t)=L$ constant and discard (3). Throughout this work, we use linear
Hookean springs so that $F(\ell_{i})=k(s-\ell_{i})$, $\ell_{i}>0$, where
$\ell_{i}(t)=x_{i+1}(t)-x_{i}(t)$ is the length of the $i$th cell, $k>0$ is
the spring constant, and $s\geq 0$ is the resting spring length [10]; we
discuss other force laws in Appendix E.

The dynamics governed by (2)–(3) describe a system in which cells mechanically
relax. Following previous work [10, 16, 14], we introduce a stochastic
mechanism that allows the cells to undergo proliferation, assuming only one
cell can divide at a time over a given interval $[t,t+\Delta t)$ for some
small duration $\Delta t$. We let the probability that the $i$th cell
proliferates be given by $G_{i}\Delta t$, where $G_{i}=G(\ell_{i})$ for some
length-dependent proliferation law $G(\ell_{i})>0$. As represented in Figure
1(c), when the $i$th cell proliferates, the cell divides into two equally-
sized daughter cells, and the boundary between the new daughter cells is
placed at the midpoint of the original cell. Throughout this work, we use a
logistic proliferation law $G(\ell_{i})=\beta[1-1/(K\ell_{i})]$ with
$\ell_{i}>1/K$, where $\beta$ is the intrinsic proliferation rate and $K$ is
the carrying capacity density; we consider other proliferation laws in
Appendix E. The implementation of the solution to these equations (2)–(3) and
the proliferation mechanism is given in the Julia package
EpithelialDynamics1D.jl; in this implementation, if $G(\ell_{i})<0$ we set
$G(\ell_{i})=0$ to be consistent with the fact that we interpret $G(\ell_{i})$
as a probability. We emphasise that, without proliferation, we need only solve
(2)–(3) once for a given initial condition in order to obtain the expected
behaviour of the discrete model, because the discrete model is deterministic
in the absence of proliferation. In contrast, incorporating proliferation
means that we need to consider several identically-prepared realisations of
the same stochastic discrete model to estimate the expected behaviour of the
discrete model for a given initial condition.

In practice, macroscopic models of populations of cells are described in terms
of cell densities rather than keeping track of the position of individual cell
boundaries. The density of the $i$th cell $(x_{i}(t),x_{i+1}(t))$ is
$1/\ell_{i}(t)$. For an interior node $x_{i}(t)$, we obtain a density
$q_{i}(t)$ by taking the inverse of the average of the cells left and right of
$x_{i}(t)$, giving

$q_{i}(t)=\frac{2}{x_{i+1}(t)-x_{i-1}(t)},\quad i=2,\ldots,n-1,$ (5)

as in Baker et al. [16]. At boundary nodes, we use

$q_{1}(t)=\frac{2}{x_{2}(t)}-\frac{2}{x_{3}(t)},\quad
q_{n}(t)=\frac{2}{x_{n}(t)-x_{n-1}(t)}-\frac{2}{x_{n}(t)-x_{n-2}(t)},$ (6)

derived by linear extrapolation of (5) to the boundary. The densities in (6)
ensure that the slope of the density curves at the boundaries, $\partial
q/\partial x$, match those in the continuum limit. We discuss the derivation
of (6) in Appendix B. In the continuum limit where the number of cells is
large and mechanical relaxation is fast, the densities evolve according to the
moving boundary problem [10, 16]

$\begin{array}[]{rcll}\partialderivative{q}{t}&=&\partialderivative{x}\left(D(q)\partialderivative{q}{x}\right)+R(q)&0<x<L(t),\,t>0,\\\\[6.0pt]
\partialderivative{q}{x}&=&0&x=0,\,t>0,\\\\[6.0pt]
\partialderivative{q}{x}&=&H(q)&x=L(t),\,t>0,\\\\[6.0pt]
q\derivative{L}{t}&=&-D(q)\partialderivative{q}{x}&x=L(t),\,t>0,\\\\[6.0pt]
\end{array}$ (7)

where $q(x,t)$ is the density at position $x$ and time $t$, $D(q)=-1/(\eta
q^{2})F^{\prime}(1/q)$, $R(q)=qG(1/q)$, $H(q)=-2qF(1/q)/[\eta D(q)]$, and
$L(t)=x_{n}(t)$ is the leading edge position with $L(0)=x_{n}(0)$. The
quantity $1/q$ in these equations can be interpreted as a continuous function
related to the length of the individual cells. The initial condition
$q(x,0)=q_{0}(x)$ is a linear interpolant of the discrete densities $q_{i}(t)$
of the cells at $t=0$. Similar to the discussion of (1), for this continuum
limit to be valid so that both $D(q)$ and $R(q)$ play a role in the continuum
model, constraints must be imposed on the discrete parameters. As discussed by
Murphy et al. [15], we require that the time scale of mechanical relaxation is
sufficiently fast relative to the time scale of proliferation. In practice
this means that for a given choice of $\beta$ we must have $k/\eta$
sufficiently large for the solution of the continuum model to match averaged
data from the discrete model. We note that, with our choices of $F$ and $G$,
the functions in (7) are given by

$D(q)=\frac{k}{\eta q^{2}},\quad R(q)=\beta q\left(1-\frac{q}{K}\right),\quad
H(q)=2q^{2}(1-qs).$ (8)

For fixed boundary problems we take $H(q)=0$ and $\mathrm{d}L/\mathrm{d}t=0$.
In Appendix C, we describe how to solve (7) numerically, as well as how to
solve the corresponding problem with fixed boundaries numerically.

##  3 Continuum-discrete comparison

Figure 2: Space-time diagrams (left column) and densities (right column) for
the four case studies from Figure 1 considered throughout the paper. The left
column shows the evolution of the discrete densities in space and time, with
(c) and (d) showing averaged results over $2500$ identically-prepared
realisations of the discrete model. In (b) and (d), the red line shows the
position of the free boundary. In the figures in the right column, the solid
curves are the discrete densities (5) and the dashed curves are solutions to
the continuum limit problem (7), and the curves are given by black, red, blue,
green, orange, and purple in the order of increasing time as indicated by the
black arrows. The times shown are (a) $t=0,1,2,3,4,5$; (b)
$t=0,5,10,25,50,100$; (c) $t=0,1,5,10,20,50$; and (d) $t=0,5,10,20,50,100$. In
(c) and (d), the shaded regions show $95\%$ confidence bands from the mean
discrete curves at each time; the curves in (a) and (b) show no shaded regions
as these models have no stochasticity.

We now consider four biologically-motivated case studies to illustrate the
performance of the continuum limit description (7). These case studies are
represented schematically in Figure 1(d)—(g). Case Studies 1 and 3, shown in
Figure 1(d) and Figure 1(f), are fixed boundary problems, where we see cells
relax mechanically towards a steady state where each cell has equal length.
Case Studies 2 and 4 are free boundary problems, where the right-most cell
boundary moves in the positive $x$-direction while all cells relax towards a
steady state where the length of each cell is given by resting spring length
$s$. Case Studies 1 and 2 have $\beta=0$ so that there is no cell
proliferation and the number of cells remains fixed during the simulations,
whereas Case Studies 3 and 4 have $\beta>0$ so that the number of cells
increases during the discrete simulations. To explore these problems, we first
consider cases where the continuum limit model is accurate, using the data
shown in Figure 2, where we show space-time diagrams and a set of averaged
density profiles for each problem in the left and right columns of Figure 2,
respectively. Case Studies 1 and 3 initially place $30$ nodes in $0\leq x\leq
5$ and $30$ nodes in $25\leq x\leq 30$, or equivalently $n=60$ with $28$ cells
in $0\leq x\leq 5$ and $28$ cells in $25\leq x\leq 30$, spacing the nodes
uniformly within each subinterval. Case Studies 2 and 4 initially place $60$
equally spaced nodes in $0\leq x\leq 5$.

Figure 3: Examples of inaccurate continuum limits for (a) Case Study 3 and (b)
Case Study 4, where both case studies use the same parameters as in Figure 2
except with $k=1/5$ rather than $k=50$. The solid curves are the discrete
densities (5) and the dashed curves are solutions to the continuum limit
problem (7). The arrows show the direction of increasing time. The density
profiles are plotted in black, red, blue, green, orange, and purple for the
respective times (a) $t=0,1,10,25,40,75$ and (b) $t=0,5,25,50,100,250$.

The problems shown in Figure 2 use parameter values such that the solution of
the continuum limit (7) is a good match to the averaged discrete density
profiles. In particular, all problems use $k=50$, $\eta=1$, $s=1/5$ and, for
Case Studies 3 and 4, $\Delta t=10^{-2}$, $K=15$, and $\beta=0.15$. The
accuracy of the continuum limit is clearly evident in the right column of
Figure 2 where, in each case, the solution of the continuum limit model is
visually indistinguishable from averaged data from the discrete model. With
proliferation, however, the continuum limit can be accurate when $k/\eta$ is
not too much larger than $\beta$, and we use Case Studies 3 and 4 to explore
this.

Figure 3 shows further continuum-discrete comparisons for Case Studies 3 and 4
where we have slowed the mechanical relaxation by taking $k=1/5$. This choice
of $k$ means that $D(q)$ and $R(q)$ are no longer on the same scale and thus
the continuum limit is no longer well defined, as explained in the discussion
of (1), meaning the continuum limit solutions are no longer accurate. In both
cases, the solution of the continuum limit model lags behind the averaged data
from the discrete model. In Appendix A, we show the $95\%$ confidence regions
for each curve in Figure 3, where we find that the solutions have much greater
variance compared to the corresponding curves in Figure 2 where $k=50$.

We are interested in developing an equation learning method for learning an
improved continuum model for problems like those in Figure 3, allowing us to
extend beyond the parameter regime where the continuum limit (7) is accurate.
We demonstrate this in Case Studies 1–4 in Section 4 where we develop such a
method.

##  4 Learning accurate continuum limit models

In this section we introduce our method for equation learning and demonstrate
the method using the four case studies from Figures 1–3. Since the equation
learning procedure is modular, adding these components into an existing
problem is straightforward. All Julia code to reproduce these results is
available at https://github.com/DanielVandH/StepwiseEQL.jl. A summary of all
the parameters used for each case study is given in Table 1.

Parameter | Case Study
---|---
1 | 2 | 3a | 3b | 4a | 4b
$k$ | $50$ | $50$ | $50$ | $1/5$ | $50$ | $1/5$
$\eta$ | $1$ | $1$ | $1$ | $1$ | $1$ | $1$
$s$ | $1/5$ | $1/5$ | $1/5$ | $1/5$ | $1/5$ | $1/5$
$\Delta t$ | — | — | $10^{-2}$ | $10^{-2}$ | $10^{-2}$ | $10^{-2}$
$\beta$ | — | — | $0.15$ | $0.15$ | $0.15$ | $0.15$
$K$ | — | — | $15$ | $15$ | $15$ | $15$
$M$ | $50$ | $150$ | $501$ | $751$ | $(25,50,100,250)$ | $(20,200,200,200)$
$t_{1}$ | $0$ | $0$ | $0$ | $0$ | $(0,0,5,10)$ | $(0,2,10,20)$
$t_{M}$ | $5$ | $15$ | $50$ | $75$ | $(10^{-1},5,10,50)$ | $(2,10,20,50)$
$n_{s}$ | — | — | $1000$ | $1000$ | $1000$ | $1000$
$n_{k}$ | — | — | $50$ | $200$ | $(25,50,100,50)$ | $(50,100,100,100)$
$\tau_{q}$ | $0.1$ | $0.35$ | $0.1$ | $0.25$ | $(0.1,0,0,0)$ | $(0,0,0,0.3)$
$\tau_{\mathrm{d}L/\mathrm{d}t}$ | — | $0.1$ | — | — | $(0,0.2,0,0)$ | $(0,0.4,0,0)$
$\tau_{t}$ | $0$ | $0$ | $0$ | $0$ | $0$ | $(0.4,0.4,0,0)$
Table 1: Parameters used for each case study. The parameters are $k$, the
spring constant; $\eta$, the drag coefficient; $s$, the resting spring length;
$\Delta t$, the proliferation duration; $\beta$, the intrinsic proliferation
rate; $K$, the carrying capacity density; $M$, the number of time points;
$t_{1}$, the initial time; $t_{M}$, the final time; $n_{s}$, the number of
identically-prepared realisations; $n_{k}$, the number of knots used for
averaging over realisations; $\tau_{q}$, which defines the $100\tau_{q}\%$ and
$100(1-\tau_{q})\%$ density quantiles; $\tau_{\mathrm{d}L/\mathrm{d}t}$, which
defines the $100\tau_{\mathrm{d}L/\mathrm{d}t}\%$ and
$100(1-\tau_{\mathrm{d}L/\mathrm{d}t})\%$ velocity quantiles; and $\tau_{t}$,
which defines the $100\tau_{t}\%$ and $100(1-\tau_{t})\%$ temporal quantiles.
Values indicated by a line are not relevant for the corresponding case study.
For Case Study 3 and 4, the label “a” refers to the accurate continuum limit
case, and “b” refers to the inaccurate continuum limit case. For Case Study 4,
some parameters are given by a set of four parameters, with the $i$th value of
this set referring to the value used when learning the $i$th mechanism; see
Section 44.4 for details.

###  4.1 Case Study 1: Fixed boundaries

Case Study 1 involves mechanical relaxation only so that there is no cell
proliferation and the boundaries are fixed, implying $R(q)=0$ and $H(q)=0$ in
(7), respectively, and the only function to learn is $D(q)$.

Our equation learning approach starts by assuming that $D(q)$ is a linear
combination of $d$ basis coefficients $\\{\theta_{1},\ldots,\theta_{d}\\}$ and
$d$ basis functions $\\{\varphi_{1},\ldots,\varphi_{d}\\}$, meaning $D(q)$ can
be represented as

$D(q)=\sum_{i=1}^{d}\theta_{i}\varphi_{i}(q).$ (9)

These basis functions could be any univariate function of $q$, for example the
basis could be
$\\{\varphi_{1},\varphi_{2},\varphi_{3}\\}=\\{1/q,1/q^{2},1/q^{3}\\}$ with
$d=3$. In this work, we impose the constraint that $D(q)\geq 0$ for
$q_{\min}\leq q\leq q_{\max}$, where $q_{\min}$ and $q_{\max}$ are the minimum
and maximum densities observed in the discrete simulations, respectively. This
constraint enforces the condition that the nonlinear diffusivity function is
positive over the density interval of interest. While it is possible to work
with some choices of nonlinear diffusivity functions for which $D(q)<0$ for
some interval of density [52, 53, 54], we wish to avoid the possibility of
having negative nonlinear diffusivity functions and our results support this
approach.

The aim is to estimate
$\bm{\theta}=(\theta_{1},\ldots,\theta_{d})^{\mkern-1.5mu\mathsf{T}}$ in (9).
We use ideas similar to the basis function approach from VandenHeuvel et al.
[42], using (9) to construct a matrix problem for $\bm{\theta}$. In
particular, let us take the PDE (7), with $R(q)=0$ and $H(q)=0$, and expand
the spatial derivative term so that we can isolate the $\theta_{k}$ terms,

$\partialderivative{q_{ij}}{t}=\sum_{k=1}^{d}\left\\{\derivative{\varphi_{k}(q_{ij})}{q}\left(\partialderivative{q_{ij}}{x}\right)^{2}+\varphi_{k}(q_{ij})\partialderivative[2]{q_{ij}}{x}\right\\}\theta_{k},$
(10)

where we let $q_{ij}$ denote the discrete density at position
$x_{ij}=x_{i}(t_{j})$ and time $t_{j}$. We note that while $q_{ij}$ is
discrete, we assume it can be approximated by a smooth function, allowing us
to define these derivatives $\partial q_{ij}/\partial t$, $\partial
q_{ij}/\partial x$, and $\partial^{2}q_{ij}/\partial x^{2}$ in (10); this
assumption is appropriate since, as shown in Figure 2, these discrete
densities can be well approximated by smooth functions. These derivatives are
estimated using finite differences, as described in Appendix D. We also
emphasise that, while (10) appears similar to results in [38, 37], the crucial
difference is that we are specifying forms for the mechanisms of the PDE
rather than the complete PDE itself; one other important difference is in how
we estimate $\bm{\theta}$, defined below and in (15). We save the solution to
the discrete problems (2)–(3) at $M$ times $0=t_{1}<t_{2}<\cdots<t_{M}$ so
that $i\in\\{1,\ldots,n\\}$ and $j\in\\{2,\ldots,M\\}$, where $n=60$ is the
number of nodes and we do not deal with data at $j=1$ since the PDE does not
apply at $t=0$. We can therefore convert (10) into a rectangular matrix
problem $\mathbf{A}\bm{\theta}=\mathbf{b}$, where the $r$th row in
$\mathbf{A}$, $r=1,\ldots,n(M-1)$, corresponding to the point $(x_{ij},t_{j})$
is given by $\mathbf{a}_{ij}\in\mathbb{R}^{1\times d}$, where

$\mathbf{a}_{ij}=\begin{bmatrix}\derivative{\varphi_{1}(q_{ij})}{q}\left(\partialderivative{q_{ij}}{x}\right)^{2}+\varphi_{1}(q_{ij})\partialderivative[2]{q_{ij}}{x},&\cdots,&\derivative{\varphi_{d}(q_{ij})}{q}\left(\partialderivative{q_{ij}}{x}\right)^{2}+\varphi_{d}(q_{ij})\partialderivative[2]{q_{ij}}{x}\end{bmatrix},$
(11)

with each element of $\mathbf{a}_{ij}$ corresponding to the contribution of
the associated basis function in (10). Thus, we obtain the system

$\mathbf{A}=\begin{bmatrix}\mathbf{a}_{12}\\\ \mathbf{a}_{22}\\\ \vdots\\\
\mathbf{a}_{nM}\end{bmatrix}\in\mathbb{R}^{n(M-1)\times
d}\quad\textnormal{and}\quad\mathbf{b}=\begin{bmatrix}\partial q_{12}/\partial
t\\\ \partial q_{22}/\partial t\\\ \vdots\\\ \partial q_{nM}/\partial
t\end{bmatrix}\in\mathbb{R}^{n(M-1)\times 1}.$ (12)

The solution of $\mathbf{A}\bm{\theta}=\mathbf{b}$, given by
$\bm{\theta}=(\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{A})^{-1}\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{b}$,
is obtained by minimising the residual
$\|\mathbf{A}\bm{\theta}-\mathbf{b}\|_{2}^{2}$, which keeps all terms present
in the learned model. We expect, however, just as in (8), that $\bm{\theta}$
is sparse so that $D(q)$ has very few terms, which makes the interpretation of
these terms feasible [37, 38]. There are several ways that we could solve
$\mathbf{A}\bm{\theta}=\mathbf{b}$ to obtain a sparse vector, such as with
sparse regression [37], but in this work we take a stepwise equation learning
approach inspired by stepwise regression [44] as this helps with both the
exposition and modularity of our approach. For this approach, we first let
$\mathcal{I}=\\{1,\ldots,d\\}$ be the set of basis function indices. We let
$\mathcal{A}_{k}$ denote the set of active coefficients at the $k$th
iteration, meaning the indices of non-zero values in $\bm{\theta}$, starting
with $\mathcal{A}_{1}=\mathcal{I}$. The set of indices of zero values in
$\bm{\theta}$, $\mathcal{I}_{k}=\mathcal{I}\setminus\mathcal{A}_{k}$, is
called the set of inactive coefficients. To obtain the next set,
$\mathcal{A}_{k+1}$, from a current set $\mathcal{A}_{k}$, we apply the
following steps:

  1. 1.

Let the vector $\bm{\theta}_{\mathcal{A}}$ denote the solution to
$\mathbf{A}\bm{\theta}=\mathbf{b}$ subject to the constraint that each
inactive coefficient $\theta_{i}$ is zero, meaning $\theta_{i}=0$ for
$i\in\mathcal{I}\setminus\mathcal{A}$ for a given set of active coefficients
$\mathcal{A}$. We compute $\bm{\theta}_{\mathcal{A}}$ by solving the reduced
problem in which the inactive columns of $\mathbf{A}$ are not included. The
vector with $\mathcal{A}=\mathcal{A}_{k}$ at step $k$ is denoted
$\bm{\theta}_{k}$. With this definition, we compute the sets

$\mathcal{M}_{k}^{+}=\left\\{\bm{\theta}_{\mathcal{A}_{k}\cup\\{i\\}}:i\notin\mathcal{A}_{k}\right\\}\quad\textnormal{and}\quad\mathcal{M}_{k}^{-}=\left\\{\bm{\theta}_{\mathcal{A}_{k}\setminus\\{i\\}}:i\in\mathcal{A}_{k}\right\\}.$
(13)

$\mathcal{M}_{k}^{+}$ is the set of all coefficient vectors $\bm{\theta}$
obtained by making each active coefficient at step $k$ inactive one at a time.
$\mathcal{M}_{k}^{-}$, is similar to $\mathcal{M}_{k+1}^{-}$ except we make
each inactive coefficient at step $k$ active one at a time. We then define
$\mathcal{M}_{k}=\\{\bm{\theta}_{k}\\}\cup\mathcal{M}_{k}^{+}\cup\mathcal{M}_{k}^{-}$,
so that $\mathcal{M}_{k}$ is the set of all coefficient vectors obtained from
activating coefficients one at a time, deactivating coefficients one at a
time, or retaining the current vector $\bm{\theta}_{k}$.

  2. 2.

Choose one of the vectors in $\mathcal{M}_{k}$ by defining a loss function
$\mathcal{L}(\bm{\theta})$:

$\underbrace{\mathcal{L}(\bm{\theta})}_{\textnormal{loss}}=\underbrace{\log\left[\frac{1}{n(M-1)}\sum_{j=2}^{M}\sum_{i=1}^{n}\left(\frac{q_{ij}-q(x_{ij},t_{j};\bm{\theta})}{q_{ij}}\right)^{2}\right]}_{\textnormal{goodness
of fit}}+\underbrace{\|\bm{\theta}\|_{0}}_{\text{model complexity}},$ (14)

where $q(x,t;\bm{\theta})$ is the solution of the PDE (7) with $R(q)=H(q)=0$
and $D(q)$ uses the coefficients $\bm{\theta}$ in (9),
$q(x_{ij},t_{j};\bm{\theta})$ is the linear interpolant of the PDE data at
$t=t_{j}$ evaluated at $x=x_{ij}$, and $\|\bm{\theta}\|_{0}$ is the number of
non-zero terms in $\bm{\theta}$. This loss function balances the goodness of
fit with model complexity. If, for some $\bm{\theta}$, $D(q)<0$ within
$q_{\min}\leq q\leq q_{\max}$, which we check by evaluating $D(q)$ at
$n_{c}=100$ equally spaced points in $q_{\min}\leq q\leq q_{\max}$, we set
$\mathcal{L}(\bm{\theta})=\infty$. With this loss function, we compute the
next coefficient vector

$\bm{\theta}_{k+1}=\operatorname*{argmin}_{\bm{\theta}\in\mathcal{M}_{k}}\mathcal{L}(\bm{\theta}).$
(15)

If $\bm{\theta}_{k+1}=\mathbf{0}$, so that all the coefficients are inactive,
we instead take the vector that attains the second-smallest loss so that a
model with no terms cannot be selected.

  3. 3.

If $\bm{\theta}_{k+1}=\bm{\theta}_{k}$, then there are no more local
improvements to be made and so the procedure stops. Otherwise, we recompute
$\mathcal{A}_{k+1}$ and $\mathcal{I}_{k+1}$ from $\bm{\theta}_{k+1}$ and
continue iterating.

The second step prevents empty models from being returned, allowing the
algorithm to more easily find an optimal model when starting with no active
coefficients. We note that Nardini et al. [32] consider a loss based on the
regression error, $\|\mathbf{A}\bm{\theta}-\mathbf{b}\|_{2}^{2}$, that has
been useful for a range of previously-considered problems [38, 37, 32]. We do
not consider the regression error in this work as we find that it typically
leads to poorer estimates for $\bm{\theta}$ compared to controlling the
density errors as we do in (15).

Let us now apply the procedure to our data from Figure 2, where we know that
the continuum limit with $D(q)=50/q^{2}$ is accurate. We use the basis
functions $\varphi_{i}=1/q^{i}$ for $i=1,2,3$ so that

$D(q)=\frac{\theta_{1}}{q}+\frac{\theta_{2}}{q^{2}}+\frac{\theta_{3}}{q^{3}},$
(16)

and we expect to learn $\bm{\theta}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$. We
save the solution to the discrete model at $M=50$ equally spaced time points
between $t_{1}=0$ and $t_{M}=5$. With this setup, and starting with all
coefficients initially active so that $\mathcal{A}_{1}=\\{1,2,3\\}$, we obtain
the results in Table 2. The first iterate gives us $\bm{\theta}_{1}$ such that
$D(q)<0$ for some range of $q$ as we show in Figure 4(a), and so we assign
$\mathcal{L}(\bm{\theta}_{1})=\infty$. To get to the next step, we remove
$\theta_{1}$, $\theta_{2}$, and $\theta_{3}$ one a time and compute the loss
for each resulting vector, and we find that removing $\theta_{3}$ leads to a
vector that gives the least loss out of those considered. We thus find
$\mathcal{A}_{2}=\\{1,2\\}$ and
$\bm{\theta}_{2}=(-1.46,47.11,0)^{\mkern-1.5mu\mathsf{T}}$. Continuing, we
find that out of the choice of removing $\theta_{1}$ or $\theta_{2}$, or
putting $\theta_{3}$ back into the model, removing $\theta_{1}$ decreases the
loss by the greatest amount, giving $\mathcal{A}_{3}=\\{2\\}$. Finally, we
find that there are no more improvements to be made, and so the algorithm
stops at $\bm{\theta}_{3}=(0,43.52,0)^{\mkern-1.5mu\mathsf{T}}$, which is
close to the continuum limit. We emphasise that this final $\bm{\theta}_{3}$
is a least squares solution with the constraint $\theta_{1}=\theta_{3}=0$,
thus there is no need to refine $\bm{\theta}_{3}$ further by eliminating
$\theta_{1}$ and $\theta_{3}$ directly in (16), as the result would be the
same. Comparing the densities from the solution of the learned PDE with
$\bm{\theta}=\bm{\theta}_{3}$ with the discrete densities in Figure 5(a), we
see that the curves are nearly visually indistinguishable near the center, but
there are some visually discernible discrepancies near the boundaries. We show
the form of $D(q)$ at each iteration in Figure 4(a), where we observe that the
first iterate captures only the higher densities, the second iterate captures
the complete range of densities, and the third iterate removes a single term
which gives no noticeable difference.

Table 2: Stepwise equation learning results for the density data for Case Study 1: Fixed boundaries using the basis expansion (16), saving the results at $M=50$ equally spaced times between $t_{1}=0$ and $t_{M}=5$ and starting with all coefficients active, $\mathcal{A}_{1}=\\{1,2,3\\}$. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | Loss
---|---|---|---|---
1 | -5.97 | 70.73 | -27.06 | $\infty$
2 | -1.46 | 47.11 | 0.00 | -4.33
3 | 0.00 | 43.52 | 0.00 | -5.18
Figure 4: Progression of $D(q)$ over each iterate for Case Study 1: Fixed
boundaries. (a) Progression from the results in Table 2 (dashed curves). (b)
As in (a), except with the results from Table 3 using matrix pruning.

To improve our learned model we introduce matrix pruning, inspired from the
data thresholding approach in VandenHeuvel et al. [42], to improve the
estimates for $\bm{\theta}$. Visual inspection of the space-time diagram in
Figure 2(a) shows that the most significant density changes occur at early
time and near to locations where $q$ changes in the initial condition, and a
significant portion of the space-time diagram involves regions where $q$ is
almost constant. These regions where $q$ has minimal change are problematic as
points which lead to a higher residual are overshadowed, affecting the least
squares problem and consequently degrading the estimates for $\bm{\theta}$
significantly, and so it is useful to only include important points in the
construction of $\mathbf{A}$. To resolve this issue, we choose to only include
points if the associated densities falls between the $10\%$ and $90\%$
quantiles for the complete set of densities, which we refer to by density
quantiles; more details on this pruning procedure are given in Appendix D.
This choice of density quantiles is made using trial and error, starting at
$0\%$ and $100\%$, respectively, and shrinking the quantile range until
suitable results are obtained. When we apply this pruning and reconstruct
$\mathbf{A}$, we obtain the improved results in Table 3 and associated
densities in Figure 5(b). Compared with Table 2, we see that the coefficient
estimates for $\bm{\theta}$ all lead to improved losses, and our final model
now has $\bm{\theta}=(0,49.83,0)^{\mkern-1.5mu\mathsf{T}}$, which is much
closer to the the continuum limit, as we see in Figure 5(b) where the solution
curves are now visually indistinguishable everywhere. Moreover, we show in
Figure 4(b) how $D(q)$ is updated at each iteration, where we see that the
learned nonlinear diffusivity functions are barely different from the expected
continuum limit result. These results demonstrate the importance of only
including the most important points in $\mathbf{A}$.

Table 3: Improved results for Case Study 1: Fixed boundaries from Table 2, now using matrix pruning so that densities outside of the $10\%$ and $90\%$ density quantiles are not included. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}$ | $\theta_{2}$ | $\theta_{3}$ | Loss
---|---|---|---|---
1 | -1.45 | 42.48 | 13.76 | -4.19
2 | 0.00 | 37.79 | 19.69 | -5.46
3 | 0.00 | 49.83 | 0.00 | -7.97
Figure 5: Stepwise equation learning results for Case Study 1: Fixed
boundaries. (a) Comparisons of the discrete density profiles (solid curves)
with those learned from PDEs obtained from the results in Table 2 (dashed
curves). (b) As in (a), except with the results from Table 3 using matrix
pruning so that densities outside of the $10\%$ and $90\%$ density quantiles
are not included. (c) Comparisons of the learned $D(q)$ from Table 2 without
pruning, Table 3 with pruning, and the continuum limit from (8). In (a)–(b),
the arrows show the direction of increasing time, and the density profiles
shown are at times $t=0,1,2,3,4,5$ in black, red, blue, green, orange, and
purple, respectively.

###  4.2 Case Study 2: Free boundaries

Case Study 2 extends Case Study 1 by allowing the right-most cell boundary to
move so that $H(q)\neq 0$. We do not consider proliferation, giving $R(q)=0$
in (7).

The equation learning procedure for this case study is similar to Case Study
1, namely we expand $D(q)$ as in (9) and constrain $D(q)\geq 0$. In addition
to learning $D(q)$, we need to learn $H(q)$ and the evolution equation
describing the free boundary. In (7), this evolution equation is given by a
conservation statement, $q\mathrm{d}L/\mathrm{d}t=-D(q)\partial q/\partial x$
with $q=q(L(t),t)$. Here we treat this moving boundary condition more
generally by introducing a function $E(t)$ so that

$q\derivative{L}{t}=-E(q)\partialderivative{q}{x}$ (17)

at $x=L(t)$ for $t>0$. While (17) could lead to local loss of conservation at
the moving boundary, our approach is to for the possibility that coefficients
in $D(q)$ and $E(q)$ differ and to explore the extent to which this is true,
or otherwise, according to our equation learning procedure. We constrain
$E(q)\geq 0$ so that (17) makes sense for our problem and we expand $D(q)$,
$H(q)$, and $E(q)$ as follows

$D(q)=\sum_{i=1}^{d}\theta_{i}^{d}\varphi_{i}^{d}(q),\quad
H(q)=\sum_{i=1}^{h}\theta_{i}^{h}\varphi_{i}^{h}(q),\quad
E(q)=\sum_{i=1}^{h}\theta_{i}^{e}\varphi_{i}^{e}(q).$ (18)

The matrix system for
$\bm{\theta}^{d}=(\theta_{1}^{d},\ldots,\theta_{d}^{d})^{\mkern-1.5mu\mathsf{T}}$
is the same as it was in Case Study 1 in (12), which we now write as
$\mathbf{A}^{d}\bm{\theta}^{d}=\mathbf{b}^{d}$ with
$\mathbf{A}^{d}\in\mathbb{R}^{n(M-1)\times d}$ and
$\mathbf{b}^{d}\in\mathbb{R}^{n(M-1)\times 1}$ given by $\mathbf{A}$ and
$\mathbf{b}$ in (12), and we can construct two other independent matrix
systems for
$\bm{\theta}^{h}=(\theta_{1}^{h},\ldots,\theta_{h}^{h})^{\mkern-1.5mu\mathsf{T}}$
and
$\bm{\theta}^{e}=(\theta_{1}^{e},\ldots,\theta_{e}^{e})^{\mkern-1.5mu\mathsf{T}}$.
To construct these matrix systems, for a given boundary point $(x_{nj},t_{j})$
we write

$\partialderivative{q_{nj}}{x}=\sum_{k=1}^{h}\theta_{k}^{h}\varphi_{k}^{h}(q_{nj}),\quad
q_{nj}\derivative{L_{j}}{t}=-\partialderivative{q_{nj}}{x}\sum_{k=1}^{e}\theta_{k}^{e}\varphi_{k}^{e}(q_{nj}),$
(19)

where $L_{j}=x_{nj}$ is the position of the leading edge at $t=t_{j}$. In (19)
we assume that $L_{j}$ can be approximated by a smooth function so that
$\mathrm{d}L_{j}/\mathrm{d}t$ can be defined. With (19) we have
$\mathbf{A}^{h}\bm{\theta}^{h}=\mathbf{b}^{h}$ and
$\mathbf{A}^{e}\bm{\theta}^{e}=\mathbf{b}^{e}$, where

$\displaystyle\displaystyle\mathbf{A}^{h}=\begin{bmatrix}\varphi_{1}^{h}(q_{12})&\cdots&\varphi_{h}^{h}(q_{12})\\\
\vdots&\ddots&\vdots\\\
\varphi_{1}^{h}(q_{nM})&\cdots&\varphi_{h}^{h}(q_{nM})\end{bmatrix},\quad\mathbf{b}^{h}=\begin{bmatrix}\partialderivative{q_{12}}{x}\\\
\vdots\\\ \partialderivative{q_{nM}}{x}\end{bmatrix}$ (20)

with $\displaystyle\mathbf{A}^{h}\in\mathbb{R}^{(M-1)\times h}$ and
$\displaystyle\mathbf{b}^{h}\in\mathbb{R}^{(M-1)\times 1}$, and

$\mathbf{A}^{e}=\begin{bmatrix}\varphi_{1}^{e}(q_{12})\partialderivative{q_{n2}}{x}&\cdots&\varphi_{e}^{e}(q_{12})\partialderivative{q_{n2}}{x}\\\
\vdots&\ddots&\vdots\\\
\varphi_{1}^{e}(q_{nM})\partialderivative{q_{nM}}{x}&\cdots&\varphi_{e}^{e}(q_{nM})\partialderivative{q_{nM}}{x}\end{bmatrix},\quad\mathbf{b}^{e}=-\begin{bmatrix}q_{n2}\derivative{L_{2}}{t}\\\
\vdots\\\ q_{nM}\derivative{L_{M}}{t}\end{bmatrix}$ (21)

with $\mathbf{A}^{e}\in\mathbb{R}^{(M-1)\times e}$ and
$\mathbf{b}^{e}\in\mathbb{R}^{(M-1)\times 1}$. Then, writing

$\mathbf{A}=\operatorname{diag}(\mathbf{A}^{d},\mathbf{A}^{h},\mathbf{A}^{e})\in\mathbb{R}^{(n+2)(M-1)\times(d+h+e)},\quad\mathbf{b}=\begin{bmatrix}\mathbf{b}^{d}\\\
\mathbf{b}^{h}\\\ \mathbf{b}^{e}\end{bmatrix}\in\mathbb{R}^{(n+2)(M-1)\times
1},$ (22)

we obtain

$\mathbf{A}\bm{\theta}=\mathbf{b},\quad\bm{\theta}=\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{h}\\\ \bm{\theta}^{e}\end{bmatrix}\in\mathbb{R}^{(d+h+e)\times
1}.$ (23)

The solution of $\mathbf{A}\bm{\theta}=\mathbf{b}$ is the combined solution of
the individual linear systems as $\mathbf{A}$ is block diagonal. Estimates for
$\bm{\theta}^{d}$, $\bm{\theta}^{h}$, and $\bm{\theta}^{e}$ are independent,
which demonstrates the modularity of our approach, where these additional
features, in particular the leading edge, are just an extra independent
component of our procedure in addition to the procedure for estimating $D(q)$.

In addition to the new matrix system $\mathbf{A}\bm{\theta}=\mathbf{b}$ in
(23), we augment the loss function (14) to incorporate information about the
location of the moving boundary. Letting $L(t;\bm{\theta})$ denote the leading
edge from the solution of the PDE (7) with parameters $\bm{\theta}$, the loss
function is

$\displaystyle\underbrace{\mathcal{L}(\bm{\theta})}_{\textnormal{loss}}$
$\displaystyle=\overbrace{\log\left[\frac{1}{n(M-1)}\sum_{j=2}^{M}\sum_{i=1}^{n}\left(\frac{q_{ij}-q\left(x_{ij},t_{j};\bm{\theta}\right)}{q_{ij}}\right)^{2}\right]}^{\textnormal{density
goodness of fit}}$
$\displaystyle+\underbrace{\log\left[\frac{1}{M-1}\sum_{j=2}^{M}\left(\frac{L_{j}-L\left(t_{j};\bm{\theta}\right)}{L_{j}}\right)^{2}\right]}_{\textnormal{leading
edge goodness of fit}}+\underbrace{\|\bm{\theta}\|_{0}}_{\textnormal{model
complexity}}.$ (24)

Let us now apply our stepwise equation learning procedure with (23) and (24).
We consider the data from Figure 2, where we know in advance that the
continuum limit with $D(q)=50/q^{2}$, $H(q)=2q^{2}-0.4q^{3}$, and
$E(q)=50/q^{2}$ is accurate. The expansions we use for $D(q)$, $H(q)$, and
$E(q)$ are given by

$\displaystyle\begin{array}[]{rcl}D(q)&=&\dfrac{\theta_{1}^{d}}{q}+\dfrac{\theta_{2}^{d}}{q^{2}}+\dfrac{\theta_{3}^{d}}{q^{3}},\\\\[6.0pt]
H(q)&=&\theta_{1}^{h}q+\theta_{2}^{h}q^{2}+\theta_{3}^{h}q^{3}+\theta_{4}^{h}q^{4}+\theta_{5}^{h}q^{5},\\\\[6.0pt]
E(q)&=&\dfrac{\theta_{1}^{e}}{q}+\dfrac{\theta_{2}^{e}}{q^{2}}+\dfrac{\theta_{3}^{e}}{q^{3}}.\end{array}$
(28)

With these expansions, we expect to learn
$\bm{\theta}^{d}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(0,2,-0.4,0,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{e}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$. We initially consider
saving the solution at $M=1000$ equally spaced times between $t_{1}=0$ and
$t_{M}=100$, and using matrix pruning so that only points whose densities fall
within the $35\%$ and $65\%$ density quantiles are included. The results with
this configuration are shown in Table 4, where we see that we are only able to
learn $H(q)=E(q)=0$ and $D(q)=25.06/q^{3}$. This outcome highlights the
importance of choosing an appropriate time interval, since Figure 2(b)
indicates that mechanical relaxation takes place over a relative short
interval which means that working with data in $0<t\leq 100$ can lead to a
poor outcome.

Table 4: Stepwise equation learning results for Case Study 2: Free boundaries, using the basis expansions (28), saving the results at $M=1000$ equally spaced times between $t_{1}=0$ and $t_{M}=100$, pruning so that densities outside of the $35\%$ and $65\%$ density quantiles are not included, and starting with all terms inactive. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}^{d}$ | $\theta_{2}^{d}$ | $\theta_{3}^{d}$ | $\theta_{1}^{h}$ | $\theta_{2}^{h}$ | $\theta_{3}^{h}$ | $\theta_{4}^{h}$ | $\theta_{5}^{h}$ | $\theta_{1}^{e}$ | $\theta_{2}^{e}$ | $\theta_{3}^{e}$ | Loss
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -1.40
2 | 0.00 | 0.00 | 25.06 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.40

We proceed by restricting our data collection to $0\leq t\leq 15$, now saving
the solution at $M=200$ equally spaced times between $t_{1}=0$ and $t_{M}=15$.
Keeping the same quantiles for the matrix pruning, the new results are shown
in Table 5 and Figure 6. We see that the densities and leading edges are
accurate for small time, but the learned mechanisms do not extrapolate as well
for $t\geq 15$, for example $L(t)$ in Figure 6(b) does not match the discrete
data. To address this issue, we can further limit the information that we
include in our matrices, looking to only include boundary points where
$\mathrm{d}L/\mathrm{d}t$ is neither too large not too small. We implement
this by excluding all points $(x_{nj},t_{j})$ from the construction of
$(\mathbf{A}^{e},\mathbf{b}^{e})$ in (21) such that
$\mathrm{d}L_{j}/\mathrm{d}t$ is outside of the $10\%$ or $90\%$ quantiles of
the vector $(\mathrm{d}L_{2}/\mathrm{d}t,\ldots,\mathrm{d}L_{M}/\mathrm{d}t)$,
called the velocity quantiles.

Table 5: Stepwise equation learning results for Case Study 2: Free boundaries, using the basis expansions (28), saving the results at $M=200$ equally spaced times between $t_{1}=0$ and $t_{M}=15$, pruning so that densities outside of the $35\%$ and $65\%$ density quantiles are not included, and starting with all terms inactive. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}^{d}$ | $\theta_{2}^{d}$ | $\theta_{3}^{d}$ | $\theta_{1}^{h}$ | $\theta_{2}^{h}$ | $\theta_{3}^{h}$ | $\theta_{4}^{h}$ | $\theta_{5}^{h}$ | $\theta_{1}^{e}$ | $\theta_{2}^{e}$ | $\theta_{3}^{e}$ | Loss
---|---|---|---|---|---|---|---|---|---|---|---|---
1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -3.37
2 | 0.00 | 0.00 | 0.00 | 0.00 | -0.03 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -2.37
3 | 0.00 | 0.00 | 0.00 | 0.00 | -0.03 | 0.00 | 0.00 | 0.00 | 8.74 | 0.00 | 0.00 | -3.68
4 | 0.00 | 47.38 | 0.00 | 0.00 | -0.03 | 0.00 | 0.00 | 0.00 | 8.74 | 0.00 | 0.00 | -4.02
5 | 0.00 | 47.38 | 0.00 | 8.41 | -1.69 | 0.00 | 0.00 | 0.00 | 8.74 | 0.00 | 0.00 | -8.14
Figure 6: Stepwise equation learning results from Table 5 for Case Study 2:
Free boundaries. (a) Comparisons of the discrete density profiles (solid
curves) with those learned from PDEs obtained from the results in Table 5
(dashed curves), plotted at the times $t=0,5,10,25,50,100$ in black, red,
blue, green, orange, and purple, respectively. The arrow shows the direction
of increasing time. (b) As in (a), except comparing the leading edges. (c)–(e)
are comparisons of the learned forms of $D(q)$, $H(q)$, and $E(q)$ with the
forms from the continuum limit (8).

Implementing thresholding on $\mathrm{d}L/\mathrm{d}t$ leads to the results
presented in Figure 7. We see that the learned densities and leading edges are
both visually indistinguishable from the discrete data. Since $H(q)$ and
$E(q)$ are only ever evaluated at $x=L(t)$, and $q(L(t),t)\approx 5$ for
$t>0$, we see that $H(q)$ and $E(q)$ only match the continuum limit at
$q\approx 5$, which means that our learned continuum limit model conserves
mass and is consistent with the traditional coarse-grained continuum limit, as
expected. We discuss in Appendix E how we can enforce $D(q)=E(q)$ to guarantee
conservation mass from the outset, however our approach in Figure 7 is more
general in the sense that our learned continuum limit is obtained without
making any a priori assumptions about the form of $E(q)$.

Figure 7: Stepwise equation learning results from Table 5 for Case Study 2:
Free boundaries, except also using matrix pruning on
$(\mathbf{A}_{3},\mathbf{b}_{3})$ so points where
$\mathrm{d}L_{j}/\mathrm{d}t$ falls outside of the $10\%$ and $90\%$ velocity
quantiles are excluded, giving $\theta_{1}^{e}=9.42$ rather than $8.74$. (a)
Comparisons of the discrete density profiles (solid curves) with those from
the learned PDE (dashed curves), plotted at the times $t=0,5,10,25,50,100$ in
black, red, blue, green, orange, and purple, respectively. The arrow shows the
direction of increasing time. (b) As in (a), except comparing the leading
edges. (c)–(e) are comparisons of the learned forms of $D(q)$, $H(q)$, and
$E(q)$ with the forms from the continuum limit (8).

###  4.3 Case Study 3: Fixed boundaries with proliferation

Case Study 3 is identical to Case Study 1 except that we incorporate cell
proliferation, implying $R(q)\neq 0$ in (7). This case is more complicated
than with mechanical relaxation only, as we have to consider how we combine
the repeated realisations to capture the average density data as well. For
this work, we average over each realisation at each time using linear
interpolants as described in Appendix D. This averaging procedure gives
$n_{k}$ points $\bar{x}_{ij}$ between $x=0$ and $x=30$ at each time $t_{j}$,
$j=1,\ldots,M$, with corresponding density value $\bar{q}_{ij}$. The
quantities $\bar{x}_{ij}$ and $\bar{q}_{ij}$ play the same role as $x_{ij}$
and $q_{ij}$ in the previous case studies.

To apply equation learning we note there is no moving boundary, giving
$H(q)=0$ in (8). We proceed by expanding $D(q)$ and $R(q)$ as follows

$D(q)=\sum_{i=1}^{d}\theta_{i}^{d}\varphi_{i}^{d}(q),\quad
R(q)=\sum_{i=1}^{r}\theta_{i}^{r}\varphi_{i}^{r}(q),$ (29)

with the aim of estimating
$\bm{\theta}^{d}=(\theta_{1}^{d},\ldots,\theta_{d}^{d})^{\mkern-1.5mu\mathsf{T}}$
and
$\bm{\theta}^{r}=(\theta_{1}^{r},\ldots,\theta_{r}^{r})^{\mkern-1.5mu\mathsf{T}}$,
again constraining $D(q)\geq 0$. We expand the PDE from (10), as in Section
4(a), and the only difference is the additional term
$\sum_{m=1}^{r}\varphi_{m}^{r}(\bar{q}_{ij})\theta_{m}^{r}$ for each point
$(\bar{x}_{ij},t_{j})$. Thus, we have the same matrix as in Section 4(a),
denoted $\mathbf{A}^{d}\in\mathbb{R}^{n_{k}(M-1)\times d}$, and a new matrix
$\mathbf{A}^{r}\in\mathbb{R}^{n_{k}(M-1)\times r}$ whose row corresponding to
the point $(\bar{x}_{ij},t_{j})$ is given by

$\displaystyle\mathbf{a}_{ij}^{r}=\begin{bmatrix}\varphi_{1}^{r}(\bar{q}_{ij})&\cdots&\varphi_{r}^{r}(\bar{q}_{ij})\end{bmatrix}\in\mathbb{R}^{1\times
r},$ (30)

so that the coefficient matrix $\mathbf{A}$ is now

$\mathbf{A}=\begin{bmatrix}\mathbf{A}^{d}&\mathbf{A}^{r}\end{bmatrix}\in\mathbb{R}^{n_{k}(M-1)\times(d+r)}.$
(31)

The corresponding entry for the point $(\bar{x}_{ij},t_{j})$ in
$\mathbf{b}\in\mathbb{R}^{n_{k}(M-1)\times 1}$ is
$\partial\bar{q}_{ij}/\partial t$. Notice that this additional term in the PDE
adds an extra block to the matrix without requiring a significant coupling
with the existing equations from the simpler problem without proliferation.
Thus, we estimate our coefficient vectors using the system

$\mathbf{A}\bm{\theta}=\mathbf{b},\quad\bm{\theta}=\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{r}\end{bmatrix}\in\mathbb{R}^{(d+r)\times 1}.$ (32)

We can take exactly the same stepwise procedure as in Section 4(a), except now
the loss function (14) uses $n_{k}$, $\bar{q}_{ij}$, and $\bar{x}_{ij}$ rather
than $n$, $q_{ij}$, and $x_{ij}$, respectively.

####  4.3.1 Accurate continuum limit

Let us now apply these ideas to our data from Figure 2, where we know that the
continuum limit with $D(q)=50/q^{2}$ and $R(q)=0.15q-0.01q^{2}$ is accurate.
The expansions we use for $D(q)$ and $R(q)$ are given by

$D(q)=\frac{\theta_{1}^{d}}{q}+\frac{\theta_{2}^{d}}{q^{2}}+\frac{\theta_{3}^{d}}{q^{3}},\quad
R(q)=\theta_{1}^{r}q+\theta_{2}^{r}q^{2}+\theta_{3}^{r}q^{3}+\theta_{4}^{r}q^{4}+\theta_{5}^{r}q^{5},$
(33)

and we expect to learn $\bm{\theta}^{d}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$ and
$\bm{\theta}^{r}=(0.15,-0.01,0,0,0)^{\mkern-1.5mu\mathsf{T}}$. We average over
$1000$ identically-prepared realisations, saving the solutions at $M=501$
equally spaced times between $t_{1}=0$ and $t_{M}=50$ with $n_{k}=50$ knots
for averaging. For this problem, and for Case Study 4 discussed later, we find
that working with $1000$ identically-prepared realisations of the stochastic
models leads to sufficiently smooth density profiles. As discussed in Appendix
F, the precise number of identically-prepared realisations is not important
provided that the number is sufficiently large; when not enough realisations
are taken, the results are inconsistent across different sets of realisations
and will fail to identify the average behaviour from the learned model. We
also use matrix pruning so that we only include points whose densities fall
within the $10\%$ and $90\%$ density quantiles, as done in Section 4(a). The
results we obtain are shown in Table 6, starting with all coefficients active.

Table 6: Stepwise equation learning results for Case Study 3: Fixed boundaries
with proliferation, where the continuum limit is accurate, using the basis
expansions (33), saving the results at $M=501$ equally spaced times between
$t_{1}=0$ and $t_{M}=50$, averaging across $1000$ realisations with $n_{k}=50$
knots, pruning so that densities outside of the $10\%$ and $90\%$ density
quantiles are not included, and starting with all diffusion and reaction
coefficients active. Coefficients highlighted in blue show the coefficient
chosen to be removed or added at the corresponding step.

Step | $\theta_{1}^{d}$ | $\theta_{2}^{d}$ | $\theta_{3}^{d}$ | $\theta_{1}^{r}$ | $\theta_{2}^{r}$ |  $\theta_{3}^{r}$ ($\times 10^{-4}$) |  $\theta_{4}^{r}$ ($\times 10^{-5}$) |  $\theta_{5}^{r}$ ($\times 10^{-7}$) | Loss
---|---|---|---|---|---|---|---|---|---
1 | -11.66 | 147.43 | -191.51 | 0.13 | -0.00 | -0.00 | $5.83$ | $-11.30$ | $\infty$
2 | -2.24 | 60.86 | 0.00 | 0.13 | -0.00 | $-5.72$ | $2.62$ | $-3.49$ | -0.71
3 | -2.25 | 60.90 | 0.00 | 0.14 | -0.01 | 0.00 | $-1.25$ | $5.95$ | -1.92
4 | 0.00 | 52.95 | 0.00 | 0.14 | -0.01 | 0.00 | $-1.36$ | $6.49$ | -3.35
5 | 0.00 | 53.02 | 0.00 | 0.15 | -0.01 | 0.00 | 0.00 | $0.32$ | -4.98
6 | 0.00 | 52.97 | 0.00 | 0.15 | -0.01 | 0.00 | 0.00 | 0.00 | -5.70

Figure 8: Stepwise equation learning results for Case Study 3: Fixed
boundaries with proliferation, where the continuum limit is accurate. (a)
Comparisons of the discrete density profiles (solid curves) with those learned
from PDEs obtained from the results in Table 6 (dashed curves), plotted at the
times $t=0,1,5,10,20,50$ in black, red, blue, green, orange, and purple,
respectively. The arrow shows the direction of increasing time. (b)–(c) are
comparisons of $D(q)$ and $R(q)$ with the forms from the continuum limit (8).

Table 6 shows that we find
$\bm{\theta}^{d}=(0,52.97,0)^{\mkern-1.5mu\mathsf{T}}$ and
$\bm{\theta}^{r}=(0.15,-0.010,0,0,0)^{\mkern-1.5mu\mathsf{T}}$, which are both
very close to the continuum limit. Figure 8 visualises these results, showing
that the PDE solutions with the learned $D(q)$ and $R(q)$ match the discrete
densities, and the mechanisms that we do learn are visually indistinguishable
with the continuum limit functions (8) as shown in Figure 8(b)–(c).

####  4.3.2 Inaccurate continuum limit

We now extend the problem so that the continuum limit is no longer accurate,
taking $k=1/5$ to be consistent with Figure 3(a). Using the same basis
expansions in (33), we save the solution at $M=751$ equally spaced times
between $t_{1}=0$ and $t_{M}=75$, averaging over $1000$ realisations with
$n_{k}=200$. We find that we need to use the $25\%$ and $75\%$ density
quantiles rather than the $10\%$ and $90\%$ density quantiles, as in the
previous example, to obtain results in this case. With this configuration, the
results we find are shown in Table 7 and Figure 9.

Table 7: Stepwise equation learning results for Case Study 3: Fixed boundaries with proliferation, where the continuum limit is inaccurate, using the basis expansions (33), saving the results at $M=751$ equally spaced times between $t_{1}=0$ and $t_{M}=75$, averaging across $1000$ realisations with $n_{k}=200$ knots, pruning so that densities outside of the $25\%$ and $75\%$ density quantiles are not included, and starting with all diffusion and reaction coefficients inactive. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}^{d}$ | $\theta_{2}^{d}$ | $\theta_{3}^{d}$ | $\theta_{1}^{r}$ | $\theta_{2}^{r}$ |  $\theta_{3}^{r}$ ($\times 10^{-4}$) |  $\theta_{4}^{r}$ ($\times 10^{-5}$) | $\theta_{5}^{r}$ | Loss
---|---|---|---|---|---|---|---|---|---
1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -0.33
2 | 0.00 | 0.00 | 0.00 | 0.02 | 0.00 | 0.00 | 0.00 | 0.00 | 0.51
3 | 0.00 | 0.00 | 0.00 | 0.11 | -0.01 | 0.00 | 0.00 | 0.00 | 0.20
4 | 0.00 | 0.11 | 0.00 | 0.11 | -0.01 | 0.00 | 0.00 | 0.00 | -0.04
5 | 0.00 | 0.12 | 0.00 | 0.13 | -0.01 | $1.59$ | 0.00 | 0.00 | -0.46
6 | 0.00 | 0.12 | 0.00 | 0.16 | -0.02 | $7.49$ | $-1.69$ | 0.00 | -1.13
Figure 9: Stepwise equation learning results for Case Study 3: Fixed
boundaries with proliferation, where the continuum limit is inaccurate. (a)
Comparisons of the discrete density profiles (solid curves) with those learned
from PDEs obtained from the results in Table 6 (dashed curves), plotted at the
times $t=0,1,10,25,40,75$ in black, red, blue, green, orange, and purple,
respectively. The arrow shows the direction of increasing time. (b)–(c) are
comparisons of $D(q)$ and $R(q)$ with the forms from the continuum limit (8).

Results in Table 7 show $\bm{\theta}^{d}=(0,0.12,0)^{\mkern-1.5mu\mathsf{T}}$,
which is reasonably close to the continuum limit with
$(0,0.2,0)^{\mkern-1.5mu\mathsf{T}}$. The reaction vector, for which the
continuum limit is $(0.15,-0.01,0,0,0)^{\mkern-1.5mu\mathsf{T}}$ so that
$R(q)$ is a quadratic, is now given by $\bm{\theta}^{r}=(0.16,-0.02,7.49\times
10^{-4},-1.69\times 10^{-5},0)^{\mkern-1.5mu\mathsf{T}}$, meaning the learned
$R(q)$ is a quartic. Figure 9 compares the averaged discrete densities with
the solution of the learned continuum limit model. Figure 9(c) compares the
learned source term with the continuum limit. While both terms are visually
indistinguishable at small densities, we see that the two source terms differ
at high densities, with the learned carrying capacity density, where $R(q)=0$,
reduced relative to the continuum limit. This is consistent with previous
results [15].

###  4.4 Case Study 4: Free boundaries with proliferation

Case Study 4 is identical to Case Study 2 except that we now introduce
proliferation into the discrete model so that $R(q)\neq 0$ in (7). First, as
in Case Study 3 and as described in Appendix D, we average our data across
each realisation from our discrete model. This averaging provides us with
points $\bar{x}_{ij}$ between $x=0$ and $x=\bar{L}_{j}$ at each time $t_{j}$,
$j=1,\ldots,M$, where $\bar{L}_{j}$ is the average leading edge at $t=t_{j}$,
with corresponding density values $\bar{q}_{ij}$, where $i=1,\ldots,n_{k}$ and
$n_{k}$ is the number of knots to use for averaging. We expand the functions
$D(q)$, $R(q)$, $H(q)$, and $E(q)$ as

$D(q)=\sum_{i=1}^{d}\theta_{i}^{d}\varphi_{i}^{d}(q),~{}R(q)=\sum_{i=1}^{r}\theta_{i}^{r}\varphi_{i}^{r}(q),~{}H(q)=\sum_{i=1}^{h}\theta_{i}^{h}\varphi_{i}^{h}(q),~{}E(q)=\sum_{i=1}^{e}\theta_{i}^{e}\varphi_{i}^{e}(q),$
(34)

again restricting $D(q),E(q)\geq 0$. The function $E(q)$ is used in the moving
boundary condition in (7), as in (17). The matrix $\mathbf{A}$ and vector
$\mathbf{b}$ are given by

$\displaystyle\mathbf{A}=\operatorname{diag}(\mathbf{A}^{dr},\mathbf{A}^{h},\mathbf{A}^{e})\in\mathbb{R}^{n_{k}(M-1)\times(d+r+h+e)},\quad\mathbf{b}=\begin{bmatrix}\mathbf{b}^{dr}\\\
\mathbf{b}^{h}\\\ \mathbf{b}^{e}\end{bmatrix}\in\mathbb{R}^{n_{k}(M-1)},$ (35)

where $\mathbf{A}^{dr}=[\mathbf{A}^{d}~{}\mathbf{A}^{r}]$ as defined in (31),
$\mathbf{A}^{h}$ and $\mathbf{A}^{e}$ are the matrices from (20) and (21),
respectively, and similarly for $\mathbf{b}^{dr}=\partial\mathbf{q}/\partial
t$, $\mathbf{b}^{h}$, and $\mathbf{b}^{e}$ from (12), (20), and (21),
respectively. Thus,

$\mathbf{A}\bm{\theta}=\mathbf{b},\quad\bm{\theta}=\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{r}\\\ \bm{\theta}^{h}\\\
\bm{\theta}^{e}\end{bmatrix}\in\mathbb{R}^{(d+r+h+e)\times 1}.$ (36)

Similar to Case Study 2, the coefficients for each mechanism are independent,
except for $\bm{\theta}^{d}$ and $\bm{\theta}^{r}$. The loss function we use
is the loss function from (24).

With this problem, it is difficult to learn all mechanisms simultaneously,
especially as mechanical relaxation and proliferation occur on different time
scales since mechanical relaxation dominates in the early part of the
simulation, whereas both proliferation and mechanical relaxation play a role
at later times. This means $D(q)$ and $R(q)$ cannot be estimated over the
entire time range as was done in Case Study 3. To address this we take a
sequential learning procedure to learn these four mechanisms using four
distinct time intervals $I^{d}$, $I^{e}$, $I^{h}$, and $I^{r}$:

  1. 1.

Fix $R(q)=H(q)=E(q)=0$ and learn $\bm{\theta}^{d}$ over $t\in I^{d}$, solving
$\mathbf{A}^{d}\bm{\theta}^{d}=\mathbf{b}^{dr}$.

  2. 2.

Fix $R(q)=H(q)=0$ and $\bm{\theta}^{d}$ and learn $\bm{\theta}^{e}$ over $t\in
I^{e}$, solving $\mathbf{A}^{e}\bm{\theta}^{e}=\mathbf{b}^{e}$.

  3. 3.

Fix $R(q)=0$, $\bm{\theta}^{d}$, and $\bm{\theta}^{e}$ and learn
$\bm{\theta}^{h}$ over $t\in I^{h}$, solving
$\mathbf{A}^{h}\bm{\theta}^{h}=\mathbf{b}^{h}$.

  4. 4.

Fix $\bm{\theta}^{d}$, $\bm{\theta}^{e}$, and $\bm{\theta}^{h}$ and learn
$\bm{\theta}^{r}$ over $t\in I^{r}$, solving
$\mathbf{A}^{r}\bm{\theta}^{r}=\mathbf{b}^{dr}-\mathbf{A}^{d}\bm{\theta}^{d}$.

In these steps, solving the system $\mathbf{A}\bm{\theta}=\mathbf{b}$ means to
apply our stepwise procedure to this system; for these problems, we start each
procedure with no active coefficients. The modularity of our approach makes
this sequential learning approach straightforward to implement. For these
steps, the interval $I^{d}$ must be over sufficiently small times so that
proliferation does not dominate, noting that fixing $R(q)=0$ will not allow us
to identify any proliferation effects when estimating the parameters. This is
less relevant for $I^{h}$ and $I^{e}$ as the estimates of $\bm{\theta}^{h}$
and $\bm{\theta}^{e}$ impact the moving boundary only.

###  4.5 Accurate continuum limit

We apply this procedure to data from Figure 2, where the continuum limit is
accurate with $D(q)=50/q^{2}$, $R(q)=0.15q-0.01q^{2}$, $H(q)=2q^{2}-0.4q^{3}$,
and $E(q)=50/q^{2}$. The expansions we use are

$\displaystyle\displaystyle\begin{array}[]{rcll}D(q)&=&\frac{\theta_{1}^{d}}{q}+\frac{\theta_{2}^{d}}{q^{2}}+\frac{\theta_{3}^{d}}{q^{3}},\\\
R(q)&=&\theta_{1}^{r}q+\theta_{2}^{r}q^{2}+\theta_{3}^{r}q^{3}+\theta_{4}^{r}q^{4}+\theta_{5}^{r}q^{5},\\\
H(q)&=&\theta_{1}^{h}q+\theta_{2}^{h}q^{2}+\theta_{3}^{h}q^{3}+\theta_{4}^{h}q^{4}+\theta_{5}^{h}q^{5},\\\
E(q)&=&\frac{\theta_{1}^{e}}{q}+\frac{\theta_{2}^{e}}{q^{2}}+\frac{\theta_{3}^{e}}{q^{3}}.\end{array}$
(41)

With these expansions, we expect to learn
$\bm{\theta}^{d}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{r}=(0.15,-0.01,0,0,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(0,2,-0.4,0,0,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{e}=(0,50,0)^{\mkern-1.5mu\mathsf{T}}$. We average the data over
$1000$ realisations. For saving the solution, the time intervals we use are
$I^{d}=[0,0.1]$, $I^{e}=[0,5]$, $I^{h}=[5,10]$, and $I^{r}=[10,50]$, with
$25$, $50$, $100$, and $250$ time points inside each time interval for saving.
For interpolating the solution to obtain the averages, we use $n_{k}=25$,
$n_{k}=50$, $n_{k}=100$, and $n_{k}=50$ over $I^{d}$, $I^{e}$, $I^{h}$, and
$I^{r}$, respectively.

To now learn the mechanisms, we apply the sequential procedure described for
learning them one at a time. For each problem, we apply pruning so that points
outside of the $10\%$ and $90\%$ density quantiles or the $20\%$ and $80\%$
velocity quantiles are not included. We find that
$\bm{\theta}^{d}=(0,49.60,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{e}=(0,49.70,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(-0.0084,0,0,-0.0011,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{r}=(0.15,-0.010,0,0,0)^{\mkern-1.5mu\mathsf{T}}$. The results
with all these learned mechanisms are shown in Figure 10. We see from the
comparisons in Figure 10(a)–(b) that the PDE results from the learned
mechanisms are nearly indistinguishable from the discrete densities. Similar
to Case Study 2, $H(q)$ only matches the continuum limit at $q(L(t),t)$. Note
also that the solutions in Figure 10(a) go up to $t=100$, despite the stepwise
procedure considering only times up to $t=50$.

Figure 10: Stepwise equation learning results for Case Study 4: Free
boundaries with proliferation, when the continuum limit is accurate, using the
learned mechanisms with
$\bm{\theta}^{d}=(0,49.60,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{e}=(0,49.70,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(-0.0084,0,0,-0.0011,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{r}=(0.15,-0.010,0,0,0)^{\mkern-1.5mu\mathsf{T}}$. (a)
Comparisons of the discrete density profiles (solid curves) with those learned
from PDEs with the given $\bm{\theta}^{d}$, $\bm{\theta}^{e}$,
$\bm{\theta}^{h}$, and $\bm{\theta}^{r}$ (dashed curves), plotted at the times
$t=0,5,10,25,50,100$ in black, red, blue, green, orange, and purple,
respectively. The arrow shows the direction of increasing time. (b) As in (a),
except comparing the leading edges. (c)–(f) are comparisons of the learned
forms of $D(q)$, $R(q)$, $H(q)$, and $E(q)$ with the forms from the continuum
limit (8).

###  4.6 Inaccurate continuum limit

We now consider data from Figure 3(b) where the continuum limit is inaccurate.
Here, $k=1/5$ and the continuum limit vectors are
$\bm{\theta}^{d}=(0,0.2,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{r}=(0.15,-0.01,0,0,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(0,2,-0.4,0,0,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{e}=(0,0.2,0)^{\mkern-1.5mu\mathsf{T}}$. Using the same
procedures and expansions as Figure 10, we average the data over $1000$
realisations. The time intervals we use are $I^{d}=[0,2]$, $I^{e}=[2,10]$,
$I^{h}=[10,20]$, and $I^{r}=[20,50]$, using $20$ time points for $I^{d}$ and
$200$ time points for $I^{e}$, $I^{h}$, and $I^{r}$. We use $n_{k}=50$ knots
for averaging the solution over $I^{d}$, and $n_{k}=100$ knots for averaging
the solution over $I^{e}$, $I^{h}$, and $I^{r}$.

To apply the equation learning procedure we prune all matrices so that points
outside of the $40\%$ and $60\%$ temporal quantiles are eliminated, where the
temporal quantiles are the quantiles of $\partial q/\partial t$ from the
averaged discrete data, and similarly for points outside of the $40\%$ and
$60\%$ velocity quantiles. We find
$\bm{\theta}^{d}=(0,0.21,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{e}=(0,0.23,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(-0.15,0,0,-0.0079,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{r}=(0.11,-0.0067,0,0,0)^{\mkern-1.5mu\mathsf{T}}$.
Interestingly, here we learn $R(q)$ is quadratic with coefficients that differ
from the continuum limit. The results with all these learned mechanisms are
shown in Figure 11. We see from the comparisons in Figure 11 that the PDE
results from the learned mechanisms are visually indistinguishable from the
discrete densities. Moreover, as in Figure 10, the learned $H(q)$ and $E(q)$
match the continuum results at $q(L(t),t)$ which confirms that the learned
continuum limit conserves mass, as expected. Note also that the solutions in
Figure 11(a) go up to $t=250$, despite the stepwise procedure considering only
times up to $t=50$, demonstrating the extrapolation power of our method.

Figure 11: Stepwise equation learning results for Case Study 4: Free
boundaries with proliferation, when the continuum limit is inaccurate, using
the learned mechanisms with
$\bm{\theta}^{d}=(0,0.21,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{e}=(0,0.23,0)^{\mkern-1.5mu\mathsf{T}}$,
$\bm{\theta}^{h}=(-0.15,0,0,-0.0079,0)^{\mkern-1.5mu\mathsf{T}}$, and
$\bm{\theta}^{r}=(0.11,-0.0067,0,0,0)^{\mkern-1.5mu\mathsf{T}}$. (a)
Comparisons of the discrete density profiles (solid curves) with those learned
from PDEs with the given $\bm{\theta}^{d}$, $\bm{\theta}^{e}$,
$\bm{\theta}^{h}$, and $\bm{\theta}^{r}$ (dashed curves), plotted at the times
$t=0,5,25,50,100,250$ in black, red, blue, green, orange, and purple,
respectively. The arrow shows the direction of increasing time. (b) As in (a),
except comparing the leading edges. (c)–(f) are comparisons of the learned
forms of $D(q)$, $R(q)$, $H(q)$, and $E(q)$ with the forms from the continuum
limit (8).

##  5 Conclusion and discussion

In this work, we presented a stepwise equation learning framework for learning
continuum descriptions of discrete models describing population biology
phenomena. Our approach provides accurate continuum approximations when
standard coarse-grained approximations are inaccurate. The framework is simple
to implement, efficient, easily parallelisable, and modular, allowing for
additional components to be added into a model with minimal changes required
to accommodate them into an existing procedure. In contrast to other
approaches, like neural networks [41] or linear regression approaches [43],
results from our procedure are interpretable in terms of the underlying
discrete process. The coefficients incorporated or removed at each stage of
our procedure give a sense of the influence each model term contributes to the
model, giving a greater interpretation of the results, highlighting an
advantage of the stepwise approach over traditional sparse regression methods
[37, 32, 38]. The learned continuum descriptions from our procedure enable the
discovery of new mechanisms and equations describing the data from the
discrete model. For example, the discovered form of $D(q)$ can be interpreted
relative to the discrete model, describing the interaction forces between
neighbouring cells. In addition, we found in Case Study 4 that, when $k=1/5$
so that the continuum limit is inaccurate, the positive root of the quadratic
form of the source term $R(q)$ is greater than the mean field carrying
capacity density $K$, as seen in Figure 11. This increase suggests that, when
the rate of mechanical relaxation is small relative to the proliferation
rates, the mean field carrying capacity density in the continuum description
can be different from that in the discrete model.

We demonstrated our approach using a series of four biologically-motivated
case studies that incrementally build on each other, studying a discrete
individual-based mechanical free boundary model of epithelial cells [16, 15,
10, 11]. In the first two case studies, we demonstrated that we can easily
rediscover the continuum limit models derived by Baker et al. [16], including
the equations describing the evolution of the free boundary. The last two case
studies demonstrate that, when the coarse-grained models are inaccurate, our
approach can learn an accurate continuum approximation. The last case study
was the most complicated, with four mechanisms needing to be learned, but the
modularity of our approach made it simple to apply a sequential procedure to
learning the mechanisms, applying the procedure to each mechanism in sequence.
Our procedure was able to recover terms that conserved mass, despite not
enforcing conservation of mass explicitly. The procedure as we have described
does have some limitations, such as assuming that the mechanisms are linear
combinations of basis functions, which could be handled more generally by
instead using nonlinear least squares [42]. The procedure may also be
sensitive to the quality of the data points included in the matrices, and thus
to the parameters used for the procedure. In Appendix F, we discuss a
parameter sensitivity study that investigates this in greater detail. In this
parameter sensitivity study, we find that the most important parameters to
choose are the pruning parameters. These parameters can be easily tuned thanks
to the efficiency of our method, modifying each parameter in sequence and
using trial and error to determine suitable parameter values.

There are many avenues for future work based on our approach. Firstly, two-
dimensional extensions of our discrete model could be considered [55, 56],
which would follow the same approach except the continuum problems would have
to be solved using a more detailed numerical approximation [57, 58, 59].
Another avenue for exploration would be to consider applying the discrete
model on a curved interface which is more realistic than considering an
epithelial sheet on a flat substrate [60, 61]. Working with heterogeneous
populations of cells, where parameters in the discrete model can vary between
individuals in the population, is also another interesting option for future
exploration [14]. Uncertainty quantification could also be considered using
bootstrapping [42] or Bayesian inference [62]. Allowing for uncertainty
quantification would also allow for noisy data sets to be modelled, unlike the
idealised, noise-free data used in this work. We emphasise that, regardless of
the approach taken for future work, we believe that our flexible stepwise
learning framework can form the basis of these potential future studies.

## References

  * [1] Maclaren OJ, Byrne HM, Fletcher AG, Maini PK. 2015 Models, measurement and inference in epithelial tissue dynamics. Mathematical Biosciences and Engineering 12, 1321–1340.
  * [2] Turner S, Sherratt JA, Painter KJ, Savill NJ. 2004 From a discrete to a continuous model of biological cell movement. Physical Review E 69, 021910.
  * [3] Alber M, Chen N, Lushnikov PM, Newman SA. 2007 Continuous macroscopic limit of a discrete stochastic model for interaction of living cells. Physical Review Letters 99, 168102.
  * [4] Zmurchok C, Bhaskar D, Edelstein-Keshet L. 2018 Coupling mechanical tension and GTPase signaling to generate cell and tissue dynamics. Physical Biology 15, 046004.
  * [5] Marée AFM, Grieneisen VA, Edelstein-Keshet L. 2012 How cells integrate complex stimuli: The effect of feedback from phosphoinositides and cell shape on cell polarization and motility. PLoS Computational Biology 8, e1002402.
  * [6] Mason J, Jack RL, Bruna M. 2023 Macroscopic behaviour in a two-species exclusion process via the method of matched asymptotics. Journal of Statistical Physics 190.
  * [7] Bruna M, Chapman SJ, Schmidtchen M. 2023 Derivation of a macroscopic model for Brownian hard needles. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 479, 20230076.
  * [8] Bruna M, Chapman SJ. 2012a Diffusion of multiple species with excluded-volume effects. The Journal of Chemical Physics 137, 204116.
  * [9] Bruna M, Chapman SJ. 2012b Excluded-volume effects in the diffusion of hard spheres. Physical Review E 85, 011103.
  * [10] Murray PJ, Edwards CM, Tindall MJ, Maini PK. 2009 From a discrete to a continuum model of cell dynamics in one dimension. Physical Review E 80, 031912.
  * [11] Murray PJ, Edwards CM, Tindall MJ, Maini PK. 2012 Classifying general nonlinear force laws in cell-based models via the continuum limit. Physical Review E 85, 021921.
  * [12] Tambyah TA, Murphy RJ, Buenzli PR, Simpson MJ. 2021 A free boundary mechanobiological model of epithelial tissues. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20200528\.
  * [13] Murphy RJ, Buenzli PR, Tambyah TA, Thompson EW, Hugo HJ, Baker RE, Simpson MJ. 2021 The role of mechanical interactions in EMT. Physical Biology 18, 046001.
  * [14] Murphy RJ, Buenzli PR, Baker RE, Simpson MJ. 2019 A one-dimensional individual-based mechanical model of cell movement in heterogeneous tissues and its coarse-grained approximation. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 475, 20180838.
  * [15] Murphy RJ, Buenzli PR, Baker RE, Simpson MJ. 2020 Mechanical cell competition in heterogeneous epithelial tissues. Bulletin of Mathematical Biology 82, 130.
  * [16] Baker RE, Parker A, Simpson MJ. 2019 A free boundary model of epithelial dynamics. Journal of Theoretical Biology 481, 61–74.
  * [17] Fozard JA, Byrne HM, Jensen OE, King JR. 2010 Continuum approximations of individual-based models for epithelial monolayers. Mathematical Medicine and Biology: A Journal of the IMA 27, 39–74.
  * [18] Supekar R, Song B, Hastewell A, Choi GPT, Mietke A, Dunkel J. 2023 Learning hydrodynamic equations for active matter from particle simulations and experiments. Proceedings of the National Academy of Sciences of the United States of America 120, e2206994120.
  * [19] Español P. 2004 Statistical mechanics of coarse-graining. Novel Methods in Soft Matter Computing 640, 69–115.
  * [20] Middleton AM, Fleck C, Grima R. 2014 A continuum approximation to an off-lattice individual-cell based model of cell migration and adhesion. Journal of Theoretical Biology 359, 220–232.
  * [21] Jeon J, Quaranta V, Cummings PT. 2010 An off-lattice hybrid discrete-continuum model of tumor growth and invasion. Biophysical Journal 98, 37–47.
  * [22] Osborne JM, Walker A, Kershaw SK, Mirams GR, Fletcher AG, Pathmanathan P, Gavaghan D, Jensen OE, Maini PK, Byrne HM. 2010 A hybrid approach to multi-scale modelling of cancer. Philosophical Transactions of the Royal Society A 368, 5013–5028.
  * [23] Buttenschön A, Edelstein-Keshet L. 2020 Bridging from single to collective cell migration: A review of models and links to experiments. PLoS Computational Biology 16, e1008411.
  * [24] Surendran A, Plank MJ, Simpson MJ. 2020 Population dynamics with spatial structure and an Allee effect. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20200501.
  * [25] Lorenzi T, Murray PJ, Ptashnyk M. 2020 From individual-based mechanical models of multicellular systems to free-boundary problems. Interfaces and Free Boundaries: Mathematical Analysis, Computation and Applications 22, 205–244.
  * [26] Romeo N, Hastewell A, Mietke A, Dunkel J. 2021 Learning developmental mode dynamics from single-cell trajectories. eLife 10, e68679.
  * [27] Pughe-Sanford JL, Quinn S, Balabanski T, Grigoriev RO. 2023 Computing chaotic time-averages from a small number of periodic orbits. (https://arxiv.org/abs/2307.09626v1).
  * [28] Vo BN, Drovandi CC, Pettitt AN, Simpson MJ. 2015 Quantifying uncertainty in parameter estimates for stochastic models of collective cell spreading using approximate Bayesian computation. Mathematical Biosciences 263, 133–142.
  * [29] Hufton PG, Lin YT, Galla T. 2019 Model reduction methods for population dynamics with fast-switching environments: Reduced master equations, stochastic differential equations, and applications. Physical Review E 99, 032122.
  * [30] Ihle T. 2016 Chapman-Enskog expansion for the Vicsek model of self-propelled particles. Journal of Statistical Mechanics: Theory and Experiment p. 083205.
  * [31] Codling EA, Plank MJ, Benhamou S. 2008 Random walk models in biology. Journal of the Royal Society Interface 5, 813–834.
  * [32] Nardini JT, Baker RE, Simpson MJ, Flores KB. 2021 Learning differential equation models from stochastic agent-based model simulations. Journal of the Royal Society Interface 18, 20200987.
  * [33] Simpson MJ, Landman KA, Hughes BD. 2010 Cell invasion with proliferation mechanisms motivated by time-lapse data. Physica A: Statistical Mechanics and its Applications 389, 3779–3790.
  * [34] Hughes BD. 1995 Random walks and random environments: Random walks. Oxford: Clarendon Press.
  * [35] Chopard B, Droz M. 1998 Cellular automata modeling of physical systems. Cambridge: Cambridge University Press.
  * [36] Evans DJ, Morriss G. 2008 Statistical mechanics of nonequilibrium liquids. Cambridge: Cambridge University Press.
  * [37] Rudy SH, Brunton SL, Proctor JL, Kutz JN. 2017 Data-driven discovery of partial differential equations. Science Advances 3, e1602614.
  * [38] Brunton SL, Proctor JL, Kutz JN. 2016 Discovering governing equations from data by sparse identification of nonlinear dynamical systems. Proceedings of the National Academy of Sciences of the United States of America 113, 3932–3937.
  * [39] Nardini JT, Lagergren JH, Hawkins-Daarud A, Curtin L, Morris B, Rutter EM, Swanson KR, Flores KB. 2020 Learning equations from biological data with limited time samples. Bulletin of Mathematical Biology 82, 119\.
  * [40] Lagergren JH, Nardini JT, Lavigne GM, Rutter EM, Flores KB. 2020a Learning partial differential equations for biological transport models from noisy spatio-temporal data. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 476, 20190800.
  * [41] Lagergren JH, Nardini JT, Baker RE, Simpson MJ, Flores KB. 2020b Biologically-informed neural networks guide mechanistic modelling from sparse experimental data. PLoS Computational Biology 16, e1008462.
  * [42] VandenHeuvel DJ, Drovandi C, Simpson MJ. 2022 Computationally efficient mechanism discovery for cell invasion with uncertainty quantification. PLoS Computational Biology 18, e1010599.
  * [43] Simpson MJ, Baker RE, Buenzli PR, Nicholson R, Maclaren OJ. 2022 Reliable and efficient parameter estimation using approximation continuum limit descriptions of stochastic models. Journal of Theoretical Biology 549, 111201.
  * [44] Yamashita T, Yamashita K, Kamimura R. 2007 A stepwise AIC method for variable selection in linear regression. Communications in Statistics - Theory and Methods 36, 2395–2403.
  * [45] Guillot C, Lecuit T. 2013 Mechanics of epithelial tissue homeostasis and morphogenesis. Science 340, 1185–1189.
  * [46] Bragulla HH, Homberger DG Structure and functions of keratin proteins in simple, stratified, keratinized and cornified epithelia. Journal of Anatomy 214, 516–559.
  * [47] Paster I, Stojadinovic O, Yin NC, Ramirez H, Nusbaum AG, Saway A, Patel SB, Khalid L, Isseroff RR, Tomic-Canic M. 2014 Epithelialization in wound healing: A comprehensive review. Advances in Wound Care 3, 445–464.
  * [48] Begnaud S, Chen T, Delacour D, Mège RM, Ladoux B. 2016 Mechanics of epithelial tissues during gap closure. Current Opinion in Cell Biology 42, 52–62.
  * [49] Paredes J, Figueiredo J, Albergaria A, Oliveira P, Carvalho J, Ribeiro AS, Caldeira J, Costa AM, Simoões-Correia J, Oliveira MJ, Pinheiro H, Pinho SS, Mateus R, Reis CA, Leite M, Fernandes MS, Schmitt F, Carneiro F, Figueiredo C, Oliveira C, Seruca R. 2012 Epithelial E- and P-cadherins: Role and clinical significance in cancer. Biochimica et Biophysica Acta 1826, 297–311.
  * [50] Hittelman WN. 2006 Genetic instability in epithelial tissues at risk for cancer. Annals of the New York Academy of Sciences 952, 1–12.
  * [51] Bezanson J, Edelman A, Karpinski S, Shah VB. 2017 Julia: A fresh approach to numerical computing. SIAM Review 59, 65–98.
  * [52] Witelski TP. 1995 Shocks in nonlinear diffusion. Applied Mathematics Letters 8, 27–32.
  * [53] Simpson MJ, Landman KA, Hughes BD, Fernando AE. 2010 A model for mesoscale patterns in motile populations. Physica A: Statistical Mechanics and its Applications 389, 1412–1424.
  * [54] Johnston ST, Baker RE, McElwain DLS, Simpson MJ. 2017 Co-operation, competition and crowding: A discrete frameworking linking allee kinetics, nonlinear diffusion, shocks and sharp-fronted travelling waves. Scientific Reports 7, 42134.
  * [55] Smith AM, Baker RE, Kay D, Maini PK. 2012 Incorporating chemical signalling factors into cell-based models of growing epithelial tissues. Journal of Mathematical Biology 65, 441–463.
  * [56] Osborne JM, Fletcher AG, Pitt-Francis JM, Maini PK, Gavaghan DJ. 2017 Comparing individual-based approaches to modelling the self-organization of multicellular tissues. PLoS Computational Biology 13, e1005387.
  * [57] Tam AKY, Simpson MJ. 2023 Pattern formation and front stability for a moving-boundary model of biological invasion and recession. Physica D: Nonlinear Phenomena 444, 133593.
  * [58] Sethian JA. 1999 Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science. Cambridge: Cambridge University Press.
  * [59] Macklin P, Lowengrub JS. 2008 A new ghost cell/level set method for moving boundary problems: Application to tumor growth. Journal of Scientific Computing 35, 266–299.
  * [60] Morris RG, Rao M. 2019 Active morphogenesis of epithelial monolayers. Physical Review E 100, 022413.
  * [61] Chang YW, Cruz-Acuña R, Tennenbaum M, Fragkopoulos AA, García AJ, Fernández-Nieves A. 2022 Quantifying epithelial cell proliferation on curved surfaces. Frontiers in Physics 10.
  * [62] Martina-Perez S, Simpson MJ, Baker RE. 2021 Bayesian uncertainty quantification for data-driven equation learning. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 477, 20210426.
  * [63] Versteeg HK, Malalasekera W. 2007 An introduction to computational fluid mechanics. Harlow: Prentice Hall 2nd edition.
  * [64] Rackauckas C, Nie Q. 2017 DifferentialEquations.jl — A performant and feature-rich ecosystem for solving differential equations in Julia. Journal of Open Research Software 5, 15.
  * [65] Hosea ME, Shampine LF. 1996 Analysis and implementation of TR-BDF2. Applied Numerical Mathematics 20, 21–37.
  * [66] Davis TA, Palamadai Natarajan E. 2010 Algorithm 907: KLU, A direct sparse solver for circuit simulation problems. ACM Transactions on Mathematical Software 37, 1–17.
  * [67] Landau HG. 1950 Heat conduction in a melting solid. The Quarterly Journal of Mechanics & Applied Mathematics 8, 81–94.
  * [68] Furzeland RM. 1980 A comparative study of numerical methods for moving boundary problems. IMA Journal of Applied Mathematics 26, 411–429.
  * [69] Amemiya T. 1985 Advanced econometrics. Cambridge, MA: Harvard University Press.
  * [70] Johnson SG. 2010 The NLopt nonlinear-optimization package. https://github.com/stevengj/nlopt. Version 2.7.1.
  * [71] Johnson SG. 2013 NLopt.jl: Package to call the NLopt nonlinear-optimization library from the Julia language. https://github.com/JuliaOpt/NLopt.jl. Version 0.6.5.

##  Appendix Appendix A Confidence bands for inaccurate continuum limits

In the paper, we show in Figure 3 a series of curves for Case Study 3 and Case
Study 4 with $k=1/5$, finding that the solution to the continuum limit is no
longer a good match to the data from the discrete model. Figure 12 shows the
confidence bands around each of these curves, showing how the uncertainty
evolves over time.

Figure 12: Complementary figure to Figure 3 in the main document, showing
inaccurate continuum limits for Case Study 3 (left column) and Case Study 4
(right column). The solid curves are the discrete densities from Equation 2.4
and the dashes curves are solutions to the continuum limit problem in Equation
2.6. The shaded regions show $95\%$ confidence bands from the mean discrete
curves at each time shown.

##  Appendix Appendix B Discrete densities at the boundaries

In the paper, we give the following formulae for computing the cell densities
from our discrete model at the boundary:

$q_{1}(t)=\dfrac{2}{x_{2}(t)-x_{1}(t)}-\dfrac{2}{x_{3}(t)-x_{1}(t)},\qquad
q_{n}(t)=\dfrac{2}{x_{n}(t)-x_{n-1}(t)}-\frac{2}{x_{n}(t)-x_{n-2}(t)},$ (42)

noting also that $x_{1}(t)=0$ in this work. In this section, we derive the
expressions for $q_{1}(t)$ and $q_{n}(t)$ and show the need for these
complicated expressions over those from Baker et al. [16], namely
$q_{1}(t)=1/(x_{2}(t)-x_{1}(t))$ and $q_{n}(t)=1/(x_{n}(t)-x_{n-1}(t))$,
through an example.

###  B.1 Derivation

We give the derivation for $q_{n}(t)$ only, as $q_{1}(t)$ is derived in the
same way. We follow the idea from Baker et al. [16], relating the cell index
$i$ to the density $q$ according to

$i(x,t)=1+\int_{0}^{x}q(y,t)\,\mathrm{d}y.$

Baker et al. [16] use $1=n-(n-1)$ together with this relationship to write

$1=\int_{x_{n-1}(t)}^{x_{n}(t)}q(y,t)\,\mathrm{d}y,$

and Baker et al. [16] then use a right endpoint rule to approximate
$q_{n}(t)$. If we instead use a trapezoidal rule, then

$1=\int_{x_{n-1}(t)}^{x_{n}(t)}q(y,t)\,\mathrm{d}y\approx\left(\dfrac{x_{n}(t)-x_{n-1}(t)}{2}\right)\left(q_{n}(t)+q_{n-1}(t)\right).$
(43)

We use this expression to solve for $q_{n}(t)$:

$\displaystyle q_{n}(t)$
$\displaystyle=\frac{2}{x_{n}(t)-x_{n-1}(t)}-q_{n-1}(t)=\frac{2}{x_{n}(t)-x_{n-1}(t)}-\frac{2}{x_{n}(t)-x_{n-2}(t)},$

which is exactly the formula in (42). We note that an alternative derivation
of this formula is to use linear extrapolation, treating the density
$1/(x_{n}(t)-x_{n-1}(t))$ as if it were placed at the cell midpoint
$(x_{n-1}(t)+x_{n}(t))/2$ rather than $x_{n}(t)$.

###  B.2 Motivation

Let us now give the motivation for why we need the modifications to the
boundary densities in (42) compared to those given in Baker et al. [16].
Consider a mechanical relaxation problem, starting with $30$ equally spaced
nodes in $0\leq x\leq 5$, taking the parameters $k=50$, $s=1/5$, $\eta=1$ and
leaving the right boundary free. Let us compare the discrete densities at
$t=2$ to those from the continuum limit, as well as estimates of the gradient
$\partial q/\partial x$ at the right boundary.

Figure 13 shows our comparisons. Focusing on the densities at the right
boundary of Figure 13(a) gives Figure 13(b), where we can see a clear
difference in the slopes of each curve. The curve obtained using the approach
of Baker et al. [16], using $q_{n}(t)=1/(x_{n}(t)-x_{n-1}(t))$, has a
different slope from the continuum limit, whereas the red curve, using
$q_{n}(t)=2/(x_{n}(t)-x_{n-1}(t))-2/(x_{n}(t)-x_{n-2}(t))$, has a slope that
is much closer to the slope of the continuum limit model at this point. These
issues become more apparent when we try to estimate $\partial q/\partial x$ at
the boundary for each time, as we would have to do in our equation learning
procedure. Shown in Figure 13(c), we see that the estimates of $\partial
q/\partial x$ that use $q_{n}(t)=1/(x_{n}(t)-x_{n-1}(t))$ do not resemble what
we expect in the continuum limit, namely $\partial q/\partial
x=H(q)=2q^{2}(1-qs)$ (using $q=q(x_{n},t)$, where $q(x,t)$ is the solution
from the continuum limit partial differential equation (PDE)). Our new
expression for $q_{n}(t)$ gives estimates for $\partial q/\partial x$ that are
much closer to $H(q)$, with $H(q)$ passing directly through the center of
these estimates across the entire time domain. Thus, our revised formulae (42)
are necessary if we want to obtain accurate estimates for the boundary
gradients.

Figure 13: Comparison of the density definitions from Baker et al. [16] to
those in (42), using data from a mechanical relaxation problem as an example.
(a) Comparing the definitions at $t=2$ together with densities from the
continuum limit PDE. The magenta rectangle shows the region that is zoomed in
on in (b). (b) Zooming in on the magenta rectangle from (a) at the right
boundary. (c) Comparing estimates of $\partial q/\partial x$ at the right
boundary using each definition along with the continuum limit boundary
condition $\partial q/\partial x=2q^{2}(1-qs)$.

##  Appendix Appendix C Numerical methods

In this section we give the details involved in solving the PDEs on the fixed
and moving domains numerically using the finite volume method [63]. We have
provided Julia packages FiniteVolumeMethod1D.jl and
MovingBoundaryProblems1D.jl to implement these methods for the fixed and
moving domains, respectively.

###  C.1 Fixed domain

We start by considering the fixed domain problem. The PDE we consider is

$\begin{array}[]{rcll}\partialderivative{q}{t}&=&\partialderivative{x}\left(D(q)\partialderivative{q}{x}\right)+R(q)&0<x<L,\,t>0,\\\\[9.0pt]
\partialderivative{q}{x}&=&0&x\in\\{0,L\\},\,t>0,\\\\[9.0pt] \end{array}$ (44)

where $L$ is the length of the domain, $D(q)$ is the nonlinear diffusivity
function, and $R(q)$ is the source term. To discretise (44), define a grid
$0=x_{1}<x_{2}<\cdots<x_{n}=L$ with $x_{i}=(n-1)\Delta x$ and $\Delta
x=L/(n-1)$. This grid enables us to define control volumes
$\Omega_{i}=[w_{i},e_{i}]$ for each $i$, where

$w_{i}=\begin{cases}x_{1}&i=1,\\\
\dfrac{1}{2}\left(x_{i-1}+x_{i}\right)&i=2,\ldots,n,\end{cases}\quad\textnormal{and}\quad
e_{i}=\begin{cases}\dfrac{1}{2}\left(x_{i}+x_{i+1}\right)&i=1,\ldots,n-1,\\\
x_{n}&i=n.\end{cases}$ (45)

The volumes of these control volumes are denoted $V_{i}=e_{i}-w_{i}$,
$i=1,\ldots,n$. We then integrate (44) over a single $\Omega_{i}$ to give

$\derivative{\bar{q}_{i}}{t}=\frac{1}{V_{i}}\left\\{D\left(q(e_{i},t)\right)\partialderivative{q(e_{i},t)}{x}-D\left(q(w_{i},t)\right)\partialderivative{q(w_{i},t)}{x}\right\\}+\bar{R}_{i},$
(46)

where

$\bar{q}_{i}=\frac{1}{V_{i}}\int_{w_{i}}^{e_{i}}q(x,t)\,\mathrm{d}x\quad\textnormal{and}\quad\bar{R}_{i}=\frac{1}{V_{i}}\int_{w_{i}}^{e_{i}}R[q(x,t)]\,\mathrm{d}x.$

To proceed, let $q_{i}=q(x_{i},t)$, $D_{i}=D(q_{i})$, $R_{i}=R(q_{i})$ and
define the following approximations:

$\begin{array}[]{rcll}\bar{q}_{i}&=&q_{i}&i=1,\ldots,n,\\\\[6.0pt]
\bar{R}_{i}&=&R_{i}&i=1,\ldots,n,\\\\[6.0pt]
D\left(q(e_{i},t)\right)&=&\frac{1}{2}\left(D_{i}+D_{i+1}\right)&i=1,\ldots,n-1,\\\\[6.0pt]
D\left(q(w_{i},t)\right)&=&\frac{1}{2}\left(D_{i-1}+D_{i}\right)&i=2,\ldots,n,\\\\[6.0pt]
\partialderivative{q(e_{i},t)}{x}&=&\frac{q_{i+1}-q_{i}}{\Delta
x}&i=1,\ldots,n-1,\\\\[6.0pt]
\partialderivative{q(w_{i},t)}{x}&=&\frac{q_{i}-q_{i-1}}{\Delta
x}&i=2,\ldots,n.\end{array}$ (47)

Using the approximations in (47), (46) becomes

$\derivative{q_{i}}{t}=\frac{1}{V_{i}}\left[\left(\frac{D_{i}+D_{i+1}}{2}\right)\left(\frac{q_{i+1}-q_{i}}{\Delta
x}\right)-\left(\frac{D_{i-1}+D_{i}}{2}\right)\left(\frac{q_{i}-q_{i-1}}{\Delta
x}\right)\right]+R_{i},$ (48)

for $i=2,\ldots,n-1$. The boundary conditions are $x=0$ and $x=L$ are
incorporated by simply setting the associated derivative term in (46) to zero,
giving

$\displaystyle\derivative{q_{1}}{t}$ $\displaystyle=\frac{1}{2V_{1}\Delta
x}\left(D_{1}+D_{2}\right)\left(q_{2}-q_{1}\right)+R_{1},$ (49)
$\displaystyle\derivative{q_{n}}{t}$ $\displaystyle=-\frac{1}{2V_{n}\Delta
x}\left(D_{n-1}+D_{n}\right)\left(q_{n}-q_{n-1}\right)+R_{n}.$ (50)

The system of ordinary differential equations (ODEs) is thus given by
(48)–(50) and defines the numerical solution to (44). In particular, letting
$\mathbf{q}^{n}=(q_{1}(t_{n}),\ldots,q_{n}(t_{n}))^{\mkern-1.5mu\mathsf{T}}$
for some time $t_{n}$, we start with
$\mathbf{q}^{0}=(q_{0}(x_{1}),\ldots,q_{0}(x_{n}))^{\mkern-1.5mu\mathsf{T}}$
using the initial condition $q(x,0)=q_{0}(x)$, and then integrate forward in
time via (48)–(50). This procedure is implemented in the Julia package
FiniteVolumeMethod1D.jl which makes use of DifferentialEquations.jl to solve
the system of ODEs with the TRBDF2(linsolve = KLUFactorization()) algorithm
[64, 65, 66].

###  C.2 Moving boundary problem

We now describe how we solve the PDEs for a moving boundary problem. The PDE
we consider is

$\begin{array}[]{rcll}\partialderivative{q}{t}&=&\partialderivative{x}\left(D(q)\partialderivative{q}{x}\right)+R(q)&0<x<L(t),\,t>0,\\\\[6.0pt]
\partialderivative{q}{x}&=&0&x=0,\,t>0,\\\\[6.0pt]
\partialderivative{q}{x}&=&H(q)&x=L(t),\,t>0,\\\\[6.0pt]
q\derivative{L}{t}&=&-E(q)\partialderivative{q}{x}&x=L(t),\,t>0.\\\\[6.0pt]
\end{array}$ (51)

We assume that $L(t)>0$ for $t\geq 0$. The discretisation starts by
transforming onto a fixed domain using the Landau transform $\xi=x/L(t)$ [67,
68, 16]. With this change of variable, (51) becomes

$\begin{array}[]{rcll}\partialderivative{q}{t}&=&\frac{\xi}{L}\derivative{L}{t}\partialderivative{q}{\xi}+\frac{1}{L^{2}}\partialderivative{\xi}\left(D(q)\partialderivative{q}{\xi}\right)+R(q)&0<\xi<1,\,t>0,\\\\[6.0pt]
\partialderivative{q}{\xi}&=&0&\xi=0,\,t>0,\\\\[6.0pt]
\partialderivative{q}{\xi}&=&LH(q)&\xi=1,\,t>0,\\\\[6.0pt]
q\derivative{L}{t}&=&-\frac{E(q)}{L}\partialderivative{q}{\xi}&\xi=1,\,t>0.\end{array}$
(52)

To now discretise (52), define $\xi_{i}=(i-1)\Delta\xi$ for $i=1,\ldots,n$,
where $\Delta\xi=1/(n-1)$, and then let

$w_{i}=\begin{cases}\xi_{1}&i=1,\\\
\dfrac{1}{2}\left(\xi_{i-1}+\xi_{i}\right)&i=2,\ldots,n,\end{cases}\quad\textnormal{and}\quad
e_{i}=\begin{cases}\dfrac{1}{2}\left(\xi_{i}+\xi_{i+1}\right)&i=1,\ldots,n-1,\\\
\xi_{n}&i=n.\end{cases}$ (53)

We then define a control volume to be the interval $\Omega_{i}=[w_{i},e_{i}]$
with volume $V_{i}=e_{i}-w_{i}$, $i=1,\ldots,n$. Next, the PDE in (52) is
integrated over this control volume to give

$\displaystyle\int_{w_{i}}^{e_{i}}\partialderivative{q}{t}\,\mathrm{d}\xi$
$\displaystyle=\frac{1}{L}\derivative{L}{t}\int_{w_{i}}^{e_{i}}\xi\partialderivative{q}{\xi}\,\mathrm{d}\xi+\int_{w_{i}}^{e_{i}}R(q)\,\mathrm{d}\xi$
$\displaystyle+\frac{1}{L^{2}}\left[D\left(q(e_{i},t)\right)\partialderivative{q(e_{i},t)}{\xi}-D\left(q(w_{i},t)\right)\partialderivative{q(w_{i},t)}{\xi}\right].$
(54)

Using integration by parts, the first integral on the right-hand side of (54)
is simply

$\int_{w_{i}}^{e_{i}}\xi\partialderivative{q}{\xi}\,\mathrm{d}\xi=e_{i}q(e_{i},t)-w_{i}q(w_{i},t)-\int_{w_{i}}^{e_{i}}q\,\mathrm{d}\xi.$

Next, define the control volume averages

$\bar{q}_{i}=\dfrac{1}{V_{i}}\int_{w_{i}}^{e_{i}}q\,\mathrm{d}\xi,\qquad\bar{R}_{i}=\dfrac{1}{V_{i}}\int_{w_{i}}^{e_{i}}R\,\mathrm{d}\xi,$

and set $q_{i}=q(\xi_{i},t)$, $D_{i}=D(q_{i})$, and $R_{i}=R(q_{i})$. With
this notation, we define the following set of approximations:

$\begin{array}[]{rcll}\bar{q}_{i}&=&q_{i}&i=1,\ldots,n,\\\\[6.0pt]
\bar{R}_{i}&=&R_{i}&i=1,\ldots,n,\\\\[6.0pt]
q(e_{i},t)&=&\frac{1}{2}\left(q_{i}+q_{i+1}\right)&i=1,\ldots,n-1,\\\\[6.0pt]
q(w_{i},t)&=&\frac{1}{2}\left(q_{i-1}+q_{i}\right)&i=2,\ldots,n,\\\\[6.0pt]
D\left(q(e_{i},t)\right)&=&\frac{1}{2}\left(D_{i}+D_{i+1}\right)&i=1,\ldots,n-1,\\\\[6.0pt]
D\left(q(w_{i},t)\right)&=&\frac{1}{2}\left(D_{i-1}+D_{i}\right)&i=2,\ldots,n,\\\\[6.0pt]
\partialderivative{q(e_{i},t)}{\xi}&=&\frac{q_{i+1}-q_{i}}{\Delta\xi}&i=1,\ldots,n-1,\\\\[6.0pt]
\partialderivative{q(w_{i},t)}{\xi}&=&\frac{q_{i}-q_{i-1}}{\Delta\xi}&i=2,\ldots,n.\end{array}$
(55)

Using these approximations, (54) becomes

$\displaystyle\displaystyle\derivative{q_{i}}{t}$
$\displaystyle\displaystyle=\frac{1}{V_{i}L}\derivative{L}{t}\left[e_{i}\left(\frac{q_{i}+q_{i+1}}{2}\right)-w_{i}\left(\frac{q_{i-1}+q_{i}}{2}\right)\right]-\frac{1}{L}\derivative{L}{t}q_{i}+R_{i}$
$\displaystyle\displaystyle+\frac{1}{V_{i}L^{2}}\left[\left(\frac{D_{i}+D_{i+1}}{2}\right)\left(\frac{q_{i+1}-q_{i}}{\Delta\xi}\right)-\left(\frac{D_{i-1}+D_{i}}{2}\right)\left(\frac{q_{i}-q_{i-1}}{\Delta\xi}\right)\right].$
(56)

The last component to handle are the boundary conditions. Since $\partial
q/\partial\xi=0$ at $\xi=0$, and since $w_{1}=\xi_{1}=0$, our discretisation
at $\xi=0$ becomes

$\derivative{q_{1}}{t}=\frac{1}{V_{1}L}\derivative{L}{t}e_{1}\left(\frac{q_{1}+q_{2}}{2}\right)-\frac{1}{L}\derivative{L}{t}q_{1}+R_{1}+\frac{1}{V_{1}L^{2}}\left(\frac{D_{1}+D_{2}}{2}\right)\left(\frac{q_{2}-q_{1}}{\Delta\xi}\right).$
(57)

The boundary condition at $\xi=1$ is $\partial q/\partial\xi=LH(q)$, thus

$\displaystyle\derivative{q_{n}}{t}$
$\displaystyle=\frac{1}{V_{n}L}\derivative{L}{t}\left[q_{n}-w_{n}\left(\frac{q_{n-1}+q_{n}}{2}\right)\right]-\frac{1}{L}\derivative{L}{t}q_{n}+R_{n}$
$\displaystyle+\frac{1}{V_{n}L^{2}}\left[D_{n}LH(q_{n})-\left(\frac{D_{n-1}+D_{n}}{2}\right)\left(\frac{q_{n}-q_{n-1}}{\Delta\xi}\right)\right].$
(58)

The remaining boundary condition is the moving boundary condition,
$q\mathrm{d}L/\mathrm{d}t=-[E(q)/L]\partial q/\partial\xi$. Since $\partial
q/\partial\xi=LH(q)$, we can write
$q_{n}\mathrm{d}L/\mathrm{d}t=-[E(q_{n})/L]LH(q_{n})=-E(q_{n})H(q_{n})$,
giving

$q_{n}\derivative{L}{t}=-E(q_{n})H(q_{n}).$ (59)

The system of ODEs (56)–(59), together with the initial conditions
$q_{i}(0)=q_{0}(\xi_{i}L(0))$ for $i=1,\ldots,n$ and $L(0)=L_{0}$, where
$q_{0}(x)$ and $L_{0}$ are the initial conditions, define our complete
discretisation. Solving these ODEs over time give values for
$q(\xi_{i},t_{j})$, for some $t_{j}$, which gets translated back in terms of
$x$ via $x_{i}=\xi_{i}L(t_{j})$. As in the fixed domain case, we solve these
ODEs using DifferentialEquations.jl together with the TRBDF2(linsolve =
KLUFactorization()) algorithm [64, 65, 66]. We provide our implementation of
this procedure in a separate Julia package, MovingBoundaryProblems1D.jl.

##  Appendix Appendix D Additional stepwise equation learning details

In this section, we give some extra details for our stepwise equation learning
procedure.

###  D.1 Discrete mechanism averaging

We start by discussing how we take multiple stochastic realisations from our
discrete cell simulations and average them into a single density function.

The discrete simulations give us $n_{s}$ identically prepared realisations
that can be averaged over to estimate the mean density curve. This average can
be estimated using a linear interpolant across each time and for each
simulation. In particular, let $n_{k}$ be the number of knots to use for the
interpolant at each time. Then, for a given time $t_{j}$, let the knots be
given by $\bar{x}_{ij}$ for $i=1,\ldots,n_{k}$. These knots are equally spaced
with $\bar{x}_{1j}=0$ and
$\bar{x}_{n_{k}j}=(1/n_{s})\sum_{\ell=1}^{n_{s}}L_{j}^{(\ell)}$, where
$L_{j}^{(\ell)}$ is the leading edge at the time $t_{j}$ from the $\ell$th
simulation. Then, letting $q^{(\ell)}(x,t_{j})$ denote the linear interpolant
of the density data at the time $t_{j}$ from the $\ell$th simulation, we
define

$\bar{q}_{ij}=\frac{1}{n_{s}}\sum_{\ell=1}^{n_{s}}q^{(\ell)}(\bar{x}_{ij},t_{j}),$
(60)

for $i=1,\ldots,n_{k}$ and $j=1,\ldots,M$. If
$q^{(\ell)}(\bar{x}_{ij},t_{j})<0$ for a given $\ell$, then we set
$q^{(\ell)}(\bar{x}_{ij},t_{j})=0$. This density data is used for computing
the system $(\mathbf{A},\mathbf{b})$ for equation learning when proliferation
is involved.

###  D.2 Derivative estimation

The equation learning system $(\mathbf{A},\mathbf{b})$ requires estimates for
the derivatives $\partial q_{ij}/\partial t$, $\partial q_{ij}/\partial x$,
$\partial^{2}q_{ij}/\partial x^{2}$, and $\mathrm{d}L_{j}/\mathrm{d}t$. To
give a formula for an estimate of these derivatives, suppose we have three
function values $\\{f_{1},f_{2},f_{3}\\}$ for some function $f(x)$ at the
points $\\{x_{1},x_{2},x_{3}\\}$, where $f_{i}=f(x_{i})$ for $i=1,2,3$. These
points do not need to be equally spaced. The Lagrange interpolating polynomial
through this data is given by

$g(x)=\frac{(x-x_{2})(x-x_{3})}{(x_{1}-x_{2})(x_{1}-x_{3})}f_{1}+\frac{(x-x_{1})(x-x_{3})}{(x_{2}-x_{1})(x_{2}-x_{3})}f_{2}+\frac{(x-x_{1})(x-x_{2})}{(x_{3}-x_{1})(x_{3}-x_{2})}f_{3},$

which can be used to estimate the derivatives via $f^{\prime}(x_{i})\approx
g^{\prime}(x_{i})$, $i=1,2,3$, and similarly for $f^{\prime\prime}(x)$. Using
this approximation, we write

$\displaystyle f^{\prime}(x_{1})$
$\displaystyle\approx\left(\frac{1}{x_{1}-x_{2}}+\frac{1}{x_{1}-x_{3}}\right)f_{1}-\frac{x_{1}-x_{3}}{(x_{1}-x_{2})(x_{2}-x_{3})}f_{2}+\frac{x_{1}-x_{2}}{(x_{1}-x_{3})(x_{2}-x_{3})}f_{3},$
(61) $\displaystyle f^{\prime}(x_{2})$
$\displaystyle\approx\frac{x_{2}-x_{3}}{(x_{1}-x_{2})(x_{1}-x_{3})}f_{1}+\left(\frac{1}{x_{2}-x_{3}}-\frac{1}{x_{1}-x_{2}}\right)f_{2}+\frac{x_{2}-x_{1}}{(x_{1}-x_{3})(x_{2}-x_{3})}f_{3},$
(62) $\displaystyle f^{\prime}(x_{3})$
$\displaystyle\approx\frac{x_{3}-x_{2}}{(x_{1}-x_{2})(x_{1}-x_{3})}f_{1}+\frac{x_{1}-x_{3}}{(x_{1}-x_{2})(x_{2}-x_{3})}f_{2}-\left(\frac{1}{x_{1}-x_{3}}+\frac{1}{x_{2}-x_{3}}\right)f_{3},$
(63) $\displaystyle f^{\prime\prime}(x_{i})$
$\displaystyle\approx\frac{2}{(x_{1}-x_{2})(x_{1}-x_{3})}f_{1}-\frac{2}{(x_{1}-x_{2})(x_{2}-x_{3})}f_{2}+\frac{2}{(x_{1}-x_{3})(x_{2}-x_{3})}f_{3},$
(64)

where (64) is valid for $i=1,2,3$.

We can use the formulae (61)–(64) to approximate our required derivatives. For
example, taking $\\{x_{1},x_{2},x_{3}\\}=\\{t_{j-1},t_{j},t_{j+1}\\}$ and
$\\{f_{1},f_{2},f_{3}\\}=\\{L_{j-1},L_{j},L_{j+1}\\}$ gives

$\derivative{L_{j}}{t}\approx\frac{L_{j+1}-L_{j-1}}{h},\quad j=2,\ldots,M-1,$
(65)

assuming the times are equally spaced with spacing $h$. The estimate for
$\mathrm{d}L_{M}/\mathrm{d}t$ is obtained by taking
$\\{x_{1},x_{2},x_{3}\\}=\\{t_{M-2},t_{M-1},t_{M}\\}$ and
$\\{f_{1},f_{2},f_{3}\\}=\\{L_{M-2},L_{M-1},L_{M}\\}$, giving

$\derivative{L_{M}}{t}\approx\frac{3L_{M}-4L_{M-1}+L_{M-2}}{2h}.$ (66)

Similarly, taking $\\{x_{1},x_{2},x_{3}\\}=\\{x_{i-1,j},x_{ij},x_{i+1,j}\\}$
and $\\{f_{1},f_{2},f_{3}\\}=\\{q_{i-1,j},q_{ij},q_{i+1,j}\\}$ gives

$\displaystyle\partialderivative[2]{q_{ij}}{x}$
$\displaystyle\approx\frac{2}{(x_{i-1,j}-x_{i,j})(x_{i-1,j}-x_{i+1,j})}q_{i-1,j}-\frac{2}{(x_{i-1,j}-x_{ij})(x_{ij}-x_{i+1,j})}q_{ij}$
$\displaystyle+\frac{2}{(x_{i-1,j}-x_{i+1,j})(x_{ij}-x_{i+1,j})},$ (67)

for $i=2,\ldots,n_{j}-1$ and $j=1,\ldots,M$, where $n_{j}$ is the number of
nodes at $t=t_{j}$. The remaining derivatives can be obtained similarly,
ensuring that the appropriate finite difference (backward, central, or
forward) is taken for the given point.

The only exception to these rules are for $\partial q/\partial x$ at the
boundaries. We find that using simple forward and backward differences there
gives better results than with (61) and (62), so we use

$\partialderivative{q_{1j}}{x}\approx\frac{q_{2j}-q_{1j}}{x_{2j}-x_{1j}},\quad\partialderivative{q_{n_{j}j}}{x}\approx\frac{q_{n_{j}j}-q_{n_{j}-1,j}}{x_{n_{j}j}-x_{n_{j}-1,j}}.$
(68)

###  D.3 Matrix pruning

We now discuss our approach to matrix pruning, wherein we discard points from
our equation learning matrix $\mathbf{A}$ that do not help to improve our
estimates for $\bm{\theta}$. The approach we take is inspired from the data
thresholding idea from VandenHeuvel et al. [42].

To start with our approach, let
$\mathbf{q}=(q_{12},\ldots,q_{n_{M}M})^{\mkern-1.5mu\mathsf{T}}$ be the vector
of all discrete densities, letting $n_{j}$ be the number of nodes at the time
$t=t_{j}$, excluding the densities from the initial condition. Then, take the
threshold tolerance $0\leq\tau_{q}<1/2$ and compute the interval
$(\mathcal{Q}_{\tau_{q}}^{\mathbf{q}},\mathcal{Q}_{1-\tau_{q}}^{\mathbf{q}})$,
where $\mathcal{Q}_{\tau}^{\mathbf{y}}$ denotes the $100\tau\%$ quantile of
the vector $\mathbf{y}$. With these intervals, we only include a row in the
matrix $\mathbf{A}$ from a given point $(x_{ij},t_{j})$ if
$\mathcal{Q}_{\tau_{q}}^{\mathbf{q}}\leq
q_{ij}\leq\mathcal{Q}_{1-\tau_{q}}^{\mathbf{q}}$.

By choosing the threshold $\tau_{q}$ appropriately, we can significantly
improve the estimates for $\bm{\theta}$ as we only include the most relevant
data for estimation, excluding all points with relatively low or high density.
Similar thresholds can be defined for the other quantities
$|\partial\mathbf{q}/\partial x|,|\partial^{2}\mathbf{q}/\partial
x^{2}|,|\partial\mathbf{q}/\partial t|$, and
$|\mathrm{d}\mathbf{L}/\mathrm{d}t|$, defining these vectors similarly to
$\mathbf{q}$, for example
$|\partial_{t}\mathbf{q}|=(|\partial_{t}q_{12}|,\ldots,|\partial_{t}q_{n_{M}M}|)^{\mkern-1.5mu\mathsf{T}}$,
with respective threshold tolerances $0\leq\tau_{\partial q/\partial
x},\tau_{\partial^{2}q/\partial x^{2}},\tau_{\partial q/\partial
t},\tau_{\mathrm{d}L/\mathrm{d}t}<1/2$.

##  Appendix Appendix E Additional examples

In this section, we give some additional case studies to further demonstrate
our method, exploring different force law and proliferation laws, and
enforcing conservation of mass together with a discussion about enforcing
equality constraints in general.

###  E.1 Enforcing conservation of mass

In the paper, we discussed at the end of Case Study 2 that it could be
possible to enforce mass conservation to fix the issue with $D(q)\neq E(q)$,
noting that mass conservation requires $D(q(L(t),t))=E(q(L(t),t))$. In this
section, we consider the results when we fix $D(q)=E(q)$ so that mass is
conserved from the outset.

This change $D(q)=E(q)$ is reasonably straightforward to implement in the
algorithm, simply replacing the boundary condition (17) so that

$q(L(t),t)\derivative{L(t)}{t}=-D\left(q(L(t),t)\right)\partialderivative{q(L(t),t)}{x}.$
(69)

This constraint $D(q)=E(q)$ also needs to be reflected in the matrix
$\mathbf{A}$. This is simple to do in this case. Previously, our matrix system
took the block diagonal form

$\begin{bmatrix}\mathbf{A}_{1}&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\mathbf{A}_{2}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{A}_{3}\end{bmatrix}\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{h}\\\
\bm{\theta}^{e}\end{bmatrix}=\begin{bmatrix}\mathbf{b}_{1}\\\
\mathbf{b}_{2}\\\ \mathbf{b}_{3}\end{bmatrix}.$ (70)

With the constraint $D(q)=E(q)$, (70) becomes

$\begin{bmatrix}\mathbf{A}_{1}&\mathbf{0}\\\ \mathbf{0}&\mathbf{A}_{2}\\\
\mathbf{A}_{3}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{h}\end{bmatrix}=\begin{bmatrix}\mathbf{b}_{1}\\\
\mathbf{b}_{2}\\\ \mathbf{b}_{3}\end{bmatrix}.$ (71)

We note that, if we wanted to enforce this constraint in Case Study 4, where
$\mathbf{A}_{1}=[\mathbf{A}^{d}~{}\mathbf{A}^{r}]$, with $\mathbf{A}^{d}$ and
$\mathbf{A}^{r}$ defined from (31), then we instead have

$\begin{bmatrix}\mathbf{A}^{d}&\mathbf{A}^{r}&\mathbf{0}\\\
\mathbf{0}&\mathbf{0}&\mathbf{A}_{2}\\\
\mathbf{A}_{3}&\mathbf{0}&\mathbf{0}\end{bmatrix}\begin{bmatrix}\bm{\theta}^{d}\\\
\bm{\theta}^{r}\\\
\bm{\theta}^{h}\end{bmatrix}=\begin{bmatrix}\mathbf{b}_{1}\\\
\mathbf{b}_{2}\\\ \mathbf{b}_{3}\end{bmatrix}.$ (72)

Let us now consider the results with mass conservation. We use the same
parameters that were used to produce the results in Figure 7. In particular,
we save the solution at $M=200$ equally spaced times between $t_{1}=0$ and
$t_{M}=15$, $\tau_{q}=0.35$, $\tau_{\mathrm{d}L/\mathrm{d}t}=0.1$, and we
start with all coefficients initially inactive. The results we obtain are
shown in Table 8 and Figure 14. We see that the form we learn for $D(q)$, and
hence for $E(q)$ also, is close to the continuum limit $50/q^{2}$, and
similarly $H(q)$ is a good match; note that $H(q)$ is only evaluated at the
boundary densities, which is approximately $5$ for $t>0$, so indeed $H(q)$
matches the continuum limit. Looking to Figure 14(a)–(b), the results are
indistinguishable from the continuum limit, which is also what we found in
Figure 7 before we enforced conservation of mass.

Table 8: Stepwise equation learning results for Case Study 2, using the basis expansions (28), saving the results at $M=200$ equally spaced times between $t_{1}=0$ and $t_{M}=15$, pruning with $\tau_{q}=0.35$ and $\tau_{\mathrm{d}L/\mathrm{d}t}=0.1$, starting with all terms inactive, and enforcing conservation of mass with $D(q)=E(q)$. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step. Step | $\theta_{1}^{d}$ | $\theta_{2}^{d}$ | $\theta_{3}^{d}$ | $\theta_{1}^{h}$ | $\theta_{2}^{h}$ | $\theta_{3}^{h}$ | $\theta_{4}^{h}$ | $\theta_{5}^{h}$ | Loss
---|---|---|---|---|---|---|---|---|---
1 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | 0.000 | -3.371
2 | 0.000 | 0.000 | 0.000 | 0.000 | -0.025 | 0.000 | 0.000 | 0.000 | -2.371
3 | 0.000 | 47.413 | 0.000 | 0.000 | -0.025 | 0.000 | 0.000 | 0.000 | -1.706
4 | 0.000 | 47.413 | 0.000 | 0.000 | 0.443 | 0.000 | 0.000 | -0.004 | -0.688
Figure 14: Stepwise equation learning results from Table 8. (a) Comparisons of
the discrete density profiles (solid curves) with those from the learned PDE
(dashed curves), plotted at the times $t=0,5,10,25,50,100$ in black, red,
blue, green, orange, and purple, respectively. (b) As in (a), except comparing
the leading edges. (c)–(e) are comparisons of the learned forms of $D(q)$,
$H(q)$, and $E(q)$ compared to the forms from the continuum limit.

####  E.1.1 Imposing linear equality constraints generically

We note that this approach to implementing the constraint $D(q)=E(q)$
requiring such a significant change to the matrix system, giving (71), and to
the boundary condition (69), might suggest that the modularity of our approach
weakens here. This does not need to be the case, and so let us briefly remark
about how constraints such as $D(q)=E(q)$, or any other linear constraints,
could be alternatively implemented in our approach seamlessly, further
demonstrating the modularity.

Suppose we take our system $\mathbf{A}\bm{\theta}=\mathbf{b}$, with
$\mathbf{A}\in\mathbb{R}^{m\times p}$, $\bm{\theta}\in\mathbb{R}^{p}$, and
$\mathbf{b}\in\mathbb{R}^{m}$, and suppose we have constraints of the form
$\mathbf{Q}^{\mkern-1.5mu\mathsf{T}}\bm{\theta}=\mathbf{c}$ where
$\mathbf{Q}\in\mathbb{R}^{p\times c}$ and $\mathbf{c}\in\mathbb{R}^{c}$, where
$c<p$ and $\mathbf{Q}$ has full rank. The constrained least squares estimator
for $\bm{\theta}$ subject to these constraints, denoted
$\hat{\bm{\theta}}^{c}$, is then given by

$\hat{\bm{\theta}}^{c}=\hat{\bm{\theta}}-\left(\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{A}\right)^{-1}\mathbf{Q}\left[\mathbf{Q}^{\mkern-1.5mu\mathsf{T}}\left(\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{A}\right)^{-1}\mathbf{Q}\right]^{-1}\left(\mathbf{Q}^{\mkern-1.5mu\mathsf{T}}\hat{\bm{\theta}}-\mathbf{c}\right),$
(73)

where
$\hat{\bm{\theta}}=(\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{A})^{-1}\mathbf{A}^{\mkern-1.5mu\mathsf{T}}\mathbf{b}$
is the unconstrained least squares estimator for
$\mathbf{A}\bm{\theta}=\mathbf{b}$ [69]. Using this formulation, imposing
$D(q)=E(q)$ is simple to enforce without changing the boundary condition or
the matrix $\mathbf{A}$, simply using $\mathbf{c}=\mathbf{0}_{3\times 1}$ and

$\mathbf{Q}=\blockarray{cc}\\\ \block{c[c]}\bm{\theta}^{d}&\mathbf{I}_{3}\\\
\bm{\theta}^{h}\mathbf{0}_{5\times 3}\\\ \bm{\theta}^{e}-\mathbf{I}_{3}\\\
=\blockarray{cccc}\\\ \block{c[ccc]}\theta_{1}^{d}&100\\\ \theta_{2}^{d}010\\\
\theta_{3}^{d}001\\\ \theta_{1}^{h}000\\\ \theta_{2}^{h}000\\\
\theta_{3}^{h}000\\\ \theta_{4}^{h}000\\\ \theta_{5}^{h}000\\\
\theta_{1}^{e}-100\\\ \theta_{2}^{e}0-10\\\ \theta_{3}^{e}00-1\\\ ,$

where $\mathbf{I}_{n}$ and $\mathbf{0}_{m\times n}$ denote the $n$-square
identity matrix and $m\times n$ zero matrix, respectively. This does not solve
the problem entirely, though, since we also have coefficients that we force to
zero throughout the stepwise procedure. These zeros constraints can also be
imposed by including additional columns of $\mathbf{Q}$. For example, if
$\theta_{1}^{h}$ and $\theta_{2}^{d}$ are inactive, then $\mathbf{Q}$ becomes

$\mathbf{Q}=\blockarray{cccc}\\\
\block{c[ccc]}\bm{\theta}^{d}&\mathbf{I}_{3}\mathbf{e}_{2}^{d}\mathbf{0}_{3\times
1}\\\ \bm{\theta}^{h}\mathbf{0}_{5\times 3}\mathbf{0}_{5\times
1}\mathbf{e}_{1}^{h}\\\ \bm{\theta}^{e}-\mathbf{I}_{3}\mathbf{0}_{3\times
1}\mathbf{0}_{3\times 1}\\\ ,$ (74)

where $\mathbf{e}_{2}^{d}=(0,1,0)^{\mkern-1.5mu\mathsf{T}}$ and
$\mathbf{e}_{1}^{h}=(1,0,0,0,0)^{\mkern-1.5mu\mathsf{T}}$. In particular, each
inactive coefficient $\theta_{i}$ corresponds to a new column with a one in
the row corresponding to that coefficient. Note that $\mathbf{Q}$ in (74) can
be further written as $\mathbf{Q}=[\mathbf{Q}_{1}~{}\mathbf{Q}_{2}]$, where
$\mathbf{Q}_{1}$ are the user-provided constraints $D(q)=E(q)$ and
$\mathbf{Q}_{2}$ are the constraints imposed by the inactive coefficients,
making it easy to incorporate constraints in this manner. Additional care is
required to ensure that there are no redundant constraints represented by
$\mathbf{Q}_{1}$ and $\mathbf{Q}_{2}$ as $\mathbf{Q}$ must be full rank. For
example, imposing $\theta_{1}^{d}=0$ and $\theta_{1}^{e}=0$ together with the
constraint $\theta_{1}^{d}=\theta_{1}^{e}$ from $D(q)=E(q)$ can be represented
using only two constraints rather than three, and the associated matrix

$\mathbf{Q}=\blockarray{cccc}\\\
\block{c[ccc]}\bm{\theta}^{d}&\mathbf{I}_{3}\mathbf{e}_{1}^{d}\mathbf{0}_{3\times
1}\\\ \bm{\theta}^{h}\mathbf{0}_{5\times 3}\mathbf{0}_{5\times
1}\mathbf{0}_{5\times 1}\\\ \bm{\theta}^{e}-\mathbf{I}_{3}\mathbf{0}_{3\times
1}\mathbf{e}_{1}^{h}\\\ ,$ (75)

where $\mathbf{e}_{1}^{h}=(1,0,0)^{\mkern-1.5mu\mathsf{T}}$, only has rank $4$
rather than the full rank $5$. This could be dealt with by finding a basis for
the column space of $\mathbf{Q}$, replacing $\mathbf{Q}$ with the
corresponding matrix of basis vectors.

To summarise this discussion, it is straightforward to implement our procedure
with the ability to enforce linear equality constraints, allowing for
additional constraints, such as conservation of mass, to be enforced. This is
easy to code without breaking the modularity of the approach and requiring a
significant change to the procedure that would be cumbersome to implement by
increasing the complexity of the corresponding code.

###  E.2 A piecewise proliferation law

In this section, we consider the problem described in Section 3.3 of Murphy et
al. [15]. This problem given by Murphy et al. [15] is used to demonstrate a
case where the solution of the continuum limit no longer gives a good match
with averaged data from the discrete model, as the value of $k$ used is too
low relative to the proliferation rate. Here, we show how our method can learn
an accurate continuum model in this case.

The example is as follows. We consider $F(\ell_{i})=k(s-\ell_{i})$ as usual,
taking $k=10^{-4}$ and $s=0$, but our proliferation law is now given by

$G(\ell_{i})=\begin{cases}0&0\leq\ell_{i}<\ell_{p},\\\
\beta&\ell_{i}\geq\ell_{p},\end{cases}$ (76)

where $\ell_{p}=0.2$ is the proliferation threshold and $\beta=10^{-2}$. We
use $\Delta t=10^{-2}$ for the proliferation events. The initial condition
places $n=41$ equally spaced nodes in $[0,10]$ so that $\ell_{i}=0.25$ at
$t=0$ for each of the $40$ cells. In Figure 15, we show a comparison of the
discrete data from this problem with the solution of the continuum limit. We
also compare the cell numbers $N(t)$, where the cell numbers from the PDE
$q(x,t)$ are obtained via $N(t)=\int_{0}^{10}q(x,t)\,\mathrm{d}x$. We see that
the densities from the solution of the continuum limit reach a capacity at
$50$ cells, while the discrete model instead reaches $80$ cells. Note that the
densities appear jagged in Figure 13 due to the combination of the averaging
procedure from Section D.1 with the variance of the densities for moderate
$t$; a better averaging method could be to build the knots at each time $t$
based on the node positions themselves, but we do not consider that here as it
does not impact the results.

Figure 15: Comparison of the solution of the piecewise proliferation law
problem with the solution of continuum limit, where
$F(\ell_{i})=k(s-\ell_{i})$ and $G(\ell_{i})=\beta$ for $\ell_{i}\geq\ell_{p}$
and $G(\ell_{i})=0$ otherwise, using $k=10^{-4}$, $s=0$, $\ell_{p}=0.2$,
$\eta=1$, $\beta=10^{-2}$, and $\Delta t=10^{-2}$. (a) The solid curves are
the discrete densities, and the dashed curves are the densities from the
solution of the continuum limit. The arrow shows the direction of increasing
time. The density profiles are shown at the times $t=0,10,50,100,250,500$ in
black, red, blue, green, orange, and purple, respectively. (b) Comparison of
the number of cells from the discrete model with that computed from the
solution of the continuum limit, using $N(t)=\int_{0}^{10}q(x,t)\,\mathrm{d}x$
for the continuum limit case. In (a)–(b), the discrete results are averaged
over $1000$ identically prepared realisations, using $n_{k}=100$ knots for the
averaging procedure described in Section D.1.

The continuum limit for this problem is

$D(q)=\frac{10^{-4}}{q^{2}}\quad\text{and}\quad
R(q)=\begin{cases}0&q>1/\ell_{p},\\\ 10^{-2}q&q\leq 1/\ell_{p}.\end{cases}$

This suggests one possible basis expansion to use for $R(q)$ in our equation
learning procedure, with the aim to learn an appropriate continuum
approximation to the results in Figure 15, could be

$R(q)=\left[\theta_{0}^{r}+\theta_{1}^{r}q+\theta_{2}^{r}q^{2}+\theta_{3}^{r}q^{3}\right]\mathbb{I}\left(q\leq\frac{1}{\ell_{p}}\right),$

where $\mathbb{I}(A)$ is the indicator function for the set $A$. We find that
this does not lead to any improved model for this problem, and so we instead
consider a polynomial model:

$R(q)=\theta_{0}^{r}+\theta_{1}^{r}q+\theta_{2}^{r}q^{2}+\theta_{3}^{r}q^{3}+\theta_{4}^{r}q^{4}+\theta_{5}^{r}q^{5}.$
(77)

For $D(q)$, this mechanism does not appear to be relevant in this example,
with the results that follow all giving visually indistinguishable regardless
of whether $D(q)=0$ or $D(q)=10^{-4}/q^{2}$. Thus, we do not bother learning
it in this case, simply fixing $D(q)=10^{-4}/q^{2}$; if we do not fix $D(q)$,
we just end up learning $D(q)=0$ in the results that follow. With (77) and
$D(q)=10^{-4}/q^{2}$, the results we obtain are shown in Table 9 and Figure
16.

Step | $\theta_{1}^{r}$ | $\theta_{2}^{r}$ | $\theta_{3}^{r}$ | $\theta_{4}^{r}$ | $\theta_{5}^{r}$ | $\theta_{6}^{r}$ | Loss
---|---|---|---|---|---|---|---
1 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -1.63
2 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | -1.53
3 | 0.077 | -0.0096 | 0.00 | 0.00 | 0.00 | 0.00 | -6.22
Table 9: Equation learning results for the piecewise proliferation law problem
in Figure 15, fixing $D(q)=10^{-4}/q^{2}$ and using the expansion of $R(q)$ in
(77). The discrete data is averaged over $1000$ identically prepared
realisations with $n_{k}=100$ knots for interpolating, and the solution is
saved every $0.1$ units of time between $t=0$ and $t=500$. Figure 16: Equation
learning results for the piecewise proliferation law problem in Figure 15,
using the results from Table 9. (a) Comparison of the averaged discrete
densities (solid curves) with the solution of the learned PDE (dashed). The
arrow shows the direction of increasing time. The arrow shows the direction of
increasing time. The density profiles are shown at the times
$t=0,10,50,100,250,500$ in black, red, blue, green, orange, and purple,
respectively.(b) Comparison of the cell numbers. (c) Comparison of the learned