Modelling the effective reproduction number of vector-borne diseases: the yellow fever outbreak in Luanda, Angola 2015–2016 as an example

Epidemic model for vector-borne diseases

We adopted the classic “susceptible-infected-removed” (SIR) structural framework to model both host and vector population dynamics (Gao et al., 2016; Tang et al., 2016; Saad-Roy, Van den Driessche & Ma, 2016; Zhao et al., 2018b; Earn et al., 2008; Brauer et al., 2016). We use S_h, I_h, R_h to denote the numbers of susceptible, infected and removed host population respectively. The S_v, I_v, R_v denote the numbers of susceptible, infected and removed vector population respectively. The susceptible host becomes infected by the “contact” with infectious vectors, eventually recovers, i.e., move to the recovered class, and remains protected from secondary infection. A similar path is also modelled in the vectors’ population. For simplicity, we ignored the incubation period, commonly denoted by E, of the infection in the epidemic model. Based on the above descriptions, the compartmental model is formulated in Eqs. (2). Figure 1 shows the schematic diagram of model (2). (2)${\displaystyle \left\{\begin{array}{c}{S}_h^{\prime }={\mu }_h{N}_h-{\beta }_{vh}\cdot \frac{S_h}{N_h}{I}_v-{\mu }_h{S}_h,\hfill \\ I_h^{\prime }={\beta }_{vh}\cdot \frac{S_h}{N_h}{I}_v-\left({\gamma }_h+{\mu }_h\right)I_h,\hfill \\ R_h^{\prime }={\gamma }_h{I}_h-{\mu }_h{R}_h,\hfill \\ S_v^{\prime }=B_v\left(t\right)-{\beta }_{hv}{S}_v\cdot \frac{I_h}{N_h}-{\mu }_v{S}_v,\hfill \\ I_v^{\prime }={\beta }_{hv}{S}_v\cdot \frac{I_h}{N_h}-\left({\gamma }_v+{\mu }_v\right)I_v,\hfill \\ R_v^{\prime }={\gamma }_v{I}_v-{\mu }_v{R}_v.\hfill \end{array}\right.}$ The model parameters are summarised in Table 1, and all parameters are assumed non-negative.

The N_h and N_v are the total numbers of the hosts and vectors respectively. We have

${\displaystyle N_h=S_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)=\text{constant};\text{and}}$ ${\displaystyle N_v\left(t\right)=S_v\left(t\right)+I_v\left(t\right)+R_v\left(t\right)\ne \text{constant}.}$ In our model, N_v = N_v(t) is time-dependent in a manner that is controlled by the mosquito birth rate B_v(t), namely $N_v^{\prime }=B_v\left(t\right)-{\mu }_v{N}_v$, which may not be necessarily equal to zero. To guarantee the biological reasonability, the values of all six classes should be non-negative, i.e., ⩾0. Then, S_h, I_h, R_h ⩽ N_h, and S_v, I_v, R_v ⩽ N_v. Although the vector-borne diseases could affect the vector’s lifespan in some occasions, e.g., tick-borne diseases (Niebylski, Peacock & Schwan, 1999), it is conventionally assumed the infected vectors can neither recover nor cause direct mortality in modelling studies, which is true for most of the mosquito-borne diseases, the most common group of vector-borne diseases. For simplicity, the disease-caused vector mortality is neglected in this modelling study. The inclusion of R_v in the model (2) seems unnatural, nevertheless letting γ_v = 0, thus R_v = 0, can resolve this conflict straightforwardly.

Most of the pathogens of the vector-borne diseases merely transmit via two paths, host-to-vector and vector-to-host, e.g., Dengue fever, yellow fever, Chikungunya fever and malaria, etc. For those can be transmitted via host-to-host path, commonly by sexual contact or blood transfusion, e.g., Zika fever, the host-to-host transmission is minor, and dominated by the aforementioned two paths. Therefore, we remark that although the model (2) takes a simple form, it is applicable to model the transmission dynamics of most vector-borne diseases. Model (2) also can be easily extended into more complex form, for instance the model system in Zhao et al. (2018b), Gao et al. (2016) and Tang et al. (2016). As long as the transmission paths remain from-host-to-vector and from-vector-to-host, the relationship between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ explored in this paper still holds.

