# Forecasting and predicting stochastic agent-based models of cell migration
with biologically-informed neural networks

John T. Nardini

###### Abstract

Collective migration, or the coordinated movement of many individuals, is an
important component of many biological processes, including wound healing,
tumorigenesis, and embryo development. Spatial agent-based models (ABMs) are
often used to model collective migration, but it is challenging to thoroughly
study these models’ behavior due to their random and computational nature.
Modelers often overcome these obstacles by coarse-graining discrete ABM rules
into continuous mean-field partial differential equation (PDE) models. These
models are advantageous because they are fast to simulate; unfortunately,
these PDE models can poorly predict ABM behavior (or even be ill-posed) at
certain parameter values. In this work, we describe how biologically-informed
neural networks (BINNs) can be used to learn BINN-guided PDE models that are
capable of accurately predicting ABM behavior. In particular, we show that
BINN-guided PDE simulations can forecast future ABM data not seen during model
training. Additionally, we demonstrate how to predict ABM data at previously-
unexplored parameter values by combining BINN-guided PDE simulations with
multivariate interpolation. We highlight these results using three separate
ABMs that consist of rules on agent pulling and/or adhesion. Surprisingly,
BINN-guided PDEs can accurately forecast and predict ABM data with a one-
compartment PDE when the mean-field PDE is ill-posed or requires two
compartments. While we focus our presentation on the biological applications,
this work is broadly applicable to studying many systems that exhibit the
collective migration of individuals.

##  1 Introduction

Cellular migration is an important biological phenomenon involved in embryo
development, wound repair, and tumorigenesis [1, 2]. While the mechanics
driving _individual_ cell migration are now well understood, we have not yet
how many cells coordinate collective _population-level_ migration. In some
contexts, collective migration results from external cell stimuli, such as
empty space, cell-cell interactions, chemical signals & gradients, etc.
Additionally, many diseases and their respective complications develop when
cells abuse these cues [3]. For example, metastatic cancers are associated
with a loss of intercellular connections within epithelial tissues [4], and a
common complication of diabetes is nonhealing wounds where high inflammatory
markers and low paracrine signaling levels prevent cells from entering and
repairing a wound [5]. There is thus a current need to establish the impacts
of each of these stimuli on collective migration, which will provide insight
into tissue homeostasis, disease progression, and effective drug therapy
development.

Mathematical modeling is a useful tool to infer how physical and chemical cues
drive collective cell migration [3, 6, 7, 8]. In particular, stochastic agent-
based models (ABMs) are a widely-used modeling framework where autonomous
agents mimic individual cells [9, 10]. ABMs are advantageous because they
capture the discrete and stochastic nature of many biological processes [7].
These models are especially influential for cell biology, where modelers can
easily implement rules governing the impacts of different stimuli on agent
actions, such as the effects of cell-cell interactions on agent migration [7,
8, 11, 12]. One can further introduce population heterogeneity into an ABM by
incorporating several agent types that perform different rules. Despite the
many benefits of ABMs, there are limitations on their use: in particular,
their model simulations are computationally intensive and time-consuming to
perform [13, 10]. These restraints prevent modelers from efficiently exploring
how model parameters alter model outputs which, in turn, make it challenging
to thoroughly understand the effects of each stimuli on collective migration.
As such, there is a need for the development of methods to address ABMs’
computational expenses by 1. forecasting future model output from limited
simulations, and 2. predicting ABM data at previously-unexplored parameter
values [13, 14, 15].

One of the most commonly-used approaches for ABM forecasting and prediction
includes coarse-graining ABM rules into a continuous differential equation
(DE) model [9, 10]. ABMs are converted into ordinary DEs (ODEs) when tracking
a time-varying output (e.g., agent density) that is spatially homogenous [13].
Partial DEs (PDEs) are suitable for describing spatially-varying ABMs [10].
Mean-field PDEs that are coarse-grained from rules on the effects of cell-cell
interaction on cell migration often take the form of parabolic PDEs with
nonlinear density-dependent diffusion rates [3, 6, 7, 12]. These equations can
be written as

$\dfrac{\partial T}{\partial t}=\nabla\cdot\big{(}\mathcal{D}(T)\nabla
T\big{)},$ (1)

where $T=T(x,t)$ denotes the total spatiotemporal agent density and
$\mathcal{D}(T)$ is the density-dependent agent diffusion rate. Such DE models
are useful surrogates for ABMs because they are cheap and efficient to
simulate and are amenable to analytical methods, which modelers can use to
precisely infer how model parameters impact their outputs. However, these
models are unable to provide insight into individual level behavior and can
lead to poor ABM predictions (which can also be ill-posed) for many parameter
values [9, 12]. For example, Baker and Simpson 2010 [9] demonstrated that the
mean-field ODE model for a population growth ABM only accurately predict ABM
data when agents proliferate slowly. Chapelle and Yates 2019 [7] coarse-
grained rules on the effects of cell pulling on agent migration into multiple
PDE models; while these models accurately predict ABM data, the authors found
their accuracy decreases with more complex model rules. Anguige and Schmeiser
2009 [6] used a stochastic space-jump model to study how cell adhesion impacts
collective migration and found that the resulting mean-field PDE model is ill-
posed for large adhesion values. Modelers may improve DE models’ predictive
capability by implementing pair-wise interactions or moment closure
approximations, but the resulting models are often complicated and may
significantly increase computation time [9, 16].

Equation learning (EQL) is a new area of research on the development and
application of algorithms to discover the dynamical systems model that best
describes a dataset [17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27]. Brunton et
al. 2016 [17] introduced what is now a widely-used least squares EQL approach
that uses sparse regression to learn a simple DE model from a user-specified
library of several candidate terms. This method has proven very successful in
recovering informative models from simulated ODE and PDE data as well as
experimental data [28]. There is a growing understanding that EQL methods can
aid the prediction of ABM data [13, 29, 30, 31]. For example, we recently
demonstrated that the least squares EQL approach learns ODE equations that
accurately describe simulated ABM data, even when the collected data is
incomplete or sparsely sampled [13]. Supekar et al. 2023 [31] coupled this
method with spectral basis representation data to discover PDE models that
capture the emergent behavior found in active matter ABMs. Another popular EQL
approach includes physics-informed neural networks (PINNs), where modelers
embed physical knowledge into the training procedure for artificial neural
networks (ANNs) [32, 33, 34, 35, 36]. Trained PINN models can predict complex,
sparse, and noisy data while also obeying known physical principles. Lagergren
et al. 2020 [26] extended the PINNs framework by replacing physics-based
mechanistic terms with function-approximating multi-layer perceptions (MLPs)
to develop the biologically-informed neural network (BINN) methodology. As a
result, BINN models can learn PDE models with nonlinear and density-dependent
modeling terms, such as the forms of Equation (1) that result from coarse-
graining various agent rules on the impacts of cell interactions on collective
migration. BINNs thus present a promising tool for ABM forecasting and
prediction, but determining how they can be used to learn predictive DE models
remains an open area of research.

In this work, we demonstrate how to combine BINNs and PDE model simulations to
forecast and predict ABM behavior. Our approach leverages BINNs’ vast data and
modeling approximation capability with the computational efficiency of PDE
models to develop a potent ABM surrogate modeling tool. In particular, we
demonstrate how to accurately forecast future ABM data at a fixed parameter
value by training a BINN model to precomputed data and then simulating the
BINN’s learned PDE model. We extend this approach to predict ABM behavior at
new parameter values by training BINNs to data from multiple parameter values
and performing multivariate interpolation over the learned modeling terms.
This prediction is performed by incorporating the interpolated modeling terms
into a simulated PDE model. Interpolation provides a straightforward approach
for learning PDE modeling terms that works well, though more complex
methodologies, such as ANNs or Gaussian processes, could be used [37].

The frequent use of ABMs to study economics, social sciences, engineering, and
many other areas has led to previous work on predictive surrogate ABM modeling
[38, 39, 40, 41]. This research is closely related to the field of designing
computer experiments, where modelers implement computationally efficient
statistical methods to emulate high-fidelity computer simulations [42, 43]. In
a typical study, modelers calibrate the chosen surrogate model to several high
fidelity computer simulations, and the calibrated surrogate model is utilized
to perform a certain task, such as identifying sensitive ABM parameters,
predicting new dynamics from the high fidelity simulation, or estimating its
parameters from data. Modelers must choose a surrogate model to use: Gaussian
processes are the most popular thanks to the influential work of [37]. The
Bayesian approximation error method is another widely-used technique [44, 45,
46], and ANNs have gained traction in recent years [15, 47]. While these
‘black-box’ model choices have proven successful in practice, they typically
ignore domain expertise on the high-fidelity simulation, which limits the
interpretability of their analyses. In this work, we implement a ‘gray-box’
approach for ABM prediction by training predictive BINN models to discover
computationally efficient surrogate PDE models [48]. Visual inspection of the
PDE modeling terms enables us to interpret how model parameters impact ABM
behavior at the population level. Our work is similar to [10], who built a
statistical model to infer the discrepancy between ABM simulations and their
coarse-grained ODE approximations; parameters with high discrepancies indicate
the assumptions underlying model coarse-graining are invalid. In this previous
study, incorporating the discrepancy model into the data’s statistical model
allowed for accurate ABM parameter estimation.

We illustrate the performance of BINNs in guiding ABM forecasting and
prediction using three separate ABMs consisting of different rules on the
impacts of cell interactions on agent migration. The first model implements
rules on cell pulling and is borrowed from [7]. The second model is a discrete
version of the space-jump model from [6] on cell adhesion. We introduce the
third model, which consists of two separate subpopulations, each of which
performs either cell pulling or adhesion rules. These models highlight how
agent interactions impact collective migrations as well as some limitations of
mean-field PDE models. Namely, the mean-field PDE model for the Adhesion model
is ill-posed for large cell adhesion values, and the mean-field model for the
Pulling & Adhesion model consists of two separate compartments for each
subpopulation; to the best of our knowledge, it is not possible to convert
this two-compartment model into a single-compartment PDE describing the total
population. Our BINN-guided approach learns a single-compartment PDE models
capable of forecasting and predicting data from all three ABMs over all
parameter values considered. We begin this work in Section 2 by presenting the
three ABMs as well as their mean-field PDE models. In Section 3, we discuss
our data analysis methods on BINNs training, multivariate interpolation, and
PDE simulation. In Section 4, we detail our results on implementing these
methods for ABM forecasting and prediction and finish this section with a
brief discussion on the computational expenses of each data analysis method.
We conclude these results and suggest areas for future work in Section 5.

##  2 Coarse-graining collectively migrating ABMs into PDE models

We present three ABMs that model how various cell-cell interactions, namely
cell pulling and adhesion, impact collective cell migration. All models are
two-dimensional cellular automata models with pulling agents that perform cell
pulling rules and/or adhesive agents that perform rules on cell adhesion. The
first model is borrowed from [7] and consists only of pulling agents; the
second model is inspired by the stochastic space jump model from [6] and
consists only of adhesive agents; to the best of our knowledge, we are the
first to study the third model, which consists of both pulling and adhesive
agents.

We detail our notation and ABM implementation in Section 2.1 and then present
the rules for each ABM and their mean-field PDE models in Section 2.2.

###  2.1 Implementation and notation details

Each model is simulated in the spatial domain $(x,y)\in[0,X]\times[0,Y]$. We
choose $X=200\text{ and }Y=40$ to represent a thin rectangle where collective
migration primarily occurs along the $x$-dimension and is not affected by the
boundary in this dimension. We represent this space with a two-dimensional
lattice with square lattice sites with length $\Delta=1$ to imitate a typical
cell length. The $(i,j)^{\text{th}}$ lattice site is centered at
$(x_{i},y_{j})$, where $x_{i}=(i-0.5)\Delta,\ i=1,\dots,X$, and
$y_{j}=(j-0.5)\Delta,\ j=1,\dots,Y.$ Each model is an exclusion process,
meaning that each agent can only occupy one lattice site at a time, and each
lattice site is occupied by at most one agent. The variables $P_{i,j}(t)$,
$H_{i,j}(t)$, and $0_{i,j}(t)$ denote the probabilities that lattice site
$(i,j)$ is occupied by a pulling agent, adhesive agent, or empty at time $t$,
respectively.

All model simulations are initialized by populating 75% of the lattice sites
in the middle 20% of columns, e.g., 75% of the lattice sites in
$\\{(x_{i},y_{j})\\}_{j=1}^{Y}$ are initially occupied for $i=80,\dots,120.$
All other columns are initially empty. Reflecting boundary conditions are used
at the boundaries of lattice to enforce a no-flux condition in the spatial
domain. Let $N^{(r)}_{P}(x_{i},t)$ and $N^{(r)}_{H}(x_{i},t)$ denote the
number of pulling and adhesive agents, respectively, in the $i^{\text{th}}$
column for $i=1,\dots,X$ from the $r^{\text{th}}$ of $R$ identically prepared
ABM simulations. To estimate the pulling and adhesive agent densities in the
$i^{\text{th}}$ column from the $r^{\text{th}}$ simulation, we compute

$P^{(r)}(x_{i},t)=\dfrac{N^{(r)}_{P}(x_{i},t)}{Y}\text{ and
}H^{(r)}(x_{i},t)=\dfrac{N^{(r)}_{H}(x_{i},t)}{Y},\ \ i=1,\dots,X,$

respectively. The total agent density is then estimated by

$T^{(r)}(x_{i},t)=P^{(r)}(x_{i},t)+H^{(r)}(x_{i},t).$

To estimate the averaged pulling, adhesive, and total agent density in the
$i^{\text{th}}$ column from $R$ identically prepared ABM simulations, we
compute:

$\displaystyle\langle P^{ABM}(x_{i},t)\rangle$
$\displaystyle=\dfrac{1}{R}\sum_{r=1}^{R}P^{(r)}(x_{i},t);$
$\displaystyle\langle H^{ABM}(x_{i},t)\rangle$
$\displaystyle=\dfrac{1}{R}\sum_{r=1}^{R}H^{(r)}(x_{i},t);\text{ and }$
$\displaystyle\langle T^{ABM}(x_{i},t)\rangle$
$\displaystyle=\dfrac{1}{R}\sum_{r=1}^{R}T^{(r)}(x_{i},t),\ \ \ \text{ for
}i=1,\dots,X.$

Figure 1: ABM rules on migration, pulling, and adhesion. a) When an agent
performs a migration event, it chooses one of the four cardinal directions to
move towards with equal probability. b) A migration event requires the lattice
site in the chosen migration direction to be empty; otherwise, the migration
event is aborted. c) Rules A-F dictate the rules on agent migration, pulling,
and adhesion. Rule A prescribes how a pulling agent (blue circle) migrates
when there is no neighboring agent. Rule B prescribes how a pulling agent
migrates and attempts to pull a neighboring pulling agent with it. Rule C
prescribes how an adhesive agent (red hexagon) migrates when there is no
neighboring agent. Rule D prescribes how a neighboring adhesive agent attempts
to adhere to a migrating adhesive agent and abort its migration event. Rule E
prescribes how a migrating pulling agent attempts to pull its neighboring
adhesive agent, while the adhesive agent attempts to adhere to the pulling
agent. Rule F prescribes how a migrating adhesive agent and neighboring
pulling agent do not interact with each other.

