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

Most ordinary differential equation (ODE) models used to describe biological or physical systems must be solved approximately using numerical methods. Perniciously, even those solvers that seem sufficiently accurate for the forward problem, i.e. for obtaining an accurate simulation, might 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 might 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 reanalyse an ODE change point 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 run-off 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.


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 loglikelihood

Log-likelihood function for an ODE model
We assume that time series data {y i } N i=1 ; y i ∈ R n are measured at time points {t i } N i=1 .These data are believed to obey some function g : R l → R n of x(t; θ) ∈ R l , where x is the solution to an ordinary differential equation: dx dt = h(t, x, θ); (1) for some function h which is informed by the relevant scientific theory and parameterized by the (potentially unknown) values θ ∈ 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: where i is a mean-zero random variable specifying the noise process.Many choices for i are possible, and the assumed distribution of this variable should be chosen carefully, as misspecification of 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: with the parameter σ representing the standard deviation of the noise process, which is also inferred from the data together with the model parameters θ.Choosing a particular noise process enables the joint probability of the observations to be expressed as a function of the parameter values θ-this is known as a likelihood.For an IID Gaussian noise process (eq.( 4)), the log-likelihood for time series data {y i } N i=1 takes the form: 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 θ most compatible with data can be found by maximising the likelihood with respect to θ (i.e., the method of maximum likelihood).Alternatively, the likelihood can be used together with a prior distribution on the parameters (p(θ)) to infer the posterior distribution according to Bayes' theorem: 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].

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 ; θ) is the true solution to the ODE at time point t i , which is unavailable and approximated by xi .The deviation between x(t i ; θ) and xi at any time point is given by the global truncation error, e(t 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 xi instead of x(t i ; θ).For brevity, we denote the accurate log-likelihood by L, and we denote the log-likelihood computed using a numerical approximation by L , which is given by Assuming Lipschitz continuity of the observation function g with Lipschitz constant K, we prove in §S1 that the difference between L and L is bounded according to: We observe an inverse relationship between σ and the bound of |L − L | 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 σ.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 σ, 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 E[L − L ] > 0: the numerical approximation of the likelihood will, on average, underestimate the true likelihood.

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: 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: where ẋ = dx/dt.

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 atol + 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 RKp(q), which use the qth order method to estimate the error (and thus control the time step), while making the actual steps using the pth 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].

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, ẋ(t = 0) = 0, and Then IID Gaussian noise was added to the solution at each of the sampled locations with σ = 0.01.Holding all other parameters fixed at their true values, the loglikelihood 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 ∆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,   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 ∆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).
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.

Forward Euler solver
One of the simplest numerical solvers for ordinary differential equations is the Forward Euler method with a uniform step size ∆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 ∆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 ∆t may alter both forward model solutions and parameter inference results.(Center) Solution for oscillator computed using a Forward Euler solver with four different choices for the time step ∆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 ∆t.

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, ẋ(t = 0) = 0, and F (t) = 1.Then, IID Gaussian noise was added to the solution at each of the sampled locations with σ = 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, ∆t = 0.1, causes substantial inaccuracy in the likelihood even though ∆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.(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.

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).

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 Thus, f 1 controls the strength of a pulse provided to the system at t = t change .
First, we consider the problem where f 1 = −1 and t 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 =  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 rtol = 10 −3 .0.2, m = 1 and an initial condition of x(t = 0) = 0, ẋ(t = 0) = 0.Then, IID Gaussian noise was added to the solution at each of the sampled locations with σ = 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 loglikelihood.
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 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 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 σ = 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.
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 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 σ = 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.

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 θ old , the Metropolis algorithm proposes new parameter values θ prop according to a proposal distribution and then accepts θ prop with a probability of: where p(θ prop ) is the prior and p(y|θ prop ) is the likelihood.To illustrate the detrimental effects of jagged errors in the likelihood, we consider a situation where θ old and θ prop have identical values under the prior and the accurately computed likelihood (this is a plausible assumption when θ old and θ prop are nearby), but we assume that the loglikelihoods 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|θ prop ) = log p(y|θ 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.

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 σ (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 σ.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 σ).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 rtol = 10 −3 , while the right column of Figure 6 shows the chains using the same solver but with more accurate tolerances of rtol = 10 −8 .At the lowest noise level considered (σ = 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 Inference was performed for the three parameters m, c, and k, as well as σ, 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 σ).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 rtol = 10 −3 .(Right) Forward simulation was performed using the RK5(4) solver with rtol = 10 −8 .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 (σ = 0.1), the three MCMC chains mix well for either tolerance choice,1 which can be explained by our bound given in eq. ( 7): that larger σ values lead to gentler variation in the log-likelihood surface and so easier exploration by inference algorithms.

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: where λ > 0 is the spreading rate of the disease, µ > 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 Assuming that S(0) ≈ 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 λ or µ); for R 0 < 1, the number of infected individuals will fall.Thus, fitting an SIR model to infection data, and estimating the spreading rate λ 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 λ 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 λ 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 ∆t = 1: 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 λ 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 ∈ {1, 2, 3} denoted the times at which λ began to (linearly) change to a new, constant value, and the time taken for these shifts was dictated by durations d i .The λ profile then has the following piecewise representation: 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 Φ 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: where where [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 Ĉt was given by: where θ = (λ 0 , λ 1 , λ 2 , λ 3 , t 1 , t 2 , t 3 , d 1 , d 2 , d 3 , µ, D, I 0 , f w , Φ w , σ) is the full vector of parameters for the differential equation model, and C t (θ) 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).

