Infection severity across scales in multi-strain immuno-epidemiological Dengue model structured by host antibody level

Abstract

Infection by distinct Dengue virus serotypes and host immunity are intricately linked. In particular, certain levels of cross-reactive antibodies in the host may actually enhance infection severity leading to Dengue hemorrhagic fever (DHF). The coupled immunological and epidemiological dynamics of Dengue calls for a multi-scale modeling approach. In this work, we formulate a within-host model which mechanistically recapitulates characteristics of antibody dependent enhancement in Dengue infection. The within-host scale is then linked to epidemiological spread by a vector–host partial differential equation model structured by host antibody level. The coupling allows for dynamic population-wide antibody levels to be tracked through primary and secondary infections by distinct Dengue strains, along with waning of cross-protective immunity after primary infection. Analysis of both the within-host and between-host systems are conducted. Stability results in the epidemic model are formulated via basic and invasion reproduction numbers as a function of immunological variables. Additionally, we develop numerical methods in order to simulate the multi-scale model and assess the influence of parameters on disease spread and DHF prevalence in the population.

Appendix: Numerical convergence rates

In this section, we provide tables showing computed rates and order of convergence for numerical experiments of the finite difference and multi-scale simulation procedure described in Sect. 4. For the numerical tests, we calculate the error in norm between computed solutions of the t, y stepping method at certain step sizes $\varDelta t,h=\varDelta y$ and reference solutions at some final time $t=T$. We utilize three different types of reference solutions: (i) the numerically approximated equilibrium given by our derived formula (31), (ii) solution of the numerical scheme with smallest step sizes $\varDelta t,{\tilde{h}}=\varDelta y$, and (iii) solution of the numerical scheme with step sizes multiplied by factor of 1/2, $\frac{\varDelta t}{2},\frac{h}{2}$. For each error calculation at step size h, $e_h$, we form a sequence by successively decreasing step size by 1/2, whereby we compute order of convergence by $\log _2(e_h/e_{h/2})$. Furthermore, we consider two different scenarios: (a) we start the initial condition where infected vectors, $I_v^1,I_v^2$ are slightly perturbed from ${\mathcal {E}}_0$ (outbreak scenario) with final time $T=50$ days, (b) we start the initial condition at numerically calculated equilibrium with final time $T=500$ days. For the former scenario (a), we do not use numerically calculated equilibrium as a reference solution since this may be far off from simulation at $t=50$.

We compute the different orders of convergence because there are several sources of error and to test different initial condition scenarios. Our method relies on distinct algorithms in addition to the finite difference scheme, such as Runge-Kutta method for within-host ODE (ode45 in MatLab), interpolation, integration and, in the case of numerical equilibrium formula, nonlinear root-finding. Each routine can produce error, which can also propagate in the form of discontinuities in recovered distribution corresponding to an influx of recovery from primary infected individuals with pre-existent antibody levels at a certain mesh points from the initial susceptible antibody distribution. In order to efficiently reduce error we utilize the trapezoidal integration when integrating with respect to initial susceptible antibody level $y_0$, but left endpoint integration for other antibody variables since there is small number of mesh points ($M_0$) for $s(\cdot , y_0)$ when compared to the range of antibody levels after infection and waning. We do also provide one numerical test with only left-endpoint integration shown in last two tables, which gives more error than trapezoidal, but has more regular order of convergence pattern. In addition, we include comparisons with a larger step size ($\varDelta y=0.06$) for scenarios (a) and the last two tables, which forces a point distribution for susceptible antibody levels, creating different error structure. Overall, from the different numerical tests, we observe convergence to certain error rates within a particular compartment and/or test scenario, ranging from orders that are sub-linear ($<1$) to larger than quadratic ($>2$). When comparing with reference solutions computed by numerical simulation at smaller step size, the order of convergence is mostly faster than linear (Tables 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16).

Table 5 Error analysis of $i^h_1(50,y,y_0)$ with initial $I_v^1(0)=I_v^2(0)=0.02$, and other components starting at ${\mathcal {E}}_0$, for step size $\varDelta y=h$ compared to reference solution $i^{{\tilde{h}}}_1$ (${\tilde{h}}=\frac{0.00375}{2}$) and $i^{h/2}_1$, respectively, in $L^1$ norm