###  2.2 ABM Rules and their mean-field PDE models

The rules governing agent pulling and adhesion from all models are visually
depicted in Figure 1, and the parameters for each rule are described in Table
1. During all rules, an agent chooses one of its four neighboring lattice site
to move into with equal probability. The migration event is aborted if the
chosen site is already occupied. Rules A, B, and E are initiated when a
pulling agent attempts to migrate. This occurs with rate $r_{m}^{pull}$,
meaning that pulling agents attempt to perform one of these rules over an
infinitesimal time interval of length $dt$ with probability $r_{m}^{pull}dt$.
Rules C, D, and F are initiated when an adhesive agent attempts to migrate,
which occurs with rate $r_{m}^{adh}$. The effective rates in Figure 1 document
the rate at which each lattice site configuration at time $t$ changes to the
updated lattice site configuration at time $t+\Delta t$. We simulate each ABM
using the Gillespie algorithm, which we provide for the Pulling & Adhesion ABM
in Algorithm S1 in appendix D.

Parameter | Description | Range
---|---|---
$r_{m}^{pull}$ | Pulling agent migration rate | $[0,\infty)$
$r_{m}^{adh}$ | Adhesive agent migration rate | $[0,\infty)$
$p_{pull}$ | Probability of successful pulling event | $[0,1]$
$p_{adh}$ | Probability of successful adhesion event | $[0,1]$
$\alpha$ | Proportion of adhesive agents | $[0,1]$
Table 1: Description of the ABM parameters involved in each ABM. The
proportion of pulling agents in each simulation is given by $1-\alpha$.

####  2.2.1 The Pulling Model

The Pulling model consists of pulling agents that migrate with rate
$r_{m}^{pull}$ and perform rules A and B from Figure 1. Suppose a pulling
agent at lattice site $(i,j)$ chooses to move rightwards into site $(i+1,j)$.
If the lattice site $(i-1,j)$ is unoccupied, then the agent performs Rule A
and moves into site $(i+1,j)$. If the lattice site $(i-1,j)$ is occupied, then
the agent attempts Rule B on agent pulling. This event succeeds with
probability $p_{pull}$, and the agent moves to site $(i+i,j)$ and pulls its
neighbor into lattice site $(i,j)$. This event fails with probability
$1-p_{pull}$, in which the agent moves into site $(i+i,j)$ but the neighbor
remains at lattice site $(i,j-1)$. These rules can be described by the
following trimolecular reaction rates:

$\displaystyle 0_{i-1,j}+P_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{r_{m}^{pull}/4}$ $\displaystyle
0_{i-1,j}+0_{i,j}+P_{i+1,j},$ (Rule A) $\displaystyle
P_{i-1,j}+P_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{p_{pull}r_{m}^{pull}/4}$ $\displaystyle
0_{i-1,j}+P_{i,j}+P_{i+1,j},$ (Rule B.1) $\displaystyle
P_{i-1,j}+P_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{(1-p_{pull})r_{m}^{pull}/4}$ $\displaystyle
P_{i-1,j}+0_{i,j}+P_{i+1,j}.$ (Rule B.2)

Equivalent reactions govern agent migration and pulling in the other three
directions.

In Appendix A.1, we show that Rules A and B can be coarse grained into the
Pulling ABM’s mean-field PDE model:

$\dfrac{\partial P}{\partial t}=\nabla\cdot\left(\mathcal{D}^{pull}(P)\nabla
P\right),\ \ \
\mathcal{D}^{pull}(P)=\dfrac{r_{m}^{pull}}{4}\left(1+3p_{pull}P^{2}\right)$
(2)

where $P=P(x,y,t)$ denotes the spatiotemporal pulling agent density. In Figure
2(a-f), we find that a simulation of Equation (2) closely matches
$P^{(1)}(x,t)$ over time for
$\bm{p}=(r_{m}^{pull},p_{pull})^{T}=(1.0,0.5)^{T}.$

Figure 2: ABM simulation snapshots and the mean-field PDE models for the
Pulling, Adesion, and Pulling & Adhesion ABMs. Blue pixels denote pulling
agents and red pixels denote adhesive agents. All ABMs were simulated on
rectangular 200$\times$40 lattices. (a-c) Snapshots of the Pulling ABM for
$r_{m}^{pull}=1.0,p_{pull}=0.5$. (d-f) The output spatiotemporal pulling agent
density (blue ‘x’ marks) is plotted against the solution of the mean-field PDE
(solid blue line) given by Equation (2). (g-i) Snapshots of the Adhesion ABM
for $r_{m}^{adh}=1.0,p_{adh}=0.5$. (j-l) The output spatiotemporal adhesive
agent density (red dots) is plotted against the solution of the mean-field PDE
(dashed red line) given by Equation (3). (m-o) Snapshots of the Pulling &
Adhesion ABM for
$r_{m}^{pull}=1.0,r_{m}^{adh}=0.25,p_{pull}=0.0.33,p_{adh}=0.33,\alpha=0.5$.
(p-r) The output spatiotemporal pulling and adhesive agent densities are
plotted against the solution of the mean-field PDE given by Supplementary
Equation (28).

####  2.2.2 The Adhesion Model

The Adhesion model consists of adhesive agents that migrate with rate
$r_{m}^{adh}$ and perform rules C and D from Figure 1. Suppose an adhesive
agent at lattice site $(i,j)$ chooses to move rightwards into site $(i+1,j)$.
If the lattice site $(i-1,j)$ is unoccupied, then the agent performs Rule C
and moves into site $(i+1,j)$. If the lattice site $(i-1,j)$ is occupied, then
the neighboring agent attempts Rule D to adhere to the migrating agent and
abort their movement. This event succeeds with probability $p_{adh}$, and
neither agent changes its location. This adhesion event fails with probability
$1-p_{adh}$, and the migratory agent moves to site $(i+i,j)$ and the neighbor
remains at lattice site $(i,j-1)$. These rules can be described by the
following trimolecular reaction rates:

$\displaystyle 0_{i-1,j}+H_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{r_{m}^{adh}/4}$ $\displaystyle
0_{i-1,j}+0_{i,j}+H_{i+1,j},$ (Rule C) $\displaystyle
H_{i-1,j}+H_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{(1-p_{adh})r_{m}^{adh}/4}$ $\displaystyle
H_{i-1,j}+0_{i,j}+H_{i+1,j}.$ (Rule D)

In Appendix A.2, we show that Rules C and D can be coarse grained into the
Adhesion ABM’s mean-field PDE model:

$\dfrac{\partial H}{\partial t}=\nabla\cdot\left(\mathcal{D}^{adh}(H)\nabla
H\right),\ \ \ \
\mathcal{D}^{adh}(H)=\dfrac{3r_{m}^{adh}}{4}\left(p_{adh}\left(H-\dfrac{2}{3}\right)^{2}+1-\dfrac{4p_{adh}}{3}\right)$
(3)

where $H=H(x,y,t)$ denotes the spatiotemporal adhesive agent density. In
Figure 2(g-l), we find that a simulation of Equation (3) closely matches
$H^{(1)}(x,t)$ over time for $\bm{p}=(r_{m}^{adh},p_{adh})^{T}=(1.0,0.5)^{T}.$

It is important to note that $\mathcal{D}^{adh}(H)$ from Equation (3) becomes
negative for some density values when $p_{adh}>0.75$. This PDE fails to
provide an ABM prediction at these parameter values because negative diffusion
is ill-posed [6].

####  2.2.3 The Pulling & Adhesion Model

The Pulling & Adhesion model consists of both pulling and adhesive agents. The
parameter $\alpha\in(0,1)$ denotes the portion of adhesive agents in the
simulation, and $(1-\alpha)$ denotes the portion of pulling agents in the
simulation. This model implements Rules A-F from Figure 1. Rules A-D are
unchanged from their descriptions in Sections 2.2.1 and 2.2.2. If a pulling
agent at lattice site $(i,j)$ chooses to move rightwards into site $(i+1,j)$
while an adhesive agent occupies site $(i-i,j)$, then Rule E dictates the
agents’ attempts to pull and adhere to each other. The migratory pulling agent
succeeds with probability $p_{pull}$ and moves to site $(i+1,j)$ while pulling
the neighboring adhesive agent into site $(i,j)$; the neighboring adhesive
agent successfully aborts the pulling agent’s migration event with probability
$p_{adh}$; both agents fail with probability $1-p_{adh}-p_{pull}$ and the
pulling agent moves to site $(i+1,j)$ while the adhesive agent remains at site
$(i-1,j)$. Based on our definition of this rule, it is not possible that both
the pulling and adhesion events succeed, so the parameters must satisfy $0\leq
p_{pull}+p_{adh}\leq 1$. Rule E can be described by the following trimolecular
reaction rate:

$\displaystyle H_{i-1,j}+P_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{p_{pull}r_{m}^{pull}/4}$ $\displaystyle
0_{i-1,j}+H_{i,j}+P_{i+1,j},$ (Rule E.1) $\displaystyle
H_{i-1,j}+P_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{(1-p_{adh}-p_{pull})r_{m}^{pull}/4}$ $\displaystyle
H_{i-1,j}+0_{i,j}+P_{i+1,j}.$ (Rule E.2)

If an adhesive agent at lattice site $(i,j)$ chooses to move rightwards into
site $(i+1,j)$ while a pulling agent occupies site $(i-i,j)$, then Rule F
dictates that the adhesive agent moves into site $(i+1,j)$ and the pulling
agent remains at site $(i-1,j)$. Rule F can be described by the following
trimolecular reaction rate:

$\displaystyle P_{i-1,j}+H_{i,j}+0_{i+1,j}$
$\displaystyle\xrightarrow{r_{m}^{adh}/4}$ $\displaystyle
P_{i-1,j}+0_{i,j}+H_{i+1,j}.$ (Rule F)

In Appendix C, we show that Rules A-F can be coarse-grained into the Pulling &
Adhesion ABM’s mean-field PDE model given by Supplementary Equation (28). This
two-compartment PDE describes the spatiotemporal densities of pulling agents,
$P(x,y,t)$, and adhesive agents, $H=H(x,y,t)$. In Figure 2(m-r), we find that
the $P$ and $H$ compartments from a simulation of Supplementary Equation (28)
closely $P^{(1)}(x,t)$ and $H^{(1)}(x,t)$, respectively, over time for
$\bm{p}=(r_{m}^{pull},r_{m}^{adh},p_{pull},p_{adh},\alpha)^{T}=(1.0,0.25,0.33,0.33,0.5)^{T}.$

To the best of our knowledge, it is not possible to convert Rules A-F into a
single-compartment PDE model describing the _total_ agent density,
$T=T(x,y,t)=H(x,y,t)+P(x,y,t)$.

##  3 Data analysis methods

We simulate all three models (the Pulling, Adhesion, and Pulling & Adhesion
ABMs) over a range of agent migration, pulling, and adhesion parameter values.
We represent the model parameters by the vector $\bm{p}$. Each model
simulation outputs 100 snapshots of agent configurations over time; from each
simulation, we generate the one-dimensional agent density along the
$x$-dimension over time. We average these densities over $R=25$ simulations to
obtain the final output ABM density, $\langle T^{ABM}(x,t;\bm{p})\rangle$. We
use the data from the first 75 timepoints as training data and the final 25
timepoints as testing data. BINN models consist of a data-approximating MLP,
$T^{MLP}(x,t)$, and a diffusion-rate-approximating MLP,
$\mathcal{D}^{MLP}(T)$. We train $T^{MLP}$ to closely approximate the ABM
training data while $T^{MLP}$ and $\mathcal{D}^{MLP}$ satisfy Equation (1).
After BINN training, the inferred $\mathcal{D}^{MLP}(T)$ function is used to
forecast and predict ABM data. To forecast ABM training and testing data, we
simulate the diffusion PDE framework using the inferred $\mathcal{D}^{MLP}(T)$
function. To predict ABM data at a new parameter value, $\bm{p}^{new}$, we
perform interpolation over several previously-inferred diffusion rate MLPs,
$\mathcal{D}^{MLP}(T;\bm{p}_{i})$ for $i=1,\dots,K_{1}$, and then simulate the
diffusion PDE framework using the resulting interpolant,
$\mathcal{D}^{interp}(T;\bm{p}^{new})$. Figure 3 gives a visual depiction of
our data analysis pipeline. The python files and notebook used for all steps
of our analysis are presented in
https://github.com/johnnardini/Forecasting_predicting_ABMs.

Figure 3:  Data analysis pipeline. 1. Simulating ABM data For a given
parameter, $\bm{p}$, we simulate the Pulling, Adhesion, or Pulling & Adhesion
ABM. Each model outputs snapshots of pulling and adhesive agent locations over
time; we summarize this data by estimating the average total agent density
along the $x$-direction for each snapshot (not shown). We perform 25 total ABM
simulations for each $\bm{p}$ and average the total spatiotemporal agent
density to obtain $\langle T^{ABM}(x,t;\bm{p})\rangle$. The first 75
timepoints are placed into an training ABM dataset, and the final 25
timepoints are placed into a testing ABM dataset. 2. Training biologically-
informed neural networks (BINNs) to ABM data. Each BINN model consists of a
data-approximating MLP, $T^{MLP}(x,t)$, and a diffusion-rate-approximating
MLP, $\mathcal{D}^{MLP}(T)$. BINN models are trained so that
$T^{MLP}(x,t)\approx\langle T^{ABM}(x,t)\rangle^{train}$ while $T^{MLP}$ and
$\mathcal{D}^{MLP}$ satisfy Equation (1). After model training, the inferred
$\mathcal{D}^{MLP}(T)$ estimates the agent diffusion rate. 3a. Forecasting ABM
data. Simulating the diffusion PDE framework with $\mathcal{D}^{MLP}$ allows
us to forecast the ABM training and testing data. 3b. Predicting new ABM data.
We predict the rate of agent diffusion at a new parameter, $\bm{p}^{new}$, by
interpolating $\mathcal{D}^{MLP}(T;\bm{p})$ over several $\bm{p}$ values to
create $\mathcal{D}^{interp}$. Simulating the diffusion PDE framework with
$\mathcal{D}^{interp}$ allows us to predict the new ABM training and testing
data.

###  3.1 Simulating ABM data

