remarkRemark hypothesisHypothesis claimClaim Learning Heterogeneous Agent
PopulationsD. Messenger, G. Wheeler, X. Liu, D. Bortz.

# Learning Anisotropic Interaction Rules from Individual Trajectories in a
Heterogeneous Cellular Population

Daniel A. Messenger Department of Applied Mathematics, University of Colorado,
Boulder, CO 80309-0526, USA<EMAIL_ADDRESS><EMAIL_ADDRESS>Graycen E. Wheeler Department of Biochemistry,
University of Colorado, Boulder, CO 80309-0526, USA. (, ).
<EMAIL_ADDRESS><EMAIL_ADDRESS>Xuedong
Liu22footnotemark: 2 David M. Bortz11footnotemark: 1

###### Abstract

Interacting particle system (IPS) models have proven to be highly successful
for describing the spatial movement of organisms. However, it has proven
challenging to infer the interaction rules directly from data. In the field of
equation discovery, the Weak form Sparse Identification of Nonlinear Dynamics
(WSINDy) methodology has been shown to be very computationally efficient for
identifying the governing equations of complex systems, even in the presence
of substantial noise. Motivated by the success of IPS models to describe the
spatial movement of organisms, we develop WSINDy for second order IPSs to
model the movement of communities of cells. Specifically, our approach learns
the directional interaction rules that govern the dynamics of a heterogeneous
population of migrating cells. Rather than aggregating cellular trajectory
data into a single best-fit model, we learn the models for each individual
cell. These models can then be efficiently classified according to the active
classes of interactions present in the model. From these classifications,
aggregated models are constructed hierarchically to simultaneously identify
different species of cells present in the population and determine best-fit
models for each species. We demonstrate the efficiency and proficiency of the
method on several test scenarios, motivated by common cell migration
experiments.

###### keywords:

Weak form, interacting particle systems, heterogeneous agents, system
identification, unsupervised learning.

##  1 Introduction

Systems of autonomous agents are ubiquitous in the natural world. Research
into their behavior has led to a plethora of proposed mathematical models,
including the agent-based “boids” model [48], ordinary differential equation
models for milling and flocking [32, 16], and nonlocal partial differential
equations [40, 60], to name a few. A general framework for rigorous analysis
of these models is by now very mature [12].

Identifying the rules of interaction between agents is necessary for
predicting and influencing the cooperative abilities of any such system,
whether composed of autonomous robots, large multi-cellular animals, single-
celled organisms, or even molecules. Methods for inferring the rules of
interaction between agents using observed trajectory data have continued to
advance since the early 2000’s. Several of the principled techniques include
force-matching [19, 27], linear regression [35, 34], mean-field formulations
[57, 36], information-theoretic tools [46], underdamped Langevin regression
[8, 9], Gaussian processes [23], and even a method based on topological rather
than metric distances [2].

These and related techniques have been successfully used to identify the
dominant drivers of collective behavior in a variety of social and biological
systems [52, 56, 1], including schools of fish [67, 62], flocks of birds [15,
45], and pedestrian traffic [68], all directly incorporating measured
trajectory data. While popular methods, such as force-matching, are useful in
identifying fields of vision and spatial statistics of interactions, they
cannot easily disentangle the combined effects of multiple forces (e.g.
attraction, repulsion, and alignment) [41, 20], let alone different
interactions between multiple species of neighbors. This limits the classes of
models they can identify and implies that new methods must be developed for
heterogeneous populations.

The field of equation discovery is a highly active area of research [30, 42,
53, 49, 10, 51, 63, 33, 14] as it offers tools to directly learn governing
differential equations. This approach is not only useful in prediction and
validation, but can be used to simultaneously identify multiple active modes
of inter-agent communication, such as repulsion, velocity alignment, and
attraction. In this work, we tackle the problem of identifying governing
equations for an interacting particle system (IPS) with multiple interacting
species. Our proposed approach is completely naive with regard to species
membership in order to specifically address problems of heterogeneity in
collective cell migration [59].

Motivated by existing hypotheses regarding the anistropy of cell-cell
interactions [50, 26, 44], we introduce our framework in the context of
directional interaction models, as defined below. Moreover, we note that the
documented significance of anisotropic interactions in general collective
systems [17, 3, 70] suggests that our approach may have wide applicability.

###  1.1 Heterogeneous Populations

Many collective populations arising in nature are inherently heterogeneous,
with the rules of interaction varying across different subsets of the
population. This is readily observable in complex mammalian populations, but
is also seen in simpler organisms, such as honey bee swarms, where bees divide
into scout and worker bee roles [55]. The advantages of heterogeneity in
collective behavior have even inspired search optimization algorithms [18,
28].

At the level of microorganisms, cells have been observed to adopt leader-like
and follower-like roles during migration [66], without the aid of a central
nervous system. Individual cell speed and persistence of motion have also been
determined to be functions of the age and size of the cell [31, 6, 5], which
may lead to heterogeneous responses to stimuli from neighboring cells. The
mechanisms which produce these heterogeneities, and the extent to which
heterogeneity is present in a given cell population, are current subjects of
debate [43, 54, 25]. Data-driven techniques may be useful in formulating
accurate mathematical models in the presence of heterogeneity.

The authors of [69] develop a highly versatile method for inferring explicit
rules of interaction in a heterogeneous population, although it is assumed
that species membership is known a priori. Several recent works have offered
methods of assessing the degree of population heterogeneity [54, 52], yet
these methods do not provide explicit mathematical models for the different
populations. In contrast, the method presented here allows one to classify the
given population into different species according to the heterogeneous
interaction rules present, and produces explicit mathematical models for each
species as a byproduct.

###  1.2 Directional Interaction Forces

It is now well-known that simple radial interaction models are incapable of
explaining many observed collective behaviors in biological settings, and that
directionally-dependent interaction rules, based on a limited field of view or
sensing angle, offer a significant advantage [17, 15, 11, 13, 21]. At the
cellular level, directional dependence of cell-cell interaction has been
proposed in the context of intracellular polarization [26], however the
cellular sensing range is not immediately obvious, since a migrating cell does
not have an obvious “field of view”. Recent works have sought to quantify the
degree to which interactions are density-dependent [7], but not which
directional modes (radial, dipolar, quadrupolar, etc.) are dominant during a
collective migration event.

In addition to providing an explanation for certain observed phenomena [70],
directional interaction rules are capable of generating spontaneous migration,
due to the total directional force between particles not being conserved in
general. In the modeling of active matter systems (such as migrating cells)
[47, 57], such symmetry breaking is commonly generated by a combination of
Brownian forcing and a self-propulsion device [65]. However, it is not clear
that self-propulsion is an appropriate mechanism for modeling cellular
movement (in comparison to fish, which are constantly swimming). Directional
forces may then be an important mechanism for symmetry breaking and
spontaneous cellular migration.

###  1.3 Weak-Form Sparse Identification of Nonlinear Dynamics

At its core, our method involves learning ordinary differential equations for
cells using available trajectory data. For this we employ the weak-form sparse
identification of nonlinear dynamics algorithm (WSINDy), which has been shown
to successfully identify governing equations from data at the levels of
ordinary different equations [38], partial differential equations [37], 1st-
order interacting particle systems [36], and even works in a small-storage
streaming scenario [39].

A significant advantage of the WSINDy method is that it identifies a single
governing equation which can be analyzed and simulated using conventional
techniques of applied mathematics. It does not involve any black-box
algorithms or mappings as would be generated in using a neural network-based
approach.111See [4] for an example where the authors first learn a neural
network model of the potential and then uses sparse identification to learn
the algebraic form of the potential.

Several alternative methods have been developed to accomplish the equation
learning task for particle systems. In particular, Lu et al. [34] develop a
method for learning general feature-dependent 2nd-order interaction rules for
heterogeneous populations, where features may include directional interaction
forces, speed dependence, and so on. The differences between this and our work
are the following. (i) We are performing the unsupervised learning task of
classifying agents by their interaction rules, whereas Lu et al.’s work
assumes knowledge of the species membership, (ii) we are interested in sparse
model representations, in particular selection of the correct modes of
interaction (e.g. attractive, repulsive, alignment, and drag force), whereas
Lu’s work assumes knowledge of both the feature-dependence and types of forces
present (e.g. for planetary systems, priori knowledge is used to rule out the
presence of an alignment force). (iii) Lastly, models are initially extracted
from single-cell trajectories. As described in the next section, rather than
aggregating data which may come from multiple cell species, we aggregate
models which are likely to describe the same species, and then use the
aggregate model to perform classification.

###  1.4 Single-Cell Learning and Model Clustering

With a possibly heterogeneous population of cell trajectories available, one
is tasked with the problem of deciding how to aggregate the data. If knowledge
of the underlying species membership is available, a more accurate model can
be inferred by pooling data from all individuals of a given species. On the
other hand, pooling data from multiple species into a single model can result
in a highly inaccurate model depending on the difference between interaction
rules in each species. In general, there exists a spectrum of possible pooling
strategies, ranging from learning few models from large subsets of the
population, to learning many models from small subsets of the populations. The
former intrinsically produces models with high bias and low variance, while
the latter produces models with low bias and high variance. Such pooling
strategies have been recently explored in [22], where it is found that
identifying a single model can be improved by pooling models learned from
subsets of the data. However, this has not been extended to classifying the
data itself into species, and finding a model for each species. Moreover, the
IPS setting offers a particular advantage on the subject of model validation,
as data can easily be assimilated into forward simulations.

In this work we investigate the extreme case of learning an individual model
${\mathcal{M}}_{i}$ for the $i$-th individual trajectory, and then clustering
the set of learned models
${\mathcal{M}}:=\\{{\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N}\\}$ according to
their identified modes of interaction. This approach is counter-intuitive
because there is no guarantee that a single cell trajectory will provide
enough information on the interaction rules of its species. To be able to
classify cells using the (potentially) insufficiently informative
trajectories, we developed an ad hoc recursive classifier which we show (in
Section 4) accurately clusters and sorts the models into species. This
approach prevents any contamination that may result from combining
trajectories of multiple species.

Once the models are clustered, an aggregate model $\overline{{\mathcal{M}}}$
is computed by averaging the models in ${\mathcal{M}}$ belonging to the most
populous cluster. The model $\overline{{\mathcal{M}}}$ is then used to
classify cells via forward simulations which are made highly efficient by
directly incorporating the data. In particular, for each trajectory in the
dataset, we use $\overline{{\mathcal{M}}}$ to simulate a new trajectory, but
with all neighbor interactions computed using the data. That is, only the new
trajectory is propagated forward in time by model $\overline{{\mathcal{M}}}$,
while the rest of the population is simply the data itself. This can then be
trivially parallelized, reducing an ${\mathcal{O}}(N^{2})$ computational cost
per timestep to $N$ cores performing ${\mathcal{O}}(N)$ updates per timestep
with no communication overhead.

We show through examples below that this hierarchical model-pooling and
validation procedure produces both correct species classification and accurate
governing equations, despite individual cell trajectory data carrying low
levels of information. For further information on the classification
algorithm, see Section 3.

###  1.5 Paper outline

In Section 2 we discuss the general form of directional interacting particle
models that will be assumed in the learning process. In Section 3 we introduce
our model selection and classification algorithm, which is composed of the
five steps: (a) learn single-cell models, (b) cluster learned models according
to force modes, (c) form an aggregate model by averaging models in the largest
cluster, (d) validate the aggregate model using data-driven forward
simulations, (e) classify cells according to performance under the aggregate
model. In Section 4 we examine the performance of the algorithm in learning
and classifying homogeneous and heterogeneous populations of one, two, and
three species. We discuss possible next directions in Section 5. Some
additional information and a summary of notation are included in the appendix.

##  2 Directional Interacting Particle Models

We use a general 2nd-order directional interaction model framework, where the
position and velocity $(x_{i},v_{i})$ of cell $i$ are governed by the
differential equations

