Troubleshooting NDSolveValue::rdelay Error in a Dengue Epidemic Model with Delay Differential Equations in Mathematica

Question

I am modeling the vector dynamics of a dengue epidemic using a system of delay differential equations (DDEs) in Wolfram Mathematica. My model includes a complex interplay of host-vector interactions, seasonality effects, and vaccination strategies.

However, I'm encountering a runtime error related to delayed time when attempting to solve the system using NDSolveValue. The specific error message is:

NDSolveValue::rdelay: Delayed time {-1. τ} = -1. τ computed at t = 0 did not evaluate to a real number.

The entire notebook can be found online at: https://www.wolframcloud.com/obj/lennertsaerens/Published/fergusonM.nb

Here is a detailed breakdown of the relevant parts of the model as I implemented it in my notebook:

Mosquito Reproduction Number and Lifecycle Parameters:

Rm = 2.69  (*Mosquito reproduction number*)
\[Omega] = 0.025 (*Larval mosquito mortality rate*)
\[CurlyEpsilon] = 0.1 (*Adult mosquito mortality rate*)
\[Alpha] = 1/19 (*Mean development time of mosquito larvae*)
(* We assign 𝑏 by fixing the required value of 𝑅𝑚,the mosquito \
reproduction number,which can be shown to be given by: *)
b = Rm / (\[Alpha]/(\[CurlyEpsilon] * (\[Alpha] + \[Omega]))) 
meq = 0.5 (*female mosquitos / person*)
K0 = meq / (\[Alpha]/\[CurlyEpsilon] (\[Alpha] (b - \
\[CurlyEpsilon])/(\[CurlyEpsilon] \[Omega]) - 
      1)) (*Mean larval mosquito carrying capacity*)
\[CapitalDelta]K = 0.3 (* Amplitude of seasonal variation in carrying \
capacity. Assigned value of 0.3 to match observed periodicity and \
amplitude of dengue epidemics. *)
POP = 1000000

K[t_] := K0*POP*(1 + \[CapitalDelta]K*Sin[2*Pi*t])

Plot[K[t], {t, 0, 10}]

Vaccination Parameters and Functions:

\[Beta]mh = 1  (* Per bite transmission probability from mosquitoes \
to humans. Should be set to match R0?? *)
\[Kappa]  = 0.6 (* Biting rate per mosquito / day *)
TD = 8 (*Mean duration of vaccine-induced protection against \
infection*)
VEMinus = 0.87 (*Maximum vaccine efficacy against infection in \
seronegative recipients*)
VEPlus = 0.4 (*Maximum vaccine efficacy against infection in \
seropositive recipients*)
A = 2 (*Age at vaccination*)
P[0] = 0.45 (*Primary infections which are symptomatic*)
P[1] = 0.8 (*Proportion of secondary infections which are symptomatic*)
P[2] = 0.12 (*Proportion of tertiary and quaternary infections which \
are symptomatic*)
P[3] = 0.12 (*Proportion of tertiary and quaternary infections which \
are symptomatic*)
P[4] = 0
P[5] = 0

S[v_, t_, a_, \[CapitalTheta]_] = 100000

(* The force of infection on humans due to serotype \
i,\[CapitalLambda]𝑖,is defined by: *)
\[CapitalLambda][i_, t_] := (\[Kappa] \[Beta]mh/M[t]) Y[t, i]
(*Print[\[CapitalLambda][i,t]]*)

(* The decay function h(𝜏) assumes exponential decay of protection \
after each vaccine dose with mean duration 𝑇: *)
h[\[Tau]_] := Piecewise[
  {{Exp[-\[Tau]/TD], \[Tau] < 0.5},
   {Exp[-(\[Tau] - 0.5)/TD], 0.5 <= \[Tau] < 1},
   {Exp[-(\[Tau] - 1)/TD], \[Tau] >= 1}}
  ]

Plot[h[\[Tau]], {\[Tau], 0, 10}, PlotLegends -> "Expressions"]

(* Heterologous vaccine-induced protection against infection decays \
over time and is described by the relative risk function 𝑓(𝜏),𝜏 being \
the time since vaccination,defined thus: *)
f[\[Tau]_, v_] := Piecewise[
  {{1, v == 0},
   {1 - VEMinus*h[\[Tau]], v == 1},
   {1 - VEPlus*h[\[Tau]], v == 2}}
  ]

Plot[{f[\[Tau], 0], f[\[Tau], 1], f[\[Tau], 2]}, {\[Tau], 0, 10}, 
 PlotLegends -> "Expressions"]

(*Incidence at time t of infection with serotype i in people with \
past exposure to serotype set theta, age a and vaccine status v*)
InfIncidence = 
 cI[v_, t_, a_, i_, \[CapitalTheta]_] := \[CapitalLambda][i, t]* 
   f[a - A, v]* S[v, t, a, \[CapitalTheta]]

(*Incidence of symptomatic disease at time t due to infection with \
serotype i in people with past exposure to serotype set theta age a \
and vaccine status v*)
SympIncidence = 
 cD[v_, t_, a_, i_, \[CapitalTheta]_] := 
  P[Length[\[CapitalTheta]] + (1 - 
       KroneckerDelta[0, v])]*\[CapitalLambda][i, t]*f[a - A, v]* 
   S[v, t, a, \[CapitalTheta]]

Mosquito Infectiousness:

\[Theta] = 2 (* Factor by which infectiousness of symptomatic \
infection exceeds that of asymptomatic infections. *)
\[Beta]hm = 1 (* Per bite transmission probability from humans to \
mosquitoes. Arbitrarily set to 1 *)
seroCombs = 
 Subsets[{1, 2, 3, 
   4}] (*Generates all possible combinations of serotypes*) 
