# A continuation technique for maximum likelihood estimators in biological
models

Tyler Cassidy  School of Mathematics, University of Leeds, Leeds, LS2 9JT, UK
<EMAIL_ADDRESS>
## Abstract

Estimating model parameters is a crucial step in mathematical modelling and
typically involves minimizing the disagreement between model predictions and
experimental data. This calibration data can change throughout a study,
particularly if modelling is performed simultaneously with the calibration
experiments, or during an on-going public health crisis as in the case of the
COVID-19 pandemic. Consequently, the optimal parameter set, or maximal
likelihood estimator (MLE), is a function of the experimental data set. Here,
we develop a numerical technique to predict the evolution of the MLE as a
function of the experimental data. We show that, when considering
perturbations from an initial data set, our approach is significantly more
computationally efficient that re-fitting model parameters while resulting in
acceptable model fits to the updated data. We use the continuation technique
to develop an explicit functional relationship between fit model parameters
and experimental data that can be used to measure the sensitivity of the MLE
to experimental data. We then leverage this inverse sensitivity analysis to
select between model fits with similar information criteria, a priori
determine the experimental measurements to which the MLE is most sensitive,
and suggest additional experiment measurements that can resolve parameter
uncertainty.

##  1 Introduction

As quantitative modeling becomes more prevalent across biology and medicine
[Altrock et al.,, 2015; Perelson,, 2002; Sanche et al.,, 2020], mathematical
models are increasingly being developed during the experimental data
collection that will inform model parameters. This cooperation facilitates the
use of mathematical modelling to inform experimental design and suggest
potential intervention strategies [Zhang et al.,, 2022; Sanche et al.,, 2020;
Cárdenas et al.,, 2022; Luo et al.,, 2022]. The COVID-19 pandemic is a
striking example of the resulting feedback loop, where mathematical models
suggest intervention strategies that influence the evolving public health
crisis before being re-calibrated to new data. [Holmdahl and Buckee,, 2020;
Thompson,, 2020; Davies et al.,, 2020].

Each updated data set requires re-calibration of the model typically through
computationally expensive optimization techniques. To reduce this
computational cost of the re-calibration step, it is common to use the
existing parameters as a starting point when performing parameter fitting to
incoming experimental data sets. This approach recycles optimization work but
does not utilize leverage the relationship between the initial and updated
experimental data set. Here, we present a computational method to incorporate
information about evolving data sets during the model validation and parameter
estimation steps.

Specifically, for given model parameters and an initial experimental data set,
we develop a method to predict the best-fit parameter set to an updated
experimental data set. Our approach can be viewed as a numerical continuation
technique [Dhooge et al.,, 2008; De Souza and Humphries,, 2019]. However,
rather than studying the dynamical properties of the mathematical model as a
function of model parameters, we consider the evolution of best-fit model
parameters as a function of the experimental data. We use the necessary
condition for a local optima to write the best-fit parameters as an implicit
function of the experimental data. Thus, we predict best-fit parameter sets
for evolving experimental data without performing any optimization. Avoiding
optimization leads to significant computational savings and we demonstrate
these gains via two examples. In both these examples, our prediction method
produces comparable model fits to randomly perturbed data sets to optimization
techniques without the computational cost of solving the inverse optimization
problem.

While our approach does lead to increased computational efficiency, the more
immediate application of our work may be in experimental design. Specifically,
we identify an explicit relationship between individual best-fit parameter
values and individual experimental data points through our continuation
approach. We can therefore quantify which experimental measurements are the
most informative for determining best-fit parameters and measure the
sensitivity of parameter estimates to perturbations in data. The role of
experimental design in model selection and parameterization has been
extensively studied [Silk et al.,, 2014; Cárdenas et al.,, 2022; Li and Vu,,
2015, 2013]. In particular, Li and Vu, [2015] studied how correlations between
best-fit model parameters can impact practical and structural identifiability
of model parameters while Silk et al., [2014]; Cárdenas et al., [2022]
explored how experimental design impacts model selection from a class of
possible mathematical models. Conversely, our contribution explicitly relates
individual experimental measurements with individual best-fit parameter
estimates. We explicitly link our continuation technique to the Fisher
information matrix commonly used in optimal experimental design [Kreutz and
Timmer,, 2009; Braniff et al., 2019b, ]. Taken together, our approach allows
the increased confidence in model parametrization from optimal experimental
design to be mapped directly to individual model parameters. Accordingly, we
can therefore design experiments to address specific uncertainties in
parameter estimates.

Furthermore, our work offers a distinct step towards understanding how robust
parameter estimates are to evolving data. Many existing computational methods
quantify confidence in parameterization; formal parameter sensitivity analyses
[Marino et al.,, 2008; Maiwald et al.,, 2016; Zi,, 2011], virtual population
approaches [Allen et al.,, 2016; Cassidy and Craig,, 2019; Jenner et al.,,
2021], or parameter identifiability analysis [Castro and de Boer,, 2020],
often via profile likelihood computation [Raue et al.,, 2009, 2014; Kreutz et
al.,, 2012], quantify how robust model predictions are to parameter variation.
In particular, these techniques view the experimental data as fixed up to
experimental noise and focus on the relationship between model parameters and
model predictions. We offer a complementary approach to existing sensitivity
analysis by explicitly studying how the best-fit parameters vary due to
changes in calibration data. As we will see, our approach encodes information
from local sensitivity analysis when calculating the functional relationship
between the best-fit parameters and the calibration data. Consequently, while
classical sensitivity analysis quantifies variability in model output due to
change in model parameters, our approach considers changes in model
parameters, and thus model predictions, as a function of the calibration data.
We demonstrate this mapping of experimental data to best-fit parameter via an
example drawn from mathematical oncology [Cassidy et al.,, 2021]. These
results, when combined with existing information criteria like the AIC or BIC
[Kass and Raftery,, 1995], allow for modellers to quantify the robustness of
best-fit parameter estimates when comparing different model fits to
experimental data.

The remainder of the article is structured as follows. We begin by defining
the optimization problem in Section 2.1. We develop the continuation method in
Section 2.2, discuss our numerical implementation in 2.3, and explore the
connection between our continuation approach and classical profile likelihood
in 3.1. We then turn to two examples from mathematical biology to illustrate
the utility of our technique in Section 3.2 before finishing with a brief
discussion.

##  2 Methods

###  2.1 Formulation of the optimization problem

Here, we introduce the framework of the underlying optimization problem. We
focus on ordinary differential equation (ODE) models representing biological
processes, as these models are common throughout mathematical biology.
However, our approach extends to partial differential equation or delay
differential equation models directly. We consider a generic ODE based model
throughout the remainder of this work.

Let the model states be given by $x(t)\in\mathbb{R}^{n}$ with model parameters
denoted by $\theta\in\Omega\subset\mathbb{R}^{p}$ where $\Omega$ is a subset
of biologically plausible parameter values. We explicitly allow the initial
condition $x(0)$ to depend explicitly on the model parameters $\theta$. Taken
together, we consider the differential equation model

$\displaystyle\frac{\textrm{d}}{\textrm{dt}}x(t)=f(x,\theta);\quad
x(0)=x_{0}(\theta)$ (1)

where $f$ is continuously differentiable in $x$ and $\theta$.

We consider calibration data $\\{\phi_{i}\\}_{i=1}^{d\times m}$ representing
$m$ measurements each taken at $d$ time points $\\{t_{i}\\}_{i=1}^{d}$. It is
possible that model species are not directly comparable against the
calibration data so we define the $m$ model observables by

$\displaystyle y_{i}(\theta)=h(x(t_{i},\theta),\theta)\in\mathbb{R}^{d\times
m}.$

In what follows, we consider $m=1$ for notational simplicity although the
analysis extends for $m\geqslant 2$.

#### Likelihood function and objective function

######  Remark 2.1.

The methods that follow do not assume a specific objective function. However,
we do assume that the objective function is twice continuously differentiable
as is commonly the case. For simplicity, we present the remainder of our
results using the common log-likelihood formulation [Stapor et al.,, 2018;
Maiwald et al.,, 2016].

The likelihood describes the probability of observing experimental data $\phi$
as a function of $\theta$ and is given by

$\displaystyle\mathcal{L}(y(x(t,\theta)),\phi)=\prod_{i=1}^{d}\frac{1}{\sqrt{2\pi\sigma^{2}_{i}}}\exp\left[-\frac{(y_{i}(\theta)-\phi^{*}_{i})^{2}}{\sigma_{i}^{2}}\right]$
(2)

