A viral load-based model for epidemic spread on spatial networks

1.   Introduction

Dengue is a vector-borne disease recognized as the major arbovirus in the world with approximately 500,000 dengue hemorrhagic fever (DHF) cases and 22,000 deaths ^[1,2,3]. Today, more than a third of the world's population lives in countries where dengue is endemic ^[4], with the dengue belt covering Central America, most of South America, sub-Saharan Africa, India, and South East Asia. Zika virus (ZIKV), a member of the Flavivirus family, was first isolated in Uganda in 1947 ^[5]. The first documented Zika outbreak occurred on Yap Island in the North Pacific in 2007 ^[6]. During 2013, a large-scale Zika outbreak was reported in French Polynesia ^[7]. After being spread to Brazil in 2015 ^[8], ZIKV was subsequently spread to other countries and territories, with more than 70 countries and territories being reported evidence of ZIKV transmission since 2007 by the end of January, 2017 ^[9].

Dengue can be transmitted to humans by one bite of the Aedes aegypti mosquitoes infected with one of four closely related dengue serotypes ^[2,10]. ZIKV is also primarily transmitted by the same mosquito species. However, ZIKV has been isolated in serum, saliva, urine, and semen ^{[11,12,13,14]}, confirming the possibility of non-vector transmission (through sexual contact). The first case of ZIKV transmission via sexual contact of female to male was reported in July, 2016 ^[15]. Subsequently, six more cases of sexual transmission of ZIKV in the U.S. and the first case in France were diagnosed ^[16]. Many mathematical modelling studies ^{[17,18,19,16,20,21,22,23,24,25,26]} have investigated the impact of sexual transmission on the epidemiology of Zika, with a conclusion that Zika sexual transmission can influence the magnitude of Zika outbreaks. In 2016, Gao et al. pioneered the modelling based estimation of the impact of sexual transmission on Zika epidemiology ^[16]. Similarly, by assessing the basic reproduction number, the study ^[25] evaluated the relative role of sexual transmission. In ^[24], a mathematical model was conducted to quantify the Zika prevalence between a source region to an import region. Furthermore, the modelling study ^[27] pointed out that the risk of sustained sexual transmission of Zika is underestimated, which strongly support that Zika should be classified as sexually transmitted infection.

Since cocirculation and particularly coinfection of dengue and Zika have been documented ^[28,29,30], it is natural to ask how Zika sexual transmission impacts the Zika transmission dynamics, and further how it affects the dengue transmission dynamics. Our main objective of this study is to address this issue through a mathematical model. Previously, we developed a mathematical model of co-infection of dengue and Zika with a particular focus on the potential impact of vaccination against dengue for Zika outbreak ^[31]. In our current work, we extend the dengue-Zika coinfection dynamic model by incorporating the natural birth and death of humans. We further expand it by including the Zika sexual transmission routes among humans. To our best knowledge, our work here is the first attempt to develop a mathematical model to address how Zika sexual transmission affects the dynamics of both Zika and dengue.

The manuscript is organized as follows. In the next section, we first develop a dengue-Zika coinfection model involving the sexual transmission of Zika among humans. In section 3, we analyse the dynamics of the proposed model, including the global dynamics of the two submodels with only one disease being considered. In section 4, we discuss the impact of Zika sexual transmission on dengue endemic through cascade effect analysis. We then investigate how Zika sexual transmission affects the transmission dynamics of both dengue and Zika through numerical simulations in section 5. Finally, we make some concluding remarks and discussion in section 6.

2.   Model formulation

Mosquito population $N_m$ is divided into four classes: susceptible ($S_m$), infected with dengue only ($I_{md}$), infected with Zika only ($I_{mz}$), and infected with dengue and Zika ($I_{mdz}$). We use a SI-type of structure for the coinfection of dengue and Zika in mosquito population. Thus the model for mosquitoes is given by

