Numerical methods and hypoexponential approximations for gamma distributed delay differential equations

Abstract

Gamma distributed delay differential equations (DDEs) arise naturally in many modelling applications. However, appropriate numerical methods for generic gamma distributed DDEs have not previously been implemented. Modellers have therefore resorted to approximating the gamma distribution with an Erlang distribution and using the linear chain technique to derive an equivalent system of ordinary differential equations (ODEs). In this work, we address the lack of appropriate numerical tools for gamma distributed DDEs in two ways. First, we develop a functional continuous Runge-Kutta (FCRK) method to numerically integrate the gamma distributed DDE without resorting to Erlang approximation. We prove the fourth-order convergence of the FCRK method and perform numerical tests to demonstrate the accuracy of the new numerical method. Nevertheless, FCRK methods for infinite delay DDEs are not widely available in existing scientific software packages. As an alternative approach to solving gamma distributed DDEs, we also derive a hypoexponential approximation of the gamma distributed DDE. This hypoexponential approach is a more accurate approximation of the true gamma distributed DDE than the common Erlang approximation but, like the Erlang approximation, can be formulated as a system of ODEs and solved numerically using standard ODE software. Using our FCRK method to provide reference solutions, we show that the common Erlang approximation may produce solutions that are qualitatively different from the underlying gamma distributed DDE. However, the proposed hypoexponential approximations do not have this limitation. Finally, we apply our hypoexponential approximations to perform statistical inference on synthetic epidemiological data to illustrate the utility of the hypoexponential approximation.

Keywords: delay differential equations; functional continuous Runge–Kutta methods; infinite delay equation; linear chain trick.

Fig. 1.

Trajectory dependence on the shape parameter is not smooth. (A) The trajectory and for calculated with the fixed (black) and smoothed (red) hypoexponential approximations in (3.6). The exact solution is not shown as it overlaps with the two curves. (B) The graph of with , using the fixed (black) and smoothed (red) parameterizations of the hypoexponential approximation, and the exact solution (blue).

Fig. 2.

Convergence plots for the linear test problem (4.1). We plot as a function of . The slope of gives the convergence rate. is the simulation of (4.1) using the fourth-order FCRK method from Section 2 with fixed step size and is the solution of the equivalent ODE (4.2). Figure (A) shows the comparison against the exact solution when , while figures (B) and (C) show the error between and for and , respectively. Figure (D) shows the solution of the DDE for each test case. The solution of the equivalent ODE (4.2) is calculated using the fourth-order RK method RK45 in Matlab with relative and absolute error tolerance of .

Fig. 3.

Convergence plots for the nonlinear test problem (4.3). We plot as a function of . The slope of gives the convergence rate. is the simulation of (4.3) using the fourth-order FCRK method from Section 2 with fixed step size and is the solution of the equivalent ODE (4.4). Panels (A—C) show the error between and for , respectively. Panel (D) shows the solution of the DDE for each test case. The solution of the equivalent ODE (4.4) is calculated using the fourth-order RK method RK45 in Matlab with relative and absolute error tolerance of .

Fig. 4.

Convergence plots for the linear gamma-distributed DDE test problem (4.5). We plot as a function of where is the simulation of (4.5) using the fourth-order FCRK method from Section 2 and is the analytical solution of the linear-distributed DDE. In B, the error reaches machine precision for A–C show the error between and for the parameter triples given by and , respectively.

Fig. 5.

Comparison of ODE approximations to the gamma-distributed DDE (1.1) using the Erlang approximation in equation (3.2) or the fixed hypoexponential approximation in (3.5). In all cases, the solution of the gamma-distributed DDE as solved using the FCRK method is in solid blue, the solution of the fixed hypoexponential two moment approximation is in dashed orange and the solution of the Erlang approximation is in purple. A–C show the solution of the linear test problem (4.1) for and , respectively. D–F show the solution of the nonlinear test problem (4.3) for and , respectively.

Fig. 6.

Comparison of ODE approximations to the gamma-distributed DDE (4.5) using the Erlang approximation in equation (3.2) or the fixed hypoexponential two moment approximation in (3.5) showing that the Erlang approximation does not have the same stability properties as the gamma-distributed DDE or the hypoexponential approximation. In all cases, the solution of the gamma-distributed DDE as solved using the FCRK method is in solid blue, the solution of the hypoexponential approximation is in dashed orange and the solution of the Erlang approximation is in purple.

Fig. 7.

Simulated epidemic data and posterior predictive checks. (A) Simulated data and the model prediction . The dark-blue band represents the credible interval, and the light-blue band the prediction interval. (B) Simulated serial intervals represented as a empirical survival function (black), and the fitted survival function for given by (blue).

Fig. 8.

Posterior marginal and joint density for the epidemic model parameters. Each dot in the joint density scatter plots represents a Monte-Carlo sample from the posterior distribution. The color of the dots indicates the density. The black vertical lines represent the ground-truth parameter values. The ground-truth parameters used for the simulation are , , , , and .