We simulate ABM data by simulating each model over many parameter values,
$\bm{p}$ (Part 1 from Figure 3). For each $\bm{p}$, we simulate $R=25$
identically prepared realizations of the ABM; each realization is completed
when time reaches $t=1000$. We estimate the total spatiotemporal agent density
from each simulation and average over all $R$ simulations to obtain $\langle
T^{ABM}(x,t;\bm{p})\rangle$. We interpolate the output data over time to the
equispaced time grid $t_{j}=(j-1)\Delta t$ with $\Delta t=10$ for
$i=1,\dots,100$ to obtain $\langle T^{ABM}(x,t)\rangle=\\{\langle
T^{ABM}(x_{i},t_{j})\rangle\\}_{i=1,\dots,X}^{j=1,\dots,100}$. We split
$\langle T^{ABM}(x,t)\rangle$ into its training and testing datasets by
setting $\langle T^{ABM}(x,t)\rangle^{train}=\\{\langle
T^{ABM}(x_{i},t_{j})\rangle\\}_{i=1,\dots,X}^{j=1,\dots,T_{f}^{train}}$ and
$\langle T^{ABM}(x,t)\rangle^{test}=\\{\langle
T^{ABM}(x_{i},t_{j})\rangle\\}_{i=1,\dots,X}^{j=T_{f}^{train}+1,\dots,T_{f}^{test}}$.
We set $T_{f}^{train}=75$ and $T_{f}^{test}=100$ to place 75% of data into the
training dataset.

###  3.2 Training Biologically-informed Neural Networks (BINNs) to ABM data

We construct BINN models that consist of two sequential multilayer perceptron
(MLP) models: $T^{MLP}(x,t)$ predicts the total agent density at the point
$(x,t)$, and $\mathcal{D}^{MLP}(T)$ predicts the agent diffusion rate at the
density value $T$ (Part 2 of Figure 3). We train these two MLPs so that
$T^{MLP}(x,t)\approx\langle T^{ABM}(x,t)\rangle^{train}$ while the two MLPs
also satisfy Equation (1) in one spatial dimension:

$\dfrac{\partial}{\partial t}T^{MLP}=\dfrac{\partial}{\partial
x}\left(\mathcal{D}^{MLP}(T^{MLP})\dfrac{\partial}{\partial x}T^{MLP}\right).$
(4)

We chose this nonlinear diffusion PDE framework for BINN model training
because both mean-field models for the Pulling and Adhesion ABMs from Section
2 obey this framework with different diffusion rates.

We further detail the BINNs model architecture in Section 3.2.1, the loss
functions in Section 3.2.2, and our training procedure in Section 3.2.3.

####  3.2.1 BINNs architecture

Following [26], we construct $T^{MLP}(x,t)$ using a fully-connected feed-
forward MLP with three hidden layers, which can be written as:

$\displaystyle z_{0}$ $\displaystyle=[x,t]$ $\displaystyle z_{1}$
$\displaystyle=\sigma\left(z_{0}W_{1}+b_{1}\right)$ $\displaystyle z_{2}$
$\displaystyle=\sigma\left(z_{1}W_{2}+b_{2}\right)$ $\displaystyle z_{3}$
$\displaystyle=\sigma\left(z_{2}W_{3}+b_{3}\right)$ $\displaystyle
T^{MLP}(x,t)$ $\displaystyle=\psi\left(z_{3}W_{4}+b_{4}\right),$ (5)

where each $z_{k}$ denotes the $k^{\text{th}}$ hidden layer; the $W_{k}$
matrices and the $b_{k}$ vectors provide the weights and biases of each hidden
layer, respectively; $\sigma$ denotes the sigmoid activation function
$\sigma(x)=1/(1+\exp{(-x)})$, and $\psi$ denotes the softplus activation
function $\psi(x)=\log(1+\exp(x))$. Each hidden layer in Equation (5) has 128
neurons, meaning that $W_{1}\in\mathbb{R}^{2\times
128};W_{2},W_{3}\in\mathbb{R}^{128\times 128};W_{4}\in\mathbb{R}^{128\times
1};b_{1},b_{2},b_{3}\in\mathbb{R}^{128};\text{ and }b_{4}\in\mathbb{R}$.

The architecture of $\mathcal{D}^{MLP}(T)$ is identical to the architecture
for $T^{MLP}$ in Equation (5), except $\mathcal{D}^{MLP}$ has a one-
dimensional input vector, $T$, instead of the two-dimensional input vector,
$[x,t]$.

####  3.2.2 Loss Function

BINNs are trained to concurrently fit the given dataset, $\langle
T^{ABM}(x,t)\rangle^{train}$, and solve Equation (4) by minimizing the
following multi-term loss function:

$\mathcal{L}_{total}=\mathcal{L}_{WLS}+\mathcal{L}_{PDE}+\mathcal{L}_{constr}.$
(6)

The $\mathcal{L}_{WLS}$ term of Equation (6) computes a weighted mean-squared
error between $T^{MLP}(x,t)$ and $\langle T^{ABM}(x,t)\rangle^{train}$:

$\mathcal{L}_{WLS}=\dfrac{1}{T_{f}^{train}X}\sum_{i=1,j=1}^{X,T_{f}^{train}}w_{i,j}\bigg{(}T^{MLP}(x_{i},t_{j})-\langle
T^{ABM}(x_{i},t_{j})\rangle\bigg{)}^{2}.$ (7)

We set $w_{i,1}=10.0$ for all values of $i$ and all other $w_{i,j}$ values to
1.0 to ensure that $T^{MLP}$ closely agrees with the ABM’s initial data. By
minimizing Equation (7), we ensure $T^{MLP}(x,t)$ closely approximates
$\langle T^{ABM}(x,t)\rangle^{train}$.

The $\mathcal{L}_{PDE}$ term of Equation (6) quantifies how closely $T^{MLP}$
and $\mathcal{D}^{MLP}$ follow Equation (4). To ensure the MLPs satisfy this
PDE framework throughout the ABM’s entire spatiotemporal domain, we uniformly
sample 10,000 points, $\\{(x_{k},t_{k})\\}_{k=1}^{10,000}$, from
$[0,X]\times[0,750]$. For notational convenience, let
$\hat{T}_{k}=T^{MLP}(x_{k},t_{k})$ and
$\hat{D}_{k}=\mathcal{D}^{MLP}\big{(}T^{MLP}(x_{k},t_{k})\big{)}$. We then
compute the mean-squared error between the left- and right-hand sides of
Equation (4) at all sampled points:

$\mathcal{L}_{PDE}=\dfrac{1}{10,000}\sum_{i=1}^{10,000}\bigg{[}\dfrac{\partial}{\partial
t}\hat{T}_{k}-\dfrac{\partial}{\partial
x}\bigg{(}\hat{D}_{k}\dfrac{\partial}{\partial
x}\hat{T}_{k}\bigg{)}\bigg{]}^{2},$ (8)

where differentiation of $T^{MLP}$ and $\mathcal{D}^{MLP}$ is performed using
automatic differentiation. Minimizing Equation (8) verifies that $T^{MLP}$ and
$\mathcal{D}^{MLP}$ together satisfy Equation (4).

The $\mathcal{L}_{constr}$ term of Equation (6) incorporates user knowledge
into BINNs training. We penalize $\mathcal{D}^{MLP}$ for outputting values
outside of the interval $[D_{\min},D_{\max}]$. We set $D_{\min}=0$ because
Equation (4) is ill-posed if $\mathcal{D}(u)<0$, and we set $D_{\max}=1.0$
because the mean-field rates of diffusion are below one for all ABM
simulations in this study. We compute this term by squaring any values of
$\hat{D_{i}}$ that are not within $[D_{\min},D_{\max}]$ and weighting these
values by $10^{10}$:

$\mathcal{L}_{constr}=\dfrac{1}{10,000}\sum_{\begin{subarray}{c}k=1\\\
\hat{D}_{k}\notin[D_{\min},D_{\max}]\end{subarray}}^{10,000}10^{10}(\hat{D_{k}})^{2}.$
(9)

This term regularizes the BINN training procedure to prevent
$\mathcal{D}^{MLP}$ from outputting unrealistic values.

####  3.2.3 BINN Training Procedure

For BINN model training, we randomly partiton the training ABM dataset into
80%/20% BINN training and BINN validation datasets. We train the BINN
parameter values (i.e., the weights and biases for $T^{MLP}$ and
$\mathcal{D}^{MLP}$) to minimize a loss function, $\mathcal{L}$, using the
gradient-based ADAM optimizer with its default hyperparameter values on the
BINN training dataset. For each new set of BINN parameters, we compute
$\mathcal{L}$ on the BINN validation dataset and save the BINN parameters if
the newly computed $\mathcal{L}$ value achieves a 1% or greater relative
improvement over the previous smallest recorded value. Following [34], we
perform training in a two-step process: in the first step, we train the BINN
to match the ABM data by optimizing $\mathcal{L}=\mathcal{L}_{WLS}$ from
Equation (7); in the second step, we train the BINN on
$\mathcal{L}=\mathcal{L}_{total}$ from Equation (6). The first training step
is performed for $10^{4}$ epochs with an early stopping criterion of $10^{3}$,
meaning that training ends early if the smallest-computed $\mathcal{L}$ value
on the validation data is unchanged for $10^{3}$ epochs. The second step is
performed for $10^{6}$ epochs with an early stopping criterion of $10^{5}$.
Each epoch is computed in minibatches of size $10^{3}$. BINN model training is
performed using the PyTorch deep learning library (version 1.7.1).

Following [26], we train five separate BINNs for each ABM dataset using
different BINN training and validation datasets because the final trained
model can be sensitive to which data is included in these two datasets. We
compute the five PDE forward simulations from these trained models and select
whichever BINN achieves the smallest mean-squared error against the ABM
training data as the final selected BINN model.

###  3.3 Forecasting ABM data

We use mean-field PDEs and BINN-guided PDEs to forecast ABM data (Part 3a of
Figure 3). Most of these PDEs are one-compartment nonlinear diffusion
equations that can be written in one spatial dimension as

$\displaystyle\dfrac{\partial T}{\partial t}$
$\displaystyle=\dfrac{\partial}{\partial x}\left(\mathcal{D}(T)\dfrac{\partial
T}{\partial x}\right),$ (10)

where $T=T(x,t)$ is the total agent density and $\mathcal{D}(T)$ is a density-
dependent rate of diffusion. Recall that for the Pulling ABM, $T(x,t)=P(x,t)$;
for the Adhesion ABM, $T(x,t)=H(x,t)$; and for the Pulling & Adhesion ABM,
$T(x,t)=P(x,t)+H(x,t)$.

For the mean-field PDE models, we simulate Equation (10) with
$\mathcal{D}(P)=\mathcal{D}^{pull}(P)$ from Equation (2) for the Pulling ABM,
and $\mathcal{D}(H)=\mathcal{D}^{adh}(H)$ from Equation (3) for the Adhesion
ABM. The mean-field PDE for the Pulling & Adhesion ABM is given by the two-
compartment PDE in Supplementary Equation (28). For BINN-guided PDE models, we
train a BINN model to $\langle T^{ABM}(x,t)\rangle^{train}$ (see Section 3.2)
and then simulate Equation (10) where $\mathcal{D}(T)$ is the
$\mathcal{D}^{MLP}(T)$ that results from BINN model training. Our
implementation method to numerically integrate Equation (10) is provided in
Appendix B.

We partition each PDE simulation,
$T(x,t)=\\{T(x_{i},t_{j})\\}_{i=1,\dots,X}^{j=1,\dots,100}$, into training and
testing datasets to match the training and testing ABM datasets:

$T(x,t)^{train}=\\{T(x_{i},t_{j})\\}_{i=1,\dots,X}^{j=1,\dots,T_{f}^{train}},\
\ \
T(x,t)^{test}=\\{T(x_{i},t_{j})\\}_{i=1,\dots,X}^{j=T_{f}^{train}+1,\dots,T_{f}^{train}}.$

We report the training mean-squared error (MSE) from each simulation as:

$\dfrac{1}{XT_{f}^{train}}\sum_{i=1,j=1}^{X,T_{f}^{train}}\left(T(x_{i},t_{j})-\langle
T^{ABM}(x_{i},t_{j}\rangle\right)^{2},$

and the testing MSE as:

$\dfrac{1}{X(T_{f}^{test}-T_{f}^{train})}\sum_{i=1,j=T_{f}^{train}+1}^{X,T_{f}^{test}}\left(T(x_{i},t_{j})-\langle
T^{ABM}(x_{i},t_{j}\rangle\right)^{2}.$

###  3.4 Predicting new ABM data

We use multivariate interpolation on previously-computed BINN-guided diffusion
rates to predict density-dependent diffusion rates for new ABM data (Part 3b
of Figure 3). We define a prior parameter collection and a new parameter
collection as

$\mathcal{P}^{prior}=\\{\bm{p}_{k}\\}_{k=1}^{K_{1}}\text{ and
}\mathcal{P}^{new}=\\{\bm{p}^{new}_{k}\\}_{k=1}^{K_{2}}.$

Our pipeline for predicting new ABM data from prior ABM simulations proceeds
as follows:

  1. 1.

Generate the prior and new ABM data collections by simulating the ABM at all
parameters from the prior and new parameter collections:

$\mathcal{T}^{prior}=\bigg{\\{}\langle
T^{ABM}(x,t;\bm{p}_{k})\rangle\bigg{\\}}_{k=1}^{K_{1}}\text{ and
}\mathcal{T}^{new}=\bigg{\\{}\langle
T^{ABM}(x,t;\bm{p}^{new}_{k})\rangle\bigg{\\}}_{k=1}^{K_{2}}.$

  2. 2.

Train a BINN model to each training ABM dataset from $\mathcal{T}^{prior}$ and
extract $\mathcal{D}^{MLP}(T;\bm{p}_{k})$ from the trained BINN model.

  3. 3.

Perform multivariate interpolation on
$\\{\mathcal{D}^{MLP}(T;\bm{p}_{k})\\}_{k=1}^{K_{1}}$ to create an
interpolant, $\mathcal{D}^{interp}(T;\bm{p})$, that matches the concatenated
vector $[T,\bm{p}_{k}]$ to the diffusion rate
$\mathcal{D}^{MLP}(T;\bm{p}_{k})$ for $k=1,\dots,K_{1}$.

  4. 4.

Predict the new ABM dataset at $\bm{p}^{new}_{k}$ by simulating Equation (10)
with $\mathcal{D}=\mathcal{D}^{interp}(T;\bm{p}^{new}_{k})$ to create
$T^{interp}(x,t;\bm{p}^{new}_{k})$. Partition
$T^{interp}(x,t;\bm{p}^{new}_{k})$ into its training and testing datasets to
match the ABM data’s training and testing datasets.

  5. 5.

Compute the training and testing MSEs between
$T^{interp}(x,t;\bm{p}^{new}_{k})$ and $\langle T^{ABM}(x,t)\rangle$ to
summarize the predictive performance of $T^{interp}(x,t;\bm{p}^{new}_{k})$ for
$k=1,\dots,K_{2}$.

We implement multi-dimensional radial basis function interpolation using Sci-
kit Learn’s (version 0.24.2) RBFInterpolator command to create
$\mathcal{D}^{interp}(T;\bm{p})$.

##  4 Results

###  4.1 PDE simulations outperform neural networks in ABM forecasting

We investigate the performance of an ANN, BINN, BINN-guided PDE model, and
mean-field PDE model in forecasting ABM data. We simulated the Pulling ABM
with $\bm{p}=(r_{m}^{pull},p_{pull})^{T}=(1.0,0.5)^{T}$ to generate the ABM
data. The ANN was trained to minimize the loss function $\mathcal{L}_{WLS}$
from Equation (7) on the training ABM dataset, whereas the BINN was trained to
minimize $\mathcal{L}_{total}$ from Equation (6). Both PDE models simulate
Equation (10). For the BINN-guided PDE, we extract $\mathcal{D}^{MLP}$ from
the trained BINN model to obtain $\mathcal{D}$; for the mean-field PDE,
$\mathcal{D}$ is given by $\mathcal{D}^{pull}$ from Equation (2).

Visual inspection suggests that all four models match the ABM training data
well (Figure 4(a-b)111the mean-field PDE is not plotted in this figure because
it is visually indistinguishable from the BINN-guided PDE.), though the
computed training MSE values reveal that the mean-field and BINN-guided PDEs
outperform the neural networks in describing this data (Table 2). The BINN,
BINN-guided PDE, and mean-field PDE all accurately forecast the testing data
(Figure 4(c)), but the two PDE models achieve smaller MSE values than the BINN
model (Table 2). The ANN’s prediction for the testing data has a protrusion
that overpredicts all data for $x>125$ (Figure 4(c) inset), which causes this
model’s computed testing MSE value to be almost an order of magnitude higher
than all others.

Figure 4: Forecasting ABM data with neural networks and PDEs. ANN and BINN models were trained to fit $\langle T^{ABM}(x,t)\rangle^{train}$ from the Pulling ABM with $\bm{p}=(r_{m}^{pull},p_{pull})^{T}=(1.0,0.5)^{T}$. These two ANNs and the mean-field and BINN-guided (BG) PDE simulations were then used to forecast (a-b) $\langle T^{ABM}(x,t)\rangle^{train}$ and (c) $\langle T^{ABM}(x,t)\rangle^{test}$. The mean-field PDE simulation is not plotted because it is visually indistinguishable from the BG PDE simulation. Model | Training MSE | Testing MSE
---|---|---
Artificial neural network | $1.17\times 10^{-4}$ | $9.36\times 10^{-4}$
Biologically-informed neural network | $9.32\times 10^{-5}$ | $1.47\times 10^{-4}$
Mean-field PDE | $7.45\times 10^{-5}$ | $1.00\times 10^{-4}$
BINN-guided PDE | $7.64\times 10^{-5}$ | $1.02\times 10^{-4}$
Table 2: Computed MSE values when forecasting $\langle
T^{ABM}(x,t)\rangle^{train}$ and $\langle T^{ABM}(x,t)\rangle^{train}$ with an
ANN, BINN, mean-field PDE, or BINN-guided PDE.

###  4.2 Forecasting ABM data with BINN-guided and mean-field PDE simulations

We investigate the performance of BINN-guided and mean-field PDE simulations
in forecasting the training and testing ABM datasets from the Pulling,
Adhesion, and Pulling & Adhesion ABMs. See Section 3.3 for implementation
details.

####  4.2.1 The BINN-guided and mean-field PDEs both accurately forecast
Pulling ABM data

The parameters for the Pulling ABM are $\bm{p}=(r_{m}^{pull},p_{pull})^{T}$.
To evaluate the BINN-guided and mean-field PDE models’ performances in
forecasting Pulling ABM data over a range of agent pulling parameter values,
we computed eleven ABM datasets by varying $p_{pull}=0.0,0.1,0.2,\dots,1.0$
while fixing $r_{m}^{pull}=1.0$. The inferred rates of agent diffusion from
both models propose that agents diffuse slower for low densities and faster
for high densities, and that larger values of $p_{pull}$ lead to increased
density-dependent diffusion rates (Figure 5(a)). The two models achieve
comparable training and testing MSE values for all values of $p_{pull}$,
though the mean-field PDE usually attains slightly smaller values (Figure
5(b)). Snapshots of both simulated PDE models against data shows that their
ABM predictions are visually indistinguishable (Supplementary Figure 12(a-c)).

Figure 5: Forecasting Pulling ABM data with the mean-field (MF) and BINN-
guided PDEs. (a) Plots of the mean-field diffusion rate,
$\mathcal{D}^{pull}(T)$, from Equation (2) and the inferred BINN diffusion
rate, $\mathcal{D}^{MLP}(T)$, for $p_{pull}=0.1,0.3,\dots,0.9$ (results not
shown for $p_{pull}=0.0,0.2,\dots,1.0$ for visual ease) while fixing
$r_{m}^{pull}=1.0$. (b) Plots of the mean-field and BINN-guided PDEs’ computed
training and testing MSE values while varying $p_{pull}$ and fixing
$r_{m}^{pull}=1.0$. (c) Plots of $\mathcal{D}^{pull}(T)$ and
$\mathcal{D}^{MLP}(T)$ for $r_{m}^{pull}=0.2,0.4,\dots,1.0$ while fixing
$p_{pull}=0.5$. (d) Plots of the mean-field and BINN-guided PDEs’ computed
training and testing MSE values while varying $r_{m}^{pull}$ and fixing
$p_{pull}=0.5$.

To evaluate both PDE models’ performances over a range of pulling agent
migration values, we computed ten Pulling ABM datasets with
$r_{m}^{pull}=0.1,0.2,\dots,1.0$ while fixing $p_{pull}=0.5$. We find close
agreement between both models’ inferred diffusion rates for values of
$r_{m}^{pull}$ (Figure 5(c)). As a result, both models achieve similar
computed training and testing MSE values (Figure 5(d)). Snapshots of both
simulated PDE models against data reveals that their ABM predictions are
visually indistinguishable (Supplementary Figure 12(d-f)).

####  4.2.2 BINN-guided PDEs accurately forecast Adhesion ABM data when the
mean-field PDE is ill-posed

The parameters for the pulling ABM are $\bm{p}=(r_{m}^{adh},p_{adh})^{T}$. To
evaluate the BINN-guided and mean-field PDE models’ performances over a range
of agent adhesion parameter values, we computed eleven ABM datasets by varying
$p_{adh}=0.0,0.1,0.2,\dots,1.0$ while fixing $r_{m}^{adh}=1.0$. The inferred
rates of agent diffusion from both models decrease with agent density for most
values of $p_{adh}$ (Figure 6(a)). When $p_{adh}=0$, BINN-guided diffusion
rate is slightly increasing and the mean-field model’s diffusion rate is
constant. The BINN-guided diffusion rates decline faster with agent density
than the corresponding mean-field diffusion rates for densities below 0.4.
Both models agree that the density-dependent rates of diffusion fall as
$p_{adh}$ increases. We computed the training and testing MSEs for both models
for all values of $p_{adh}$ (Figure 6(b)) and partition the results as follows
:

  * •

When $\bm{p_{adh}<0.5}$: both models achieve similar training MSE values near
$7\times 10^{-5}$ and testing MSE values around $10^{-4}$.

  * •

When $\bm{0.5\leq p_{adh}\leq 0.75}$: the mean-field PDE model’s training and
testing MSE values increase as $p_{adh}$ increases, with a maximum computed
value above $3\times 10^{-4}$. The BINN-guided PDE model’s training and
testing MSE values remain near $7\times 10^{-5}$ and $10^{-4}$, respectively.

  * •

When $\bm{p_{adh}>0.75}$: the mean-field PDE model is ill-posed and cannot
forecast this ABM data. The BINN-guided PDE model’s computed training and
testing MSE values increase as $p_{adh}$ increases, with a maximum computed
value of $2\times 10^{-4}$.

Close inspection of snapshots from both PDE model simulations against ABM data
from $p_{adh}=0.7$ reveals that the mean-field PDE model slightly overpredicts
the data at high densities above 0.5 and low densities below 0.1, whereas the
BINN-guided PDE closely matches the data (Supplementary Figure 13(c) inset).

Figure 6: Forecasting Adhesion ABM data with the mean-field (MF) and BINN-
guided PDEs. (a) Plots of the mean-field diffusion rate,
$\mathcal{D}^{adh}(T)$, from Equation (3) and the inferred BINN diffusion
rate, $\mathcal{D}^{MLP}(T)$, for $p_{adh}=0.1,0.3,\dots,0.9$ (results not
shown for $p_{adh}=0.0,0.2,\dots,1.0$ for visual ease) while fixing
$r_{m}^{adh}=1.0$. (b) Plots of the mean-field and BINN-guided PDEs’ computed
training and testing MSE values while varying $p_{adh}$ and fixing
$r_{m}^{adh}=1.0$. (c) Plots of $\mathcal{D}^{adh}(T)$ and
$\mathcal{D}^{MLP}(T)$ for $r_{m}^{adh}=0.2,0.4,\dots,1.0$ while fixing
$p_{adh}=0.5$. (d) Plots of the mean-field and BINN-guided PDEs’ computed
training and testing MSE values while varying $r_{m}^{adh}$ and fixing
$p_{adh}=0.5$.

To evaluate both PDE models’ performances over a range of adhesive agent
migration values, we computed ten ABM datasets with
$r_{m}^{adh}=0.1,0.2,\dots,1.0$ while fixing $p_{adh}=0.5$. Both PDE models
similarly propose that agent density-dependent diffusion rates decrease for
larger agent density values and that these rates increase for larger values of
$r_{m}^{adh}$ (Figure 6(c)). As a result, both PDEs achieve similar computed
training and testing MSE values for most values of $r_{m}^{adh}$ (Figure
6(d)). When $r_{m}^{adh}=0.1$, however, the BINN-guided PDE’s testing MSE
value is close to $10^{-4}$, whereas the mean-field PDE attains a much lower
testing MSE value near $6\times 10^{-5}$. Despite these differences, the two
model simulations appear similar at these parameter values (Supplementary
Figure 13(d-f)).

####  4.2.3 BINN-guided PDEs accurately forecast Pulling & Adhesion ABM data
with a one-compartment model

The parameters for the Pulling & Adhesion ABM are
$\bm{p}=(r_{m}^{pull},r_{m}^{adh},p_{pull},p_{adh},\alpha)^{T}$. We evaluate
the performance of the BINN-guided and mean-field DE models in forecasting
data from the Pulling & Adhesion ABM. We created 48 ABM datasets by fixing the
base parameter values at $\bm{p}_{base}=(1.0,0.25,0.33,0.33,0.5)^{T}$ and then
varying one parameter at a time over several values. We vary
$r_{m}^{pull}=0.5,0.6,\dots,1.5$; $r_{m}^{adh}=0.0,0.1,\dots,1.0$;
$p_{pull}=0.1,0.2,\dots,0.6,0.67$; $p_{adh}=0.1,0.2,\dots,0.6,0.67$; and
$\alpha=0.0,0.1,\dots,1.0$. These parameter values were chosen to always
satisfy $p_{pull}+p_{adh}\leq 1.$

The BINN models’ inferred diffusion rates are often U-shaped with large
diffusion values at low and high agent densities and smaller values at
intermediate densities (Figure 7). This U-shape tends to increase for larger
values of $r_{m}^{pull},r_{m}^{adh},\text{ and }p_{pull}$ and decrease for
larger values of $p_{adh}$ and $\alpha$. The inferred diffusion rates appear
most sensitive to changes in the $\alpha$ parameter: at $\alpha=0.0$, it
strictly increases with agent density and attains an average value of 0.289;
at $\alpha=1.0$, it is strictly decreasing and has an average value of 0.051.
The inferred diffusion rate is also sensitive to the $r_{m}^{adh}$ and
$r_{m}^{pull}$ parameters: varying $r_{m}^{adh}$ primarily alters the BINN
diffusion rate at intermediate agent density values, whereas varying
$r_{m}^{pull}$ changes the BINN diffusion rate at low and high agent densitiy
values.

Figure 7: The inferred BINN diffusion rates for Pulling & Adhesion ABM data.
Plots of the inferred BINN diffusion rate, $\mathcal{D}^{MLP}(T)$, when
varying (a) $r_{m}^{pull}$, (b) $r_{m}^{adh}$, (c) $p_{pull}$, (d) $p_{adh}$,
(e) $\alpha$.

Recall that the BINN-guided PDE computes a single compartment to forecast the
total agent density, $T(x,t)$, whereas the mean-field PDE computes two
compartments forecasting the Pulling and Adhesive agent densities, $P(x,t)$
and $H(x,t)$, respectively. We forecast the total agent density with the mean-
field PDE by setting $T(x,t)=P(x,t)+H(x,t)$. The BINN-guided and mean-field
PDE models achieve similar training MSE values for most parameter values that
we considered (Figure 8). The mean-field model’s testing MSE values are often
smaller than the BINN-guided testing MSE values, though the BINN-guided PDE
also achieves small testing MSE values. For example, both PDE simulations
accurately predict ABM data when $p_{adh}$ is set to $0.4$, but visualizing
both PDE simulations shows that the mean-field PDE better matches the elbow of
the data than the BINN-guided PDE (Supplementary Figure 14(a-c)). The BINN-
guided PDE outperforms the mean-field PDE in forecasting data for small values
of $r_{m}^{adh}$: plotting both PDE simulations against data from
$r_{m}^{adh}=0.1$ shows that the mean-field PDE underpredicts the largest
agent density values, while the BINN-guided PDE accurately matches this data
(Supplementary Figure 14(d-f)).

Figure 8: Forecasting Pulling & Adhesion ABM data with the mean-field and
BINN-guided PDEs. Plots of the mean-field and BINN-guided PDEs’ computed
training and testing values while varying (a) $r_{m}^{pull}$, (b)
$r_{m}^{adh}$, (c) $p_{pull}$, (d) $p_{adh}$, (e) $\alpha$.

###  4.3 Predicting ABM data at new parameter values

ABM simulations can be computationally expensive when the model includes
complex rules or consists of many agents. This computational bottleneck makes
it challenging to investigate ABM behavior at many parameter values. We now
examine how performing multivariate interpolation on several BINN-inferred
diffusion rates can aid the prediction of previously-unseen ABM data at new
parameter values (see Section 3.4 for implementation details).

We predict new data from the Adhesion and Pulling & Adhesion ABMs in this
section. We do not include the Pulling ABM in this work because the mean-field
PDE model accurately forecasted ABM data for all parameter values that we
considered in Section 4.2.1.

####  4.3.1 Predicting Adhesion ABM data

The parameters for the Adhesion ABM are $\bm{p}=(r_{m}^{adh},p_{adh})^{T}$. We
perform ABM data prediction for $p_{adh}\geq 0.5$ in this section because we
found that the mean-field PDE model accurately forecasted ABM data for
$p_{adh}\leq 0.5$ in Section 4.2.2.

We first perform interpolation over the $p_{adh}$ parameter while fixing
$r_{m}^{adh}$. The prior data collection consists of six ABM datasets
generated by varying $p_{adh}=0.5,0.6,0.7,\dots,1.0$ while fixing
$r_{m}^{adh}=1.0$; the new data collection consists of five ABM datasets
generated by varying $p_{adh}=0.55,0.65,0.75,0.85,\text{ and }0.95$ while
fixing $r_{m}^{adh}=1.0$. We performed multivariate interpolation over the six
inferred $\mathcal{D}^{MLP}(T;\bm{p})$ terms from the prior data collection to
generate $\mathcal{D}^{interp}(T;\bm{p})$. We use this interpolant to predict
the diffusion rates for all parameters from the new data collection (9(a)).
All inferred diffusion rates decrease with agent density and tends to fall
with larger $p_{adh}$ values. Most of the computed training and testing MSE
values on the new data collection are comparable to their counterparts from
the prior data collection, except the testing MSE at $p_{adh}=0.95$ exceeds
$5\times 10^{-4}$ while the testing MSEs at $p_{adh}=0.9\text{ and }1.0$ do
not exceed $2.5\times 10^{-4}$ (Figure 9(b)). Visual inspection of the
simulated PDE prediction against ABM data at $p_{adh}=0.95$ reveals that it
matches the data well but slightly mispredicts the data’s heel at later time
points (Supplementary Figure 15(a-c)).

Figure 9: Predicting Adhesion ABM data with BINN-guided PDEs and multivariate
interpolation for new $p_{adh}$ values. The parameters for the Adhesion ABM
are given by $\bm{p}=(r_{m}^{adh},p_{adh})^{T}.$ Here, we vary $p_{adh}$ while
fixing $r_{m}^{adh}=1.0$. The prior data collection consists of
$p_{adh}=0.5,0.6,\dots,1.0$ and the new data collection consists of
$p_{adh}=0.55,0.65,\dots,0.95$ (a) Plots of the learned
$\mathcal{D}^{MLP}(T;\bm{p})$ diffusion rates for the prior data collection.
We performed multivariate interpolation on these rates to obtain
$\mathcal{D}^{interp}(T;\bm{p})$, which we plot for the new data collection.
(b) Plots of the BINN-guided PDEs’ computed training and testing values on the
prior data collection, and the interpolated PDE’s training and testing values
on the new data collection.

We next perform interpolation over the $r_{m}^{adh}$ and $p_{adh}$ parameters.
The prior data collection consists of 18 ABM datasets generated by varying
$r_{m}^{adh}=0.1,0.5,1.0$ and $p_{adh}=0.5,0.6,\dots,1.0$; the new data
collection consists of ten ABM datasets generated from a latin hypercube
sampling of $(r_{m}^{adh},p_{adh})\in[0.1,1.0]\times[0.5,1.0]$ (Figure 10(a)
and Supplementary Table 4). We performed multivariate interpolation over each
$\mathcal{D}^{MLP}(T;\bm{p})$ from the prior data collection to generate
$\mathcal{D}^{interp}(T;\bm{p})$. The predicted diffusion rates for the new
data collection decrease with agent density, rise for larger $r_{m}^{adh}$
values, and decrease faster for larger $p_{adh}$ values (Figure 10(b)). We
order the parameters from the new data collection by increasing training MSE
values; four of the five lowest training MSE values result from the five
smallest $p_{adh}$ values, and four of the five highest MSE values result from
the five highest $p_{adh}$ values (Figure 10(c)). The four lowest training and
testing MSE values are all below $110^{-4}$, the eight lowest are all below
$2\times 10^{-4}$, and the highest testing MSE value reaches $1.6\times
10^{-3}$. Visual inspection of the simulated PDE prediction with the highest
testing MSE value reveals that this simulation mispredicts the data’s heel but
otherwise matches the ABM data well (Supplementary Figure 16(a-c)). Visual
inspection of the simulated PDE prediction with the third-highest MSE value
shows that this simulations accurately matches the ABM data (Supplementary
Figure 16(d-f)).

Figure 10: Predicting Adhesion ABM data with BINN-guided PDEs and multivariate
interpolation for new $r_{m}^{adh}$ and $p_{adh}$ values. The parameters for
the Adhesion ABM are given by $\bm{p}=(r_{m}^{adh},p_{adh})^{T}.$ Here, we
vary both parameters. (a) The prior data collection consists of
$r_{m}^{adh}=0.1,0.5,1.0$ and $p_{adh}=0.5,0.6,\dots,1.0$ and the new data
collection consists of a Latin hypercube (LHC) sampling of
$\bm{p}\in[0.1,1.0]\times[0.5,1.0]$ with 10 samples. (b) We performed
multivariate interpolation on the $\mathcal{D}^{MLP}(T;\bm{p})$ rates on the
prior data collection to obtain $\mathcal{D}^{interp}(T;\bm{p})$. We plot
three illustrative $\mathcal{D}^{interp}(T;\bm{p})$ values from the new data
collection. (c) Plots of the interpolated PDE’s training and testing values on
the new data collection.

####  4.3.2 Predicting Adhesion & Pulling ABM data

The parameters for the Pulling & Adhesion ABM are
$\bm{p}=(r_{m}^{pull},r_{m}^{adh},p_{pull},p_{adh},\alpha)^{T}$. We perform
ABM data prediction over a large range of parameter values to determine if the
one-compartment BINN-guided PDE simulations can predict this ABM’s data that
results from two interacting subpopulations.

We perform multivariate interpolation over the $p_{pull},p_{adh},$ and
$\alpha$ parameters while fixing $r_{m}^{pull}=1.0$ and $r_{m}^{adh}=0.25$.
The prior and new data collections consist of 40 and 20 ABM parameter
combinations, respectively, that were generated from Latin hypercube samplings
of $(p_{pull},p_{adh},\alpha)\in[0,0.67]\times[0,0.67]\times[0,1]$ (Figure
11(a) and Supplementary Tables 5 and 6). We chose samplings where
$p_{pull}+p_{adh}\leq 1.0$ for all samples. The computed training and testing
MSE values for the new parameter collection suggest all simulated PDE
predictions accurately match the ABM data at those parameters (Figure 11(b)).
Of the 20 computed testing MSE values in the new data collection, four are
below $1.0\times 10^{-4}$, 16 are below $2.0\times 10^{-4}$, and all are below
$5.0\times 10^{-4}$. The highest and third highest testing MSE value results
from $(p_{pull},p_{adh},\alpha)=(0.218,0.553,0.675)$ and
$(0.251,0.486,0.975)$, respectively. Visually inspecting the simulated PDE
predictions from these parameter values against ABM data reveals that both
match the data well, though the worst prediction overpredicts the highest ABM
density values (Supplementary Figure 17).

Figure 11: Predicting Pulling & Adhesion ABM data with BINN-guided PDEs and
multivariate interpolation for new $p_{pull},p_{adh}$, and $\alpha$ values.
The parameters for the Adhesion ABM are given by
$\bm{p}=(r_{m}^{adh},r_{m}^{pull},p_{adh},p_{pull},\alpha)^{T}.$ Here, we vary
$p_{pull},p_{adh},$ and $\alpha$ while fixing $r_{m}^{pull}=1.0$ and
$r_{m}^{adh}=0.25$. (a) The prior data consists of a Latin hypercube (LHC)
sampling of $(p_{pull},p_{adh},\alpha)\in[0,0.67]\times[0,0.67]\times[0,1]$
with 40 samples and the new data consists of a LHC sampling of the same domain
with 20 samples. (b) Plots of the interpolated PDE’s training and testing
values on the new data.

###  4.4 Comparing the computational expense of each modeling approach

We finish with a discussion on the computational expense of all approaches
discussed in this work (Table 3 and Supplementary Figure 18). We recorded the
computed wall times to simulate each ABM, train each BINN model, and simulate
each PDE in Section 4.2. Averaging across all ABMs suggests that the average
ABM dataset took 40.0 minutes to generate with a standard deviation of 15.6
minutes. The average mean-field PDE model simulations for the Pulling ABM and
the Adhesion ABM took 0.6 and 0.5 seconds to complete, respectively, which are
about 4,000 and 4,500 times faster than the average ABM simulations time. The
average mean-field PDE model simulation time for the Pulling & Adhesion ABM
was4.7 seconds, which is 542 times faster than the average ABM simulation
time. Training a BINN model is the most time-consuming task with an average
time of 11.2 hours across all ABMs with an average standard deviation of 4.32
hours. The average BINN-guided PDE simulation takes 82.9 seconds with a
standard deviation of 77.12 seconds, which is approximately 28 times faster
than simulating the ABM.

ABM Name | ABM simulation | MF PDE simulation | BINN Training | BG PDE simulation
---|---|---|---|---
Adhesion | 37.5 (15.4) minutes | 0.5 (0.15) seconds | 10.6 (4.44) hours | 16.9 (23.65) seconds
Pulling | 39.9 (15.8) minutes | 0.6 (0.20) seconds | 10.0 (3.99) hours | 164.8 (156.9) seconds
Pulling & Adhesion | 42.5 (15.52) minutes | 4.7 (1.20) seconds | 13.1 (4.54) hours | 66.9 (50.81) seconds
Table 3: Computational expenses of each modeling approach. The mean wall time
computations (standard deviation in parentheses) for ABM simulations, BINN
training, mean-field PDE simulations, and BINN-guided PDE simulations for all
three ABMs.

##  5 Discussion and Future Work

The primary aim of this work was to introduce how BINNs can be used to learn
PDE models capable of forecasting and predicting ABM data. We used three
stochastic ABMs that consist of rules on how agent pulling and adhesion impact
agent migration to illustrate this approach. These models capture some of the
key cellular interactions for tumor invasion, wound healing, and development
[3, 12]. It is challenging to predict how the parameters that characterize
these processes will impact the output model behaviors, however, due to the
models’ high computational expenses. Modelers frequently address this
limitation by coarse-graining ABM rules into computationally-efficient PDE
models. Unfortunately, the resulting models can give misleading ABM
predictions and may be ill-posed for certain parameter values [6, 9]. Here, we
demonstrated that training BINNs to ABM data allows us to learn BINN-guided
PDE models that are capable of both 1.) forecasting future ABM data at a fixed
parameter value and 2.) predicting ABM data from new parameter values. We use
BINNs to forecast ABM data by simulating Equation (10) with the trained BINN
model’s inferred diffusion rate. ABM prediction requires previously-computed
ABM data from several parameter values. We train a BINN model to each dataset,
and then perform multivariate interpolation over their inferred diffusion
rates. We predict the output ABM data at new parameters by simulating Equation
(10) with the resulting interpolant as the diffusion rate.

