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

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–5], exclusion processes [6–9], mechanical models of epithelial tissues [10–17], hydrodynamics [18,19] and a variety of other types of individual-based models [1,20–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]. By 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 [31,33,34]. 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 $\mathrm{\Delta }$ with probability $P_{\text{m}}\in [0,1]$ per time step of duration $\tau$. The stochastic proliferation process involves undergoing proliferation with probability $P_{\text{p}}\in [0,1]$ per unit time step. The continuum limit description of this kind of discrete process can be written as [33]

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 $\mathrm{\partial }q/\mathrm{\partial }t=\mathcal{N}(q,\mathcal{D},\mathit{\theta })$, where $q$ is the population density, $\mathcal{N}$ is some nonlinear function parametrized by $\mathit{\theta }$, $\mathcal{D}$ is a collection of differential operators and $\mathit{\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 $\mathit{\theta }$, such as sparse regression or nonlinear optimization [32,37–39,41,42], representing $\mathcal{N}$ as a library of functions [37–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 $\mathit{\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 parallelizable, allowing for rapid exploration of results with different discrete parameters and different forms of $\mathcal{N}$ for a given dataset. 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 emphasize 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 $\mathit{\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 [10,17]. Epithelial tissues are biological tissue composed of cells, organized 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–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

Following Murray et al. [10] and Baker et al. [16], we suppose that we have a set of nodes $x_1,\dots ,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,\dots ,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 \text{d}{x}_i(t)/\text{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]

• Download figure
• Open in new tab
• Download PowerPoint

The dynamics governed by (2.1)–(2.2) describe a system in which cells mechanically relax. Following previous work [10,14,16], 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+\mathrm{\Delta }t)$ for some small duration $\mathrm{\Delta }t$. We let the probability that the $i$th cell proliferates be given by $G_i\mathrm{\Delta }t$, where $G_i=G({\ell }_i)$ for some length-dependent proliferation law $G({\ell }_i)>0$. As represented in figure 1c, 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 electronic supplementary material, §S5. The implementation of the solution to these equations (2.1)–(2.2) 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 emphasize that, without proliferation, we need only solve (2.1)–(2.2) 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. By contrast, incorporating proliferation means that we need to consider several identically prepared realizations 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 density of the cells left and right of $x_i(t)$, giving

3. Continuum-discrete comparison

We now consider four biologically motivated case studies to illustrate the performance of the continuum limit description (2.6). These case studies are represented schematically in figure 1d–g. Case Studies 1 and 3, shown in figure 1d,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 in the left column of figure 2 and a set of averaged density profiles for each problem in the right column of figure 2. Case Studies 1 and 3 initially place $30$ nodes in $0\le x\le 5$ and $30$ nodes in $25\le x\le 30$, or equivalently $n=60$ with 28 cells in $0\le x\le 5$ and $28$ cells in $25\le x\le 30$, spacing the nodes uniformly within each subinterval. Case Studies 2 and 4 initially place $60$ equally spaced nodes in $0\le x\le 5$.

• Download figure
• Open in new tab
• Download PowerPoint

The problems shown in figure 2 use parameter values such that the solution of the continuum limit (2.6) 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, $\mathrm{\Delta }t={10}^{-2}$, $K=15$ and $\beta =0.15$. The accuracy of the continuum limit is clearly evident in the density profiles in 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.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 electronic supplementary material, §S1, we show the $95\mathrm{\%}$ confidence regions for each curve in figure 3, where we find that the solutions have much greater variance compared with the corresponding curves in figure 2 where $k=50$.

• Download figure
• Open in new tab
• Download PowerPoint

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 (2.6) is accurate. We demonstrate this in Case Studies 1–4 in §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.

(a) 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 (2.6), 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,\dots ,{\theta }_d\}$ and $d$ basis functions $\{{\phi }_1,\dots ,{\phi }_d\}$, meaning $D(q)$ can be represented as

$$D(q)=\sum _{i=1}^d{\theta }_i{\phi }_i(q).$$4.1

These basis functions could be any univariate function of $q$, for example the basis could be $\{{\phi }_1,{\phi }_2,{\phi }_3\}=\{1/q,1/q^2,1/q^3\}$ with $d=3$. In this work, we impose the constraint that $D(q)\ge 0$ for $q_{min}\le q\le 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–54], we wish to avoid the possibility of having negative nonlinear diffusivity functions and our results support this approach.

The aim is to estimate $\mathit{\theta }={({\theta }_1,\dots ,{\theta }_d)}^{\mathsf{T}}$ in (4.1). We use ideas similar to the basis function approach from VandenHeuvel et al. [42], using (4.1) to construct a matrix problem for $\mathit{\theta }$. In particular, let us take the PDE (2.6), with $R(q)=0$ and $H(q)=0$, and expand the spatial derivative term so that we can isolate the ${\theta }_k$ terms,

$$\frac{\mathrm{\partial }{q}_{ij}}{\mathrm{\partial }t}=\sum _{k=1}^d\{\frac{\mathrm{d}{\phi }_k(q_{ij})}{\mathrm{d}q}{\left(\frac{\mathrm{\partial }{q}_{ij}}{\mathrm{\partial }x}\right)}^2+{\phi }_k(q_{ij})\frac{{\mathrm{\partial }}^2{q}_{ij}}{\mathrm{\partial }{x}^2}\}{\theta }_k,$$4.2

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 $\mathrm{\partial }{q}_{ij}/\mathrm{\partial }t$, $\mathrm{\partial }{q}_{ij}/\mathrm{\partial }x$ and ${\mathrm{\partial }}^2{q}_{ij}/\mathrm{\partial }{x}^2$ in (4.2); 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 electronic supplementary material, §S4. We also emphasize that, while (4.2) appears similar to results in [37,38], 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 $\mathit{\theta }$, defined below and in (4.7). We save the solution to the discrete problems (2.1)–(2.2) at $M$ times $0=t_1<t_2<\cdots <t_M$ so that $i\in \{1,\dots ,n\}$ and $j\in \{2,\dots ,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 (4.2) into a rectangular matrix problem $\mathbf{\text{A}}\mathit{\theta }=\mathbf{\text{b}}$, where the $r$th row in $\mathbf{\text{A}}$, $r=1,\dots ,n(M-1)$, corresponding to the point $(x_{ij},t_j)$ is given by ${\mathbf{\text{a}}}_{ij}\in {\mathbb{R}}^{1\times d}$, where

$${\mathbf{\text{a}}}_{ij}=\left[\begin{array}{ccc}{\displaystyle \frac{\mathrm{d}{\phi }_1(q_{ij})}{\mathrm{d}q}}{\left({\displaystyle \frac{\mathrm{\partial }{q}_{ij}}{\mathrm{\partial }x}}\right)}^2+{\phi }_1(q_{ij}){\displaystyle \frac{{\mathrm{\partial }}^2{q}_{ij}}{\mathrm{\partial }{x}^2}},& \dots ,& {\displaystyle \frac{\mathrm{d}{\phi }_d(q_{ij})}{\mathrm{d}q}}{\left({\displaystyle \frac{\mathrm{\partial }{q}_{ij}}{\mathrm{\partial }x}}\right)}^2+{\phi }_d(q_{ij}){\displaystyle \frac{{\mathrm{\partial }}^2{q}_{ij}}{\mathrm{\partial }{x}^2}}\end{array}\right],$$4.3

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

$$\mathbf{\text{A}}=\left[\begin{array}{c}{\mathbf{\text{a}}}_{12}\\ {\mathbf{\text{a}}}_{22}\\ ⋮\\ {\mathbf{\text{a}}}_{nM}\end{array}\right]\in {\mathbb{R}}^{n(M-1)\times d}\text{and}\mathbf{\text{b}}=\left[\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }{q}_{12}}{\mathrm{\partial }t}}\\ {\displaystyle \frac{\mathrm{\partial }{q}_{22}}{\mathrm{\partial }t}}\\ ⋮\\ {\displaystyle \frac{\mathrm{\partial }{q}_{nM}}{\mathrm{\partial }t}}\end{array}\right]\in {\mathbb{R}}^{n(M-1)\times 1}.$$4.4

