Stability and Hopf bifurcation of a two species malaria model with time delays

We present a mathematical model of the transmission dynamics of two species of malaria with time lags. The model is equally applicable to two strains of a malaria species. The reproduction numbers of the two species are obtained and used as threshold parameters to study the stability and bifurcations of the equilibria of the model. We find that the model has a disease free equilibrium, which is a global attractor when the reproduction number of each species is less than one. Further, we observe that the non-disease free equilibrium of the model contains stability switches and Hopf bifurcations take place when the delays exceed the critical values.


Introduction
Malaria is a mosquito-borne disease that resulted in about 438, 000 deaths globally in 2015 (WHO, 2016). It is caused by different species of the Plasmodium parasite; the most common being Plasmodium falciparum, Plasmodium vivax, Plasmodium ovale and Plasmodium malariae. The life cycle begins when an infected female Anopheles mosquito bites a human host and injects the malaria parasites in the form of sporozites into the blood stream. The sporozites travel to the human liver where they grow, multiply and re-enter the blood stream as merozoites. A proportion of the merozoites replicates asexually within red blood cells, in the process destroying them and causing sickness while others develop into gametocytes. The incubation period within the human host is about 7-90 days depending on the species. The life cycle continues when some of the gametocytes are picked up by a female mosquito during a blood meal. The parasites grow and multiply within the mosquito. After 10-21 days, depending on the species, the parasites now in the form of sporozoites are ready for transmission to another human at the next bite, thus continuing the cycle.
Models for the transmission dynamics of the malaria parasite continue to evolve (Teboh-Ewungkem, Ngwa, & Ngonghala, 2013). Malaria models that incorporate time lags to account for the incubation periods within humans or mosquitoes have been considered in the literature but do pale in comparison to models that assume instantaneous effects.
Recently, Agyingi, Ngwa, and Wiandt (2016) proposed a model for studying the dynamics of multiple species and strains of malaria. The susceptible-infectious (SI) model for the n species/strains was developed using 2(n+1) ordinary differential equations. The immediate limitation of the model is that it does not account for the incubation periods within the host and vector. While this can be remedied by simply including latent compartments, the resulting system suffers from the downfall of containing many parameters which cannot be determined experimentally. In this paper we upgrade the SI model in Agyingi et al. (2016) to incorporate time lags which serve as the incubation periods, for two species of the parasite. We remark here that the derived model can also be used to study two strains of a single species. The introduction of time lags are known to excite instabilities in dynamical models (MacDonald, 1989). We analytically determine the critical values for the onset of bifurcations and affirm them with numerical computations. Further, we deduce competitive exclusion for the model as a consequence of determining the equilibria of the system and also confirm it numerically.

Mathematical model
Maintaining all the assumptions made in the model by Agyingi et al. (2016), we add delay parameters to both host and vectors and obtain the governing equations of the model for two species as follows: where all the parameters are assumed to be positive. The variables in the above equation are defined as follows: S h is the susceptible human population; S m is the susceptible mosquito population; and I hk are the infectious human populations for k = 1, 2; I mk are the infectious mosquito populations for k = 1, 2. N h is the total human population and N m is the total mosquitoes population. We define the parameters as: θ is the birth/death rate for humans; ν is the natural birth/death rate for mosquitoes; γ k and β k are the human recovery and death rates for species k; a is the average biting rate of the mosquitoes; ρ j is the probability that a mosquito bite will lead to an infection in a human; and σ j is the probability that mosquitoes will be infected as a result of taking blood meals from infectious humans. For k = 1, 2, τ k and τ * k are the delay terms within host and vectors, respectively. All analysis from this point forward will dwell on the infective equations only. To get rid of other variables, we assume the total population sizes for humans (N h ) and mosquitoes (N m ) are constant, normalize each of them to one and define new variables: u h = S h /N h , u m = S m /N m , x k = I hk /N h , and y k = I mk /N m , for k = 1, 2. The infective equations of the system are now given bẏ The reproduction number for species (or strain) k, denoted by R k 0 , is the average number of secondary infections that we expect to be produced by introducing a single infected human of species (or strain) k into a naive population. The introduction of an infected human will lead to the infection of α k cases in mosquitos, each of which will lead to σ k cases in humans. R k 0 is computed as the product of the quantities α k and δ k (see Agyingi et al., 2016), that is, We set the basic reproduction number, R 0 , of the model to be that of the dominant species, that is, R 0 = max{R 1 0 , R 2 0 }.

