Analyzing the control of dengue by releasing Wolbachia-infected male

Abstract: To date, an innovative strategy to control dengue is to release Wolbachia-infected male mosquitoes into wild areas to sterilize wild female mosquito vectors by cytoplasmic incompatibility (CI). To investigate the efficacy of Wolbachia in blocking dengue virus transmission, we develop a deterministic mathematical model of human and mosquito populations in which one dengue serotype circulates. The delay differential equation model captures the respective extrinsic and intrinsic incubation periods (EIP and IIP) in the mosquito and human, as well as the maturation delay between mating and emergence of adult mosquitoes, which have received relatively little attention. We analyze the existence and stability of disease-free equilibria, and obtain a sufficient and necessary condition on the existence of the disease-endemic equilibrium. We also determine two threshold values of the release ratio θ, denoted by θ∗ 1 and θ ∗ 2 with θ ∗ 1 > θ ∗ 2. When θ > θ ∗ 1, the mosquito population will be eradicated eventually. When θ∗ 2 < θ < θ ∗ 1, a complete mosquito eradication becomes impossible, but virus eradication is ensured at the meantime. When θ < θ∗ 2, the disease-endemic equilibrium emerges that allows dengue virus to circulate between humans and mosquitoes. We carry out sensitivity analysis of the threshold values in terms of the model parameters, and simulate several possible control strategies with different release ratios, which confirm the public awareness that reducing mosquito bites and killing adult mosquitoes are the most effective strategy to control the epidemic. Our model provides new insights on the effectiveness of Wolbachia in reducing dengue at a population level.


Introduction
Dengue is a viral disease transmitted mainly by Aedes mosquitoes, including Aedes aegypti and Aedes albopictus. Infected human beings carrying dengue viruses may get high fever, headaches, and skin rash, which may progress to dengue shock syndrome or dengue hemorrhagic fever. The annually reported dengue cases increased sharply from 2.2 million in 2010 to 3.2 million in 2015 [1]. More than 500,000 people with severe dengue are hospitalized annually, and the case fatality percentage is about 2.5%. Due to the lack of commercially available vaccines for dengue control, traditional methods have been focused on vector control by heavy applications of insecticides and environmental management [2]. However, these programs have not prevented the spread of these diseases due to the rapid development of insecticide resistance [3] and the continual creation of ubiquitous breeding sites.
An innovational mosquito control method utilizes Wolbachia, an intracellular bacterium present in about 60% of insect species [4], including some mosquitoes. Wolbachia can induce cytoplasmic incompatibility (CI) in mosquitoes, which results in early embryonic death from matings between Wolbachiainfected males and females that are either uninfected or harbor a different Wolbachia strain [5,6]. In the world's largest "mosquito factory" located in Guangzhou, China, 20 million Wolbachia-infected males are produced each week. These mosquitoes have been released in several urban and suburb areas since 2015, and killed more than 95% mosquitoes on the Shazai island [7]. Motivated by the success and the challenge in field trials, the study on the complex Wolbachia spread dynamics has become a hot research topic. Various mathematical models have been developed, including models of ordinary differential equations [8][9][10], delay differential equations [9,[11][12][13], impulsive differential equations [14], stochastic equations [15], and reaction-diffusion equations [16]. These models are mainly devoted to analyzing the threshold dynamics for population replacement with Wolbachia-infected female mosquitoes or CI-driven population suppression, which only included mosquito population into the models. However, it is usually a formidable task for complete population replacement or eradication in large-scale field trials, and we can only concede mosquito replacement/eradication to virus eradication by restricting the mosquito densities below the epidemic risk threshold.
To assess the efficacy of blocking dengue virus transmission by Wolbachia, a deterministic mathematical model of human and mosquito populations interfered by the circulation of a single dengue serotype was developed in the framework of SIER model [17]. This important study has inspired further development of compartmental models to analyze the transmission dynamics of dengue [18][19][20][21][22]. In this paper, we extend their efforts by incorporating the Wolbachia-infected male mosquito release into a compartmental model, where the infected males are maintained at a fixed ratio to the adult female mosquito population. We divide the mosquito population into three compartments: susceptible, exposed and infectious, and divide the human population into four compartments: susceptible, exposed, infectious and recovered. In particular, we include the respective extrinsic and intrinsic incubation periods (EIP and IIP) in the mosquito and human populations in our model. These periods have been shown to be crucial in clinical diagnosis, outbreak investigation, and dengue control [23], but have received relatively little attention [9,18]. Further, the release of Wolbachia-infected male mosquitoes is not effective immediately in sterilizing wild females due to the maturation delay between mating and emergence of adults [24]. To characterize the effect of EIP, IIP and the maturation delay on dengue transmission, we introduce three delays into our model.
The model treats IIP in humans, EIP in mosquitoes, and the maturation delay as three delays. We aim to find the threshold values of the release ratio for mosquito or virus eradication under the proportional release policy. By analyzing the existence and stability of disease-free equilibria, we obtain the sufficient and necessary condition on the existence of the disease-endemic equilibrium. Two threshold values of the release ratio θ, denoted by θ * 1 and θ * 2 with θ * 1 > θ * 2 are explicitly expressed.
When θ > θ * 1 , the mosquito population will be eradicated eventually. If it fails for mosquito eradication but θ * 2 < θ < θ * 1 , virus eradication is ensured together with the persistence of susceptible mosquitoes. When θ < θ * 2 , the disease-endemic equilibrium emerges that allows dengue virus to circulate between humans and mosquitoes through mosquito bites. Sensitivity analysis of the threshold values in terms of the model parameters, and numerical simulations on several possible control strategies with different release ratios confirm the public awareness that reducing mosquito bites and killing adult mosquitoes are the most effective strategy to control the epidemic. Our model will provide new insights on the effectiveness of Wolbachia in reducing dengue at a population level.

