Calculation of final size for vector-transmitted epidemic model

1.   Introduction

Threatening vector-transmitted epidemics such as Malaria, Zika virus infection and Dengue infection are mediated by mosquitoes. Demand for practical intervention to prevent vector-transmitted infectious disease has been increasing in non-tropical countries where an increase of incidences has been reported possibly due to elevated global transportation and global warming. Quantitative indicator for epidemics is indispensable for estimating the impact of epidemics.

Mathematical models of vector-transmitted infectious disease have been applied not only to understand qualitative behavior but also to define quantitative indicators of an epidemic process ^{[1,2,5,6,5,6,7,8]}. The number of sub-population experiencing infection during an epidemic process, referred to as final size, is a useful indicator to estimate the impact of epidemics. The final size of a susceptible population can be numerically computed for the classical Kermack McKendrick equation ^[9,10,11,12] (see Appendix 4 for a brief description and derivation). In a practical setting, final size calculation has been used to estimate possible damage quantitatively and to prepare possible prevention.

Despite the usefulness of final size calculation, application of final size calculation is limited to a specific class of epidemic models such as the classical Kermack McKendrick model. Although some extensions to increase the usability of final size calculation have been proposed, the majority of practical applications focus only on the classical Kermack McKendrick model. To our best knowledge, however, final size equations for vector transmission epidemic models have not been reported.

In this paper, we consider a vector-transmitted epidemic model and derive a final size equation. For comparison, we derive final size equations for several epidemic models. Our focus includes epidemic models with mass-action type/standard incidence transmission rate or vector transmission. The organization of the present paper is as follows. In the next section, we propose a vector-transmitted epidemic model with standard incidence rate. For the main model, quasi-steady state approximation is applied to obtain a simpler model, enabling further mathematical analysis and derivation of a final size equation. One of significant findings is the presence of a threshold curve which determines the existence of two bistable distinct final sizes when infection induced death exists.

2.   Model formulation

We consider the following system of differential equations:

where variables in the human compartments $S$, $I$ and $R$ represent the number of susceptible, infective and recovered individuals, respectively. Note that the total population $N = S+I+R$ satisfies $\frac{dN}{dt} = 0$: the population is closed. Parameters in the human compartment $\beta_1$ and $\gamma$ denote transmission coefficient from infective mosquito to susceptible human, and recovery rate. Similarly, variables $V$ and $W$ represent the number of susceptible and infective mosquitoes. Parameters $g$, $\beta_2$ and $\mu$ denote constant reproduction rate, transmission coefficient from infective human to susceptible mosquito, and death rate of mosquitoes, respectively. In model (2.1), infection from mosquitoes to humans is given by $\beta_1 \frac{S}{N} W$. Similarly, infection from humans to mosquitoes is given by $\beta_2 \frac{I}{N} V$. Let $d$ denote the average death rate of humans and $\varepsilon \mathrel{\mathop:} = \frac{d}{\mu}$. Hereafter we assume that the death rate of mosquitoes is higher than that of humans, that is, $\varepsilon = \frac{d}{\mu} \ll 1$. This assumption leads to apply a quasi-steady state approximation (QSSA) to model (2.1). QSSA is a convenient approximation which is widely used to reduce the dimension of a dynamical system when a variable representing fast dynamics is included. By QSSA, variables representing fast dynamics can be ignored from the main system (a concrete example can be found in ^[13]). In our case, the fourth and fifth equations of (2.1) are rewritten as

By solving the fourth equation with respect to $V(t)$,

Futhermore, by substituting this into the fifth equation,

Hence we obtain that

By substituting (2.3) into the first and second equations of (2.1), we obtain that

where the force of infection is $\beta_1 \frac{\beta_2 gI}{{\mu}^2N+\beta_2 \mu I}$. Model (2.4) can be interpreted as a vector transmission epidemic model with a nonlinear incidence rate. In the subsequent section, we consider two cases whether infectives can recover or not (i.e., die). Model (2.4) corresponds to a vector-transmitted epidemic model with recovery. Hereafter, we refer to model (2.4) as recovery model. Dynamics of model (2.4) with $\beta_1 = \beta_2 = 0.01$, $g = 1000$, $\mu = \frac{1}{30}$ and $\gamma = 0.2$ is summarized in the caption of Figure 1.

3.   Comparison of final size equations

3.1.   Recovery model

In this subsection, we derive the final size equation for recovery model (2.4). It follows from the first and second equations of (2.4) that

By solving the first equation of (2.4) with respect to $I$, we obtain that

Then (3.1) is rewritten as

where

Hence $\left|-\frac{\mu N}{\beta_1 g} \frac{d}{dt} \ln{S}\right| < 1$. By applying a power series expansion, we obtain that

Hence

By integrating both sides of (3.4) from 0 to $\infty$ with respect to $t$, we obtain the following final size equation:

where $S(0)+I(0) = N$. Note that $I(\infty)$ must satisfy $I(\infty) = 0$. In fact, $R(t) \to \infty$ as $t \to \infty$ when $I(\infty) > 0$, which contradicts to $S(t)+I(t)+R(t) = N$ (constant).

We show that the final size equation for the classical Kermack McKendrick equation can be obtained as a first order approximation of (3.5) for sufficiently large $N$. In fact, if $N$ is large, then

Then final size equation (3.5) is reduced to

Note that this corresponds to the final size equation for the classical Kermack McKendrick model (see Appendix 4 for details). First order approximation can be a good estimate for a particular set of parameters. Even for small total population size such as $N = 200$, Figure 2 shows a good correspondence between numerical solutions of (3.6) obtained by the Newton method and numerical simulation results of recovery model (2.4) for sufficiently large time.

3.2.   Non-recovery model

In this subsection, we consider a situation that no recovered individuals exist. In other words, we assume that $\frac{dR}{dt} = 0$ and $R(0) = 0$ in recovery model (2.4). Then we obtain the following system of differential equations:

For model (3.7), parameter $\gamma$ can be interpreted as infection induced death rate. Hereafter model (3.7) is referred to as non-recovery model. Since $\frac{dS}{dt}+\frac{dI}{dt}+\frac{dR}{dt} = -\gamma I$ and $R(t) = 0 \ (t \geq 0)$, the total population $S(t)+I(t)+R(t) = S(t)+I(t) = N$ depends on time and is monotonically decreasing with respect to time.

The first equation of (3.7) is reduced to the following quadratic equation with respect to $I$:

Define $c_0$, $c_1$ and $c_2$ by

Since $dS/dt < 0$, $c_0 < 0$ and $c_2 < 0$. Furthermore we can show that $c_1 > 0$. In fact, for any $S, I > 0$,

This implies that there exists a possibility for (3.8) to have two positive roots if and only if $c_1^2-4c_0c_2 > 0$. Let us denote a positive root of (3.8) by $I^+(S(t), dS(t)/dt)$. By integrating both sides of the first equation of (3.7) from $0$ to $\infty$ with respect to $t$, we obtain the final size equation for model (3.7):

Two panels of Figure 3 show numerical simulation results for non-recovery model (3.7) with $\beta_1 = \beta_2 = 0.01$, $g = 1000$, $\mu = \frac{1}{30}$, $\gamma = 0.2$. Interestingly, the final size differs between two panels. More specifically, $S(\infty) = 0$ on the left panel while $S(\infty) > 0$ on the right panel. By implementing several numerical simulations, we also find that there exists a threshold curve which separates a whole region into two sub-regions. In one region, all trajectories of model (3.7) converge to a trivial equilibrium point at which $S(\infty) = 0$, while each trajectory in another sub-region converges to a different point on the $S$-axis at which infectives disappears (see two panels of Figure 4). Although analytical characterization for the existence of the threshold curve has not been obtained, hereafter we perform extensive numerical simulations to draw and to figure out the property of a threshold curve.

Let us define $\bar{S}$ by

Then we show that the basic reproduction number defined for an equilibrium of non-recovery model (3.7) corresponds to $1$ if $S(0) = \bar{S}$. In fact, it follows from the second equation of model (3.7) that for $I(0) \approx 0$,

Infectives can increase at the initial phase if

Therefore the basic reproduction number is given by

The threshold curve, denoted by $I = f(S)$, is defined as a curve satisfying $0 = f(\bar{S})$. In other words, the intersection of the threshold curve and $S$-axis is point $(\bar{S}, 0)$. For convenience, we define $\tilde{S}$ with which $I$-nullcline takes its maximum value:

Let $L_1$ (or $L_2$) denote a boundary which defines the interface of $A \cup B$ and $C \cup D$ (or $D$ and $E$), respectively (see $L_1$ and $L_2$ in Figure 4):

Sub-regions $A$-$E$ are defined by $L_1$, $L_2$ or the inner and outer of $I$-nullcline as follows:

To prove that solutions converge to $S(\infty) > 0$ or 0 depending on the initial values, we classify $\mathbb{R}_+^2$ into five regions (see the right panel of Figure 4). We introduce the following propositions:

Proposition 1. $\frac{dS}{dt} < 0$ for any $t < \infty$.

Proposition 2. Any trajectories cross $I$-nullcline horizontally along $S$-axis if and only if initial condition is taken in region $B$, $C$ or $D$.

Proposition 3. Threshold curve $I = f(S)$ entirely separates region $E$ and the others.