Stability analysis
In this section we compute and study the behavior of the equilibrium points of the proposed model. For purpose of simplicity and clarity in subsequent analysis, we let the delays within humans to be equal for all species and the delays in the vector to be identical for all species, that is, τ 1 = τ 2 = τ and τ * 1 = τ * 2 = τ * . Further, to eliminate the susceptible variables from (2), we let u h = 1 − x 1 − x 2 and u m = 1 − y 1 − y 2 . System (2) for k = 1, 2 becomeṡ or equivalently,ẋ Remark 1: It is obvious from system (4) that the origin is an equilibrium point of the model. This equilibrium indicates the absence of malaria, and is designated as the disease free equilibrium (DFE). In addition to the DFE, system (4) has two more equilibrium points as stated in the result below (Agyingi et al., 2016). Lemma 1: If R k 0 > 1 for k = 1, 2, then system (4) possesses two non-DFEs given by We remark here that the non-DFEs stated above are positive and therefore biologically relevant only if R k 0 = α k δ k > 1, otherwise the equilibria are negative. In order to study the behavior of the equilibria, we compute the characteristic equation of the delay system (4) as where the Jacobian matrices J 0 , J τ , and J τ * are defined, respectively as Theorem 1: The disease-free equilibrium of system (4) is locally asymptotically stable if R k 0 < 1 for k = 1, 2 and unstable if there exists an R k 0 > 1 for some k. Proof: At the DFE, the characteristic Equation (5) In case of τ = τ * = 0, we know (Agyingi et al., 2016) that the disease-free equilibrium is locally asymptotically stable when R 1 0 < 1 and R 2 0 < 1. The functions f (λ; τ , τ * ) and g(λ; τ , τ * ) are both continuous in λ, τ , and τ * , thus by the implicit function theorem their roots λ f (τ , τ * ) and λ g (τ , τ * ) are continuous functions of τ and τ * . For τ = τ * = 0 the roots are on the left side Re(z) < 0 of the complex plane (whenever R 1 0 < 1 and R 2 0 < 1), thus if for some τ > 0, τ * > 0 values any of the roots have positive real part, there exist some τ > 0 and τ The real and imaginary parts of the above equation are, respectively, given by These are the same as cos (ω(τ + τ * ))c 1 α 1 νδ 1 = −ω 2 + c 1 ν and sin ω(τ + τ * )c 1 α 1 νδ 1 = −(c 1 + ν)ω, thus, by adding the squares of the left and right sides of the equations, we obtain that This implies that no ω can satisfy this equation, as the left side is nonnegative and the right side is negative. The computation for g(λ; τ , τ * ) is analogous. In case any of the R k 0 is larger than 1, then there exist roots with positive real part and the disease-free equilibrium is unstable.
In the case when the reproduction number of each species is less than 1, we have the following stronger result. Theorem 2: The disease-free equilibrium of system (4) is globally asymptotically stable if R k 0 < 1 for k = 1, 2. Proof: Recall that (4) is derived from (2), thus it is sufficient to prove global asymptotic stability of the DFE for (2). Since S h ≤ N h and S m ≤ N m , we obtain u h ≤ 1 and u m ≤ 1 in (2); thus we have the differential inequalitẏ This is the same asẋ Multiplying (8) by e c k t , we obtain thaṫ for k = 1, 2. Integrating the above inequality on (0, t), we obtain Rearranging this equation, we have that for k = 1, 2 Now this implies that Thus we get that for k = 1, 2. The analogous computation for (9) gives that lim sup thus combining these two inequalities we obtain that for k = 1, 2 and since R k 0 < 1, we get that lim sup t x k (t) = 0; similarly, lim sup t y k (t) = 0. This concludes the proof.

Remark 2:
Recalling the non-DFEs as we see that the model immediately exhibits competitive exclusion, since only one species is represented in each of the equilibria. The analysis for the behaviors of the non-DFEs (x * 1 , x * 2 , y * 1 , y * 2 ) and (x * * 1 , x * * 2 , y * * 1 , y * * 2 ) are similar and therefore we will treat just the former.