We highlighted the strong forecasting and prediction capabilities of this
approach using multiple simulations from the Pulling, Adhesion, and Pulling &
Adhesion ABMs. For the Pulling ABM, the trained BINNs learn diffusion rates
that are similar to the mean-field diffusion rates for all parameter values;
as such, both models perform similarly in forecasting future ABM data. Due to
the similar computed results for all parameter values, we suggest that the
mean-field PDE model can be used (in lieu of performing multivariate
interpolation over several BINN diffusion rates) to predict new Pulling ABM
data. For the Adhesion ABM, the BINN-guided PDEs accurately forecast ABM data,
even for large adhesion values where the mean-field PDE is ill-posed. We can
perform multivariate interpolation to accurately predict new Adhesion ABM data
for large adhesion values. For the Pulling & Adhesion ABM, both the BINN-
guided and mean-field PDE simulations accurately forecasts the total ABM data.
The BINN-guided PDE achieves this with a one-compartment model, whereas the
mean-field PDE must compute two compartments (one for each agent type). We
were able to perform multivariate interpolation to accurately predict new
Pulling & Adhesion ABM data when varying parameters that alter the proportions
of agent types in the simulation, and the rates of agent pulling and adhesion.

A limitation of our approach for ABM forecasting and prediction is the
computational expense of BINN model training. The average BINN training
procedure in this study took 11.2 hours, which is about 17 times longer than
the average ABM data generation time of 40 minutes. Once a BINN model has been
trained, however, the average BINN-guided PDE simulation took 83 seconds,
which is roughly 28 times faster than the average time to generate an ABM
dataset. One possible source of these long BINN training times is our chosen
BINN model architecture, which consists of over over 50,000 parameters to
train. Kaplarevi-Malii et al. [33] proposed a genetic algorithm to identify
the optimal model archictecture for PINN models. In future work, we plan to
implement this algorithm to identify simpler BINN model architectures that can
be efficiently trained to learn predictive PDE models for ABMs.

Future work may include modeling more complex ABMs using the methods proposed
in this study. We focused on simple pulling rules that involve two agents, but
Chappelle and Yates [7] considered rules where multiple agents interact during
each pulling event. They found that the predictive accuracy of the coarse-
grained PDE models decreases as more agents become involved in pulling
processes. VandenHeuvel et al. [49] studied ABMs involving mechanical
relaxation to mimic epithelial tissues, but the coarse-grained PDEs are only
valid for fast mechanical relaxation rates. The coarse-grained PDE models for
the migration rules from these ABMs obey Equation (10) in one spatial
dimension; we could use BINNs to learn PDEs that better forecast and predict
ABM data for these challenging model rules. Agent proliferation is an
important component of biological invasion in wound healing and cancer growth
that we did not consider in this work. Many previous studies have shown that
coarse-grained PDE models fail to accurately describe ABM simulations when
agent proliferate quickly [9]. We could easily extend our methods to predict
ABM data from models where agents migrate and proliferate by adding a
population growth term into Equation (10) during BINN training.

An additional area of future work includes further BINN model development as a
means to deepen our ABM analysis. We studied population heterogeneity in this
work using an ABM with two agent types. We used a one-compartment PDE during
BINN training to demonstrate BINNs’ ability to forecast and predict complex
data with simple PDE models. Alternatively, we could update the BINN modeling
framework to predict the density of each agent type (as well as their
diffusion and/or growth terms) with a multi-compartment PDE framework. We
could infer how model parameters affect the modeling terms for each type of
agent. Incorporating the Bayesian inference techniques introduced in [34]
would allow us to develop Bayesian BINNs capable of performing uncertainty
quantification for each modeling term. There are thus many ways to extend the
methods proposed in this work, which will allow us to study more complex ABMs
and improve the learned models’ predictive capabilities.

Data Availability statement: All code for this work is publicly available at
https://github.com/johnnardini/Forecasting_predicting_ABMs.