Effect of time step on the forward solution
We first study the effect of assuming ∆t = 1 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 λ, µ, D I 0 , and σ) 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 λ) were then used to generate forward simulations according to the full model without weekly   13)-( 16)), with one change point, and pre-specified values of λ 0 and λ 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 λ values) are considered, which each yield two sets of simulations: one corresponding to Forward Euler with ∆t = 1 day (as in [Dehning et al., 2020]) and another with ∆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 λ 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 λ(t), particularly during the (crucial) growth phase of the epidemic.

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 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 (λ, µ).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 λ value (cf. Figure 7), meaning a larger λ 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 ∆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 ∆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 rainfallrunoff models used in hydrology [Schoups andVrugt, 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: 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 σ.
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 rtol = 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 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  18)-( 22).For each parameter, the solid line indicates the likelihood calculated using an RK3(2) adaptive solver with rtol = atol = 10 −3 , while the dashed line indicates the likelihood calculated using the CVODE adaptive solver with rtol = atol = 10 −7 .estimated posterior medians.Likelihood surfaces are plotted for two adaptive step size solvers: the RK3(2) solver from SciPy with rtol = 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.

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

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.

Conflict of interest declaration
We declare that we have no competing interests.Using this in eq. ( 24),

Funding
As discussed further in the main text, this bound indicates an inverse relationship between σ and |L − L | 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 loglikelihood 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 ∼ N (g(x(t i ; θ)), σ).For brevity, as before, let g i = g(x(t i ; θ)) and ĝi = g(x i ), and define a i = g i − ĝi .Then, we can write y i = g i + i , where i ∼ N (0, σ).Starting from §S1, eq. ( 23), we have: Noting that the first term in eq. ( 26) is the sum of independent random variables each with mean 0 and variance (a 2 i /σ 2 ), we obtain the result: We . 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 σ is smaller.Also, in the case that g is the identity function,

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 [?].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: 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 change is the time when the stimulus changes in value.
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 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, ẋ(t = 0) = 0.Then, IID Gaussian noise was added to the solution at each of the sampled locations with σ = 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   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 .
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.

Figure 1 :
Figure1: 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 ∆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).

Figure 2 :
Figure2: 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 ∆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 ∆t.

Figure 3 :
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.

Figure 4 :
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 rtol = 10 −3 .

Figure 5 :
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 rtol = 10 −3 .

Figure 6 :
Figure6: Effect of noise on MCMC convergence.Data were generated according to the same specifications as for Figure4, with f 1 = −1, and the indicated values of σ.Inference was performed for the three parameters m, c, and k, as well as σ, 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 σ).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 rtol = 10 −3 .(Right) Forward simulation was performed using the RK5(4) solver with rtol = 10 −8 .

Figure 7 :
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 λ(t), and the bottom panel shows the number of daily cases resulting from these trajectories for each choice of the time step ∆t.

Figure 8 :
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 ∆t = 1 and ∆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 λ (eqs.(13)-(17)) and two different values for the ODE solver time step, ∆t.In both panels, lines indicate the posterior median and shaded regions indicate the central 95% of the posterior.

Figure 9 :
Figure9: Rainfall run-off model: inference using accurate and inaccurate adaptive solvers.Here, we plot the likelihood surface for each parameter for the rainfallrunoff model defined by eqs.(18)-(22).For each parameter, the solid line indicates the likelihood calculated using an RK3(2) adaptive solver with rtol = atol = 10 −3 , while the dashed line indicates the likelihood calculated using the CVODE adaptive solver with rtol = atol = 10 −7 .
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.
t i )| 2 , and so the average size of the underestimation decreases as the global truncation error decreases.

Figure S1 :
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 .