Applying stiff integrators for ODEs and DDEs to problems with distributed delays

TL;DR

The paper tackles numerical solution of models with distributed delays by approximating the delay kernel with a sum of exponentials and converting the integro-differential problem into an augmented system of ODEs/DDEs. This enables applying established stiff integrators like Radau5 and its delay-enabled variant Radar5 without redesigning time-stepping schemes. A key contribution is an efficient linear-algebra strategy that reduces the augmented-systems solve to a rank-one update of the original Jacobian, yielding near-cubic cost in the original dimension plus linear cost in the augmentation size. The authors provide practical kernel-approximation techniques (Beylkin–Monzón and Braess–Hackbusch), detail parameter selection for gamma and Pareto kernels, and demonstrate substantial accuracy and speedups on pharmacokinetic/pharmacodynamic tests and a chemotherapy-model application. The approach broadens the applicability of powerful stiff solvers to Volterra-type and distributed-delay problems, with clear implications for pharmacology and related fields.

Abstract

There exist excellent codes for an efficient numerical treatment of stiff and differential-algebraic problems. Let us mention {\sc Radau5} which is based on the $3$-stage Radau IIA collocation method, and its extension to problems with discrete delays {\sc Radar5}. The aim of the present work is to present a technique that permits a direct application of these codes to problems having a right-hand side with an additional distributed delay term (which is a special case of an integro-differential equation). Models with distributed delays are of increasing importance in pharmacodynamics and pharmacokinetics for the study of the interaction between drugs and the body. The main idea is to approximate the distribution kernel of the integral term by a sum of exponential functions or by a quasi-polynomial expansion, and then to transform the distributed (integral) delay term into a set of ordinary differential equations. This set is typically stiff and, for some distribution kernels (e.g., Pareto distribution), it contains discrete delay terms with constant delay. The original equations augmented by this set of ordinary differential equations can have a very large dimension, and a careful treatment of the solution of the arising linear systems is necessary. The use of the codes {\sc Radau5} and {\sc Radar5} is illustrated at three examples (two test equations and one problem taken from pharmacodynamics). The driver programs for these examples are publicly available from the homepages of the authors.

Figures (3)

• Figure 4.1: Relative error $| t^{-\alpha } - T(t,h) |/ t^{-\alpha}$ as a function of the step size $h$. The different curves correspond to $\alpha = 0.1, 0.2, \ldots , 0.9$ (bold curve for $\alpha = 0.5$). \newlabelfig:err-trap
• Figure 4.2: Relative error of the the truncation errors as a function of $M$ and $N$, respectively. \newlabelfig:err-trunc
• Figure 6.1: Solution components $y(t)$ and $w(t)$ of the system \ref{['example3eh']} for the parameter values of Table \ref{['tab:31eh']}: upper row (problem 1) and lower row (problem 2). \newlabelfig:solution