By Propositions 1-3, we obtain the following observations. If $(S(0), I(0)) \in {A}$, then any trajectories remain to stay in ${A}$, and tend to trivial equilibrium $(0, 0)$. If $(S(0), I(0)) \in {B}$, then any trajectories cross the left part of $I$-nullcline within a finite time, and finally enter ${A}$. Similarly, if $(S(0), I(0)) \in {C}$, then any trajectories cross $L_1$ within a finite time, and finally enter ${B}$. If $(S(0), I(0)) \in {D}$, then any trajectories enter either ${C}$ or ${A}$ without crossing $L_2$. Finally, if $(S(0), I(0)) \in {E}$, then any trajectories remain in ${E}$, and tend to a non-trivial equilibrium point $\left(\frac{\beta_1 \beta_2 g}{\gamma {\mu}^2}+c, 0 \right)$, where $c > 0$ is a constant and depends on $(S(0), I(0))$.

3.2.1.   Explicit final size for less fatal infection model

Although the final size equation for non-recovery model (3.7) is obtained as (3.10), it is difficult to calculate a final size. Here we impose a feasible assumption for non-recovery model (3.7), and derive an explicit final size for the simpler model. More precisely, we assume that infection induced death rate $\gamma$ is small: $\gamma \ll 1$. Then for sufficiently large $t$, $S(t)+I(t) \simeq N$. This leads to the following approximation:

Whenever $I(t)$ is small enough, we apply a linear approximation to (3.15).

Define $\beta$ by $\beta: = \frac{\beta_1\beta_2g}{\mu^2N}$. Then non-recovery model (3.7) is reduced to

Let us refer to model (3.17) as less fatal infection model. Note that model (3.17) has been studied in ^[6] in a different context. Two panels in Figure 5 show numerical simulation results. The basic reproduction number defined for less fatal infection model (3.17) is $\mathcal{R}_0 = \frac{\beta}{\gamma}$, which corresponds to the classical Kermack McKendrick model (see Appendix 4 for details). Furthermore, $I-$nullcline is a straight line which has slope $\frac{\beta}{\gamma}-1$ through the origin. In other words, $\mathcal{R}_{0} > 1$ if and only if the slope is positive (see Figure 6). From the first and second equations of (3.17),

By dividing both sides of (3.18) by $S+I$, we obtain that

By integrating both sides of (3.19) from 0 to $\infty$ with respect to $t$, we have

Note that an equilibrium state of (3.17) is given by the solution of the following system of equations:

Hereafter we show that $I(\infty) = 0$ by contradiction. Suppose that $I(\infty) > 0$. By dividing both sides of the second equation of (3.20) by $I(\infty)$,

It follows from the first equation of (3.20) that $\gamma = 0$, which contradicts to the assumption $\gamma > 0$. Thus, $I(\infty) = 0$. Then we obtain that

for $S(\infty) > 0$. In summary, the explicit form of final size $S(\infty)$ is given as follows:

We also show that (3.21) can hold only if $\beta < \gamma$. In fact, if $\beta \geq \gamma$, then

which is a contradiction. Since $\beta < \gamma$,

Hence we have

Moreover, by taking contraposition,

In conclusion, the explicit form of the final size $S(\infty)$ is given by

This indicates that both susceptible and infective individuals vanish if $\mathcal{R}_{0} \geq 1$. By contrast, if $\mathcal{R}_{0} < 1$, only a fraction of individuals vanishes (see Figure 6). Note that this result is consistent with numerical computation results in subsection 3.2.

4.   Conclusion

In this paper, we formulate a vector-transmitted epidemic model which describes interactions among susceptible/infected human individuals and susceptible/infected mosquitoes. By assuming a fast turnover rate of mosquito life-cycle, quasi-steady state approximation was applied to model (2.1) to obtain a simpler model. Explicit forms of final size equations were obtained for different epidemic models. We derived important indicators which characterize epidemics for each epidemic models as summarized in Table 1. Numerical computation of final size for non-recovery model (3.7) exhibits a qualitatively distinct threshold phenomenon which is not observed for the classical Kermack McKendrick model: there exists a threshold curve which separates the whole region into two sub-regions. Interestingly, final sizes differ between these two sub-regions. In other words, the initial number of infectious individuals may be crucial for the epidemic outcome of vector-transmitted infection if infection induced death occurs.

Acknowledgements

This research was partly supported by JST PRESTO Grant Number JPMJPR16E9, the Japan Society for the Promotion of Science (JSPS) through the "Grant-in-Aid (C) 16K05265 (to S.N.)" and "Grant-in-Aid 26400211 (to Y.T.)". The authors are grateful to anonymous referees whose valuable suggestions helped improve the quality of this paper.

Conflict of interest

The authors declare no competing interests.

Appendix

Size equation for the classical Kermack McKendrick model

The classical Kermack and McKendrick model is given by

Definition for variables and parameters of model (4.1) are written in Section 2. Note that the total population $N = S+I+R$ satisfies $\frac{dN}{dt} = 0$. The basic reproduction number for model (4.1) is given by

Note that

Moreover, the final size equation for model (4.1) is given by

where $S(0)+I(0) = N$. Figure 7 shows several trajectories with different initial conditions. Note that $S(\infty) > 0$. In other words, there exists no trajectory which converges to a trivial equilibrium point $(S(\infty), I(\infty)) = (0, 0)$.