## References

  * [1] Michalina Janiszewska, Marina C. Primi, and Tina Izard.  Cell adhesion in cancer: beyond the migration of single cells.  Journal of Biological Chemistry, 295(8):2495–2505, February 2020\.
  * [2] Katheryn E. Rothenberg, Yujun Chen, Jocelyn A. McDonald, and Rodrigo Fernandez-Gonzalez.  Rap1 coordinates cell-cell adhesion and cytoskeletal reorganization to drive collective cell migration in vivo.  Current Biology, 33(13):2587–2601.e5, July 2023.
  * [3] John T. Nardini, Douglas A. Chapnick, Xuedong Liu, and David M. Bortz.  Modeling keratinocyte wound healing: cell-cell adhesions promote sustained migration.  Journal of Theoretical Biology, 400:103–117, July 2016.
  * [4] Peter Friedl and Darren Gilmour.  Collective cell migration in morphogenesis, regeneration and cancer.  Nature Reviews Molecular Cell Biology, 10(7):445–457, July 2009\.
  * [5] R. A. F. Clark and P. Henson.  The molecular and cellular biology of wound repair.  Plenum Press, New York, second edition, 1995.
  * [6] Kathleen Anguige and Christian Schmeiser.  A one-dimensional model of cell diffusion and aggregation, incorporating volume filling and cell-to-cell adhesion.  Journal of Mathematical Biology, 58(3):395, March 2009.
  * [7] George Chappelle and Christian A. Yates.  Pulling in models of cell migration.  Physical Review E, 99(6):062413, June 2019.
  * [8] Shahzeb Raja Noureen, Jennifer P. Owen, Richard L. Mort, and Christian A. Yates.  Swapping in lattice-based cell migration models.  Physical Review E, 107(4):044402, April 2023.
  * [9] Ruth E. Baker and Matthew J. Simpson.  Correcting mean-field approximations for birth-death-movement processes.  Physical Review E, 82(4):041905, October 2010.
  * [10] Matthew J. Simpson, Ruth E. Baker, Pascal R. Buenzli, Ruanui Nicholson, and Oliver J. Maclaren.  Reliable and efficient parameter estimation using approximate continuum limit descriptions of stochastic models.  Journal of Theoretical Biology, 549:111201, September 2022.
  * [11] Robert J. H. Ross, Christian. A. Yates, and Ruth E. Baker.  Inference of cell–cell interactions from population density characteristics and cell trajectories on static and growing domains.  Mathematical Biosciences, 264:108–118, June 2015.
  * [12] Robin N. Thompson, Christian A. Yates, and Ruth E. Baker.  Modelling cell migration and adhesion during development.  Bulletin of Mathematical Biology, 74(12):2793–2809, December 2012\.
  * [13] John T. Nardini, Ruth E. Baker, Matthew J. Simpson, and Kevin B. Flores.  Learning differential equation models from stochastic agent-based model simulations.  Journal of The Royal Society Interface, 18(176):rsif.2020.0987, 20200987, March 2021.
  * [14] Le-Minh Kieu, Nicolas Malleson, and Alison Heppenstall.  Dealing with uncertainty in agent-based models for short-term predictions.  Royal Society Open Science, 7(1):191074, January 2020.
  * [15] Dale Larie, Gary An, and R. Chase Cockrell.  The use of artificial neural networks to forecast the behavior of agent-based models of pathophysiology: an example utilizing an agent-based model of sepsis.  Frontiers in Physiology, 12:716434, October 2021.
  * [16] Matthew J. Simpson, Jesse A. Sharp, and Ruth E. Baker.  Distinguishing between mean-field, moment dynamics and stochastic descriptions of birth–death–movement processes.  Physica A: Statistical Mechanics and its Applications, 395:236–246, February 2014.
  * [17] Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Discovering governing equations from data by sparse identification of nonlinear dynamical systems.  Proceedings of the National Academy of Sciences, 113(15):3932–3937, April 2016.
  * [18] E. Kaiser, J. Nathan Kutz, and Steven L. Brunton.  Sparse identification of nonlinear dynamics for model predictive control in the low-data limit.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 474(2219):20180335, November 2018.
  * [19] Samuel Rudy, Alessandro Alla, Steven L. Brunton, and J. Nathan Kutz.  Data-driven identification of parametric partial differential equations.  SIAM Journal on Applied Dynamical Systems, 18(2):643–660, January 2019.  Publisher: Society for Industrial and Applied Mathematics.
  * [20] Kathleen Champion, Bethany Lusch, J. Nathan Kutz, and Steven L. Brunton.  Data-driven discovery of coordinates and governing equations.  Proceedings of the National Academy of Sciences, 116(45):22445–22451, November 2019.  Publisher: National Academy of Sciences Section: Physical Sciences.
  * [21] Niall M. Mangan, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Inferring biological networks by sparse identification of nonlinear dynamics.  IEEE Transactions on Molecular, Biological and Multi-Scale Communications, 2(1):52–63, June 2016.  Conference Name: IEEE Transactions on Molecular, Biological and Multi-Scale Communications.
  * [22] Niall M. Mangan, J. Nathan Kutz, Steven L. Brunton, and Joshua L. Proctor.  Model selection for dynamical systems via sparse regression and information criteria.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2204):20170009, August 2017.
  * [23] Daniel A. Messenger and David M. Bortz.  Weak SINDy: galerkin-based data-driven model selection.  Multiscale Modeling & Simulation, 19(3):1474–1497, January 2021\.
  * [24] Daniel A. Messenger and David M. Bortz.  Weak SINDy for partial differential equations.  Journal of Computational Physics, 443:110525, October 2021.
  * [25] John H. Lagergren, John T. Nardini, G. Michael Lavigne, Erica M. Rutter, and Kevin B. Flores.  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(2234):20190800, February 2020.
  * [26] John H. Lagergren, John T. Nardini, Ruth E. Baker, Matthew J. Simpson, and Kevin B. Flores.  Biologically-informed neural networks guide mechanistic modeling from sparse experimental data.  PLOS Computational Biology, 16(12):e1008462, December 2020.  Publisher: Public Library of Science.
  * [27] John T. Nardini, John H. Lagergren, Andrea Hawkins-Daarud, Lee Curtin, Bethan Morris, Erica M. Rutter, Kristin R. Swanson, and Kevin B. Flores.  Learning equations from biological data with limited time samples.  Bulletin of Mathematical Biology, 82(9):119, September 2020.
  * [28] Samuel H. Rudy, Steven L. Brunton, Joshua L. Proctor, and J. Nathan Kutz.  Data-driven discovery of partial differential equations.  Science Advances, 3(4):e1602614, April 2017.  Publisher: American Association for the Advancement of Science Section: Research Article.
  * [29] Daniel A. Messenger, Graycen E. Wheeler, Xuedong Liu, and David M. Bortz.  Learning anisotropic interaction rules from individual trajectories in a heterogeneous cellular population.  Journal of The Royal Society Interface, 19(195):20220412, October 2022.
  * [30] Daniel A. Messenger and David M. Bortz.  Learning mean-field equations from particle data using WSINDy.  Physica D: Nonlinear Phenomena, 439:133406, November 2022.
  * [31] Rohit Supekar, Boya Song, Alasdair Hastewell, Gary P. T. Choi, Alexander Mietke, and Jörn Dunkel.  Learning hydrodynamic equations for active matter from particle simulations and experiments.  Proceedings of the National Academy of Sciences, 120(7):e2206994120, February 2023.
  * [32] Shengze Cai, Zhicheng Wang, Sifan Wang, Paris Perdikaris, and George Em Karniadakis.  Physics-informed neural networks for heat transfer Problems.  Journal of Heat Transfer, 143(6):060801, June 2021.
  * [33] Ana Kaplarevic-Malisic, Branka Andrijevic, Filip Bojovic, Srdan Nikolic, Lazar Krstic, Boban Stojanovic, and Milos Ivanovic.  Identifying optimal architectures of physics-informed neural networks by evolutionary strategy.  Applied Soft Computing, page 110646, July 2023.
  * [34] Kevin Linka, Amelie Schäfer, Xuhui Meng, Zongren Zou, George Em Karniadakis, and Ellen Kuhl.  Bayesian physics informed neural networks for real-world nonlinear dynamical systems.  Computer Methods in Applied Mechanics and Engineering, 402:115346, December 2022.
  * [35] M. Raissi, P. Perdikaris, and G.E. Karniadakis.  Physics-informed neural networks: a deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations.  Journal of Computational Physics, 378:686–707, February 2019.
  * [36] Yeonjong Shin, Jerome Darbon, and George Em Karniadakis.  On the convergence of physics informed neural networks for linear second-order elliptic and parabolic type PDEs.  Communications in Computational Physics, 28(5):2042–2074, June 2020\.  arXiv: 2004.01806.
  * [37] Marc C. Kennedy and Anthony O’Hagan.  Bayesian calibration of computer models.  Journal of the Royal Statistical Society Series B: Statistical Methodology, 63(3):425–464, September 2001.
  * [38] David L. Banks and Mevin B. Hooten.  Statistical challenges in agent-based modeling.  The American Statistician, 75(3):235–242, July 2021.
  * [39] Sokratis Papadopoulos and Elie Azar.  Integrating building performance simulation in agent-based modeling using regression surrogate models: A novel human-in-the-loop energy modeling approach.  Energy and Buildings, 128:214–223, September 2016.
  * [40] Guus Ten Broeke, George Van Voorn, Arend Ligtenberg, and Jaap Molenaar.  The use of surrogate models to analyse agent-based models.  Journal of Artificial Societies and Social Simulation, 24(2):3, 2021\.
  * [41] Yi Zhang, Zhe Li, and Yongchao Zhang.  Validation and calibration of an agent-based model: a surrogate approach.  Discrete Dynamics in Nature and Society, 2020:1–9, January 2020\.
  * [42] Sushant S. Garud, Iftekhar A. Karimi, and Markus Kraft.  Design of computer experiments: a review.  Computers & Chemical Engineering, 106:71–95, November 2017.
  * [43] Ralph C. Smith.  Uncertainty quantification: theory, implementation, and applications.  Society for Industrial and Applied Mathematics, Philadelphia, PA, January 2013.
  * [44] Jari Kaipio and Ville Kolehmainen.  Approximate marginalization over modelling errors and uncertainties in inverse problems.  In Paul Damien, Petros Dellaportas, Nicholas G. Polson, and David A. Stephens, editors, Bayesian Theory and Applications, pages 644–672. Oxford University Press, January 2013.
  * [45] Jari P. Kaipio and Erkki Somersalo.  Statistical and computational inverse problems, volume 160 of Applied Mathematical Sciences.  Springer New York, New York, NY, 2005.
  * [46] Jari Kaipio and Erkki Somersalo.  Statistical inverse problems: discretization, model reduction and inverse crimes.  Journal of Computational and Applied Mathematics, 198(2):493–504, January 2007.
  * [47] Claudio Angione, Eric Silverman, and Elisabeth Yaneske.  Using machine learning as a surrogate model for agent-based simulations.  PLOS ONE, 17(2):e0263150, February 2022.
  * [48] Abdul Afram and Farrokh Janabi-Sharifi.  Black-box modeling of residential HVAC system and comparison of gray-box and black-box modeling methods.  Energy and Buildings, 94:121–149, May 2015.
  * [49] Daniel J. VandenHeuvel, Pascal R. Buenzli, and Matthew J. Simpson.  Pushing coarse-grained models beyond the continuum limit using equation learning, August 2023.  arXiv:2308.11086 [math, q-bio].
  * [50] Alexander Kurganov and Eitan Tadmor.  New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations.  Journal of Computational Physics, 160(1):241–282, May 2000.
  * [51] Linda Petzold.  Automatic selection of methods for solving stiff and nonstiff systems of ordinary differential equations.  SIAM Journal on Scientific and Statistical Computing, 4(1):136–148, March 1983.

##  Appendix A Coarse-graining ABMs into PDE models

We will coarse-grain the Pulling, Adhesion, and Pulling & Adhesion ABMs into
their mean-field PDE models. Each ABM consists of a combination of Rules A-F
from Figure 1. Each rule updates the occupancies of three consecutive lattice
sites, such as $\\{(i,j-1),(i,j),(i,j+1)\\}$. Recall from Section 2 that
$0_{i,j}(t)$, $P_{i,j}(t)$, and $H_{i,j}(t)$ denote the probabilities that the
individual lattice site $(i,j)$ is unoccupied, occupied by a pulling agent, or
occupied by an adhesive agent, respectively, at time $t$. To convert each rule
into a PDE model, we invoke the _mean-field assumption_ , which supposes that
all lattice site occupancies are independent of each other. This assumption
simplifies model coarse-graining by allowing us to replace the joint
probability of three lattice site occupancies with the product of the three
individual lattice site occupancy probabilities. For example, under the mean-
field assumption, we can write the probability that lattice sites
$(i,j-1),(i,j),\text{ and }(i,j+1)$ are all occupied by pulling agents at time
$t$ as $P_{i,j-1}(t)P_{i,j}(t)P_{i,j+1}(t)$; otherwise, we must consider the
joint occupancy probability for this triplet of lattice sites. Mean-field DE
models can poorly predict ABM behavior when the mean-field assumption is
violated during ABM simulations, see [9, 10, 13] for further details.

###  A.1 Coarse-graining the Pulling ABM

The Pulling ABM is composed of Rules A and B from Figure 1 and Section 2.2.1.
We begin coarse-graining this ABM into a PDE model by writing the master
equation governing how $P_{i,j}(t)$ changes according to these rules:

$\dfrac{\partial P_{i,j}(t)}{\partial
t}=K^{\ref{eq:ruleA}}+K^{\ref{eq:ruleB1}}+K^{\ref{eq:ruleB2}}.$ (11)

Rule A specifies how pulling agents migrate into an empty lattice site with
rate $r_{m}^{pull}/4$ when there is no neighboring agent in the lattice site
opposite the direction of migration. This rate is divided by four because the
agent randomly chooses to attempt to migrate into one of its four neighboring
lattice sites. We write this rule in the master equation as:

$\displaystyle K^{\ref{eq:ruleA}}=$
$\displaystyle-\dfrac{2r_{m}^{pull}}{4}\left[0_{i,j-1}(t)P_{i,j}(t)0_{i,j+1}(t)+0_{i-1,j}(t)P_{i,j}(t)0_{i+1,j}(t)\right]$
$\displaystyle+\dfrac{r_{m}^{pull}}{4}\left[0_{i,j-2}(t)P_{i,j-1}(t)0_{i,j}(t)+0_{i,j}P_{i,j+1}0_{i,j+2}+0_{i-2,j}P_{i-1,j}0_{i,j}+0_{i,j}P_{i+1,j}0_{i+2,j}\right],$
(12)

where the first line describes how a pulling agent moves out of lattice site
$(i,j)$, and the second line describes how a pulling agent moves into lattice
site $(i,j)$.

Rule B.1 specifies how a pulling agent migrates into an empty neighboring
lattice site and pulls its neighbor along with it, which occurs with
probability $p_{pull}$. We write this rule in the master equation as:

$\displaystyle K^{\ref{eq:ruleB1}}=-\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
P_{i,j}(t)P_{i,j+1}(t)0_{i,j+2}(t)+0_{i,j-2}(t)P_{i,j-1}(t)P_{i,j}(t)+$
$\displaystyle
P_{i,j}(t)P_{i+1,j}(t)0_{i+2,j}(t)+0_{i-2,j}(t)P_{i-1,j}(t)P_{i,j}(t)\bigg{]}$
$\displaystyle\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$ $\displaystyle
P_{i,j-2}(t)P_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)P_{i,j+1}(t)P_{i,j+2}(t)+$
$\displaystyle
P_{i-2,j}(t)P_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)P_{i+1,j}(t)P_{i+2,j}(t)\bigg{]}.$
(13)

Rule B.2 specifies how a pulling agent migrates into an empty neighboring
lattice site and fails to pull its neighbor along with it, which occurs with
probability $1-p_{pull}$. We write this rule in the master equation as:

$\displaystyle
K^{\ref{eq:ruleB2}}=-\dfrac{(1-p_{pull})r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
P_{i,j-1}(t)P_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)P_{i,j}(t)P_{i,j+1}(t)+$
$\displaystyle
P_{i-1,j}(t)P_{i,j}(t)0_{i+1,j+1}(t)+0_{i,j-1}(t)P_{i,j}(t)P_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{(1-p_{pull})r_{m}^{pull}}{4}\bigg{[}$ $\displaystyle
P_{i,j-2}(t)P_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)P_{i,j+1}(t)P_{i,j+2}(t)+$
$\displaystyle
P_{i-2,j}(t)P_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)P_{i+1,j}(t)P_{i+2,j}(t)\bigg{]}.$
(14)

To obtain the resulting PDE model for the Pulling ABM, we substitute Equations
(12), (13), and (14) into Equation (11) and set $0_{i,j}=1-P_{i,j}.$ We
replace each term with its Taylor expansion, up to second order:

$\displaystyle P_{i\pm m,j}(t)$ $\displaystyle=P_{i,j}(t)\pm
m\Delta(P_{i,j}(t))_{x}+\dfrac{m\Delta^{2}}{2}(P_{i,j}(t))_{xx}+\mathcal{O}(\Delta^{3}),$
$\displaystyle m=-2,-1,0,1,2;$ $\displaystyle P_{i,j\pm n}(t)$
$\displaystyle=P_{i,j}(t)\pm
n\Delta(P_{i,j}(t))_{y}+\dfrac{n\Delta^{2}}{2}(P_{i,j}(t))_{yy}+\mathcal{O}(\Delta^{3}),$
$\displaystyle n=-2,-1,0,1,2;$ (15)

where subscripts denote differentiation with respect the the shown variable,
and $\Delta$ is the length of each lattice site. As shown in the Mathematica
notebook Pulling_model_coarse_graining.nb, taking the limit of the resulting
expression as $\Delta\rightarrow 0$ leads to the mean-field PDE model for the
Pulling ABM:

$\dfrac{\partial P}{\partial
t}=\nabla\cdot\left(\dfrac{r_{m}^{pull}}{4}\left(1+3p_{pull}P^{2}\right)\nabla
P\right),$ (16)

where $P=P_{i,j}(t)$.

###  A.2 Coarse-graining the Adhesion ABM

The Adhesion ABM is composed of Rules C and D from Figure 1 and Section 2.2.2.
We begin coarse-graining this ABM into a PDE model by writing the master
equation governing how $H_{i,j}(t)$ changes according to these rules:

$\dfrac{\partial H_{i,j}(t)}{\partial
t}=K^{\ref{eq:ruleC}}+K^{\ref{eq:ruleD}}.$ (17)

Rule C specifies how adhesive agents migrate into an empty lattice site with
rate $r_{m}^{adh}/4$ when there is no neighboring agent in the lattice site
opposite the direction of migration. We write this rule in the master equation
as:

$\displaystyle K^{\ref{eq:ruleC}}=-\dfrac{2r_{m}^{adh}}{4}\bigg{[}$
$\displaystyle
0_{i,j-1}(t)H_{i,j}(t)0_{i,j+1}(t)+0_{i-1,j}(t)H_{i,j}(t)0_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{r_{m}^{adh}}{4}\bigg{[}$ $\displaystyle
0_{i,j-2}(t)H_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)H_{i,j+1}(t)0_{i,j+2}(t)+$
$\displaystyle
0_{i-2,j}(t)H_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)H_{i+1,j}(t)0_{i+2,j}(t)\bigg{]},$
(18)

where the first line describes how an adhesive agent moves out of lattice site
$(i,j)$, and the second and third lines describe how an adhesive agent moves
into lattice site $(i,j)$.

Rule D specifies how adhesive agents migrate into an empty neighboring lattice
site when a neighboring adhesive agent is in the lattice site opposite the
direction of migration. The neighboring adhesive agent attempts to adhere to
the migratory agent and abort the migration event. The adhesion event succeeds
with probability $p_{adh}$, and neither agent changes its position. The
adhesion event fails with probability $1-p_{adh}$, and the migratory agent
shifts into the previously-empty lattice site while the neighboring agent
remains in its previous lattice site. We write this rule in the master
equation as:

$\displaystyle K^{\ref{eq:ruleD}}=-\dfrac{(1-p_{adh})r_{m}^{adh}}{4}\bigg{[}$
$\displaystyle
H_{i,j-1}(t)H_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)H_{i,j}(t)H_{i,j+1}(t)+$
$\displaystyle
H_{i-1,j}(t)H_{i,j}(t)0_{i+1,j}(t)+0_{i-1,j}(t)H_{i,j}(t)H_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{(1-p_{adh})r_{m}^{adh}}{4}\bigg{[}$ $\displaystyle
H_{i,j-2}(t)H_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)H_{i,j+1}(t)H_{i,j+2}(t)+$
$\displaystyle
H_{i-2,j}(t)H_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)H_{i+1,j}(t)H_{i+2,j}(t)\bigg{]}.$
(19)

To obtain the resulting PDE model for the Adhesion ABM, we substitute
Equations (18) and (19) into Equation (17) and set $0_{i,j}=1-H_{i,j}$. We
replace each term with its Taylor expansion, up to second order:

$\displaystyle H_{i\pm m,j}(t)$ $\displaystyle=H_{i,j}(t)\pm
m\Delta(H_{i,j}(t))_{x}+\dfrac{m\Delta^{2}}{2}(H_{i,j}(t))_{xx}+\mathcal{O}(\Delta^{3}),$
$\displaystyle m=-2,-1,0,1,2;$ $\displaystyle H_{i,j\pm n}(t)$
$\displaystyle=H_{i,j}(t)\pm
n\Delta(H_{i,j}(t))_{y}+\dfrac{n\Delta^{2}}{2}(H_{i,j}(t))_{yy}+\mathcal{O}(\Delta^{3}),$
$\displaystyle n=-2,-1,0,1,2.$ (20)

As shown in the Mathematica notebook Adhesion_model_coarse_graining.nb, taking
the limit of the resulting expression as $\Delta\rightarrow 0$ leads to the
mean-field PDE model for the Adhesion ABM:

$\dfrac{\partial H}{\partial
t}=\nabla\cdot\left(\dfrac{r_{m}^{adh}}{4}\left(3p_{adh}\left(H-\dfrac{2}{3}\right)^{2}+1-\dfrac{4p_{adh}}{3}\right)\nabla
H\right)$ (21)

where $H=H_{i,j}(t)$.

###  A.3 Coarse-graining the Pulling & Adhesion ABM

The Pulling & Adhesion ABM is composed of Rules A to F from Figure 1 and
Sections 2.2.1-2.2.3. We begin coarse-graining this ABM into a PDE model by
writing the master system of equations governing how both $P_{i,j}(t)$ and
$H_{i,j}(t)$ change according to these rules:

$\displaystyle\dfrac{\partial P_{i,j}(t)}{\partial t}$
$\displaystyle=K^{\ref{eq:ruleA}}+K^{\ref{eq:ruleB1}}+K^{\ref{eq:ruleB2}}+K^{\ref{eq:rulee1}}_{P}+K^{\ref{eq:rulee2}}$
(22) $\displaystyle\dfrac{\partial H_{i,j}(t)}{\partial t}$
$\displaystyle=K^{\ref{eq:ruleC}}+K^{\ref{eq:ruleD}}+K^{\ref{eq:rulee1}}_{H}+K^{\ref{eq:rulef}},$
(23)

where $K^{\ref{eq:rulee1}}_{P}$ denotes how $P_{i,j}(t)$ is affected by Rule
E.1 and $K^{\ref{eq:rulee1}}_{H}$ denotes how $H_{i,j}(t)$ is affected by Rule
E.1. All other rules affect either $P_{i,j}(t)$ or $H_{i,j}(t)$, but not both.
Rules A-D are described in Sections A.1 and A.2, and we do not restate them
here.

Rule E specifies how a pulling agent migrates into an empty neighboring
lattice site when a neighboring adhesive agent is present in the lattice site
opposite the direction of migration. In Rule E.1, the pulling agent
successfully pulls the adhesive agent as it migrates, which occurs with
probability $p_{pull}$. In this scenario, the pulling agent shifts into the
previously-empty lattice site and the adhesive agent moves into the site
previously occupied by the pulling agent. We write this rule in the master
equation for $P_{i,j}(t)$ as:

$\displaystyle
K^{\ref{eq:rulee1}}_{P}=-\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
H_{i,j-1}(t)P_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)P_{i,j}(t)H_{i,j+1}(t)+$
$\displaystyle
H_{i-1,j}(t)P_{i,j}(t)0_{i+1,j}(t)+0_{i-1,j}(t)P_{i,j}(t)H_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$ $\displaystyle
H_{i,j-2}(t)P_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)P_{i,j+1}(t)H_{i,j+2}(t)+$
$\displaystyle
H_{i-2,j}(t)P_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)P_{i+1,j}(t)H_{i+2,j}(t)\bigg{]},$
(24)

and in the master equation for $H_{i,j}(t)$ as:

$\displaystyle
K^{\ref{eq:rulee1}}_{H}=-\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
0_{i,j-2}(t)P_{i,j-1}(t)H_{i,j}(t)+H_{i,j}(t)P_{i,j+1}(t)0_{i,j+2}(t)+$
$\displaystyle
0_{i-2,j}(t)P_{i-1,j}(t)H_{i,j}(t)+H_{i,j}(t)P_{i+1,j}(t)0_{i+2,j}(t)\bigg{]}$
$\displaystyle+\dfrac{p_{pull}r_{m}^{pull}}{4}\bigg{[}$ $\displaystyle
H_{i,j-1}(t)P_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)P_{i,j}(t)H_{i,j+1}(t)+$
$\displaystyle
H_{i-1,j}(t)P_{i,j}(t)0_{i+1,j}(t)+0_{i-1,j}(t)P_{i,j}(t)H_{i+1,j}(t)\bigg{]}.$
(25)

The neighboring adhesive agent successfully adheres to the migrating pulling
agent and aborts its migration event with probability $p_{adh}$. Neither
$P_{i,j}(t)$ or $H_{i,j}(t)$ changes in this scenario as no agents change
their locations in response to the adhesion event. In Rule E.2, the adhesive
agent fails to adhere to the pulling agent and the pulling agent fails to pull
the adhesive agent, which occurs with probability $1-p_{adh}-p_{pull}$. In
this scenario, the pulling agent shifts into the previously-empty lattice site
while the neighboring adhesive agent remains in its previous lattice site. We
write this rule in the master equation as:

$\displaystyle
K^{\ref{eq:rulee2}}=-\dfrac{(1-p_{adh}-p_{pull})r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
H_{i,j-1}(t)P_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)P_{i,j}(t)H_{i,j+1}(t)+$
$\displaystyle
H_{i-1,j}(t)P_{i,j}(t)0_{i+1,j}(t)++0_{i-1,j}(t)P_{i,j}(t)H_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{(1-p_{adh}-p_{pull})r_{m}^{pull}}{4}\bigg{[}$
$\displaystyle
H_{i,j-2}(t)P_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)P_{i,j+1}(t)H_{i,j+2}(t)+$
$\displaystyle
H_{i-2,j}(t)P_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)P_{i+1,j}(t)H_{i+2,j}(t)\bigg{]}.$
(26)

Rule F specifies how adhesive agents migrate into an empty neighboring lattice
site when a neighboring pulling agent is in the lattice site opposite the
direction of migration. The two agents do not interact with each other in this
scenario. As such, the adhesive agent migrates into the empty lattice site
with rate $r_{m}^{adh}/4$. We write this rule in the master equation as:

$\displaystyle K^{\ref{eq:rulef}}=-\dfrac{r_{m}^{adh}}{4}\bigg{[}$
$\displaystyle
P_{i,j-1}(t)H_{i,j}(t)0_{i,j+1}(t)+0_{i,j-1}(t)H_{i,j}(t)P_{i,j+1}(t)+$
$\displaystyle
P_{i-1,j}(t)H_{i,j}(t)0_{i+1,j}(t)+0_{i-1,j}(t)H_{i,j}(t)P_{i+1,j}(t)\bigg{]}$
$\displaystyle+\dfrac{r_{m}^{adh}}{4}\bigg{[}$ $\displaystyle
P_{i,j-2}(t)H_{i,j-1}(t)0_{i,j}(t)+0_{i,j}(t)H_{i,j+1}(t)P_{i,j+2}(t)+$
$\displaystyle
P_{i-2,j}(t)H_{i-1,j}(t)0_{i,j}(t)+0_{i,j}(t)H_{i+1,j}(t)P_{i+2,j}(t)\bigg{]}.$
(27)

To obtain the resulting system of differential equations for the Pulling &
Adhesion ABM, we substitute Equations (12), (13), (14), (18), (19), (24),
(25), (26), and (27) into Equation (23) and set $0_{i,j}=1-T_{i,j}$, where
$T_{i,j}=P_{i,j}+H_{i,j}$. We replace each term with its Taylor expansion up
to second order from Equations (15) and (20). As shown in the Mathematica
notebook Pulling-Adhesion_coarse_graining.nb, taking the limit of the
resulting expression as $\Delta\rightarrow 0$ leads to the mean-field system
of PDEs for the Pulling & Adhesion ABM:

$\displaystyle\dfrac{\partial P}{\partial t}=$
$\displaystyle\dfrac{r_{m}^{pull}}{4}\nabla\cdot\bigg{(}(1-T)\nabla P+P\nabla
T\bigg{)}$
$\displaystyle+p_{adh}\dfrac{r_{m}^{pull}}{4}\nabla\cdot\bigg{(}-3P(1-T)\nabla
H-H(1-T)\nabla P-HP\nabla T\bigg{)}$
$\displaystyle+p_{pull}\dfrac{r_{m}^{pull}}{4}\nabla\cdot\bigg{(}3P^{2}\nabla
T\bigg{)}$ $\displaystyle\dfrac{\partial H}{\partial t}=$
$\displaystyle\dfrac{r_{m}^{adh}}{4}\nabla\cdot\bigg{(}(1-T)\nabla H+H\nabla
T\bigg{)}$
$\displaystyle+p_{adh}\dfrac{r_{m}^{adh}}{4}\nabla\cdot\bigg{(}-4(1-T)H\nabla
H-H^{2}\nabla T\bigg{)}$
$\displaystyle+p_{pull}\dfrac{r_{m}^{pull}}{4}\nabla\cdot\bigg{(}-(1-T)H\nabla
P+(1-T)P\nabla H+3HP\nabla T\bigg{)},$ (28)

where $P=P_{i,j}(t),H=H_{i,j}(t),\text{ and }T=T_{i,j}(t)$.

##  Appendix B Numerical integration of PDEs

When simulating Equation (10), we populate the middle 20% of the spatial
dimension with 75% confluence and zero confluence everywhere else to match the
initial ABM configurations and implement no-flux boundary conditions:

$\displaystyle T(x,0)$ $\displaystyle=\begin{cases}0.75,&80\leq x\leq 120\\\
0,&\text{otherwise},\end{cases},$ $\displaystyle\dfrac{\partial T}{\partial
x}(0,t)$ $\displaystyle=\dfrac{\partial u}{\partial x}(X,t)=0.$ (29)

Before integration, we discretize the spatial domain as $x_{i}=i\Delta x$ with
$i=0,...,199$ and $\Delta x=1.0$. For ease of notation, let
$T_{i}(t)=T(x_{i},t)$ and $\mathcal{D}_{i}(t)=\mathcal{D}(T_{i}(t))$. We then
use the method of lines approach to integrate Equation (10). To discretize the
right hand side of Equation (10), we let

$\dfrac{\partial T_{i}(t)}{\partial x}\left(\mathcal{D}_{i}(t)\dfrac{\partial
T_{i}(t)}{\partial
x}\right)\approx\dfrac{P_{i+\nicefrac{{1}}{{2}}}(t)-P_{i-\nicefrac{{1}}{{2}}}(t)}{\Delta
x},$

where $P_{i\pm\nicefrac{{1}}{{2}}}(t)$ denotes the right or left flux through
location $x_{i}$, respectively. Following [50], we approximate these fluxes by

$\displaystyle P_{i+\nicefrac{{1}}{{2}}}(t)$
$\displaystyle=\dfrac{1}{2}\left(\mathcal{D}_{i}(t)\dfrac{T_{i+1}(t)-T_{i}(t)}{\Delta
x}+\mathcal{D}_{i+1}(t)\dfrac{T_{i+1}(t)-T_{i}(t)}{\Delta x}\right)$
$\displaystyle P_{i-\nicefrac{{1}}{{2}}}(t)$
$\displaystyle=\dfrac{1}{2}\left(\mathcal{D}_{i-1}(t)\dfrac{T_{i}(t)-T_{i-1}(t)}{\Delta
x}+\mathcal{D}_{i}(t)\dfrac{T_{i}(t)-T_{i-1}(t)}{\Delta x}\right).$ (30)

To implement the no-flux boundary conditions, we incorporate the ghost points
$x_{-1}$ and $x_{200}$ that enforce $u_{-1}(t)=u_{1}(t)$ and
$u_{198}(t)=u_{200}(t)$ into Equation (30). We integrate Equation (10) using
the odeint command from Scipy’s integration package (version 1.8.0), which
implements the Livermore Solver for Differential Equations (LSODA) method
[51].

##  Appendix C Supplementary figures

Figure 12: Forecasting Pulling ABM data with mean-field (MF) and BINN-guided PDE models. The mean-field and BINN-guided PDE simulations are used to forecast Pulling ABM data for (a-c) $r_{m}^{pull}=1.0,p_{pull}=0.8$ (d-f) $r_{m}^{pull}=0.9,p_{pull}=0.5$. Figure 13: Forecasting Adhesion ABM data with mean-field and BINN-guided PDE models. The mean-field and BINN-guided PDE simulations are used to forecast Adhesion ABM data for (a-c) $r_{m}^{adh}=1.0,p_{adh}=0.7$ (d-f) $r_{m}^{adh}=0.1,p_{adh}=0.5$. Figure 14: Forecasting Pulling & Adhesion ABM data with mean-field (MF) and BINN-guided PDE models. The mean-field and BINN-guided PDE simulations are used to forecast Pulling & Adhesion ABM data for the base parameter values ($r_{m}^{pull}=1.0,r_{m}^{adh}=0.25,p_{pull}=0.33,p_{adh}=0.33$, and $\alpha=0.5$), except (a-c) $p_{adh}=0.4$ (d-f) $r_{m}^{adh}=0.1.$ Figure 15: Predicting Adhesion ABM data with the interpolated PDE model. The interpolated PDE model predicts Adhesion ABM data for (a-c) $r_{m}^{adh}=1.0$ and $p_{adh}=0.95$. Sample | $\bm{p}=(r_{m}^{adh},\ {p_{adh}})^{T}$
---|---
1 | (0.145, 0.825)T
2 | (0.505, 0.575)T
3 | (0.415, 0.725)T
4 | (0.865, 0.525)T
5 | (0.955, 0.625)T
6 | (0.235, 0.775)T
7 | (0.685, 0.675)T
8 | (0.325, 0.875)T
9 | (0.775, 0.925)T
10 | (0.595, 0.975)T
Table 4: Latin hypercube sampling for the Adhesion ABM. The samples from the new parameter dataset for the Adhesion ABM when varying $r_{m}^{adh}$ and $p_{adh}$. The samples are ordered by increasing testing MSE values (see Figure 10(c)). Figure 16: Predicting Adhesion ABM data with the interpolated PDE model. The interpolated PDE model predicts Adhesion ABM data for (a-c) $r_{m}^{adh}=0.595$ and $p_{adh}=0.975$ and (d-f) $r_{m}^{adh}=0.325$ and $p_{adh}=0.875$. Sample | $\bm{p}=\ (r_{m}^{pull},\ r_{m}^{adh},\ p_{pull},\ p_{adh},\ \alpha)^{T}$
---|---
1 |  (1.0, 0.25, 0.394, 0.578, 0.912)T
2 |  (1.0, 0.25, 0.293, 0.528, 0.938)T
3 |  (1.0, 0.25, 0.008, 0.226, 0.988)T
4 |  (1.0, 0.25, 0.511, 0.477, 0.862)T
5 |  (1.0, 0.25, 0.41, 0.109, 0.962)T
6 |  (1.0, 0.25, 0.075, 0.595, 0.888)T
7 |  (1.0, 0.25, 0.042, 0.544, 0.838)T
8 |  (1.0, 0.25, 0.327, 0.059, 0.712)T
9 |  (1.0, 0.25, 0.444, 0.31, 0.662)T
10 |  (1.0, 0.25, 0.209, 0.209, 0.612)T
11 |  (1.0, 0.25, 0.126, 0.41, 0.762)T
12 |  (1.0, 0.25, 0.193, 0.042, 0.588)T
13 |  (1.0, 0.25, 0.059, 0.561, 0.462)T
14 |  (1.0, 0.25, 0.243, 0.26, 0.788)T
15 |  (1.0, 0.25, 0.427, 0.494, 0.512)T
16 |  (1.0, 0.25, 0.595, 0.327, 0.812)T
17 |  (1.0, 0.25, 0.025, 0.461, 0.388)T
18 |  (1.0, 0.25, 0.377, 0.176, 0.488)T
19 |  (1.0, 0.25, 0.226, 0.645, 0.538)T
20 |  (1.0, 0.25, 0.528, 0.126, 0.688)T
21 |  (1.0, 0.25, 0.561, 0.075, 0.562)T
22 |  (1.0, 0.25, 0.142, 0.193, 0.362)T
23 |  (1.0, 0.25, 0.31, 0.092, 0.738)T
24 |  (1.0, 0.25, 0.176, 0.662, 0.412)T
25 |  (1.0, 0.25, 0.645, 0.008, 0.638)T
26 |  (1.0, 0.25, 0.343, 0.293, 0.312)T
27 |  (1.0, 0.25, 0.092, 0.611, 0.238)T
28 |  (1.0, 0.25, 0.109, 0.628, 0.012)T
29 |  (1.0, 0.25, 0.159, 0.343, 0.212)T
30 |  (1.0, 0.25, 0.26, 0.142, 0.188)T
31 |  (1.0, 0.25, 0.36, 0.377, 0.262)T
32 |  (1.0, 0.25, 0.276, 0.36, 0.038)T
33 |  (1.0, 0.25, 0.578, 0.243, 0.288)T
34 |  (1.0, 0.25, 0.628, 0.159, 0.062)T
35 |  (1.0, 0.25, 0.477, 0.511, 0.138)T
36 |  (1.0, 0.25, 0.611, 0.276, 0.338)T
37 |  (1.0, 0.25, 0.461, 0.444, 0.162)T
38 |  (1.0, 0.25, 0.544, 0.427, 0.112)T
39 |  (1.0, 0.25, 0.494, 0.394, 0.088)T
40 |  (1.0, 0.25, 0.662, 0.025, 0.438)T
Table 5: Latin hypercube sampling for the Pulling & Adhesion ABM. The samples from the prior parameter dataset for the Pulling & Adhesion ABM when varying $p_{pull}$, $p_{adh}$, and $\alpha$. The samples are ordered by increasing testing MSE values. Sample | $\bm{p}=\ (r_{m}^{pull},\ r_{m}^{adh},\ p_{pull},\ p_{adh},\ \alpha)^{T}$
---|---
1 | (1.0, 0.25, 0.285, 0.519, 0.775)T
2 | (1.0, 0.25, 0.419, 0.352, 0.875)T
3 | (1.0, 0.25, 0.486, 0.117, 0.525)T
4 | (1.0, 0.25, 0.553, 0.285, 0.375)T
5 | (1.0, 0.25, 0.385, 0.586, 0.475)T
6 | (1.0, 0.25, 0.586, 0.184, 0.175)T
7 | (1.0, 0.25, 0.62, 0.151, 0.325)T
8 | (1.0, 0.25, 0.184, 0.084, 0.625)T
9 | (1.0, 0.25, 0.352, 0.385, 0.925)T
10 | (1.0, 0.25, 0.653, 0.05, 0.275)T
11 | (1.0, 0.25, 0.151, 0.653, 0.075)T
12 | (1.0, 0.25, 0.452, 0.251, 0.125)T
13 | (1.0, 0.25, 0.084, 0.218, 0.225)T
14 | (1.0, 0.25, 0.318, 0.62, 0.725)T
15 | (1.0, 0.25, 0.519, 0.017, 0.825)T
16 | (1.0, 0.25, 0.117, 0.419, 0.425)T
17 | (1.0, 0.25, 0.251, 0.486, 0.975)T
18 | (1.0, 0.25, 0.017, 0.452, 0.025)T
19 | (1.0, 0.25, 0.05, 0.318, 0.575)T
20 | (1.0, 0.25, 0.218, 0.553, 0.675)T
Table 6: Latin hypercube sampling for the Pulling & Adhesion ABM. The samples
from the new parameter dataset for the Pulling & Adhesion ABM when varying
$p_{pull}$, $p_{adh}$, and $\alpha$. The samples are ordered by increasing
testing MSE values (see Figure 11(b)). Figure 17: Predicting Pulling &
Adhesion ABM data with the interpolated PDE model. The interpolated PDE model
predicts Adhesion ABM data for $r_{m}^{pull}=1.0$, $r_{m}^{adh}=0.25$, and
(a-c) $p_{pull}=0.218$, $p_{adh}=0.553$, and $\alpha=0.675$ (d-f)
$p_{pull}=0.251$, $p_{adh}=0.486$, and $\alpha=0.975$. Figure 18:
Computational expenses of each modeling approach. Violin plots represent the
distribution of wall time computations for ABM simulations, BINN training,
mean-field PDE simulations, and BINN-guided PDE simulations for the (a)
Pulling ABM, (b) Adhesion ABM, and (c) Pulling & Adhesion ABM.

##  Appendix D Gillespie algorithm

Create an $X\times Y$ lattice with user-specified placement of agents

Set $t=0$

Set maximum simulation time $t_{\text{end}}$

Set $P(t)$ and $H(t)$ equal to the number of Pulling and Adhesive agents on
the lattice, respectively

while  _$t <t_{\text{end}}$_ do

Calculate the following random variables, uniformly distributed on
$[0,1]:\gamma_{1},\gamma_{2}$

Calculate the propensity function $a(t)=r_{m}^{pull}P(t)+r_{m}^{adh}H(t)$

Calculate time step $\tau=-\ln(\gamma_{1})/a(t)$

$t=t+\tau$

$R=a(t)\gamma_{2}$

if  _$R <r_{m}^{pull}P(t)$_ then

Perform Pulling agent migration (Algorthm S2)

else if  _$R <r_{m}^{pull}P(t)+r_{m}^{adh}H(t)$_ then

Perform Adhesive agent migration (Algorthm S3)

end while

Algorithm 1 Gillespie algorithm for the Pulling & Adhesion ABM

Randomly choose a pulling agent and determine its lattice site index,
$\vec{x}=(i,j)^{T}$

Choose one of the four cardinal migration directions,
$\vec{dx}=(dx,dy)^{T}\in\\{(1,0)^{T},(-1,0)^{T},(0,1)^{T},(0,-1)^{T}\\}$, with
equal probability, $1/4$. The neighboring direction is given by
$\hat{dx}=-\vec{dx}$

if  _$\vec{x}+\vec{dx}$ is empty_ then

if  _$\vec{x}+\hat{dx}$ is empty_ then

/* Rule A */

Move the chosen pulling agent to lattice site $\vec{x}+\vec{dx}$

else if  _$\vec{x}+\hat{dx}$ is occupied by a Pulling agent_ then

/* Rule B */

Calculate the random variable, $\gamma_{3}$, uniformly distributed on $[0,1]$

if  _$\gamma_{3}\leq p_{pull}$_ then

Move the chosen pulling agent to lattice site $\vec{x}+\vec{dx}$

Move the neighboring agent to lattice site $\vec{x}$

else if  _$\gamma_{3} >p_{pull}$_ then

Move the chosen pulling agent to lattice site $\vec{x}+\vec{dx}$

else if  _$\vec{x}+\hat{dx}$ is occupied by an Adhesive agent_ then

/* Rule E */

Calculate the random variable, $\gamma_{3}$, uniformly distributed on $[0,1]$

if  _$\gamma_{3}\leq p_{pull}$_ then

Move the chosen pulling agent to lattice site $\vec{x}+\vec{dx}$

Move the neighboring agent to lattice site $\vec{x}$

else if  _$\gamma_{3}\leq p_{pull}+1-p_{adh}$ _ then

Move the chosen pulling agent to lattice site $\vec{x}+\vec{dx}$

Algorithm 2 Pulling Agent migration

Randomly choose an adhesive agent and determine its lattice site index,
$\vec{x}=(i,j)^{T}$

Choose one of the four cardinal migration directions,
$\vec{dx}=(dx,dy)^{T}\in\\{(1,0)^{T},(-1,0)^{T},(0,1)^{T},(0,-1)^{T}\\}$, with
equal probability, $1/4$. The neighboring direction is given by
$\hat{dx}=-\vec{dx}$

if  _$\vec{x}+\vec{dx}$ is empty_ then

if  _$\vec{x}+\hat{dx}$ is empty_ then

/* Rule C */

Move the chosen adhesive agent to lattice site $\vec{x}+\vec{dx}$

else if  _$\vec{x}+\hat{dx}$ is occupied by an adhesive agent_ then

/* Rule D */

Calculate the random variable, $\gamma_{3}$, uniformly distributed on $[0,1]$

if  _$\gamma_{3}\leq(1-p_{adh})$ _ then

Move the chosen adhesive agent to lattice site $\vec{x}+\vec{dx}$

else if  _$\vec{x}+\hat{dx}$ is occupied by a Pulling agent_ then

/* Rule F */

Move the chosen adhesive agent to lattice site $\vec{x}+\vec{dx}$

Algorithm 3 Adhesive agent migration