The experimental error at each measurement point, $\sigma_{i}$, can be
estimated as an additional model parameter or fixed to a known value. Here, we
follow Sharp et al., [2022] and take $\sigma_{i}$ fixed at a known constant
value, although it is possible to include $\sigma_{i}$ in the vector of
unknown parameters $\theta$. The maximum likelihood estimator (MLE)
$\theta^{*}$, and thus best-fit model parameters for the given experimental
data $\phi$, is defined by the solution of the inverse problem

$\displaystyle\theta^{*}=\textrm{argmax}_{\theta\in\Omega}\mathcal{L}(\theta,\phi^{*}).$

As the differential equations defining $y(x(t,\theta))$ rarely have explicit
solutions, the likelihood (2) is difficult to evaluate analytically. It is
therefore standard to minimize the negative log-likelihood
$G(\theta,\phi)=-\log\left(\mathcal{L}(y(x(t,\theta)),\phi^{*})\right)$ given
by

$\displaystyle
G(\theta,\phi)=\displaystyle\sum_{i=1}^{d}\log\left(\sqrt{2\pi\sigma_{i}^{2}}\right)+\frac{(y_{i}(\theta)-\phi_{i}^{*})^{2}}{\sigma_{i}^{2}}.$
(3)

Under the assumption that $\sigma_{i}=\sigma$ is fixed, the error term
$\log\left(\sqrt{2\pi\sigma^{2}}\right)$ and denominator of $G(\theta,\phi)$
are constant and do not influence the solution of the optimization problem.
The maximum likelihood estimator $\theta^{*}$ is the parameter set that
minimizes $G(\theta,\phi^{*})$. A number of computational techniques exist to
minimize $G(\theta,\phi)$ and thus calculate $\theta^{*}$. These optimization
techniques typically require simulating the mathematical model (1) at each
optimization step. Further complicating the optimization, $G(\theta,\phi^{*})$
is often non-convex with multiple local minima.

###  2.2 Continuation of maximal likelihood estimator

In (3), we explicitly write the objective function $G$ as a function of the
model parameters $\theta$ and the experimental data $\phi$. Accordingly, the
MLE $\theta^{*}$ is an implicit function of the experimental data $\phi$
defined as the solution of the optimization problem

$\displaystyle\theta^{*}(\phi)=\textrm{argmax}_{\theta\in\Omega}\mathcal{L}(\theta,\phi).$
(4)

Model fitting is increasingly performed concurrently with experiments [Luo et
al.,, 2022] or obtained from an evolving real-world scenario, as in epidemic
modelling [Sanche et al.,, 2020]. In both of these cases, the experimental
data is evolving and should not be considered as known and constant.
Accordingly, we are interested in the MLE as a function of the experimental
data $\phi$. Most existing optimization techniques consider the experimental
data fixed and omit this dependence. Here, we develop a continuation type
technique to compute the evolution of $\theta^{*}$ numerically as a function
of $\phi$ from an initial solution of the optimization problem. Ultimately, we
calculate the evolution of $\theta^{*}(\phi)$ as the calibration data varies
to generate a curve of potential MLEs in $(\phi,\theta^{*})$ space using a
numerical continuation technique.

Numerical continuation methods compute branches of implicitly defined curves.
A standard application of these continuation type techniques in mathematical
biology is numerical bifurcation analysis [Dhooge et al.,, 2008; Sanche et
al.,, 2022]. In their most common form, numerical bifurcation techniques
compute equilibrium systems of a non-linear dynamical system as a function of
model parameters but can be used to detect much richer dynamical behaviour [De
Souza and Humphries,, 2019]. Often, these continuation techniques leverage
“predictor-corrector” algorithms. Predictor-corrector approaches use the
implicit function theorem to predict the solution to the corresponding non-
linear system of equations. Then, the predicted solution is used as a starting
value to explicitly calculate the solution of the system of equations during
the corrector step. Here, we develop a similar “prediction-correction”
strategy to predict the behaviour of the solution $\theta^{*}(\phi)$ of the
inverse problem (4) as a function of the data $\phi$. We focus on the
“predictor” step, as the corrector step, if necessary, can utilize existing
numerical optimization techniques to calculate the MLE.

As the log-likelihood (3) is continuously differentiable, local optimal must
satisfy

$\displaystyle\mathrm{D}_{\theta}G(\theta^{*},\phi)=0,$ (5)

so we necessarily have

$\displaystyle\theta^{*}(\phi)\in\\{\theta\in\Omega|\mathrm{D}_{\theta}G(\theta^{*},\phi)=0\\}.$

However, unlike the implicit equation used to determine equilibria of a
dynamical system and used in continuation techniques for numerical bifurcation
analysis, the optimality condition (5) is a necessary, but not sufficient,
condition for $\theta^{*}$ to be a MLE. Models that are not structurally
identifiable [Raue et al.,, 2014] have manifolds in parameter space on which
this optimality constraint holds but are not necessarily MLEs. We discuss the
relationship between our approach and profile likelihood classifications of
structural identifiability in Section 3.1.

Now, let $\theta^{*}_{0}$ be the MLE for calibration data $\phi_{0}$. Further,
let the Hessian $\mathrm{D}^{2}_{\theta}G(\theta,\phi)$ be invertible at
$(\theta_{0}^{*},\phi_{0})\in\mathbb{R}^{p}\times\mathbb{R}^{d}$ and consider
the function

$\displaystyle\mathrm{D}_{\theta}G(\theta^{*},\phi):\mathbb{R}^{p}\times\mathbb{R}^{d}\to\mathbb{R}^{p}.$

Then, the implicit function theorem ensures the existence of a function
$\Psi(\phi)$ such that

$\displaystyle\mathrm{D}_{\theta}G(\Psi(\phi),\phi)=0$

in a neighbourhood of $\phi_{0}$ with $\Psi(\phi_{0})=\theta^{*}(\phi_{0})$.
It is natural to consider $\Psi(\phi)$ as the predicted MLE $\theta^{*}(\phi)$
for $\phi$ in a neighbourhood of $\phi_{0}$.

The implicit function theorem ensures that $\Psi$ exists but computing
$\Psi(\phi)$ analytically is functionally impossible. However, the implicit
function $\Psi(\phi)$ is continuously differentiable and we expand $\Psi$ as a
function of the calibration data $\phi$ using Taylor series

$\displaystyle\Psi(\phi+\Delta\phi)=\Psi(\phi)+\mathrm{D}\Psi(\phi)\Delta\phi+\mathcal{O}(\Delta\phi^{2}).$
(6)

where $\phi+\Delta\phi$ is the updated calibration data. Then, to predict
$\Psi$ starting from a known solution $\Psi(\phi)=\theta^{*}$ we calculate
$\mathrm{D}\Psi(\phi)$. The implicit function theorem implies that

$\displaystyle\mathrm{D}\Psi=-\left[\mathrm{D}_{\theta}^{2}G(\Psi(\phi),\phi)\right]^{-1}\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi).$

We thus use $\mathrm{D}\Psi$ to evaluate (6) and thus perform the continuation
step.

###  2.3 Numerical Implementation

We now show how to use the objective function (3) to calculate finite
difference approximations to the derivatives included in (6). As before, we
assume that we are given a point
$(\theta_{0}^{*},\phi_{0})\in\mathbb{R}^{p}\times\mathbb{R}^{d}$ such that

$\displaystyle\theta_{0}^{*}=\textrm{argmin}_{\theta\in\Omega}G(\theta,\phi_{0}).$

For $\theta_{n}$ denoting the $n$-th parameter, we calculate

$\displaystyle\frac{\partial
G(\theta,\phi)}{\partial\theta_{n}}=\displaystyle\sum_{i=1}^{d}2\left(y_{i}(\theta)-\phi_{i}\right)\frac{\partial
y_{i}(\theta)}{\partial\theta_{n}}$

and so

$\displaystyle\left[\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)\right]_{(n,i)}=-2\frac{\partial
y_{i}(\theta)}{\partial\theta_{n}}.$ (7)

The derivatives $\partial_{\theta_{n}}y_{i}(\theta)$ can be calculated through
finite difference schemes [Zi,, 2011]

$\displaystyle\frac{\partial
y_{i}(\theta)}{\partial\theta_{n}}=\frac{y_{i}(\theta+\Delta\theta_{n})-y_{i}(\theta-\Delta\theta_{n})}{2\Delta\theta_{n}}+\mathcal{O}\left((\Delta\theta_{n})^{2}\right),$

where $\Delta\theta_{n}$ is a small perturbation in only the $n$-th parameter.
In practice, it is standard to take $\Delta\theta_{n}$ to be some small
percentage of the initial parameter $\theta_{n}$ [Li et al.,, 2011]. In this
case, computing $\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$ requires $2p$
model simulations where $p$ is the number of model parameters. We note that
$\partial_{\theta_{n}}y_{i}(\theta)$ is commonly used to perform local
sensitivity analysis and that more accurate finite difference approximations,
such as centered differences, can be used to calculate
$\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$.

