# Understanding the impact of numerical solvers on inference for differential
equation models

Richard Creswell  Department of Computer Science, University of Oxford
Katherine M. Shepherd  Department of Computer Science, University of Oxford
Ben Lambert  Department of Mathematics, University of Exeter  Gary R. Mirams
School of Mathematical Sciences, University of Nottingham  Chon Lok Lei
Institute of Translational Medicine, Faculty of Health Sciences, University of
Macau; Department of Biomedical Sciences, Faculty of Health Sciences,
University of Macau  Simon Tavener  Department of Mathematics, Colorado State
University  Martin Robinson  Department of Computer Science, University of
Oxford  David J. Gavaghan  Department of Computer Science, University of
Oxford

###### Abstract

Most ordinary differential equation (ODE) models used to describe biological
or physical systems must be solved approximately using numerical methods.
Perniciously, even those solvers which seem sufficiently accurate for the
_forward problem_ , i.e., for obtaining an accurate simulation, may not be
sufficiently accurate for the _inverse problem_ , i.e., for inferring the
model parameters from data. We show that for both fixed step and adaptive step
ODE solvers, solving the forward problem with insufficient accuracy can
distort likelihood surfaces, which may become jagged, causing inference
algorithms to get stuck in local “phantom” optima. We demonstrate that biases
in inference arising from numerical approximation of ODEs are potentially most
severe in systems involving low noise and rapid nonlinear dynamics. We
reanalyze an ODE changepoint model previously fit to the COVID-19 outbreak in
Germany and show the effect of the step size on simulation and inference
results. We then fit a more complicated rainfall-runoff model to hydrological
data and illustrate the importance of tuning solver tolerances to avoid
distorted likelihood surfaces. Our results indicate that when performing
inference for ODE model parameters, adaptive step size solver tolerances must
be set cautiously and likelihood surfaces should be inspected for
characteristic signs of numerical issues.

Keywords: ordinary differential equations, inference, Bayesian statistics,
truncation error, compartmental models, hydrological modelling

##  1 Introduction

Many scientific phenomena involve time-varying signals or outputs. These
phenomena are often believed to obey some parametric model, whose parameter
values are _a priori_ unknown but can be inferred from observed data. In this
paper, we present and analyze the key challenges that arise in parameter
inference when models involve ordinary differential equations (ODEs). ODEs are
used throughout the biological and physical sciences to express dynamic
processes; a few examples amongst the myriad of their application areas
include epidemiology [van der Vegt et al., 2022], hydrology [Kavetski et al.,
2003], cardiac electrophysiology [Whittaker et al., 2020], and population
dynamics [Shertzer et al., 2002].

In general, the differential equations used in scientific applications cannot
be solved analytically. However, a wide range of computational methods have
been developed to obtain a numerical approximation of their solutions.
(Solving the differential equation for a particular value of the parameters is
known as the _forward problem_.) While numerical algorithms to solve the
forward problem introduce error, the properties of this error are generally
well understood and can be controlled. In solvers using a fixed time step
(discussed in §3.2), the error can be reduced by decreasing the size of the
time step [Gautschi, 1997]. In solvers in which the time step is set
adaptively (discussed in §3.3), the error is typically controlled through
user-specified relative and absolute tolerances on the local truncation error
(the error in the solution introduced by a single time step of the solver)
[Dormand and Prince, 1980].

Our focus in this paper is on the interplay between the numerical
approximations inherent in the forward problem and the _inverse problem_ ,
which consists of learning values of the parameters that are compatible with
an observed time series. Widely used approaches to the inverse problem include
optimization of an objective function which measures the quality of fit
between the model and the data (e.g., maximum likelihood), or Bayesian
approaches which generate samples from the posterior distribution of the
parameters (e.g., Markov chain Monte Carlo (MCMC)). These approaches to the
inverse problem require the forward problem to be solved at multiple different
parameter values. The errors in each numerical solution of the forward
problem, even when individually small, are liable to introduce significant
bias to inference results.

The rest of this paper is organised as follows. In §2, we present the widely
used independent and identically distributed Gaussian noise log-likelihood
function for fitting ODE models and derive a bound on the error in this log-
likelihood arising from the use of an approximate solution to the ODE. On the
basis of this bound, and results presented subsequently, we argue that the
biases in inference results arising from numerical solvers are likely to be
most severe in systems which have low noise and rapid nonlinear dynamics. In
§3, we study two broad classes of ODE solvers: those involving a fixed time
step, and those involving a time step set adaptively to control the error on
the solution. For both classes of solver, we illustrate the effects that
solver inaccuracy may have on inference, and illustrate this using synthetic
data. Additionally, in §S3, we study how smoothing approximations can reduce
the influence of numerical error on computation of the likelihood. Finally, in
§4 and §5 we consider inference of ODE models for real data series. in §4 we
reanalyze an ODE model of disease transmission fit to the COVID-19 outbreak in
Germany and show that, when using a solver with a fixed time step, the choice
of time step can alter inference and simulation results, and in §5 we fit a
rainfall-runoff model to hydrological data to illustrate the pitfalls of
performing parameter inference using an adaptive step size solver with
insufficient local tolerances.

##  2 Effects of numerical error on computation of the log-likelihood

###  2.1 Log-likelihood function for an ODE model

We assume that time series data $\\{y_{i}\\}_{i=1}^{N}$;
$y_{i}\in\mathbb{R}^{n}$ are measured at time points $\\{t_{i}\\}_{i=1}^{N}$.
These data are believed to obey some function
$g:\mathbb{R}^{l}\rightarrow\mathbb{R}^{n}$ of $x(t;\theta)\in\mathbb{R}^{l}$,
where $x$ is the solution to an ordinary differential equation:

$\displaystyle\frac{dx}{dt}$ $\displaystyle=h(t,x,\theta);$ (1) $\displaystyle
x(0)$ $\displaystyle=x_{0},$ (2)

for some function $h$ which is informed by the relevant scientific theory and
parameterized by the (potentially unknown) values $\theta\in\mathbb{R}^{m}$.
(In some inference problems, the initial value $x_{0}$ is also unknown and
inferred from the data.) Eq. (1) has been written as a _first order_ equation
(i.e., only involving the first order derivative of $x$); higher order
differential equations may be rewritten as systems of first order equations.

Deterministic forward models never fully explain the variation in real
observations. To make the forward model a feasible explanation of the data, an
additional stochastic component representing unmodelled elements (for example,
processes involved in the measurement of the signal) is included. This
stochastic component is often included in additive form, yielding the proposed
data generating process:

$y_{i}=g(x(t_{i};\theta))+\epsilon_{i},$ (3)

where $\epsilon_{i}$ is a mean-zero random variable specifying the noise
process. Many choices for $\epsilon_{i}$ are possible, and the assumed
distribution of this variable should be chosen carefully, as misspecification
of $\epsilon_{i}$ may cause inaccurate inference results
[Lambert et al., 2023]. A standard choice, which we use for synthetic data
generation in this paper, is the independent and identically distributed (IID)
Gaussian:

$\epsilon_{i}\overset{\text{IID}}{\sim}N(0,\sigma),$ (4)

with the parameter $\sigma$ representing the standard deviation of the noise
process, which is also inferred from the data together with the model
parameters $\theta$. Choosing a particular noise process enables the joint
probability of the observations to be expressed as a function of the parameter
values $\theta$—this is known as a likelihood. For an IID Gaussian noise
process (eq. (4)), the log-likelihood for time series data
$\\{y_{i}\\}_{i=1}^{N}$ takes the form:

$\log
p(y_{1},\dots,y_{N}|\theta,\sigma)=-\frac{N}{2}\log(2\pi)-\frac{N}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-g(x(t_{i};\theta)))^{2}.$
(5)

The likelihood expresses the quality of the fit between the model output and
the data, with higher values of the likelihood indicating a superior fit.
Thus, the values of $\theta$ most compatible with data can be found by
maximising the likelihood with respect to $\theta$ (i.e., the method of
maximum likelihood). Alternatively, the likelihood can be used together with a
prior distribution on the parameters ($p(\theta)$) to infer the posterior
distribution according to Bayes’ theorem:

$p(\theta|y_{1},\dots,y_{N})=\frac{p(y_{1},\dots,y_{N}|\theta)p(\theta)}{p(y)}.$

In this paper, we consider the typical case where the posterior cannot easily
be expressed in closed form but can be approximated using Markov chain Monte
Carlo (MCMC) sampling methods [Gelman et al., 2013].

###  2.2 Error in the log-likelihood arising from approximation of the forward
solution

The data are assumed to obey the IID Gaussian log-likelihood, eq. (5). We
assume that $x(t_{i};\theta)$ is the _true_ solution to the ODE at time point
$t_{i}$, which is unavailable and approximated by $\hat{x}_{i}$. The deviation
between $x(t_{i};\theta)$ and $\hat{x}_{i}$ at any time point is given by the
global truncation error, $e(t_{i})$:

$e(t_{i})=x(t_{i};\theta)-\hat{x}_{i}.$

In general, $e(t_{i})$ is unknown, although, for particular numerical solvers,
its magnitude may be bounded by some function of the step size or some other
quantity which can be used to tune the accuracy of the solver.

The log-likelihood available to the inference algorithm takes the same form as
eq. (5), but computed using the numerical approximation $\hat{x}_{i}$ instead
of $x(t_{i};\theta)$. For brevity, we denote the accurate log-likelihood by
$\mathcal{L}$, and we denote the log-likelihood computed using a numerical
approximation by $\mathcal{L}^{\prime}$, which is given by

$\mathcal{L}^{\prime}=-\frac{N}{2}\log(2\pi)-\frac{N}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-g(\hat{x}_{i}))^{2}.$
(6)

Assuming Lipschitz continuity of the observation function $g$ with Lipschitz
constant $K$, we prove in §S1 that the difference between
$\mathcal{L}^{\prime}$ and $\mathcal{L}$ is bounded according to:

$|\mathcal{L}-\mathcal{L}^{\prime}|\leq\sum_{i=1}^{N}\left(\frac{K^{2}}{2\sigma^{2}}|e(t_{i})|^{2}+\frac{K}{\sigma^{2}}|e(t_{i})||y_{i}-g(x(t_{i};\theta))|\right).$
(7)

We observe an inverse relationship between $\sigma$ and the bound of
$|\mathcal{L}-\mathcal{L}^{\prime}|$ when $e(t_{i})$ is held constant. Thus,
when a solver is tuned to yield a particular global truncation error
$e(t_{i})$, we expect the absolute bias in the log-likelihood to be more
severe at smaller values of $\sigma$. Eq. (7) also indicates that more severe
biases may be expected when $N$ is larger—i.e., there are more observations.
Additionally, at a fixed level of $\sigma$, we expect the bias in the log-
likelihood to decrease as the global truncation error is decreased. In §S2, we
additionally derive the distribution of the error in the likelihood, and show
that $\mathbb{E}[\mathcal{L}-\mathcal{L}^{\prime}]>0$: the numerical
approximation of the likelihood will, on average, underestimate the true
likelihood.

##  3 Effects of ODE solvers on inference

To study the interplay between ODE solvers and inference, we introduce the
following differential equation problem which describes an oscillatory system
with damping and forcing:

$m\frac{d^{2}x}{dt^{2}}+c\frac{dx}{dt}+kx=F(t).$

The model has three parameters which will be treated as unknown: $(m,c,k)$. In
classical mechanics, these represent the mass, damping coefficient, and spring
constant respectively. $F(t)$ represents the forcing function or stimulus, and
in this paper takes a variety of forms throughout our results. This damped and
forced oscillator is described by a second order differential equation; to
apply ODE solvers straightforwardly, we rewrite it as a first order
differential equation of two state variables:

$\frac{d}{dt}\begin{pmatrix}x\\\
\dot{x}\end{pmatrix}=\begin{pmatrix}\dot{x}\\\
\frac{F(t)}{m}-\frac{c}{m}\dot{x}-\frac{k}{m}x\end{pmatrix},$ (8)

where $\dot{x}=dx/dt$.

###  3.1 Fixed step and adaptive step ODE solvers

A wide range of numerical algorithms have been developed to obtain approximate
solutions to initial value problems (IVPs) of the form given in eq. (1). These
algorithms typically work by computing an approximate solution on a grid of
time points (in general, distinct from the time points where the data are
located) and then using an interpolation algorithm to obtain the solution at
intermediate time points.

Most simply, the grid of solver time points can be prespecified in advance (we
refer to such methods as fixed time step solvers). However, in general, it is
inefficient to use the same time step throughout the entire time range on
which the ODE is being solved, particularly when solved repeatedly over a
range of parameters. Solvers can employ large time steps in regions where the
solution and its gradients change gradually without causing much error in the
solution; however, in regions where the derivative changes rapidly, small time
steps are required to maintain a low error. This motivated the development of
ODE solvers which adjust the step size throughout the time domain over which
the ODE is solved. While fixed step solvers are still commonly used, adaptive
step solvers are standard in high performance computing and are widely
implemented in software libraries for ODE solving.

When using an adaptive step size solver, the user does not specify a step
size, but rather a local error tolerance. The algorithm then selects a time-
varying sequence of step sizes such that the local error in the solution falls
below the specified tolerance. The total number of time steps used by the
solver thus depends on the selected tolerance and the properties of the
solution. Typically, an interpolation scheme is then used to obtain the
solution at intermediate time values. Tolerances can be expressed either as an
absolute value or relative to the magnitude of the solution. In many
implementations, both are available to the user: for example, the SciPy
library allows the user to specify both an absolute tolerance atol and a
relative tolerance rtol, and chooses step sizes such that the magnitude of the
local truncation error on the solution $x$ does not exceed
$\texttt{atol}+\texttt{rtol}|x|$ [Virtanen et al., 2020]. For the results
presented in this paper, we fix atol to a value of $10^{-9}$ and tune rtol to
control the accuracy of the solver. Adaptive step sizes have been implemented
for a wide variety of ODE solver algorithms. Here, we focus on Runge-Kutta
methods of the form RK$p(q)$, which use the $q$th order method to estimate the
error (and thus control the time step), while making the actual steps using
the $p$th order method [Dormand and Prince, 1980]. Runge-Kutta methods are not
described in detail here for brevity—they are widely used and details can be
found in many standard texts (for example, [Gautschi, 1997]). We rely on the
SciPy adaptive time step Runge-Kutta implementation, which employs a quartic
interpolation polynomial for RK5(4) and a cubic Hermite interpolation
polynomial for RK3(2) [Virtanen et al., 2020].

####  3.1.1 Typical log-likelihood surface shapes

We now consider the influence of the two numerical solution methods for
parameter inference. Because fixed step solvers use the same grid throughout
parameter space, while adaptive step solvers may employ different grids at
different parameter values, these two classes of solvers differ in the
characteristics of the error that they may introduce into the likelihood
function.

We illustrate this by computing the likelihood surface for the $k$ parameter
in the oscillator problem, eq. (8). $75$ evenly spaced data points were
generated from and including $t=0$ to $t=50$ from the model with an accurate
solver (the RK5(4) solver with relative tolerance set to $10^{-8}$), using
true parameter values $k=1$, $c=0.2$, $m=1$, initial conditions of $x(t=0)=0$,
$\dot{x}(t=0)=0$, and

$F(t)=\begin{cases}1,&t<25,\\\ 0.9,&t\geq 25.\end{cases}$

Then IID Gaussian noise was added to the solution at each of the sampled
locations with $\sigma=0.01$. Holding all other parameters fixed at their true
values, the log-likelihood was calculated for a range of values of $k$, using
three different ODE solvers. First, the RK5(4) solver with relative tolerance
set to $10^{-8}$ was used to compute the accurate (‘True’) likelihood. Next,
the Forward Euler solver with a fixed time step of $\Delta t=0.01$ was used.
Finally, we used the RK5(4) solver, but with its relative tolerance tuned so
the observed magnitude in the error in the log-likelihood at the true
parameter values was equal to that produced by the Forward Euler solver (for
this problem, this resulted in relative tolerance tuned to $0.00944$). These
results are shown in Figure 1.

At the true parameter value, both solvers result in a slight underestimation
of the log-likelihood. Across the parameter range considered, the fixed time
step solver results in a log-likelihood which is shifted relative to the true
one, but retains the smooth, unimodal shape. However, the adaptive step solver
results in a log-likelihood surface which in addition to being shifted
exhibits jagged, discontinuous fluctuations. In the remainder of §3, we
examine these two phenomena in more detail.

Figure 1: Comparison of log-likelihood surfaces calculated using fixed step
and adaptive step solvers. Log-likelihood for the parameter $k$ calculated
from data generated from the oscillator model eq. (8), with all other
parameters held at their true values. The log-likelihood was calculated from
eq. (5) using an adaptive step RK5(4) solver with relative tolerance set to
$10^{-8}$ (True), a Forward Euler solver with a fixed time step $\Delta
t=0.01$, and an adaptive step RK5(4) solver with tolerance tuned such that at
the true parameter values (vertical line) it introduces the same magnitude of
error in the log-likelihood as the fixed step Forward Euler solver
(corresponding to a relative tolerance of $0.00944$).

###  3.2 Fixed time step solvers

####  3.2.1 Forward Euler solver

One of the simplest numerical solvers for ordinary differential equations is
the Forward Euler method with a uniform step size $\Delta t$. This solver is
easily implemented and thus has achieved wide usage despite its simplicity and
typically mediocre performance.