(2.1)
$\begin{dcases}\ddot{x}_{i}=\frac{1}{N_{\text{tot}}}\sum_{j=1}^{N_{\text{tot}}}f_{\text{a-r}}(|x_{i}-x_{j}|,\theta_{ij})(x_{i}-x_{j})\\\
\qquad+\frac{1}{N_{\text{tot}}}\sum_{j=1}^{N_{\text{tot}}}f_{\text{align}}(|x_{i}-x_{j}|,\theta_{ij})(v_{i}-v_{j})\\\
\qquad+\frac{1}{N_{\text{tot}}}\sum_{j=1}^{N_{\text{tot}}}f_{\text{drag}}(|v_{i}|,\theta_{ij})v_{i}+\eta_{i}.\end{dcases}$

Here, $\theta_{ij}$ is the angle between $v_{i}$ and $x_{j}-x_{i}$ (see the
diagram in Figure 2.1) and $(\eta_{i})_{i\in[N_{tot}]}$ are white noise
processes. The attractive-repulsive force $f_{\text{a-r}}$, alignment force
$f_{\text{align}}$, and the drag force $f_{\text{drag}}$ define the rules by
which cell $i$ communicates with the rest of the population. Our primary
objective is to identify a set of interaction rules
$\\{(f_{\text{a-r}},f_{\text{align}},f_{\text{drag}})_{\ell}\\}_{1\leq\ell\leq
S}$, one for each of the $S$ species present in the population.

$\displaystyle\theta_{i}{}_{j}$$\displaystyle x_{i}$$\displaystyle
x_{j}$$\displaystyle v_{i}$$\displaystyle v_{j}$$\displaystyle x_{j}-x_{i}$
Figure 2.1: Diagram of social interactions depending on angle $\theta_{ij}$
between cell $i$’s velocity and cell $j$’s position relative to $i$.

###  2.1 Directionality $\theta_{ij}$

As mentioned above, directional variation in the interaction forces between
cells can arise from various factors, including intracellular polarization, or
heterogeneous distribution of membrane-bound receptors, asymmetry in the
protrusion/contraction of lamellopodia as the cell crawls on the substrate,
and so on. In the current study, we assume that each of these effects is
unobservable, hence we model the aggregate directional effect using the angles
$\theta_{ij}$, depicted in Figure 2.1. Dependence on angle $\theta_{ij}$
allows for interactions between cell $i$ and cell $j$ to vary depending on the
direction of motion222Note that asymmetric interactions
$\theta_{ij}\neq\theta_{ji}$ lead to symmetry breaking and spontaneous cell
migration from an initially motionless state.. Put another way, in the
reference frame of cell $i$, the polar coordinates $r_{ij}=|x_{i}-x_{j}|$ and
$\theta_{ij}$ allow one to represent any interaction force that varies over
the 2D plane.

In this study we restrict the angular dependence to
$\\{1,\cos(\theta_{ij}),\cos(2\theta_{ij})\\}$, which allows for a combination
of radial, dipolar, and quadrupolar interactions (see Figure 4.1 for examples
of dipolar (right) and quadrupolar (left) forces used in this study). Higher-
order directionality can usually be assumed to be negligible, however
extension to higher modes is straightforward.

###  2.2 Attractive-repulsive force $f_{\text{a-r}}$

The interaction force $f_{\text{a-r}}$ acts along the vector from cell $i$ to
cell $j$ and captures short-range repulsion and long-range signaling. Many IPS
models include only an attractive-repulsive force, due to its extensive
pattern-forming capabilities [61, 24].

We impose the following natural constraints on $f_{\text{a-r}}$:

(2.2) $f_{\text{a-r}}(r,\theta)\begin{dcases}\geq 0,&0\leq r<r_{nf}\\\ \leq
0,&r\geq r_{ff}\\\
\in\text{span}\\{1,\cos(\theta),\cos(2\theta)\\},&\text{every $r$
fixed,}\end{dcases}$

where $r_{nf}$ is the near-field threshold and $r_{ff}$ is the far-field
threshold, i.e., a large distance.

The first inequality enforces that $f_{\text{a-r}}$ is near-field repulsive,
which must be true by volume exclusion. In practice we define the $r_{nf}$
implicitly by

$\mathbb{P}\left(|x_{i}-x_{j}|<r_{nf}\right)=p_{nf},$

where in this work we set $p_{nf}=0.001$, and the dataset is used to compute
the probability. This states that the near-field region is defined by short-
distance interactions which account for less than $0.1\%$ of observed
neighbor-neighbor distances.

The second equality enforces long-range decay, as well as model stability.
Decay is natural since interactions can be expected to be small outside of
some large distance $r_{ff}$. We enforce that interactions are attractive at
large distances (allowing for decay as well), so that in simulation the
particles do not spread to infinity. We set $r_{ff}=1$ throughout, although
$r_{ff}$ can easily be chosen from the data (e.g. $r_{ff}=50r_{nf}$
corresponds to an effective interaction range of 50 cell radii). (See Appendix
A for resulting values of $r_{nf}$ and $r_{ff}$ and other hyperparameters for
examples below).

The third set inclusion simply reiterates the assumptions on directionality
described above.

###  2.3 Alignment force $f_{\text{align}}$

The alignment force $f_{\text{align}}$ captures cells’ tendency to match the
velocity of neighboring cells. There are many theories as to how this arises
physically [58, 50]. Perhaps protrusions from cells inform the cell about the
bulk direction of motion, which would be a very local effect. However,
alignment models which have been proven to lead to flocking depend on
sufficiently long-range alignment [16].

We impose the following constraints on $f_{\text{align}}$:

(2.3) $f_{\text{align}}(r,\theta)\begin{dcases}\leq 0,&0\leq r\\\
\in\text{span}\\{1,\cos(\theta),\cos(2\theta)\\},&\text{every $r$
fixed}\end{dcases}$

The first inequality enforces that $f_{\text{align}}$ is non-positive, which
is necessary for the constant velocity state $v_{i}=v_{j}=v$ to be a stable
configuration. If not, small perturbations away from $v_{i}=v_{j}$ result in
cells accelerating away from each other, which is a redundant force given that
cells can be pushed away from each other through $f_{\text{a-r}}$ (it is also
not hard to see that $f_{\text{align}}>0$ is unphysical). The second
constraint restricts the alignment force to be a combination of radial,
dipolar, or quadrupolar modes, similar to $f_{\text{a-r}}$.

###  2.4 Drag force $f_{\text{drag}}$

The drag force $f_{\text{drag}}$ captures energy expenditure due to general
resistance to motion (resulting e.g. from substrate roughness), however we
allow an angular dependence on $\theta_{ij}$ to capture possible decreases or
increases in drag depending on local neighbor distribution. For this we impose
the following constraints:

(2.4) $f_{\text{drag}}(s,\theta)\begin{dcases}\leq 0,&0\leq s<\infty\\\
\in\text{span}\\{1,\cos(\theta)\\},&\text{every $s$ fixed}\end{dcases}$

where $s$ indicates the speed of the cell. The force $f_{\text{drag}}$ is
chosen to be negative so that cells do not have a “self-propulsion” device. As
mentioned previously, many models of active matter include self-propulsion as
a partial mechanism for symmetry breaking and general non-equilibrium effects.
To reiterate, we do not expect cells to have a self-propulsion device, in
fact, we wish to learn how migration occurs spontaneously, incited by pairwise
interactions. In addition, a positive drag force leads to populations
spreading outside of the range of meaningful interactions. In this way,
negative drag is computationally beneficial, as it leads to improved model
stability.

##  3 Algorithm

Figure 3.1: Classification pipeline for cells from heterogeneous populations.
(a) An ensemble of models
${\mathcal{M}}=\\{{\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N}\\}$ is learned,
each from an individual trajectory; (b) ${\mathcal{M}}$ is partitioned into
clusters ${\mathcal{C}}=\\{{\mathcal{C}}_{1},\dots,{\mathcal{C}}_{r}\\}$
according to active forces in each model; (c) models in the largest cluster
$\overline{{\mathcal{C}}}$ are averaged together, giving
$\overline{{\mathcal{M}}}$; (d) $\overline{{\mathcal{M}}}$ is validated along
each individual trajectory; (e) validation errors are classified, producing an
identified species $\overline{{\mathcal{S}}}$ (cyan checkmarks) and the
remaining cells (red X’s) are returned to step (c) to be clustered again.
Steps (b)-(e) repeat until all cells are classified. Note that the number and
members of model clusters ${\mathcal{C}}$ and resulting aggregate model
$\overline{{\mathcal{M}}}$ will change each iteration depending on the
identity of remaining unlabeled cells.

Our algorithm involves first learning an ensemble of directional force models
${\mathcal{M}}=\\{{\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N}\\}$, that is, one
model for each of the $N$ focal cells selected for learning. We next group
models into clusters
${\mathcal{C}}:=\\{{\mathcal{C}}_{1},\dots,{\mathcal{C}}_{r}\\}$ and compute
an aggregrate model $\overline{{\mathcal{M}}}$ from the largest cluster,
denoted by $\overline{{\mathcal{C}}}$. A new species
$\overline{{\mathcal{S}}}$ is then identified, with membership in
$\overline{{\mathcal{S}}}$ determined by the accuracy of data-driven forward
simulations of model $\overline{{\mathcal{M}}}$. Cells in the species
$\overline{{\mathcal{S}}}$ are then removed from the population and the
remaining cells are returned to the clustering phase. More explicitly, the
algorithm is composed of the following steps, which are visualized in Figure
3.1.

  1. (a)

Identify individual cell models
${\mathcal{M}}=\\{{\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N}\\}$ using the
WSINDy algorithm

  2. (b)

Partition learned models ${\mathcal{M}}$ into clusters
$\\{{\mathcal{C}}_{1},\dots,{\mathcal{C}}_{r}\\}$ according to active force
modes

  3. (c)

Identify the largest cluster $\overline{{\mathcal{C}}}$ and average together
the models in $\overline{{\mathcal{C}}}$ to arrive at a single model
$\overline{{\mathcal{M}}}$

  4. (d)

Simulate the model $\overline{{\mathcal{M}}}$ separately for each remaining
unlabeled cell, using the data to calculate neighbor interactions

  5. (e)

Classify cells based on simulation error under $\overline{{\mathcal{M}}}$ and
label the lowest-error class as the new species $\overline{{\mathcal{S}}}$

  6. (f)

Remove cells in $\overline{{\mathcal{S}}}$ from the population and return to
step (b)

To summarize, the result is a set of $S$ models and species
$\\{(\overline{{\mathcal{M}}}_{\ell},\overline{{\mathcal{S}}}_{\ell})\\}_{\ell=1}^{S}$,
where each model $\overline{{\mathcal{M}}}_{\ell}$ is constructed from an
average of individual models within a cluster. We use the notation333See Table
5 for a complete list of notations used. of ${\mathcal{M}}_{i}$ to be the
model for the $i$-th cell, ${\mathcal{C}}_{j}$ to be the $j$-th cluster of
models, and $\overline{{\mathcal{S}}}_{\ell}$ to be the set of cells
identified as the $\ell$-th species and which obey model
$\overline{{\mathcal{M}}}_{\ell}$.

We now give more detail on steps (a)-(f), including stopping criteria, leaving
the more technical aspects to the appendix.

###  3.1 Learning single-cell models

The first step in the algorithm is to learn an ensemble of single-cell models,
one for each of the $N$ focal cell selected from the $N_{\text{tot}}$ total
cells tracked during the experiment444For the artificial data used here, the
set of focal cells is the entire population ($N=N_{tot}$), however we discuss
in the conclusion various situations in which a subset of $N<N_{\text{tot}}$
focal cells should be chosen.. By “single-cell” model, we mean that the scope
of each model is limited to learning only the dynamics of its focal cell,
however data from the remaining tracked cells are incorporated to learn the
interaction forces.

#### WSINDy

The main ingredient in learning single-cell models in ${\mathcal{M}}$ is the
WSINDy algorithm, together with careful choices for the bases used to
represent the three main forces $f_{\text{a-r}}$, $f_{\text{align}}$, and
$f_{\text{drag}}$. Each cell $i$ is identified by a position and velocity
$(x_{i}(t),v_{i}(t))$ which we assume is well-approximated by a 2nd-order
model of the form (2.1). The dynamics take the general form

(3.1) $\ddot{x}_{i}(t)=F_{i}(X(t),V(t)),$

