Discrete-time model for malaria transmission with constant releases of sterile mosquitoes

ABSTRACT In this study, we first formulate a baseline discrete-time mathematical model for malaria transmission where the survival function of mosquitoes is of Beverton–Holt type. We then introduce sterile mosquitoes to the baseline model to explore the transmission dynamics with sterile mosquitoes. We derive formulas for the reproductive number of infection and determine the existence and uniqueness of endemic fixed points as well, for the models with or without sterile mosquitoes. We then study the impact of the releases of sterile mosquitoes on the disease transmissions by investigating the effects of varying the release rates of the sterile mosquitoes. We use a numerical example to illustrate our results for all cases and finally give brief discussions of our findings.


Introduction
Mosquito-borne diseases, such as malaria, are a big concern for the public health worldwide. Malaria is a leading cause of death in many developing countries. There are 3.2 billion people, half of the world's population, living in areas at risk of malaria transmission in 106 countries and territories. In 2016, estimated 216 million cases of malaria occurred worldwide and 445,000 people died, mostly children in the African Region. About 1700 cases of malaria are diagnosed in the United States each year [13,23,30].
Malaria is transmitted indirectly not from a human to a human but through vectors as blood-feeding mosquitoes. The effective way to prevent the transmission of such diseases is to control mosquitoes and among the biological methods in control of mosquitoes is the sterile insect technique (SIT) [7], in which the natural reproductive process of the target population is disrupted. By chemical or physical methods, male mosquitoes are genetically modified to be sterile despite being sexually active. These sterile male mosquitoes are then released into the environment to mate with the wild female mosquitoes. A wild female mosquito that mates with a sterile male mosquito will either not reproduce or produce eggs that do not hatch. Repeated releases of genetically modified mosquitoes or the releases of a significantly large number of sterile mosquitoes may eventually wipe out a wild mosquito population, although it is, in practice, often more useful to consider controlling the population rather than eradicating it [5,8,29].
Mathematical models have played an important role in providing insights into the transmission dynamics of diseases and in influencing the decision making for preventing the diseases. There are many models in the literature for the study of vector-borne diseases and models incorporating sterile mosquitoes, are also formulated for the disease transmission dynamics [1,9,17,18,20,28]. However, most of the models are of continuous time, based on differential equations, to take advantage of the rich theoretical foundation of the well-developed dynamical systems theory and bifurcation theory. As time steps are sufficiently small or population sizes are sufficiently large, these models are valid. The results from those continuous-time models have shown their significance in helping understand the transmission dynamics and make useful strategies for control and prevention of the disease.
On the other hand, the timescales for the dynamics of the human and mosquito populations are significantly different. Then, discrete-time models, based on difference equations, become equally important and even more appropriate [3,10,11,14]. Moreover, discrete-time models also have many strengths in the theory of population dynamics and mathematical epidemiology [2,10,11,27,31].
In this paper, we first formulate a discrete-time susceptible-exposed-infective-recovered (SEIR) compartmental model for humans and a susceptible-exposed-infective (SEI) compartmental model for mosquitoes, based on the model in [21], as our baseline model in Section 2. We derive a formula for the reproductive number R 0 of infection and study the existence of endemic fixed points. We then incorporate the sterile mosquitoes into the baseline model such that the susceptible mosquitoes consist of wild and sterile mosquitoes in Section 3. We consider the case where the release rate of sterile mosquitoes is constant and derive a formula for the reproductive number. We investigate the impact of the releases of sterile mosquitoes on the malaria transmission based on the reproductive number and the endemic fixed points, and give a numerical example to verify our theoretical results in Section 4. Brief discussions of our findings are provided in Section 5.