Forward Euler has been used for inference in some recent high-profile
epidemiological research where $\Delta t$ was set to a value comparable to the
time scale of the behavior of the system (e.g., [Birrell et al., 2021, Dehning
et al., 2020]). Whether these applications are representative of the use of
Forward Euler more generally is unclear, but our results in §4 indicate that
such choices of $\Delta t$ may alter both forward model solutions and
parameter inference results.

####  3.2.2 Inference for the damped, driven oscillator using Forward Euler

We now exemplify the impact of using Forward Euler with insufficiently small
time steps on inference by using synthetic noisy data generated from the
(accurate) solution of eq. (8). $25$ evenly spaced data points were generated
from and including $t=0$ to $t=5$ from the model with an accurate solver (the
RK5(4) solver with relative tolerance set to $10^{-8}$), using true parameter
values $k=1$, $c=0.2$, $m=1$, an initial condition of $x(t=0)=0$,
$\dot{x}(t=0)=0$, and $F(t)=1$. Then, IID Gaussian noise was added to the
solution at each of the sampled locations with $\sigma=0.1$. Holding all other
parameters fixed at their true values, the log-likelihood was calculated for a
range of values of $k$, using the Forward Euler solver with various time
steps.

Figure 2 shows the impact of using Forward Euler on the likelihood surface.
The results show the typical effect of a fixed step solver with insufficiently
small time steps: the likelihood surface maintains a smooth shape, but it is
shifted relative to its true location. The longest time step considered in
this study, $\Delta t=0.1$, causes substantial inaccuracy in the likelihood
even though $\Delta t=0.1$ is small compared to the time scale of the dynamics
of the system and the system with $F(t)=1$ contains no discontinuities or
other challenging features.

As the step size is refined, the log-likelihood curves converge. This suggests
a diagnostic technique which could be incorporated into inference algorithms:
once the optimal parameter values have been determined, the log-likelihood
should be evaluated at those parameter values with the step size on the solver
slightly adjusted; if the solver is sufficiently accurate, the value of the
log-likelihood should not be a strong function of the step size.

Figure 2: Damped oscillator inference using Forward Euler. (Left) Synthetic
data for the damped driven oscillator. The curved line indicates the accurate
solution to the ODE with these parameters, while the points indicate the noisy
data. (Center) Solution for oscillator computed using a Forward Euler solver
with four different choices for the time step $\Delta t$. (Right) Log-
likelihood for the parameter $k$ calculated from the noisy data, with all
other parameters held at their true values. The log-likelihood was calculated
from eq. (5) using a Forward Euler solver with four different choices for the
time step $\Delta t$.

###  3.3 Adaptive step size solvers

Adaptive step size solvers enable increased efficiency in obtaining accurate
solutions to ODEs. However, when used in inference problems, they can convert
a smooth likelihood surface into a rough one, characterized by rapid (and
entirely phantom) changes in the likelihood which interfere with inference
algorithms. These inaccuracies in the likelihood can be observed even at
tolerances in the solution error where further refinements do not visibly
influence the solution. For example, in cardiac electrophysiology, jagged
parameter likelihoods have been observed with adaptive step size ODE solvers
with tolerances as low $10^{-7}$ [Johnstone, 2018, Mirams, 2018]. Next, we
investigate the origin of the jagged likelihoods using synthetic data from the
oscillator model described in eq. (8).

####  3.3.1 Inference for the damped, driven oscillator using an adaptive step
size solver

We first study the effects of adaptive time step solvers on inference using
the model system that was introduced at the beginning of §3 (eq. (8)). Here,
we set the input stimulus according to

$F(t)=\begin{cases}1,&t<t_{\text{change}},\\\ f_{1},&t\geq
t_{\text{change}}.\end{cases}$ (9)

Thus, $f_{1}$ controls the strength of a pulse provided to the system at
$t=t_{\text{change}}$.

Figure 3: Damped oscillator inference using adaptive time step Runge-Kutta.
(Left) Synthetic data for the damped driven oscillator. The curved line
indicates the accurate solution to the ODE with these parameters, while the
points indicate the noisy data. (Center) Solution for oscillator computed
using an RK5(4) solver with three different choices for the relative tolerance
(indicated by tol in the legend). (Right) Log-likelihood for the parameter $k$
calculated from the noisy data, with all other parameters held at their true
values. The log-likelihood was calculated from eq. (5) using an RK5(4) solver
with three different choices for the tolerance.

First, we consider the problem where $f_{1}=-1$ and $t_{\text{change}}=2.5$
for different choices of the RK5(4) solver tolerance. $25$ evenly spaced data
points were generated from and including $t=0$ to $t=5$ from the model with an
accurate solver (the RK5(4) solver with relative tolerance set to $10^{-8}$),
using true parameter values $k=1$, $c=0.2$, $m=1$ and an initial condition of
$x(t=0)=0$, $\dot{x}(t=0)=0$. Then, IID Gaussian noise was added to the
solution at each of the sampled locations with $\sigma=0.1$. Holding all other
parameters fixed at their true values, the log-likelihood was calculated for a
range of values of $k$, using the RK5(4) solver with various tolerances. These
results are shown in Figure 3. At insufficient tolerances, the log-likelihood
surface exhibits significant erroneous jaggedness. Notably, visual changes
between the forward simulations are minor even at tolerances which cause
drastic differences in the log-likelihood.

Figure 4: Damped oscillator model: forward simulations and inference using an
adaptive solver. Time series data and parameter likelihood surfaces are shown
for four values of $f_{1}$ in the oscillator problem: eqs. (8) & (9). For each
value of $f_{1}$, the top plot shows the accurate ODE solution (line) and
noisy synthetic data (points) generated from it. The bottom plot panels show
the corresponding log-likelihood surface for $k$ over an interval centred on
the true value, $k=1$, while all other parameters are held at their true
values. For generating the likelihood surfaces, an RK5(4) solver was used with
$\texttt{rtol}=10^{-3}$.

Next, we fix the adaptive solver tolerance and study how introducing more
rapid changes in the system’s behavior affects the log-likelihood surface. In
Figure 4, we fix $t_{\text{change}}=25$ and consider four different values of
$f_{1}$ and plot the likelihood surface for the model parameter $m$ calculated
according to an RK5(4) solver with $\texttt{rtol}=10^{-3}$. For each value of
$f_{1}$, 75 evenly spaced data points were generated on the interval from and
including $t=0$ to $t=50$, using parameter values $k=1$, $c=0.2$, and $m=1$.
IID Gaussian noise was added to the solution at each of the sampled locations
with $\sigma=0.01$. The likelihood was then calculated over a range of values
of $k$, with all other parameters held at their correct values. For $f_{1}=1$,
the stimulus $F(t)$ is constant over time, and the likelihood surface appears
smooth. However, as $f_{1}$ is adjusted so the stimulus is a stronger pulse,
the likelihood becomes jagged with large deviations away from the true
likelihood surface. (This is an example of a challenging RHS which could be
made more tractable for inference using smoothing approximations, which we
analyze in §S3.) Overall, these results indicate that the more rapid the
changes in a system’s behavior, generally the tighter solver tolerances are
required to solve the inverse problem.

Figure 5: Damped oscillator model: likelihood discontinuities caused by
variation in the number of adaptive steps. The log-likelihood surface for the
parameter $5$ in the oscillator problem (black solid line) and the number of
time points used by the adaptive step size ODE solver in the calculation of
each value of the likelihood (blue dashed line) are shown. An RK5(4) solver
was used with $\texttt{rtol}=10^{-3}$.

A fundamental point to note is that these inaccuracies arise because different
values of the parameters represent different forward problems, and the solver
selects a different sequence of step sizes for each. When the solution
contains regions of rapid change, differences in the positions of the solver
time steps, and, particularly, the inevitably discontinuous jumps in the total
number of time steps used by the solver, cause errors in the likelihood. This
phenomenon is investigated more closely in Figure 5. For this study, the
oscillator model eq. (8) was again used. 50 evenly spaced data points were
generated on the time interval from and including $t=0$ to $t=10$, with
$t_{\text{change}}=5$ and $f_{1}=-5$, using parameter values $m=1$, $c=0.2$,
and $k=1$. IID Gaussian noise was added to the solution at each of the sampled
locations with $\sigma=0.01$. The likelihood for $k$ was calculated as before
and is plotted in Figure 5. In this case, the figure is restricted to a very
narrow range of $k$ values, and the total number of time points selected by
the adaptive solver for the calculation of the likelihood at each value of $k$
is overlaid on the plot. Here, the large jumps in the likelihood correspond to
the addition or removal of a solver time point. Smaller spikes and jaggedness
where the total number of solver time points is constant correspond to
shifting of the solver time points.

####  3.3.2 Effect of jaggedness on inference algorithms