Calculating the Hessian $\mathrm{D}_{\theta}^{2}G(\theta,\phi)$ via finite
differences is simple to implement but computationally expensive due to the
number of objective function evaluations. However, the Hessian, or the
observed Fisher Information, is commonly used throughout parameter
optimization algorithms and other techniques such as profile likelihood
calculations, estimates of the likelihood function, and classical sensitivity
anaylsis, which has led to recent advances in the development of
computationally efficient techniques to calculate
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$ [Stapor et al.,, 2018] and the ability
to recycle these calculations to avoid computational cost.

In the following examples, we use a finite difference scheme to calculate
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$. We calculate the diagonal elements of
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$ using forward second order differences
and the off-diagonal terms by

$\displaystyle\frac{\partial
G(\theta,\phi)}{\partial\theta_{i}\partial\theta_{j}}$
$\displaystyle=\left(\frac{1}{4(\Delta\theta_{i})(\Delta\theta_{j})}\right)\left[G(\theta+\Delta\theta_{i}+\Delta\theta_{j},\phi)-G(\theta+\Delta\theta_{i}-\Delta\theta_{j},\phi)\right.$
$\displaystyle{}\quad\left.+G(\theta-\Delta\theta_{i}+\Delta\theta_{j},\phi)+G(\theta-\Delta\theta_{i}-\Delta\theta_{j},\phi)\right]+\mathcal{O}\left((\Delta\theta_{i})^{2},(\Delta\theta_{j})^{2}\right).$

Thus, our computation of the Hessian requires $2p(p+1)$ objective function
evaluations, although, as mentioned, more efficient implementations are
available. In fact, many gradient-based optimization techniques approximate
the Hessian $D^{2}_{\theta,\theta}G(\theta,\phi)$ at each iteration [MATLAB,,
2017]. For example, both fmincon and fminunc in [MATLAB,, 2017] calculate
$D^{2}_{\theta,\theta}G(\theta,\phi)$ at each step and print the pre-computed
Hessian as an output of the optimizer. It is therefore possible, and
efficient, to recycle this calculation when calculating an update to
$\theta_{0}^{*}$ using (5).

All told, this numerical implementation requires $2p(p+2)$ model simulations
to evaluate (5). This computational cost is certainly not optimal but does
benefit from re-using calculations performed in local sensitivity analysis and
the optimization step. Finally, while we have written (5) with the inverse of
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$, it is computationally more
appropriate to solve the linear system of equations

$\displaystyle\mathrm{D}_{\theta}^{2}G(\theta,\phi)\mathrm{D}\Psi=-\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$

for the unknown $\mathrm{D}\Psi$.

Code to implement this continuation technique is available at
https://github.com/ttcassid/MLE$\\_$Continuation.

##  3 Results

###  3.1 Relationship with existing techniques

There are a number of existing techniques to study the relationship between
model parameters and data. While our continuation technique focuses on the
relationship between the MLE and the calibration data, it has many ties to
these existing techniques. We therefore discuss how this continuation method
relates to parameter identifiability as assessed by the profile likelihood;
local sensitivity analysis; and experimental design, with a focus on using the
explicit relationship between data and the MLE to suggest additional
experimental measurements.

#### Parameter identifiability

Thus far, we have explicitly written the MLE estimator as a function of the
experimental data used to fit a model. Our approach is intrinsically related
to parameter identifability analysis. Identifiability analysis attempts to
determine if available experimental observations are capable to uniquely
determine model parameters. Accordingly, the practical identifiability of a
mathematical model depends on available experimental data. The profile
likelihood, given by

$\displaystyle
PLE_{\theta_{i}}(c)=\min_{\theta_{i}=c,\theta\in\mathbb{R}^{p}}G(\theta,\phi),$

and introduced by Raue et al., [2009], is a projection of the likelihood
function onto the model parameter $\theta_{i}=c$. The profile likelihood
illustrates the behaviour of the likelihood function as the parameter
$\theta_{i}$ is fixed away from the optimal value $\theta^{*}_{i}$. The shape
of $PLE_{\theta_{i}}(c)$ illustrates the confidence interval of the parameter
estimate $\theta_{i}^{*}$ for given experimental data. Formally, Raue et al.,
[2009] define these confidence intervals by

$\displaystyle\textrm{C.I.}(\theta_{i},\alpha)=\\{c|PLE_{\theta_{i}}(c)-PLE_{\theta_{i}}(\theta_{i}^{*})<\Delta_{\alpha}\\}$

where $\Delta_{\alpha}=\chi^{2}(\alpha,df)$ is the $\chi^{2}$ distribution at
significance level $\alpha$ and $df$ degrees of freedom [Raue et al.,, 2009].
A parameter is practically identifiable in the sense of Raue et al., [2009]
with confidence level $\alpha$ if $\textrm{C.I.}(\theta_{i},\alpha)$ is
bounded in parameter space for given experimental data. Conversely, a non-
identifiable parameter has a profile likelihood that does not increase past
the threshold $\Delta_{\alpha}$.

The profile likelihood is intrinsically linked to the available experimental
data $\phi_{i}$. We view the PLE as a function of both the parameter
$\theta_{i}$ and the experimental data $\phi$

$\displaystyle
PLE_{\theta_{i}}(c,\phi)=\min_{\theta_{i}=c,\theta\in\mathbb{R}^{p}}G(\theta,\phi).$

For practically unidentifiable models, it is natural to ask what perturbations
to the experimental data could render the model practically identifiable. Raue
et al., [2009] use the profile likelihood of a model parameter to suggest
additional experiments to resolve practical non-identifiability. They simulate
the model for parameter values along $PLE_{\theta_{i}}$ to suggest additional
experimental measurements at times $t_{s,i}$, where $t_{s,i}$ represents the
$i-$th simulated measurement time. In our framework, we define

$\displaystyle\theta^{*}|_{\theta_{i}=c}(\phi)=\textrm{argmin}_{\theta_{i}=c,\theta\in\mathbb{R}^{p}}G(\theta,\phi),$

so that

$\displaystyle
PLE_{\theta_{i}}(c,\phi)=G(\theta^{*}|_{\theta_{i}=c}(\phi),\phi).$

We note that the definition of $\theta^{*}|_{\theta_{i}=c}(\phi)$ is precisely
that of $\theta^{*}(\phi)$ with the added constraint that $\theta_{i}=c$. We
can calculate $\mathrm{D}_{\phi}\theta^{*}|_{\theta_{i}=c}$ as a function of
the experimental data $\phi$ in precisely the same manner as described
previously. Consequently, our continuation approach can complement the
experimental design approach suggested by Raue et al., [2009] by incorporating
the sensitivity of the MLE to perturbations in the (simulated or experimental)
calibration data.

#### Sensitivity analysis

Local sensitivity analysis quantifies how small perturbations of the best-fit
parameters impact model output [Zi,, 2011]. A standard approach to local
sensitivity analysis is using the finite difference approximation of

$\displaystyle S_{n}(t)=\frac{\partial
y(\theta)}{\partial\theta_{n}}=\frac{h(t_{i},\theta+\Delta\theta_{n})-h(t_{i},\theta-\Delta\theta_{n})}{\Delta\theta_{n}}+\mathcal{O}\left(\Delta\theta_{n}\right)$

to identify which parameter values strongly impact model projections. When
$|S_{n}|$ is small, the model output is considered to be insensitive to
$\theta_{n}$. The $n$-th row of
$\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$ is precisely $S_{n}(t_{i})$
for $t_{i}$ corresponding to calibration data measurements. When implementing
(5), the magnitude of the continuation step $\mathrm{D}\Psi(\phi)\Delta\phi$
in the direction of $\theta_{n}$ is scaled by $S_{n}$. This scaling encodes
the local sensitivity of model predictions to variations in parameters in the
prediction of $\Psi(\phi)$. Consequently, our continuation method naturally
includes the information gained from local sensitivity analysis.

#### Experimental design

In our derivation of $\mathrm{D}\Psi$, we assumed that the Hessian matrix
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$ was invertible. The Hessian gives the
curvature of the loglikelihood and is known as the observed Fisher information
matrix $\mathcal{I}_{obs}$. The observed Fisher information is a local
measurement in data space. Conversely, the expected Fisher information
considers the entirety of data space for fixed model parameters $\theta$. The
expected Fisher information is obtained by taking the expectation of
$\mathrm{D}_{\theta}^{2}G(\theta,\phi)$ over all possible experimental
measurements $\phi$ and is defined via

$\displaystyle\mathcal{I}=\mathbb{E}\left[\mathrm{D}_{\theta}^{2}G(\theta,\phi)\right].$