Lemma 2: The characteristic equation of the non-DFE
and the proof follows directly.
Theorem 3: If R 1 0 < R 2 0 , then the non-DFE (x * 1 , x * 2 , y * 1 , y * 2 ) of system (4) can lose stability only if an eigenvalue of the second factor, g(λ; τ , τ * ), in the characteristic equation crosses over to the right side of the complex plane.
Similarly to earlier computations, consider the square sum of the left and right side of the appropriate equations for the cosine and sine; we obtain This gives us that thus if R 2 0 = α 2 δ 2 > α 1 δ 1 = R 1 0 , then the roots of the factor f (λ; τ , τ * ) do not cross over the imaginary axis since no ω satisfies the above equation. Therefore the loss of stability for the non-DFE (x * 1 , x * 2 , y * 1 , y * 2 ) has to come from the second factor, g(λ; τ , τ * ). To further study the stability of the non disease-free equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ), we rewrite the second factor in the product of the characteristic equation stated in Lemma 2 as To get a better insight into the structure of g(λ; τ , τ * ), we consider five different cases for the values of the delays τ and τ * . Case 1: If τ = 0 and τ * = 0 in g(λ; τ , τ * ), we obtain the no delay system and the equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) is asymptotically stable (Agyingi et al., 2016). Case 2: If τ > 0 and τ * = 0 in g(λ; τ , τ * ), then we get where F 1 = c 2 + ν + A 1 and F 2 = c 2 (ν + A 1 ). Let λ = iω be a root of the characteristic Equation (11). We obtain the following real and imaginary equations, respectively: and ωA cos ωτ − B sin ωτ = −ωF 1 , (13) where B = A(ν + A 1 ) − c 2 ν. Squaring and adding the Equations (12) and (13), we get If F 2 − B < 0, then Equation (14) has a single positive root given aŝ (12) and (13), we get the critical values of the delay τ aŝ

Substitutingω in
When τ =τ j , the characteristic Equation (11) has a pair of purely imaginary roots ±iω.
Next we investigate whether the root crosses over into the positive half of the complex plane. Differentiating the characteristic Equation (11) with respect to τ , we get At the root λ = iω, the real part of the above equation is given as It is evident from the preceding equation that if F 2 − B < 0 then Re[dλ/dτ ] −1 > 0 at τ =τ 0 . Therefore the root crosses the imaginary axis into the positive half of the complex plane. We summarize this case with the following result. Lemma 3: If τ * = 0 and the condition F 2 −B < 0 holds, then the equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) is, (i) asymptotically stable for τ ∈ [0,τ 0 ), (ii) unstable when τ >τ 0 and (iii) undergoes a Hopf bifurcation when τ =τ 0 .
Case 4: If τ = τ * > 0 in g(λ; τ , τ * ), then we get Multiplying Equation (19) through by e λτ we get where (20), then we obtain the corresponding real and imaginary parts of the resulting complex equation as where K 1 = c 2 ν − ω 2 + D 3 , K 2 = c 2 ν − ω 2 − D 3 , and K = ω(c 2 + ν). Solving for cos ωτ and sin ωτ from the Equation (21), we get Upon squaring and adding the terms in (20), we obtain The preceding equation is a polynomial of degree 8 in the variable ω: If c 2 ν + D 3 − D 2 < 0 in the above equation then h(0) < 0 and also as ω → ∞, h → ∞, thus there exists at least one positive root. For each positive root ω we can compute the value of the delay using (22), and choose the minimum as the critical delay value. Similarly to the previous cases, we can establish that the root crosses the imaginary axis into the positive half of the complex plane and therefore we obtain the following result. Lemma 5: If τ = τ * > 0 and the condition c 2 ν + D 3 − D 2 < 0 holds, then the equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) is, (i) asymptotically stable for τ ∈ [0, τ 0 ), (ii) unstable when τ > τ 0 and (iii) undergoes a Hopf bifurcation when τ = τ 0 .
Case 5: If τ > 0 and τ * > 0 in g(λ; τ , τ * ). To study this case, we keep one of the delays, for example τ * , fixed in its stable interval, assume that τ > 0, and follow the same approach illustrated in Cases 2-4. Due to the complexity of the algebra involved, we skip all the details and propose the following result. Proposition 1: For some fixed value of τ * = 0 and under appropriate conditions, there exists a critical value τ =τ > 0 such that the equilibrium (x * 1 , y * 1 , x * 2 , y * 2 ) is, (i) asymptotically stable for τ ∈ [0,τ ), (ii) unstable when τ >τ and (iii) undergoes a Hopf bifurcation when τ =τ .

