WGMwgm QEqe

1]organization=Department of Mathematics,addressline=University of California,
Riverside, state=CA, country=US

2]organization=Department of Mathematics,addressline=Lafayette College,
city=Easton, state=PA, country=US

3]organization=Department of Mathematics,addressline=University of Oxford,
city=Oxford, country=UK

4]organization=Interdisciplinary Center for Quantitative Modeling in Biology,
University of California Riverside, country=USA

[3]Corresponding author

# Designing experimental conditions to use the Lotka-Volterra model to infer
tumor cell line interaction types

Heyrim Cho  Allison L. Lewis  Kathleen M. Storey  Helen M. Byrne  [  [  [  [

###### Abstract

The Lotka-Volterra model is widely used to model interactions between two
species. Here, we generate synthetic data mimicking competitive, mutualistic
and antagonistic interactions between two tumor cell lines, and then use the
Lotka-Volterra model to infer the interaction type. Structural identifiability
of the Lotka-Volterra model is confirmed, and practical identifiability is
assessed for three experimental designs: (a) use of a single data set, with a
mixture of both cell lines observed over time, (b) a sequential design where
growth rates and carrying capacities are estimated using data from experiments
in which each cell line is grown in isolation, and then interaction parameters
are estimated from an experiment involving a mixture of both cell lines, and
(c) a parallel experimental design where all model parameters are fitted to
data from two mixtures (containing both cell lines but with different initial
ratios) simultaneously. Each design is tested on data generated from the
Lotka-Volterra model with noise added, to determine efficacy in an ideal
sense. In addition to assessing each design for practical identifiability, we
investigate how the predictive power of the model—i.e., its ability to fit
data for initial ratios other than those to which it was calibrated—is
affected by the choice of experimental design. The parallel calibration
procedure is found to be optimal and is further tested on in silico data
generated from a spatially-resolved cellular automaton model, which accounts
for oxygen consumption and allows for variation in the intensity level of the
interaction between the two cell lines. We use this study to highlight the
care that must be taken when interpreting parameter estimates for the
spatially-averaged Lotka-Volterra model when it is calibrated against data
produced by the spatially-resolved cellular automaton model, since baseline
competition for space and resources in the CA model may contribute to a
discrepancy between the type of interaction used to generate the CA data and
the type of interaction inferred by the LV model.

###### keywords:

Lotka-Volterra model parameter identifiability experimental design high- to
low-fidelity model calibration

##  1 Introduction

In many cancers, cellular heterogeneity plays a significant role in resistance
to treatment and tumor recurrence [32, 43]. Thus, it is crucial to gain a
better understanding of how a tumor’s growth dynamics and its response to
treatment influence—and are influenced by—cellular heterogeneity and the
interactions between different cell populations, in order to make more
accurate predictions about treatment responses. Mathematical modeling is an
important tool that can be used to gain insight into the mechanisms that drive
tumor heterogeneity and to make patient-specific predictions. Ordinary
differential equation models are especially useful when fitting to sparse
data, and the efficacy of such models has been compared in multiple studies
[5, 1, 13, 26, 4, 17, 40]. However, to ensure practical utility of such models
for making predictions that are accurate and contain limited uncertainty,
unique parameter identifiability—in the sense that the same model output
cannot be produced by two different sets of parameter values—must be achieved.
A thorough discussion of the importance of parameter identifiability in
systems biology models is provided in [48]. Further, Simpson et al. (2022)
performed an identifiability analysis on a series of population growth models,
of the type often used to simulate tumor growth [40]. Here, we investigate the
parameter identifiability of the Lotka-Volterra model, a two-compartment
ordinary differential equation model, in order to assess its ability to
distinguish different types of interactions between two tumor cell lines.

The Lotka-Volterra (LV) model is widely used to model interactions between two
populations [28, 47]. Its structure is sufficiently flexible to describe
competitive, mutualistic, and antagonistic interactions. In particular, the
competitive LV model has been widely used to study competition between
different subpopulations of cancer cells, the host tissue, and other cell
populations [16, 11, 15, 14, 22]. However, recent studies suggest that
interactions between different cancer populations may be more complex than
competition, and may contribute to tumor invasion and treatment response [34,
7, 49, 6, 33, 42]. Since tumor cell lines can interact in multiple ways, with
consequential implications for tumor progression and recurrence, we seek to
determine whether it is possible to infer interaction type when fitting the LV
model to data from two interacting tumor cell lines. In Paczkowski et al.
(2021), a preliminary investigation was performed by fitting the LV model to
dynamic time course data from tumor spheroids cultured from prostate cancer
cell lines [34]. Previously, a range of techniques has been used to estimate
parameter values in the LV model from experimental data [39, 24, 27, 31]. We
extend this work by first considering the structural identifiability of the LV
model given perfect, noise-free data about the volumes of two cell lines. We
then assess its practical identifiability using synthetic data generated from
the LV model with noise added, and from a spatially-resolved cellular
automaton (CA) model. Our CA model incorporates various types of cellular
interaction and adds to the growing body of literature in which the growth of
tumor spheroids is simulated using different types of agent-based models
(e.g., on-lattice models in [29, 30, 44, 36] and off-lattice models in [21,
45, 23, 18]).

As part of the practical identifiability investigation, we compare the ability
of different experimental design structures to fit the LV model and to infer
the type of interaction between the two cell lines given dynamic tumor data on
the volumes (and/or proportions) of the two cell types. We consider three
experimental design structures: (a) individual calibration, where we fit to a
single data set which contains dynamic volume information for both cell lines;
(b) sequential calibration of the cell-line specific intrinsic growth rates
and carrying capacities using dynamic data from each cell line, and subsequent
calibration of the interaction parameters using time course data from an
experiment in which the two cell lines are co-cultured; and (c) parallel
calibration of all model parameters using data from two mixture experiments,
which differ in terms of their initial conditions.

We use synthetic data from the LV model to identify the optimal experimental
design scheme, to eliminate issues relating to model discrepancy. We then
calibrate the model to synthetic data generated from the CA model. The CA
model accounts for spatial heterogeneity, oxygen consumption, and allows the
intensity of the imposed cell-cell interactions to be varied. This enables us
to study the robustness of the calibration procedure on a data set that better
approximates in vitro experimental data for multicellular tumor spheroids.

The results from our structural and practical identifiability analyses provide
insight into the type of experimental design and data needed to obtain
identifiable parameter estimates and infer the interaction type between two
cell lines. Our investigations suggest that where volumes data about both cell
lines is available for multiple sets of initial conditions, the LV model can
be fitted to the data and used to make predictions that may inform treatment
decisions in the clinic. At the same time, care is needed when interpreting
parameter estimates of the (spatially-averaged) LV model when it is fit to
spatially-resolved data.

In Section 2, we describe the LV model and the methods that we use to assess
its structural and practical identifiability. In Section 3 we demonstrate the
structural identifiability of the LV model, and assess its practical
identifiability under the three proposed calibration procedures using
synthetic data generated by the LV model with noise added. We then test the
parallel calibration procedure on data generated from the CA model. We also
discuss the limitations of using a spatially-averaged model, such as the LV
model, to make inferences about spatially-resolved data. In Section 4, we
summarize our findings and describe possible directions for future work.
Further details on the CA model are described in Appendix A, and the
sequential calibration results of fitting the LV model to the CA data are
presented in Appendix B.

##  2 Methods

###  2.1 Lotka-Volterra Model and Synthetic Data

We use the Lotka-Volterra (LV) model [28, 47] to describe the growth of a
heterogeneous tumor spheroid comprised of two cancer cell types, which we term
Type-$S$ and Type-$R$. In isolation, each population is assumed to undergo
logistic growth, with growth rates $r_{S}$ and $r_{R}$ and carrying capacities
$K_{S}$ and $K_{R}$ respectively. Interactions between the two populations are
determined by the signs and magnitudes of the parameters $\gamma_{S}$ and
$\gamma_{R}$. Written in terms of these parameters, the governing equations
read:

$\displaystyle\frac{dS}{dt}$ $\displaystyle=$ $\displaystyle
r_{S}S\left(1-\frac{S}{K_{S}}-\frac{\gamma_{R}R}{K_{S}}\right),$ (1)
$\displaystyle\frac{dR}{dt}$ $\displaystyle=$ $\displaystyle
r_{R}R\left(1-\frac{R}{K_{R}}-\frac{\gamma_{S}S}{K_{R}}\right),$ (2)

where $S(t)$ and $R(t)$ represent the volume (mm3) of Type-$S$ and Type-$R$
cancer cells at time $t$, respectively. We close Equations (1)-(2) by assuming
that the initial volumes of the Type-$S$ and Type-$R$ cells are known, so
that:

$S(0)=S_{0}\quad\mbox{and}\quad R(0)=R_{0}.$ (3)

The parameter set consists of
$\theta=\\{r_{S},r_{R},K_{S},K_{R},\gamma_{S},\gamma_{R}\\}$, where $r_{S}$,
$r_{R}$, $K_{S}$, and $K_{R}$ are assumed to be non-negative, while
$\gamma_{S}$ and $\gamma_{R}$ may be positive, negative, or zero. The signs of
$\gamma_{R}$ and $\gamma_{S}$ define the type of interaction between the two
cell lines: it is competitive if $\gamma_{R}$ and $\gamma_{S}$ are positive,
mutualistic if they are both negative, and antagonistic if they are of
opposite signs (see Table 1 for a summary).

| $\gamma_{S}>0$ | $\gamma_{S}=0$ | $\gamma_{S}<0$
---|---|---|---
$\gamma_{R}>0$ | Competitive |  | R antagonizes S
$\gamma_{R}=0$ |  | Neutral |
$\gamma_{R}<0$ | S antagonizes R |  | Mutual
Table 1: A summary of how the signs of the parameters $\gamma_{S}$ and
$\gamma_{R}$ in Equations (1)-(2) determine the type of interaction between
the Type-$S$ and Type-$R$ cells.

In Section 3.2, we assess the practical identifiability of the LV model
parameters by fitting it to data from the LV model. We generate 29 data sets:
two pure experiments (i.e., homogeneous tumor data comprising only Type-$S$ or
Type-$R$ cells), and nine mixture experiments with initial ratios 1:9, 2:8,…,
9:1 of Type-$S$:Type-$R$ cells, for competitive, mutual, and antagonistic
interaction types (assuming, without loss of generality, that Type-$R$ cells
antagonize Type-$S$ cells). Following [34], we fix $r_{S}=r_{R}=0.36$ and
$K_{S}=K_{R}=0.85$, and vary $\gamma_{S}$ and $\gamma_{R}$ to simulate the
different interaction types; Table 2 summarizes the parameter values used to
generate the data sets.

Interaction Type | $r_{S}$ | $r_{R}$ | $K_{S}$ | $K_{R}$ | $\gamma_{S}$ | $\gamma_{R}$
---|---|---|---|---|---|---
Competitive | 0.36 | 0.36 | 0.85 | 0.85 | 0.5 | 0.5
Mutual | 0.36 | 0.36 | 0.85 | 0.85 | -0.5 | -0.5
Antagonistic | 0.36 | 0.36 | 0.85 | 0.85 | -0.5 | 0.5
Table 2: Parameter values used to generate synthetic data from the Lotka-
Volterra model.

All simulated tumors are initialized with a total tumor volume of 0.02 mm3;
thus, the nine sets of initial conditions are given by
$(S_{0}^{j},R_{0}^{j})=(0.002j,0.02-0.002j)$ for $j=1,2...,9$. Volume
measurements for both cell lines are generated at $t=7,14,21,\dots,56$ days.
For a given initial ratio, and at each data point—denoted $(S_{i},R_{i})$ for
$i=1,\dots,8$—we add 5% noise: i.e.,
$S_{\text{noise}}(t_{i})=S(t_{i})(1+0.05\xi_{S,i})$ and
$R_{\text{noise}}(t_{i})=R(t_{i})(1+0.05\xi_{R,i})$, where $\xi_{S}$ and
$\xi_{R}$ are uniform random variables on $[-1,1]$. Sample synthetic data for
the 1:9, 5:5, and 9:1 initial ratios for the three interaction types are shown
in Figure 1. Since the intrinsic growth rates and carrying capacities for both
cell lines are identical, the simulations for the competitive and mutual cases
predict long-time coexistence of the two cell lines at equal volumes. In the
antagonistic cases, the dominance of the Type-$R$ cells results in a higher
predicted volume for the Type-$R$ population, regardless of the initial
conditions.

(a)

(b)

(c)

Figure 1: Synthetic data from the Lotka-Volterra model, generated with the
parameter values listed in Table 2 and initial Type-$S$:Type-$R$ ratios of
1:9, 5:5, and 9:1, for (a) competitive, (b) mutual, and (c) antagonistic
interaction types.

###  2.2 Cellular Automaton Model and Synthetic Data

As a test of the ability of the LV model to fit in vitro tumor spheroid data,
we extend a cellular automaton (CA) model from [34, 13, 12]. Cells are
arranged on a two-dimensional lattice, representing a two-dimensional cross-
section of a three-dimensional tumor spheroid in vitro. Each cell can be
classified as proliferating, quiescent, or necrotic. Cell state is determined
by the local oxygen concentration, which is modeled using a reaction-diffusion
equation. As in the LV model, we consider two types of tumor cells, Type-$S$
and Type-$R$, which consume oxygen, proliferate, and die at rates which depend
on the local $O_{2}$ concentration. As in the LV data generation discussed in
Section 2.1, we further assume that the two cell types are identical except
for how they interact with each other, although this could later be modified
to represent differences in other characteristics, e.g. radiosensitive versus
radioresistant.

In the CA model, we consider four interaction types: neutral, competitive,
mutualistic, and antagonistic (where, without loss of generality, the Type-$R$
population antagonizes the Type-$S$ population, and the Type-$S$ population
promotes the Type-$R$ population). The interaction type affects both the
oxygen consumption and proliferation rates. For example, in the antagonistic
case, Type-$S$ cells surrounded by a large number of Type-$R$ cells consume
less oxygen, and as a result, divide at a slower rate; likewise, Type-$R$
cells increase their rates of oxygen consumption and division when surrounded
by Type-$S$ cells.

The model includes a parameter, $I$, which determines the intensity of the
cellular interactions, with $I=1$ for neutral interactions (all cell cycle and
oxygen consumption rates fixed at their baseline values), and $I>1$ for
stronger interactions (the cell cycle and oxygen consumption rates altered to
reflect the interaction type). With the exception of the sensitivity analysis
presented in Figure 15, we fix $I=4$ for competitive, mutual, and antagonistic
interactions. Further details about the CA model are included in A. Cellular
Automaton Model.

We generate in silico data from the CA model for the four interaction types,
and different values of the initial proportion of Type-$S$ cells, as for the
LV data. In all cases, we initialize the cells in a circular region with a
radius of 9 cells, corresponding to a total tumor volume of approximately 0.02
mm3. The cell types at each initial site are randomly chosen, according to the
specified initial conditions. For each choice of initial conditions and
parameter values, we simulate 10 model replicates and record the average
volumes of Type-$S$ cells and Type-$R$ cells at $t=7,14,21,\dots,56$, to
reduce the impact of stochastic outliers on the calibration process.

Figure 2(a) shows how the spatial distributions of Type-$S$ and Type-$R$ cells
change over time for typical simulations, with an initial ratio of
Type-$S$:Type-$R$ cells of 9:1. This figure compares the distribution of cells
for the competitive, mutual, and antagonistic interaction types. Figure 2(b)
shows how the total numbers of Type-$S$ cells and Type-$R$ cells change over
time for the particular simulations shown in part (a) and the averaged
replicates. We observe that in the competitive case, the Type-$S$ cells
quickly become dominant, due to the large proportion of Type-$S$ cells at
$t=0$. In the mutual case, the Type-$R$ population catches up with the
Type-$S$ population size by day 10, and both cell types are able to co-exist
throughout the tumor. In the antagonistic case, the Type-$R$ cells have a
selective advantage, so they overtake the Type-$S$ cells by day 10 and proceed
to dominate the spheroid. In contrast to the data generated by the LV model in
Section 2.1, asymptotic coexistence of the two cell lines is predicted only
for the mutual case. The parameter values used to generate the simulated tumor
spheroids are given in Table 5.

Figure 2: CA simulation outputs, showing (a) two-dimensional cross sections of
spheroids at Days 3, 10, 20, and 40, for competitive, mutual, and antagonistic
interaction types, and (b) mean trajectories, with error bars indicating the
standard deviation, from 10 replicates of the CA. The yellow and purple curves
correspond to the specific simulations shown in part (a). In all cases, the
simulations were initialized with a 9:1 ratio of Type-$S$ to Type-$R$ cells,
with parameters as listed in Table 5.

###  2.3 Model Calibration

Throughout the investigation, we use Bayesian inference to perform parameter
estimation. In more detail, we use a Metropolis Hastings algorithm to generate
Markov parameter chains whose stationary distributions represent the posterior
distributions of the parameters. These posterior densities can be used to
estimate optimal values of the parameters when the LV model is fitted to the
data, denoted by
$\hat{\theta}=\\{\hat{r}_{S},\hat{r}_{R},\hat{K}_{S},\hat{K}_{R},\hat{\gamma}_{S},\hat{\gamma}_{R}\\}$,
and reflect the uncertainty in the estimates. We employ non-informative
uniform prior distributions of $\mathcal{U}(0,1)$ for
$\\{r_{S},r_{R},K_{S},K_{R}\\}$, and $\mathcal{U}(-3,3)$ for the interaction
parameters $\gamma_{S}$ and $\gamma_{R}$. The lower bounds for $r_{S}$,
$r_{R}$, $K_{S}$, and $K_{R}$ are chosen to ensure positivity of these
parameter estimates, while the ranges for $\gamma_{S}$ and $\gamma_{R}$ are
chosen to allow inference of all interaction types (mutual, competitive, and
antagonistic). The optimal parameter estimates are identified by maximizing
the likelihood function; if we assume that errors are independent and normally
distributed as $\varepsilon\sim\mathcal{N}(0,\sigma^{2})$, then the likelihood
function is of the form

$\mathcal{L}(y|\theta)=\frac{1}{(2\pi\sigma^{2})^{n/2}}\exp{\left(\frac{1}{2\sigma^{2}}\sum\limits_{i=1}^{n}[y_{i}-f_{i}(\theta)]^{2}\right)},$

where $y_{i}$ represents the $i$th observed data point for $i=1,\dots,n$, and
$f_{i}(\theta)$ represents the quantity of interest for parameter set $\theta$
at time $t_{i}$. With this formulation, maximizing the likelihood function is
equivalent to minimizing the sum-of-squares error between the observed data
and the model predictions. The Markov chain of accepted candidates is
constructed using the proposal function
$J(\theta^{*}|\theta^{k-1})=\mathcal{N}(\theta^{k-1},C)$, where $C$ represents
the estimated covariance matrix for the parameter set $\theta$. That is, each
parameter candidate is sampled from a multivariate normal distribution
centered at the previous accepted candidate. If the new candidate,
$\theta^{*}$, improves the likelihood function (i.e., yields a smaller sum-of-
squares error than the previous candidate), we set $\theta^{k}=\theta^{*}$,
accepting this candidate into our Markov chain. Else, we accept the candidate
with probability $1-\alpha$, where $\alpha$ is defined as the ratio of the
likelihoods using the current candidate as the numerator and the previous
candidate as the denominator—see [41] for details—and reject otherwise. The
occasional acceptance of a parameter candidate whose likelihood is inferior to
those already sampled promotes well-mixing of the Markov parameter chains and
full exploration of the admissible parameter space, permitting quicker
convergence to the stationary distribution and ensuring that the chain does
not stall in a local optimum. Upon convergence, these Markov chains can be
used to visualize the parameter posterior densities. Full details regarding
the implementation of the Metropolis Hastings algorithm are provided at [41,
20]. For all calibration procedures, we use the MCMC toolbox for MATLAB with
the Delayed Rejection Adaptive Metropolis variation of the Metropolis Hastings
algorithm [25].

###  2.4 Identifiability Analysis

#### Structural Identifiability

Before undertaking parameter estimation, it is important to establish whether
the model is a priori structurally identifiable, i.e., is not over-
parametrized. At the input level, we say that a parameter is not structurally
identifiable if perturbing the parameter can still result in the same model
output through compensation via changes in other parameters [48, 9, 41]; else,
it is structurally identifiable. A model is only structurally identifiable if
all of its parameters are structurally identifiable; essentially, we require a
1-to-1 mapping between the parameter set and model output. Several algorithms
for assessing structural identifiability have been proposed, including methods
that employ generating series, power series expansions, differential algebra,
and differential geometry. For further details, we refer the interested reader
to [46, 3, 10, 37].

Previous studies have used algebraic methods and Lie derivatives to
investigate the structural identifiability of the LV model [38, 19]. In this
investigation, structural identifiability of the LV model (1)-(2) will be
assessed using the Taylor series expansion method as outlined in [8, 35].

We denote by $y(t;\theta)$ the time-dependent observable quantity (or
quantities) evaluated at parameter set $\theta$. Assuming that $y(t;\theta)$
is analytic in a neighborhood of the initial conditions $t=t_{0}$, we can
express $y(t;\theta)$ and each of its time derivatives in terms of the model
parameters and initial conditions using the Taylor series expansion

$\displaystyle y(t;\theta)=y(t_{0};\theta)+y^{(1)}(t_{0};\theta)\
t+y^{(2)}(t_{0};\theta)\ \frac{t^{2}}{2!}+\dots+y^{(i)}(t_{0};\theta)\
\frac{t^{i}}{i!}+\dots,$

where $y^{(i)}(t_{0};\theta)$ denotes the $i$th time derivative of the
observable quantity evaluated at the initial condition $t_{0}$ with parameter
set $\theta$.

By considering the coefficients of the Taylor series expansion, we can restate
the problem of structural identifiability as a system of algebraic equations
relating our unknown parameter set $\theta$ to the observable quantity
$y(t;\theta)$ and its time derivatives at $t=t_{0}$, where the number of
equations to be considered depends upon the number of parameters to be
estimated. If these equations admit a unique solution for $\theta$ in terms of
$y(t;\theta)$ and its derivatives at $t=t_{0}$, we declare the model to be
globally structurally identifiable. If no such solution exists, we conclude
that the model is either non-identifiable, or locally structurally
identifiable, in the sense that the identifiability of the model parameters
depends on the region of parameter space in which it is evaluated.

Establishing structural identifiability of nonlinear systems is a challenging
task, whose complexity increases with the dimension of the parameter space.
The GenSSI (Generating Series for testing Structural Identifiability) MATLAB
toolbox enables non-experts to perform such analyses [9]. It employs the
generating series approach and produces identifiability tableaux, which can be
used to determine how to handle parameters classified as non-identifiable [2].
When assessing structural identifiability of the LV model in Section 3.1, we
will confirm the results of our Taylor series identifiability analysis using
the GenSSI toolbox.

#### Practical Identifiability

The focus of practical, or a posteriori, identifiability is to determine
whether model parameters can be inferred from potentially noisy data via model
calibration. Potential issues with practical identifiability include:
difficulty inferring parameter values due to measurement errors, model
discrepancy (when a model is unable to accurately represent the underlying
system), and an experimental design that is unable to naturally “excite”—i.e.,
activate—certain parameters. Each issue may preclude a posteriori parameter
identification, even for structurally identifiable models. These issues depend
on the quality, design, and availability of the data, rather than the model
structure.

For our purposes, we say that a parameter is practically identifiable if the
posterior distribution obtained from the Metropolis Hastings algorithm is
unimodal, exhibiting a clear and unique optimum. Otherwise, the parameter is
deemed practically non-identifiable, in the sense that multiple values of the
input may yield the same model output. In this way, non-identifiability of a
parameter may manifest as a posterior distribution that is
multimodal—suggesting multiple local maxima—or as a distribution that is
relatively unchanged from its non-informative prior distribution, indicating
that the parameter is uninformed by the data. A model is pronounced
practically identifiable given a particular data set only if all model
parameters are practically identifiable.

###  2.5 Experimental Design Options

For our practical identifiability investigation, we compare three methods to
calibrate the LV model to data. First, we consider whether the parameters in
the LV model may be uniquely identified using data from a single
experiment—containing time course volumetric data about both cell lines—with
the initial ratio of the two cell lines specified, as described in Section
2.1. We term this method “individual calibration”.

Our second method is termed “sequential calibration.” Here, we first estimate
$(r_{S},K_{S})$ and $(r_{R},K_{R})$, the intrinsic growth rates and carrying
capacities for each cell line, using data describing the dynamics of one cell
line grown in isolation. Subsequently, these parameters are fixed at their
estimated values $\\{\hat{r}_{S},\hat{r}_{R},\hat{K}_{S},\hat{K}_{R}\\}$, and
the interaction parameters $\gamma_{S}$ and $\gamma_{R}$ are inferred from a
third data set containing a mixture of the two cell lines. Thus, this method
requires three data sets—two relating to homogeneous growth of each cell line
and one to heterogeneous growth of a mixture.

Our final calibration procedure is referred to as “parallel calibration,” as
all six parameters are estimated simultaneously using two mixture data sets
which differ only in terms of their initial conditions. For example, we might
attempt to infer the complete parameter set
$\\{r_{S},r_{R},K_{S},K_{R},\gamma_{S},\gamma_{R}\\}$ by fitting the LV model
to growth data initialized with ratios of 1:9 and 9:1 for the Type-$S$ and
Type-$R$ cells.

A summary of the three calibration methods is presented in Figure 3. All
methods are assessed in Section 3.2, using synthetic data generated from the
LV model. In Section 3.3, we conduct further tests, applying the parallel
procedure to data generated from the CA model.

Figure 3: Summary of the three calibration procedures outlined in Section 2.5.

###  2.6 Assessment Metrics

We use three metrics—termed $E_{1}$, $E_{2}$, and $E_{3}$—to assess the
interaction inference, accuracy, and predictive power of the resulting model
fits. Metric $E_{1}$ is a label which indicates the interaction type inferred
by the LV model, based upon the signs of the interaction parameters:

$\displaystyle E_{1}=\left\\{\begin{array}[]{cc}C,&\text{ if }\
\hat{\gamma}_{S}>0,\ \hat{\gamma}_{R}>0\\\ M,&\text{ if }\
\hat{\gamma}_{S}<0,\ \hat{\gamma}_{R}<0\\\ A,&\text{ if }\
\hat{\gamma}_{S}<0,\ \hat{\gamma}_{R}>0.\end{array}\right.$ (7)

where $C$, $M$, and $A$ indicate “competitive”, “mutual”, and “antagonistic”
relationships, respectively. We note that the interaction type inferred by the
LV model may not match the type of interaction used to generate the data.

Metric $E_{2}$ measures the accuracy and predictive power of the inferred
model fit. That is, in addition to measuring the error of the model fit to the
training data set used for calibration, we quantify the accuracy with which
the model reproduces the dynamics of other heterogeneous tumors which differ
from the training data set in terms of their initial conditions, using testing
data sets generated from the other eight initial ratio combinations. We denote
the observed data by $d_{i}^{j}=(s_{i}^{j},r_{i}^{j}),$ where $s_{i}^{j}$ and
$r_{i}^{j}$ represent the volume measurements of Type-$S$ and Type-$R$ cells,
respectively, with initial ratio $j:(10-j)$, for $j=1,\dots,9$ at time $t_{i}$
for $i=1,\dots,N$. We denote the solution to the LV model at time $t_{i}$,
parameterized with parameters $\theta$ and initial conditions $[S(0),R(0)]$,
by $y(t_{i};\theta,[S(0),R(0)])=(S(t_{i}),R(t_{i}))$. Metric $E_{2}$ is then
defined to be the log of the sum-of-squares error between the training and
testing data sets and the predictions of the LV model for all nine initial
conditions. Thus, we have

$\displaystyle
E_{2}=\ln\left(\sum_{i=1}^{N}\sum_{j=1}^{9}||d_{i}^{j}-y(t_{i};\hat{\theta},[S_{0}^{j},R_{0}^{j}])||_{2}^{2}\right),$
(8)

where $\hat{\theta}$ represents the inferred parameter set resulting from the
calibration of the LV model to the training data set.

Our third metric measures the ability of the LV model and proposed calibration
procedure to recover the true values of the parameters used to generate the
data. As such, it is used only in Section 3.2 where we use the LV model to
generate the training and testing data. If the true parameter values are
denoted by
$\tilde{\theta}=\\{\tilde{r}_{S},\tilde{r}_{R},\tilde{K}_{S},\tilde{K}_{R},\tilde{\gamma}_{S},\tilde{\gamma}_{R}\\}$
and the estimated parameters by
$\hat{\theta}=[\hat{r}_{S},\hat{r}_{R},\hat{K}_{S},\hat{K}_{R},\hat{\gamma}_{S},\hat{\gamma}_{R}]$,
then metric $E_{3}$ is defined as

$\displaystyle E_{3}=\|\tilde{\theta}-\hat{\theta}\|_{1}.$ (9)

Accordingly, small values of $E_{3}$ indicate that the recovered parameter
estimates are close to the true parameter values.

##  3 Results

###  3.1 Structural Identifiability

We illustrate that the LV model is structurally identifiable. Using the Taylor
series expansion method discussed in Section 2.4, we evaluate
$y(t;\theta)=\\{S,R\\}$ and its successive time derivatives in terms of the
model parameters and initial conditions at time $t=0$, assuming that the
observables $S(t)$ and $R(t)$ are analytic in a neighborhood of $t=0^{+}$.
Recall, these derivatives represent the coefficients of the Taylor series
expansion (up to a constant multiple) of the observable. If we introduce the
following notation,

$a_{0}=S(0^{+}),\ \ a_{1}=S^{\prime}(0^{+}),\ \
a_{2}=S^{\prime\prime}(0^{+}),\ \ a_{3}=S^{\prime\prime\prime}(0^{+}),$
$b_{0}=R(0^{+}),\ \ b_{1}=R^{\prime}(0^{+}),\ \
b_{2}=R^{\prime\prime}(0^{+}),\ \ b_{3}=R^{\prime\prime\prime}(0^{+}),$

then it is straightforward to show that Equations (1)-(2) supply the following
identities:

$\displaystyle\frac{a_{1}}{a_{0}}$ $\displaystyle=$ $\displaystyle
r_{S}-a_{0}\left(\frac{r_{S}}{K_{S}}\right)-b_{0}\left(\frac{r_{S}\gamma_{R}}{K_{S}}\right),$
(10) $\displaystyle\frac{b_{1}}{b_{0}}$ $\displaystyle=$ $\displaystyle
r_{R}-b_{0}\left(\frac{r_{R}}{K_{R}}\right)-a_{0}\left(\frac{r_{R}\gamma_{S}}{K_{R}}\right).$
(11)

We view Equations (10)-(11) as two algebraic equations that relate the unknown
model parameters to measurable quantities $a_{i}$ and $b_{i}$
($i=0,1,2,\dots$). Since we seek to identify six parameters, at least four
more equations are needed. We obtain these additional equations by
differentiating with respect to time and substituting for $a_{i}$ and $b_{i}$.
In this way, we obtain the following equations:

$\displaystyle\frac{a_{2}}{a_{1}}$ $\displaystyle=$ $\displaystyle
r_{S}-2a_{0}\left(\frac{r_{S}}{K_{S}}\right)-\left(\frac{a_{0}b_{1}}{a_{1}}+b_{0}\right)\left(\frac{r_{S}\gamma_{R}}{K_{S}}\right)$
(12) $\displaystyle\frac{b_{2}}{b_{1}}$ $\displaystyle=$ $\displaystyle
r_{R}-2b_{0}\left(\frac{r_{R}}{K_{R}}\right)-\left(a_{0}+\frac{a_{1}b_{0}}{b_{1}}\right)\left(\frac{r_{R}\gamma_{S}}{K_{R}}\right)$
(13) $\displaystyle\frac{a_{3}}{a_{2}}$ $\displaystyle=$ $\displaystyle
r_{S}-2\left(a_{0}+\frac{a_{1}^{2}}{a_{2}}\right)\left(\frac{r_{S}}{K_{S}}\right)-\left(\frac{a_{0}b_{2}}{a_{2}}+\frac{2a_{1}b_{1}}{a_{2}}+b_{0}\right)\left(\frac{r_{S}\gamma_{R}}{K_{S}}\right)$
(14) $\displaystyle\frac{b_{3}}{b_{2}}$ $\displaystyle=$ $\displaystyle
r_{R}-2\left(b_{0}+\frac{b_{1}^{2}}{b_{2}}\right)\left(\frac{r_{R}}{K_{R}}\right)-\left(a_{0}+\frac{2a_{1}b_{1}}{b_{2}}+\frac{a_{2}b_{0}}{b_{2}}\right)\left(\frac{r_{R}\gamma_{S}}{K_{R}}\right)$
(15)

If we make the substitution

$\displaystyle p_{1}=r_{S},\ \ \ \ p_{2}=r_{R},\ \ \ \
p_{3}=\frac{r_{S}}{K_{S}},\ \ \ \ p_{4}=\frac{r_{R}}{K_{R}},\ \ \ \
p_{5}=\frac{r_{s}\gamma_{R}}{K_{S}},\ \ \ \
p_{6}=\frac{r_{R}\gamma_{S}}{K_{R}},$ (16)

then it is straightforward to show that Equations (10)-(15) admit a unique
solution for the six parameters $\\{p_{1},p_{2},p_{3},p_{4},p_{5},p_{6}\\}$
provided that all initial conditions are known:

$\displaystyle p_{1}$ $\displaystyle=$
$\displaystyle-\frac{a_{0}^{3}a_{2}b_{2}-a_{0}^{3}a_{3}b_{1}-2a_{0}^{2}a_{1}^{2}b_{2}+4a_{0}^{2}a_{1}a_{2}b_{1}+a_{0}^{2}a_{1}a_{3}b_{0}-a_{0}^{2}a_{2}^{2}b_{0}-2a_{0}a_{1}^{3}b_{1}-2a_{0}a_{1}^{2}a_{2}b_{0}+2a_{1}^{4}b_{0}}{a_{0}^{3}(a_{1}b_{2}-a_{2}b_{1})},$
$\displaystyle p_{2}$ $\displaystyle=$
$\displaystyle\frac{a_{0}b_{0}^{2}b_{1}b_{3}-a_{0}b_{0}^{2}b_{2}^{2}-2a_{0}b_{0}b_{1}^{2}b_{2}+2a_{0}b_{1}^{4}-a_{1}b_{0}^{3}b_{3}+4a_{1}b_{0}^{2}b_{1}b_{2}-2a_{1}b_{0}b_{1}^{3}+a_{2}b_{0}^{3}b_{2}-2a_{2}b_{0}^{2}b_{1}^{2}}{b_{0}^{3}(a_{1}b_{2}-a_{2}b_{1})},$
$\displaystyle p_{3}$ $\displaystyle=$
$\displaystyle-\frac{a_{0}^{2}a_{2}b_{2}-a_{0}^{2}a_{3}b_{1}-a_{0}a_{1}^{2}b_{2}+3a_{0}a_{1}a_{2}b_{1}-2a_{1}^{3}b_{1}}{a_{0}^{3}(a_{1}b_{2}-a_{2}b_{1}},$
$\displaystyle p_{4}$ $\displaystyle=$
$\displaystyle-\frac{a_{1}b_{0}^{2}b_{3}-3a_{1}b_{0}b_{1}b_{2}+2a_{1}b_{1}^{3}-a_{2}b_{0}^{2}b_{2}+a_{2}b_{0}b_{1}^{2}}{b_{0}^{3}(a_{1}b_{2}-a_{2}b_{1}},$
$\displaystyle p_{5}$ $\displaystyle=$
$\displaystyle-\frac{a_{0}^{2}a_{1}a_{3}-a_{0}^{2}a_{2}^{2}-2a_{0}a_{1}^{2}a_{2}+2a_{1}^{4}}{a_{0}^{3}(a_{1}b_{2}-a_{2}b_{1})},$
$\displaystyle p_{6}$ $\displaystyle=$
$\displaystyle\frac{b_{0}^{2}b_{1}b_{3}-b_{0}^{2}b_{2}^{2}-2b_{0}b_{1}^{2}b_{2}+2b_{1}^{4}}{b_{0}^{3}(a_{1}b_{2}-a_{2}b_{1})}.$

The original parameter set
$\\{r_{S},r_{R},K_{S},K_{R},\gamma_{S},\gamma_{R}\\}$ can be recovered using
Equation (16). Since Equations (10)-(15) admit a unique solution for the full
parameter set, we conclude that the model is structurally identifiable given
observables $y(t;\theta)=\\{S,R\\}$ and known initial conditions
$\\{a_{0},a_{1},a_{2},a_{3},b_{0},b_{1},b_{2},b_{3}\\}$.

The structural identifiability analysis performed above requires knowledge of
the Type-$S$ and Type-$R$ observables; it assumes that, in addition to data on
total tumor volumes, the proportions of the two cell lines are also known. In
practice, it may not be feasible to collect such data in a clinical setting.
Thus, it is natural to ask whether the system remains structurally
identifiable if the total tumor volume, $V=S+R$, is the only observable
quantity. It is possible to show that a unique parameterization exists in this
case, but only if the initial volumes of the two cell lines, $S(0)$ and
$R(0)$, are known. Given these initial conditions, structural identifiability
can be established by proceeding to 6th order in the Taylor series expansion
for $\frac{dV}{dt}=\frac{dS}{dt}+\frac{dR}{dt}$ (results not shown).

To verify the results from the Taylor series method, we use the GenSSI toolbox
for both cases: (1) assuming observables $y(t;\theta)=\\{S,R\\}$ with
available information about $S(0)$, $R(0)$, and their derivatives, and (2)
assuming only a single observable $y(t;\theta)=\\{V=S+R\\}$ with initial
estimates of both cell lines, $S(0)$ and $R(0)$. For the first case, GenSSI
confirms that the model is globally structurally identifiable. For the second
case, with a single observable, $V=S+R$, GenSSI confirms that the model is, at
least, locally structurally identifiable.

###  3.2 Practical Identifiability Using Synthetic LV Data

Having established that the LV model parameters are structurally identifiable,
we focus now on assessing practical identifiability. We begin by attempting to
fit the LV model to noisy data generated from the LV model using the three
methods discussed in Section 2.5.

#### Individual Calibration

The results presented in Figure 4, show how the metrics $E_{1}$, $E_{2}$, and
$E_{3}$—as defined in Section 2.6—change as the initial proportion of Type-$S$
cells varies, when we use the individual calibration procedure. Figure 4(a)
shows that, in most cases, the model fit to the data predicts the correct
interaction type, as indicated by the signs of the inferred interaction
parameters. The only exception is the competitive case with an initial ratio
of Type-$S$ to Type-$R$ cells of 5:5, for which an antagonistic relationship
is predicted. However, the posterior parameter densities for this case reveal
that the carrying capacity and interaction parameters are uninformed by the
data, suggesting that the chosen parameter estimates are non-unique and that
the parameters are unidentifiable (see Figure 5).

(a)

(b)

(c)

Figure 4: (Individual Calibration) (a) Metric $E_{1}$ in Eq. (7), stating the
inferred interaction type (green indicates a match between the interaction
type used to generate the data and that inferred by the model), (b) metric
$E_{2}$ in Eq. (8), measuring the log of sum of squares error across all nine
data sets, and (c) metric $E_{3}$ in Eq. (9), measuring the distance from the
true parameters, for the competitive, mutual, and antagonistic cases with
calibration using a single data set.
Figure 5: (Individual Calibration) Posterior parameter densities computed by
Metropolis Hastings algorithm for the 5:5 initial ratio of Type-$S$ to
Type-$R$ cells.

Results for the prediction error metric $E_{2}$ are presented in Figure 4(b).
Errors in the model fits are generally highest for the mutual interaction data
sets, although even the worst case scenario yields a reasonable visual fit—see
Figure 6 for a comparison of the case with the smallest error (competitive
interaction, with an initial ratio of 1:9) and the case with the largest error
(mutual interaction, with an initial ratio of 4:6). Results for the parameter
error metric $E_{3}$ are presented in Figure 4(c). In general, high values of
$E_{2}$ correspond to parameter estimates that differ markedly from those used
to generate the data, as indicated by large values of $E_{3}$. While it is
unsurprising that the errors are large when the original parameter values
cannot be recovered, these results indicate which experimental data sets are
difficult to fit.

(a) (b)

Figure 6: (Individual Calibration) (a) Competitive case model fits using the
parameter values obtained from the 1:9 experiment, and (b) Mutual case model
fits using the parameter values obtained from the 4:6 experiment.

#### Sequential Calibration

The $E_{1}-E_{3}$ metrics for the sequential calibration procedure are
reported in Figure 7(a)-(c), respectively. As a reminder, these fits are
generated using three data sets: a pure data set containing only Type-$S$
cells, a pure data set containing only Type-$R$ cells, and a mixture data set
with a known value of the initial ratio of Type-$R$ to Type-$S$ cells (the
$x$-axes in Figure 7 indicate the value of this initial ratio). In all cases,
the interaction type is correctly inferred. We observe that while the errors
in the model fits for the mutual data sets are comparable to those for the
individual procedure, the errors for the competitive and antagonistic cases
are markedly reduced. Further, the difference between the true and estimated
parameter values are smaller for the sequential procedure than for the
individual procedure across all three interaction types.

(a)

(b)

(c)

Figure 7: (Sequential Calibration) (a) Metric $E_{1}$ in Eq. (7), stating the
inferred interaction type (green indicates a match between the interaction
type used to generate the data and that inferred by the model), (b) metric
$E_{2}$ in Eq. (8), measuring the log of sum of squares error across all nine
data sets, and (c) metric $E_{3}$ in Eq. (9), measuring the distance from the
true parameters, for the competitive, mutual, and antagonistic cases with
calibration using a single data set.

#### Parallel Calibration

For the parallel scheme, all six parameters are estimated simultaneously using
two different mixture data sets. The results for metrics $E_{1}-E_{3}$ are
displayed in Figure 8(a)-(c) for all three interaction types. The $x$\- and
$y-$axes indicate the initial proportions of Type-$S$ cells employed in the
two mixtures. As indicated by metric $E_{1}$, the inferred interaction matches
the interaction type used to generate the data in all cases. In general,
results for the $E_{2}$ and $E_{3}$ metrics are comparable to those from the
sequential procedure, although certain mixture combinations yield better
results than others. In particular, without knowing a priori the interaction
type of the observed data, calibrating the model parameters using two
experimental mixtures at the extremes—i.e., 1:9 and 9:1—would be a reasonable
strategy for ensuring reasonable errors in model fits and estimated parameter
values that are close to their true values.

(a) (b) (c)

Figure 8: (Parallel Calibration) (a) Inferred interaction type (metric $E_{1}$
in Eq. (7), left), (b) errors in model calibration across all nine data sets
(metric $E_{2}$ in Eq. (8), middle) and (c) distance from true parameter set
(metric $E_{3}$ in Eq. (9), right) for competitive (top), mutual (middle), and
antagonistic (bottom) interaction types.

#### Comparison of Calibration Procedures

To further compare the three calibration methods, we present sample parameter
posterior densities in Figure 9, using the 1:9 ratio mixture for the
individual calibration, the pure data sets and 1:9 ratio mixture for the
sequential method, and the 1:9 and 9:1 mixtures for the parallel calibration.
We focus on the antagonistic case. Our first observation is that the posterior
densities from the individual calibration are highly uninformed by the data,
particularly those for the carrying capacity and interaction parameters
(further reflecting the trend noted in Figure 5). This is indicative of
parameter non-identifiability, in the sense that parameters may not be
uniquely inferred from the data. On the contrary, the sequential and parallel
posterior densities are well-informed by the data, as illustrated by the
narrower distributions and clear maximum a posteriori (or MAP) point
estimates. It is perhaps unsurprising that the posterior densities for the
sequential procedure are the narrowest, since it employs three data sets
whereas the parallel procedure uses only two. Further, the parameters are
estimated two-at-a-time in the sequential procedure, reducing the total
variability at any given step.

Figure 9: Comparison of posterior densities for three calibration methods
(individual using 1:9, sequential using two pure data sets and 1:9, and
parallel using 1:9 and 9:1 mixtures) for the antagonistic interaction type.
True values used to generate the data are depicted with a dashed black line.

Even so, the sequential procedure has drawbacks. While the sequential
procedure ensures the practical identifiability of the LV model parameters
when fitting to data generated from the LV model, the accuracy of inferred
parameters may deteriorate when data is error-prone. In particular, using a
homogeneous data set with large noise level can bias the intrinsic growth
rates and carrying capacities for one or both cell lines; this bias must then
be compensated for in the later inference of the interaction parameters. To
illustrate this drawback, we fit the LV model to 100 sets of synthetic LV data
for four different noise levels and the antagonistic interaction type. Figure
10 summarizes how the inferred values of the interaction parameters,
$\hat{\gamma}_{S}$ and $\hat{\gamma}_{R}$, change as the noise level varies on
the set $\\{5\%,10\%,20\%,50\%\\}$ and the initial ratio of Type-$S$ to
Type-$R$ cells varies on the set $\\{$1:9, 5:5, 9:1$\\}$. The accuracy of the
estimated interaction parameter $\hat{\gamma}_{S}$ depends on the initial
proportion of Type-$S$ cells. When the initial proportion of Type-$S$ cells is
large (e.g., 90%), the inferred interaction parameter values are more accurate
and robust despite the noise, even for the largest noise level of 50%.
However, when the initial proportion of Type-$S$ cells is small (e.g., 10%),
the variance of the estimated $\hat{\gamma}_{S}$ is significantly larger and
the accuracy deteriorates as the noise level increases. Thus, an experimental
design using small initial proportions of Type-$S$ is not an appropriate
design with which to accurately infer parameter values in cases where Type-$R$
cells antagonize Type-$S$ cells when performing a sequential calibration. When
a large proportion of the antagonizing population (e.g. Type-$R$ in case of
Type-$R$ antagonizing Type-$S$) is mixed with a small amount of the
antagonized population, the antagonizing population tends to dominate quickly.
The resulting dynamics are visually similar to those for a competitive
interaction with small Type-$S$ initial proportion; thus, the inferred
interaction type often manifests as competitive, with $\hat{\gamma}_{S}>0$. As
the initial proportion of Type-$S$ cells increases, the antagonizing effect of
Type-$R$ is better recognized, since the Type-$R$ ultimately dominates despite
its small initial proportion. This is a major drawback of sequentially
calibrating the model, since the interaction type may be inconsistently
inferred across different values of the initial ratio of the Type-$S$ and
Type-$R$ cells if the measurement noise is large.

Figure 10: (Sequential Calibration) Comparing fitted values of the interaction
parameters $\gamma_{S}$ and $\gamma_{R}$ when Type-$R$ cells antagonize
Type-$S$ cells, using 100 sets each of synthetic data generated from the LV
model with noise intensity levels of 5%, 10%, 20%, and 50%. The true values of
$\gamma_{S}=-0.5$ and $\gamma_{R}=0.5$ are marked as red dashed lines.

In contrast, the parallel calibration procedure is less sensitive to
measurement noise, since all parameters are estimated simultaneously and two
initial mixture ratios are employed. Interaction parameter estimates for the
parallel model fits to 100 sets of synthetic LV data with initial ratios 1:9
and 9:1 are shown in Figure 11. For the first three noise levels—5%, 10%, and
20%—all parameter estimates correctly predict an antagonistic relationship.
Only for the largest noise levels, 50%, do we observe positive estimates for
$\gamma_{S}$ in nine of the 100 cases, so that the interaction manifests as
competitive according to the LV model definition. The remaining 91 interaction
types are inferred as antagonistic. The increased consistency in the inferred
interaction type is a benefit of using the parallel procedure. Additionally,
when comparing the sequential results in Figure 10 to the parallel results in
Figure 11, we observe that the variability in the parameter estimates across
the 100 fits is much smaller for the parallel calibration, since compensation
for the bias in the error-prone data can be absorbed across all six parameters
instead of isolated to the interaction parameters alone. This phenomenon can
also be observed in Figure 9, where the point estimates associated with the
parallel procedure are closer to the true parameter values than those for the
sequential method ($E_{3}=0.03724$ for the parallel procedure using 1:9 and
9:1 mixtures, and $E_{3}=0.07628$ for the sequential method using the 1:9
mixture).

Figure 11: (Parallel Calibration) Comparing fitted values of the interaction
parameters $\gamma_{S}$ and $\gamma_{R}$ when Type-$R$ cells antagonize
Type-$S$ cells, using 100 sets each of synthetic data generated from the LV
model with 1:9 and 9:1 mixtures and noise intensity levels of 5%, 10%, 20%,
and 50%. The true values of $\gamma_{S}=-0.5$ and $\gamma_{R}=0.5$ are marked
as red dashed lines.

Thus, the choice of optimal method depends on the question of interest. If the
goal is to minimize the uncertainty in the parameter estimates, a sequential
experimental design may be preferable due to the narrower posterior densities
obtained from the use of three data sets, although biases in the values of
parameters estimated early in the sequential design may limit its efficacy
depending on the chosen initial ratio. If minimizing experimental costs is
important, then the parallel procedure is preferred, as it requires two data
sets rather than three. Given that obtaining clinical or in vitro data for
tumor growth measurements can be technically challenging and expensive, and
given that the parallel procedure is more robust to noise and provides point
estimates which are closest to the true parameter values, we focus on this
method for the rest of the investigation.

###  3.3 Practical Identifiability Using Synthetic Cellular Automaton Data

As the parallel procedure was chosen as the preferred calibration methodology
in Section 3.2, henceforth we focus on results for the parallel calibration to
the CA data. Results of the sequential calibration to synthetic CA data are
included in B. Sequential Calibration to Synthetic CA Data for comparison.

We begin by fitting the LV model to CA simulations generated using a neutral
interaction, as described in Section 2.2 and A. Cellular Automaton Model. The
results for the $E_{1}$ and $E_{2}$ metrics are shown in Figure 12. Notably,
nearly all LV model fits infer a competitive interaction type, with the
exception of one inferred antagonistic relationship. This discrepancy is a
drawback of employing a spatially-averaged model (Lotka-Volterra) to fit data
generated from a spatially-resolved model (cellular automaton). In the
spatially-resolved CA model, any imposed interaction mechanism is enforced
atop a baseline competition for space and resources (e.g., oxygen and
nutrients), which cannot be separately accounted for in the spatially-averaged
LV model defined in Equations (1)-(2). Thus, caution is needed when
interpreting the results of the interaction inference for the remainder of
this section.

(a)

(b)

Figure 12: (Parallel Calibration) (a) Metric $E_{1}$ in Eq. (7) and (b) metric
$E_{2}$ in Eq. (8) for the calibration of the Lotka-Volterra model to
synthetic CA data for the neutral interaction type.

Metric $E_{1}$ for the competitive, mutual, and antagonistic interaction types
is illustrated in Figure 13. Recall, when calibrating with the parallel
procedure using data generated from the LV model, the correct interaction type
was inferred in all cases. Here, with calibration to the CA data, the signs of
the inferred interaction parameters in the LV model may not match those in the
CA model, as expected given the results for the neutral case presented in
Figure 12. Although all competitive CA simulations are represented by
competitive interactions in the LV model fits, such a match is only observed
for mutual and antagonistic simulations if the initial proportions of the two
mixtures are at the far extremes (e.g., 1:9 and 9:1). Further, CA simulations
with mutual interactions are often represented by antagonistic interactions
when fit to the LV model. Similarly, CA simulations of antagonistic cases are
often represented by competitive interactions in the LV fits.

(a)

(b)

(c)

Figure 13: (Parallel Calibration) Metric $E_{1}$ in Eq. (7) for the
calibration of the Lotka-Volterra model to synthetic CA data for (a)
competitive, (b) mutual, and (c) antagonistic interaction types, using CA
intensity parameter $I=4$. Green indicates mixture combinations for which
inferred interaction type in the LV model matches the interaction type used to
generate the CA data.

To further quantify the bias toward a competitive interaction in the LV model,
we present the estimated values of the interaction parameters for the four
interaction types in Figure 14, demonstrating their dependence on both the
interaction type imposed in the CA and the first mixture ratio employed. For
each first mixture ratio choice, we average the parameter estimates for the
eight remaining second mixture ratios for easier visualization. Consistent
with Figure 12, the parameter estimates for $\gamma_{S}$ and $\gamma_{R}$ in
the neutral case tend to positive values, such that the LV model infers a
competitive interaction; in fact, the estimated parameter values for the
neutral and competitive cases are similar in magnitude. Additionally, it is
rare for both interaction parameters in LV model fits to CA simulations with
mutual interactions to be less than zero simultaneously, as was expected in
our earlier idealistic calibration to spatially-averaged data. Thus, LV model
fits to mutual CA simulations often infer antagonistic interactions. In
general, applying a spatially-averaged model to fit spatially-resolved data
can lead to a misinterpretation of the type of interaction present, as one is
unable to distinguish between the effects of the interaction mechanism and the
effects of the baseline competition for space and resources; both effects are
absorbed into the values of the interaction parameters $\gamma_{S}$ and
$\gamma_{R}$, and cannot be separately identified.

(a) (b)

Figure 14: (Parallel Calibration) Interaction parameter estimates for (a)
$\gamma_{S}$ and (b) $\gamma_{R}$ according to first mixture Type-$S$:Type-$R$
ratio and interaction type, using CA intensity parameter $I=4$. Each data
point represents the average interaction parameter over the eight remaining
second mixture ratios.

Based on this result, we hypothesize that interaction inference in the LV
model may depend upon the level of intensity of the interaction mechanism
imposed between the two cell lines in the CA model. Figure 15 illustrates the
dependence of metric $E_{1}$ on the value of the intensity parameter, $I$,
used to generate the CA data. For $I=3$ in the antagonistic case, we observe
that there is no combination of mixtures for which the interaction inferred
from the LV model matches the interaction type of the CA simulations. As the
intensity level increases, we observe that using two extreme mixture
combinations enables correct inference of the interaction type used to
generate the CA data. We conclude that the interaction inferred by the LV
model will only match that used in the CA simulation if the intensity level in
the CA model is sufficiently strong.

Intensity $I=3$

t  Intensity $I=4$ Intensity $I=5$
(a)
(b)
(c)

Figure 15: (Parallel Calibration) Comparing LV model interaction inference
across different CA interaction intensity levels using parallel calibration
(left to right: $I=3$, $I=4$, and $I=5$) for (a) competitive, (b) mutual, and
(c) antagonistic. Green indicates mixture combinations for which inferred
interaction type in the LV model matches the interaction type used to generate
the CA data.

Figure 16 details the $E_{2}$ error metric for the competitive, mutual, and
antagonistic interaction types. In general, the errors in the model fits are
much larger—particularly in the competitive and antagonistic cases—as compared
to the idealistic calibration in Figure 8. For demonstration, we plot the
model fits and associated credible intervals for two of the competitive
scenarios in Figure 17. The first set, using the 5:5 and 7:3 mixtures, is the
best case scenario in terms of the $E_{2}$ metric, with $E_{2}=-2.477.$ The
second set, using the 7:3 and 8:2 mixtures, is the worst case scenario, with
$E_{2}=2.309.$ In both cases, the model fits are reasonable for initial ratios
close to one of the combinations used for calibration, but for other values of
the initial ratio the fits (and uncertainties) are poor. Again, we attribute
this phenomenon to model discrepancy: i.e., our use of a spatially-averaged
model to fit spatially-resolved data. Comparing this figure to Figure 6, where
the data was generated from the LV model, reveals a major discrepancy in the
behavior of the data. Generation of competitive interaction data from the LV
model results in asymptotic coexistence between the two cell lines. However,
the data shown in Figure 17 tends toward elimination of one cell line over
time; the additional competition for space and resources inhibits coexistence.
The inability of the LV model to account for this added competition is an
issue of model discrepancy, leading to larger model errors and uncertainties
when calibrating to the synthetic CA data.

(a)

(b)

(c)

Figure 16: (Parallel Calibration) Metric $E_{2}$ in Eq. (8) for the
calibration of the Lotka-Volterra model to synthetic CA data for (a)
competitive, (b) mutual, and (c) antagonistic interaction types, using CA
intensity parameter $I=4$.

(a) (b)

Figure 17: (Parallel Calibration) Competitive model fits and credible
intervals using (a) parameter values obtained from the 5:5 and 7:3
experiments, and (b) parameter values obtained from the 7:3 and 8:2
experiments, with synthetic data generated from the CA model with intensity
parameter $I=4$.

##  4 Discussion

In this work, we investigated whether the type of interaction between two cell
lines in a tumor spheroid (competitive, mutualistic, or antagonistic) can be
inferred using the Lotka-Volterra model. We verified the structural
identifiability of the system, and considered the practical identifiability of
the model for three different experimental designs and calibration procedures.
The first procedure tested whether the parameters of the Lotka-Volterra model
could be uniquely inferred using a single data set containing information
about both cell line volumes over time. The second design was a sequential
calibration procedure, using two homogeneous cell-line experiments to
calibrate the growth rate and capacity parameters and a mixture data set to
calibrate the interaction parameters. The final scheme was a parallel
calibration of all parameters using two mixture data sets.

By testing the three experimental design procedures on data generated from the
Lotka-Volterra model, we determined that both the sequential and parallel
designs yield identifiable parameter sets, while the posterior densities
resulting from the individual calibration scheme were not well-informed by the
data. Though the sequential procedure produced the narrowest posteriors, the
parallel procedure was shown to be more robust to noise and recovered point
estimates closest to the true parameters. Thus, we favored the parallel
procedure since it could be employed with one fewer data set, an important
consideration for applications with budgetary restrictions on experiments.

To assess the robustness of our parallel calibration procedure on a more
realistic data set, we developed a cellular automaton model that tracks two
cell types exhibiting these distinct types of interactions with parameters
estimated from a prostate cancer cell line. The imposed interaction type
dictated the speed at which cells divide and their corresponding oxygen
consumption, when surrounded by a sufficient number of cells of the opposite
type. We used this model to generate in silico data, to which we calibrated
the Lotka-Volterra model using the proposed parallel calibration design. The
use of the spatially-averaged Lotka-Volterra model to fit data generated from
the spatially-resolved cellular automaton model uncovered a model discrepancy
issue; namely, the effects of both the primary interaction mechanism between
the two cell lines and the baseline competition for space and resources were
both absorbed into the Lotka-Volterra interaction parameters $\gamma_{S}$ and
$\gamma_{R}$, with the result that the inferred interaction type by the Lotka-
Volterra model would only match that used to generate the CA data if the
intensity level of the interaction—as specified in the CA data generation—was
strong enough to overcome the underlying competitive nature of the spatially-
resolved data. With a strong enough intensity level, this match generally
occurs for all interaction types when two asymmetric initial proportions are
used for calibration, but is not guaranteed for other mixture combinations.
Thus, researchers are cautioned to treat their interpretation of the Lotka-
Volterra interaction parameters with care, as interaction types inferred from
spatially-averaged data may not accurately reflect the interaction type of a
spatially-resolved experiment.

This work serves as a preliminary investigation into the identifiability of
the Lotka-Volterra model using in silico data, with special emphasis given to
distinguishing the interaction type between the two cell lines. A number of
possible directions are available from here. First, it remains to be seen
whether there exists a correlation between the initial spatial distribution of
the two cell lines and the type of interactions that result. A natural follow-
up would be to expand the Lotka-Volterra into PDE form to account for spatial
variation, and to determine whether there exists a relationship between
spatial distribution of cell types and the model’s ability to infer the
interaction mechanism. This investigation could be performed using both the
current version of the cellular automaton model described in Section 3.3 and
A. Cellular Automaton Model, which mimics two-dimensional growth such as that
observed in in vitro laboratory settings, but could also be tested using a
scaled-up version of the CA model in three-dimensions, a closer representation
of in vivo tumor composition. Eventually, we plan to verify these results
using experimental data for two different cell lines.

Additionally, the analysis in this study is a first step towards the
incorporation of radiotherapy treatment, and potentially other types of
treatment. Specifically, we are interested in studying two cell types with
differing intrinsic characteristics and different levels of radiosensitivity
(e.g., Type-$S$ for “sensitive” and Type-$R$ for “resistant”), to determine
whether these inherent differences in the two cell lines have an impact on the
inferred interaction type when the Lotka-Volterra model is fit to experimental
data.

## References

  * Altrock et al. [2015] Altrock, P.M., Liu, L.L., Michor, F., 2015.  The mathematics of cancer: integrating quantitative models.  Nature Reviews Cancer 15, 730–745.  doi:10.1038/nrc4029.
  * Balsa-Canto et al. [2010] Balsa-Canto, E., Alonso, A.A., Banga, J.R., 2010.  An iterative identification procedure for dynamic modeling of biochemical networks.  BMS Systems Biology 4.  doi:doi: 10.1186/1752-0509-4-11.
  * Bates et al. [2019] Bates, D.J., Hauenstein, J.D., Meshkat, N., 2019.  Identifiability and numerical algebraic geometry.  PLoS One 14.  doi:doi: 10.1371/journal.pone.0226299.
  * Benzekry et al. [2014] Benzekry, S., Lamont, C., Beheshti, A., Tracz, A., Ebos, J.M., Hlatky, L., Hahnfeldt, P., 2014.  Classical mathematical models for description and prediction of experimental tumor growth.  PLoS computational biology 10, e1003800.
  * Byrne [2010] Byrne, H.M., 2010.  Dissecting cancer through mathematics: From the cell to the animal model.  Nature Reviews Cancer 10, 221–230.  doi:10.1038/nrc2808.
  * Caswell and Swanton [2017] Caswell, D.R., Swanton, C., 2017\.  The role of tumour heterogeneity and clonal cooperativity in metastasis, immune evasion and clinical outcome.  BMC medicine 15, 1–9.
  * Chapman et al. [2014] Chapman, A., del Ama, L.F., Ferguson, J., Kamarashev, J., Wellbrock, C., Hurlstone, A., 2014\.  Heterogeneous tumor subpopulations cooperate to drive invasion.  Cell reports 8, 688–695.
  * Chappell et al. [1990] Chappell, M.J., Godfrey, K.R., Vajda, S., 1990.  Global identifiability of the parameters of nonlinear systems with specified inputs: a comparison of methods.  Mathematical Biosciences 102, 41–73.
  * Chis et al. [2011a] Chis, O.T., Banga, J.R., Balsa-Canto, E., 2011a.  Genssi: a software toolbox for structural identifiability analysis of biological models.  Bioinformatics 27, 2610–2611.
  * Chis et al. [2011b] Chis, O.T., Banga, J.R., Balsa-Canto, E., 2011b.  Structural identifiability of systems biology models: a critical comparison of methods.  PLoS One 6.  doi:doi: 10.1371/journal.pone.0027755.
  * Cho and Levy [2020] Cho, H., Levy, D., 2020.  The impact of competition between cancer cells and healthy cells on optimal drug delivery.  Mathematical Modelling of Natural Phenomena 15, 42.
  * Cho et al. [2020] Cho, H., Lewis, A.L., Storey, K.M., 2020.  Bayesian information-theoretic calibration of radiotherapy sensitivity parameters for informing effective scanning protocols in cancer.  Journal of Clinical Medicine 9.  doi:10.3390/jcm9103208.
  * Cho et al. [2021] Cho, H., Lewis, A.L., Storey, K.M., Jennings, R., Shtylla, B., Reynolds, A.M., Byrne, H.M., 2021.  A framework for performing data-driven modeling of tumor growth with radiotherapy treatment.  Springer Special Issue: Using Mathematics to Understand Biological Complexity, Women in Mathematical Biology , 179–216.
  * Cunningham et al. [2018] Cunningham, J.J., Brown, J.S., Gatenby, R.A., Staňková, K., 2018\.  Optimal control to develop therapeutic strategies for metastatic castrate resistant prostate cancer.  Journal of theoretical biology 459, 67–78.
  * Diego et al. [2013] Diego, D., Calvo, G.F., Pérez-García, V.M., 2013.  Modeling the connection between primary and metastatic tumors.  Journal of mathematical biology 67, 657–692.
  * Gatenby and Gawlinski [1996] Gatenby, R.A., Gawlinski, E.T., 1996\.  A reaction-diffusion model of cancer invasion.  Cancer research 56, 5745–5753.
  * Gerlee [2013] Gerlee, P., 2013.  The model muddle: in search of tumor growth laws.  Cancer research 73, 2407–2411.
  * Ghaffarizadeh et al. [2018] Ghaffarizadeh, A., Heiland, R., Friedman, S.H., Mumenthaler, S.M., Macklin, P., 2018.  Physicell: an open source physics-based cell simulator for 3-d multicellular systems.  PLoS computational biology 14, e1005991.
  * Greene et al. [2019] Greene, J.M., Gevertz, J.L., Sontag, E.D., 2019.  Mathematical approach to differentiate spontaneous and induced evolution to drug resistance during cancer treatment.  JCO Clinical Cancer Informatics doi:doi: 10.1200/CCI.18.00087.
  * Haario et al. [2006] Haario, H., Laine, M., Mira, A., Saksman, E., 2006\.  DRAM: Efficient adaptive MCMC.  Statistics and Computing 16, 339–354.  doi:10.1007/s11222-006-9438-0.
  * Hyun and Macklin [2013] Hyun, A.Z., Macklin, P., 2013\.  Improved patient-specific calibration for agent-based cancer modeling.  Journal of theoretical biology 317, 422.
  * Isea and Lonngren [2015] Isea, R., Lonngren, K.E., 2015\.  A mathematical model of cancer under radiotherapy.  International Journal of Public Health Research 3, 340–344.
  * Kim et al. [2014] Kim, Y., Kang, H., Lawler, S., Eladdadi, A., Kim, P., Mallet, D., 2014.  Mathematical models of tumor-immune system dynamics.
  * Kloppers and Greeff [2013] Kloppers, P., Greeff, J., 2013\.  Lotka–volterra model parameter estimation using experiential data.  Applied Mathematics and Computation 224, 817–825.  doi:https://doi.org/10.1016/j.amc.2013.08.093.
  * Laine [2018] Laine, M., 2018.  MCMC toolbox for Matlab.  URL: https://mjlaine.github.io/mcmcstat/.
  * Laleh et al. [2022] Laleh, N.G., Loeffler, C.M.L., Grajek, J., Staňková, K., Pearson, A.T., Muti, H.S., Trautwein, C., Enderling, H., Poleszczuk, J., Kather, J.N., 2022\.  Classical mathematical models for prediction of response to chemotherapy and immunotherapy.  PLoS Comp Biol 18, e1009822.
  * Liepe et al. [2013] Liepe, J., Filippi, S., Komorowski, M., Stumpf, M., 2013\.  Maximizing the Information Content of Experiments in Systems Biology.  PLoS Computational Biology 9, 1–13.  doi:10.1371/journal.pcbi.1002888.
  * Lotka [1925] Lotka, A.J., 1925.  Elements of Physical Biology.  Williams & Wilkins, Baltimore.
  * Lowengrub et al. [2009] Lowengrub, J.S., Frieboes, H.B., Jin, F., Chuang, Y.L., Li, X., Macklin, P., Wise, S.M., Cristini, V., 2009\.  Nonlinear modelling of cancer: bridging the gap between cells and tumours.  Nonlinearity 23, R1.
  * Macklin et al. [2010] Macklin, P., Edgerton, M., Cristini, V., Cristini, V., Lowengrub, J., 2010.  Multiscale modeling of cancer.  Cristini and Lowengrub , 216–244.
  * Marasco et al. [2016] Marasco, A., Picucci, A., Romano, A., 2016.  Market share dynamics using lotka–volterra models.  Technological Forecasting and Social Change 105, 49–62.  doi:https://doi.org/10.1016/j.techfore.2016.01.017.
  * McGranahan and Swanton [2017] McGranahan, N., Swanton, C., 2017\.  Clonal heterogeneity and tumor evolution: Past, present, and the future.  Cell 168, 613–628.  doi:10.1016/j.cell.2017.01.018.
  * Noble et al. [2021] Noble, R.J., Walther, V., Roumestand, C., Hochberg, M.E., Hibner, U., Lassus, P., 2021\.  Paracrine behaviors arbitrate parasite-like interactions between tumor subclones.  Frontiers in ecology and evolution 9\.
  * Paczkowski et al. [2021] Paczkowski, M., Kretzschmar, W.W., Markelc, B., Liu, S.K., Kunz-Schughart, L.A., Harris, A.L., Partridge, M., Byrne, H.M., Kannan, P., 2021.  Reciprocal interactions between tumour cell populations enhance growth and reduce radiation sensitivity in prostate cancer.  Communications Biology 4, 1–13.
  * Pohjanpalo [1978] Pohjanpalo, H., 1978.  System identifiability based on the power series expansion of the solution.  Mathematical Biosciences 41, 21–33.
  * Poleszczuk et al. [2016] Poleszczuk, J., Macklin, P., Enderling, H., 2016.  Agent-based modeling of cancer stem cell driven solid tumor growth.  Methods Mol Biol , 335–346doi:doi: 10.1007/7651_2016_346.
  * Raue et al. [2014] Raue, A., Karlsson, J., Saccomani, M.P., Jirstrand, M., Timmer, J., 2014.  Comparison of approaches for parameter identifiability analysis of biology systems.  Bioinformatics 30, 1440–1448.  doi:doi: 10.1093/bioinformatics/btu006.
  * Remien and Ridenhour [2021] Remien, C., Ridenhour, B., 2021\.  Structural identifiability of the generalized lotka-volterra model for microbiome studies.  Royal Society Open Science 8.  doi:10.1098/rsos.201378.
  * Sim et al. [2012] Sim, A., Filippi, S., Stumpf, M.P.H., 2012.  Information geometry and sequential monte carlo.  URL: https://arxiv.org/abs/1212.0764, doi:10.48550/ARXIV.1212.0764.
  * Simpson et al. [2022] Simpson, M.J., Browning, A.P., Warne, D.J., Maclaren, O.J., Baker, R.E., 2022.  Parameter identifiability and model selection for sigmoid population growth models.  J Theor Biol. 535.  doi:doi: 10.1016/j.jtbi.2021.110998.
  * Smith [2014] Smith, R., 2014.  Uncertainty Quantification: Theory, Implementation, and Applications.  SIAM, Philadelphia, PA.
  * Susswein et al. [2022] Susswein, Z., Sengupta, S., Clarke, R., Bansal, S., 2022\.  Borrowing ecological theory to infer interactions between sensitive and resistant breast cancer cell populations.  bioRxiv .
  * Tabassum and Polyak [2015] Tabassum, D., Polyak, K., 2015\.  Tumorigenesis: it takes a village.  Nat Rev Cancer 15, 473–483.  doi:10.1038/nrc3971.
  * Van Liedekerke et al. [2015] Van Liedekerke, P., Palm, M., Jagiella, N., Drasdo, D., 2015\.  Simulating tissue mechanics with agent-based models: concepts, perspectives and some novel results.  Computational particle mechanics 2, 401–444.
  * Venugopalan et al. [2014] Venugopalan, G., Camarillo, D.B., Webster, K.D., Reber, C.D., Sethian, J.A., Weaver, V.M., Fletcher, D.A., El-Samad, H., Rycroft, C.H., 2014.  Multicellular architecture of malignant breast epithelia influences mechanics.  PLoS One 9, e101955.
  * Villaverde [2019] Villaverde, A.F., 2019.  Observability and structural identifiability of nonlinear biological systems.  Complexity , 1–12doi:doi: 10.1155/2019/8497093.
  * Volterra [1926] Volterra, V., 1926.  Fluctuations in the abundance of a species considered mathematically.  Nature 118, 558–560.
  * Wieland et al. [2021] Wieland, F.G., Hauber, A.L., Rosenblatt, M., Tönsing, C., Timmer, J., 2021.  On structural and practical identifiability.  Current Opinion in Systems Biology 25, 60–69.
  * Wu et al. [2010] Wu, M., Pastor-Pareja, J.C., Xu, T., 2010.  Interaction between rasv12 and scribbled clones induces tumour growth and invasion.  Nature 463, 545–548.

## Appendix

## A. Cellular Automaton Model

We adapt our CA model from the models in [34, 13, 12], by allowing for
different types of interactions between the two cell types. In the CA model,
the interaction types affect the proliferation and oxygen consumption rates.
Let $\tau$ denote the discrete time step length in the model simulation, and
let $\tau_{cycle}(\mathbf{x},t)$ denote the cell cycle counter for a
proliferating cell at site $\mathbf{x}$ at time $t$. When
$\tau_{cycle}(\mathbf{x},t)\leq 0$, the cell at site $\mathbf{x}$ divides.
During each time step, the reduction in the cycle counter for a given cell and
its oxygen consumption depend on the composition of its neighborhood. For a
proliferating Type-$S$ cell or Type-$R$ cell at site $\mathbf{x}$, its cell
cycle counter reduces by $\Delta\tau_{S}(I)$ or $\Delta\tau_{R}(I)$,
respectively, and it consumes oxygen at rate $\kappa_{S}(I)$ or
$\kappa_{R}(I)$, respectively (mol cm-3 s-1). The values of
$\Delta\tau_{y}(I)$ and $\kappa_{y}(I)$, for $y\in(S,R)$, depend on the local
neighborhood of each cell, on the interaction type that we are modeling, and
on the intensity parameter $I$, which quantifies the intensity of the
interaction between the cell populations. Note that $I\geq 1$, where $I=1$
yields the neutral interaction case in the CA model, and higher values of $I$
correspond to stronger interactions between the cell types. In Section 3.3 of
the main text, we discuss the sensitivity of the Lotka-Volterra interaction
parameter calibration to the CA interaction parameter $I$.

We consider four different interaction types between cell populations:
neutral, competitive, mutualistic, antagonistic (with the Type-$R$ population
antagonizing the Type-$S$ population, and the Type-$S$ population promoting
the Type-$R$ population). In the competitive case, cells surrounded by a large
number of cells of the other type consume less oxygen, and, as a result,
divide more slowly. In the mutual case, cells surrounded by a large number of
cells of the other type consume more oxygen and divide more quickly. In the
antagonistic case, Type-$S$ cells surrounded by large numbers of Type-$R$
cells consume less oxygen and divide more slowly, while Type-$R$ cells
surrounded by large numbers of Type-$S$ cells consume more oxygen and divide
more quickly. In order to summarize the specific details of these effects, we
first define relevant terminology. Let $N(\mathbf{x},t)$ denote the total
number of neighbors in the two-dimensional Moore neighborhood of site
$\mathbf{x}$ at time $t$, with $0\leq N(\mathbf{x},t)\leq 8$. Let
$N_{S}(\mathbf{x},t)$ and $N_{R}(\mathbf{x},t)$ be the number of Type-$S$
cells and Type-$R$ cells, respectively, in the neighborhood of site
$\mathbf{x}$ at time $t$, with $0\leq N_{y}(\mathbf{x},t)\leq
N(\mathbf{x},t)\leq 8$ for $y\in(S,R)$. We denote by $T_{1}$ and $T_{2}$ two
threshold values for the neighborhood size, with $0\leq T_{1}<T_{2}\leq 8$,
which are used to determine the values of $\Delta\tau_{y}(I)$ and
$\kappa_{y}(I)$.

Table 3 summarizes the value for the cell cycle counter reduction for a
Type-$S$ cell, $\Delta\tau_{S}(I)$, dependent upon the total number of
neighbors in its Moore neighborhood, the number of Type-$R$ neighbors, and the
interaction type used in the simulation. We note that there are baseline
levels of cell cycle counter reduction, dependent on the total number of
neighbors, to simulate contact inhibition of proliferation. For interaction
types with $I>1$, additional changes to the cell cycle counter reduction
depend on the number of neighbors of the other type. Similarly to the cell
cycle counter reduction, Table 4 displays the value for the corresponding
oxygen consumption rate, $\kappa_{S}(I)$. Almost all values in Tables 3 and 4
are identical for a Type-$R$ cell, when considering the number of Type-$S$
neighbors in its neighborhood, except the two values highlighted in blue, for
the antagonistic case. In this case, $\Delta\tau_{R}(I)=\frac{\tau I}{2}$ and
$\kappa_{R}(I)=\frac{\kappa_{P}I}{2}$ when $T_{1}<N(\mathbf{x},t)\leq T_{2}$
and $N_{S}(\mathbf{x},t)>T_{1}$, and $\Delta\tau_{R}(I)=\frac{\tau I}{4}$ and
$\kappa_{R}(I)=\frac{\kappa_{P}I}{4}$ when $N(\mathbf{x},t)>T_{2}$ and
$N_{S}(\mathbf{x},t)>T_{2}$.

| Neutral | Competitive | Mutual | Antagonistic
---|---|---|---|---
$N(\mathbf{x},t)\leq T_{1}$ | $\tau$ | $\tau$ | $\tau$ | $\tau$
$T_{1}<N(\mathbf{x},t)\leq T_{2}$ and $N_{R}(\mathbf{x},t)\leq T_{1}$ | $\frac{\tau}{2}$ | $\frac{\tau}{2}$ | $\frac{\tau}{2}$ | $\frac{\tau}{2}$
$T_{1}<N(\mathbf{x},t)\leq T_{2}$ and $N_{R}(\mathbf{x},t)>T_{1}$ | $\frac{\tau}{2}$ | $\frac{\tau}{2I}$ | $\frac{\tau I}{2}$ | $\frac{\tau}{2I}$
$N(\mathbf{x},t)>T_{2}$ and $N_{R}(\mathbf{x},t)\leq T_{2}$ | $\frac{\tau}{4}$ | $\frac{\tau}{4}$ | $\frac{\tau}{4}$ | $\frac{\tau}{4}$
$N(\mathbf{x},t)>T_{2}$ and $N_{R}(\mathbf{x},t)>T_{2}$ | $\frac{\tau}{4}$ | $\frac{\tau}{4I}$ | $\frac{\tau I}{4}$ | $\frac{\tau}{4I}$
Table 3: The value for the cell cycle counter reduction, $\Delta\tau_{S}(I)$, for each neighborhood scenario and interaction type. Alternatively, when considering a resistant cell at site x and the number of sensitive neighbors, the only values that differ for $\Delta\tau_{R}(I)$ are those highlighted in blue. | Neutral | Competitive | Mutual | Antagonistic
---|---|---|---|---
$N(\mathbf{x},t)\leq T_{1}$ | $\kappa_{P}$ | $\kappa_{P}$ | $\kappa_{P}$ | $\kappa_{P}$
$T_{1}<N(\mathbf{x},t)\leq T_{2}$ and $N_{R}(\mathbf{x},t)\leq T_{1}$ | $\frac{\kappa_{P}}{2}$ | $\frac{\kappa_{P}}{2}$ | $\frac{\kappa_{P}}{2}$ | $\frac{\kappa_{P}}{2}$
$T_{1}<N(\mathbf{x},t)\leq T_{2}$ and $N_{R}(\mathbf{x},t)>T_{1}$ | $\frac{\kappa_{P}}{2}$ | $\frac{\kappa_{P}}{2I}$ | $\frac{\kappa_{P}I}{2}$ | $\frac{\kappa_{P}}{2I}$
$N(\mathbf{x},t)>T_{2}$ and $N_{R}(\mathbf{x},t)\leq T_{2}$ | $\frac{\kappa_{P}}{4}$ | $\frac{\kappa_{P}}{4}$ | $\frac{\kappa_{P}}{4}$ | $\frac{\kappa_{P}}{4}$
$N(\mathbf{x},t)>T_{2}$ and $N_{R}(\mathbf{x},t)>T_{2}$ | $\frac{\kappa_{P}}{4}$ | $\frac{\kappa_{P}}{4I}$ | $\frac{\kappa_{P}I}{4}$ | $\frac{\kappa_{P}}{4I}$
Table 4: The value for the oxygen consumption rate, $\kappa_{S}(I)$, for each
neighborhood scenario and interaction type. Alternatively, when considering a
resistant cell at site x and the number of sensitive neighbors, the only
values that differ for $\kappa_{R}(I)$ are those highlighted in blue.

As part of the CA model, we incorporate oxygen, the single growth-rate-
limiting nutrient, using a reaction-diffusion equation. In particular, the
change in oxygen concentration $c(\textbf{x},t)$ (mol cm-3) at location x at
time $t$ is described by:

$\frac{\partial c(\textbf{x},t)}{\partial
t}=D\nabla^{2}c(\textbf{x},t)-\Gamma(\textbf{x},t),$

where $D$ is the oxygen diffusion coefficient (cm2 s-1), and
$\Gamma(\textbf{x},t)$ is the oxygen consumption rate (mol cm-3 s-1). In the
neutral interaction case, $\Gamma(\textbf{x},t)$ is defined as follows:

$\displaystyle\Gamma(\textbf{x},t)=\begin{cases}\kappa_{S}(\textbf{x},t)&\text{
if }\textbf{x}\text{ is occupied by a proliferating Type-$S$ cell}\\\
\kappa_{R}(\textbf{x},t)&\text{ if }\textbf{x}\text{ is occupied by a
proliferating Type-$R$ cell}\\\ \kappa_{Q}&\text{ if }\textbf{x}\text{ is
occupied by a quiescent cell}\\\ 0&\text{ otherwise.}\end{cases}$

Recall $\kappa_{S}(\textbf{x},t)$ and $\kappa_{R}(\textbf{x},t)$ are
summarized in Table 4. The parameter $\kappa_{Q}>0$ denotes the constant rate,
in mol cm-3 s-1, at which quiescent cells consume oxygen. We note that the
concentration $c(\textbf{x},t)$ will not become negative because cells exposed
to oxygen below the threshold $c_{N}$ will become necrotic, and necrotic cells
do not consume oxygen. We also use the following initial and boundary
conditions to simulate oxygen diffusion from the boundaries of a square Petri
dish into the culture medium:

$\displaystyle c(x,y,0)=c_{\infty},$ $\displaystyle
c(0,y,t)=c(L,y,t)=c(x,0,t)=c(x,L,t)=c_{\infty},$

where $L$ is the domain length, and $c_{\infty}$ is the background O2
concentration.

Local oxygen levels determine the cell state at each lattice site. Sites can
contain proliferating cells, quiescent cells, or necrotic cells, or they can
be empty. The threshold oxygen levels $c_{Q}$ and $c_{N}$ determine the state
at each site, as follows:

  * •

If $c_{Q}<c(\textbf{x},t)\leq c_{\infty}$, then the cell at site x is
proliferating.

  * •

If $c_{N}<c(\textbf{x},t)\leq c_{Q}$, then the cell at site x is quiescent.

  * •

If $0\leq c(\textbf{x},t)\leq c_{N}$, then the cell at site x is necrotic.

Table 5 summarizes the parameter values that we used to simulate the CA model
throughout this work.

Param. | Description | Value | Units
---|---|---|---
$l$ | Cell size | 0.0018 | cm
$L$ | Domain length | 0.36 | cm
$\bar{\tau}_{cycle}$ | Mean (SD) cell cycle time | 18.3 (1.4) | h
$c_{\infty}$ | Background O2 concentration | $2.8\times 10^{-7}$ | mol cm-3
$D$ | O2 diffusion constant | $1.8\times 10^{-5}$ | cm2s-1
$c_{Q}$ | O2 concentration threshold for proliferating cells | 1.82$\times 10^{-7}$ | mol cm-3
$c_{N}$ | O2 concentration threshold for quiescent cells | 1.68$\times 10^{-7}$ | mol cm-3
$\kappa_{P}$ | O2 consumption rate of proliferating cells | 1.0$\times 10^{-8}$ | mol cm-3s-1
$\kappa_{Q}$ | O2 consumption rate of quiescent cells | 5.0$\times 10^{-9}$ | mol cm-3s-1
$p_{NR}$ | Rate of lysis of necrotic cells | 0.01 | hr-1
$I$ | Interaction intensity level | Varies | –
$T_{1}$ | Neighborhood threshold 1 for cell cycle reduction | 4 | cells
$T_{2}$ | Neighborhood threshold 2 for cell cycle reduction | 7 | cells
Table 5: A summary of the parameters used in the CA model and their default
values. Parameter values are estimated using experimental data from the
prostate cancer cell line, PC3, in [34].

We generate the in silico data from the CA model with different initial
proportions of Type-$S$ cells, varying this initial proportion across all 0.1
increments between 0 and 1. We generate data exhibiting each of the four
interaction types: neutral, competitive, mutual, and antagonistic. In all
cases, we initialize the cells in a circular region of radius 9 cells, in the
center of a domain of size 200 cells $\times$ 200 cells. The type of cell at
each initial site is randomly chosen, according to the specified initial
distribution. All cells are initially proliferating cells, each with an
initial cell cycle counter $\tau_{cycle}(\mathbf{x},0)$ chosen as a random
integer between 0 and the mean cell cycle time, $\bar{\tau}_{cycle}$.

## B. Sequential Calibration to Synthetic CA Data

Here, we present the sequential calibration results when fitting the LV model
to synthetic CA data, for comparison to the results from the parallel
procedure detailed in Section 3.3. Metric $E_{1}$, shown in Figure 18(a),
reveals that for antagonistic CA simulations, there is inconsistency in the
type of interaction inferred by the LV model. Specifically, for Type-$S$
initial proportions from 0.1 to 0.6, the LV model predicts a competitive
interaction, while for larger initial proportions of Type-$S$, an antagonistic
relationship is inferred. Metric $E_{2}$ is shown in Figure 18(b); in general,
the errors are larger than those computed for the parallel calibration scheme
in Figure 8, further motivating our choice to use the parallel procedure in
Section 3.3. The parameter estimates for each of the interaction types and
initial ratios are reported in Figure 18(c). We note that the interaction
parameters associated with mutualistic interactions appear to vanish to zero,
but they are actually slightly negative in value. As discussed in Section 3.3,
the tendency of the interaction parameters to skew toward positive values is a
result of the baseline competition for space and resources in the spatially-
resolved CA data. As such, the LV model often predicts a competitive
relationship for antagonistic CA simulations. In particular, if the Type-$R$
cells antagonize the Type-$S$ cells and a large number of Type-$R$ cells are
mixed with a small number of Type-$S$ cells, then the type of interaction
tends to be misclassified. Since the initial fraction of the Type-$S$ cells is
small, it is likely that the Type-$R$ cells rapidly become dominant and the
Type-$S$ cells are rapidly eliminated. Thus, when the antagonized population
Type-$S$ comprise only 10% of the initial tumor volume, the interaction
parameters are estimated to have the same signs as for the competitive
interaction case. However, as the initial proportion of Type-$S$ cells
increases, the antagonizing effect of the Type-$R$ cells is discernible, since
the Type-$R$ ultimately dominates despite its small initial portion.
Specifically, we find that the interaction parameter $\gamma_{S}$ is negative
when the initial proportion of Type-$S$ cells exceeds 0.6.

(a)

(b)

(c)

Figure 18:  (Sequential Calibration) (a) Metric $E_{1}$ in Eq. (7) for the
sequential calibration of the LV model to synthetic CA data. Green indicates
mixtures for which the inferred interaction type in the LV model matches the
interaction type used to generate the CA data. (b) Metric $E_{2}$ in Eq. (8),
measuring the log sum-of-squares error across all nine model fits. (c) The
values of fitted parameters $\hat{\gamma}_{S}$ and $\hat{\gamma}_{R}$ for each
interaction type across different initial proportions of Type-$S$ cells.