The jagged spikes appearing in the likelihood surface as a result of
insufficiently accurate adaptive step size solvers plague computational
inference algorithms. A common approach to Bayesian inference is to use the
Metropolis MCMC algorithm, or variants of it [Gelman et al., 2013]. This
algorithm generates a sequence of parameter values via a Markov chain whose
stationary distribution is the posterior distribution of the parameters. Given
the most recent parameter values in the chain $\theta^{\text{old}}$, the
Metropolis algorithm proposes new parameter values $\theta^{\text{prop}}$
according to a proposal distribution and then accepts $\theta^{\text{prop}}$
with a probability of:

$r=\min\left(1,\frac{p(\theta^{\text{prop}})}{p(\theta^{\text{old}})}\frac{p(y|\theta^{\text{prop}})}{p(y|\theta^{\text{old}})}\right),$

where $p(\theta^{\text{prop}})$ is the prior and $p(y|\theta^{\text{prop}})$
is the likelihood. To illustrate the detrimental effects of jagged errors in
the likelihood, we consider a situation where $\theta^{\text{old}}$ and
$\theta^{\text{prop}}$ have identical values under the prior and the
accurately computed likelihood (this is a plausible assumption when
$\theta^{\text{old}}$ and $\theta^{\text{prop}}$ are nearby), but we assume
that the log-likelihoods at these two parameter values computed using the
numerical approximation differ by some factor $c$ driven by numerical error in
the adaptive step size solver (i.e., $\log p(y|\theta^{\text{prop}})=\log
p(y|\theta^{\text{old}})-c$, for $c>0$). This assumption of a jump in computed
likelihood values at nearby parameter values is analogous to the spikes
appearing in the log-likelihood in our results in Figures 4 and 5.

Under these assumptions, $\log r=-c$ or $r=\exp(-c)$. For a value of $c=10$
(smaller than many of the magnitudes of spikes observed in our results), the
probability of accepting the proposal is less than 1 in 20,000. Even a
relatively small jump of magnitude $c=3$ will be traversed by the sampler with
a probability of only about 5%. Although these computations are based on
simplistic assumptions, they suggest that even minor warping of the log-
likelihood may severely restrict the ability of a Metropolis-Hastings sampler
(or similar inference algorithm) to traverse the parameter space efficiently.

###  3.4 The impact of observation error magnitude on inference and sampling
performance

In this section, we empirically study the effects of different levels of
observation noise on inference. We performed Bayesian inference using MCMC for
the oscillator problem with varying levels of noise in the data. We considered
two values of $\sigma$ ($0.01$ and $0.1$) to generate the data, fixed
$f_{1}=-1$, and otherwise generated data exactly as described for Figure 4. We
set a uniform prior on $[0.1,1.5]$ for the three model parameters $m$, $c$,
and $k$, and a uniform prior on $[0,1]$ for the $\sigma$. Three MCMC chains
were run, initialized at random samples from the prior (with the same MCMC
starting point being used for both choices of the true $\sigma$). 1500
iterations of MCMC were performed using the Haario-Bardenet adaptive
covariance algorithm as implemented in PINTS to sample from the posterior
[Haario et al., 2001, Clerx et al., 2018]. The MCMC chains for the $m$
parameter are plotted in the left column of Figure 6 using the RK5(4) solver
with $\texttt{rtol}=10^{-3}$, while the right column of Figure 6 shows the
chains using the same solver but with more accurate tolerances of
$\texttt{rtol}=10^{-8}$.

At the lowest noise level considered ($\sigma=0.01$), the three MCMC chains
using the less accurate solver move towards the true value of the parameter
but fail to mix with each other. Instead, each chain remains stuck in a narrow
region of parameter space near the true parameter value, corresponding to the
phantom local maxima in the likelihood surface observed in Figure 4. Reducing
the solver tolerance to $10^{-8}$ was, however, sufficient to ensure chain
mixing, indicating that the lack of convergence was purely an artefact of
using an inaccurate solver. At the highest level of noise considered here
($\sigma=0.1$), the three MCMC chains mix well for either tolerance
choice,111We note that, for this level of noise, the centers of the sampling
distributions are shifted slightly away from the true parameter value because
the noise limits our ability to estimate this parameter. which can be
explained by our bound given in eq. (7): that larger $\sigma$ values lead to
gentler variation in the log-likelihood surface and so easier exploration by
inference algorithms.

Figure 6: Effect of noise on MCMC convergence. Data were generated according
to the same specifications as for Figure 4, with $f_{1}=-1$, and the indicated
values of $\sigma$. Inference was performed for the three parameters $m$, $c$,
and $k$, as well as $\sigma$, via adaptive covariance MCMC [Haario et al.,
2001] with three independent chains initialized at random samples from the
prior (uniform on $[0.1,1.5]$ for the model parameters, and uniform on $[0,1]$
for $\sigma$). 1500 MCMC iterations were performed. The plots show the three
chains for the $m$ parameter. (Left) Forward simulation was performed using
the RK5(4) solver with $\texttt{rtol}=10^{-3}$. (Right) Forward simulation was
performed using the RK5(4) solver with $\texttt{rtol}=10^{-8}$.

##  4 Fixed step solvers applied to an SIR change point model of the spread of
COVID-19 in Germany

A widely used class of differential equation models in epidemiology are
compartmental models, which divide the population into a number of
compartments representing different diseased or non-diseased states and
specify the rates at which individuals move from one compartment to another
[van der Vegt et al., 2022]. A simple yet commonly used example is the SIR
model (susceptible-infected-recovered) [Weiss, 2013]. This model keeps track
of the number of susceptible individuals $S(t)$ (those who can be infected
with the disease), infected individuals $I(t)$ (those who are currently
infectious with the disease), and recovered individuals $R(t)$ (those who have
recovered from the disease and are assumed immune). Neglecting births and
deaths, the model is expressed by the following system of differential
equations:

$\frac{dS}{dt}=-\lambda\frac{SI}{N}$ (10)
$\frac{dI}{dt}=\lambda\frac{SI}{N}-\mu I$ (11) $\frac{dR}{dt}=\mu I,$ (12)

where $\lambda>0$ is the spreading rate of the disease, $\mu>0$ is the
recovery rate, and $N>0$ is the total size of the population. The system
additionally requires the specification of initial conditions for each
compartment ($S(0)$, $I(0)$, $R(0$)). $I(0)$ must exceed zero for an infection
to spread.

The qualitative behaviour of the SIR model can be determined by the basic
reproduction number, $R_{0}$, where

$R_{0}=\frac{\lambda}{\mu}.$

Assuming that $S(0)\approx N$ and $I(0)>0$, when $R_{0}>1$, the number of
infected individuals will grow, and the epidemic will eventually infect the
entire population (barring a change in $\lambda$ or $\mu$); for $R_{0}<1$, the
number of infected individuals will fall. Thus, fitting an SIR model to
infection data, and estimating the spreading rate $\lambda$ and reproduction
number $R_{0}$, are important steps in understanding and predicting the
progression of an epidemic.

An extension to the standard SIR model has $\lambda$ vary over time, allowing
the model to capture changes in the spread of a disease caused by behavioural
changes or government policy. In the aftermath of the outbreak of COVID-19 in
Europe in early 2020, an SIR model allowing changes in $\lambda$ through time
was used in a high profile paper which attempted to capture the impact of
major public health policy interventions on COVID-19 transmission in Germany
[Dehning et al., 2020]. The authors used the model eqs. (10)–(12), discretised
with a one day time step, equivalent to a Forward Euler solver with $\Delta
t=1$:

$S_{t}=S_{t-1}-\lambda(t)\Delta t\frac{S_{t-1}I_{t-1}}{N}$ (13)
$I_{t}=I_{t-1}+\lambda(t)\Delta t\frac{S_{t-1}I_{t-1}}{N}-\mu\Delta tI_{t-1}$
(14) $R_{t}=R_{t-1}+\mu\Delta tI_{t-1}.$ (15)

The initial condition was given by an unknown parameter $I_{0}=I(0)$. The
system was closed with $R(0)=0$ and $S(0)=N-I_{0}$. The spreading rate
$\lambda$ was assumed to be a continuous function of time and was allowed to
shift at three time points, whose locations were estimated from the data.
Specifically, these three time points, $t_{i},i\in\\{1,2,3\\}$ denoted the
times at which $\lambda$ began to (linearly) change to a new, constant value,
and the time taken for these shifts was dictated by durations $d_{i}$. The
$\lambda$ profile then has the following piecewise representation:

$\lambda(t)=\begin{cases}\lambda_{0},&t<t_{1},\\\
\lambda_{0}+\frac{\lambda_{1}-\lambda_{0}}{d_{1}}(t-t_{1}),&t_{1}\leq
t<t_{1}+d_{1},\\\ \lambda_{1},&t_{1}+d_{1}\leq t<t_{2},\\\
\lambda_{1}+\frac{\lambda_{2}-\lambda_{1}}{d_{2}}(t-t_{2}),&t_{2}\leq
t<t_{2}+d_{2},\\\ \lambda_{2},&t_{2}+d_{2}\leq t<t_{3},\\\
\lambda_{2},+\frac{\lambda_{3}-\lambda_{2}}{d_{3}}(t-t_{3})&t_{3}\leq
t<t_{3}+d_{3},\\\ \lambda_{3},&t_{3}+d_{3}\leq t.\end{cases}$

Additional features of the model included a reporting delay and a weekly
modulation. The reporting delay was characterised by a single parameter $D$
indicating the number of days between the time at which new infections occur
and the time at which they are reported. The modulation accounts for the
weekly periodicity evident in the data and is characterised by two parameters
$f_{w}$ and $\Phi_{w}$. This significant periodicity likely arises from
processes involved in the reporting of COVID-19 cases and deaths [Gallagher et
al., 2023]. Specifically, cases $C_{t}$ are modelled by:

$C_{t}=(1-f(t))I^{\text{new}}_{t-D},$ (16)

where

$f(t)=(1-f_{w})\left(1-\left|\sin\left(\frac{\pi}{7}t-\frac{1}{2}\Phi_{w}\right)\right|\right),$
(17)

where $I^{\text{new}}_{t}=S_{t-1}-S_{t}$. [Dehning et al., 2020] assumed a
Student-t distribution with four degrees of freedom and multiplicative noise
for the likelihood, such that the likelihood for observed cases $\hat{C_{t}}$
was given by:

$p(\hat{C_{t}}|\theta,\sigma)=\text{Student-t}_{\nu=4}(\text{mean}=C_{t}(\theta),\text{scale}=\sigma\sqrt{C_{t}(\theta)}),$

where
$\theta=(\lambda_{0},\lambda_{1},\lambda_{2},\lambda_{3},t_{1},t_{2},t_{3},d_{1},d_{2},d_{3},\mu,D,I_{0},f_{w},\Phi_{w},\sigma)$
is the full vector of parameters for the differential equation model, and
$C_{t}(\theta)$ is the deterministic solution which may be computed using a
range of different time steps. The prior distributions for the parameters are
given in Table S1 (supplementary information).

###  4.1 Effect of time step on the forward solution

We first study the effect of assuming $\Delta t=1\text{ day}$ on forward
simulations of the model. We set up the forward simulations using the same
settings that [Dehning et al., 2020] used to generate their Figure 2. The
parameters of an SIR model without change points or weekly modulation (i.e., a
single value of $\lambda$, $\mu$, $D$ $I_{0}$, and $\sigma$) were inferred
from an early period of the German daily reported COVID-19 cases, from 2 March
2020 to 15 March 2020. The posterior median values of these parameters
(excepting $\lambda$) were then used to generate forward simulations according
to the full model without weekly modulation (eqs. (13)–(16)), with one change
point, and pre-specified values of $\lambda_{0}$ and $\lambda_{1}$.

As in [Dehning et al., 2020], the first set of simulations considered how
different levels of social restrictions could influence the course of disease
transmission, as measured by cases. Three levels of social restrictions
(assumed to be captured by different $\lambda$ values) are considered, which
each yield two sets of simulations: one corresponding to Forward Euler with
$\Delta t=1\text{ day}$ (as in [Dehning et al., 2020]) and another with
$\Delta t=0.1$ days. We choose $0.1$ days as the more accurate comparator
method, as further refinement of the step size yields little change in forward
solutions or inference results but is increasingly costly to run. The results
of this are shown in Figure 7A. Our second set of simulations, shown in Figure
7B, considered only our “strong” social distancing scenario and explored three
different times at which the change in $\lambda$ might occur (e.g., if a
public health intervention were implemented at different times). These
simulations illustrate how using a large time step generally leads to a
substantial underestimation of case counts for a given choice of $\lambda(t)$,
particularly during the (crucial) growth phase of the epidemic.

Figure 7: COVID-19 model: forward simulations using Forward Euler. In both (a)
and (b), the top panel shows three different pre-specified trajectories of
$\lambda(t)$, and the bottom panel shows the number of daily cases resulting
from these trajectories for each choice of the time step $\Delta t$. Figure 8:
COVID-19 model: inference using Forward Euler (Left) Real data and model fits
for the number of daily COVID-19 cases in Germany over the period 2 March 2020
to 21 April 2020. Note that the model fits for $\Delta t=1$ and $\Delta t=0.1$
overlap almost completely. (Right) Inferred basic reproduction number over
time for the Germany COVID-19 data, using the SIR model with change points in
$\lambda$ (eqs. (13)–(17)) and two different values for the ODE solver time
step, $\Delta t$. In both panels, lines indicate the posterior median and
shaded regions indicate the central $95\%$ of the posterior.

###  4.2 Effect of time step on the posterior distributions

We also studied the effect of the time step on parameter inference for the
full model (eqs. (13)–(17)) using the German daily cases data from 2 March
2020 to 21 April 2020 as was done in [Dehning et al., 2020]. Inference was
performed using the PyMC3 No-U-Turn MCMC Sampler (NUTS) [Salvatier et al.,
2016, Gelman et al., 2013] using the model developed by [Dehning et al.,
2020], modified to allow the $0.1$ day step size. To initialize the chains,
automatic differentiation variational inference [Kucukelbir et al., 2017] as
implemented in PyMC3 [Salvatier et al., 2016] was performed to generate an
approximate posterior (which, however, does not capture correlations between
the parameters). Four MCMC chains were then initialized by sampling from this
approximation of the posterior. The chains were run for $500$ iterations of
NUTS, with the first $100$ discarded as burn-in, and convergence assessed by
requiring that $\hat{R}<1.05$ [Gelman et al., 2013]. These results are shown
in Figure 8.

Both models achieve a near identical visual fit to the data, using the median
values of the recovered parameters. However, the parameter estimates of the
two models differed. We focus on the posterior distribution for the basic
reproduction number $R_{0}$, which is calculated using the MCMC samples of the
joint posterior for $(\lambda,\mu)$. The one day time step results in
overestimation of $R_{0}$ (by approximately 10% relative to the $0.1$ day time
step) during the early stages of the epidemic (i.e., before the first change
point). This is because, during the growth phase of the epidemic, the larger
time step results in slower growth for a given $\lambda$ value (cf. Figure 7),
meaning a larger $\lambda$ value is estimated to compensate. During the later
stages of the epidemic, the values of $R_{0}$ are more similar between the two
models. Additionally, the change point locations are not much affected by the
choice of time step (though, this is expected as the change points have fairly
informative priors).

Our results indicate that while the discrete version of the SIR change point
model using $\Delta t=1$ appears visually to obtain a good fit to German
COVID-19 data, the growth parameters of the discrete model using this time
step vary markedly from those recovered using $\Delta t=0.1$, and thus care
should be taken in the deployment of such discrete models and the reporting of
their results.

##  5 Numerical errors arising in rainfall-runoff models of river streamflow
data

In this section we use real data from the French Broad River at Asheville,
North Carolina to investigate the impact of adaptive solvers in performing
inference for rainfall-runoff models used in hydrology [Schoups and Vrugt,
2010, Schoups et al., 2010].

Rainfall-runoff models divide the flow of water through a river basin into
several spatially grouped components representing different hydrological
processes. The model we consider here is governed by a system of five ODEs:

$\frac{dS_{i}}{dt}=\text{Precip}(t)-\text{InterceptEvap}(t)-\text{EffectPrecip}(t)$
(18)
$\frac{dS_{u}}{dt}=\text{EffectPrecip}(t)-\text{UnsatEvap}(t)-\text{Percolation}(t)-\text{Runoff}(t)$
(19) $\frac{dS_{s}}{dt}=\text{Percolation}(t)-\text{SlowStream}(t)$ (20)
$\frac{dS_{f}}{dt}=\text{Runoff}(t)-\text{FastStream}(t)$ (21)
$\frac{dz}{dt}=\text{SlowStream}(t)+\text{FastStream}(t),$ (22)

Each term in this equation is defined in Table S2 (supplementary information),
and the seven unknown parameters of the model and their prior distributions
are defined in Table S3 (supplementary information). The data consist of daily
streamflow measurements ($dz/dt$), and the authors assume an IID Gaussian
likelihood with unknown standard deviation $\sigma$.

Figure 9: Rainfall run-off model: inference using accurate and inaccurate
adaptive solvers. Here, we plot the likelihood surface for each parameter for
the rainfall-runoff model defined by eqs. (18)–(22). For each parameter, the
solid line indicates the likelihood calculated using an RK3(2) adaptive solver
with $\texttt{rtol}=\texttt{atol}=10^{-3}$, while the dashed line indicates
the likelihood calculated using the CVODE adaptive solver with
$\texttt{rtol}=\texttt{atol}=10^{-7}$.