Human population is stratified into: susceptible to dengue and Zika ($S$), dengue infected but susceptible to Zika ($I_d$), Zika infected but susceptible to dengue ($I_z$), dengue and Zika co-infected ($I_{dz}$), recovered from dengue and susceptible to Zika ($R_d$), recovered from Zika and susceptible to dengue ($R_z$), Zika infected but immune to dengue ($J_d^z$), dengue infected but immune to Zika ($J_z^d$), recovered from both dengue and Zika ($R_{dz}$). $N_h = S+I_d+I_z+I_{dz}+R_d+R_z+J_d^z+J_z^d+R_{dz}$ denotes the total number of humans.

We assume a SIR-type model for the co-infection of dengue and Zika among humans, which gives

Here, $\beta_s$ is the sexual transmission coefficient from humans infected with Zika only ($I_z$) to humans who are susceptible to both viruses ($S$), $a_{s1}$ and $a_{s2}$ denote the relative human-to-human sexual transmissibility of humans infected with both virus ($I_{dz}$) and humans with Zika infection but dengue immune ($J_d^z$), respectively, compared with the humans infected with Zika only ($I_z$). Similarly, $\beta_{1s}$ is the sexual transmission coefficient from humans with Zika only ($I_z$) to human with dengue only ($I_d$), with $a_1$ and $a_2$ being the relative human-to-human transmissibility corresponding to $I_{dz}$ and $J_d^z$, respectively. $\beta_{rs}$ represents the sexual transmission coefficient from humans with Zika only ($I_z$) to humans recovered from dengue but susceptible to Zika ($R_d$), with $a_{r1}$ and $a_{r2}$ being the relative human-to-human transmissibility of $I_{dz}$ and $J_d^z$, respectively. $\Lambda_h$ is the constant birth rate of humans, $\mu_h$ is the mortality rate of humans while $1/\mu_h$ can be seen as the life span of sexual activity for humans ^[19]. The transmission diagram is shown in Figure 1. The definitions of the variables are listed in Table 1 while the definitions and values for the other parameters of model (2.1)–(2.2) can refer to Table 2.

3.   Dynamics of disease-free equilibria

In this section, we mainly discuss the existence and stability of the disease-free equilibria for system (2.1)–(2.2). Define

Obviously, $\mathcal{D}$ is a positively invariant and attracting region in $\mathbb{R}_+^{13}$ for system (2.1)–(2.2).

It is easy to see that system (2.1)–(2.2) has a disease-free equilibrium, which is given

Using the next generation matrix method introduced in papers ^[37,38], the basic reproduction number for system (2.1)–(2.2) is calculated as $\mathcal{R}_0 = \max\{\mathcal{R}_d, \mathcal{R}_z^s\}$ (see Appendix A for more details), where

are the dengue basic reproduction number and the Zika basic reproduction number, respectively. Here

are the basic reproduction numbers of sexual transmission and vectorial transmission for Zika, respectively. In summary, we have established

Theorem 3.1. The disease-free equilibrium $E_0$ is locally asymptotically stable if $\mathcal{R}_0 < 1$ and unstable if $\mathcal{R}_0 > 1$.

Next, we will have an insight into the dynamics of the disease-free equilibrium if only one disease circulates among humans and mosquitoes. Firstly, we assume that only dengue occurs, then system (2.1)–(2.2) can be simplified to the following subsystem

For system (3.3), $\mathcal{D}_1 = \{(S, I_d, S_m, I_{md})\in \mathbb{R}_+^4|0 < S+I_d\leq\frac{\Lambda_h}{\mu_h}, 0 < S_m+I_{md}\leq\frac{\Lambda_m}{\mu_m}\}$ is a positively invariant and attracting region in $\mathbb{R}_+^4$. Correspondingly, the disease-free equilibrium reduces to $E_0^d = (\frac{\Lambda_h}{\mu_h}, 0, \frac{\Lambda_m}{\mu_m}, 0)$, which is locally asymptotically stable when $\mathcal{R}_d < 1$ with $\mathcal{R}_d$ being given in equation (3.1). Note that $\mathcal{R}_d$ is also the basic reproduction number of model (3.3) ^[32,39]. In fact, $\mathcal{R}_d < 1$ can ensure that $E_0^d$ is globally stable, that is,