where $(X(t),V(t))\in\mathbb{R}^{2dN_{\text{tot}}}$ denotes the entire
population of positions and velocities in the colony at time $t$. We then
assume that we have available a dataset of positions
$\mathbf{X}=(\mathbf{x}_{1}(t_{k}),\dots,\mathbf{x}_{N_{\text{tot}}}(t_{k}))_{k=1}^{L}$
sampled from the system $X$ at $L$ timepoints. Our goal is to identify $F_{i}$
using $\mathbf{X}$.

The SINDy approach involves representing $F_{i}$ as a sparse linear
combination of basis elements $\Theta(X,V):=(f_{j}(X,V))_{1\leq j\leq J}$,
such that

$F_{i,d^{\prime}}(X,V)=\sum_{j=1}^{J}\mathbf{w}^{\star}_{i,j}f_{j,d^{\prime}}(X,V),$

where subscript $d^{\prime}$ indicates the spatial coordinate
($d^{\prime}\in\\{1,2\\}$ in this study). The basis $\Theta$ is chosen by the
user and determines the accuracy of the learned model as well as the
conditioning of the WSINDy algorithm.

The available cell position data $\mathbf{X}$ is used to approximate
velocities $\mathbf{V}:=\dot{\mathbf{X}}\approx\dot{X}$ and accelerations
$\ddot{\mathbf{X}}\approx\ddot{X}$, using e.g. finite differences, leading
from (3.1) to the data-driven linear system

(3.2)
$\ddot{\mathbf{x}}_{i}\approx\Theta(\mathbf{X},\mathbf{V})\mathbf{w}^{\star}_{i}.$

With some abuse of notation, we denote by $\Theta(\mathbf{X},\mathbf{V})$ the
matrix that results from evaluating the basis $\Theta(X,V)$ at the time-series
data $(\mathbf{X},\mathbf{V})$. The entries are
$\Theta(\mathbf{X},\mathbf{V})_{k+(d^{\prime}-1),j}=f_{j,d^{\prime}}(\mathbf{X}(t_{k}),\mathbf{V}(t_{k}))$.

The data $\mathbf{X}$ is often corrupted by measurement noise, which leads to
inaccurate computations of derivatives $\mathbf{V}$. For the current setting,
the standard SINDy approach just outlined requires second-order derivatives
$\ddot{\mathbf{X}}$, which are even less accurate to compute from noisy data.
To prevent some of the corruption from noise, we can use the weak form, which
leads to WSINDy. Returning to equation (3.1), we convert to the weak form by
multiplying by a test function $\phi(t)$ and integrating in time,

(3.3)
$\left\langle\phi,\ddot{x}_{i}\right\rangle:=\left\langle\phi,F_{i}(X,V)\right\rangle.$

where the inner product $\left\langle\cdot,\cdot\right\rangle$ denotes the
time intergral

$\left\langle f,g\right\rangle=\int_{-\infty}^{\infty}f(t)g(t)dt.$

Choosing $\phi$ to be twice differentiable and zero outside of some interval
$(a,b)$, we then integrate by parts twice on the left-hand side to arrive at

$\langle\ddot{\phi},x_{i}\rangle=\left\langle\phi,F_{i}(X,V)\right\rangle,$

so that the second derivative has been removed from $x_{i}$ and placed on
$\phi$. Choosing a basis of test functions $\Phi:=(\phi_{q})_{1\leq q\leq Q}$,
we build the weak-form linear system

(3.4) $\mathbf{b}^{(i)}\approx\mathbf{G}^{(i)}{\widehat{\mathbf{w}}}^{(i)}$

where
$\mathbf{b}^{(i)}_{q+(d^{\prime}-1)}=\langle\ddot{\phi}_{q},\mathbf{x}_{i,d^{\prime}}\rangle$
and
$\mathbf{G}^{(i)}_{q+(d^{\prime}-1),j}=\left\langle\phi_{q},f_{j,d^{\prime}}(\mathbf{X},\mathbf{V})\right\rangle$.

As well as choosing $\Theta$, in order to compute
$(\mathbf{G}^{(i)},\mathbf{b}^{(i)})$, we need to compute555In some cases we
can use weak-form integration to eliminate this step but we do not pursue this
in the current work. $\mathbf{V}$ from the position data $\mathbf{X}$, choose
a test function basis $\Phi$, and discretize integrals appearing in the linear
system. For simplicity, we compute $\mathbf{V}$ using 2nd-order centered
finite difference, although a number of methods exist for numerical
differentiation from data [29, 64]. For integration, we use the trapezoidal
rule, and we use test functions of the form

(3.5) $\phi_{q}(t)=\max\left(1-\left(\frac{t-t_{q}}{m\Delta t}\right)^{2},\
0\right)^{p},$

for shape parameters $m$ and $p$, and timestamps $t_{q}$ in the range of the
available time series. We refer to [38, 37] for methods of choosing
$(m,p,t_{q})$ and for accuracy and robustness results concerning the
trapezoidal rule combined with this particular class of test functions. In
this work, we use the changepoint algorithm in [37] with $\tau=10^{-10}$ and
$\widehat{\tau}=3$, leading to $m\in\\{31,\dots,38\\}$ and $p\in\\{8,9\\}$
(Table 6 lists these values used in the examples below). Since the time series
below are relative short ($L=200$ or $L=400$), we use all available $t_{q}$,
i.e. $(t_{q})_{q=1}^{Q}=(m\Delta t,\dots,(L-m-1)\Delta t)$ so that $Q=L-2m$.

#### Regression

We solve the linear system (3.4) for coefficients
${\widehat{\mathbf{w}}}^{(i)}$ by approximately solving the following
constrained sparse regression problem:

(3.6) ${\widehat{\mathbf{w}}}^{(i)}=\operatorname*{arg\,min}_{\mathbf{w}\
s.t.\
\mathbf{C}\mathbf{w}\leq\mathbf{d}}\left\\{\left\|{\mathbf{G}^{(i)}\mathbf{w}-\mathbf{b}^{(i)}}\right\|_{2}^{2}+\lambda^{2}\left\|{\mathbf{w}}\right\|_{0}\right\\}.$

The linear inequality constraint $\mathbf{C}\mathbf{w}\leq\mathbf{d}$ encodes
the constraints listed in (2.2), (2.3), and (2.4) on the forces on
$f_{\text{a-r}}$, $f_{\text{align}}$, and $f_{\text{drag}}$, and $\lambda$ is
the sparsity threshold. We employ the modified sequential thresholding
algorithm from [37, 36], with least-squares iterations replaced by solving the
associated linearly-constrained quadratic program666Details on this
implementation, in particlar $\mathbf{C}$ and $\mathbf{d}$, are given in
Appendix C..

Since the coefficients ${\widehat{\mathbf{w}}}^{(i)}$ have no a priori
absolute magnitude, we threshold only on the magnitudes of the given term
relative to the response vector $\mathbf{b}^{(i)}$, namely, we define the
thresholding operator $H_{\lambda}(\mathbf{w})$ by

(3.7)
$H_{\lambda}(\mathbf{w}_{j})=\begin{dcases}0,&\frac{\left\|{\mathbf{G}^{(i)}_{j}\mathbf{w}_{j}}\right\|}{\left\|{\mathbf{b}^{(i)}}\right\|}\notin[\lambda,\
\lambda^{-1}]\\\ \mathbf{w}_{j},&\text{otherwise}.\end{dcases}$

The constrained sequential thresholding algorithm then consists of alternate
between (I) solving (3.6) with $\lambda=0$ for $\widetilde{\mathbf{w}}$ with
the additional constraint
$\text{supp}\left(\widetilde{\mathbf{w}}\right)\subset\text{supp}\left(\mathbf{w}_{\ell}\right)$,
and (II) setting $\mathbf{w}^{(\ell+1)}=H_{\lambda}(\widetilde{\mathbf{w}})$.
A sweep over $\lambda$ values in the range $10^{-4}$ to $10^{0}$ is performed
according to [38, 36] to choose an appropriate threshold $\lambda$.

#### Trial Basis Functions

In the case of the directional force model (2.1), we require three bases
${\mathcal{F}}_{\text{a-r}}$, ${\mathcal{F}}_{\text{align}}$, and
${\mathcal{F}}_{\text{drag}}$ for the three proposed forces $f_{\text{a-r}}$,
$f_{\text{align}}$, and $f_{\text{drag}}$. We seek a sparse model, and so
choose global basis functions, rather than a model composed of a large sum
over basis functions that are spatially localized.

For the attractive-repulsive basis ${\mathcal{F}}_{\text{a-r}}$ we choose
products of cosines and scaled and weighted Laguerre polynomials,

${\mathcal{F}}_{\text{a-r}}=\\{\cos(n\theta)p_{\ell}(\alpha
r)e^{-\frac{\alpha}{2}r}\\}_{n=0,\ell=0}^{2,17}$

for $\ell$th degree Laguerre polynomial $p_{\ell}$. The scale $\alpha$ is
chosen from $r_{\max}$, the maximum observed distance between cells, such that
$e^{-\frac{\alpha}{2}r_{\max}}=\epsilon_{\text{mach}}\approx e^{-36}$. We set
$\alpha=36$ in all cases below since $r_{\max}\approx 2$.

The pattern of attractive and repulsive regions of the force $f_{\text{a-r}}$
is not known a priori, hence the Laguerre basis offers flexibility. The choice
of weighted Laguerre polynomials (with weight $\omega(r)=e^{-r/2}$) is guided
by the orthogonality relation

$\int_{0}^{\infty}p_{m}(r)p_{n}(r)\omega^{2}(r)dr=\delta_{mn}.$

We find ${\mathcal{F}}_{\text{a-r}}$ leads to a well-conditioned matrix
$\mathbf{G}$ despite orthogonality not holding with respect to the data
distribution.

For the alignment force we choose a basis of shifted cosines and exponential
functions,

${\mathcal{F}}_{\text{align}}=\\{(1+\cos(n\theta))e^{-2^{\ell}r}\\}_{n=0,\ell=-2}^{2,5}$

which is informed by the fact that $f_{\text{align}}$ must be negative. This
is easily controlled with ${\mathcal{F}}_{\text{align}}$ by simply enforcing
that the coefficients ${\widehat{\mathbf{w}}}_{\text{align}}$ be negative. For
the same reason, we choose the drag force from a basis of monomials and
cosines,

${\mathcal{F}}_{\text{drag}}=\\{(1+\cos(n\theta))|v|^{\ell}\\}_{n=0,\ell=0}^{1,4}$

as this can also be easily controlled to yield an overall negative
$f_{\text{drag}}$ force by constraining only the basis elements
${\widehat{\mathbf{w}}}_{\text{drag}}$. Moreover, monomials capture the
physical assumption that resistance to motion should increase with speed.

###  3.2 Cluster

Once the ensemble of models ${\mathcal{M}}$ has been learned, each model is
validated on a small neighborhood of cells, and models with poor performance
are replaced by a neighboring model with better performance, if one exists.
This step can be seen as a form of “cross-pollination” that greatly increases
the chance of finding accurate models and eliminating inaccurate models.
Details on this step can be found in Appendix D.

Models are then clustered according to the force modes present. Specifically,
using the bases above, we can expand each force according to distinct
directional modes:

$\displaystyle f_{\text{a-r}}(r,\theta)$
$\displaystyle=f^{(0)}_{\text{a-r}}(r)+\cos(\theta)f^{(1)}_{\text{a-r}}(r)+\cos(2\theta)f^{(2)}_{\text{a-r}}(r)$
$\displaystyle f_{\text{align}}(r,\theta)$
$\displaystyle=f^{(0)}_{\text{align}}(r)+\cos(\theta)f^{(1)}_{\text{align}}(r)+\cos(2\theta)f^{(2)}_{\text{align}}(r)$
$\displaystyle f_{\text{drag}}(|v|,\theta)$
$\displaystyle=f^{(0)}_{\text{drag}}(|v|)+\cos(\theta)f^{(1)}_{\text{drag}}(|v|)$

This leads to eight possible force modes, which we order as follows:

$\left\\{f^{(0)}_{\text{a-r}},f^{(1)}_{\text{a-r}},f^{(2)}_{\text{a-r}},f^{(0)}_{\text{align}},f^{(1)}_{\text{align}},f^{(2)}_{\text{align}},f^{(0)}_{\text{drag}},f^{(1)}_{\text{drag}}\right\\}.$

We associate the sparsity pattern of the force modes with the set of all 8-bit
codes, giving a total of $2^{8}=256$ possible model clusters. Models are
partitioned into clusters
${\mathcal{C}}=\\{{\mathcal{C}}_{1},\dots,{\mathcal{C}}_{r}\\}$ based on their
associated codes. For example, species A listed in Table 1 below is associated
with the code $10111010$ indicating that $f^{(0)}_{\text{a-r}}$,
$f^{(2)}_{\text{a-r}}$, $f^{(0)}_{\text{align}}$, $f^{(1)}_{\text{align}}$,
and $f^{(0)}_{\text{drag}}$ are present in the model.

There are several other options for model replacement and clustering, include
clustering based on the sparsity pattern of ${\widehat{\mathbf{w}}}$, or
simply on the presence or absence of each of the three forces
$f_{\text{a-r}}$, $f_{\text{align}}$, $f_{\text{drag}}$. The former
significantly increases the number of possible clusters, while the latter
leads to just 8 possible clusters. Our choice reflects the desire to
disentangle directionalities of the governing forces without introducing a
strong dependence on the bases used to approximate each force.

###  3.3 Aggregate

Having formed the model clusters ${\mathcal{C}}$, let
$\overline{{\mathcal{C}}}$ be the cluster with the most members. We then
compute $\overline{{\mathcal{M}}}$ by averaging over the coefficients from
models in $\overline{{\mathcal{C}}}$, that is, we compute

(3.8)
$\overline{\mathbf{w}}=\frac{1}{|\overline{{\mathcal{C}}}|}\sum_{i\in\overline{{\mathcal{C}}}}{\widehat{\mathbf{w}}}^{(i)},$

where $|\overline{{\mathcal{C}}}|$ denotes the number of elements in
$\overline{{\mathcal{C}}}$. This ensemble average preserves the force modes
that identify $\overline{{\mathcal{C}}}$. As explored in [22], in some cases
it may be more appropriate to use the coefficient median, or take a weighted
average. We leave these for future work.

###  3.4 Validate

To validate the aggregate model $\overline{{\mathcal{M}}}$, we perform forward
simulations over the remaining unclassified cells in a highly parallelizable
way that utilizes the experimental data to efficiently march forward in time.

Let $N^{\prime}\leq N$ be the number of remaining unclassified cells. For each
$i=1,\dots,N^{\prime}$, we simulate a new trajectory
$\\{(\overline{x}_{i}(t_{k}),\overline{v}_{i}(t_{k}))\\}_{k=1}^{L}$ using the
averaged model $\overline{{\mathcal{M}}}$ with the experimental initial
conditions
$(\overline{x}_{i}(0),\overline{v}_{i}(0))=(\mathbf{x}_{i}(0),\mathbf{v}_{i}(0))$.
We march forward in time according to the forward Euler update:

(3.9) $\displaystyle\overline{v}_{i}(t_{k+1})$
$\displaystyle=\overline{v}_{i}(t_{k})+\Delta
t\overline{{\mathcal{M}}}(\overline{x}_{i}(t_{k}),\overline{v}_{i}(t_{k}),\mathbf{X}^{\prime
i}(t_{k}),\mathbf{V}^{\prime i}(t_{k}))$ (3.10)
$\displaystyle\overline{x}_{i}(t_{k+1})$
$\displaystyle=\overline{x}_{i}(t_{k})+\Delta t\overline{v}_{i}(t_{k})$

where $(\mathbf{X}^{\prime i}(t_{k}),\mathbf{V}^{\prime i}(t_{k}))$ indicates
$(\mathbf{X}(t_{k}),\mathbf{V}(t_{k}))$ with the $i$-th cell removed777By some
abuse of notation, $\overline{{\mathcal{M}}}(x,v,X,V)$ is used to denote the
instantaneous force on a particle $(x,v)$ from neighboring cells $(X,V)$ under
model $\overline{{\mathcal{M}}}$. Since the time resolution of the data
$\Delta t$ is assumed to be coarse, we perform the simulation on a finer grid
with timestep $\Delta t^{\prime}=2^{-5}\Delta t$, and use piecewise cubic
hermite interpolation to generate positions and velocities of neighbors
$(\mathbf{X}^{\prime i},\mathbf{V}^{\prime i})$ at intermediate times. We
stress that we do not update the neighbor cells using the model
$\overline{{\mathcal{M}}}$, which would be much more costly, we merely use
neighbor positions and velocities from the data to compute interactions that
govern the motion of cell $i$. The resulting trajectories
$\\{(\overline{x}_{i},\overline{v}_{i})\\}_{i=1}^{N^{\prime}}$ can then be
trivially computed in parallel888See Table 7 for simulation walltimes..

We then define the validation error for cells $i=1,\dots,N^{\prime}$ as the
relative velocity difference

(3.11) $\Delta
V_{i}:=\sqrt{\frac{\sum_{k=1}^{L^{\prime}}\left\|{\overline{v}_{i}(t_{k})-\mathbf{v}_{i}(t_{k})}\right\|_{2}^{2}}{\sum_{k=1}^{L^{\prime}}\left\|{\mathbf{v}_{i}(t_{k})}\right\|_{2}^{2}}},$

where $L^{\prime}\leq L$ is a subset of the time series over which the
simulation is expected to remain accurate999In particular, for chaotic systems
the trajectories cannot be expected to remain close for large times, however
the correct model will be initially accurate.. In this work we choose
$L^{\prime}=0.25L$. With $L=200$ timesteps in the examples below, this results
in a comparison with the data over the first 50 timesteps at the original
scale $\Delta t$, or equivalently 1600 timesteps on the finer scale $\Delta
t^{\prime}$.

###  3.5 Classify

Let $VE$ be the set of validation errors, $VE=\\{\Delta V_{1},\dots,\Delta
V_{N^{\prime}}\\}$. An empirical observation used in this work is that when
$\overline{{\mathcal{M}}}$ approximates well an underlying model for a true
species, the log-transformed validation errors $\log_{10}(VE)$ are fit well by
a Gaussian mixture model (GMM) with two mixtures (see Figures 4.4,4.7,4.9). We
thus use a 2-mixture GMM to classify the remaining cells. Cells are granted
membership into the mixture that yields the highest posterior probability of
generating its log-validation error, and the class with lowest mean error is
labeled as a species. This can thought of as a sequential binary
classification scheme.

For example, in each plot of Figure 4.4, a representative GMM resulting from a
two-species population, the left-most mixture corresponds to low validation
errors under the model $\overline{{\mathcal{M}}}$ and is classified as a
species $\overline{{\mathcal{S}}}$ (in this case, species C in Table 1 is
identified). The cells in $\overline{{\mathcal{S}}}$ are subsequently removed
from the population, and cells in the right-most mixture are returned to the
cluster phase (b).

###  3.6 Stopping conditions

Steps (b)-(e) are repeated until one the following conditions is met.

  1. (1)

Less than 3 cells remain: $N^{\prime}\leq 2$

  2. (2)

More than $99\%$ of remaining cells have less than $5\%$ validation error:
$\\#\\{VE<0.05\\}\geq 0.99N^{\prime}$

  3. (3)

More than 10 species have been identified: $S>10$

The first case is an obvious criterion to prevent infinite looping over
outlier cells for which there is not enough information to learn an adequate
model. The second condition skips the GMM fitting process when all cells have
sufficiently low error101010We choose this cutoff at $5\%$ but note that this
parameter may need tuning., and directly assigns all unlabeled cells to a new
species. This is necessary to account for the cases of high-accuracy recovery,
where it is observed that the $VE$ is no longer approximately log-normal. For
example, this is observed in the rightmost plot of Figure 4.3. Finally, for
$N$ very large, it may be necessary to restrict the total number of species,
which is encapsulated in the third condition, although we did not observe the
number of iterations exceeding 5 in any trials with $N\leq 1000$ and $S\leq 3$
true species.

##  4 Results: Artificial Cells

|
---|---
Figure 4.1: Forces used to generate artificial data, motivated by experiment.
Left: quadrupolar attractive-repulsive force $f^{\star}_{\text{a-r}}$. Right:
dipolar alignment force $f^{\star}_{\text{align}}$. Note that all artificial
species are prescribed the same drag force
$f^{\star}_{\text{drag}}(s,\theta)=-5s$, linear in cell speed $s$.

(4.1)
$\boxed{\begin{array}[]{rl}f_{\text{a-r}}^{\star}(r,\theta)&:=(15+10\cos(2\theta))\left(e^{-20r}-0.25e^{-10r}\right)\\\
f_{\text{align}}^{\star}(r,\theta)&:=-(8+8\cos(\theta))e^{-8r}\\\
f_{\text{drag}}^{\star}(s,\theta)&:=-5s\end{array}}$

| $f^{\star}_{\text{a-r}}$ | $f^{\star}_{\text{align}}$ | $f^{\star}_{\text{drag}}$
---|---|---|---
Species A | $\checkmark$ | $\checkmark$ | $\checkmark$
Species B | $\checkmark$ |  | $\checkmark$
Species C |  | $\checkmark$ | $\checkmark$

Table 1: Species delineation by active force modes.

We examine artificial cell cultures with combinations of 1-3 distinct cell
types, denoted by species A, species B, and species C. Each species has a
unique combination of the forces depicted in Figure 4.1, with the combinations
and force definitions specified in Table 1. The forces include a quadrupolar
attractive-repulsive force $f^{\star}_{\text{a-r}}$, a dipolar alignment force
$f^{\star}_{\text{align}}$, and a monopolar drag force
$f_{\text{drag}}^{\star}$ which is linear in its speed argument. As we will
see below, species A and species B share the force $f_{\text{a-r}}$ and hence
result in similar dynamics, which presents a challenge to identification. It
turns out that using a longer time series results in correct classification.

We let $\mathbf{X}_{P}$ denote a simulation with individuals from species $P$,
for example, $\mathbf{X}_{A}$ is a simulation with only individuals from
species A and $\mathbf{X}_{A,B}$ is a simulation with a mixed population of
species A and species B. Each simulation has 1000 individuals and an even
number of members in each species (up to rounding). More details on the
simulations can be found in the appendix.

We are particularly interested in three traits of our learning algorithm: (1)
was the classification successful? (2) are the learned forces close to the
true forces? (3) are simulated trajectories using the learned model close to
the original trajectories? To assess (1) we report the classification success
CS$(i)$ for $i\in\\{A,B,C\\}$ as the fraction of individual from species $i$
that ended up in the cluster in question, where clusters are listed as subrows
(rows not separated by horizontal lines) within each row in Tables 2, 3, and
4, in the order they were identified in step (e) of the algorithm. For
example, in row 2 of Table 3, two clusters are identified from the two-species
data $\mathbf{X}_{A,C}$, with the first cluster containing $100\%$ of the
species C cells, indicated by CS$(C)=1.000$, and the second cluster containing
$100\%$ of the species A cells, indicated by CS$(A)=1.000$, with no outliers.

To assess the accuracy of learned forces with respect to the ground truth
forces, for each of the three forces we report the relative $L^{2}$ error over
$r\in[0,2]$ (since $r_{\max}\approx 2$ for all examples here) and
$\theta\in\\{0,\pi/4,\pi/2,3\pi/4,\pi\\}$, denoting this by $\Delta
f_{\text{a-r}}$, $\Delta f_{\text{align}}$, and $\Delta f_{\text{drag}}$.

Lastly, we assess the difference in learned and true trajectories using the
average validation error $\Delta
V=\frac{1}{|\overline{{\mathcal{S}}}|}\sum_{i\in\overline{{\mathcal{S}}}}\Delta
V_{i}$, where $\Delta V_{i}$ is computed from model $\overline{{\mathcal{M}}}$
associated with identified species $\overline{{\mathcal{S}}}$ using (3.11).

###  4.1 Homogeneous Populations