Previous work has shown that using large time steps with such hydrological
models can bias inferences [Kavetski et al., 2003]. We show that using an
adaptive step size method (as suggested by [Schoups et al., 2010]) can also
cause inaccurate inference results, unless the error is tightly controlled.

Using an accurate ODE solver (the CVODE multistep solver from the SUNDIALS
library [Hindmarsh et al., 2005] with $\texttt{rtol}=\texttt{atol}=10^{-7}$),
we obtained the posterior distributions for the seven parameters of the model,
using USGS data for the streamflow at Asheville, North Carolina (USGS station
03451500) over a 200 day period starting 1 January 1960. Sampling was
performed using the Dream multi-chain MCMC algorithm as implemented in PINTS
[Vrugt et al., 2009, Clerx et al., 2018], using 6 chains with each initialized
by a sample from the prior (Table S3, supplementary information). $25000$ MCMC
iterations were performed, and convergence of the chains was assessed by
requiring that $\hat{R}<1.05$ [Gelman et al., 2013]. In Figure 9, we plot the
likelihood surfaces of the parameters for slices through parameter space near
the estimated posterior medians. Likelihood surfaces are plotted for two
adaptive step size solvers: the RK3(2) solver from SciPy with
$\texttt{rtol}=\texttt{atol}=10^{-3}$, and the CVODE solver as described
above. For all parameters, the $10^{-3}$ tolerance solver causes highly jagged
likelihoods, of sufficient magnitude to interfere with inference via MCMC or
maximum likelihood estimation. This is in accordance with our earlier results
using the oscillator model in §3, as rapid changes in the RHS cause spurious
jaggedness in the computed likelihood. The likelihoods calculated using the
more accurate solver have similar broadscale shapes but are smooth enough for
accurate inference to be performed.

##  6 Discussion

Inaccurate solution of ODEs through either fixed time step or adaptive solvers
can lead to biased inferences which are generally exacerbated when there is
low observation noise. For adaptive solvers, these biases may manifest through
the presence of phantom jaggedness in the likelihood surface, which can wreak
havoc for inference algorithms attempting to navigate the surface. They may
also lead researchers to falsely conclude that a model is unidentifiable,
when, for example, the chains in an MCMC run fail to mix. They may then choose
to modify the model in arbitrary ways when, in fact, all that was required to
render inference soluble was a reduction in solver tolerances. Tolerances
which seem good enough for forward simulation are likely insufficient for
solving inverse problems. For example, a default relative tolerance of
$10^{-3}$ was insufficient for both the synthetic data and real data studied
in §3.3 and §5. When using an ODE solver library to perform inference, default
settings may well not suffice and, ideally, the solver tolerance should be set
by inspection of the likelihood surface.

Unless there is a bifurcation in system behavior at points in parameter space,
likelihood surfaces should not have abrupt discontinuities. So, the presence
of such changes may well be an artefact of using an adaptive ODE solver with
insufficient tolerances. MCMC and optimisation algorithms could be augmented
by monitoring for such jumps and warning the user should they occur.

ODEs involving discontinuous RHS functions are known to be particularly
challenging to solve accurately. Our results indicate that RHS functions
involving rapid changes over time, such as those involving discontinuities,
also curse computational inference when adaptive ODE solvers are used.
However, our results in §S3 also indicate that errors in the likelihood
arising from discontinuous RHS functions can be ameliorated through the use of
simple smoothing approximations—a potentially more computationally efficient
alternative to increasing tolerances. We argue that in many scientific systems
such smoothing approximations are additionally more realistic descriptions of
the phenomena being modelled. Although the appropriate degree of smoothing may
be difficult to determine in general, for certain systems, the level of
smoothing can be tuned based on knowledge of the process being modelled.

Much of the work on error control for ODE solvers has focused merely on the
accuracy of the forward problem. The accuracies of widely used ODE solvers are
typically tuned via their step sizes or local truncation error tolerances, but
these are not the most relevant quantities for inference—instead, it is the
error in the log-likelihood which must be controlled. ODE solvers which
control the error on the log-likelihood directly would avoid much of thie
problems highlighted in this paper, and we suggest this as a fruitful resarch
direction.

##  7 End section

###  7.1 Data, code and materials

The code to perform the computer experiments presented in this paper was
written in Python 3.7 and is available in an open source Python library at
https://github.com/rccreswell/ode_inference. To run the COVID-19 simulations,
we adapted the software library developed by [Dehning et al., 2020]. The
version of the code including our modifications is available at
https://github.com/rccreswell/covid19_inference_forecast.

###  7.2 Authors’ contributions

R.C.: conceptualization, formal analysis, investigation, methodology,
software, validation, visualization, writing—original draft, writing—review
and editing; K.M.S.: formal analysis, investigation, methodology, software,
visualization, writing—original draft, writing—review and editing; B.L.:
conceptualization, formal analysis, investigation, methodology, project
administration, supervision, validation, writing—original draft,
writing—review and editing; G.R.M: conceptualization, investigation,
methodology, writing—review and editing; C.L.L.: conceptualization,
investigation, methodology, supervision, writing—original draft,
writing—review and editing; S.T.: formal analysis, investigation, methodology,
writing—original draft, writing—review and editing; M.R.: conceptualization,
investigation, supervision, methodology, writing—original draft,
writing—review and editing; D.J.G.: conceptualization, formal analysis,
investigation, methodology, project administration, supervision,
writing—original draft, writing—review and editing.

###  7.3 Conflict of interest declaration

We declare that we have no competing interests.

###  7.4 Funding

R.C. acknowledges support from the EPSRC via a doctoral training partnership
studentship in the Department of Computer Science at the University of Oxford.
K.M.S. and D.J.G. acknowledge funding from the EPSRC CDT in Sustainable
Approaches to Biomedical Science: Responsible and Reproducible Research -
SABS:R3 (EP/S024093/1). C.L.L. acknowledges support from the Science and
Technology Development Fund, Macao SAR (FDCT) (reference number 0048/2022/A)
and support from the University of Macau via a UM Macao Fellowship.