Theorem 3.2. If $\mathcal{R}_d < 1$, the disease-free equilibrium $E_0^d$ of system (3.3) is globally asymptotically stable in $\mathcal{D}_1$.

Note that, the similar models of dengue have been studied in ^[40,32]. We can use the same methods in these two studies to prove Theorem 3.2 (see Appendix B for details).

Performing the similar process by assuming that only Zika circulates among the humans and mosquitoes, system (2.1)–(2.2) becomes

The positively invariant and attracting region should now be defined as $\mathcal{D}_2 = \{(S, I_z, S_m, I_{mz})\in \mathbb{R}_+^4|0 < S+I_z\leq\frac{\Lambda_h}{\mu_h}, 0 < S_m+I_{mz}\leq\frac{\Lambda_m}{\mu_m}\}$ in $\mathbb{R}_+^4$. Also, system (3.4) has a disease-free equilibrium $E_0^z = (\frac{\Lambda_h}{\mu_h}, 0, \frac{\Lambda_m}{\mu_m}, 0)$. It is locally asymptotically stable if $\mathcal{R}_z^s < 1$ and unstable if $R_z^s > 1$, here, $\mathcal{R}_z^s$ is the basic reproduction number of system (3.4) with the formula being given in equation (3.1). The global dynamics of $E_0^z$ can be concluded as in the following.

Theorem 3.3. If $\mathcal{R}_z^s < 1$, the disease-free equilibrium $E_0^z$ of system (3.4) is globally asymptotically stable in $\mathcal{D}_2$.

Proof. Similar to Theorem 3.2, if $\mathcal{R}_z^s < 1$, there must exist an $\varepsilon_2 > 0$ such that $c\eta_z\frac{\mu_h\Lambda_m}{\mu_m(\gamma_z+\mu_h)\Lambda_h} \frac{1}{1-\frac{\beta_s}{\gamma_z+\mu_h}} < \varepsilon_2 < \frac{\mu_m}{c\beta_z}$. We introduce a function $L_2 = \varepsilon_2 I_z+I_{mz}$, which satisfies $L_2\geq 0$ along the solution of system (3.4) with $L_2 = 0$ if and only if both $I_z$ and $I_{mz}$ are zero. The derivative of $L_2$ along the solution of system (3.4) satisfies

since $\frac{S}{N_h}\leq1$ and $S_m\leq\frac{\Lambda_m}{\mu_m}$ hold true in $\mathcal{D}_2$ and the assumptions $\mathcal{R}_s < 1$, $c\eta_z\frac{\mu_h\Lambda_m}{\mu_m(\gamma_z+\mu_h)\Lambda_h} \frac{1}{1-\frac{\beta_s}{\gamma_z+\mu_h}} < \varepsilon_2 < \frac{\mu_m}{c\beta_z}$ are satisfied. Further, $L_2' = 0$ if and only if $I_z = 0$ and $I_{mz} = 0$. Therefore, according to the LaSalle's Invariance Principle, we have that $E_0^z$ is globally asymptotically stable in $\mathcal{D}_2$.

4.   Cascade effect of Zika transmission on dengue endemic

In this section, we will illustrate how the Zika sexual transmission can first lead to Zika endemic, and then lead to dengue endemic through cascade effect analysis. The endemic equilibrium of Zika only model $E_z(\hat{S}, \hat{I}_z, \hat{R}_z, \hat{S}_m, \hat{I}_{mz})$ satisfies the following conditions:

Solving the above equations, there are

where $\hat{I}_z$ is the positive root of the following equation

with $A = \frac{\mu_m\Lambda_h}{c\eta_z\mu_h} > 0, \ B = \frac{\Lambda_h}{\gamma_z+\mu_h} > 0$ and $\mathcal{R}_s, \ \mathcal{R}_z$ being given in (3.2).

Denote $\Delta = [A(\mathcal{R}_s+{\mathcal{R}_z}^2)+B(1-\mathcal{R}_s)]^2-4AB\mathcal{R}_s(1-\mathcal{R}_s-{\mathcal{R}_z}^2)$, then