Many existing experimental design methods leverage the expected Fisher
information matrix to minimize the covariance in model parameter estimates via
the Cramér-Rao inequality. These experimental design techniques typically
maximize some aspect, often the determinant, of the Fisher information matrix
as a function of possible data to select the most informative calibration data
set [Kreutz and Timmer,, 2009]. From a geometric perspective, maximizing the
determinant of the Fisher information matrix corresponds to minimizing the
volume of the confidence ellipsoid engendered from the covariance matrix
[Braniff et al., 2019b, ].

In particular, Braniff et al., 2019a  considered the case of bistable gene
regulatory networks where the fold bifurcation and unstable manifold between
stable equilibria complicates experimental design and parameter estimation.
Sharp et al., [2022] considered an information-geometry perspective to propose
the expected Fisher information matrix and resulting Riemannian manifold as a
guide for data collection. As is often the case, both Sharp et al., [2022] and
Braniff et al., 2019a  used the expected Fisher information, which considers
all possible calibration data via the expectation over $\phi$. Here, we show
how our approach complements the classical Fisher information approach to
experimental design, albeit through a local measurement, in $(\theta,\phi)$
space. We recall that

$\displaystyle\mathrm{D}\Psi\Delta\phi=-\left[\mathcal{I}_{obs}\right]^{-1}\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)\Delta\phi,$

so if $\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$ were the identity, then
$\mathrm{D}\Psi$ would correspond to the Fisher information approach to
measuring uncertainty in MLE.

In the calculation of $\mathrm{D}\Psi\Delta\phi$, the matrix
$\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$ maps perturbations in the
calibration data $\Delta\phi$ through the curvature of the loglikelihood to
changes in the MLE. Consequently,
$\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)$ acts as a change of basis
matrix from the space of calibration data to parameter space. Simply,
$\mathrm{D}_{\theta,\phi}^{2}G(\Psi(\phi),\phi)\Delta\phi$ scales changes in
the calibration data to the confidence ellipsoid in parameter space obtained
from $\left[\mathcal{I}_{obs}\right]^{-1}$. Geometrically, if
$\mathrm{D}^{2}_{\theta}G$ has eigenvalues $\lambda_{i}$ with corresponding
eigenvectors $\nu_{i}$, then choosing $\Delta\phi$ such that
$\nu_{i}=\mathrm{D}^{2}_{\theta}G\Delta\phi$ translates perturbations in
calibration data to the corresponding eigenspace of the covariance matrix.

For example, the $i-$th column of $\mathrm{D}\Psi$ maps perturbations of the
$i-$th data point to changes in the MLE. Specifically, the sum

$\displaystyle\frac{\Delta\theta^{*}}{\Delta\phi_{k}}=\displaystyle\sum_{k=1}^{p}|\mathrm{D}\Psi_{k,j}|$

measures the sensitivity of the MLE $\theta^{*}$ to perturbations in the
$k-$th data point. Thus,

$\displaystyle\|\mathrm{D}\Psi\|_{1}=\displaystyle\max_{k=1,2,...,p}\frac{\Delta\theta^{*}}{\Delta\phi_{k}}$

and the most informative data point satisfies

$\displaystyle
l=\displaystyle\textrm{argmax}_{k=1,2,...,p}\displaystyle\frac{\Delta\theta^{*}}{\Delta\phi_{k}},$

where informative is understood as the data point inducing the largest
sensitivity in the MLE. As an extreme example, if

$\displaystyle\frac{\Delta\theta^{*}}{\Delta\phi_{n}}=0,$

then perturbations in $\phi_{n}$ do not impact the MLE estimate, which implies
complete insensitivity of the model fit to $\phi_{k}$. This example
corresponds to $\Delta\phi$ belonging to the kernel of the matrix
$\mathrm{D}^{2}_{\theta,\phi}G$ since we have assumed that
$\mathrm{D}^{2}_{\theta}G$ is invertible.

We can therefore utilize our analysis to identify which additional
experimental measurements could increase confidence in model parameterization.
Consider $k$ additional measurements
$\\{\phi_{s,i}=y_{s,i}(\theta^{*})\\}_{i=1}^{k}$ taken directly from the model
simulation at times $\\{t_{s,i}\\}_{i=1}^{k}$ where the subscript $s$
indicates simulated data. Including $\\{\phi_{s,i}\\}$ in the objective
function (3) does not change the MLE or objective value function as these
simulated data exactly match the model values. However,
$\|\mathrm{D}\Psi(\phi+\Delta\phi_{s,i})\|$ quantifies the sensitivity of the
MLE to variability in the $k$ simulated measurements. Accordingly, the
measurement that maximizes $\|\mathrm{D}\Psi(\phi+\Delta\phi_{s,i})\|$ for a
fixed perturbation size $\Delta$ is a good candidate for an additional
experimental measurement to decrease parameter uncertainty.

###  3.2 Examples

The continuation framework derived earlier is applicable to a large variety of
models throughout in the mathematical biology literature. To demonstrate the
utility of the continuation method, we consider two examples from distinct
fields and model formulations. First, we consider a mathematical model of
phenotypic heterogeneity in non-small cell lung cancer (NSCLC) [Cassidy et
al.,, 2021]. This model is given by a system of two non-local, structured PDEs
representing the density of drug-sensitive and drug-tolerant NSCLC cells. The
PDE model is equivalent to a system of integral equations following the
introduction of two auxiliary variables which can be further reduced to a
system of ODEs (see [Cassidy et al.,, 2021] for details). The parameters of
the ODE model were fit to in vitro NSCLC data taken from growth experiments in
treated and untreated media [Cassidy et al.,, 2021].

We also consider a classical model of HIV-1 viral dynamics. This model has
been used extensively to understand viral dynamics data [Perelson,, 2002] and
the identifiability of model parameters was considered by Wu et al., [2008].
In that work, Wu et al., [2008] used simulated data to validate their
identifiability results; we follow Wu et al., [2008] and use simulated data to
illustrate our approach.

#### A PDE model of phenotypic switching in mathematical oncology

Non-genetic phenotypic heterogeneity has been increasingly studied as a driver
of treatment resistance in solid cancers [Goldman et al.,, 2015]. A number of
mathematical models have been derived to study the emergence of phenotypic
plasticity in cancer cell lines [Gunnarsson et al.,, 2020; Jolly et al.,,
2018; Sahoo et al.,, 2021; Craig et al.,, 2019]. We consider the Cassidy et
al., [2021] model that tracks the density of NSCLC cells with a drug-sensitive
($A(t,a)$) or drug-tolerant ($B(t,a)$) phenotype at time $t$ and age $a$. The
total number of cells of each phenotype is given by

$\bar{A}(t)=\int_{0}^{\infty}A(t,a)\mathrm{d}a\quad\textrm{and}\quad\bar{B}(t)=\int_{0}^{\infty}B(t,a)\mathrm{d}a.$
(8)

The total number of NSCLC cells is given by $N(t)=\bar{A}(t)+\bar{B}(t)$.
Cassidy et al., [2021] considered logistic growth with an Allee effect,
wherein cooperation between cells of the same phenotype can lead to increased
growth rates, given by

$\displaystyle R_{A}(\bar{A}(t),\bar{B}(t))$
$\displaystyle=r_{A}\left(1-\frac{\bar{A}(t)+\bar{B}(t)}{K}\right)\quad\textrm{and}$
$\displaystyle R_{B}(\bar{A}(t),\bar{B}(t))$
$\displaystyle=r_{B}\left(1-\frac{\bar{A}(t)+\bar{B}(t)}{K}\right)f_{n}(\bar{A}(t),\bar{B}(t)).$
(9)

where $r_{A}$ and $r_{B}$ are phenotype specific growth rates, the carrying
capacity is $K$, and the strength of the Allee effect is

$\displaystyle
f_{n}(\bar{A}(t),\bar{B}(t))=1+\left(\frac{r_{A}-r_{B}}{r_{B}}\right)\left(\frac{\bar{B}(t)^{n}}{\bar{A}(t)^{n}+\bar{B}(t)^{n}}\right).$

Finally, drug-tolerant and drug-sensitive cells have phenotype-specific death
rates $d_{B}$ and

$d_{A}=\left\\{\begin{array}[]{cc}d_{A}&\textrm{If untreated}\\\
d_{A}^{max}&\textrm{During treatment.}\end{array}\right.$

$A(t,a)$ and $B(t,a)$ satisfy the age structured PDEs

$\left.\begin{aligned}
\partial_{t}A(t,a)+\partial_{a}A(t,a)&=-[d_{A}+R_{A}(\bar{A}(t),\bar{B}(t))]A(t,a)\\\
\partial_{t}B(t,a)+\partial_{a}B(t,a)&=-[d_{B}+R_{B}(\bar{A}(t),\bar{B}(t))]B(t,a)\end{aligned}\right\\}$
(10)