vacStats = {0, 1, 2} (*All possible vaccination statuses*)
TIP = 6 (* Intrinsic incubation period *)
TInf = 4 (* infectious period in humans *)
maxAge = 120 (* Maximum age for humans *)

innerFunction[i_, t_, \[Tau]_, a_, \[CapitalTheta]_, v_] := 
 If[Not[MemberQ[\[CapitalTheta], i]], 
  cI[v, t - \[Tau], a, i, \[CapitalTheta]] + (\[Theta] - 1)*
    cD[v, t - \[Tau], a, i, \[CapitalTheta]], 0]

(*The force of infection on mosquitoes due to serotype \
i,\[CapitalPsi]𝑖,is defined by:*)
\[CapitalPsi][i_, t_] := \[Kappa]*\[Beta]hm/POP*
  Integrate[
   Sum[
    innerFunction[i, t, \[Tau], a, \[CapitalTheta], v],
    {\[CapitalTheta], seroCombs},
    {v, vacStats}
    ],
   {a, 0, maxAge},
   {\[Tau], TIP, TIP + TInf}
   ]

Differential Equations:

\[Omega] = 0.025 (*Larval mosquito mortality rate*)
\[Eta] = 1/10 (*Mean extrinsic incubation period*)

eqLarvae = 
 D[L[t], t] == 
  b*M[t] - \[Alpha]*L[t] - \[Omega]*L[t] (1 + L[t]/K[t]) (*Larvae*)
eqAdults = 
 D[Am[t], t] == \[Alpha]*L[t] - 
   Sum[\[CapitalPsi][i, t]*Am[t], {i, 4}] - \[CurlyEpsilon]*
    Am[t] (*Uninfected adult mosquitos*)
eqInfected = 
 Flatten[Table[
   D[H[t, i, j], 
     t] == (KroneckerDelta[1, j]* \[CapitalPsi][i, t]*
       Am[t]) + (4 \[Eta]*(1 - KroneckerDelta[1, j])*
       H[t, i, j - 1]) - (4 \[Eta] + \[CurlyEpsilon])*H[t, i, j], {i, 
    4}, {j, 4}]] (* Infected but not yet infectious mosquitos *)
eqInfectious = 
 Flatten[Table[
   D[Y[t, i], t] == 4 \[Eta]*H[t, i, 4] - \[CurlyEpsilon]*Y[t, i], {i,
     4}]] (* Infectious mosquitos *)

M[t_] := 
 Am[t] + Sum[H[t, i, j], {i, 4}, {j, 4}] + 
  Sum[Y[t, i], {i, 4}] (*Total adult mosquitos, 
  everything except larvae*)

histories = {
  L[t /; t <= 0] == 0,
  Am[t /; t <= 0] == 100000,
  H[t /; t <= 0, i, j] == 0 /. {i -> #, j -> #2} & @@@ 
   Tuples[Range[4], 2],
  Y[t /; t <= 0, i] == 1000 /. i -> # & /@ Range[4]
  }

Solving the Differential Equations:

solution = NDSolveValue[
  Flatten[{eqLarvae, eqAdults, eqInfected, eqInfectious, histories}],
  Flatten[{L[t], Am[t], H[t, 1, 1], Y[t, 1]}],
  {t, 0, 200},
  Method -> {"EquationSimplification" -> "Residual"}
  ]

Upon execution, I receive the error message related to the delayed time. I suspect the issue is with the handling of the delay term in the DDEs, but I'm not certain how to proceed with the corrections. The delay terms are integral to the model, representing the intrinsic incubation period and infectious period in humans which are necessary for calculating the force of infection on mosquitoes.

Could someone with experience in Mathematica's handling of DDEs provide insight into what might be causing this error and how to resolve it? Suggestions for rewriting the delay terms or adjusting the NDSolveValue implementation would be greatly appreciated.

In my attempt to model the dynamics of a dengue epidemic using delay differential equations (DDEs) in Wolfram Mathematica, I encountered a runtime error related to the handling of delayed terms when using NDSolveValue. To address this, I meticulously added initial history conditions for each function where a delay term appears, as required by Mathematica's DDE solver.

I expected that by properly initializing these histories, Mathematica would be able to handle the delayed terms and solve the system of equations. However, despite this initialization, I still received an error message indicating that the delayed time did not evaluate to a real number at t = 0. This was surprising, as I believed that the initial histories would provide the necessary information for all time points required by the delay terms.

To further investigate the issue, I modified the system by removing the delay terms altogether. After this alteration, the system worked as intended, suggesting that the error is specifically associated with the implementation of the delay. The successful results from this non-delayed system are documented at the bottom of my online Mathematica notebook, which indicates that the rest of the model is correctly specified and that the issue is isolated to the handling of delays.

Answer 1

At the time of posting this question, the NDSolve implementation in the latest version of Mathematica (version 14) did not support delay differential equations with an infinite number of delays. Because of this limitation, changing the integral that represents the incubation period in humans to a sum fixed the issue. Many thanks to Chris for suggesting this change in the comments of the question.

New version of the force of infection on mosquitos:

Ψ[i_, t_] := κ*βhm/POP*
  Integrate[
   Sum[
    innerFunction[i, t, τ, a, Θ, v],
    {τ, TIP, TIP + TInf},
    {Θ, seroCombs},
    {v, vacStats}
    ],
   {a, 0, maxAge}
   ]

With the rest of the code unchanged, I am now able to solve the system of equations describing the vector dynamics of a dengue epidemic.