Through some straightforward analysis, we find that $f(I_z) = 0$ has only one positive root given by

if and only if $1-\mathcal{R}_s-{\mathcal{R}_z}^2 < 0$ which is equivalent to $\mathcal{R}_z^s > 1$. That is, the endemic equilibrium of Zika only model exits if and only if $\mathcal{R}_z^s > 1$. Therefore, Zika can become endemic by including the sexual transmission because $\mathcal{R}_z^s$ is increasing as $\beta_s$ increases, for the area that the vectorial transmission only can not (i.e. $\mathcal{R}_z < 1$).

We then consider the transmission dynamics of dengue after Zika reaching its endemic state. That is, we assume that only dengue circulates among humans and mosquitoes and ignore all the Zika transmission, but we set the initial conditions as the Zika endemic state. Then system (2.1)–(2.2) can be reduced to the following model:

Note that, now the disease-free equilibrium of dengue should be $E_d^*(\hat{S}, 0, \hat{I}_z, 0, \hat{R}_z, 0, 0, \hat{S}_m, 0, \hat{I}_{mz})$.

In order to investigate the stability of $E_d^*$, we also use the next generation operator. Thus, we just need to consider the compartments $I_d$, $J_z^d$ and $I_{md}$ with

Therefore, we have

Through easy calculation, we then get

with

Then by solving $|\lambda I-FV^{-1}| = 0$, there is

Therefore, we can define a reproduction number of dengue with Zika being endemic which is called invasion reproduction number ^[41] as

Thus, the disease-free equilibrium of dengue $E_d^*$ is stable if $\mathcal{R}_d^* < 1$ and unstable if $\mathcal{R}_d^* > 1$. In other words, dengue can further become endemic after Zika being endemic if $\mathcal{R}_d^* > 1$.

Then we discuss how the Zika sexual transmission will affect dengue transmission dynamics when $\mathcal{R}_d < 1$. It is easy to see that if $\mathcal{R}_z^s < 1$, there can be $\hat{S} = \Lambda_h/\mu_h$, $\hat{S}_m = \Lambda_m/\mu_m$, and $\hat{I}_z = \hat{I}_{mz} = \hat{R}_z = 0$, hence $\mathcal{R}_d^* = \mathcal{R}_d$. This means, Zika sexual transmission can not lead to dengue to be endemic before invoking a Zika endemic. Further, if we set $\beta_{rd} = \beta_{1d} = \beta_d$, $\eta_{1d} = \eta_{jd} = \eta_{1jd} = \eta_d$ and $\gamma_d = \gamma_z^d$, there is also $\mathcal{R}_d^* = \mathcal{R}_d < 1$ with $\hat{S}+\hat{I}_z+\hat{R}_z = N_h = \Lambda_h/\mu_h$ and $\hat{S}_m+\hat{I}_{mz} = \Lambda_m/\mu_m$. Therefore, under this situation, the sexual transmission of Zika can not lead dengue to be endemic as well.

It is worth mentioning that the phenomenon of antibody dependent enhancement (ADE) usually occurs among the different serotypes of dengue, and between dengue and Zika. There are studies^{[42,43,44,45]} showing that previous exposure to one dengue serotype can increase the risk of the infection by a second serotype. The similar results that plasma to dengue can enhance the infection of Zika have been reported in the clinical studies ^[46,47,48]. Furthermore, ADE impact between dengue and Zika is actually bidirectional ^[49]. Therefore, we can assume that there is the ADE of Zika for dengue infection, that is we assume that $\beta_{rd} > \beta_d$. Also, we let $\beta_{1d} = \beta_d$, $\eta_{1d} = \eta_{jd} = \eta_{1jd} = \eta_d$ and $\gamma_d = \gamma_z^d$. Then, we have that $\mathcal{R}_d^*$ is increasing as $\hat{R}_z$ increases. Furthermore, $\hat{R}_z$ is increasing as $\beta_s$ increases (see Appendix C for more details). That is, $\mathcal{R}_d^*$ is increasing as $\beta_s$ increases. This meas that Zika sexual transmission can first induce a Zika endemic, and then can a dengue endemic if there is ADE of Zika for dengue infection.