with boundary conditions corresponding to cellular reproduction given by

$\left.\begin{aligned}
A(t,0)&=2\int_{0}^{\infty}\left[R_{A}(\bar{A}(t),\bar{B}(t))\beta_{AA}(a)A(t,a)+f_{n}(\bar{A}(t),\bar{B}(t))R_{B}(\bar{A}(t),\bar{B}(t))\beta_{BA}(a)B(t,a)\right]\mathrm{d}a\\\
B(t,0)&=2\int_{0}^{\infty}\left[R_{A}(\bar{A}(t),\bar{B}(t))\beta_{AB}(a)A(t,a)+f_{n}(\bar{A}(t),\bar{B}(t))R_{B}(\bar{A}(t),\bar{B}(t))\beta_{BB}(a)B(t,a)\right]\mathrm{d}a.\end{aligned}\right\\}$
(11)

The functions $\beta_{ij}$ represent the probability of a reproducing mother
cell with age $a$ and phenotype $i$ giving birth to a daughter cell with
phenotype $j$. The probability of phenotypic inheritance is given by

$\beta_{ii}(a)=P_{ii}^{*}+(P_{ii}^{max}-P_{ii}^{*})\exp\left(-\sigma_{i}a\right),$

where $\sigma_{i}$ represents the decay rate of intracellular signalling
factors that modulate how ageing impacts the probability of daughter cells
retaining the mother cells phenotype, and

$\beta_{AB}(a)=1-\beta_{AA}(a)\quad\textrm{and}\quad\beta_{BA}(a)=1-\beta_{BB}(a).$

Further details, including a derivation of the initial conditions of (10),
model analysis, and reduction of the phenotype switching mode (10) to a system
of ODEs can be found in Cassidy et al., [2021].

The model (10) was fit to in vitro experimental data corresponding to NSCLC
cell population growth in untreated and treated environments where treatment
is applied from day $3$ onwards. The calibration data is 4 data points
$\\{\phi_{i}\\}_{i=1}^{4}$ collected at time $t_{i}=0,2,4,6$ days in the
control experiment, and two additional data points $\\{\phi_{i}\\}_{i=5}^{6}$
collected on days $t_{i}=4,6$ days during the treated experiment. As anti-
cancer treatment is applied from day 3 on-wards and decreases the cancer cell
population, we necessarily have $\phi_{5}\leqslant\phi_{3}$ and
$\phi_{6}\leqslant\phi_{4}$. We denote the experimental data used to
parametrize the model by $\\{\phi_{i}^{0}\\}_{i=1}^{6}$. The model output
corresponding to the experimental measurements is thus

$\displaystyle y_{i}(\theta)=N(t_{i},\theta),$

and the objective function is the standard sum of squares error given by

$\displaystyle
G_{pheno}(\theta,\phi)=\sqrt{\displaystyle\sum_{i=1}^{6}\left(\log_{10}(N(t_{i},\theta)-\log_{10}(\phi_{i})\right)^{2}}.$

Cassidy et al., [2021] fit model parameters
$[r_{A},r_{B},d_{A}=d_{B},d_{A}^{max}]$ to treated and untreated experimental
data simultaneously for a number of cell lines. The MLE found by Cassidy et
al., [2021] corresponds to
$\theta^{*}(\phi^{0})=[0.4827,0.3498,0.7025,0.4198]$.

We perturbed the experimental data collected by Craig et al., [2019] with
increasing amounts of Gaussian noise. We created 10 perturbed data sets
$\\{\phi_{i}^{j}\\}_{i=1}^{6}$ where the index $j=1,2,...,10,$ denotes the
$j$-th perturbed data set and the normally distributed noise with $\mu=0$,
$\sigma^{2}=1$, and scaled such that

$\displaystyle\|\log_{10}(\phi_{i}^{j})-\log_{10}(\phi_{i}^{*})\|=(0.05+jh_{step}\times\left(0.75-2\times
0.05\right)\left(\frac{2}{10(11)}\right)\|\log_{10}(\phi_{i}^{0})\|$

where $h_{step}=0.65/55$ was chosen such that
$\|\log_{10}(\phi_{i}^{10})-\log_{10}(\phi_{i}^{0})\|=0.75\|\log_{10}(\phi_{i}^{0})\|.$

We enforce that this randomly perturbed data satisfies
$\phi_{5}\leqslant\phi_{3}$ and $\phi_{6}\leqslant\phi_{4}$. For each
perturbed data set $\\{\phi_{i}^{j}\\}$, we used the continuation method
described in Section 2.2 to calculate

$\displaystyle\Psi(\phi^{j})=\theta^{*}(\phi^{j-1})+\mathrm{D}\Psi(\phi^{j-1})\Delta\phi+\mathcal{O}(\Delta\phi^{2}).$
(12)

The naive approach to calculate the MLE $\theta^{*}(\phi^{j})$ for updated
data $\phi^{j}$ would be to use the MLE from the previous data,
$\theta^{*}(\phi^{j-1})$, as an initial starting guess for the parameter
fitting step. Hence, to illustrate the utility of our continuation technique,
we calculated $\Psi(\phi^{j})$ using (12) and then calculated
$G_{pheno}(\Psi(\phi^{j}),\phi^{j})$. We also calculated the true MLE
$\theta^{*}(\phi^{j})$ using the Matlab algorithm fmincon from the starting
guesses $\Psi(\phi^{j})$ and $\theta^{*}(\phi^{j-1})$. In Figure 1 A), we show
the objective function value evaluated at the updated data $\phi^{j}$ and
three parameter sets : the naive starting point, $\theta^{*}(\phi^{j-1})$; the
predicted MLE, $\Psi(\phi^{j})$; and the true MLE, $\theta^{*}(\phi^{j})$. We
note that the non-monotonic profile of the objective function $G_{pheno}$ in
Figure 1 A) is to be expected as we are adding noise to experimental data.
This noise may perturb the existing data away from dynamics that can be well-
described by the mathematical model. Accordingly, the important information
from Figure 1 A) is the comparison

$\displaystyle G_{pheno}(\theta^{*}(\phi^{i}),\phi^{i})\leqslant
G_{pheno}(\Psi(\phi^{i}),\phi^{i})<G_{pheno}(\theta^{*}(\phi^{i-1}),\phi^{i}),$

which demonstrates the accuracy of the continuation step (5) in driving a
relative decrease in $G_{pheno}$.

Further, in Figure 1 B), we show the cumulative number of objective function
evaluations when calculating $\theta^{*}(\phi^{j})$ for $j=1,2,...,10$ when
starting the optimization from $\theta^{*}(\phi^{j-1})$ and $\Psi(\phi^{j})$.
The total number of function evaluations used is lower when starting the
optimization from the predicted MLE $\Psi(\phi^{j})$ than when starting from
$\theta^{*}(\phi^{j-1}$. More strikingly, the predicted MLE
$G(\Psi(\phi^{j}),\phi^{j})$ is comparable against
$G(\theta^{*}(\phi^{j},\phi^{j})$ in Figure 1 A) and there is computational
benefit to only calculating the predicted MLE $\Psi(\phi^{j})$ rather than re-
fitting the parameters. Taken together, the results shown in Figure 1
demonstrate the accuracy and computation efficiency gained by calculating
$\Psi(\phi^{j})$.

Figure 1:  Comparison between MLE estimates obtained using the naive and
continuation approaches. Panel A shows a comparison of the objective function
value for the naive and continuation guesses as well as the true minimal
objective function value as a function of the perturbation of the experimental
data from the initial data. Panel B shows a comparison of the number of
objective value evaluations required to obtain the minimal value when starting
from the naive or predicted MLE with the number of function evaluations
required to calculate $\Psi(\theta^{i})$.

We now demonstrate how to utilize the continuation frame work to identify
additional time points to increase confidence in model parameters. We focus on
the treated environment and consider additional time points
$t_{s,i}=3.1,3.2,3.3,3.4,3.5,5,7$ days with corresponding simulated
measurements $\\{\phi_{i,s}\\}_{i=1}^{7}=N(t_{s,i}).$ We perturb each of these
simulated measurements by a fixed amount, $\Delta\phi=\pm 0.3N(3.1)$, to give
14 additional, perturbed measurements. We appended each of these 14
measurements to the experimental data and predicted the MLE to these appended
data sets.

We calculated the relative change in the MLE for each model parameter and each
of the 14 appended data sets. We note that each of the simulated data point
occurs following the beginning of therapy. The immediate decrease observed in
$N(t)$ following the beginning of treatment is due to the death of sensitive
cells following treatment administration and controlled by the parameter
$d_{a}^{max}$. From the biological interpretation of the parameters, we expect
$d_{a}^{max}$ to be highly sensitive to perturbations in these data points.

As expected, $d_{a}^{max}$ was the most sensitive model parameter to
perturbations of the simulated data and we show the percent relative change in
$d_{a}^{max}$ from the unperturbed data in Figure 2 B). As expected, the
maximal death rate of sensitive cells increased when the simulated data point
was decreased from the true value and decreased when the simulated data point
was increased.