Experiment | $\Delta f_{\text{a-r}}$ | $\Delta f_{\text{align}}$ | $\Delta f_{\text{drag}}$ | CS($A$) | CS($B$) | CS($C$) | $\Delta V$
---|---|---|---|---|---|---|---
$\mathbf{X}_{A}$ | $0.0211$ | $0.0384$ | $0.0382$ | 1.000 | — | — | 0.0100
$\mathbf{X}_{B}$ | $0.0098$ | — | $0.0125$ | — | 0.997 | — | 0.0076
$\mathbf{X}_{C}$ | — | $0.0011$ | $0.0016$ | — | — | 1.000 | 0.0016
Table 2: Performance of model learning and classification algorithm of homogeneous populations. |  |
---|---|---
Figure 4.2: Examples of learned and original trajectories from homogeneous
populations. Left to right: $\mathbf{X}_{A}$, $\mathbf{X}_{B}$,
$\mathbf{X}_{C}$.

As an initial benchmark we detect single-species populations from homogeneous
data. While simpler than the heterogeneous case, this is a nontrivial task due
to the variability of single-cell trajectories and local environments within
the population. Our method successfully identifies the models for species $A$,
$B$, and $C$ from homogeneous simulations, achieving less than $1\%$ mean
validation errors in each case, and less than $4\%$ relative force errors
$\Delta f$ (Table 2). In simulation $\mathbf{X}_{B}$, three cells are
identified as outliers (appearing in the right tail of Figure 4.3 (middle)),
and all other cells in $\mathbf{X}_{A}$, $\mathbf{X}_{B}$, and
$\mathbf{X}_{C}$ are correctly classified. A comparison between original and
learned trajectories is depicted in Figure 4.2, with learned trajectories
overlapping original trajectories in each case.

|  |
---|---|---
Figure 4.3: Distribution of log-validation errors for homogeneous cell
experiments $\mathbf{X}_{A}$, $\mathbf{X}_{B}$, $\mathbf{X}_{C}$.
Distributions for $\mathbf{X}_{A}$ and $\mathbf{X}_{B}$ are fit well by a
single Gaussian, indicating a single species is present. The distribution for
$\mathbf{X}_{C}$ has a non-Gaussian tail, although all errors are below $1\%$,
indicating that the candidate model fits the population up to the specified
error tolerance.

###  4.2 Two-Species Populations

Experiment | $\Delta f_{\text{a-r}}$ | $\Delta f_{\text{align}}$ | $\Delta f_{\text{drag}}$ | CS($A$) | CS($B$) | CS($C$) | $\Delta V$
---|---|---|---|---|---|---|---
$\mathbf{X}_{A,C}$ |  | —
---
$0.0226$
| $0.0011$
---
$0.0941$
| $0.0002$
---
$0.0472$
| 0
---
1.000
| —
---
—
| 1.000
---
0
| 0.0005
---
0.0200
$\mathbf{X}_{B,C}$ |  | —
---
$0.0328$
| $0.0077$
---
—
| $0.0051$
---
$0.0461$
| —
---
—
| 0
---
1.000
| 1.000
---
0
| 0.0023
---
0.0339
$\mathbf{X}_{A,B}$ |  | $0.0341$
---
$0.0075$
$0.4780$
| —
---
—
$0.3660$
| $0.0133$
---
$0.0538$
—
| 0.102
---
0.084
0.814
| 0.628
---
0.340
0.012
| —
---
—
—
| 0.1201
---
0.0362
0.3945
$\mathbf{X}_{A,B}$(long) |  | $0.0018$
---
$0.0034$
$0.0071$
$0.0034$
| —
---
0.0023
—
0.0051
| 0.0042
---
0.0067
0.0044
0.0144
| 0.002
---
0.994
0
0.004
| 0.978
---
0
0.023
0
| —
---
—
—
—
| 0.0045
---
0.0070
0.0568
0.0199
Table 3: Performance of model learning and classification algorithm for two-
species populations. $\mathbf{X}_{A,B}$(long) is simply the continuation of
$\mathbf{X}_{A,B}$ to twice the time horizon, and significantly improves
classification over $\mathbf{X}_{A,B}$.

Next we examine the ability of the learning algorithm to detect two-species
populations along with accurate aggregate models. Figure 4.4 displays two
representative Gaussian mixture fits to the log-validation errors for
$\mathbf{X}_{A,C}$ (left) and $\mathbf{X}_{B,C}$ (right). In both cases, the
log-errors are well-approximated by Gaussian mixtures with wide separations
between mixtures. This allows for complete classification in both cases, as
indicated by CS$(A)$, CS$(B)$, and CS$(C)$ in rows 2 and 3 of Table 3. Force
differences $\Delta f$ are less than $5\%$ in each case, with trajectory
validation errors less than $3.5\%$. In particular, species C achieves less
that $0.3\%$ validation error, which is due to the true force
$f_{\text{align}}^{\star}$ existing in the learning library, whereas
$f_{\text{a-r}}^{\star}$ is approximated using a truncated series expansion,
resulting in larger errors. See Figures 4.5 and 4.6 for comparison between
original and learned trajectories.

For experiment $\mathbf{X}_{A,B}$, initially the method is incapable of
correctly classifying cells into species A and species B. Three clusters are
identified with sub-optimal models (Table 3 row 4). Accurate classification is
achieved by running the algorithm with a longer experiment
$\mathbf{X}_{A,B}$(long) (Table 3 row 5) which is the continuation of
$\mathbf{X}_{A,B}$ for twice the total time points, at the same temporal
resolution. An initial cluster is identified containing $97.8\%$ of species
$B$ along with $0.2\%$ (a single cell) of the existing species $A$ cells,
followed by a second cluster with $99.4\%$ of the species $A$ cells and no
cells from species $B$. The last two clusters correctly partition the
remaining cells (12 in total), again finding accurate models, allowing for
recombination with the first two cluster during post-processing.

Figure 4.8 shows a comparison between original and learned trajectories for
$\mathbf{X}_{A,B}$(long) and Figure 4.7 depicts representative Gaussian
mixture models. In particular, Figure 4.7 (left) shows increased overlap
between the two Gaussian mixtures in the first iteration, compared to Figure
4.4, however model performance is still sufficiently different as to classify
$\approx 98\%$ of cells correctly.

|
---|---
Figure 4.4: Distribution of log-validation errors for heterogeneous cell experiments $\mathbf{X}_{A,C}$ and $\mathbf{X}_{B,C}$. In each case, species $C$ is identified in the first iteration, and a clear separation between the two species allows for accurate clustering using Gaussian mixture models. |  |
---|---|---
Figure 4.5: Example trajectories from experiment $\mathbf{X}_{A,C}$. Cells with true color labels are depicted on the left, but are passed into the algorithm unlabeled. The algorithm then classifies the population into different species and returns accurate models for each species. Classified cells from species $C$ (middle) and species $A$ (right) are highlighted showing excellent agreement between data and simulation. |  |
---|---|---
Figure 4.6: Example trajectories from experiment $\mathbf{X}_{B,C}$. As in Figure 4.5, cells with true color labels are depicted on the left. Classified cells from species $C$ (middle) and species $B$ (right) are highlighted showing excellent agreement original data and output of the learned models. |
---|---
Figure 4.7: Log-validation errors for two-species population $\mathbf{X}_{A,B}$(long). Strong similarities between the two species present an initial challenge to identification, which is overcome by taking a longer time series. The initial Gaussian mixture model (left) identifies a majority species B cluster. In the second iteration (right), a cluster with all Species A cells is identified, and a small group of cells remains which is then partitioned correctly (see row 5, columns CS(A) and CS(B) of Table 3). |  |
---|---|---
Figure 4.8: Example trajectories from experiment $\mathbf{X}_{A,B}$(long).
Cells with true labels are depicted on the left and classified cells from
species $B$ (middle) are species $A$ (right) are depicted with model output
overlapping the input data in each case.

###  4.3 Three-Species Population

As a final test we identify species from the three-species experiment
$\mathbf{X}_{A,B,C}$. Similar to the case $\mathbf{X}_{A,B}$, we see
improvements with a longer time-series $\mathbf{X}_{A,B,C}$(long). For the
initial experiment $\mathbf{X}_{A,B,C}$, species C is completely identified in
the first cluster (Table 4 row 2), and in the second cluster $98.8\%$ of
species $B$ cells are identified along with $9.1\%$ of species $A$ cells,
leading to a fairly inaccurate model ($\Delta V\approx 0.06$). The subsequent
clusters divide the remaining species $A$ and $B$ cells.

Doubling the time series with $\mathbf{X}_{A,B,C}$(long), we find the majority
of each species residing in its own cluster (Table 4 row 3). Cluster 1
contains all of the species C cells, cluster 2 consists of 94% of species B
cells and 1.8% of species A, and cluster 4 consists of 95% of species A.
Moreover, the aggregrate models for each of these clusters result in
validation errors under $1\%$.

Clusters 4 and 5 of $\mathbf{X}_{A,B,C}$(long) are each small mixtures of
species A and species B, however the learned forces in each case are
nevertheless accurate. We find that $\Delta f\leq 5\%$ (row 3 of Table 4,
subrows 4 and 5) while the validation error is fairly high, with $\Delta V\geq
6.6\%$. This indicates that the cells in clusters 4 and 5 have trajectories
that are particularly sensitive to perturbations. Given the complexity of the
dynamics (one can observe sharp turns taken by cells in the bottom two plots
of Figure 4.10), trajectories cannot be expected to remain close for all time,
and in this case the validation error (3.11) may be too strong a
metric111111We discuss this further in the conclusion but leave a complete
investigation to future work.. It is remarkable that the aggregate models for
clusters 1, 2, and 3 produce accurate learned trajectories for the whole time
series.

Experiment | $\Delta f_{\text{a-r}}$ | $\Delta f_{\text{align}}$ | $\Delta f_{\text{drag}}$ | CS($A$) | CS($B$) | CS($C$) | $\Delta V$
---|---|---|---|---|---|---|---
$\mathbf{X}_{A,B,C}$ |  | —
---
$0.0287$
0.0243
0.0022
0.0478
| 0.0005
---
—
0.2321
0.0001
—
| 0.0020
---
0.0191
—
0.0050
0.0118
| 0
---
0.091
0.542
0.332
0.009
| 0
---
0.988
0.006
0
0.006
| 1.000
---
0
0
0
0
| 0.0016
---
0.0620
0.5010
0.0406
0.3462
$\mathbf{X}_{A,B,C}$(long) |  | —
---
$0.0117$
0.0031
0.0090
0.0082
| 0.0009
---
—
0.0034
—
0.0303
| 0.0043
---
0.0048
0.0041
0.0011
0.0404
| 0
---
0.018
0.949
0.003
0.040
| 0
---
0.940
0
0.06
0
| 1.000
---
0
0
0
0
| 0.0011
---
0.0095
0.0094
0.0663
0.0951
Table 4: Performance of model learning and classification algorithm of three-species populations. $\mathbf{X}_{A,B,C}$(long) is simply the continuation of $\mathbf{X}_{A,B,C}$ to twice the time horizon. |  |
---|---|---
Figure 4.9: Gaussian mixture models for classifying the three-species experiment $\mathbf{X}_{A,B,C}$(long) (see Table 4 row 3 for details). We see an initial complete separation of species $C$ (left), followed by a mixed cluster containing $94\%$ of the species B cells and $1.8\%$ of the species $A$ cells (middle). The next iteration classifies a majority species A cluster (right). Clusters 4 and 5 are effectively outliers and contain the remaining 31/1000 cells. |
---|---
|
Figure 4.10: Example trajectories simulated using learned models from
$\mathbf{X}_{A,B,C}$(long). Top left: example cells with true labels. Top
right, bottom left, bottom right: example trajectories from clusters 1-3 (see
Table 4 row 3 for details). Note in particular that learned models for
clusters 2 and 3 are able to capture sharp turns in the true dynamics.

##  5 Discussion