5.   Numerical simulations

In this section, through numerical analyses, we investigate how sexual transmission of Zika affects the dynamics of both dengue and Zika. First, in order to determine the most important parameters for impacting the Zika basic reproduction number, we explored the parameter space by performing an uncertainty analysis using a Latin hypercube sampling method. We used a partial rank correlation coefficients (PRCCs) ^[50,51,52] to examine the sensitivity analysis for the Zika basic reproduction number with respect to all the parameters involved in $\mathcal{R}_z^s$. In the absence of data to inform distribution functions, we chose a uniform distribution for all input parameters. The PRCC values for $\mathcal{R}_z^s$ is shown in Figure 2. Figure 2 shows that the first three parameters with most impact on $\mathcal{R}_z^s$ are the mosquito mortality rate $\mu_m$, the biting rate $c$ and the mosquito recruitment rate $\Lambda_m$, and the sexual transmission coefficient is positively correlated to Zika basic reproduction number with a relatively large PRCCs.

As mentioned in the last section, we can assume that there is the ADE of Zika for dengue infection, and we use a parameter $\kappa$ to describe the multiplication factor of the susceptibility to dengue induced by ADE. Thus, we set $\beta_{rd} = \kappa\beta_d$, $\beta_{rzd} = \kappa p \beta_{dz}$. And we assume no ADE of dengue for Zika infection with $\beta_{rz} = \beta_{z}$, $\beta_{rdz} = (1-p)\beta_{dz}$. In order to perform the simulations, we mainly let the parameters $\beta_s$ and $\kappa$ change and fixed the other parameter values. Table 1 has given an overview of the setup of the basic parameters. The other mosquito-to-humans transmission probabilities are assumed as: $\beta_{1d} = \beta_d$, $\beta_{1zd} = p\beta_{dz}$, $\beta_{1z} = \beta_z$ and $\beta_{1dz} = (1-p)\beta_{dz}$, and the rest parameter values are fixed as:

By fixing $\kappa = 1$, we plotted the solutions of system (2.1)–(2.2) with various values of $\beta_s$, as shown in Figure 3. From Figure 3(A) and (C), we find that if there is no sexual transmission of Zika (i.e. $\beta_s = 0$), the basic reproduction number values are $\mathcal{R}_d = \mathcal{R}_z^s = 0.916 < 1$ and the disease-free equilibrium $E_0$ is stable. When the sexual transmission of Zika is included, it can make the Zika only endemic equilibrium $E_z$ become stable by increasing the Zika basic reproduction number above the threshold 1, as shown in Figure 3(B) and (D). The similar results can also be obtained from the simulation of Zika only model (i.e. model (3.4)), as shown in Figure 4. Furthermore, it follows from Figure 4 that the number of humans infected with Zika can stabilize at a higher level as $\beta_s$ increases. This means that the endemic level of Zika is increasing as Zika sexual transmission coefficient increases. In Figure 3 we fixed $\kappa = 1$, that is, we assumed that there is no ADE of Zika for dengue infection. Under this scenario, sexual transmission of Zika does not have a significant effect on the dynamics of dengue while it always dies out finally.

Further, by letting $\kappa = 3$, we plotted the solutions of system (2.1)–(2.2) with various values of Zika sexual transmission coefficient $\beta_s$ in Figure 5. It follows from Figure 5 that the disease-free equilibrium $E_0$ is globally stable when $\beta_s = 0$. However, if we include the Zika sexual transmission by letting $\beta_s = 0.8$, then system (2.1)–(2.2) has a stable endemic equilibrium. Note that, when $\beta_s = 0.8$, $\mathcal{R}_d$ is also equal to 0.916 which is less than the threshold 1. However, as we can see from Theorem 3.2, dengue always dies out whenever $\mathcal{R}_d < 1$ if no Zika circulates within this area. This means that, due to the cocirculation of dengue and Zika and ADE of Zika for dengue infection, for an area where vector transmission only is not enough to make dengue and Zika be endemic, Zika sexual transmission can not only lead Zika to be endemic, but also make dengue become endemic even though the dengue basic reproduction number is less than unit.