The treatment sensitive population rapidly shrinks during therapy. The
stabilization and rebound of the population during therapy is due to the
expansion of the drug resistant population. This stabilization occurs once the
drug sensitive population has been maximally suppressed which due to the drug
effect. The most informative simulated data point, as measured by the
magnitude of the relative change in the parameter $d_{a}^{max}$, was at time
$t_{i,s}=3.4$. At $t=3.4$, drug sensitive cells are no longer dominant due to
drug pressure. The depth of the population response to treatment, as measured
by $N(3.4),$ is thus highly sensitive to death rate of these drug sensitive
cells under treatment. In Figure 2 A), we show the simulated experimental
measurements and predicted model dynamics for the most informative time point.
The predicted model simulations capture the perturbed data point while
retaining good fits to the true experimental data.

Figure 2: Evaluating additional time points to identify $d_{a}^{max}$ in an in
vitro model of NSCLC. Panel A shows the a selection of predicted model
dynamics when fit to experimental data with a single additional time point
$\phi_{i,s}^{*}$ that is perturbed by a $\Delta\phi$ from the true simulated
value. For figure clarity, model trajectories corresponding to the
perturbation of $\\{\phi_{4,s}\\}$ is shown. Panel B shows a tornado plot of
the predicted relative change in the best-fit parameter $d_{a}^{max}$ for each
additional simulated data point $\\{\phi_{i,s}\\}_{i=1}^{7}$. The left side of
the tornado plot, in blue, shows the relative change when the perturbed value
$\phi_{i,s}=\phi_{i,s}^{*}+\Delta\phi$ is larger than the simulated value
$\phi_{i,s}^{*}$. The right-hand side, in orange, shows the relative change in
$d_{a}^{max}$ when $\phi_{i,s}=\phi_{i,s}^{*}+\Delta\phi$ is smaller than the
simulated value $\phi_{i,s}^{*}$.

#### Parameter continuation in a viral dynamics model

The standard viral dynamics model has been extensively used to understand the
dynamics of viral infection in HIV-1 [Perelson,, 2002]. The model tracks the
concentration of uninfected target cells, $T(t)$, infected cells $I(t)$, and
free infectious virus $V(t)$. Here, we follow Wu et al., [2008] and consider a
model of HIV-1 dynamics where the target cells are CD$4^{+}$ T-cells. These
cells are produced at a constant rate $\lambda$ and cleared linearly at rate
$d$. Infection occurs at a rate $\beta$ following contact between a target
cell and infectious viral particle and infected cells are cleared at rate
$\delta$. Upon lysis, infected cells release $N$ viral particles into the
circulation and free virus is cleared at a constant rate $c$. The viral
dynamics model is given by

$\left.\begin{aligned} \frac{\textrm{d}}{\textrm{dt}}T(t)&=\lambda-\beta
T(t)V(t)-dT(t)\\\ \frac{\textrm{d}}{\textrm{dt}}I(t)&=\beta T(t)V(t)-\delta
I(t)\\\ \frac{\textrm{d}}{\textrm{dt}}V(t)&=\delta
NI(t)-cV(t).\end{aligned}\right\\}$ (13)

It is common to set $p=\delta N$ so the final equation for $V(t)$ becomes

$\displaystyle\frac{\textrm{d}}{\textrm{dt}}V(t)$ $\displaystyle=pI(t)-cV(t),$

and the system (13) is equipped with initial conditions
$T(0)=T_{0},I(0)=I_{0},$ and $V(0)=V_{0}$. In typical clinical studies,
temporal data is only collected for circulating free virus so the model output
corresponding to the calibration measurements is

$\displaystyle y_{i}(\theta)=\log_{10}(V(t_{i},\theta)),$

where using $log_{10}$ measurements of viral load is standard in HIV studies.

During antiretroviral therapy (ART), the viral load may fall below the limit
of detection of standard assays. While there are a number of techniques to
account for this censored data, we do not consider data collected during ART,
so the objective function is given by the sum of squares error

$\displaystyle
G_{HIV}(\theta,\phi)=\sqrt{\displaystyle\sum_{i=1}^{n}\left(\log_{10}(V(t_{i},\theta)-\log_{10}(\phi_{i})\right)^{2}}.$
(14)

Wu et al., [2008] characterized the identifiability of this model using a
higher order derivative method. They found that, if the initial conditions of
the model $T_{0},I_{0},$ and $V_{0}$ are known, then all six model parameters
$\theta=\\{\beta,d,\delta,c,N,\lambda,\\}$ are identifiable. To illustrate
their results, they fixed $\theta=\\{(2\times 10^{-5},0.15,0.55,5.5,900,80\\}$
and simulated the ODE model (13). They sampled the simulated viral load at
$37$ distinct time points and added noise $\varepsilon_{i}$ sampled from a
Gaussian distribution with $\mu=0$ and $\sigma^{2}=1$ [Wu et al.,, 2008].

In Section 3.2, we demonstrated the effectiveness of our continuation
technique by focusing on objective value function and computational efficiency
in calculating the MLE. Here, we illustrate how model dynamics evolve during
MLE continuation. We follow Wu et al., [2008] but consider a smaller subset of
calibration data collected at time $t_{i}=\\{0.4,1,8,14,20,36,46,58\\}$. We
add noise $\varepsilon_{i}^{0}$ sampled from a Gaussian distribution with
$\mu=0$ and $\sigma^{2}=0.15$ so the initial calibration data is

$\displaystyle\phi_{i}^{0}=\log_{10}(V(t_{i},\theta))+\varepsilon_{i}^{0}.$

We first fit the model to the simulated data $\phi_{i}^{0}$ to obtain an
initial MLE. We then generate 4 additional viral load time courses
$\\{\phi_{i}^{j}\\}_{j=1}^{10}$ by

$\displaystyle\phi_{i}^{j}=\phi_{i}^{0}+h_{step}|\varepsilon_{i}^{j}|$

for $\varepsilon_{i}^{j}$ sampled from a Gaussian distribution with $\mu=0$
and $\sigma^{2}=1$ and $h_{step}=\pm 0.1,\pm 0.2$. This collection of $4$ data
sets could feasibly represent experimental data measured from an increasingly
large sample drawn from a population of HIV-1 positive individuals with
population viral dynamic parameters given by $\theta=\\{(2\times
10^{-5},0.15,0.55,5.5,900,80\\}$. Here, we test the ability of our
continuation technique to predict reasonable viral dynamic curves without
refitting the data.

In Figure 3 A), we compute the predicted $\Psi(\phi^{j})$ and plot the
predicted model dynamics obtained from $\Psi(\phi^{j})$ against the perturbed
data $\phi^{j}$. In Figure 3 B), we show the fit model predictions to the
perturbed data. In each case, the viral dynamics show comparable model
predictions for the fit and predicted model parameters demonstrating that our
continuation method can successfully predict reasonable model simulations. In
fact, the Bayesian Information Criteria [Kass and Raftery,, 1995] indicates no
significant differences between the predicted and true MLE for all 4 data
sets. However, Figure 3 C) shows the significant computational improvement
obtained by only calculating the continuation step rather than fitting all
model parameters at each step. The predicted model dynamics track the true
viral load trajectory.

Figure 3:  Comparison of predicted model fits to randomly perturbed data.
Panels A and B show model trajectories obtained using predicted model
parameters to the simulated experimental data perturbed by
$\phi_{i}^{j}=\phi_{i}^{0}+h_{step}|\varepsilon_{i}^{j}|$. Panel A shows the
predicted model fits to the experimental data while B shows the model fits to
data resulting from the true MLE. Panel C shows the number of objective value
evaluations required to predict the MLE using this continuation technique or
fit the model parameters to the perturbed data using the known parameters as a
starting guess.

It is common to find numerous local minima of (14) when fitting (13) to
simulated data. As measured by the value of the log-likelihood function or
information criteria, these local minima can produce comparable fits to a
given data set despite different dynamics. We perturbed the initial data set
$\phi_{0}$ by

$\displaystyle\log(\phi_{i}^{1})=\log(\phi_{i}^{0})+0.8\varepsilon_{i}$