We have introduced a method for performing a combined classification and model
selection task relevant to heterogeneous systems of autonomous agents.
Specifically, we have shown that learning an ensemble of interacting particle
models (one for each agent) allows iterative classification of agents into
species according to their forward simulation accuracy. This is surprising due
to the limited information carried in a single trajectory. Fortunately, the
validation errors empirically approximate a log-normal distribution, hence
Gaussian mixture model classification arises as the appropriate tool for
identifying species membership.

Computational feasibility of this approach is grounded in the parallel nature
of both the learning algorithm and the simulation component. Learning each
single-trajectory model is cheap, with linear systems of size $m\times n$ with
$m$ and $n$ not exceeding several hundreds (see Table 7 for walltimes). In the
simulation step, each trajectory is validated separately and in parallel, and
neighbor interactions are computed using the measurement data itself,
resulting in ${\mathcal{O}}(N)$ pairwise interaction computations per timestep
instead of an ${\mathcal{O}}(N^{2})$ full simulation121212We are not aware of
other simulation approaches in the context of interacting particles that use
the measurement data for direct calculation of the forces..

Several aspects of this approach deserve a more technical analysis, which may
lead to improvements. We aim in future work to undertake more rigorous study
of each component of the algorithm as outlined below.

#### Learning single-cell models

  1. (I)

Information content: It would be beneficial to quantify the information
content in each cell trajectory, possibly eliminating trajectories that do not
provide sufficient information. We saw in the experiments $\mathbf{X}_{A,B}$
and $\mathbf{X}_{A,B,C}$ that increasing the length of the trajectory leads to
better classification when two species exhibit similar dynamics. It may be
possible to use existing techniques, such as force matching, to identify
highly-informative cells from forces magnitudes, neighbor distributions, etc.
Model replacement, as outlined in Appendix D, is an initial step in this
direction.

  2. (II)

Noisy trajectories: While we expect measurement noise to be filtered out by
the weak form to similar extent demonstrated on ODE systems [38], we leave
examination of the robustness to both intrinsic (e.g. Brownian) and extrinsic
(e.g. measurement) noise to future work.

  3. (III)

Model library: We chose the force bases and constraints to reflect physical
properties, namely short-distance repulsion, long-distance decay, negative
alignment, and negative drag. Directional modes enforce bilateral symmetry,
and are low order (monopole, dipole, quadrupole). These can easily be adapted
to incorporate other known information, moreover the bases themselves may be
adapted to the data (note this is partially done, using the neighbor distance
distribution $\rho_{rr}$ to restrict the range of interactions). One major
assumption here is that there is no propulsion force, that energy is increased
only through anisotropic interactions with neighbors. It would be interesting
to examine whether this assumption holds true.

  4. (IV)

Regression approach: We employ modified sequential thresholding, which looks
for an overall sparse solution, although we threshold only on the term
magnitude $\left\|{\mathbf{G}_{j}\mathbf{w}_{j}}\right\|$ and not that raw
coefficient $\mathbf{w}_{j}$. This in particular allows
$f_{\text{a-r}},f_{\text{align}},f_{\text{drag}}$ to have equal opportunity to
enter the model despite different scales and bases used. The effect is that
the resulting model is sparse in the force modes (as described in Section
3.2), while each force mode may have many components (in fact $f_{\text{a-r}}$
is usually not sparse on a given directional mode). It may be more appropriate
to use a group sparsity-enforcing method, such as constrained group LASSO. In
general, this lies at the intersection of approximation and selection, where
sparse selection is required to select the correct modes, however the force
content in each mode requires approximation. Explorations of the appropriate
balance between selection and approximation would be valuable.

#### Cluster & Aggregate

  1. (I)

Here we cluster models based on the directional modes present. This can easily
be extended to the full pattern of nonzero elements in the model vector
${\widehat{\mathbf{w}}}$, although this depends on the number of resulting
clusters.

  2. (II)

We have used a simple uniformly-weighted average (3.8) to aggregate models,
however the use of model replacement (see Appendix D) implies that the average
is implicitly weighted according to model generalizability. Results may be
improved if other criteria (e.g. information criteria) are incorporated into
the weighted average, or if the median is taken instead of the mean.

  3. (III)

In the examples above, the aggregate model is used as the final model for the
give class. Instead, one could further refine the model by performing an
additional regression combining all data from the identified species.

#### Validate & Classify

In practice, the validation errors can easily be checked to satisfy log-
normalcy a posteriori. If this is not satisfied, it may not be straightforward
to cluster based on the validation error. In particular, chaotic trajectories
cannot be expected to achieve a low validation error, in which case another
metric is needed. In this case, it is reasonable to require that trajectories
to be long enough to compute statistics. We aim to investigate the
requirements for performing classification with chaotic interacting particles
in future work.

##  Appendix A Notation, Algorithm Hyperparameters, and Walltimes

Symbol | Definition
---|---
$f_{\text{a-r}}$ | attractive-repulsive force
$f_{\text{align}}$ | alignment force
$f_{\text{drag}}$ | drag force
$(x_{i},v_{i})$ | position and velocity of cell $i$
$\theta_{ij}$ | angle between $v_{i}$ and $x_{i}-x_{j}$
$N$ | number of focal cells selected for learning
$N_{\text{tot}}$ | total number of cells in the population
$\eta_{i}$ | noise in force law for cell $i$
$r_{nf}$ | near-field radius, below which $f_{\text{a-r}}$ is repulsive
$r_{ff}$ | far-field radius, above which $f_{\text{a-r}}$ is attractive
$F_{i}$ | total force on particle $i$
$(X,V)$ | entire particle system position and velocity
$(\mathbf{X},\mathbf{V})$ | particle position and velocity time-series data
$\Theta(X,V)=(f_{j}(X,V))_{1\leq j\leq J}$ | basis of force functions for learning
$(\Theta(\mathbf{X},\mathbf{V}),\ddot{\mathbf{X}})$ | linear system for learning using SINDy
$(\mathbf{G},\mathbf{b})$ | linear system for learning using WSINDy
$\mathbf{w}^{\star}$ | true model coefficients
$\\{\phi_{q}\\}_{q\in Q}$ | test functions to compute weak time derivatives
$(m,p,t_{q})$ | test function hyperparameters
${\mathcal{F}}_{\text{a-r}}$, ${\mathcal{F}}_{\text{align}}$, and ${\mathcal{F}}_{\text{drag}}$ | trial bases
$p_{\ell}$ | Laguerre polynomials
$r_{\max}$ | maximum inter-particle distance observed in $\mathbf{X}$
${\mathcal{M}}$ | set of single-cell models
$f_{\text{force}}^{(i)}$ | directional force modes (force $\in\\{\text{a-r},\text{align},\text{drag}\\}$)
${\mathcal{C}}$ | set of model clusters based on force modes
$\overline{{\mathcal{C}}}$ | cluster with the most members
$\overline{\mathbf{w}}$ | coefficients obtained from averaging the models in cluster $\overline{{\mathcal{C}}}$
$\overline{{\mathcal{M}}}$ | model associated with coefficients $\overline{\mathbf{w}}$
$\overline{{\mathcal{S}}}$ | species identified as obeying the model $\overline{{\mathcal{M}}}$
$S$ | Number of species identified
$(\overline{x}_{i},\overline{v}_{i})$ | simulated cells using $\overline{{\mathcal{M}}}$
$L$ | number of time steps in data
$\Delta t$ | time step of data
$\Delta t^{\prime}$ | time step for simulation
$\overline{{\mathcal{M}}}(x,v,X,V)$ | force on $(x,v)$ resulting from neighbors $(X,V)$ using model $\overline{{\mathcal{M}}}$
$(\mathbf{X}^{\prime i},\mathbf{V}^{\prime i})$ | set of particle positions and velocities with $i$th cell removed
$\Delta V_{i}$ | validation error of cell $i$
$VE$ | set of validation errors
Table 5: Summary of notation used througout.

In Table 5 we include all notation used throughout the article. In Table 6 we
list several hyperparameter choices for the examples above which were not
covered in Sections 3-4. Finally, we include walltimes for the main components
of the algorithm in Table 7, recorded in MATLAB using an AMD Ryzen 7 pro 4750u
processor.

Experiment | $m$ | $p$ | $r_{nf}$ | $r_{ff}$
---|---|---|---|---
$\mathbf{X}_{A}$ | $38$ | $8$ | $0.0497$ | 1
$\mathbf{X}_{B}$ | $38$ | $8$ | $0.0365$ | 1
$\mathbf{X}_{C}$ | $32$ | $9$ | $0.0219$ | 1
$\mathbf{X}_{A,C}$ | $35$ | $9$ | $0.0494$ | 1
$\mathbf{X}_{B,C}$ | $35$ | $9$ | $0.0499$ | 1
$\mathbf{X}_{A,B}$ | $38$ | $8$ | $0.0365$ | 1
$\mathbf{X}_{A,B}$(long) | $31$ | $9$ | $0.0219$ | 1
$\mathbf{X}_{A,B,C}$ | $38$ | $8$ | $0.0569$ | 1
$\mathbf{X}_{A,B,C}$(long) | $31$ | $9$ | $0.0253$ | 1
Table 6: Hyperparameters used for each example. $\mathbf{x}_{i}\to{\mathcal{M}}_{i}$ | ${\mathcal{M}}\to(\overline{{\mathcal{C}}},\overline{{\mathcal{M}}})$ | $(\mathbf{x}_{i},\overline{{\mathcal{M}}})\to\Delta V_{i}$ | $VE\to\overline{{\mathcal{S}}}$
---|---|---|---
5-10 sec. | 10-30 sec. | 10-30 sec. |  $<1$ sec.
Table 7: Wall times for main components of one iteration of the algorithm,
recorded for the $\mathbf{X}_{A,B,C}$(long) experiment ($N_{\text{tot}}=1000$
cells and $L=400$ timepoints). Each component contains an inner iteration
which may be trivially parallelized. The reported times are for one step of
the respective inner iteration. Specifically, it takes 5-10 seconds to learn
each model ${\mathcal{M}}_{i}$, while computation time for the full set of
models ${\mathcal{M}}=\\{{\mathcal{M}}_{1},\dots,{\mathcal{M}}_{N}\\}$ depends
on the number of available CPUs. To form the model clusters
${\mathcal{C}}=\\{{\mathcal{C}}_{1},\dots\mathcal{C}_{r}\\}$, find the largest
cluster $\overline{{\mathcal{C}}}$, and compute the averaged model
$\overline{{\mathcal{M}}}$, forward simulations of each model
${\mathcal{M}}_{i}$ (see Section D) are performed which take 10-30 seconds
(depending on the complexity of ${\mathcal{M}}_{i}$). Similarly, computation
of each $\Delta V_{i}$ requires one forward simulation (10-30 seconds).
Lastly, to identify $\overline{{\mathcal{S}}}$, it takes $<1$ second to
perform GMM classification of $\log_{10}(VE)$, and the results of 20 rounds of
classification are averaged. If full parallelization is available, the
walltime is less than two minutes per outer iteration, or $<2S$ minutes in
total, where $S$ is the number of identified species. For comparison, the cost
of generating the data for $\mathbf{X}_{A,B,C}$(long) takes 5-6 hours.

##  Appendix B Simulation details

All example data was generated from the exact models using Forward Euler with
a timestep of $\Delta t^{*}\approx 0.00042$ up until a final time of $T\approx
26$. The time series was then coarsened to a resolution of $\Delta t=0.13$,
resulting in a total of $L=200$ time points for learning. Experiments
$\mathbf{X}_{A,B}$(long) and $\mathbf{X}_{A,B,C}$(long) are simply extensions
of $\mathbf{X}_{A,B}$ or $\mathbf{X}_{A,B,C}$ by an addition 200 timepoints at
the same resolution. Initial positions were generated using Latin hypercube
sampling in a box of side length $2$, while initial velocities were drawn from
a Gaussian distribution with mean 0 and covariance $0.0025\mathbf{I}_{2}$.

##  Appendix C Constrained sparse regression

The constrained sequential thresholding algorithm requires solving at each
thresholding iteration $\ell$ a linearly constrained quadratic program of the
form