Basic reproduction number

The basic reproduction number, ${\mathcal{R}}_0$, is the expected number of secondary cases produced by one typical infection joining in a completely susceptible population during its infectious period (Heffernan, Smith & Wahl, 2005). When ${\mathcal{R}}_0<1$, the disease would die out in the long run. While if ${\mathcal{R}}_0>1$, the disease would spread among the population and may cause a pandemic.

In epidemiology, the next generation matrix is a method used to derive the basic reproduction number, ${\mathcal{R}}_0$, for a compartmental model of the spread of infectious diseases. According to Van den Driessche & Watmough (2002), Diekmann, Heesterbeek & Roberts (2009) and Diekmann, Heesterbeek & Metz (1990), a systematic procedure to calculate the ${\mathcal{R}}_0$ by solving the dominant eigenvalue, i.e., the eigenvalue with the largest real part, of the next generation matrix, G, at the disease-free equilibrium (DFE). Mathematically, it can be shown easily that the DFE exists and is stable for our epidemic model (2). For the next generation matrix G = FV⁻¹, the F is the new infection (or transmission) matrix and V is the infection transfer (or transition) matrix. The entry of ith row and jth column of matrix F is denoted by F_i,j, and $F_{i,j}=\frac{\partial {\mathcal{F}}_i}{\partial x_j}$ where ${\mathcal{F}}_i$ is the ith equation of $\mathcal{F}$ and x_j is the jth variable of the vector of infected classes. For instance, in Eqs. (2), the I_h and I_v are the infected classes. Similarly, the entry of ith row and jth column of matrix V is denoted by V_i,j, and $V_{i,j}=\frac{\partial {\mathcal{V}}_i}{\partial x_j}$ where ${\mathcal{V}}_i$ is the ith equation of $\mathcal{V}$ and x_j is the jth variable of the vector of infected classes. The vector $\mathcal{F}$ is the transmission rates’ vector quantity, i.e., the changing rates from infected to non-infected classes, and vector $\mathcal{V}$ is the transition rates’ vector quantity, i.e., the changing rates among infected classes. The F is the Jacobian matrix of $\mathcal{F}$, and V is the Jacobian matrix of $\mathcal{V}$.

For the compartmental model in Eqs. (2), the I_h and I_v are the infected classes, which should be included in the vectors ($\mathcal{F}$ and $\mathcal{V}$) of infected classes. Then, we have $\mathcal{F}=\left(\begin{array}{c}\hfill {\beta }_{vh}\cdot \frac{S_h}{N_h}{I}_v\hfill \\ \hfill {\beta }_{hv}{S}_v\cdot \frac{I_h}{N_h}\hfill \end{array}\right)$, and $\mathcal{V}=\left(\begin{array}{c}\hfill \left({\gamma }_h+{\mu }_h\right)I_h\hfill \\ \hfill \left({\gamma }_v+{\mu }_v\right)I_v\hfill \end{array}\right)$. Hence, $\mathbf{F}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill {\beta }_{vh}\hfill \\ \hfill {\beta }_{hv}\cdot \frac{N_v}{N_h}\hfill & \hfill 0\hfill \end{array}\right)$, and $\mathbf{V}=\left(\begin{array}{cc}\hfill {\gamma }_h+{\mu }_h\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\gamma }_v+{\mu }_v\hfill \end{array}\right)$. The next generation matrix, G, is given as follows. ${\displaystyle \mathbf{G}=\mathbf{F}{\mathbf{V}}^{-1}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill \frac{{\beta }_{vh}}{{\gamma }_v+{\mu }_v}\hfill \\ \hfill m\cdot \frac{{\beta }_{hv}}{{\gamma }_h+{\mu }_h}\hfill & \hfill 0\hfill \end{array}\right),}$ where the term m is the vector to host ratio, which is defined by $m=\frac{N_v}{N_h}$. By solving the dominant eigenvalue of G (Van den Driessche & Watmough, 2002; Diekmann, Heesterbeek & Roberts, 2009; Cushing & Diekmann, 2016; Cushing & Diekmann, 2016), we derive the ${\mathcal{R}}_0$ of the model (2) in Eq. (3). (3)${\displaystyle {\mathcal{R}}_0=\sqrt{\left(\frac{m{\beta }_{hv}}{{\gamma }_h+{\mu }_h}\right)\cdot \left(\frac{{\beta }_{vh}}{{\gamma }_v+{\mu }_v}\right)}=\sqrt{m\cdot \frac{{\beta }_{hv}{\beta }_{vh}}{\left({\gamma }_h+{\mu }_h\right)\left({\gamma }_v+{\mu }_v\right)}}.}$ In Eq. (3), the first ratio under the square root, i.e., $\frac{m{\beta }_{hv}}{{\gamma }_h+{\mu }_h}$, represents the number of vector infections caused by one infected host, and the second, i.e., $\frac{{\beta }_{vh}}{{\gamma }_v+{\mu }_v}$, represents the number of host infections caused by one infected vector. The square root represents the geometric mean that takes the average number of secondary host (or vector) infections produced by a single infected host (or vector) (Van den Driessche, 2017).