for $\varepsilon_{i}$ sampled from a Gaussian distribution with $\mu=0$ and
$\sigma^{2}=1$. We fit this perturbed data from 10 distinct initial guesses
using fmincon [MATLAB,, 2017]. These 10 starting initial guesses converged to
two local minima. We denote the corresponding parameter estimates by
$\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ and plot the resulting model
trajectories in Fig 4. These fits are indistinguishable by BIC and both appear
to accurately describe the viral load data. Consequently, it is not obvious
which of $\hat{\theta}_{1}$ and $\hat{\theta}_{2}$ best describe the data.

However, it is reasonable to expect that the MLE should be robust to small
perturbations of the calibration data. We measure the robustness of each of
these minima by calculating $\|\mathrm{D}\Psi(\phi^{1})\|$ at
$\hat{\theta}_{1}$ and $\hat{\theta}_{2}$. A smaller norm
$\|\mathrm{D}\Psi(\phi^{1})\|$ implies less sensitivity of the MLE to
perturbations of the calibration data. For the example shown in Fig 4, there
is a 16 fold difference in sensitivity to calibration data. In this way,
$\mathrm{D}\Psi$ can be used to distinguish between otherwise similar fits. We
suggest that, when choosing between multiple fits with similar BIC values, the
parameter estimate with the smaller sensitivity to the data is a more robust,
and thus preferential, fit.

Figure 4:  Comparison of two potential fits to randomly perturbed viral
dynamics models. Model trajectories obtained from two local minima from
fitting 10 initial guesses to viral load data shown in black. Both
trajectories accurately describe the viral load dynamics as evidenced by a
small difference in BIC. However, the parameter estimate corresponding to the
oscillatory trajectory is much more robust as measured by
$\|\mathrm{D}\Psi(\phi^{1})\|$.

##  4 Discussion

Parameter fitting is crucial step when using mathematical models to predict
novel treatment strategies, extrapolate from clinical trials, identify new
drug targets or schedules, or propose non-pharmaceutical interventions [Brady
and Enderling,, 2019; Cassidy et al.,, 2020; Cassidy and Craig,, 2019].
However, parameter fitting can be difficult and computationally expensive. A
large variety of fitting techniques have therefore been developed to calibrate
model predictions against data [Toni et al.,, 2009; Horbelt et al.,, 2002;
Kreutz et al.,, 2013; Lauß et al.,, 2018]. Moreover, mathematical modeling is
increasingly applied to understand emerging data and make real-time
predictions. In this case, as new data emerges, the model parameters must be
refit with potential computational cost. Here, we developed a continuation
type technique to quantify how updates to experimental data will impact the
MLE and predict the evolution of the MLE as a function of the experimental
data used to calibrate the model.

We used the implicit function theorem to calculate the trajectory of the MLE
through parameter space. As the implicit function theorem only guarantees the
existence of a differentiable trajectory $\Psi$ through calibration
data–parameter space, we utilized the first order Taylor expansion $\Psi$ to
extrapolate the evolution of the MLE due to changes in experimental data. We
showed how this calculation is intrinsically linked to local sensitivity
analysis and the curvature of the objective function. In two examples drawn
from mathematical biology, we showed how this continuation technique can
predict acceptable model fits to experimental data while significantly
reducing computational overhead. In fact, in most applications, our
continuation technique requires no dedicated computational overhead as the
Hessian of the objective function is calculated at each step when using common
optimization algorithms, such as fmincon [MATLAB,, 2017], and local
sensitivity analysis is a standard step in model fitting.

Perhaps more importantly that gains in computational efficiency, our approach
explicitly identifies relationships between individual experimental
measurements and parameter estimates. Our approach addresses similar questions
to local sensitivity analysis from a distinct perspective. Rather than using
simulations to understand how small perturbations in model parameters from the
best-fit parameters change model outputs as in standard sensitivity analysis,
we quantify how changes in the training data impact the best-fit parameters
and measure the sensitivity of the best-fit parameters to variations in this
calibration data. As we showed in Section 3.2, this perspective can be used to
suggest additional experimental measurements to increase confidence in model
parameterization. Further, we showed how to use $D\Psi$ to understand which
experimental measurements are most informative for model parameterizations and
identify redundant measurements that do not provide additional information for
parameter estimation.

Our technique is a type of local analysis that explores the functional
dependence of the MLE on experimental data starting from a pre-identified MLE.
Specifically, we assume that the Hessian of the objective function is
invertible at the MLE and our results are necessarily local in parameter space
as we are extrapolating from a pre-identified MLE. Nevertheless, our examples
show the utility of our continuation approach for even large perturbations of
the experimental data.

Despite these limitations, we developed a continuation-type technique to
predict the functional dependence of a MLE on the experimental data used to
train a mathematical model. While we have focused on applications in
mathematical biology, our approach is immediately portable to other domains.
As our method is independent of the number of data points, our approach could
be particularly useful in big-data applications. Ultimately, our results offer
a unified approach to quantify the relationship between training data and
best-fit model parameters and to leverage this understanding to suggest
additional experiments to increase confidence in model parameterization.

## Data access statement

The code and data underlying the results in this manuscript is available at
https://github.com/ttcassid/MLE$\\_$Continuation.