Results and discussion
We present in this section numerical computations that demonstrate the analysis provided in the previous section for two species and two strains of the disease. We simulate P. falciparum and P. vivax for the species, and for the strains we consider resistant and nonresistant P. falciparum. The values of the parameters used in our simulations were taken from Agyingi et al. (2016), Pongsumpun andTang (2009), Esteva, Gumel, andde Leon (2009) and Aron and May (1982). A concise table of these values is given in Agyingi et al. (2016). In all our simulations, we assumed that there were no disease induced deaths. All simulations were made for the same initial conditions, x 1 (0) = .1 and x 2 (0) = .1, as small perturbations did not alter the behavior of the equilibria presented below.
We begin by simulating P. falciparum and P. vivax, for delay values of 7 and 10 days in humans and mosquitoes, respectively. The results are given in Figure 1. In Figure 1(a), the computed reproduction number for each species was less than one and as expected we see    that the DFE is globally asymptotically stable since both species go extinct with time. In Figure 1(b), we double the biting rate of mosquitoes while retaining every other parameter, and observed that the computed reproduction numbers were both bigger than one. This case demonstrate the asymptotic stability of a non-DFE. It also demonstrates competitive exclusion, a direct consequence of the nature of the non-DFEs as stated in Remark 2. The surviving species P. falciparum had the larger reproduction number.
In the next set of results, we investigate the behavior of the non-DFE analysed in Theorem 3-Proposition 1. We omit the case where τ = τ * = 0 as it has been discussed in Agyingi et al. (2016).
The results in In Figure 2 are for the case where τ * = 0 and the value of τ > 0 is varied. For this simulation, we investigated P. vivax and a resistant form of P. falciparum. We deduce from Figure 2(a) and (b) that the critical value of the delay τ =τ 0 that excites instability is between 28 and 29 days, that is, 28 <τ 0 < 29. Initial oscillations that appear before the critical valueτ 0 is reached, converge to the non-DFE rendering the equilibrium asymptotically stable as seen in Figure 2(a). Oscillations that start afterτ 0 only increase in amplitude making the non-DFE unstable as portrayed in Figure 2(b). Therefore the non-DFE undergoes a Hopf bifurcation at the critical delay valueτ 0 .
On investigating the case τ = 0 and the value of τ * > 0 varied, we did not see any instability. This does not contradict the result stated in Lemma 4 because the realistic incubation period of the parasite within the mosquito is at most 21 days. Thus, the maximum value of τ * = 21 was not sufficient to cause instability in the model. The next result investigates the case in which τ = τ * = 0. We see a transition from stable to unstable as indicated in Figure 3. The instability occurs when the delay value is very close to the maximum value of the parasite's incubation period in mosquitoes. From Figure 3(a) and (b), we see that the critical delay value τ 0 lies within the interval 20 < τ 0 < 21.
Finally, we consider the case summarized in Proposition 1, that is, keeping the delay τ * fixed at some non-zero value, and we vary the value of the delay τ . The results for this case are given in Figure 4. The simulations are for two strains of P. falciparum; one resistant and the other non-resistant. The value of the delay τ * was fixed at τ * = 10. On varying the value of the delay τ , we straddled a critical value τ =τ 0 on the interval 38 <τ 0 < 39 as depicted in Figure 4(a) and (b). This result confirms the existence of a Hopf bifurcation. Illustrated in Figure 5 is a phase plane portrait for the behavior of the resistant strain. We see the onset of oscillations in Figure 5(a) which eventually converge to the equilibrium point. However, in Figure 5(b) the oscillations grow in amplitude rendering the system unstable. We remark here that resistant strain survives while the non-resistant strain dies out, affirming competitive exclusion.

Conclusion
In this paper we have presented a deterministic model for malaria transmission possessing time delays that account for the incubation period of the parasite within humans and mosquitoes. The model is formulated so that it can be used to study two species or two strains of a single species of malaria. Three equilibrium points are computed for the model, one being the DFE. We show both analytically and numerically that the DFE is globally stable when each species has a reproduction number smaller than one. We investigated the non-disease free equilibria theoretically and confirmed numerically that there exist critical delay values delineating stability and instability. Further, Hopf bifurcations occur at these critical delay values. The theoretical structure of the non-disease free equilibria suggests competitive exclusion in favor of the species or strain with the higher reproduction number.

Disclosure statement
No potential conflict of interest was reported by the authors.