Model formulation
We divide adult mosquitoes into subgroups of susceptible (S M ), exposed (E M ) and infectious (I M ). Let N M = S M + E M + I M be the total number of mosquitoes. The human population is divided into four subpopulations: susceptible (S H ), exposed (infected but not infectious, E H ), infectious (I H ), and recovered (R H ). The total number of humans is denoted by N H = S H + E H + I H + R H . Without the bothering of the dengue virus, i.e., E M = I M = E H = I H = R H = 0, we assume that mosquito and human populations follow the logistic growth [25] satisfying where Let τ H be the intrinsic incubation period (IIP) in humans between infection and the onset of infectiousness. IIP is an important determinant of dengue transmission dynamics, which varies from 3 to 14 days. Exposed humans become infectious at the rate With further assumption that infected but not infectious humans have the same death rate as that of susceptible humans, E H follows Let γ be the recovery rate of humans. Then for I H and R H , we arrive at where we disregard the negligible dengue mortality in humans [26] and set an identical mortality rate for humans.
With the release of Wolbachia-infected male mosquitoes, denoted by R(t) at time t, the birth rate of mosquitoes is reduced from r M to r M (1 − Probability of complete CI occurrence).
Under random mating and equal sex determination, the probability of CI occurrence at time t is the proportion of Wolbachia-infected male mosquitoes among all male mosquitoes, i.e., R(t)/(R(t) + N M (t)). However, there is a delay between the release of Wolbachia-infected males and the reduction of the wild mosquitoes which is caused by the maturation delay between mating and emergence of adult mosquitoes [24], denoted by τ e . Hence, the number of susceptible mosquitoes that survive the maturation period is The susceptible mosquitoes become exposed at a rate of [bβ HM S M (I H /N H )], where β HM is the fraction of virus transmission from infectious humans to susceptible mosquitoes through blood meals. Taking these considerations into account, we have Exposed mosquitoes spend the extrinsic incubation period (EIP), denoted by τ M , which is the viral incubation period between the time when a female mosquito takes a viraemia blood meal from an infectious human and the time when that mosquito becomes infectious [27], typically 8∼12 days. EIP has been frequently recognized as a crucial component of dengue virus transmission dynamics [23]. In view of this point, the rate that exposed mosquitoes become infectious per unit time is Then the whole set of equations in our model ends with for exposed and infectious mosquitoes, respectively. Our purpose is to find the threshold values in terms of the release ratio for mosquito or virus eradication under the proportional release policy, where the Wolbachia-infected male mosquitoes is maintained in a fixed proportion to the adult female mosquito population, i.e., (2.8) The explicit expressions of the threshold values are obtained as follows: We get the following theorem. Theorem 1. If θ > θ * 1 , then eradication of mosquitoes occurs. If θ * 2 < θ < θ * 1 , then eradication of virus occurs. If θ < θ * 2 , then there exists a unique disease-endemic equilibrium of system (2.1)-(2.7). The proof of Theorem 1 is embodied in the analysis of the existence and stability of equilibrium points of system (2.1)-(2.7) in the following two sections, and we omit it here.