(C.1)
$\mathbf{w}^{(\ell+1)}=\operatorname*{arg\,min}_{\begin{subarray}{c}\mathbf{w}\text{
s.t. }\mathbf{C}\mathbf{w}\leq\mathbf{d}\\\
\text{supp}\left(\mathbf{w}\right)\subset{\mathcal{I}}^{(\ell)}\end{subarray}}\left\|{\mathbf{G}\mathbf{w}-\mathbf{b}}\right\|_{2}^{2}$

where ${\mathcal{I}}^{(\ell)}$ is the set of coefficients of
$\mathbf{w}^{(\ell)}$ satisfying (3.7). The constraint system
$\mathbf{C}\mathbf{w}\leq\mathbf{d}$ has the following four components.

  1. 1.

$f_{\text{a-r}}\geq 0$ when $0\leq r<r_{nf}$: for the near-field repulsion of
$f_{\text{a-r}}$ we discretize the region $\\{(r,\theta)\ :\ 0\leq r<r_{nf},\
\theta\in[0,2\pi)\\}$ choosing 5 equally-spaced points in $r$ from $10^{-6}$
to $r_{nf}$ and 5 equally-spaced points in $\theta$ from $0$ to $\pi$.
Evaluating each of the basis function for $f_{\text{a-r}}$ at this grid
results in a constraint system
$\mathbf{C}_{\text{a-r},nf}\mathbf{w}_{\text{a-r}}\leq\textbf{0}$ of dimension
$25\times J_{\text{a-r}}$ where $J_{\text{a-r}}$ is the number of basis
functions used to approximate $f_{\text{a-r}}$ and $\mathbf{w}_{\text{a-r}}$
is the restriction of $\mathbf{w}$ to coefficients of $f_{\text{a-r}}$.

  2. 2.

$f_{\text{a-r}}\leq 0$ when $r\geq r_{ff}$: similarly for the far-field
region, we choose $10$ equally-spaced points in $r$ from $r_{ff}$ to
$r_{\max}$, where $r_{\max}$ is the maximum observed neighbor-neighbor
distance in the simulation, and $\theta$ over the same points as previous.
This results in a constraint system
$\mathbf{C}_{\text{a-r},ff}\mathbf{w}_{\text{a-r}}\leq\textbf{0}$ of
dimensions $50\times J_{\text{a-r}}$.

  3. 3.

$f_{\text{align}}\leq 0$: since the basis is positive, the constrain system
here is simply
$\mathbf{I}_{J_{\text{align}}}\mathbf{w}_{\text{align}}\leq\textbf{0}$, where
$\mathbf{I}_{n}$ indicates the identity on $\mathbb{R}^{n}$.

  4. 4.

$f_{\text{drag}}\leq 0$: similarly the basis is positive, so the constraint
system is
$\mathbf{I}_{J_{\text{drag}}}\mathbf{w}_{\text{drag}}\leq\textbf{0}$.

Altogether we get $\mathbf{d}=\textbf{0}$ and

$\mathbf{C}=\begin{bmatrix}\mathbf{C}_{\text{a-r},nf}&0&0\\\
\mathbf{C}_{\text{a-r},ff}&0&0\\\ 0&\mathbf{I}_{J_{\text{align}}}&0\\\
0&0&\mathbf{I}_{J_{\text{drag}}}\end{bmatrix}$

We use MATLAB’s quadprog with constraint tolerance $10^{-10}$ and maximum
iterations set to 1000. Note that since $\mathbf{d}=\textbf{0}$, we do not
lose feasibility during the thresholding step, which is possible in general.

##  Appendix D Model replacement

Once the initial batch of $N$ models ${\mathcal{M}}$ is learned, we simulate
each model ${\mathcal{M}}_{i}$ as outlined in 3.4 on $K$ different validation
cells selected from the data, where $K=20$ throughout. For a given model
${\mathcal{M}}_{i}$, we select these $K$ validation cells by finding cells in
the population that match well certain statistics of cell $i$. In particular,
we define the following distributions:

$\displaystyle\rho^{(i)}_{rr}(r)$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\mathbb{P}_{x\sim
X^{\prime}}\left(\|x_{i}(t)-x(t)\|<r\right)dt$
$\displaystyle\rho^{(i)}_{vv}(s)$
$\displaystyle=\frac{1}{T}\int_{0}^{T}\mathbb{P}_{v\sim
V^{\prime}}\left(\|v_{i}(t)-v(t)\|<s\right)dt$
$\displaystyle\rho^{(i)}_{v}(s)$
$\displaystyle=\mathbb{P}\left(\|v_{i}\|<s\right)$

where $(X^{\prime},V^{\prime})$ denote the remainder of the cell population
excluding cell $i$. Respectively, these denote the distribution of distances
from cell $i$ to all other cells, the distribution of velocity differences
between cell $i$ and all other cells, and the distribution of speeds that cell
$i$ experiences131313These statistics are likely to correspond to the
information content that cell $i$ carries about its own forces
$f_{\text{a-r}}$, $f_{\text{align}}$, and $f_{\text{drag}}$, given the force
dependencies. We approximate these distributions from the data using
histograms with 50 bins.

For each cell $i$ we compute the Kullback-Leibler (KL) divergence between its
distributions $\rho^{(i)}_{rr}$, $\rho^{(i)}_{vv}$, $\rho^{(i)}_{v}$ and those
of the rest of the population141414To limit the computational overhead, we
only compute these KL divergences to cell $i$’s nearest 200 neighbors in the
Euclidean sense., where the KL divergence between densities $\rho$ and $\nu$
is given by

${\mathcal{D}}_{KL}(\rho|\nu)=-\int\rho(x)\log\left(\frac{\nu(x)}{\rho(x)}\right)\,dx.$

The $K$ validation cells used to validate model $i$ are the $K$ cells with
smallest cost $\mathscr{L}$, defined by

$\mathscr{L}:={\mathcal{D}}_{KL}(\rho^{(i)}_{rr}|\rho^{(j)}_{rr})^{2}+{\mathcal{D}}_{KL}(\rho^{(i)}_{vv}|\rho^{(j)}_{vv})^{2}+{\mathcal{D}}_{KL}(\rho^{(i)}_{v}|\rho^{(j)}_{v})^{2}.$

Let the validation error $\Delta V_{i\to j}$ be defined as in (3.11), but
indicating ${\mathcal{M}}_{i}$ used to validate cell $j$ (i.e. using the
initial conditions of cell $j$). We replace ${\mathcal{M}}_{i}$ with
${\mathcal{M}}_{j}$ if the following three conditions are met:

  1. (1)

$\Delta V_{i\to i}>\Delta V_{j\to i}$

  2. (2)

$\Delta V_{i\to j}>\Delta V_{j\to j}$

  3. (3)

$\max\\{\Delta V_{j\to i},\Delta V_{j\to j}\\}<\text{tol}$

where we set $\text{tol}=0.25$ in this work. In words, ${\mathcal{M}}_{j}$
replaces ${\mathcal{M}}_{i}$ if (1) ${\mathcal{M}}_{j}$ performs better than
${\mathcal{M}}_{i}$ on cell $i$, (2) ${\mathcal{M}}_{j}$ performs better than
${\mathcal{M}}_{i}$ on cell $j$, and (3) ${\mathcal{M}}_{j}$ achieves a
reasonably low error (defined by tol) on both cell $i$ and cell $j$. (Note
that cell $i$ and cell $j$ are required to be mutual validation cells for a
model replacement to occur). Furthermore, if ${\mathcal{M}}_{j}$ replaces
${\mathcal{M}}_{i}$, and another model ${\mathcal{M}}_{k}$ replaces
${\mathcal{M}}_{j}$, we reassign ${\mathcal{M}}_{i}$ to ${\mathcal{M}}_{k}$ as
well, even if cells $i$ and $k$ are not mutual validation cells.

We find this approach to be crucial to increasing accuracy of the learned
models, as it transfers successful learning of few cells with highly
informative trajectories to cells with less informative trajectories. As with
all validation steps of our algorithm, this approach would be infeasible if
not for fast data-driven forward simulations.

##  Appendix E Gaussian Mixture Model (GMM) classification

Since the GMM fitting is performed using the expectation-maximization
algorithm with random conditions, we perform the GMM fitting for 20 trials and
identify $\overline{{\mathcal{S}}}$ as the cells that in more than half of the
trials appear in the mixture with lowest error.

At some stage in the algorithm, all the remaining cells will be homogeneous.
In this case, a 2-species Gaussian is the wrong model. To account for this, we
do an initial fit to a single Gaussian and compute its Bayesian Information
Criterion (BIC). We accept the 2-mixture GMM if the average BIC of all 20
trial GMM fits with 2 mixtures is lower than that of the single Gaussian.

#### Data Access

All software used to generate the results in this work is available at this
repository: https://github.com/MathBioCU/WSINDy_CellCluster.git

#### Funding

This research was supported in part by the NSF/NIH Joint DMS/NIGMS
Mathematical Biology Initiative grant R01GM126559, in part by the NSF
Mathematical Biology MODULUS grant 2054085, and in part by the NSF Computing
and Communications Foundations grant 1815983. This work also utilized
resources from the University of Colorado Boulder Research Computing Group,
which is supported by the National Science Foundation (awards ACI-1532235 and
ACI-1532236), the University of Colorado Boulder, and Colorado State
University.

#### Acknowledgements

The authors wish to thank Prof. Vanja Dukić (Department of Applied
Mathematics, University of Colorado, Boulder) for insightful comments about
statistical aspects of this work as well as the mathematical form of the
directional interaction kernel.