## References

  * Allen et al., [2016] Allen, R. J., Rieger, T. R., and Musante, C. J. (2016).  Efficient Generation and Selection of Virtual Populations in Quantitative Systems Pharmacology Models.  CPT: Pharmacometrics & Systems Pharmacology, 5(3):140–146.
  * Altrock et al., [2015] Altrock, P. M., Liu, L. L., and Michor, F. (2015).  The mathematics of cancer: integrating quantitative models.  Nature Reviews Cancer, 15(12):730–745.
  * Brady and Enderling, [2019] Brady, R. and Enderling, H. (2019).  Mathematical Models of Cancer: When to Predict Novel Therapies, and When Not to.  Bulletin of Mathematical Biology, 81(10):3722–3731.
  * [4] Braniff, N., Richards, A., and Ingalls, B. (2019a).  Optimal Experimental Design for a Bistable Gene Regulatory Network.  IFAC-PapersOnLine, 52(26):255–261.
  * [5] Braniff, N., Scott, M., and Ingalls, B. (2019b).  Component Characterization in a Growth-Dependent Physiological Context: Optimal Experimental Design.  Processes, 7(1):52.
  * Cárdenas et al., [2022] Cárdenas, S. D., Reznik, C. J., Ranaweera, R., Song, F., Chung, C. H., Fertig, E. J., and Gevertz, J. L. (2022).  Model-informed experimental design recommendations for distinguishing intrinsic and acquired targeted therapeutic resistance in head and neck cancer.  npj Systems Biology and Applications, 8(1):32.
  * Cassidy and Craig, [2019] Cassidy, T. and Craig, M. (2019).  Determinants of combination GM-CSF immunotherapy and oncolytic virotherapy success identified through in silico treatment personalization.  PLOS Computational Biology, 15(11):e1007495.
  * Cassidy et al., [2020] Cassidy, T., Humphries, A. R., Craig, M., and Mackey, M. C. (2020).  Characterizing Chemotherapy-Induced Neutropenia and Monocytopenia Through Mathematical Modelling.  Bulletin of Mathematical Biology, 82(8):104.
  * Cassidy et al., [2021] Cassidy, T., Nichol, D., Robertson-Tessi, M., Craig, M., and Anderson, A. R. A. (2021).  The role of memory in non-genetic inheritance and its impact on cancer treatment resistance.  PLOS Computational Biology, 17(8):e1009348.
  * Castro and de Boer, [2020] Castro, M. and de Boer, R. J. (2020).  Testing structural identifiability by a simple scaling method.  PLOS Computational Biology, 16(11):e1008248.
  * Craig et al., [2019] Craig, M., Kaveh, K., Woosley, A., Brown, A. S., Goldman, D., Eton, E., Mehta, R. M., Dhawan, A., Arai, K., Rahman, M. M., Chen, S., Nowak, M. A., and Goldman, A. (2019).  Cooperative adaptation to therapy (CAT) confers resistance in heterogeneous non-small cell lung cancer.  PLOS Computational Biology, 15(8):e1007278.
  * Davies et al., [2020] Davies, N. G., Kucharski, A. J., Eggo, R. M., Gimma, A., Edmunds, W. J., Jombart, T., O’Reilly, K., Endo, A., Hellewell, J., Nightingale, E. S., Quilty, B. J., Jarvis, C. I., Russell, T. W., Klepac, P., Bosse, N. I., Funk, S., Abbott, S., Medley, G. F., Gibbs, H., Pearson, C. A. B., Flasche, S., Jit, M., Clifford, S., Prem, K., Diamond, C., Emery, J., Deol, A. K., Procter, S. R., van Zandvoort, K., Sun, Y. F., Munday, J. D., Rosello, A., Auzenbergs, M., Knight, G., Houben, R. M. G. J., and Liu, Y. (2020).  Effects of non-pharmaceutical interventions on COVID-19 cases, deaths, and demand for hospital services in the UK: a modelling study.  The Lancet Public Health, 5(7):e375–e385.
  * De Souza and Humphries, [2019] De Souza, D. C. and Humphries, A. R. (2019).  Dynamics of a Mathematical Hematopoietic Stem-Cell Population Model.  SIAM Journal on Applied Dynamical Systems, 18(2):808–852.
  * Dhooge et al., [2008] Dhooge, A., Govaerts, W., Kuznetsov, Y. A., Meijer, H. G., and Sautois, B. (2008).  New features of the software MatCont for bifurcation analysis of dynamical systems.  Mathematical and Computer Modelling of Dynamical Systems, 14(2):147–175.
  * Goldman et al., [2015] Goldman, A., Majumder, B., Dhawan, A., Ravi, S., Goldman, D., Kohandel, M., Majumder, P. K., and Sengupta, S. (2015).  Temporally sequenced anticancer drugs overcome adaptive resistance by targeting a vulnerable chemotherapy-induced phenotypic transition.  Nature Communications, 6(1):6139.
  * Gunnarsson et al., [2020] Gunnarsson, E. B., De, S., Leder, K., and Foo, J. (2020).  Understanding the role of phenotypic switching in cancer drug resistance.  Journal of Theoretical Biology, 490:110162.
  * Holmdahl and Buckee, [2020] Holmdahl, I. and Buckee, C. (2020).  Wrong but Useful – What Covid-19 Epidemiologic Models Can and Cannot Tell Us.  New England Journal of Medicine, 383(4):303–305.
  * Horbelt et al., [2002] Horbelt, W., Timmer, J., and Voss, H. (2002).  Parameter estimation in nonlinear delayed feedback systems from noisy data.  Physics Letters A, 299(5-6):513–521.
  * Jenner et al., [2021] Jenner, A. L., Cassidy, T., Belaid, K., Bourgeois-Daigneault, M.-C., and Craig, M. (2021).  In silico trials predict that combination strategies for enhancing vesicular stomatitis oncolytic virus are determined by tumor aggressivity.  Journal for ImmunoTherapy of Cancer, 9(2):e001387.
  * Jolly et al., [2018] Jolly, M. K., Kulkarni, P., Weninger, K., Orban, J., and Levine, H. (2018).  Phenotypic Plasticity, Bet-Hedging, and Androgen Independence in Prostate Cancer: Role of Non-Genetic Heterogeneity.  Frontiers in Oncology, 8(MAR):1–12.
  * Kass and Raftery, [1995] Kass, R. E. and Raftery, A. E. (1995).  Bayes Factors.  Journal of the American Statistical Association, 90(430):773–795.
  * Kreutz et al., [2013] Kreutz, C., Raue, A., Kaschek, D., and Timmer, J. (2013).  Profile likelihood in systems biology.  FEBS Journal, 280(11):2564–2571.
  * Kreutz et al., [2012] Kreutz, C., Raue, A., and Timmer, J. (2012).  Likelihood based observability analysis and confidence intervals for predictions of dynamic models.  BMC Systems Biology, 6.
  * Kreutz and Timmer, [2009] Kreutz, C. and Timmer, J. (2009).  Systems biology: Experimental design.  FEBS Journal, 276(4):923–942.
  * Lauß et al., [2018] Lauß, T., Oberpeilsteiner, S., Steiner, W., and Nachbagauer, K. (2018).  The discrete adjoint method for parameter identification in multibody system dynamics.  Multibody System Dynamics, 42(4):397–410.
  * Li and Vu, [2013] Li, P. and Vu, Q. D. (2013).  Identification of parameter correlations for parameter estimation in dynamic biological models.  BMC Systems Biology, 7(1):91.
  * Li and Vu, [2015] Li, P. and Vu, Q. D. (2015).  A simple method for identifying parameter correlations in partially observed linear dynamic models.  BMC Systems Biology, 9(1):92.
  * Li et al., [2011] Li, Y., O’Dell, S., Walker, L. M., Wu, X., Guenaga, J., Feng, Y., Schmidt, S. D., McKee, K., Louder, M. K., Ledgerwood, J. E., Graham, B. S., Haynes, B. F., Burton, D. R., Wyatt, R. T., and Mascola, J. R. (2011).  Mechanism of Neutralization by the Broadly Neutralizing HIV-1 Monoclonal Antibody VRC01.  Journal of Virology, 85(17):8954–8967.
  * Luo et al., [2022] Luo, M. C., Nikolopoulou, E., and Gevertz, J. L. (2022).  From Fitting the Average to Fitting the Individual: A Cautionary Tale for Mathematical Modelers.  Frontiers in Oncology, 12(April):8–11.
  * Maiwald et al., [2016] Maiwald, T., Hass, H., Steiert, B., Vanlier, J., Engesser, R., Raue, A., Kipkeew, F., Bock, H. H., Kaschek, D., Kreutz, C., and Timmer, J. (2016).  Driving the model to its limit: Profile likelihood based model reduction.  PLoS ONE, 11(9):1–18.
  * Marino et al., [2008] Marino, S., Hogue, I., Ray, C., and Kirschner, D. (2008).  A methodology for performing global uncertainty and sensitivity analysis in systems biology.  Journal of Theoretical Biology, 254:178–196.
  * MATLAB, [2017] MATLAB (2017).  R2017a.  The MathWorks Inc., Natick, Massachusetts.
  * Perelson, [2002] Perelson, A. S. (2002).  Modelling viral and immune system dynamics.  Nature Reviews Immunology, 2(1):28–36.
  * Raue et al., [2014] Raue, A., Karlsson, J., Saccomani, M. P., Jirstrand, M., and Timmer, J. (2014).  Comparison of approaches for parameter identifiability analysis of biological systems.  Bioinformatics, 30(10):1440–1448.
  * Raue et al., [2009] Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U., and Timmer, J. (2009).  Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood.  Bioinformatics, 25(15):1923–1929.
  * Sahoo et al., [2021] Sahoo, S., Mishra, A., Kaur, H., Hari, K., Muralidharan, S., Mandal, S., and Kumar Jolly, M. (2021).  A mechanistic model captures the emergence and implications of non-genetic heterogeneity and reversible drug resistance in ER+ breast cancer cells.  bioRxiv, page 2021.03.14.435359.
  * Sanche et al., [2022] Sanche, S., Cassidy, T., Chu, P., Perelson, A. S., Ribeiro, R. M., and Ke, R. (2022).  A simple model of COVID-19 explains disease severity and the effect of treatments.  Scientific Reports, 12(1):14210.
  * Sanche et al., [2020] Sanche, S., Lin, Y. T., Xu, C., Romero-Severson, E., Hengartner, N., and Ke, R. (2020).  High Contagiousness and Rapid Spread of Severe Acute Respiratory Syndrome Coronavirus 2.  Emerging Infectious Diseases, 26(7):1470–1477.
  * Sharp et al., [2022] Sharp, J. A., Browning, A. P., Burrage, K., and Simpson, M. J. (2022).  Parameter estimation and uncertainty quantification using information geometry.  Journal of The Royal Society Interface, 19(189).
  * Silk et al., [2014] Silk, D., Kirk, P. D. W., Barnes, C. P., Toni, T., and Stumpf, M. P. H. (2014).  Model Selection in Systems Biology Depends on Experimental Design.  PLoS Computational Biology, 10(6):e1003650.
  * Stapor et al., [2018] Stapor, P., Fröhlich, F., and Hasenauer, J. (2018).  Optimization and profile calculation of ODE models using second order adjoint sensitivity analysis.  Bioinformatics, 34(13):i151–i159.
  * Thompson, [2020] Thompson, R. N. (2020).  Epidemiological models are important tools for guiding COVID-19 interventions.  BMC Medicine, 18(1):152.
  * Toni et al., [2009] Toni, T., Welch, D., Strelkowa, N., Ipsen, A., and Stumpf, M. P. H. (2009).  Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems.  Journal of the Royal Society, Interface, 6(31):187–202.
  * Wu et al., [2008] Wu, H., Zhu, H., Miao, H., and Perelson, A. S. (2008).  Parameter Identifiability and Estimation of HIV/AIDS Dynamic Models.  Bulletin of Mathematical Biology, 70(3):785–799.
  * Zhang et al., [2022] Zhang, J., Cunningham, J., Brown, J., and Gatenby, R. (2022).  Evolution-based mathematical models significantly prolong response to abiraterone in metastatic castrate-resistant prostate cancer and identify strategies to further improve outcomes.  eLife, 11:1–105.
  * Zi, [2011] Zi, Z. (2011).  Sensitivity analysis approaches applied to systems biology models.  IET Systems Biology, 5(6):336–346.