使用延迟微分方程修复登革热流行模型的 Wolfram Mathematica 中的延迟时间错误“NDSolveValue::rdelay”

问题描述 投票：0回答：1

我正在使用 Wolfram Mathematica 中的延迟微分方程 (DDE) 系统对登革热流行病的矢量动力学进行建模。我的模型包括宿主-载体相互作用、季节性效应和疫苗接种策略之间复杂的相互作用。

但是，在尝试使用 

NDSolveValue
求解系统时，我遇到了与延迟时间相关的运行时错误。具体错误信息为：

NDSolveValue::rdelay：延迟时间{-1。 τ} = -1。在 t = 0 时计算的 τ 未计算为实数。

整个笔记本可以在线找到：https://www.wolframcloud.com/obj/lennertsaerens/Published/fergusonM.nb

以下是我在笔记本中实现的模型相关部分的详细分解：

蚊子繁殖数量和生命周期参数：

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}]

疫苗接种参数及功能：

\[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]]

蚊子传染性：

\[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}
   ]

微分方程：

\[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]
  }

求解微分方程：

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"}
  ]

执行后，我收到与延迟时间相关的错误消息。我怀疑问题出在 DDE 中延迟项的处理上，但我不确定如何进行更正。延迟项是模型的组成部分，代表人类的固有潜伏期和感染期，这是计算蚊子感染力所必需的。

有 Mathematica 处理 DDE 经验的人能否深入了解可能导致此错误的原因以及如何解决它？如果您提出重写延迟条款或调整 

NDSolveValue
实施的建议，我们将不胜感激。

在尝试使用 Wolfram Mathematica 中的延迟微分方程 (DDE) 对登革热流行的动态进行建模时，我在使用 

NDSolveValue
时遇到了与延迟项处理相关的运行时错误。为了解决这个问题，我按照 Mathematica 的 DDE 求解器的要求，为每个出现延迟项的函数精心添加了初始历史条件。

我期望通过正确初始化这些历史记录，Mathematica 将能够处理延迟项并求解方程组。然而，尽管进行了这样的初始化，我仍然收到一条错误消息，表明延迟时间并未计算为 

t = 0
处的实数。这是令人惊讶的，因为我相信初始历史将为延迟项所需的所有时间点提供必要的信息。

为了进一步调查该问题，我修改了系统，完全删除了延迟项。进行此更改后，系统按预期工作，这表明该错误与延迟的实施具体相关。这个非延迟系统的成功结果记录在我的在线 Mathematica 笔记本的底部，这表明模型的其余部分已正确指定，并且该问题与延迟的处理无关。

debugging simulation wolfram-mathematica differential-equations
1个回答
0
投票

在发布此问题时，最新版本 Mathematica（版本 14）中的 

NDSolve
实现不支持具有无限数量延迟的延迟微分方程。由于这一限制，将代表人类潜伏期的积分更改为总和可以解决该问题。非常感谢 Chris 在 Mathematica Stack 交换的问题评论中提出了这一更改。

新版蚊子感染力：

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

在其余代码不变的情况下，我现在能够求解描述登革热流行病矢量动力学的方程组。

最新问题
© www.soinside.com 2019 - 2024. All rights reserved.