## References

  * [Birrell et al., 2021] Birrell, P., Blake, J., Van Leeuwen, E., Gent, N., and De Angelis, D. (2021).  Real-time nowcasting and forecasting of covid-19 dynamics in england: the first wave.  Philosophical Transactions of the Royal Society B, 376(1829):20200279.
  * [Clerx et al., 2018] Clerx, M., Robinson, M., Lambert, B., Lei, C. L., Ghosh, S., Mirams, G. R., and Gavaghan, D. J. (2018).  Probabilistic inference on noisy time series (pints).  arXiv preprint arXiv:1812.07388.
  * [Dehning et al., 2020] Dehning, J., Zierenberg, J., Spitzner, F. P., Wibral, M., Neto, J. P., Wilczek, M., and Priesemann, V. (2020).  Inferring change points in the spread of covid-19 reveals the effectiveness of interventions.  Science, 369(6500).
  * [Dormand and Prince, 1980] Dormand, J. R. and Prince, P. J. (1980).  A family of embedded Runge-Kutta formulae.  Journal of Computational and Applied Mathematics, 6(1):19–26.
  * [Gallagher et al., 2023] Gallagher, K., Creswell, R., Gavaghan, D., and Lambert, B. (2023).  Identification and attribution of weekly periodic biases in epidemiological time series data.  medRxiv.
  * [Gautschi, 1997] Gautschi, W. (1997).  Numerical analysis.  Springer Science & Business Media.
  * [Gelman et al., 2013] Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian data analysis.  CRC press.
  * [Haario et al., 2001] Haario, H., Saksman, E., and Tamminen, J. (2001).  An adaptive metropolis algorithm.  Bernoulli, pages 223–242.
  * [Hindmarsh et al., 2005] Hindmarsh, A. C., Brown, P. N., Grant, K. E., Lee, S. L., Serban, R., Shumaker, D. E., and Woodward, C. S. (2005).  Sundials: Suite of nonlinear and differential/algebraic equation solvers.  ACM Transactions on Mathematical Software (TOMS), 31(3):363–396.
  * [Johnstone, 2018] Johnstone, R. H. (2018).  Uncertainty characterisation in action potential modelling for cardiac drug safety.  PhD thesis, University of Oxford.
  * [Kavetski et al., 2003] Kavetski, D., Kuczera, G., and Franks, S. W. (2003).  Semidistributed hydrological modeling: A “saturation path” perspective on topmodel and vic.  Water Resources Research, 39(9).
  * [Kucukelbir et al., 2017] Kucukelbir, A., Tran, D., Ranganath, R., Gelman, A., and Blei, D. M. (2017).  Automatic differentiation variational inference.  Journal of Machine Learning Research.
  * [Lambert et al., 2023] Lambert, B., Lei, C. L., Robinson, M., Clerx, M., Creswell, R., Ghosh, S., Tavener, S., and Gavaghan, D. (2023).  Autocorrelated measurement processes and inference for ordinary differential equation models of biological systems.  Journal of the Royal Society Interface, 20.
  * [Mirams, 2018] Mirams, G. (2018).  Numerical errors from ODE solvers can mess up optimisation and inference very easily.  Mathematical Matters of the Heart.
  * [Salvatier et al., 2016] Salvatier, J., Wiecki, T. V., and Fonnesbeck, C. (2016).  Probabilistic programming in Python using PyMC3.  PeerJ Computer Science, 2:e55.
  * [Schoups et al., 2010] Schoups, G., Vrugt, J., Fenicia, F., and Van De Giesen, N. (2010).  Corruption of accuracy and efficiency of Markov chain Monte Carlo simulation by inaccurate numerical implementation of conceptual hydrologic models.  Water Resources Research, 46(10).
  * [Schoups and Vrugt, 2010] Schoups, G. and Vrugt, J. A. (2010).  A formal likelihood function for parameter and predictive inference of hydrologic models with correlated, heteroscedastic, and non-Gaussian errors.  Water Resources Research, 46(10).
  * [Shertzer et al., 2002] Shertzer, K. W., Ellner, S. P., Fussmann, G. F., and Hairston Jr, N. G. (2002).  Predator-prey cycles in an aquatic microcosm: testing hypotheses of mechanism.  Journal of Animal Ecology, pages 802–815.
  * [van der Vegt et al., 2022] van der Vegt, S. A., Dai, L., Bouros, I., Farm, H. J., Creswell, R., Dimdore-Miles, O., Cazimoglu, I., Bajaj, S., Hopkins, L., Seiferth, D., Cooper, F., Lok, C. L., Gavaghan, D., and Lambert, B. (2022).  Learning transmission dynamics modelling of COVID-19 using comomodels.  Mathematical Biosciences, 349:108824.
  * [Virtanen et al., 2020] Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E., Peterson, P., Weckesser, W., Bright, J., van der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, E., Carey, C. J., Polat, İ., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors (2020).  SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python.  Nature Methods, 17:261–272.
  * [Vrugt et al., 2009] Vrugt, J. A., Ter Braak, C., Diks, C., Robinson, B. A., Hyman, J. M., and Higdon, D. (2009).  Accelerating Markov chain Monte Carlo simulation by differential evolution with self-adaptive randomized subspace sampling.  International Journal of Nonlinear Sciences and Numerical Simulation, 10(3):273–290.
  * [Weiss, 2013] Weiss, H. H. (2013).  The SIR model and the foundations of public health.  Materials Mathematics, pages 1–17.
  * [Whittaker et al., 2020] Whittaker, D. G., Clerx, M., Lei, C. L., Christini, D. J., and Mirams, G. R. (2020).  Calibration of ionic and cellular cardiac electrophysiology models.  Wiley Interdisciplinary Reviews: Systems Biology and Medicine, 12(4):e1482.

##  S1 Proof of bound on the error in the log-likelihood

As stated in the main text, we denote the log-likelihood computed using the
numerical approximation to the forward map by:

$\mathcal{L}^{\prime}=-\frac{N}{2}\log(2\pi)-\frac{N}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-g(\hat{x}_{i}))^{2},$

while the log-likelihood computed using the hypothetical, true solution to the
ODE is:

$\mathcal{L}=-\frac{N}{2}\log(2\pi)-\frac{N}{2}\log(\sigma^{2})-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-g(x(t_{i};\theta)))^{2},$

where $x(t;\theta)$ is the true solution to the ODE,
$\\{\hat{x}_{i}\\}_{i=1}^{N}$ is the approximate solution to the ODE at the
data time points $\\{t_{i}\\}_{i=1}^{N}$, and $g$ is the observation operator.

For brevity, let $g_{i}=g(x(t_{i};\theta))$ and $\hat{g}_{i}=g(\hat{x}_{i})$.
We have:

$\displaystyle\mathcal{L}-\mathcal{L}^{\prime}$
$\displaystyle=-\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-g_{i})^{2}+\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(y_{i}-\hat{g}_{i})^{2}$
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left[-g_{i}^{2}+\hat{g}_{i}^{2}+2y_{i}(g_{i}-\hat{g}_{i})\right]$
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}(g_{i}-\hat{g}_{i})\left(2(y_{i}-g_{i})+(g_{i}-\hat{g}_{i})\right)$
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left((g_{i}-\hat{g}_{i})^{2}+2(g_{i}-\hat{g}_{i})(y_{i}-g_{i})\right).$
(23)

Taking the absolute value,

$\displaystyle|\mathcal{L}-\mathcal{L}^{\prime}|$
$\displaystyle=\left|\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left((g_{i}-\hat{g}_{i})^{2}+2(g_{i}-\hat{g}_{i})(y_{i}-g_{i})\right)\right|$
$\displaystyle\leq\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left|(g_{i}-\hat{g}_{i})^{2}+2(g_{i}-\hat{g}_{i})(y_{i}-g_{i})\right|$
$\displaystyle\leq\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left((g_{i}-\hat{g}_{i})^{2}+2\left|(g_{i}-\hat{g}_{i})(y_{i}-g_{i})\right|\right)$
(24)

To proceed further, we impose the assumption of Lipschitz continuity of the
observation function $g$ with Lipschitz constant $K$, i.e.,
$|g(x_{1})-g(x_{2})|\leq K|x_{1}-x_{2}|$ for all
$x_{1},x_{2}\in\mathbb{R}^{l}$. We thus bound:

$\displaystyle|g_{i}-\hat{g}_{i}|=|g(x(t_{i};\theta))-g(\hat{x_{i}})|$
$\displaystyle\leq K|x(t_{i};\theta)-\hat{x_{i}}|$
$\displaystyle=K|e(t_{i})|.$ (25)

Using this in eq. (24),

$\displaystyle|\mathcal{L}-\mathcal{L}^{\prime}|$
$\displaystyle\leq\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left(K^{2}|e(t_{i})|^{2}+2K|e(t_{i})||y_{i}-g_{i}|\right).$

As discussed further in the main text, this bound indicates an inverse
relationship between $\sigma$ and $|\mathcal{L}-\mathcal{L}^{\prime}|$ when
$e(t_{i})$ is held constant.

##  S2 Distribution of the error in the log-likelihood

In this section, rather than deriving a bound on the _absolute value_ of the
error in the log-likelihood arising from numerical errors in the forward
solution, we study the _distribution_ of the difference between the true and
numerically approximated log-likelihoods.

We assume that the $y_{i}$ are distributed according to the specified IID
Gaussian likelihood at the same parameter values at which the likelihood is
being evaluated, so we have $y_{i}\sim N(g(x(t_{i};\theta)),\sigma)$. For
brevity, as before, let $g_{i}=g(x(t_{i};\theta))$ and
$\hat{g}_{i}=g(\hat{x}_{i})$, and define $a_{i}=g_{i}-\hat{g}_{i}$. Then, we
can write $y_{i}=g_{i}+\epsilon_{i}$, where $\epsilon_{i}\sim N(0,\sigma)$.
Starting from §S1, eq. (23), we have:

$\displaystyle\mathcal{L}-\mathcal{L}^{\prime}$
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left(a_{i}^{2}+2a_{i}(y_{i}-g_{i})\right)$
$\displaystyle=\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}\left(a_{i}^{2}+2a_{i}\epsilon_{i}\right)\text{
using }y_{i}=g_{i}+\epsilon_{i}$
$\displaystyle=\frac{1}{\sigma^{2}}\sum_{i=1}^{N}a_{i}\epsilon_{i}+\frac{1}{2\sigma^{2}}\sum_{i=1}^{N}a_{i}^{2}.$
(26)

Noting that the first term in eq. (26) is the sum of independent random
variables each with mean $0$ and variance $(a_{i}^{2}/\sigma^{2})$, we obtain
the result:

$\displaystyle\mathcal{L}-\mathcal{L}^{\prime}\sim
N\left(\frac{\sum_{i=1}^{N}a_{i}^{2}}{2\sigma^{2}},\frac{\sqrt{\sum_{i=1}^{N}a_{i}^{2}}}{\sigma}\right).$

We have
$\mathbb{E}[\mathcal{L}-\mathcal{L}^{\prime}]=\frac{\sum_{i=1}^{N}a_{i}^{2}}{2\sigma^{2}}$.
Therefore, on average the numerical approximation of the log-likelihood
underestimates the true log-likelihood (when both are computed at the same
parameter values which generated the data). The average size of this
underestimation is greater when $\sigma$ is smaller. Also, in the case that
$g$ is the identity function,
$\sum_{i=1}^{N}a_{i}^{2}=\sum_{i=1}^{N}|e(t_{i})|^{2}$, and so the average
size of the underestimation decreases as the global truncation error
decreases.