Effective reproduction number

During an epidemic, the susceptible individuals (S_h or S_v) are gradually consumed, become infected and finally removed from the disease transmission cycle. The effective reproduction number, ${\mathcal{R}}_{\text{eff}}$, is the expected number of secondary cases produced by one typical infection joining in a population during its infectious period (Van den Driessche & Watmough, 2002). The ${\mathcal{R}}_{\text{eff}}$ is time-varying, which is denoted by ${\mathcal{R}}_{\text{eff}}\left(t\right)$ and sometimes ${\mathcal{R}}_t$ for the discretized situation (Ali, Kadi & Ferguson, 2013; Fraser, 2007; Cori et al., 2013), and ${\mathcal{R}}_{\text{eff}}$ quantifies the instantaneous transmissibility of the disease. By applying the approach in ‘Basic Reproduction Number’, the next generation matrix, G, is given by ${\displaystyle \mathbf{G}=\mathbf{F}{\mathbf{V}}^{-1}=\left(\begin{array}{cc}\hfill 0\hfill & \hfill {\beta }_{vh}\frac{S_h}{N_h}\hfill \\ \hfill {\beta }_{hv}\cdot \frac{S_v}{N_h}\hfill & \hfill 0\hfill \end{array}\right)\times \left(\begin{array}{cc}\hfill {\left({\gamma }_h+{\mu }_h\right)}^{-1}\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {\left({\gamma }_v+{\mu }_v\right)}^{-1}\hfill \end{array}\right)=\left(\begin{array}{cc}\hfill 0\hfill & \hfill \frac{S_h}{N_h}\cdot \frac{{\beta }_{vh}}{{\gamma }_v+{\mu }_v}\hfill \\ \hfill m\cdot \frac{S_v}{N_v}\cdot \frac{{\beta }_{hv}}{{\gamma }_h+{\mu }_h}\hfill & \hfill 0\hfill \end{array}\right).}$ Then, we derive ${\mathcal{R}}_{\text{eff}}$ from G. ${\displaystyle {\mathcal{R}}_{\text{eff}}=\sqrt{m\cdot \frac{{\beta }_{hv}{\beta }_{vh}}{\left({\gamma }_h+{\mu }_h\right)\left({\gamma }_v+{\mu }_v\right)}\cdot \frac{S_h}{N_h}\cdot \frac{S_v}{N_v}}={\mathcal{R}}_0\sqrt{\frac{S_h}{N_h}\cdot \frac{S_v}{N_v}},}$ where ${\mathcal{R}}_0$ is given in Eq. (3). We further define the susceptibilities of hosts (${\mathcal{S}}_{\text{h}}$) and vectors (${\mathcal{S}}_{\text{v}}$) by (4)${\displaystyle {\mathcal{S}}_{\text{h}}=\frac{S_h}{N_h},\text{and}{\mathcal{S}}_{\text{v}}=\frac{S_v}{N_v}.}$ In the epidemic model (2), we have S_h, I_h, R_h ⩽ N_h, and S_v, I_v, R_v ⩽ N_v. Therefore, $0⩽{\mathcal{S}}_{\text{h}},{\mathcal{S}}_{\text{v}}⩽1$. Henceforth, the effective reproduction number is in Eq. (5). (5)${\displaystyle {\mathcal{R}}_{\text{eff}}={\mathcal{R}}_0\sqrt{{\mathcal{S}}_{\text{h}}{\mathcal{S}}_{\text{v}}}.}$ Since $0⩽{\mathcal{S}}_{\text{h}},{\mathcal{S}}_{\text{v}}⩽1$, we have ${\mathcal{R}}_{\text{eff}}⩽{\mathcal{R}}_0$. Furthermore, (i) for most of the vector-borne diseases, the vector’s lifespan is much shorter than host’s lifespan, e.g., mosquito’s lifespan is around 7 to 60 days and tick’s lifespan is around 3 months to 3 years, and (ii) the infected vectors do not recover. These two facts will lead to an outcome in model (2) that class R_v = 0 and class I_v is extremely small and almost zero. Therefore, for simplicity, we set ${\mathcal{S}}_{\text{v}}=100\text{\%}$ for the remaining parts of this study.