## References

  * [1] H. Akaike.  A new look at the statistical model identification.  IEEE Transactions on Automatic Control, 19(6):716–723, December 1974.
  * [2] M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic.  Interaction ruling animal collective behavior depends on topological rather than metric distance: Evidence from a field study.  Proceedings of the National Academy of Sciences, 105(4):1232–1237, 2008.
  * [3] Sara Bernardi, Raluca Eftimie, and Kevin J. Painter.  Leadership Through Influence: What Mechanisms Allow Leaders to Steer a Swarm?  Bull Math Biol, 83(6):69, June 2021.
  * [4] Harish S. Bhat.  Learning and Interpreting Potentials for Classical Hamiltonian Systems.  In Peggy Cellier and Kurt Driessens, editors, Machine Learning and Knowledge Discovery in Databases, volume 1167, pages 217–228. Springer International Publishing, Cham, 2020.
  * [5] Manfred Boehm and Elizabeth G Nabel.  Cell cycle and cell migration: new pieces to the puzzle, 2001.
  * [6] C Bonneton, J-B Sibarita, and J-P Thiery.  Relationship between cell migration and cell cycle during the initiation of epithelial to fibroblastoid transition.  Cell motility and the cytoskeleton, 43(4):288–295, 1999.
  * [7] Alexander P Browning, Wang Jin, Michael J Plank, and Matthew J Simpson.  Identifying density-dependent interactions in collective cell behaviour.  Journal of The Royal Society Interface, 17(165):20200143, 2020.
  * [8] David B Brückner, Nicolas Arlt, Alexandra Fink, Pierre Ronceray, Joachim O Rädler, and Chase P Broedersz.  Learning the dynamics of cell–cell interactions in confined cell migration.  Proceedings of the National Academy of Sciences, 118(7), 2021.
  * [9] David B Brückner, Pierre Ronceray, and Chase P Broedersz.  Inferring the dynamics of underdamped stochastic systems.  Physical review letters, 125(5):058103, 2020.
  * [10] 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, 2016.
  * [11] Anna Q Cai, Kerry A Landman, and Barry D Hughes.  Modelling directional guidance and motility regulation in cell migration.  Bulletin of mathematical biology, 68(1):25–52, 2006.
  * [12] José A Carrillo and Young-Pil Choi.  Mean-field limits: from particle descriptions to macroscopic equations.  Archive for Rational Mechanics and Analysis, 241(3):1529–1573, 2021\.
  * [13] José A Carrillo, Massimo Fornasier, Giuseppe Toscani, and Francesco Vecil.  Particle, kinetic, and hydrodynamic models of swarming.  In Mathematical modeling of collective behavior in socio-economic and life sciences, pages 297–336. Springer, 2010.
  * [14] 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, 2019.
  * [15] Duxin Chen, Bowen Xu, Tao Zhu, Tao Zhou, and Hai-Tao Zhang.  Anisotropic interaction rules in circular motions of pigeon flocks: an empirical study based on sparse bayesian learning.  Physical Review E, 96(2):022411, 2017.
  * [16] Felipe Cucker and Steve Smale.  Emergent behavior in flocks.  IEEE Transactions on automatic control, 52(5):852–862, 2007.
  * [17] Raluca Eftimie, Gerda de Vries, and Mark A Lewis.  Complex spatial group patterns result from different animal communication mechanisms.  Proceedings of the National Academy of Sciences, 104(17):6974–6979, 2007.
  * [18] Andries P Engelbrecht.  Heterogeneous particle swarm optimization.  In International Conference on Swarm Intelligence, pages 191–202. Springer, 2010.
  * [19] Anders Eriksson, Martin Nilsson Jacobi, Johan Nyström, and Kolbjørn Tunstrøm.  Determining interaction rules in animal swarms.  Behavioral Ecology, 21(5):1106–1111, 2010.
  * [20] R Escobedo, V Lecheval, V Papaspyros, F Bonnet, F Mondada, Clément Sire, and Guy Theraulaz.  A data-driven method for reconstructing and modelling social interactions in moving animal groups.  Philosophical Transactions of the Royal Society B, 375(1807):20190380, 2020.
  * [21] Joep HM Evers, Razvan C Fetecau, and Lenya Ryzhik.  Anisotropic interactions in a first-order aggregation model.  Nonlinearity, 28(8):2847, 2015.
  * [22] U. Fasel, J. N. Kutz, B. W. Brunton, and S. L. Brunton.  Ensemble-SINDy: Robust sparse model discovery in the low-data, high-noise limit, with active learning and control.  Proc. R. Soc. Math. Phys. Eng. Sci., 478(2260):20210904, April 2022\.
  * [23] Jinchao Feng, Yunxiang Ren, and Sui Tang.  Data-driven discovery of interacting particle systems using gaussian processes.  arXiv preprint arXiv:2106.02735, 2021.
  * [24] Razvan C Fetecau, Yanghong Huang, and Theodore Kolokolnikov.  Swarm dynamics and equilibria for a nonlocal aggregation model.  Nonlinearity, 24(10):2681, 2011.
  * [25] Anna Haeger, Katarina Wolf, Mirjam M. Zegers, and Peter Friedl.  Collective cell migration: Guidance principles and hierarchies.  Trends Cell Biol., 25(9):556–566, September 2015.
  * [26] William R Holmes, JinSeok Park, Andre Levchenko, and Leah Edelstein-Keshet.  A mathematical model coupling polarity signaling to cell adhesion explains diverse cell migration patterns.  PLoS computational biology, 13(5):e1005524, 2017.
  * [27] Yael Katz, Kolbjørn Tunstrøm, Christos C Ioannou, Cristián Huepe, and Iain D Couzin.  Inferring the structure and dynamics of interactions in schooling fish.  Proceedings of the National Academy of Sciences, 108(46):18720–18725, 2011.
  * [28] Daniela Kengyel, Heiko Hamann, Payam Zahadat, Gerald Radspieler, Franz Wotawa, and Thomas Schmickl.  Potential of heterogeneity in collective behaviors: A case study on heterogeneous swarms.  In International conference on principles and practice of multi-agent systems, pages 201–217. Springer, 2015.
  * [29] Ian Knowles and Robert J Renka.  Methods for numerical differentiation of noisy data.  Electronic Journal of Differential Equations, 21:235–246, 2014\.
  * [30] 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.  Proc. R. Soc. A., 476(2234):20190800, February 2020.
  * [31] Tian Lan, Meng Yu, Weisheng Chen, Jun Yin, Hsiang-Tsun Chang, Shan Tang, Ye Zhao, Spyros Svoronos, Samuel WK Wong, and Yiider Tseng.  Decomposition of cell activities revealing the role of the cell cycle in driving biofunctional heterogeneity.  Scientific reports, 11(1):1–15, 2021.
  * [32] Herbert Levine, Wouter-Jan Rappel, and Inon Cohen.  Self-organization in systems of self-propelled particles.  Physical Review E, 63(1):017101, 2000.
  * [33] Zichao Long, Yiping Lu, Xianzhong Ma, and Bin Dong.  Pde-net: Learning pdes from data.  In International Conference on Machine Learning, pages 3208–3216. PMLR, 2018.
  * [34] Fei Lu, Mauro Maggioni, and Sui Tang.  Learning interaction kernels in heterogeneous systems of agents from multiple trajectories.  J. Mach. Learn. Res., 22:32–1, 2021.
  * [35] Ryan Lukeman, Yue-Xian Li, and Leah Edelstein-Keshet.  Inferring individual rules from collective behavior.  Proceedings of the National Academy of Sciences, 107(28):12576–12580, 2010.
  * [36] Daniel A Messenger and David M Bortz.  Learning mean-field equations from particle data using wsindy.  arXiv preprint arXiv:2110.07756, 2021.
  * [37] Daniel A Messenger and David M Bortz.  Weak SINDy for partial differential equations.  Journal of Computational Physics, page 110525, 2021.
  * [38] Daniel A Messenger and David M Bortz.  Weak SINDy: Galerkin-based data-driven model selection.  Multiscale Modeling & Simulation, 19(3):1474–1497, 2021.
  * [39] Daniel A. Messenger, Emiliano Dall’Anese, and David M. Bortz.  Online Weak-form Sparse Identification of Partial Differential Equations.  arXiv:2203.03979, (submitted), March 2022.
  * [40] Alexander Mogilner and Leah Edelstein-Keshet.  A non-local model for a swarm.  Journal of mathematical biology, 38(6):534–570, 1999.
  * [41] Rajnesh K Mudaliar and Timothy M Schaerf.  Examination of an averaging method for estimating repulsion and attraction interactions in moving groups.  Plos one, 15(12):e0243631, 2020.
  * [42] 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):20200987, 2021.
  * [43] John T Nardini, Douglas A Chapnick, Xuedong Liu, and David M Bortz.  Modeling keratinocyte wound healing dynamics: Cell–cell adhesion promotes sustained collective migration.  Journal of theoretical biology, 400:103–117, 2016.
  * [44] Jacob Notbohm, Shiladitya Banerjee, Kazage JC Utuje, Bomi Gweon, Hwanseok Jang, Yongdoo Park, Jennifer Shin, James P Butler, Jeffrey J Fredberg, and M Cristina Marchetti.  Cellular contraction and polarization drive collective cellular motion.  Biophysical journal, 110(12):2729–2738, 2016.
  * [45] Aditya A Paranjape, Soon-Jo Chung, Kyunam Kim, and David Hyunchul Shim.  Robotic herding of a flock of birds using an unmanned aerial vehicle.  IEEE Transactions on Robotics, 34(4):901–915, 2018.
  * [46] KR Pilkiewicz, BH Lemasson, MA Rowland, A Hein, J Sun, A Berdahl, ML Mayo, J Moehlis, M Porfiri, E Fernández-Juricic, et al.  Decoding collective communications using information theory tools.  Journal of the Royal Society Interface, 17(164):20190563, 2020.
  * [47] Sriram Ramaswamy.  The mechanics and statistics of active matter.  Annu. Rev. Condens. Matter Phys., 1(1):323–345, 2010.
  * [48] Craig W Reynolds.  Flocks, herds and schools: A distributed behavioral model.  In Proceedings of the 14th annual conference on Computer graphics and interactive techniques, pages 25–34, 1987.
  * [49] 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, 2017.
  * [50] Debarati Sarkar, Gerhard Gompper, and Jens Elgeti.  A minimal model for structure, dynamics, and tension of monolayered cell colonies.  Communications Physics, 4(1):1–8, 2021.
  * [51] Hayden Schaeffer.  Learning partial differential equations via data discovery and sparse optimization.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 473(2197):20160446, 2017.
  * [52] T. M. Schaerf, J. E. Herbert-Read, and A. J. W. Ward.  A statistical method for identifying different rules of interaction between individuals in moving animal groups.  J. R. Soc. Interface, 18(176):rsif.2020.0925, 20200925, March 2021\.
  * [53] Michael Schmidt and Hod Lipson.  Distilling free-form natural laws from experimental data.  science, 324(5923):81–85, 2009.
  * [54] Linus J. Schumacher, Philip K. Maini, and Ruth E. Baker.  Semblance of Heterogeneity in Collective Cell Migration.  Cell Systems, 5(2):119–127.e1, August 2017.
  * [55] Thomas D Seeley, P Kirk Visscher, and Kevin M Passino.  Group decision making in honey bee swarms: When 10,000 bees go house hunting, how do they cooperatively choose their new nesting site?  American scientist, 94(3):220–229, 2006.
  * [56] David J. T. Sumpter, Alex Szorkovszky, Alexander Kotrschal, Niclas Kolm, and James E. Herbert-Read.  Using activity and sociability to characterize collective motion.  Phil. Trans. R. Soc. B, 373(1746):20170015, May 2018.
  * [57] Rohit Supekar, Boya Song, Alasdair Hastewell, Alexander Mietke, and Jörn Dunkel.  Learning hydrodynamic equations for active matter from particle simulations and experiments.  arXiv preprint arXiv:2101.06568, 2021.
  * [58] Dhananjay T Tambe, C Corey Hardin, Thomas E Angelini, Kavitha Rajendran, Chan Young Park, Xavier Serra-Picamal, Enhua H Zhou, Muhammad H Zaman, James P Butler, David A Weitz, et al.  Collective cell guidance by cooperative intercellular forces.  Nature materials, 10(6):469–475, 2011.
  * [59] Dean G Tang.  Understanding cancer stem cell heterogeneity and plasticity.  Cell research, 22(3):457–472, 2012.
  * [60] John Toner and Yuhai Tu.  Flocks, herds, and schools: A quantitative theory of flocking.  Physical review E, 58(4):4828, 1998.
  * [61] Chad M Topaz and Andrea L Bertozzi.  Swarming patterns in a two-dimensional kinematic model for biological groups.  SIAM Journal on Applied Mathematics, 65(1):152–174, 2004.
  * [62] Kolbjørn Tunstrøm, Yael Katz, Christos C. Ioannou, Cristián Huepe, Matthew J. Lutz, and Iain D. Couzin.  Collective States, Multistability and Transitional Behavior in Schooling Fish.  PLoS Comput Biol, 9(2):e1002915, February 2013.
  * [63] Silviu-Marian Udrescu and Max Tegmark.  Ai feynman: A physics-inspired method for symbolic regression.  Science Advances, 6(16):eaay2631, 2020.
  * [64] Floris Van Breugel, J Nathan Kutz, and Bingni W Brunton.  Numerical differentiation of noisy data: A unifying multi-objective optimization framework.  IEEE Access, 8:196865–196877, 2020.
  * [65] Tamás Vicsek and Anna Zafeiris.  Collective motion.  Physics reports, 517(3-4):71–140, 2012.
  * [66] Medhavi Vishwakarma, Jacopo Di Russo, Dimitri Probst, Ulrich S Schwarz, Tamal Das, and Joachim P Spatz.  Mechanical interactions among followers determine the emergence of leaders in migrating epithelial cell collectives.  Nature communications, 9(1):1–12, 2018.
  * [67] Ashley J. W. Ward, Timothy M. Schaerf, James E. Herbert-Read, Lesley Morrell, David J. T. Sumpter, and Mike M. Webster.  Local interactions and global properties of wild, free-ranging stickleback shoals.  R. Soc. open sci., 4(7):170043, July 2017.
  * [68] William H. Warren.  Collective Motion in Human Crowds.  Curr Dir Psychol Sci, 27(4):232–240, August 2018.
  * [69] Ming Zhong, Jason Miller, and Mauro Maggioni.  Data-driven discovery of emergent behaviors in collective dynamics.  Physica D: Nonlinear Phenomena, 411:132542, 2020.
  * [70] Cole Zmurchok and Gerda de Vries.  Direction-dependent interaction rules enrich pattern formation in an individual-based model of collective behavior.  PloS one, 13(6):e0198550, 2018.