##  S3 Smoothing forcing terms to reduce numerical errors in the likelihood

As indicated by our results in Figures 4 and 5, discontinuities in the right-
hand side (RHS) of an ODE can lead to substantial errors in the likelihood
when adaptive step size solvers are used. In general, errors in the likelihood
arising from numerical errors in the solution can be reduced by refining the
tolerance of the adaptive solver. However, when the RHS suffers from a
discontinuity, the required solver tolerance to obtain an acceptable
likelihood surface may employ a prohibitively large grid of solver points.
Several approaches to remove discontinuities from the RHS have been developed
to enable more accurate forward simulations, including smoothing
approximations and solving the ODE separately within regions where the RHS is
continuous [stewart2011dynamics]. These techniques may be particularly
advantageous when performing inference. In this section, we study the effects
on the computation of the log-likelihood of one of these approaches, which is
to smooth discontinuities in the RHS of the ODE. Smoothing is often a
particularly appropriate assumption for biological models, where a continuous
rather than an instant change may in fact more realistically represent the
true behaviour of the system. For example, in epidemiology, interventions
(such as the introduction of a vaccination campaign) may be naively
represented by discontinuous step functions in the RHS of a compartmental
epidemiological model; however, a function smoothly moving between two values
(corresponding to the intervention reaching its full effect gradually over an
appropriate period of time) is both more realistic and more tractable for
numerical solvers for the forward problem [van der Vegt et al., 2022].

The hyberbolic tangent function ($\tanh$) is a useful smooth approximation to
a step function. In the forced oscillator problem, we can use $\tanh$ to
approximate the step function stimulus, eq. (9) (main text), with $f_{1}=-1$
according to:

$F^{\text{smooth}}(t)=-\tanh\left(\frac{t-t_{\text{change}}}{a}\right)$ (27)

where $a$ is a tuning parameter controlling the level of smoothing, with
larger values of $a$ leading to a more gradual change in the stimulus, and
$t_{\text{change}}$ is the time when the stimulus changes in value.

Figure S1: Effect of tanh-smoothing on likelihood surface. (a.) Synthetic data
for the damped driven oscillator. The curved line indicates the accurate
solution to the ODE with these parameters, while the points indicate the noisy
data. (b.) The three considered forms of the stimulus. $a=0$ indicates the
unsmoothed stimulus (eq. (9), main text), while the positive values of $a$
indicate the tanh-smoothed stimulus according to eq. (27). (c.) Solution for
oscillator computed using an RK5(4) solver with relative tolerance $10^{-3}$,
with three different forms of the stimulus, at the true parameter values. (d.)
Log-likelihood for the parameter $k$ calculated from the noisy data, with all
other parameters held at their true values. The log-likelihood was calculated
from eq. (5) (main text) using an RK5(4) solver with relative tolerance
$10^{-3}$.

To examine the effect of the smoothing approximation on inference, we computed
the likelihood surface for the $k$ parameter in the forced oscillator model
using a variety of choices for the smoothing parameter, with results shown in
Figure S1. Using $f_{1}=-1$ and $t_{\text{change}}=2.5$, $25$ evenly spaced
data points were generated from and including $t=0$ to $t=5$ from the model
with an accurate solver (the RK5(4) solver with relative tolerance set to
$10^{-8}$), using true parameter values $k=1$, $c=0.2$, $m=1$ and an initial
condition of $x(t=0)=0$, $\dot{x}(t=0)=0$. Then, IID Gaussian noise was added
to the solution at each of the sampled locations with $\sigma=0.1$. Holding
all other parameters fixed at their true values, the log-likelihood was
calculated for a range of values of $k$, using the RK5(4) solver with relative
tolerance tuned to $10^{-3}$. The likelihood was computed using both the
original step function stimulus eq. (9) (main text) (indicated in Figure S1 by
$a=0$), as well as the smooth approximation eq. (27) with two different
choices of $a>0$. Without smoothing, we observe significant jagged biases in
the likelihood, as expected due to the insufficient solver tolerance. However,
with smoothing, a smooth, tractable likelihood surface is obtained despite the
mediocre solver tolerance. This is despite the fact that all forward
simulations are visually very similar. This is in accordance with our results
in Figure 3 (main text), where even visibly small changes in the forward
solution may hide the fact that there lurks substantial distortions of the
likelihood surface.

##  S4 Supplementary information

Parameter | Prior
---|---
$\lambda_{0}$ | $\text{log normal}(\log(0.4),0.5)$
$\lambda_{1}$ | $\text{log normal}(\log(0.2),0.5)$
$\lambda_{2}$ | $\text{log normal}(\log(0.125),0.5)$
$\lambda_{3}$ | $\text{log normal}(\log(0.0625),0.5)$
$t_{1}$ | $N(2020\text{ March }9,3\text{ days})$
$t_{2}$ | $N(2020\text{ March }16,1\text{ day})$
$t_{3}$ | $N(2020\text{ March }23,1\text{ day})$
$d_{i}$ | $\text{log normal}(\log(3),0.3)$
$\mu$ | $\text{log normal}(\log(0.0625),0.2)$
$D$ | $\text{log normal}(\log(8),0.2)$
$I_{0}$ | $\text{half Cauchy}(100)$
$f_{w}$ | $\text{beta}(0.7,0.17)$
$\Phi_{w}$ | $\text{Von-Mises}(0,0.01)$
$\sigma$ | $\text{half Cauchy}(10)$
Table S1: Prior distributions for parameters in the SIR changepoint model.

Term |  Definition |  Description
---|---|---
$S_{i}$ |  Interception storage |  Water which strikes vegetative surfaces.
$S_{u}$ |  Unsaturated storage |  Storage of water in the soil above the water table.
$S_{s}$ |  Slow reservoir |  Water moving to the river via percolation.
$S_{f}$ |  Fast reservoir |  Water moving to the river via surface runoff.
$z$ |  River discharge |  Water flowed out of the river at the measuring location.
$f(S,a)$ |  $\frac{1-e^{-aS}}{1-e^{-a}}$ |  Nonlinear flux function.
$\text{Precip}(t)$ |  Precipitation |  Areal precipitation in the river basin, provided as input to the model.
$\text{Evap}(t)$ |  Evaporation |  Evaporation from the river basin, provided as input to the model.
$\text{InterceptEvap}(t)$ |  $\text{Evap}(t)f(S_{i}/I_{\text{max}},\alpha_{i})$ |  Evaporation from interception.
$\text{EffectPrecip}(t)$ |  $\text{Precip}(t)f(S_{i}/I_{\text{max}},-\alpha_{i})$ |  Effective precipitation which reaches unsaturated storage.
$\text{UnsatEvap}(t)$ |  $\max(0,\text{Evap}(t)-\text{InterceptEvap}(t))f(S_{u}/S_{\text{u,max}},\alpha_{e})$ |  Evaporation from unsaturated storage.
$\text{Percolation}(t)$ |  $Q_{\text{s,max}}f(S_{u}/S_{\text{u,max}},\alpha_{s})$ |  Trickling of water through the ground.
$\text{Runoff}(t)$ |  $\text{EffectPrecip}(t)f(S_{u}/S_{\text{u,max}},\alpha_{f})$ |  Flow of water on the surface.
$\text{SlowStream}(t)$ |  $S_{s}/K_{s}$ |  Slow component of the river flow.
$\text{FastStream}(t)$ |  $S_{f}/K_{f}$ |  Fast component of the river flow.
Table S2: Description of the terms which appear in the rainfall-runoff model.

Parameter |  Definition |  Prior
---|---|---
$I_{\text{max}}$ |  Maximum interception |  $\text{Uniform}(0,10)$
$S_{\text{u,max}}$ |  Unsaturated storage capacity |  $\text{Uniform}(10,1000)$
$Q_{\text{s,max}}$ |  Maximum percolation |  $\text{Uniform}(0,100)$
$\alpha_{e}$ |  Evaporation flux shape |  $\text{Uniform}(0,100)$
$\alpha_{f}$ |  Runoff flux shape |  $\text{Uniform}(-10,10)$
$K_{s}$ |  Slow reservoir time constant |  $\text{Uniform}(0,150)$
$K_{F}$ |  Fast reservoir time constant |  $\text{Uniform}(0,10)$
$\alpha_{s}=0$ |  Percolation flux shape |  -
$\alpha_{i}=50$ |  Interception flux shape |  -
$\sigma$ |  Noise standard deviation |  Uniform(0, 10)
Table S3: Description of the seven unknown parameters of the model, and the
two parameters with fixed values.