Baseline model for malaria transmission
To study the impact of interactive wild and sterile mosquitoes on the malaria transmission, based on the discrete-time malaria model in [21], we formulate a new discrete-time model as our baseline model in this section, where the survivability of mosquitoes assumes the Beverton-Holt type of nonlinearity instead of the Ricker type in [21].
Since the mosquitoes lifespan is shorter than their infective period, we assume that mosquitoes are unable to recover after infection. Then we divide the mosquito population into groups of susceptible, exposed or incubating, and infective individuals, and denote their numbers, at generation n, by S v (n), E v (n) and I v (n), respectively [6,21].
For the human population, we divide it into groups of susceptible, exposed or incubating, infective and recovered individuals. We let S h (n) be the number of susceptible humans, E h (n) the number of exposed or incubating humans who are infected but not infectious yet, I h (n) the number of infective humans who are infected and also infectious, R h (n) the number of humans who are recovered from infection but partly lose their immunity [19, [24][25][26], and N h (n) = S h (n) + E h (n) + I h (n) + R h (n), the total number of humans at time n. Then, the model equations for humans and mosquitoes described in Figure 1 are given by where the parameters are given in Table 1. We assume that the infection has no effects on the birth and that the intraspecific competition only affects the survivability of the susceptible newborns. We further assume that the survivability of the susceptible newborns, , has the Beverton-Holt-II type of nonlinearity as To determine the infection rates λ v (n) and λ h (n), we let r be the average number of bites that a single mosquito makes on all human hosts. Then the total number of bites made by susceptible mosquitoes is rS v , among which, a proportion of I h /N h goes to the infected hosts, causing r · S v · (I h /N h ) · β h = β h r(I h /N h ) · S v new infections for mosquitoes, where β h is the transmission probability per bite to a susceptible mosquito. The constant input flow of the susceptible humans λ h (n): The human infection rate with 0 ≤ λ h (n) ≤ 1 α h k : k = 1, . . . , 8. Survival probabilities of susceptible, exposed or incubating, infective and recovered humans, respectively, which are all constants γ h : The developing rate of incubating humans becoming infective such that 1/γ h is the incubating period θ h : The rate of partial immunity loss, which is a constant λ v (n): The mosquito infection rate with 0 ≤ λ v (n) ≤ 1 γ v : The rate of incubating individuals becoming infective, which is a constant with 0 ≤ γ v ≤ 1 b v : The per capita birth rate of mosquitoes α v k : k = 1, . . . , 5. Survival probabilities of susceptible, incubating and infective mosquitoes, respectively, which are all constants The maximum survival probability for newborns with 0 The density-dependant factor γ v : The developing rate of incubating mosquitoes becoming infective. 1/γ v is the incubating period β h : The transmission probability per bite to a susceptible mosquito from an infective human β v : The transmission probability per bite to a susceptible human from an infective mosquito r: The number of average bites by a single mosquito on all human hosts Thus the infection rate for mosquitoes is given by where we assume β h r < 1.
Similarly, let r 1 be the average number of bites a human host receives from all mosquitoes. Then the total number of bites received by susceptible human hosts is r 1 S h , among which, a proportion of I v /N v is from infected mosquitoes, causing r 1 S h · (I v /N v ) · β v = β v r 1 (I v /N v )S h new infections for hosts, where β v is the transmission probability per bite to a susceptible human. Here, however, we need to note that the number of infected mosquitoes can be sufficiently larger than the total number of humans, that is, Then, it may lead to the human infection rate greater than 1 even when β v r < 1, and thus S h (n + 1) < 0 from (2.1), if we directly define λ h (n) in the same way as for λ v (n). Thus we consider the following factor for the infection rate for humans: Moreover, the total number of bites that all mosquitoes make is the same as the total number of bites all human hosts receive. We then need the following balance constraint: and this balance constraint and (2.4) immediately lead to Since the number of infected mosquitoes can be sufficiently larger than the total number of humans, which may lead to the factor L(n) much greater than 1, we assume that the infection rate for humans is given by where G is a positive function of L, satisfying the following conditions [10]:

The reproductive number
We derive a formula for the reproductive number by investigating the local stability of the infection-free fixed point. The Jacobian matrix of system (2.1) at the infection-free fixed point has the form where F is the fertility matrix given by T is the transition matrix given by and matrix C has the form where Then, the next generation matrix of system (2.1) [4,12,15,16] is given by where , , which leads to the eigenvalues as ρ 1 = 0 and ρ 2 = √ q 13 q 32 . According to [16], the infection-free fixed point is locally asymptotically stable if and only if the eigenvalues ρ 1,2 < 1. Thus we define the net reproductive number of infection as R 0 := ρ 2 = √ q 13 q 32 . Hence, with S v 0 and S h 0 given in (2.7). We have the following results.

The endemic fixed point
We next explore the existence of endemic fixed points of system (2.1).
The components of an endemic fixed point need to satisfy the following equations: where Then we have where . (2.14) Substituting (2.11)-(2.13) into (2.3) and (2.5), respectively, we obtain Hence, there exists an endemic fixed point if and only if there is a positive solution to equation

The interactive transmission model with sterile mosquitoes
Now suppose that sterile mosquitoes are released into a wild mosquito population and we let B(n) be the number of the sterile mosquitoes released at time n. Since sterile mosquitoes do not reproduce, there is no maturation process from larvae to adults for sterile mosquitoes. Hence, the number of sterile mosquitoes at n is just the number of released sterile mosquitoes, and the size of total mosquitoes is N v (n) + B(n). We assume sterile mosquitoes are constantly released such that B(n) := b > 0 is a positive constant. After the sterile mosquitoes are released, the mating interaction between the wild and sterile mosquitoes takes place. The number of susceptible offspring produced per mating is Then the model for the interactive mosquitoes and the human population can be described as where the human part remains the same as described in system (2.1). Assume that all of the survival probabilities of mosquitoes α v i , i = 1, 2, 3, 4, 5, in system (3.1) are the same, denoted by α v . Then the interactive dynamics between the total wild mosquitoes and the constant sterile mosquitoes are governed by the following equation: It follows from [22] that the results for (3.2) can be stated as follows.  Notice that if b > b c , for system (3.2), there exists no positive fixed point and the trivial fixed point N vb 0 = 0 is the only fixed point, which is globally asymptotically stable. Correspondingly, for system (3.1), all wild mosquitoes eventually go extinct so that there will be no malaria transmission. Thus we only consider the case of b < b c , and we assume the initial size of the wild mosquitoes is greater than N vb(−) , as shown in Figure 2.

The reproductive number and disease spread
We derive a formula for the reproductive number of infection after the sterile mosquitoes are released into the wild mosquito population with b < b c and the initial size in (N vb(−) , ∞) as follows. Consider the following infection-free fixed point of system (3.1): It follows from (2.7) that S hb 0 = S h 0 and from (3.1) that We then define function P(b) := (S vb Clearly, there exists no, one or two positive roots of (3 . Therefore, in the case of b < b c , there exist two infection-free fixed points. The Jacobian matrix J b evaluated at an infection-free fixed point has the following form: where F b is the fertility matrix F b = F, T b is the transition matrix T b = T, with F and T given in (2.8), and It follows from (3.6) that which implies that, for the two positive roots S vb that is, On the other hand, it follows from c b 11 < 1 for S vb 0 (+) that all eigenvalues of matrix C are inside the unit circle, and thus the local stability of this infection-free fixed point is determined by matrix F b + T b . It follows then from Section 2.1 that we can define the reproductive number of infection for system (3.1) as It follows from (3.6) that and we have . That is to say, the composed function R b 0 is monotone decreasing with respect to b. Hence, there is a unique thresholdb such that In fact, the thresholdb can be explicitly solved as follows.
Thus, the infection-free fixed . Thus all wild mosquitoes will be wiped out if b ≥ b c and henceb < b c . We summarize the results as follows.

Theorem 3.1: Assume sterile mosquitoes are released into the wild mosquito population constantly with the number of releases b > 0. Define the two threshold values of releases b c andb in (3.3) and (3.9), respectively. Then we have the following results.
• If b > b c , there exists no positive fixed point for the interactive mosquitoes system (3.2) and the only trivial fixed point N vb 0 = 0 of system (3.2) is globally asymptotically stable. Correspondingly, for system (3.1), all wild mosquitoes are wiped out and there will be no infection.
• If b = b c , the unique positive fixed point of system (3.2) is unstable and the trivial fixed point N vb 0 = 0 of system (3.2) is globally asymptotically stable. Correspondingly, for system (3.1), all wild mosquitoes are wiped out as well and hence there is no infection.
• Ifb < b < b c , the sterile and wild mosquitoes coexist, but the reproductive number R b 0 < 1 and the infection-free fixed point of system (3

.1) associated with the locally asymptotically stable positive fixed point N vb(+) given in (3.4) of system (3.2) is locally asymptotically stable. Thus the infection will eventually go extinct.
• If b <b such that R b 0 > 1, then the infection-free fixed point of system (3.1) is unstable. The disease spreads when the initial size of wild mosquito is greater than N vb(−) given in (3.4).

Endemic fixed point
Similarly as in Section 2.2, we determine the existence of endemic fixed points for system (3.1) as follows.
The components of an endemic fixed point for the human part are the same as those in Section 2.2 for the baseline model. The components for the wild mosquitoes at an endemic fixed point satisfy the following system: which leads to , where ω v 1 and ω v 2 are defined in (2.14). Then we obtain .
(3.12) Substituting (3.11) and (2.15) into L(t) in (2.5) yields It follows from (3.12) that . Then Thus if R b 0 > 1, there exists a positive solution λ * to H b (λ h ) = 0, that is, to λ * = G(L(λ * )). Then there exists an endemic fixed point if R b 0 > 1. Next, we prove the uniqueness of this endemic fixed point. We rewrite the equation G(L(x)) − x for convenience and take the first derivative with respect to variable x. Then we have where G (L(x)) > 0. It follows from (3.13) that and we write If k < 0, there exists a unique positive solution x * satisfying l 1 ( We then prove that L (x) is monotone decreasing for x < x * , which is equivalent to It follows from (3.13) that Hence there exists at most one positive solution x 0 satisfying l 2 (x 0 ) = 0 from Descartes' rule of sign.
In order to compare x * with x 0 , we assume x * > x 0 . Then, for any x ∈ (x 0 , x * ), L (x) is increasing and L (x) > 0 but L (x * ) = 0, which is a contradiction. Thus, x * < x 0 , and then L (x) < 0 for x < x * . It follows from (3.17) where G (L(x)) < 0, that H b (x) < 0 for x < x * . Therefore, function H b (x) is concave down for x < x * and then decreasing for x > x * with H b (0) > 0, which implies that the endemic fixed point is unique. If k > 0, then A h 1 > 0 and thus L (x) > 0 for all x. If the constant part of l 2 (x) from (3.16) is negative, l 2 (x) < 0 for all x ∈ (0, 1] and thus L (x) < 0, which implies that H b (x) < 0 from (3.17). Since function H b (x) is concave down for all x ∈ (0, 1], the endemic fixed point is unique.
If the constant part of l 2 (x) from (3.16) is positive, there exists a unique solution x 1 such that L ( where there exists a unique solutionx such that H b (x) = 0. Therefore, function H b (x) is monotone increasing for x <x and monotone decreasing for x >x. Thus the endemic fixed point is unique.
In summary, we have shown that when R b 0 > 1, the infection-free fixed point is unstable and there exists a unique endemic fixed point.
In the following, we prove that if R b 0 < 1, there exists no endemic fixed point. If k < 0, it follows from the analysis above that H b (x) is concave down if x < x * and monotone decreasing if x > x * . Thus there exists no intersection with x-axis since If k > 0, it follows from (3.16) that H b (x) is concave down with H b (0) < 0 if the constant part of l 2 (x) is negative, and hence there exists no endemic fixed point. On the other hand, if the constant part of l 2 (x) is positive, is monotone decreasing for all x ∈ (0, 1] and there exists no endemic fixed point.

Impact of releases of sterile mosquitoes
To study the impact of the releases of sterile mosquitoes on the malaria transmission dynamics, first, it is clear that if b ≥b, there is no infection. We then consider the interval (0,b), whereb is given in (3.9). For each b ∈ (0,b), the corresponding reproductive number R b 0 > 1 and there exists a unique endemic fixed point associated with λ * to equation Then it follows from (3.13) that L(b) = L(R b 0 (b), λ * (b)) since both R b 0 and λ * can be regarded as functions of variable b. Thus (4.1) and by taking the derivative of equation H b (λ(b)) = 0 with respect to b on both sides, Hence λ * (b) < 0, which implies that λ * is monotone decreasing with respect to b.
Moreover, it follows from (2.11) that , and then Thus component I h (b) of the unique endemic fixed point is monotone decreasing with respect to b. In other words, we can reduce the population of the infected humans by increasing the amount of sterile mosquitoes. Furthermore, it follows from (3.11) and (3.12) that .
Taking the derivative with respect to b, we obtain In addition, it follows from (2.15) that and N vb+ (b) < 0 from (3.8). Consequently, I v (b) < 0, which implies that component I v (b) of the unique endemic fixed point is also monotone decreasing with respect to b. As a result of increasing the releases of the sterile mosquitoes, we can possibly wipe out all of the wild mosquitoes eventually. Therefore, for b ∈ (0,b), even though R b 0 (b) > 1 such that the disease spreads and goes to a positive steady state for n → ∞, we can still reduce the components of the infected humans and mosquitoes by increasing the releases of the sterile mosquitoes to make the transmission under control.
We provide an example below to demonstrate our findings.

Example 4.1:
We use the following parameters for the transmission model:  where N vb(−) is unstable and N vb(+) is locally asymptotically stable with initial value larger than N vb (−) . N vb(+) is a decreasing function with respect to release value b.
However, the threshold value for the disease spread isb = 21.2330 such that the reproductive number R b 0 (b) < 1 if b = 30 > 21.2330 and thus the disease dies out eventually, which is shown in Figure 4 (right). If b < 21.2330, the reproductive number R b 0 (b) > 1, which makes the infection-free fixed point unstable and thus the disease spreads. Moreover, is monotone decreasing with respect to b, as shown in Figure 4 (left).
is a positive and decreasing function as shown in Figure 5 (left). Corresponding to λ h (b), there exists a unique endemic fixed point whose components I v (b) and I h (b) are monotone decreasing with respect to b as shown in Figure 5 (right), which indicates that increasing the releases of sterile mosquitoes reduces the disease spread.
To show the dynamical impact of the releases of the sterile mosquitoes, we also present the solutions of the disease transmission system in Figure 6. When no sterile mosquitoes are released, the reproductive number R 0 = 1.3039 > 1 and the disease spreads as shown in Figure 6

Concluding remarks
To have a better understanding of the impact of releasing sterile mosquitoes on malaria transmission, we first formulated a simple compartmental SEIR model for the malaria transmission, the survivability of the susceptible newborns has the Beverton-Holt-II type of nonlinearity, in (2.1) as our baseline model in Section 2. We derived a formula for the reproductive number of infection R 0 in (2.9) and showed that the infection-free fixed point of the baseline model is asymptotically stable if R 0 < 1 and unstable if R 0 > 1. We also showed that if R 0 > 1, there exists an endemic fixed point for the baseline model. We then introduced sterile mosquitoes into our baseline model and formulated the interactive model in Section 3. The human components are the same as those in (2.1) but the mosquito components are given in (3.1). We only consider the case of constant releases of sterile mosquitoes. We derived a formula for the reproductive number R b 0 , presented in (3.7), and showed R b 0 < R 0 for any release amount b. We proved that the infection-free fixed point is asymptotically stable if R b 0 < 1 and is unstable if R b 0 > 1. Using the rate of releases b as a variable, we determined threshold value b c for the existence of a unique endemic fixed point for the interactive wild and sterile mosquitoes and threshold valueb that determines whether R b 0 < 1 or R b 0 > 1, that is, whether the disease dies out or spreads. We studied the impact of the releases of sterile mosquitoes on the transmission dynamics in Section 4 by investigating the variation of the reproductive number, R b 0 (b), and the infected components, I h (b) and I v (b), induced by λ h (b), as b varies, respectively. Based on the threshold values b c andb, we provided Example 4.1 to confirm and demonstrate our findings. If the rate of releases b is greater than threshold b c , the trivial fixed point for the interactive wild and sterile mosquitoes is locally asymptotically stable. All wild mosquitoes are wiped out and thus there is no infection. On the other hand, if b < b c , while the wild mosquitoes cannot be wiped out and the two types of mosquitoes coexist, the disease can still go extinct when there are sufficient sterile mosquitoes withb < b < b c , which leads to R b 0 < 1. In the case even if we are unable to release enough sterile mosquitoes with b <b, which leads to R b 0 > 1 and thus the disease spreads when the initial wild mosquitoes is larger than N vb(−) given in (3.4), the infected components I h (b) and I v (b) are monotone decreasing with respect to b; that is, the infection can be reduced. All of the results are demonstrated by numerical simulations in Example 4.1.
In conclusion, based on the discrete-time interactive malaria model we formulated and the analysis we accomplished in this study, we showed that the releases of sterile mosquitoes have significant impacts on controlling the transmission of malaria. It will be useful and helpful in considering effective strategies for control and prevention of mosquito-borne diseases as well.