The solution of $\mathbf{\text{A}}\mathit{\theta }=\mathbf{\text{b}}$, given by $\mathit{\theta }={({\mathbf{\text{A}}}^{\mathsf{T}}\mathbf{\text{A}})}^{-1}{\mathbf{\text{A}}}^{\mathsf{T}}\mathbf{\text{b}}$, is obtained by minimizing the residual $||\mathbf{\text{A}}\mathit{\theta }-\mathbf{\text{b}}|{|}_2^2$, which keeps all terms present in the learned model. We expect, however, just as in (2.7), that $\mathit{\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{\text{A}}\mathit{\theta }=\mathbf{\text{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,\dots ,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 $\mathit{\theta }$, starting with ${\mathcal{A}}_1=\mathcal{I}$. The set of indices of zero values in $\mathit{\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:

(i)	Let the vector ${\mathit{\theta }}_{\mathcal{A}}$ denote the solution to $\mathbf{\text{A}}\mathit{\theta }=\mathbf{\text{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 ${\mathit{\theta }}_{\mathcal{A}}$ by solving the reduced problem in which the inactive columns of $\mathbf{\text{A}}$ are not included. The vector with $\mathcal{A}={\mathcal{A}}_k$ at step $k$ is denoted ${\mathit{\theta }}_k$. With this definition, we compute the sets $${\mathcal{M}}_k^{+}=\{{\mathit{\theta }}_{{\mathcal{A}}_k\cup \{i\}}:i\notin {\mathcal{A}}_k\}\text{and}{\mathcal{M}}_k^{-}=\{{\mathit{\theta }}_{{\mathcal{A}}_k\setminus \{i\}}:i\in {\mathcal{A}}_k\}.$$4.5 ${\mathcal{M}}_k^{+}$ is the set of all coefficient vectors $\mathit{\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=\{{\mathit{\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 ${\mathit{\theta }}_k$.
(ii)	Choose one of the vectors in ${\mathcal{M}}_k$ by defining a loss function $\mathcal{L}(\mathit{\theta })$, $$\underset{\text{loss}}{\underbrace{\mathcal{L}(\mathit{\theta })}}=\underset{\text{goodness of fit}}{\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;\mathit{\theta })}{q_{ij}}\right)}^2\right]}}+\underset{\text{model complexity}}{\underbrace{||\mathit{\theta }|{|}_0}},$$4.6 where $q(x,t;\mathit{\theta })$ is the solution of the PDE (2.6) with $R(q)=H(q)=0$ and $D(q)$ uses the coefficients $\mathit{\theta }$ in (4.1), $q(x_{ij},t_j;\mathit{\theta })$ is the linear interpolant of the PDE data at $t=t_j$ evaluated at $x=x_{ij}$, and $||\mathit{\theta }|{|}_0$ is the number of non-zero terms in $\mathit{\theta }$. This loss function balances the goodness of fit with model complexity. If, for some $\mathit{\theta }$, $D(q)<0$ within $q_{min}\le q\le q_{max}$, which we check by evaluating $D(q)$ at $n_c=100$ equally spaced points in $q_{min}\le q\le q_{max}$, we set $\mathcal{L}(\mathit{\theta })=\mathrm{\infty }$. With this loss function, we compute the next coefficient vector $${\mathit{\theta }}_{k+1}=\underset{\mathit{\theta }\in {\mathcal{M}}_k}{argmin}\mathcal{L}(\mathit{\theta }).$$4.7 If ${\mathit{\theta }}_{k+1}=\mathbf{\text{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.
(iii)	If ${\mathit{\theta }}_{k+1}={\mathit{\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 ${\mathit{\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{\text{A}}\mathit{\theta }-\mathbf{\text{b}}|{|}_2^2$, that has been useful for a range of previously considered problems [32,37,38]. We do not consider the regression error in this work as we find that it typically leads to poorer estimates for $\mathit{\theta }$ compared with controlling the density errors as we do in (4.7).

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 ${\phi }_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},$$4.8

and we expect to learn $\mathit{\theta }={(0,50,0)}^{\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 set-up, and starting with all coefficients initially active so that ${\mathcal{A}}_1=\{1,2,3\}$, we obtain the results in table 2. The first iteration gives us ${\mathit{\theta }}_1$ such that $D(q)<0$ for some range of $q$ as we show in figure 4a, and so we assign $\mathcal{L}({\mathit{\theta }}_1)=\mathrm{\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 ${\mathit{\theta }}_2={(-1.46,47.11,0)}^{\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 ${\mathit{\theta }}_3={(0,43.52,0)}^{\mathsf{T}}$, which is close to the continuum limit. We emphasize that this final ${\mathit{\theta }}_3$ is a least-squares solution with the constraint ${\theta }_1={\theta }_3=0$, thus there is no need to refine ${\mathit{\theta }}_3$ further by eliminating ${\theta }_1$ and ${\theta }_3$ directly in (4.8), as the result would be the same. Comparing the densities from the solution of the learned PDE with $\mathit{\theta }={\mathit{\theta }}_3$ with the discrete densities in figure 5a, we see that the curves are nearly visually indistinguishable near the centre, but there are some visually discernible discrepancies near the boundaries. We show the form of $D(q)$ at each iteration in figure 4a, where we observe that the first iteration captures only the higher densities, the second iteration captures the complete range of densities, and the third iteration removes a single term which gives no noticeable difference.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

Table 2. Stepwise equation learning results for the density data for Case Study 1: fixed boundaries using the basis expansion (4.8), 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.

Collapse

step	${\theta }_1$	${\theta }_2$	${\theta }_3$	loss
1	−5.97	70.73	−27.06	$\mathrm{\infty }$
2	−1.46	47.11	0.00	−4.33
3	0.00	43.52	0.00	−5.18

Table 3. Improved results for Case Study 1: fixed boundaries from table 2, now using matrix pruning so that densities outside of the $10\mathrm{\%}$ and $90\mathrm{\%}$ density quantiles are not included. Coefficients highlighted in blue show the coefficient chosen to be removed or added at the corresponding step.

Collapse

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$

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 $\mathit{\theta }$. Visual inspection of the space–time diagram in figure 2a 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 $\mathit{\theta }$ significantly, and so it is useful to only include important points in the construction of $\mathbf{\text{A}}$. To resolve this issue, we choose to only include points if the associated densities fall between the $10\mathrm{\%}$ and $90\mathrm{\%}$ quantiles for the complete set of densities, which we refer to by density quantiles; more details on this pruning procedure are given in electronic supplementary material, §S4. This choice of density quantiles is made using trial and error, starting at $0\mathrm{\%}$ and $100\mathrm{\%}$, respectively, and shrinking the quantile range until suitable results are obtained. When we apply this pruning and reconstruct $\mathbf{\text{A}}$, we obtain the improved results in table 3 and associated densities in figure 5b. Compared with table 2, we see that the coefficient estimates for $\mathit{\theta }$ all lead to improved losses, and our final model now has $\mathit{\theta }={(0,49.83,0)}^{\mathsf{T}}$, which is much closer to the continuum limit, as we see in figure 5b) where the solution curves are now visually indistinguishable everywhere. Moreover, we show in figure 4b 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{\text{A}}$.

(b) 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)\ne 0$. We do not consider proliferation, giving $R(q)=0$ in (2.6).

The equation learning procedure for this case study is similar to Case Study 1, namely we expand $D(q)$ as in (4.1) and constrain $D(q)\ge 0$. In addition to learning $D(q)$, we need to learn $H(q)$ and the evolution equation describing the free boundary. In (2.6), this evolution equation is given by a conservation statement, $q\text{d}L/\text{d}t=-D(q)\mathrm{\partial }q/\mathrm{\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\frac{\mathrm{d}L}{\mathrm{d}t}=-E(q)\frac{\mathrm{\partial }q}{\mathrm{\partial }x},$$4.9

at $x=L(t)$ for $t>0$. While (4.9) could lead to local loss of conservation at the moving boundary, our approach is to allow 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)\ge 0$ so that (4.9) 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{\phi }_i^d(q),H(q)=\sum _{i=1}^h{\theta }_i^h{\phi }_i^h(q),\hspace{0.17em}E(q)=\sum _{i=1}^h{\theta }_i^e{\phi }_i^e(q).$$4.10

The matrix system for ${\mathit{\theta }}^d={({\theta }_1^d,\dots ,{\theta }_d^d)}^{\mathsf{T}}$ is the same as it was in Case Study 1 in (4.4), which we now write as ${\mathbf{\text{A}}}^d{\mathit{\theta }}^d={\mathbf{\text{b}}}^d$ with ${\mathbf{\text{A}}}^d\in {\mathbb{R}}^{n(M-1)\times d}$ and ${\mathbf{\text{b}}}^d\in {\mathbb{R}}^{n(M-1)\times 1}$ given by $\mathbf{\text{A}}$ and $\mathbf{\text{b}}$ in (4.4), and we can construct two other independent matrix systems for ${\mathit{\theta }}^h={({\theta }_1^h,\dots ,{\theta }_h^h)}^{\mathsf{T}}$ and ${\mathit{\theta }}^e={({\theta }_1^e,\dots ,{\theta }_e^e)}^{\mathsf{T}}$. To construct these matrix systems, for a given boundary point $(x_{nj},t_j),$ we write

$$\frac{\mathrm{\partial }{q}_{nj}}{\mathrm{\partial }x}=\sum _{k=1}^h{\theta }_k^h{\phi }_k^h(q_{nj})\mathrm{and}{q}_{nj}\frac{\mathrm{d}{L}_j}{\mathrm{d}t}=-\frac{\mathrm{\partial }{q}_{nj}}{\mathrm{\partial }x}\sum _{k=1}^e{\theta }_k^e{\phi }_k^e(q_{nj}),$$4.11

where $L_j=x_{nj}$ is the position of the leading edge at $t=t_j$. In (4.11), we assume that $L_j$ can be approximated by a smooth function so that $\text{d}{L}_j/\text{d}t$ can be defined. With (4.11) we have ${\mathbf{\text{A}}}^h{\mathit{\theta }}^h={\mathbf{\text{b}}}^h$ and ${\mathbf{\text{A}}}^e{\mathit{\theta }}^e={\mathbf{\text{b}}}^e$, where

$${\mathbf{\text{A}}}^h=\left[\begin{array}{ccc}{\phi }_1^h(q_{12})& \cdots & {\phi }_h^h(q_{12})\\ ⋮& \ddots & ⋮\\ {\phi }_1^h(q_{nM})& \cdots & {\phi }_h^h(q_{nM})\end{array}\right],{\mathbf{\text{b}}}^h=\left[\begin{array}{c}{\displaystyle \frac{\mathrm{\partial }{q}_{12}}{\mathrm{\partial }x}}\\ ⋮\\ {\displaystyle \frac{\mathrm{\partial }{q}_{nM}}{\mathrm{\partial }x}}\end{array}\right]$$4.12

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

$${\mathbf{\text{A}}}^e=\left[\begin{array}{ccc}{\phi }_1^e(q_{12}){\displaystyle \frac{\mathrm{\partial }{q}_{n2}}{\mathrm{\partial }x}}& \cdots & {\phi }_e^e(q_{12}){\displaystyle \frac{\mathrm{\partial }{q}_{n2}}{\mathrm{\partial }x}}\\ ⋮& \ddots & ⋮\\ {\phi }_1^e(q_{nM}){\displaystyle \frac{\mathrm{\partial }{q}_{nM}}{\mathrm{\partial }x}}& \cdots & {\phi }_e^e(q_{nM}){\displaystyle \frac{\mathrm{\partial }{q}_{nM}}{\mathrm{\partial }x}}\end{array}\right],{\mathbf{\text{b}}}^e=-\left[\begin{array}{c}{q}_{n2}{\displaystyle \frac{\mathrm{d}{L}_2}{\mathrm{d}t}}\\ ⋮\\ q_{nM}{\displaystyle \frac{\mathrm{d}{L}_M}{\mathrm{d}t}}\end{array}\right]$$4.13

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

$$\mathbf{\text{A}}=\mathrm{diag}({\mathbf{\text{A}}}^d,{\mathbf{\text{A}}}^h,{\mathbf{\text{A}}}^e)\in {\mathbb{R}}^{(n+2)(M-1)\times (d+h+e)},\mathbf{\text{b}}=\left[\begin{array}{c}{\mathbf{\text{b}}}^d\\ {\mathbf{\text{b}}}^h\\ {\mathbf{\text{b}}}^e\end{array}\right]\in {\mathbb{R}}^{(n+2)(M-1)\times 1},$$4.14

we obtain

$$\mathbf{\text{A}}\mathit{\theta }=\mathbf{\text{b}},\mathit{\theta }=\left[\begin{array}{c}{\mathit{\theta }}^d\\ {\mathit{\theta }}^h\\ {\mathit{\theta }}^e\end{array}\right]\in {\mathbb{R}}^{(d+h+e)\times 1}.$$4.15

The solution of $\mathbf{\text{A}}\mathit{\theta }=\mathbf{\text{b}}$ is the combined solution of the individual linear systems as $\mathbf{\text{A}}$ is block diagonal. Estimates for ${\mathit{\theta }}^d$, ${\mathit{\theta }}^h$ and ${\mathit{\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{\text{A}}\mathit{\theta }=\mathbf{\text{b}}$ in (4.15), we augment the loss function (4.6) to incorporate information about the location of the moving boundary. Letting $L(t;\mathit{\theta })$ denote the leading edge from the solution of the PDE (2.6) with parameters $\mathit{\theta }$, the loss function is

$$\begin{array}{rl}\underset{\text{loss}}{\underbrace{\mathcal{L}(\mathit{\theta })}}& \hspace{0.17em}=\stackrel{\text{density goodness of fit}}{\overbrace{\log \left[\frac{1}{n(M-1)}\sum _{j=2}^M\sum _{i=1}^n{\left(\frac{q_{ij}-q(x_{ij},t_j;\mathit{\theta })}{q_{ij}}\right)}^2\right]}}\\ & \hspace{0.17em}+\underset{\text{leading edge goodness of fit}}{\underbrace{\log \left[\frac{1}{M-1}\sum _{j=2}^M{\left(\frac{L_j-L(t_j;\mathit{\theta })}{L_j}\right)}^2\right]}}+\underset{\text{model complexity}}{\underbrace{||\mathit{\theta }|{|}_0}}.\end{array}$$4.16

Let us now apply our stepwise equation learning procedure with (4.15) and (4.16). 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

$$\begin{array}{rl}& D(q)={\displaystyle \frac{{\theta }_1^d}{q}}+{\displaystyle \frac{{\theta }_2^d}{q^2}}+{\displaystyle \frac{{\theta }_3^d}{q^3}},\\ & 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\\ \mathrm{and}& E(q)={\displaystyle \frac{{\theta }_1^e}{q}}+{\displaystyle \frac{{\theta }_2^e}{q^2}}+{\displaystyle \frac{{\theta }_3^e}{q^3}}.\end{array}\}$$4.17

With these expansions, we expect to learn ${\mathit{\theta }}^d={(0,50,0)}^{\mathsf{T}}$, ${\mathit{\theta }}^h={(0,2,-0.4,0,0)}^{\mathsf{T}}$ and ${\mathit{\theta }}^e={(0,50,0)}^{\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\mathrm{\%}$ and $65\mathrm{\%}$ 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 2b indicates that mechanical relaxation takes place over a relative short interval, which means that working with data in $0<t\le 100$ can lead to a poor outcome.

• Download figure
• Open in new tab
• Download PowerPoint

Table 4. Stepwise equation learning results for Case Study 2: free boundaries, using the basis expansions (4.17), 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\mathrm{\%}$ and $65\mathrm{\%}$ 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.

Collapse

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	$-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\le t\le 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\ge 15$, for example $L(t)$ in figure 6b 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 $\text{d}L/\text{d}t$ is neither too large nor too small. We implement this by excluding all points $(x_{nj},t_j)$ from the construction of $({\mathbf{\text{A}}}^e,{\mathbf{\text{b}}}^e)$ in (4.13) such that $\text{d}{L}_j/\text{d}t$ is outside of the $10\mathrm{\%}$ or $90\mathrm{\%}$ quantiles of the vector $(\text{d}{L}_2/\text{d}t,\dots ,\text{d}{L}_M/\text{d}t)$, called the velocityquantiles.

• Download figure
• Open in new tab
• Download PowerPoint

Table 5. Stepwise equation learning results for Case Study 2: free boundaries, using the basis expansions (4.17), 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\mathrm{\%}$ and $65\mathrm{\%}$ 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.

Collapse

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	−3.37
2	0.00	0.00	0.00	0.00	−0.03	0.00	0.00	0.00		0.00	0.00	−2.37
3	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.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

Implementing thresholding on $\text{d}L/\text{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 electronic supplementary material, §S5, 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)$.