Table 6 Error analysis of $i^h_{12}(50,y,y_0)$ with initial $I_v^1(0)=I_v^2(0)=0.02$, and other components starting at ${\mathcal {E}}_0$, for step size $\varDelta y=h$ compared to reference solution $i_{12}^{{\tilde{h}}}$ (${\tilde{h}}=\frac{0.00375}{2}$) and $i^{h/2}_{12}$, respectively, in $L^1$ norm

Table 7 Error analysis of $r^h_1(50,y)$ with initial $I_v^1(0)=I_v^2(0)=0.02$, and other components starting at ${\mathcal {E}}_0$, for step size $\varDelta y=h$ compared to reference solution $r^{{\tilde{h}}}_1$ (${\tilde{h}}=\frac{0.00375}{2}$) and $r^{h/2}_1$, respectively, in $L^1$ norm

Table 8 Error analysis of $I_v^{1,h}(50)$ with initial $I_v^1(0)=I_v^2(0)=0.02$, and other components starting at ${\mathcal {E}}_0$, for step size $\varDelta y=h$ compared to reference solution $I_v^{{\tilde{h}}}$ (${\tilde{h}}=\frac{0.00375}{2}$) and $I_v^{h/2}$, respectively, in $L^1$ norm

Table 9 Error analysis of $s^h(50,y)$ with initial $I_v^1(0)=I_v^2(0)=0.02$, and other components starting at ${\mathcal {E}}_0$, for step size $\varDelta y=h$ compared to reference solution $s^{{\tilde{h}}}$ (${\tilde{h}}=\frac{0.00375}{2}$) and $s^{h/2}$, respectively, in $L^1$ norm

Table 10 Error analysis of $i^h_1(500,y,y_0)$ with initial condition set at numerically calculated equilibrium ${\bar{i}}_1(y,y_0)$ for step sizes $\varDelta y=h$ compared to ${\bar{i}}_1$, reference solution $i^{{\tilde{h}}}_1$ (${\tilde{h}}=.00375$) and $i^{h/2}_1$, respectively, in $L^1$ norm

Table 11 Error analysis of $i^h_{12}(500,y,y_0)$ with initial condition set at numerically calculated equilibrium ${\bar{i}}_{12}(y,y_0)$ for step sizes $\varDelta y=h$ compared to ${\bar{i}}_{12}$, reference solution $i^{{\tilde{h}}}_{12}$ (${\tilde{h}}=.00375$) and $i^{h/2}_{12}$, respectively, in $L^1$ norm

Table 12 Error analysis of $r^h_1(500,y)$ with initial condition set at numerically calculated equilibrium ${\bar{r}}_1(y)$ for step sizes $\varDelta y=h$ compared to ${\bar{r}}_1$, reference solution $r^{{\tilde{h}}}_1$ (${\tilde{h}}=.00375$) and $r^{h/2}_1$, respectively, in $L^1$ norm

Table 13 Error analysis of $I_v^{1,h}(500)$ with initial condition set at numerically calculated equilibrium ${\bar{I}}_v$ for step sizes $\varDelta y=h$ compared to ${\bar{I}}_v$, reference solution $I_v^{{\tilde{h}}}$ (${\tilde{h}}=.00375$) and $I^{h/2}_v$, respectively

Table 14 Error analysis of $s^h(500,y)$ with initial condition set at numerically calculated equilibrium ${\bar{s}}(y)$ for step sizes $\varDelta y=h$ compared to ${\bar{s}}$, reference solution $s^{{\tilde{h}}}$ (${\tilde{h}}=.00375$) and $s^{h/2}$, respectively, in $L^1$ norm

Table 15 Error analysis with initial condition set at numerically calculated equilibrium ${\mathcal {E}}_c$ for step sizes $\varDelta y=h$ compared to ${\mathcal {E}}_c$, using left-end point approximation integration.

Table 16 Error analysis with initial condition set at numerically calculated equilibrium ${\mathcal {E}}_c$ for step sizes $\varDelta y=h$ compared to ${\mathcal {E}}_c$, using left-end point approximation integration