Importantly, the relationship between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ in Eq. (5) holds whenever the transmission path of the vector-borne pathogens is from a host to a host via a vector.

An example of the yellow fever epidemic in Luanda 2015–2016

To illustrate the relationship between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ in Eq. (5), we compared different calculations or estimations of the effective reproduction number. We adopted the recent yellow fever (YF) outbreak in Luanda, Angola from 2015 to 2016 as a case study (Zhao et al., 2018b). The YF cases time series in Luanda and the local vaccination coverage were obtained from the situation reports released by the African Health Observatory (AHO) (WHO, 2017). Similar to the WHO (WHO, 2017) and previous literature (Kraemer et al., 2017; Zhao et al., 2018b), both probable and confirmed cases were grouped together, and were considered as the “YF cases” for further analyses (Fig. 2A).

We intend to demonstrate that the relationship between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ in Eq. (5) is more proper for pinpoint the epidemic control threshold of the vector-borne diseases than that in Eq. (1). To do so, we treat the instantaneous reproduction number, ${\mathcal{R}}_t$, calculated by the serial interval (SI) approach as the true effective reproduction number. The calculation of ${\mathcal{R}}_t$ is introduced in detail in ‘Instantaneously Reproduction Number Estimation by Renewable Equation’. After we find ${\mathcal{R}}_t$ series, we compare the calculations of the ${\mathcal{R}}_{\text{eff}}$s by using Eq. (1) or Eq. (5) to the ${\mathcal{R}}_{\text{eff}}$. The focus of this part is to compare the two relationships between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ described in Eqs. (1) and (5), and to identify which one is more proper for pinpoint the epidemic control threshold of the vector-borne diseases.

Note that although the formulation of ${\mathcal{R}}_0$ in model (2) is different from that in the YF epidemic model in Zhao et al. (2018b), the relationship between ${\mathcal{R}}_0$ and ${\mathcal{R}}_{\text{eff}}$ explored in Eq. (5) holds regardless of the complicity of the epidemic model.

The reproduction number, ${\mathcal{R}}_0\left(t\right)$ and ${\mathcal{R}}_{\text{eff}}\left(t\right)$, reconstruction approach proposed in Zhao et al. (2018b) is concisely introduced in ‘Reconstruction of the Reproduction Numbers from Compartmental Model’. The instantaneous (effective) reproduction number, $\mathcal{R}\left(t\right)$ or ${\mathcal{R}}_t$ estimation is described in details in ‘Instantaneously Reproduction Number Estimation by Renewable Equation’. We compared the two forms of effective reproduction number, ${\mathcal{R}}_{\text{eff}}\left(t\right)$, as in Eqs. (1) and (5) with the estimation by using the renewable equation in Eq. (7). To summary, Table 2 listed the relevant notations in this section as well as the remains of this study.