6.   Discussion and conclusion

Recently, there is increasing evidence confirming that sexual transmission of Zika occurs among humans. Several studies have investigated how Zika sexual transmission affects the dynamics of the spread of Zika via mathematical modelling. In this work, we developed a dengue-Zika coinfection model to study the impact of Zika sexual transmission on the transmission dynamics of both dengue and Zika when the two diseases cocirculate within a same area.

The basic reproduction number for both dengue and Zika, $\mathcal{R}_d$ and $\mathcal{R}_z^s$, is determined, and $\mathcal{R}_0 = \max\{\mathcal{R}_d, \mathcal{R}_z^s\} < 1$ is the necessary and sufficient condition for the local stability of the disease-free equilibrium $E_0$. Theoretically, we have proved that the disease-free equilibrium for the dengue only model ($E_0^d$) or Zika only model ($E_0^z$) is globally stable if and only if their basic reproduction number is less than 1. Correspondingly, if the Zika (or dengue) basic reproduction number exceeds the threshold 1, then a Zika (or dengue) only endemic equilibrium for system (2.1)–(2.2) appears.

Through cascade effect analysis and numerical simulations, we try to have a further insight into the impact of Zika sexual transmission on the dynamics of dengue and Zika. The cascade effect analysis shows that Zika sexual transmission can not lead dengue to be endemic before making Zika be endemic, but it can first make Zika be endemic, and then lead dengue to be endemic if there is ADE of Zika for dengue infection. Furthermore, our sensitivity analysis shows that the first three parameters with significant impact on the Zika basic reproduction number are the Zika sexual transmission coefficient, the mortality rate of mosquitoes, and the Zika recovery rate of humans. $\mathcal{R}_z^s$ increases if the sexual transmission route is included. As a result, the sexual transmission of Zika can lead Zika to be endemic in an area where Zika will die out if there is only the vector-borne transmission. We further found that the epidemic level of Zika increases as Zika transmission coefficient increases.

As mentioned above, when $\mathcal{R}_d < 1$, the disease-free equilibrium for the dengue only model is globally stable. That is, if the dengue basic reproduction number is less than 1, dengue always goes to extinction if no Zika circulates in this area. In Figure 3, we showed that Zika sexual transmission may not affect the dynamics of dengue significantly when there isn't ADE of Zika for dengue infection. However, if we assume that there is ADE of Zika for dengue infection, increasing Zika sexual transmission coefficient can make dengue and Zika be endemic simultaneously, as shown in Figure 5. Note that, here the dengue basic reproduction number is equal to 0.916 which is less than 1. In the current study, we take $1/\mu_h$ as the lifespan of sexual activity for humans as we focus on the effect of sexual transmission on the dynamics of both dengue and Zika in sexual active group. The ignorance of sexual inaction group will disregard the contribution of sexual inactive people to the vector-borne transmission and hence may underestimate disease infections. To estimate the contribution of sexual inactive group on the vector-borne transmission, we plotted the solutions of model (3.4) in Figure 6 by changing the death rate and the recruitment rate of humans and fixing all the other parameters. In Figure 6, we set the lifespan of human as 35 years (taking $1/\mu_h = 35$ years as the life span of sexual activity for humans) and 70 years (taking $1/\mu_h = 70$ years as the natural death rate), respectively. Correspondingly, we assume that the recruitment rate of the total population (including the sexual active and inactive groups) is 1.3 or 1.6 times to the recruitment rate of the sexual active group. We obtain that the values of $I_z$ at the endemic equilibrium are $20.16$ for $\mu_h = 1/(35\times 365), \Lambda_h = 2.2$, $22.43$ for $\mu_h = 1/(70\times 365), \Lambda_h = 2.2\times 1.3$, and $25.85$ for $\mu_h = 1/(70\times 365), \Lambda_h = 2.2\times 1.3$. Thus, we can estimate the contribution of the sexual inactive group to the vector-borne transmission is about 10 percent with 1.3 times recruitment rate and 22 percent for 1.6 times recruitment rate. There need further study to address this issue in more details.