Existence of equilibria
To determine the steady-state solutions of system (2.1)-(2.7), we set the right sides of (2.1)-(2.7) to zero and ignore the time lags to arrive at Adding (3.1a) to (3.1d) together, we have provided that µ H < r H holds. Similarly, adding (3.1e) to (3.1g) together, we have Define the net reproductive numbers for humans and mosquitoes respectively by Besides the zero equilibrium, system (2.1)-(2.7) admits two disease-free equilibrium points Remark 1. The equilibrium point E * 01 corresponds to mosquito eradication as well as the infectiousfree state for humans. The equilibrium point E * 02 biologically corresponds to the coexistence of mosquitoes and humans, without the infection of dengue virus, which more closely fits the actual situation. In the following discussion, we always assume that R H 0 > 1. To find the disease-endemic equilibrium point of system (2.1)-(2.7), denoted by we need to solve the algebraic equations (3.1a)-(3.1g). If we define the basic reproduction number by then we obtain the existence of the disease-endemic equilibrium as follows.
From (3.1b) and (3.1c), we respectively have Hence we get the relation between E * H and I * H as i.e., Recall that at any steady state, we have Combining (3.1a) and (3.1c), one has Hence, we get the relation between E * M and I * M which reads as Again, notice that at any steady state, we have Then from (3.1e) and (3.1g), we arrive at From (3.5) to (3.9), to get the explicit expression of the disease-endemic equilibrium point, we only need to solve for I * H and I * M . To this end, combining (3.1c) and (3.7), we have Similarly, from (3.1g) and (3.9), we have Tedious but direct computations from (3.10)-(3.11) offers a linear relation between I * M and I * H as Plugging (3.12) into (3.10) and (3.11) yields the unique I * H and I * M as (3.14) It is easy to see from ( (3.9) requires again N * M > 0. From (3.9) and (3.14), with N * H > 0, we have Hence, the disease-endemic equilibrium point exists if and only if i.e., R 0 > 1. With the conclusion of Theorem 2, we complete the proof.

Stability of equilibria
Because we are dealing with a system of delay differential equations, the characteristic equation has an infinite number of roots satisfying where I is the identity matrix and the matrices J, J τ H , J τ M and J τ e have entries that are the partial derivatives of the right sides of (2.1)-(2.7) with respect to, respectively, The characteristic matrix at E * 01 is (4.5) We firstly notice that when R H 0 > 1, µ H − r H < 0. To prove the conclusion, we only need to prove that all roots of − µ M − λ + e −(λ+µ M )τ e r M 1 + θ = 0 (4.6) have negative real parts if and only if R M 0 < 1. Recall Theorem 4.7 in Smith [28] which states that λ = a + be −λτ has no roots with non-negative real parts if a + b < 0 and b ≥ a, irrelevant of the value of τ > 0. By this theorem, we see that when i.e., R M 0 < 1, all roots of (4.6) have negative real parts. Next we prove that the condition (4.7) is also sufficient to exclude the possibility for (4.6) to have roots with non-negative real parts. To proceed, assume that λ = z 1 + iz 2 with z 1 ≥ 0 is a solution of (4.6). Then r M e −µ M τ e 1 + θ e −τ e z 1 cos(τ e z 2 ) = µ M + z 1 , (4.8) It is easy to see that if (z 1 , z 2 ) satisfies (4.8)-(4.9), so does (z 1 , −z 2 ). Without loss of generality, we assume that z 2 ≥ 0. Taking square of (4.8) and (4.9) and adding them together, we have (4.10) If z 2 = 0, then z 1 > 0 due to (4.7), and from (4.10), we have a contradiction to (4.7). If z 2 > 0, from (4.10), we have also a contradiction to (4.7). This completes the proof.
Remark 2. Theorem 4 implies that if the proportion of Wolbachia-infected males to wild mosquito population, θ, satisfies θ > θ * 1 , then wild mosquitoes will be eradicated eventually, which has proved the first conclusion of Theorem 1.

Sensitivity analysis of R 0
The basic reproduction number R 0 is a function of the parameters listed in Table 2. To estimate the response of R 0 to different parameters, following the procedure in [20], we define the relative sensitivity index of R 0 with respect to each parameter p in Table 2 by where p * is taken as the baseline value in Table 2 which also yields R * 0 = 1.5869. The relative sensitivity indices are shown in Figure 1 which are ranked in their absolute values. It turns out that avoiding mosquito bites through physical and chemical means or a combination of both is the most direct and effective method to reduce the transmission of dengue virus, which has comparable performance of elevating the death rate of mosquitoes, µ M . Compared with the parameter µ M , the birth rate of mosquitoes, r M , has much less contribution to R 0 . The elevation of R 0 due to the increase of µ H can be almost equally offset by the decrease of r H . Parameters β MH , β HM , the ratio of the carrying capacities of mosquitoes to humans, r MH , and the recovery rate of humans γ, contribute equally to R 0 . Among the three delay parameters, the extrinsic incubation period of mosquitoes, τ M , is the most sensitive parameter, while the intrinsic incubation period of humans, τ H , is the least sensitive one. Unexpectedly, the release ratio, θ, only plays a minor role in controlling the basic reproduction number. The most likely reason is that the baseline value of θ is set as 1, much smaller than the release ratio in field trials. Coincidentally, further computation shows that when the release ratio is increased from 5 to 6, the basic reproduction number is decreased from 1.0447 to 0.9092(< 1). This observation is in line with our claim that for mosquito eradication or virus eradication, the optimal release ratio of Wolbachia-infected male mosquitoes to wild males is about 5 to 1 [9,33] Table 2.

5.2.
Dynamics of θ * 1 and θ * 2 Theorem 1 shows that θ * 1 is the threshold value for mosquito eradication. If mosquito eradication fails, then θ * 2 is the threshold value for virus eradication. Dynamics of θ * 1 and θ * 2 are shown in Figure  2, which are in terms of the designated parameter while letting other parameters in system (2.1)-(2.7) be fixed as the baseline values in Table 2. The threshold values θ * 1 and θ * 2 increase as r M increase; see Figure 2(a), which decrease as µ M or τ e decrease; see Figure 2(b) and (c). The threshold value θ * 1 is fixed at 12.7072 when r M , µ M and τ e are fixed. In terms of b, β HM , β MH , and r MH , the threshold θ * 2 presents a quasi-logistic growth mode; see Figure 2(d), (e), and (f). The threshold θ * 2 is a linear decreasing function of γ or τ M ; see Figure 2(g) and (i). When τ H lies between 3 and 14, it only brings a negligible change of θ * 2 ; see Figure 2(h) which is consistent with the observation in Figure 1 that the relative sensitivity index of τ H is very close to 0.

Dynamics of mosquitoes and/or viruses under different release ratios
Given the baseline values in Table 2 Figure 3 shows that if we take the release ratio θ = 13 > θ * 1 , then both the virus and mosquito will eventually be eradicated. When the production of Wolbachia-infected male mosquitoes fails to meet θ * 1 , we can only concede mosquito eradication to virus eradication.  Figure 3. The release ratio greater than θ * 1 is capable of wiping both virus and mosquito. Figure 4 shows that the release ratio θ = 6 successfully clears the virus, together with the persistence of mosquito populations. However, further lower of the release ratio to 2 which is less than the threshold value θ * 2 can neither clear virus nor eradicate mosquitoes. To see this, we initiate system  When θ = 2, Figure 5 shows that the number of infectious humans oscillates in the vicinity of their steady-state I  Figure 5. The release ratio less than θ * 2 can neither clear virus nor eradicate mosquitoes.

Conclusions
Dengue fever is one of the most common mosquito-borne viral diseases. Due to the lack of commercially available vaccines and efficient clinical cures, traditional methods have been focused on vector control by heavy applications of insecticides and environmental management. However, these programs have not prevented the spread of these diseases due to the rapid development of insecticide resistance and the continual creation of ubiquitous larval breeding sites. One novel dengue control method involves the intracellular bacterium Wolbachia, whose infection in Aedes aegypti or Aedes albopictus, the major mosquito vector of dengue virus, can greatly reduce the virus replication in mosquitoes. Wolbachia can also induce cytoplasmic incompatibility (CI) in mosquitoes, which results in early embryonic death from matings between Wolbachia-infected males and females that are either uninfected or harbor a different Wolbachia strain.
CI mechanism drives Wolbachia-infected male mosquitoes as a new weapon to suppress or eradicate wild females. However, complete population eradication in large-scale field trials is usually formidable, and we can only concede mosquito eradication to virus eradication by restricting the mosquito densities below the epidemic critical threshold. Motivated by the success and the challenge in field trials, we developed a deterministic mathematical model of human and mosquito populations interfered by the circulation of a single dengue serotype in the framework of SIER model to assess the efficacy of blocking dengue virus transmission by Wolbachia. We extended the SIER model by incorporating the Wolbachia-infected male mosquito release, where the infected males are maintained at a fixed ratio to the adult female mosquito population. Furthermore, the extrinsic incubation period in mosquito (EIP), the intrinsic incubation period in human (IIP), and the maturation delay between mating and emergence of adult mosquitoes were embedded as three delays in our model.
The threshold values of the release ratio θ for mosquito or virus eradication were found by seeking the sufficient and necessary condition on the existence of the disease-endemic equilibrium. Two explicit expressions on threshold values of θ, denoted by θ * 1 and θ * 2 with θ * 1 > θ * 2 , were obtained. When θ > θ * 1 , the mosquito population will be eradicated eventually. If it fails for mosquito eradication but θ * 2 < θ < θ * 1 , virus eradication is ensured together with the persistence of susceptible mosquitoes. When θ < θ * 2 , the emergence of the disease-endemic equilibrium makes dengue virus circulation between humans and mosquitoes possible. Sensitivity analysis of the threshold values showed that avoiding mosquito bites through physical and chemical means is the most direct and effective method to reduce the transmission of dengue virus, which has comparable performance of elevating the death rate of mosquitoes. Results from numerical simulations also confirmed our previous claim that for mosquito eradication or virus eradication, the optimal release ratio of Wolbachia-infected male mosquitoes to wild males is about 5 to 1.