Our model captures some important results that Zika sexual transmission can not only lead Zika to be endemic, but also further make dengue endemic with an assumption that there is ADE of Zika for dengue infection for the area both the diseases will die out with vectorial transmission only. This means that, due to the coinfection and ADE, increasing the transmission probability (or the endemic level) of one disease can lead to the persistent of the other disease and increase its endemic level. Note that, we set $\beta_{dz} = \beta_d = \beta_z$ in Figure 5. However, due to the competition on two pathogens within a host, the transmission probability of the co-infection class may be less than the probability of the classes infected of one disease. Thus, we have done further numerical experiments by setting $\beta_{dz} = 0.01 < \beta_d = \beta_z = 0.15$ and fixed all the other parameter values as those in Figure 5, and we get the similar conclusion to those shown in Figure 5. This implies that the choice of parameter $\beta_{dz}$ may not be the unique reason for inducing the interesting result. Hence a deep study will be needed to identify the key factors that promote the transmission of both pathogens when they are co-circulating within the same area, and we leave this for future study. For another aspect, ADE between dengue and Zika is usually bidirectional, and our model is easy to be modified considering this bidirectional effect. Another non-vector borne transmission, vertical transmission, has also been reported ^[53]. How it will affect the dynamics of both dengue and Zika is still unclear. Addressing this issue would require a further indepth research, and will be the goal of future work. In our model, we did not include the exposed classes for both humans and mosquitos, which can have a significant effect on the transmission dynamics of vector-borne diseases ^[16,54,55]. It was known that ignoring the extrinsic incubation period in mosquitoes may lead to overestimation of the infection risk. Hence our conclusion obtained here may overestimate the disease infection, but we hope the approaches we used for investigation of coinfection are able to be applied more generally. The higher dimensional model or delay differential equations will be formulated to model the coinfection of dengue and Zika in future study.

Acknowledgments

This project was supported by the National Natural Science Foundation of China (NSFC, 11631012, 11571273(YX)), and by the Canada Research Chair Program (230720 (JW)), the Natural Sciences and Engineering Research Council of Canada (105588-2011 (JW)).

Conflict of interest

The authors declare no conflict of interest.

Supplementary

Appendix A. Calculation of $\mathcal{R}_0$ for system (2.1)–(2.2)

The basic reproduction number is determined using the next generation operator. Considering the equations of $I_d, I_z, I_{dz}, J_d^z, J_z^d, I_{md}, I_{mz}, I_{mdz}$, we have

and

Then, we obtain

with $\Theta = \dfrac{\mu_h\Lambda_m}{\mu_m\Lambda_h}, \; C_1 = c(1-q), \; C_2 = c(1-p)$, and

with $v_{33} = \gamma_{dz}^d+\gamma_{dz}^{z}+\mu_h$. Thus we have

with

Therefore,

where

By simple calculating, the corresponding characteristic equation is

from which we have that the dominant eigenvalues of $FV^{-1}$ are

Therefore, the basic reproduction number for system (2.1)–(2.2) is the spectral radius of $FV^{-1}$, that is $\mathcal{R}_0 = \rho(FV^{-1}) = \max\{\mathcal{R}_d, \mathcal{R}_z^s\}$, where

is the basic reproduction number of dengue, and

is the basic reproduction number of Zika with

being the basic reproduction numbers of sexual transmission and vectorial transmission for Zika, respectively.

Appendix B. Proof of Theorem 3.2.

Proof. If $\mathcal{R}_d < 1$, which means $\frac{c\beta_d}{\mu_m}\frac{c\eta_d\mu_h\Lambda_m}{\mu_m(\gamma_d+\mu_h)\Lambda_h} < 1$, there must exist an $\varepsilon_1 > 0$ satisfying $\frac{c\eta_d\mu_h\Lambda_m}{\mu_m(\gamma_d+\mu_h)\Lambda_h} < \varepsilon_1 < \frac{\mu_m}{c\beta_d}$. Define a function $L_1 = \varepsilon_1 I_d+I_{md}$. It satisfies that $L_1\geq 0$ along the solution of system (3.3) with $L_1 = 0$ if and only if $I_d$ and $I_{md}$ are both equal to zero. The derivative of $L_1$ along the solution of system (3.3) satisfies

due to the fact that $\frac{S}{N_h}\leq1$, $S_m\leq\frac{\Lambda_m}{\mu_m}$ in $\mathcal{D}_1$ and the assumption $\frac{c\eta_d\mu_h\Lambda_m}{\mu_m(\gamma_d+\mu_h)\Lambda_h} < \varepsilon_1 < \frac{\mu_m}{c\beta_d}$. Further, we have $L_1' = 0$ if and only if $I_d = 0$ and $I_{md} = 0$. Thus according to the LaSalle's Invariance Principle, we obtain that $E_0^d$ is globally asymptotically stable in $\mathcal{D}_1$.

Appendix C. $\hat{S}$ and $\hat{S}_m$ are decreasing, and $\hat{I}_z$, $\hat{R}_z$ and $\hat{I}_{mz}$ are increasing as $\beta_s$ increases if $\mathcal{R}_z^s > 1$.

Regarding $\hat{I}_z$ as a function of $\beta_s$, we can prove that it is increasing as $\beta_s$ increases whenever it exists. It follows from the formula of $\hat{I}_z$ that only $\mathcal{R}_s$ includes the parameter $\beta_s$, and it is a monotonous increasing function of $\beta_s$. Thus we just need to verify the monotonicity of $\hat{I}_z$ with respect to $\mathcal{R}_s$ instead. Calculating the derivative of $\hat{I}_z$ with respect to $\mathcal{R}_s$, there is

Definitely there are $A{\mathcal{R}_z}^2+B > 0$ and $\frac{1}{\sqrt{\Delta}} > 0$, further if $(B-A{\mathcal{R}_z}^2)(A\mathcal{R}_s+B\mathcal{R}_s+A{\mathcal{R}_z}^2-B)- 4AB\mathcal{R}_s{\mathcal{R}_z}^2\geq0$, there must be $\frac{d\hat{I}_z}{d\mathcal{R}_s} > 0$, which means that $\hat{I}_z$ is an increasing function of $\mathcal{R}_s$. Otherwise, when $(B-A{\mathcal{R}_z}^2)(A\mathcal{R}_s+B\mathcal{R}_s+A{\mathcal{R}_z}^2-B)-4AB\mathcal{R}_s{\mathcal{R}_z}^2 < 0$, we can verify that $A{\mathcal{R}_z}^2+B+\frac{1}{\sqrt{\Delta}} [(B-A{\mathcal{R}_z}^2)(A\mathcal{R}_s+B\mathcal{R}_s+ A{\mathcal{R}_z}^2-B)-4AB\mathcal{R}_s{\mathcal{R}_z}^2] > 0$ due to $(A{\mathcal{R}_z}^2+B)^2\Delta-\left[(B-A{\mathcal{R}_z}^2)(A\mathcal{R}_s+B\mathcal{R}_s+A{\mathcal{R}_z}^2-B)- 4AB\mathcal{R}_s{\mathcal{R}_z}^2\right]^2 = (A+B)^2{\mathcal{R}_s}^2 > 0.$ As a conclusion, there is always that $\hat{I}_z$ is increasing as $\mathcal{R}_s$ increases, hence increasing as $\beta_s$ increases. Then, it follows from the relationships of $\hat{S}$, $\hat{R}_z$, $\hat{I}_{mz}$ and $\hat{S}_m$ that $\hat{S}$ and $\hat{S}_m$ are decreasing as $\beta_s$ increases, and $\hat{R}_z$ and $\hat{I}_{mz}$ are increasing as $\beta_s$ increases.