# Dynamical Analysis on a Malaria Model with Relapse Preventive Treatment and Saturated Fumigation

Authors: Dipo Aldila
Journal: Computational and Mathematical Methods in Medicine
 (2022)

---

## Abstract

Malaria has produced health issues in many parts of the world. One of the reason is due to the recurrence phenomenon, which can happen years after the main infection has appeared in the human body. Furthermore, the fumigation intervention, which has become a major worry in several regions of the world, has yielded unsatisfactory results, as seen by the high number of cases reported each year in several African countries. We present a novel mathematical model that integrates tafenoquine treatments to prevent relapse in the human population and saturation fumigation to control mosquito populations in this study. The endemic threshold, also known as the basic reproduction number, is calculated analytically, as is the existence and local stability of the equilibrium points. Through careful investigation, we discovered that the malaria-free equilibrium is locally asymptotically stable if the basic reproduction number is less than one and unstable if it is greater than one. According to the sensitivity analysis, the utilization of tafenoquine treatment is inversely proportional to the basic reproduction number. Although our model never exhibits a backward bifurcation at the basic reproduction number equal to one, we have demonstrated that it is possible; when the basic reproduction number is greater than one, two stable malaria-endemic equilibrium can exist. As a result, when the basic reproduction number is more than one, the final state will be determined by the initial condition of the population. As a result, enormous temporal fumigation can shift the stability of our malaria model from a big endemic size to a smaller endemic size, which is more advantageous in terms of the malaria prevention strategy. Despite the fact that this is not a case study, the numerical results presented in this article are intended to support any theoretical analysis of current malaria eradication tactics in the field.


---

## Body

## 1. Introduction

Malaria is a vector-borne disease caused by the bite of a female mosquito that has been infected with Plasmodium. Of more than 100 species, only five Plasmodium species cause malaria, namely, Plasmodium vivax, Plasmodium malaria, Plasmodium falciparum, Plasmodium knowlesi, and Plasmodium ovale [1]. When this Plasmodium has entered the human bloodstream, it will attack several vital organs in the human body, especially the liver and red blood cells [2]. People who have been infected with malaria will show a variety of symptoms, including chills, fever, and headache, which can even result in death in most cases in the pediatric population.Until now, there have been many interventions launched by governments in various countries in the world to tackle the spread of malaria. These interventions include the use of vaccines, treatment, use of insecticide-treated bed nets (ITN), and vector control with fumigation and larvicides [3]. Among these mentioned interventions, vector control with fumigation is considered as the most promising and easiest policy to implement [1]. However, several problems in its implementation arise, such as the tendency of mosquitoes to become resistant to fumigants when the intervention is not well controlled [4] or the problem of limited implementation costs. In some cases in the field, the high fumigation intensity needs to be reduced when infected people are too high. This is due to the difficulty of implementation in the field when intervention costs must also be allocated to treatment for infected individuals in the hospitals.Vaccines for malaria have become one of the main concerns of governments in many parts of the world and the World Health Organization [5]. In 2021, the R21/Matrix-M vaccine has been investigated to be the second malaria vaccine, and it is stated that it has reached the minimum efficacy limit required by WHO, which is a minimum efficacy of 75% [6]. This type of vaccine has an efficacy level of 77% to reduce the chance of successful infection in humans due to an infected female Anopheles mosquito bite. In addition to vaccination, treatment interventions are also needed to prevent the severity or incidence of relapse in malaria patients. Until now, it was stated that primaquine was the primary drug used to avoid relapse in individuals infected with malaria. However, because this drug has to be taken on a regular basis (every 14 days), it has resulted in many treatments not being successful [7]. Therefore, MMV and GlaxoSmithKline (GSK) collaborated to develop a new malaria drug known as tafenoquine, which is a single dose treatment for preventing relapse in malaria-infected individuals [8].The complexity of malaria has attracted the attention of many researchers to take part in efforts to understand the mechanism of spread and the best intervention for malaria. This is due to the complexity of its infection mechanism, such as recurrence phenomena (relapse, reinfection, and recrudescence), to the problem of the most appropriate intervention. Among these researches, mathematical modeling would play an essential role. Many authors have introduced mathematical models to understand how malaria spreads among human and mosquito populations. The first mathematical model for malaria was introduced by Ross in the early 19th century [9], where he introduced the mechanism of malaria transmission involving mosquito and human populations in his model. Ross’ research was then continued by Macdonald [10] where he introduced the concept of basic reproduction number in his model. Since then, many mathematical models have been introduced by researchers to understand various important factors in the spread of malaria. Authors in [11] proposed a malaria model considering immunological memory which boost protection of reinfection phenomenon. Two host types in malaria transmission are discussed by author in [12]. Furthermore, a two-age class model for malaria transmission is discussed in [13]. A periodic biting rate of malaria mosquitoes is modelled by author in [14]. They used Floquet theory to analyze the stability of their model. Recently, author in [15] proposed a malaria model with optimal control on saturated treatment rate. Another new strategy of transmission blocking drugs for malaria is modelled by Wu and Hu in [16]. They found that increasing the transmission blocking drugs is a more pronounced effect compared to treatment intervention. Another important factors have been discussed through mathematical models such as vector-bias effect [15, 17], relapse [18, 19], reinfection [20, 21], fumigation [15, 22], temperature and seasonality [23–25], impact of Wolbachia [26], and coinfection [27]. However, the best that we know, there is no mathematical model considering the impact of potential new treatment (tafenoquine) into their model.In this paper, we introduce and investigate a new mathematical model on malaria transmission. In this model, we divide the human population into five epidemiological classes based on their health status while the mosquito population into two epidemiological classes. Several vital factors were introduced into our model: the effect of a potential new treatment for malaria to prevent relapse; vector-bias phenomena where mosquitoes are more attracted to bite the infected individuals; and fumigation intervention which depends on the number of infected individuals at timet. We perform our mathematical analysis to show the existence of a forward bifurcation and forward hysteresis phenomena on our model, which allows the possibility of existence of three different endemic equilibrium, where two of them is locally stable. Based on this phenomenon, we show from a numerical simulation that a massive fumigation intervention in a limited time window can change the dynamic direction of the system from a large endemic equilibrium to a smaller endemic point. We also show a sensitivity analysis to determine the most influential parameter to our model.This paper is organized as follows. We formulate our model in Section2. The stability of the malaria-free equilibrium point and the form of the respected basic reproduction number are shown in Section 3. In Section 4, we analyze the existence of the malaria-endemic equilibrium point. In addition, we show the possible forward hysteresis from our model in this section. Some numerical experiments on the proposed model are shown in Section 5. Finally, some relevant conclusions are given in the last section.


## 2. Mathematical Model Formulation

The proposed dynamic model for malaria transmission in this article is inspired by our previous work in [15], by taking into account two important factors. The first factor involved in our new model is the involvement of a malaria relapse prevention drug intervention (for example, tafenoquine [28]). The second factor involves fumigation intervention which is not a monotonous function. We assume that when the number of humans infected with malaria is approaching the outbreak, the intervention given can be quite large. However, when the number of infected people continues to grow, efforts for fumigation will be reduced because of the difficulty of intervention during the outbreak.This model divides the human population based on their health status and whether they received any malaria treatment or not. Hence, let human population be divided into five epidemiological classes: susceptibleS, latent E, infected I, exposed treated T, and recovered R. On the other hand, we classify mosquito population only into two classes: susceptible U and infected W. The latent individual is an individual who has already been exposed to malaria and has Plasmodium in their lever. If an individual in E gets treated with tafenoquine to prevent relapse, then they will be classified into the class of T. Only infected individual (I) can transmit the Plasmodium into the susceptible mosquito. Based on this assumption, we have the total human population which is given by
(1)Nh=S+E+T+I+R,and the total mosquito population is given by
(2)Nv=U+W.The malaria model is governed by the following system of ordinary differential equations:(3)dSdt=Λh−ΠhN,W−μhS+ξR,dEdt=ΠhN,W−u1+η+μhE,dTdt=u1E−1−pδ+pκ+μhT,dIdt=1−pδT+ηE−γ+μhI,dRdt=pκT+γI−μh+ξR,dUdt=Λv−ΠvN,U−μv+ΨI,u2U,dWdt=ΠvN,U−μv+ΨI,u2W,where ΠhN,W and ΠvN,U are the infection rate in human and mosquito population, respectively, while ΨI,u2 presents the fumigation effectiveness factors.The per capita of birth on humans and mosquitoes is denoted byΛh and Λv, respectively. The natural death rate of humans and mosquitoes is given by μh and μv, respectively. Furthermore, parameters u1 and u2 present medical treatment intervention with tafenoquine and vector control with fumigation, respectively. Letpbe the proportion of exposed individuals who get tafenoquine and succeeded in avoiding relapse afterκ−1duration of treatment. On the other hand, we assume that the 1−p proportion of individuals in T failed in treatment. Hence, we have 1−pδT as the transition from T to I due to treatment failure, where δ−1 is the incubation period of Plasmodium with the effect of tafenoquine. We denote that the recovery rate from malaria is γ, while ξ−1 is the duration of temporal immunity.We construct the force of infection in humanΠhN,W as follows. Let b be the average bite per mosquito per day. In our model, we take into account the preference of mosquito to be more attracted to bite infected human, rather than noninfected human. This phenomenon is commonly known as “vector-bias” phenomenon [29]. Based on this “vector-biased” assumption, the probability of a mosquito encountering a susceptible human is given by S/S+E+T+αI+R, where α>1 is the vector-bias parameter. Hence, total bite of all mosquito per day is given by bWS/S+E+T+αI+R. Assuming νh as the probability that the bite of infected mosquito succeeded in infecting susceptible human, then bνhWS/S+E+T+αI+R present the total of susceptible human who get infected by malaria per time. Since b and νh are constant parameters with a dimension of bite/day and 1/bite×mosquito, respectively, we assume βh≔bνh. Therefore, we have that
(4)ΠhN,W=βhWSS+E+T+αI+R.Using a similar approach, letνv be the probability of successful infection in mosquitoes; the force of infection on mosquitoes is given by
(5)ΠvN,U=βvUαIS+E+T+αI+R,whereβv≔bνv with a dimension of bite/day and 1/bite×human for b and νv, respectively.Now, we construct our fumigation termΨI,u2. We assume that the fumigation intervention depends on the number of infected individuals. Indicators of the endemic of malaria in the field can not be seen from the number of infected mosquitoes, but it can be identified by the high number of infected individuals which is reported in the media. When the number of infected individuals increases, then the intensity of fumigation will increase. However, when the number of infected individuals increases more significantly, then the effectiveness of fumigation will decrease since the policymaker may concentrate more on the number of infected individuals in the hospital, which makes them overwhelmed to control vector population in the field. Hence, we assume that ΨI,u2 should have the following properties:
(i)
When the number of the infected individual is zero, then the fumigation intervention is zero. Hence, we haveΨ0,u2=0(ii)
The fumigation intervention increases at the beginning when the number of infected individual start to increase but will decrease when the number of infected individual is sufficiently large. Hence, we have that∂ΨI,u2/∂I>0 for I∈0,Icritical and ∂ΨI,u2/∂I≤0 for I∈Icritical,∞. Note that Icritical denote the critical number of I when the policymaker is already overwhelmed to conduct an effective fumigation intervention in the fieldBased on the above assumption, we model our fumigation intervention as(6)ΨI,u2=u2Ia+I2,where a>0 is the saturated coefficient.According to the mentioned assumptions on the infection and fumigation functions, system (3) now is read as
(7)dSdt=Λh−βhWSS+E+T+αI+R−μhS+ξR,dEdt=βhWSS+E+T+αI+R−u1+η+μhE,dTdt=u1E−1−pδ+pκ+μhT,dIdt=1−pδT+ηE−γ+μhI,dRdt=pκT+γI−μh+ξR,dUdt=Λv−βvUαIS+E+T+αI+R−μv+u2Ia+I2U,dWdt=βvUαIS+E+T+αI+R−μv+u2Ia+I2W,with a nonnegative initial conditions given at timet=0. Figure 1 depicts the flow chart of our malaria transmission model. Biological interpretation and the unity of all parameters in system (7) are given in Table 1.Figure 1
Transmission diagram of malaria model in (7).Table 1
Biological interpretation of parameters in system (7).
ParDescriptionDimensionValueRef.ΛhNumber of newborn in human per dayHuman/Day1000/65×365[15]ΛvNumber of newborn in mosquitoes per dayMosquitoes/Day1000/21[15]βhInfection rate of mosquito to human1/Mosquito×day0.022[30, 31]βvInfection rate of human to mosquito1/Human×day0.24[30, 31]αVector-bias coefficient—4[32]u1Rate of treatment with tafenoquine1/Day[0,1]Assumedu2Vector control with fumigation1/Day[0,1]AssumedμhNatural death rate of human1/Day1/65×365[15]μvNatural death rate of mosquito1/Day1/21[30]ηNatural incubation rate1/day0.0833[18]pProportion of treated individual who succeeds in treatment—0.8AssumedδIncubation rate due to use of tafenoquine1/Day0.016AssumedκRecovery rate tafenoquine treatment1/Day0.0166AssumedγRecovery rate1/Day0.0035[30]aSaturation coefficient of fumigationHuman10AssumedξWaning rate of temporal immunity1/Day0.005[33]Let system (7) have an initial condition in the following set:
(8)D=S,E,T,I,R,U,W∈ℝ+7S,U>0,E,T,I,R,W≥0.To describe the feasible solution of system (7) and its biological interpretation, we have the following theorem.Theorem 1.
For initial values in (8), malaria model in system (7) has a unique solution and remains in D for all time t≥0.Proof.
Please see AppendixA for the proof.


## 3. Malaria-Free Equilibrium and the Basic Reproduction Number

The first equilibrium point of our model is the malaria-free equilibrium point. This equilibrium present a situation where all nonsusceptible population do not exist in the equilibrium condition. For this reason, letE=0, T=0, I=0, R=0, and W=0, and then, malaria-free equilibrium MFE is obtained by the following subsystem:
(9)dSdt=Λh−μhS,dUdt=Λv−μvU.Taking the right hand side of the above system, it follows that the malaria-free equilibrium of system (7) is given by
(10)MFE=S∗,E∗,T∗,I∗,R∗,U∗,W∗=Λhμh,0,0,0,0,Λvμv,0.To conduct further analysis on the qualitative behaviour of our model, it is important to determine the related basic reproduction number of our proposed model. In many epidemiological models, basic reproduction number holds an important role in determining that the diseases die out or exist in the population [34–38]. Basic reproduction number is defined as the expected number of secondary cases caused by one primary case during infection period in a completely susceptible population [39, 40]. The basic reproduction number is calculated using the next-generation matrix approach [41]. From system (7), we have that the infected compartments consist of E, T, I, and W. The Jacobian matrix of subsystem of infected compartment on system (7) evaluated in MFE can be written as F+V, where
(11)F=000βh0000000000βvΛvαμhμvΛh0,V=−u1−η−μh000u1−1−pδ−κp−μh00η1−pδ−γ−μh0000−μv,whereF and V present the transmission and transition terms. Using formula in [41], we have the next-generation matrix NGM of system (7) which is given by
(12)NGM=−ETFV−1E=0βhμvαΛvβvμhδηp+δpu1−ηκp−δη−δu1−ημhμvΛhu1+η+μhδp−κp−δ−μhγ+μh0,where ET is the transpose of E, with
(13)E=10000001.Note that each column ofF can be spanned by each column of E. Hence, the basic reproduction number of system (7) is given by
(14)R0=βvΛvαμhδη+u11−p+ημh+pκβhΛhμv2δ1−p+μh+pκu1+η+μhγ+μh.More example on the method of next-generation matrix method to determine the basic reproduction number in various epidemiological models can be seen in [42–44]. The above expression can be rewritten as a multiplication between four important component on malaria transmission on system (7) as follows.
(15)R0=C1×C2×C3×C4,whereC1=βh/u1+η+μh present the number of new latent infected human per infection period of E, C2=αβv/μv present the number of new infected mosquitoes per infection period of W, C3=Nv/Nh present the ratio of mosquitoes and human, and C4=η+u11/1+pκ+μh/1−pδ present the impact of tafenoquine intervention.According to Theorem 2 in [45], we have the following theorem regarding the local stability criteria of the malaria-free equilibrium of system (7).Theorem 2.
The malaria-free equilibrium of system (7) is locally asymptotically stable if R0<1 and unstable if R0>1.


### 3.1. Sensitivity Analysis on the Basic Reproduction Number

In many mathematical epidemiology models, understanding the impact of key parameters in determining the size of the basic reproduction number is essential to find the best optimal strategy. Therefore, we study the normalized sensitivity analysis of the basic reproduction number using the following formula [30]:
(16)ΓρR0=∂R0∂p×pR0,where ρ is any key parameter in malaria model in system (7). In our paper, we are only interested in the following parameters: βh, βv, α, u1, u2, η, δ, p, κ, γ, ξ, and a. Furthermore, we do not show the partial derivative of these parameters since it has a long expressions. Using parameter values as in Table 1, u1=0.2, u2=0, and p=0.8; the normalized sensitivity of R0 is given in Table 2 and visualized in Figure 2.Table 2
Normalized sensitivity indices ofR0 with respect to βh, βv, α, u1, u2, η, δ, p, κ, γ, ξ, and a.
Par (ρ)ΓρR0Par (ρ)ΓρR0Par (ρ)ΓρR0βh0.5βv0.5α0.5u1-0.191u20η0.191δ0.129p−0.6477κ-0.120γ-0.494ξ0a0Figure 2
Histogram of normalized sensitivity analysis ofR0.The normalized indices in Table2 are a nondimensional value, which present the percentage change of R0 for each increase value of parameter ρ for 1%. For an example, since ΓpR0=−0.6477, then increasing probability of individuals in T to succeed in treatment for 10% will reduce R0 for 6.477%. On the other hand, since ΓβhR0=0.5, then increasing βh for 10% will increase R0 for 5%. From Figure 2, we can see that βh, βv, α, η, and δ are proportional to R0. Increasing these mentioned values will increase R0. On the other hand, parameters p, γ, u1, and κ are inversely proportional to R0. Therefore, increasing the value of p, γ, u1, and κ will reduce R0. In addition, we can see that fumigation u2, rate of loss of immunity ξ, and saturated parameter a do not affect R0. Figure 2 shows the most to the less influential parameter to R0 in descending order, from left to the right.Figure3 shows the level set of R0 with respect to u1, α, and p. From Figure 3(a), we can see that increasing the value of p reduces R0. It means that more people succeed due to treatment with tafenoquine; then, the possibility to achieve malaria-free equilibrium is bigger. Same interpretation to the rate of treatment u1. We can see that more intense intervention of tafenoquine will reduce R0. In addition, we can see clearly that better quality of tafenoquine will reduce the burden of intervention in providing tafenoquine treatment to achieve malaria-free conditions. The effect of vector-bias on the success of tafenoquine intervention to reduce R0 can be seen in Figure 3(b). We can see that more bias the mosquito to be more preferring infected human will increase the R0, which makes the intervention of tafenoquine should be given more intense to reduce the value of R0.Figure 3
Level set ofR0 with respect to p, α, and u1.
(a)
R0p,u1(b)
R0α,u1


## 3.1. Sensitivity Analysis on the Basic Reproduction Number

In many mathematical epidemiology models, understanding the impact of key parameters in determining the size of the basic reproduction number is essential to find the best optimal strategy. Therefore, we study the normalized sensitivity analysis of the basic reproduction number using the following formula [30]:
(16)ΓρR0=∂R0∂p×pR0,where ρ is any key parameter in malaria model in system (7). In our paper, we are only interested in the following parameters: βh, βv, α, u1, u2, η, δ, p, κ, γ, ξ, and a. Furthermore, we do not show the partial derivative of these parameters since it has a long expressions. Using parameter values as in Table 1, u1=0.2, u2=0, and p=0.8; the normalized sensitivity of R0 is given in Table 2 and visualized in Figure 2.Table 2
Normalized sensitivity indices ofR0 with respect to βh, βv, α, u1, u2, η, δ, p, κ, γ, ξ, and a.
Par (ρ)ΓρR0Par (ρ)ΓρR0Par (ρ)ΓρR0βh0.5βv0.5α0.5u1-0.191u20η0.191δ0.129p−0.6477κ-0.120γ-0.494ξ0a0Figure 2
Histogram of normalized sensitivity analysis ofR0.The normalized indices in Table2 are a nondimensional value, which present the percentage change of R0 for each increase value of parameter ρ for 1%. For an example, since ΓpR0=−0.6477, then increasing probability of individuals in T to succeed in treatment for 10% will reduce R0 for 6.477%. On the other hand, since ΓβhR0=0.5, then increasing βh for 10% will increase R0 for 5%. From Figure 2, we can see that βh, βv, α, η, and δ are proportional to R0. Increasing these mentioned values will increase R0. On the other hand, parameters p, γ, u1, and κ are inversely proportional to R0. Therefore, increasing the value of p, γ, u1, and κ will reduce R0. In addition, we can see that fumigation u2, rate of loss of immunity ξ, and saturated parameter a do not affect R0. Figure 2 shows the most to the less influential parameter to R0 in descending order, from left to the right.Figure3 shows the level set of R0 with respect to u1, α, and p. From Figure 3(a), we can see that increasing the value of p reduces R0. It means that more people succeed due to treatment with tafenoquine; then, the possibility to achieve malaria-free equilibrium is bigger. Same interpretation to the rate of treatment u1. We can see that more intense intervention of tafenoquine will reduce R0. In addition, we can see clearly that better quality of tafenoquine will reduce the burden of intervention in providing tafenoquine treatment to achieve malaria-free conditions. The effect of vector-bias on the success of tafenoquine intervention to reduce R0 can be seen in Figure 3(b). We can see that more bias the mosquito to be more preferring infected human will increase the R0, which makes the intervention of tafenoquine should be given more intense to reduce the value of R0.Figure 3
Level set ofR0 with respect to p, α, and u1.
(a)
R0p,u1(b)
R0α,u1


## 4. The Malaria-Endemic Equilibrium

### 4.1. Existence of Malaria-Endemic Equilibrium

The malaria-endemic equilibrium of system (7) is given by
(17)MEE=S†,E†,T†,I†,R†,U†,W†,where
(18)S†=Λhμh−E†−I†−T†−R†,E†=I†δμh+γ1−p+γ+μhμh+pκδu1+η1−p+ημh+pκ,T†=I†u1γ+μhδu1+η1−p+ημh+pκ,R†=δγ+μh1−p+γκpη+μh+μhηγ+pκu1δu1+η1−p+ημh+pκξ+μh,U†=ΛvS†+E†+αI†+T†+R†a+I†2Σi=03ci,W†=Λva+I†2u2I†+μva+I†2−U†,withc0=αμvS†+E†+T†+R†, c1=aαβv+μv+u2S†+E†+T†+R†, c2=αu2+μvS†+E†+T†+R†, and c3=αβv+μv. Note that I† is taken from the positive root of the following polynomial:
(19)GΩ,I=∑j=16kjIj=0,whereΩ is the set of parameter in system (7), and
(20)k6=−μh2μvξ+μhδ1−p+μh+pκα−1u1+η+μhαβv+μvα−1γ+μh,k0=Λhμv2δ1−p+μh+pκu1+η+μhγ+μhR02−1,whilek5, k4, k3, k2, and k1 have a complex form to be written in this article. It can be seen that whenever I†>0, then E†, T†, R†, and U† are also positive. On the other hand, S† is always positive since Nh≤Λh/μh (see the proof of Theorem 1). On the other hand, since
(21)W†=Λva+I†2u2I†+μva+I†2−U†<Λvμv−U†and Nv≤Λv/μv (see the proof of Theorem 1), then we can guarantee that W† is also positive.From the expression of polynomial in (19), k6 is always negative since α>1, k0>0⇐R0>1, while another coefficient is difficult to be determined, whether it was positive or negative. Hence, using the Descartes rules of sign [46], there exists at least one positive root of polynomial (19) whenever R0>1. According to this result and the expression of MEE, we have the following result.Theorem 3.
System (7) has at least one malaria-endemic equilibrium point if R0>1.Since polynomial in (19) is a six-degree polynomial, it is possible that system (7) have more than one malaria-endemic equilibrium point. We use Descartes rules of sign [46] to analyze the maximum possibility of positive root of polynomial in (19). The result is given in Table 3 for the case when R0>1, and Table 4 for the case when R0<1.Table 3
Possible number of positive roots of polynomialGΩ,I, when R0>1.
Casek6k5k4k3k2k1k0Possible positive roots1—++++++12—++++—+1 or 33—+++—++1 or 34—+++——+1 or 35—++—+++1 or 36—++—+—+1, 3, or 57—++——++1 or 38—++———+1 or 39—+—++++1 or 310—+—++—+1, 3, or 511—+—+—++1, 3, or 512—+—+——+1, 3, or 513—+——+++1 or 314—+——+—+1, 3, or 515—+———++1 or 316—+————+1 or 317——+++++118——+++—+1 or 319——++—++1 or 320——++——+1 or 321——+—+++1 or 322——+—+—+1, 3, or 523——+——++1 or 324——+———+1 or 325———++++126———++—+1 or 327———+—++1 or 328———+——+1 or 329————+++130————+—+1 or 331—————++132——————+1Table 4
Possible number of positive roots of polynomialGΩ,I, when R0<1.
Casek6k5k4k3k2k1k0Possible positive roots33—+++++—0 or 234—++++——0 or 235—+++—+—0, 2, or 436—+++———0 or 237—++—++—0, 2, or 438—++—+——0, 2, or 439—++——+—0, 2, or 440—++————0 or 241—+—+++—0, 2, or 442—+—++——0, 2, or 443—+—+—+—0, 2, 4, or 644—+—+———0, 2, or 445—+——++—0, 2, or 446—+——+——0, 2, or 447—+———+—0, 2, or 448—+—————0 or 249——++++—0 or 250——+++——0 or 251——++—+—0, 2, or 452——++———0 or 253——+—++—0, 2, or 454——+—+——0, 2, or 455——+——+—0, 2, or 456——+————0 or 257———+++—0 or 258———++——0 or 259———+—+—0, 2, or 460———+———0 or 261————++—0 or 262————+——0 or 263—————+—0 or 264———————0From Table2, we can confirm the result in Theorem 3 that we always have at least one malaria-endemic equilibrium when R0>1. If R0>1, then we always have an odd number possibility of the positive root of polynomial (19), i.e., 1, 3, or 5 positive roots. On the other hand, malaria-endemic equilibrium is possible to vanish only when R0<1. However, we still possible to have 2, 4, or 6 positive roots of polynomial (19) when R0<1.


### 4.2. Bifurcation Analysis

In this section, we perform the bifurcation analysis of our proposed malaria model in system (7). To do this analysis, we use the well-known Castillo-Song bifurcation theorem [47] (please see [48–51] for more examples on the use of this theorem on epidemiological models). First, for numerical calculation purposes, let us redefine our proposed system (7) as follows:
(22)f1≔Λh−βhx6x1x1+x2+x3+αx4+x5−μhx1+ξx5,f2≔βhx6x1x1+x2+x3+αx4+x5−u1+η+μhx2,f3≔u1x2−1−pδ+pκ+μhx3,f4≔1−pδx3+ηx2−γ+μhx4,f5≔pκx3+γx4−μh+ξx5,f6≔Λv−βvx6αx4x1+x2+x3+αx4+x5−μv+u2x4a+x42x6,f7≔βvx6αx4x1+x2+x3+αx4+x5−μv+u2x4a+x42x7,where xi for i=1,2,⋯7 present S, E, T, I, R, U, and W, respectively, Next, we determine our bifurcation parameter to replace R0. By solving R0=1 with respect to βh, we obtain the bifurcation parameter, namely, βh=β∗, as follows:
(23)βh=β∗=p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv.The linearization ofMFE of system (22) at βh=β∗ is given by
(24)JMFE≔−μh000ξ0c170c220000c270u1c3300000η1−pδ−γ−μh00000κpγ−ξ−μh00000c640−μv0000βvΛvαμhμvΛh00−μv,with
(25)c17=−p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv,c22=−u1−η−μh,c27=p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv,c33=−1−pδ−κp−μh,c64=−βvΛvαμhμvΛh−u2Λvaμv.The Jacobian matrixJMFE has a simple zero eigenvalue, and the other three eigenvalues are explicitly negative −μh,−μh,−μh+ξ, while the other three is coming from the root of the following polynomial:
(26)Pλ=c3λ3+c2λ2+c1λ+c0=0,where
(27)c3=Λhμvu1+μv+3μh+η+γ+pκ+δ1−p,c2=δ1−pη+γ+u1+μv+2μh3μh2+μh2κp+2η+2γ+3μv+2u1+⋯+μvκp+η+γ+u1+γκp+η+u1+pκη+u1,c1=μv23μh2+μh21−pδ+2κp+η+γ+u1+δ1−pη+γ+u1+pκη+γ+u1+η+u1+μh1−pδ+pκ+μhγ+μhΛh+Λvαηβ∗βvμh,c0=Λhμv2γ+μhu1+η+μh1−pδ+pκ+μh+βh∗βvμhαΛvδ1−pη+u1+pκη+ημh.Since1−p>0, then ci for i=0,1,2,3 are positive. Since all the coefficients of Pλ are positive, then all other three eigenvalues of JMFE are negative. Therefore, we can continue using the center manifold theory to analyze the bifurcation phenomena at R0=1. Next, we use the Castillo-Chavez and Song theorem [47] to analyze the bifurcation phenomena of system (7) at R0=1.First, we calculate the right and left eigenvector ofJMFE with respect to the zero eigenvalue. The right eigenvector is given by w=w1,w2,w3,w4,w5,w6,w7T, with
(28)w1=1γ+μhξ+μhu1−μh3+p−1δ−κp−ξ−η−γ−u1μh2+ξ+η+γ+u1p−1δ−κξ+η+γ+u1p+−ξ−η−u1γ−ξu1+ημh+p−1ξ+η+u1γ+ξu1+ηδ−κξ+η+u1γ+ηξp−γξu1,w2=−δp−κp−δ−μhu1,w3=1,w4=−δηp+δpu1−ηκp−δη−δu1−ημhγ+μhu1,w5=−δ+κp+δu1−ηδ−κp−δ−μhγ+κpμhu1γ+μhξ+μhu1,w6=−aαβvμh+u2Λh−u1+ηp−1δ+ηκp+μhw3Λvu1γ+μhaΛhμv2,w7=−Λvαβvμhη+u1p−1δ−ηκp+μhΛhμv2γ+μhu1.On the other hand, the left eigenvector is given byv=v1,v2,v3,v4,v5,v6,v7 where
(29)v1=0,v2=η+u1p−1δ−ηκp+μhη+u1+μhδp−1,v3=1,v4=p−1δ−κp−μhp−1δ,v5=0,v6=0,v7=γ+μhΛhμvp−1δ−κp−μhΛvαβvμhδp−1.It is obvious thatv1=v5=v6=0. Furthermore, f3 and f4 are one degree functions. Thus, we only need to consider the second-order partial derivative of f2 and f7. By algebraic computation, we obtain the following second-order partial derivatives which have nonzero values after substituting the MFE.
(30)∂2f2∂x2∂x7=∂2f2∂x7∂x2=−βhμhΛh,∂2f2∂x3∂x7=∂2f2∂x7∂x3=−βhμhΛh,∂2f2∂x4∂x7=∂2f2∂x7∂x4=−βhμhαΛh,∂2f2∂x5∂x7=∂2f2∂x7∂x5=−βhμhΛh,∂2f7∂x1∂x4=∂2f7∂x4∂x1=−βvΛvαμh2μvΛh2,∂2f7∂x2∂x4=∂2f7∂x4∂x2=−βvΛvαμh2μvΛh2,∂2f7∂x3∂x4=∂2f7∂x4∂x3=−βvΛvαμh2μvΛh2,∂2f7∂x5∂x4=∂2f7∂x4∂x5=−βvΛvαμh2μvΛh2,∂2f7∂x4∂x6=∂2f7∂x6∂x4=βvαμhΛh,∂2f7∂x4∂x7=∂2f7∂x7∂x4=−u2a,∂2f7∂x4∂x4=−2βvΛvα2μh2μvΛh2.For the bifurcation indicators, we calculateA for system (22) which is expressed by
(31)A=v2∑i,j=17wi∂2f2∂xi∂xj+v7∑i=17wi∂2f7∂xi∂xj.We can confirm thatA is always negative (please see the expression of A in Appendix D). Meanwhile, B is given by
(32)B=v2∑i=17wi∂2f2∂xi∂βh=u1+ηp−1δ−ηκp+μh2μhβvαΛvu1+η+μhδ1−pΛhμv2γ+μhu1.Since all parameters are positive, and1−p>0, then we have that B>0. According to Castillo-Chavez and Song theorem [47], since the quantity of A is negative and B is positive, then system (22)) indicates a forward bifurcation at R0=1. We state the result in the following theorem.Theorem 4.
System (7) always exhibits a forward bifurcation at R0=1.


### 4.3. Numerical Experiments on Theorem4

In this section, we show the numerical interpretation of Theorem4. The first numerical experiment is for the bifurcation diagram of system (7), which is given in Figure 4. We use parameter values as mentioned in Table 1, except that it states differently. With this set of parameter values, we have R0=1 when βh=0.0004079. For the case of a=400, u1=0.2, and u2=0, the bifurcation diagram is shown in Figure 4(a). It can be seen that the forward bifurcation phenomenon appears, which indicates there always exists a unique endemic equilibrium point when R0>1, and no endemic equilibrium when R0<1. Furthermore, we can see that the malaria-endemic equilibrium is always stable (solid red) when R0>1. The autonomous simulation for various initial conditions is shown in Figure 5. We use Runge-Kutta adaptive step size method in MATLAB to run the autonomous simulation in this article [52] (please see [53] for further detail on the method and its algorithm). It can be seen that when R0=0.8<1, then the solution from all different initial conditions tends to the malaria-free equilibrium point (Figure 5). On the other hand, when R0>1, then all trajectories tend to the malaria-endemic equilibrium (Figure 6).Figure 4
Type of bifurcation phenomena of system (7). The red figure presents I† in MEEE, the blue curve is I∗ in MFE, and the magenta curve presents the basic reproduction number as a function of βh. The solid and dotted curve present stable and unstable equilibrium point, respectively.
(a)
Forward bifurcation phenomena of system (7) when u2=0 and a=400(b)
Forward hysteresis phenomena of system (7) when u2=0.3 and a=10Figure 5
Autonomous simulation of Figure4(a): trajectories of infected compartments for many different initial conditions toward MFE when R0=0.98<1. The left figure is simulation for the first 300 days, while the right figure is simulation for days 95000 to 100000.Figure 6
Autonomous simulation of Figure4(a): trajectories of infected compartments for many different initial conditions toward MEE when R0=4.902>1. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.The autonomous simulation of system (7) when forward hysteresis (Figure 4(b)) appears is given in Figures 7 and 8. The numerical results is using the same parameter values as in Figure 4(b). We only conduct two cases for this scenario, namely, when R0>1 but close to one in which only one stable malaria-endemic appears (Figure 7) and when two stable malaria-endemic equilibrium appears (Figure 8) when R0>1, but not too far from 1. In the first case, as shown in Figure 7, we can see that all trajectories from all different initial conditions tend to the same malaria-endemic equilibrium. However, when hysteresis starts to appear, which causes two stable malaria-endemic equilibrium, the solutions will tend to two different stable malaria-endemic equilibrium points, depending on their initial conditions. We can see that when the initial condition is close enough to the bigger malaria-endemic equilibrium (blue curve), then the solution tends to the bigger size of malaria-endemic equilibrium. The same thing happens when the initial value of infection is small enough, and then, the solution leads to the smallest stable malaria-endemic equilibrium. These simulation results indicate that fumigation may trigger the existence of multiple stable malaria-endemic equilibrium for some value when R0>1. Figure 9 confirms the statement. It can be seen that an increase in fumigation rate increases the interval when multiple stable malaria-endemic equilibrium appears.Figure 7
Autonomous simulation of Figure4(b): trajectories of infected compartments for many different initial conditions toward single MEE when R0=1.225>1. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.Figure 8
Autonomous simulation of Figure4(b): trajectories of infected compartments for many different initial conditions toward two stable MEE when R0=3.6769>1, depending on the initial conditions. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.Figure 9
The bifurcation diagram of system (7) depends on the value of fumigation rate. We use the same parameter values for each curve, except u2=0 (black), u2=0.1 (magenta), u2=0.2 (green), u2=0.3 (red), and u2=0.4 (cyan)


## 4.1. Existence of Malaria-Endemic Equilibrium

The malaria-endemic equilibrium of system (7) is given by
(17)MEE=S†,E†,T†,I†,R†,U†,W†,where
(18)S†=Λhμh−E†−I†−T†−R†,E†=I†δμh+γ1−p+γ+μhμh+pκδu1+η1−p+ημh+pκ,T†=I†u1γ+μhδu1+η1−p+ημh+pκ,R†=δγ+μh1−p+γκpη+μh+μhηγ+pκu1δu1+η1−p+ημh+pκξ+μh,U†=ΛvS†+E†+αI†+T†+R†a+I†2Σi=03ci,W†=Λva+I†2u2I†+μva+I†2−U†,withc0=αμvS†+E†+T†+R†, c1=aαβv+μv+u2S†+E†+T†+R†, c2=αu2+μvS†+E†+T†+R†, and c3=αβv+μv. Note that I† is taken from the positive root of the following polynomial:
(19)GΩ,I=∑j=16kjIj=0,whereΩ is the set of parameter in system (7), and
(20)k6=−μh2μvξ+μhδ1−p+μh+pκα−1u1+η+μhαβv+μvα−1γ+μh,k0=Λhμv2δ1−p+μh+pκu1+η+μhγ+μhR02−1,whilek5, k4, k3, k2, and k1 have a complex form to be written in this article. It can be seen that whenever I†>0, then E†, T†, R†, and U† are also positive. On the other hand, S† is always positive since Nh≤Λh/μh (see the proof of Theorem 1). On the other hand, since
(21)W†=Λva+I†2u2I†+μva+I†2−U†<Λvμv−U†and Nv≤Λv/μv (see the proof of Theorem 1), then we can guarantee that W† is also positive.From the expression of polynomial in (19), k6 is always negative since α>1, k0>0⇐R0>1, while another coefficient is difficult to be determined, whether it was positive or negative. Hence, using the Descartes rules of sign [46], there exists at least one positive root of polynomial (19) whenever R0>1. According to this result and the expression of MEE, we have the following result.Theorem 3.
System (7) has at least one malaria-endemic equilibrium point if R0>1.Since polynomial in (19) is a six-degree polynomial, it is possible that system (7) have more than one malaria-endemic equilibrium point. We use Descartes rules of sign [46] to analyze the maximum possibility of positive root of polynomial in (19). The result is given in Table 3 for the case when R0>1, and Table 4 for the case when R0<1.Table 3
Possible number of positive roots of polynomialGΩ,I, when R0>1.
Casek6k5k4k3k2k1k0Possible positive roots1—++++++12—++++—+1 or 33—+++—++1 or 34—+++——+1 or 35—++—+++1 or 36—++—+—+1, 3, or 57—++——++1 or 38—++———+1 or 39—+—++++1 or 310—+—++—+1, 3, or 511—+—+—++1, 3, or 512—+—+——+1, 3, or 513—+——+++1 or 314—+——+—+1, 3, or 515—+———++1 or 316—+————+1 or 317——+++++118——+++—+1 or 319——++—++1 or 320——++——+1 or 321——+—+++1 or 322——+—+—+1, 3, or 523——+——++1 or 324——+———+1 or 325———++++126———++—+1 or 327———+—++1 or 328———+——+1 or 329————+++130————+—+1 or 331—————++132——————+1Table 4
Possible number of positive roots of polynomialGΩ,I, when R0<1.
Casek6k5k4k3k2k1k0Possible positive roots33—+++++—0 or 234—++++——0 or 235—+++—+—0, 2, or 436—+++———0 or 237—++—++—0, 2, or 438—++—+——0, 2, or 439—++——+—0, 2, or 440—++————0 or 241—+—+++—0, 2, or 442—+—++——0, 2, or 443—+—+—+—0, 2, 4, or 644—+—+———0, 2, or 445—+——++—0, 2, or 446—+——+——0, 2, or 447—+———+—0, 2, or 448—+—————0 or 249——++++—0 or 250——+++——0 or 251——++—+—0, 2, or 452——++———0 or 253——+—++—0, 2, or 454——+—+——0, 2, or 455——+——+—0, 2, or 456——+————0 or 257———+++—0 or 258———++——0 or 259———+—+—0, 2, or 460———+———0 or 261————++—0 or 262————+——0 or 263—————+—0 or 264———————0From Table2, we can confirm the result in Theorem 3 that we always have at least one malaria-endemic equilibrium when R0>1. If R0>1, then we always have an odd number possibility of the positive root of polynomial (19), i.e., 1, 3, or 5 positive roots. On the other hand, malaria-endemic equilibrium is possible to vanish only when R0<1. However, we still possible to have 2, 4, or 6 positive roots of polynomial (19) when R0<1.


## 4.2. Bifurcation Analysis

In this section, we perform the bifurcation analysis of our proposed malaria model in system (7). To do this analysis, we use the well-known Castillo-Song bifurcation theorem [47] (please see [48–51] for more examples on the use of this theorem on epidemiological models). First, for numerical calculation purposes, let us redefine our proposed system (7) as follows:
(22)f1≔Λh−βhx6x1x1+x2+x3+αx4+x5−μhx1+ξx5,f2≔βhx6x1x1+x2+x3+αx4+x5−u1+η+μhx2,f3≔u1x2−1−pδ+pκ+μhx3,f4≔1−pδx3+ηx2−γ+μhx4,f5≔pκx3+γx4−μh+ξx5,f6≔Λv−βvx6αx4x1+x2+x3+αx4+x5−μv+u2x4a+x42x6,f7≔βvx6αx4x1+x2+x3+αx4+x5−μv+u2x4a+x42x7,where xi for i=1,2,⋯7 present S, E, T, I, R, U, and W, respectively, Next, we determine our bifurcation parameter to replace R0. By solving R0=1 with respect to βh, we obtain the bifurcation parameter, namely, βh=β∗, as follows:
(23)βh=β∗=p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv.The linearization ofMFE of system (22) at βh=β∗ is given by
(24)JMFE≔−μh000ξ0c170c220000c270u1c3300000η1−pδ−γ−μh00000κpγ−ξ−μh00000c640−μv0000βvΛvαμhμvΛh00−μv,with
(25)c17=−p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv,c22=−u1−η−μh,c27=p−1δ−κp−μhγ+μhu1+η+μhμv2Λhu1+ηp−1δ−ηκp+μhμhβvαΛv,c33=−1−pδ−κp−μh,c64=−βvΛvαμhμvΛh−u2Λvaμv.The Jacobian matrixJMFE has a simple zero eigenvalue, and the other three eigenvalues are explicitly negative −μh,−μh,−μh+ξ, while the other three is coming from the root of the following polynomial:
(26)Pλ=c3λ3+c2λ2+c1λ+c0=0,where
(27)c3=Λhμvu1+μv+3μh+η+γ+pκ+δ1−p,c2=δ1−pη+γ+u1+μv+2μh3μh2+μh2κp+2η+2γ+3μv+2u1+⋯+μvκp+η+γ+u1+γκp+η+u1+pκη+u1,c1=μv23μh2+μh21−pδ+2κp+η+γ+u1+δ1−pη+γ+u1+pκη+γ+u1+η+u1+μh1−pδ+pκ+μhγ+μhΛh+Λvαηβ∗βvμh,c0=Λhμv2γ+μhu1+η+μh1−pδ+pκ+μh+βh∗βvμhαΛvδ1−pη+u1+pκη+ημh.Since1−p>0, then ci for i=0,1,2,3 are positive. Since all the coefficients of Pλ are positive, then all other three eigenvalues of JMFE are negative. Therefore, we can continue using the center manifold theory to analyze the bifurcation phenomena at R0=1. Next, we use the Castillo-Chavez and Song theorem [47] to analyze the bifurcation phenomena of system (7) at R0=1.First, we calculate the right and left eigenvector ofJMFE with respect to the zero eigenvalue. The right eigenvector is given by w=w1,w2,w3,w4,w5,w6,w7T, with
(28)w1=1γ+μhξ+μhu1−μh3+p−1δ−κp−ξ−η−γ−u1μh2+ξ+η+γ+u1p−1δ−κξ+η+γ+u1p+−ξ−η−u1γ−ξu1+ημh+p−1ξ+η+u1γ+ξu1+ηδ−κξ+η+u1γ+ηξp−γξu1,w2=−δp−κp−δ−μhu1,w3=1,w4=−δηp+δpu1−ηκp−δη−δu1−ημhγ+μhu1,w5=−δ+κp+δu1−ηδ−κp−δ−μhγ+κpμhu1γ+μhξ+μhu1,w6=−aαβvμh+u2Λh−u1+ηp−1δ+ηκp+μhw3Λvu1γ+μhaΛhμv2,w7=−Λvαβvμhη+u1p−1δ−ηκp+μhΛhμv2γ+μhu1.On the other hand, the left eigenvector is given byv=v1,v2,v3,v4,v5,v6,v7 where
(29)v1=0,v2=η+u1p−1δ−ηκp+μhη+u1+μhδp−1,v3=1,v4=p−1δ−κp−μhp−1δ,v5=0,v6=0,v7=γ+μhΛhμvp−1δ−κp−μhΛvαβvμhδp−1.It is obvious thatv1=v5=v6=0. Furthermore, f3 and f4 are one degree functions. Thus, we only need to consider the second-order partial derivative of f2 and f7. By algebraic computation, we obtain the following second-order partial derivatives which have nonzero values after substituting the MFE.
(30)∂2f2∂x2∂x7=∂2f2∂x7∂x2=−βhμhΛh,∂2f2∂x3∂x7=∂2f2∂x7∂x3=−βhμhΛh,∂2f2∂x4∂x7=∂2f2∂x7∂x4=−βhμhαΛh,∂2f2∂x5∂x7=∂2f2∂x7∂x5=−βhμhΛh,∂2f7∂x1∂x4=∂2f7∂x4∂x1=−βvΛvαμh2μvΛh2,∂2f7∂x2∂x4=∂2f7∂x4∂x2=−βvΛvαμh2μvΛh2,∂2f7∂x3∂x4=∂2f7∂x4∂x3=−βvΛvαμh2μvΛh2,∂2f7∂x5∂x4=∂2f7∂x4∂x5=−βvΛvαμh2μvΛh2,∂2f7∂x4∂x6=∂2f7∂x6∂x4=βvαμhΛh,∂2f7∂x4∂x7=∂2f7∂x7∂x4=−u2a,∂2f7∂x4∂x4=−2βvΛvα2μh2μvΛh2.For the bifurcation indicators, we calculateA for system (22) which is expressed by
(31)A=v2∑i,j=17wi∂2f2∂xi∂xj+v7∑i=17wi∂2f7∂xi∂xj.We can confirm thatA is always negative (please see the expression of A in Appendix D). Meanwhile, B is given by
(32)B=v2∑i=17wi∂2f2∂xi∂βh=u1+ηp−1δ−ηκp+μh2μhβvαΛvu1+η+μhδ1−pΛhμv2γ+μhu1.Since all parameters are positive, and1−p>0, then we have that B>0. According to Castillo-Chavez and Song theorem [47], since the quantity of A is negative and B is positive, then system (22)) indicates a forward bifurcation at R0=1. We state the result in the following theorem.Theorem 4.
System (7) always exhibits a forward bifurcation at R0=1.


## 4.3. Numerical Experiments on Theorem4

In this section, we show the numerical interpretation of Theorem4. The first numerical experiment is for the bifurcation diagram of system (7), which is given in Figure 4. We use parameter values as mentioned in Table 1, except that it states differently. With this set of parameter values, we have R0=1 when βh=0.0004079. For the case of a=400, u1=0.2, and u2=0, the bifurcation diagram is shown in Figure 4(a). It can be seen that the forward bifurcation phenomenon appears, which indicates there always exists a unique endemic equilibrium point when R0>1, and no endemic equilibrium when R0<1. Furthermore, we can see that the malaria-endemic equilibrium is always stable (solid red) when R0>1. The autonomous simulation for various initial conditions is shown in Figure 5. We use Runge-Kutta adaptive step size method in MATLAB to run the autonomous simulation in this article [52] (please see [53] for further detail on the method and its algorithm). It can be seen that when R0=0.8<1, then the solution from all different initial conditions tends to the malaria-free equilibrium point (Figure 5). On the other hand, when R0>1, then all trajectories tend to the malaria-endemic equilibrium (Figure 6).Figure 4
Type of bifurcation phenomena of system (7). The red figure presents I† in MEEE, the blue curve is I∗ in MFE, and the magenta curve presents the basic reproduction number as a function of βh. The solid and dotted curve present stable and unstable equilibrium point, respectively.
(a)
Forward bifurcation phenomena of system (7) when u2=0 and a=400(b)
Forward hysteresis phenomena of system (7) when u2=0.3 and a=10Figure 5
Autonomous simulation of Figure4(a): trajectories of infected compartments for many different initial conditions toward MFE when R0=0.98<1. The left figure is simulation for the first 300 days, while the right figure is simulation for days 95000 to 100000.Figure 6
Autonomous simulation of Figure4(a): trajectories of infected compartments for many different initial conditions toward MEE when R0=4.902>1. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.The autonomous simulation of system (7) when forward hysteresis (Figure 4(b)) appears is given in Figures 7 and 8. The numerical results is using the same parameter values as in Figure 4(b). We only conduct two cases for this scenario, namely, when R0>1 but close to one in which only one stable malaria-endemic appears (Figure 7) and when two stable malaria-endemic equilibrium appears (Figure 8) when R0>1, but not too far from 1. In the first case, as shown in Figure 7, we can see that all trajectories from all different initial conditions tend to the same malaria-endemic equilibrium. However, when hysteresis starts to appear, which causes two stable malaria-endemic equilibrium, the solutions will tend to two different stable malaria-endemic equilibrium points, depending on their initial conditions. We can see that when the initial condition is close enough to the bigger malaria-endemic equilibrium (blue curve), then the solution tends to the bigger size of malaria-endemic equilibrium. The same thing happens when the initial value of infection is small enough, and then, the solution leads to the smallest stable malaria-endemic equilibrium. These simulation results indicate that fumigation may trigger the existence of multiple stable malaria-endemic equilibrium for some value when R0>1. Figure 9 confirms the statement. It can be seen that an increase in fumigation rate increases the interval when multiple stable malaria-endemic equilibrium appears.Figure 7
Autonomous simulation of Figure4(b): trajectories of infected compartments for many different initial conditions toward single MEE when R0=1.225>1. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.Figure 8
Autonomous simulation of Figure4(b): trajectories of infected compartments for many different initial conditions toward two stable MEE when R0=3.6769>1, depending on the initial conditions. The left figure is simulation for 300 days, while the right figure is simulation for days 95000 to 100000.Figure 9
The bifurcation diagram of system (7) depends on the value of fumigation rate. We use the same parameter values for each curve, except u2=0 (black), u2=0.1 (magenta), u2=0.2 (green), u2=0.3 (red), and u2=0.4 (cyan)


## 5. Autonomous Simulation

From the previous mathematical analysis, we found that our proposed malaria model always exhibits a forward bifurcation atR0=1. These results indicate that the basic reproduction number becomes the only endemic indicator on our proposed model. However, our model may show a multiple stable endemic equilibrium when R0>1. This phenomenon is called a forward hysteresis [54]. We found that this phenomenon was affected by the intensity of fumigation u2 and the level of population awareness a. Furthermore, our sensitivity analysis indicates how important is the intervention of tafenoquine to prevent the occurrence of relapse and fumigation to control the number of Anopheles mosquitoes in the environment. To visualize our mentioned results, we perform several numerical simulations on our autonomous simulations for several scenarios.


### 5.1. Effect of Vector-Bias

In malaria transmission, vector-bias has an important role in determining the endemic condition of the population [29]. The larger the vector-bias values, the more mosquito attracted to hunt infected humans for their meal. Figure 10 depicts the dynamic of the solution of our malaria model in (7) for several values of vector-bias parameter. We use the same parameter values as in Table 2, except u1=0.1, u2=0.05, p=0.8, and varying α from 1 to 5. With these parameters, R0 is always larger than 1, which makes the solution of system (7) tends to the malaria-endemic equilibrium. We can see that an increased value of the vector-bias parameter at the malaria-endemic equilibrium situation will increase the total population in infected humans but reduce the size of the infected mosquito population. This means that the more Anopheles mosquitoes attracted to bite infected humans than healthy humans can negatively impact the human population, where the endemic size can increase. Therefore, efforts to control the mosquito population are essential in this situation.Figure 10
Simulations showing the effect of vector-bias parameter (α) on the total of infected human (left) and mosquitoes (right). We use same parameter values as in Table 2, except u1=0.1, u2=0.05, p=0.8, and α varying: α=5 (black), α=4 (red), α=3 (green), α=2 (blue), and α=1 (cyan).
(a)(b)


### 5.2. Effect of Fumigation Saturation Parameter

The first autonomous simulation was conducted to show the impact of the fumigation saturation parametera. As we mentioned before, a smaller value a indicates a more prepared community to the increasing number of infected individuals. From the expression of R0 in (15), it can be seen that a does not appear in R0. Hence, we conclude that a does not impact the size of R0. However, as we have shown in Figure 11, a smaller value of a reduces the size of total infected humans and mosquitoes in the malaria-endemic equilibrium point. Therefore, it can be concluded that although the level of community readiness to carry out fumigation does not affect the final state of population (endemic or not), it is clear that the higher the community readiness (the smaller the value of a), then the smaller the total size of the infected population in malaria-endemic equilibrium.Figure 11
Simulations showing the effect of fumigation saturation parameter (a) on the total of infected human (left) and mosquitoes (right). We use same parameter values as in Table 2, except u1=0.1, u2=0.2, p=0.8, and a varying: a=20 (blue), a=50 (green), a=100 (red), a=200 (cyan), and a=300 (black). With this set of parameter, we have that R0=1.23.
(a)(b)


### 5.3. Effect of Different Fumigation Strategy

As we have mentioned in sensitivity analysis onR0, we find that fumigation does not affect the size of R0, but it can reduce the size of malaria-endemic equilibrium when fumigation intervention increases, as shown in Figure 12.Figure 12
Effect of fumigation on the endemic size of total infected human (left) and mosquitoes (right). Three different values ofu2 are given: 0 (red), 0.5 (blue), and 1 (green) do not change the value of R0, which is always 1.77.
(a)(b)Now, we conduct our simulation with three different scenarios, based on the measured fumigation control depending on the implementation time. In the 1st and the 2nd scenarios, we choose u2 to be changed depending on the time interval, using the following step function:
(33)u21stscenario=0.3,t≤500,0.9,500≤t≤2000,0.3,2000≤t≤10000,u22ndscenario=0.3,t≤10,0.9,10≤t≤1510,0.3,1510≤t≤10000,while the 3rd scenario when u2=0.3 for all time t∈0,10000. The result is given in Figure 13. We can see from Figure 13 that when there exist two stable malaria-endemic equilibrium points, then proper fumigation intervention may change the direction of stability of the system, which in our numerical experiment is from the large endemic size into small endemic size. When the improvement of fumigation is given too early (2nd scenario), then after the fumigation intervention loosened back into 0.3, then the dynamic of total infected human goes back to the large endemic size. On the other hand, when the intervention is given several times after the first implementation (1st scenario), then the dynamic of total infected humans is continuously going to the small endemic equilibrium. Based on this, it is necessary to consider the time for implementing an appropriate increase in fumigation intervention so that the solution dynamics can be directed to a smaller endemic point if the bistability phenomenon appears.Figure 13
Simulations on the total infected human showing the effect of early (green), late (red), and no change (blue) of fumigation intervention.


## 5.1. Effect of Vector-Bias

In malaria transmission, vector-bias has an important role in determining the endemic condition of the population [29]. The larger the vector-bias values, the more mosquito attracted to hunt infected humans for their meal. Figure 10 depicts the dynamic of the solution of our malaria model in (7) for several values of vector-bias parameter. We use the same parameter values as in Table 2, except u1=0.1, u2=0.05, p=0.8, and varying α from 1 to 5. With these parameters, R0 is always larger than 1, which makes the solution of system (7) tends to the malaria-endemic equilibrium. We can see that an increased value of the vector-bias parameter at the malaria-endemic equilibrium situation will increase the total population in infected humans but reduce the size of the infected mosquito population. This means that the more Anopheles mosquitoes attracted to bite infected humans than healthy humans can negatively impact the human population, where the endemic size can increase. Therefore, efforts to control the mosquito population are essential in this situation.Figure 10
Simulations showing the effect of vector-bias parameter (α) on the total of infected human (left) and mosquitoes (right). We use same parameter values as in Table 2, except u1=0.1, u2=0.05, p=0.8, and α varying: α=5 (black), α=4 (red), α=3 (green), α=2 (blue), and α=1 (cyan).
(a)(b)


## 5.2. Effect of Fumigation Saturation Parameter

The first autonomous simulation was conducted to show the impact of the fumigation saturation parametera. As we mentioned before, a smaller value a indicates a more prepared community to the increasing number of infected individuals. From the expression of R0 in (15), it can be seen that a does not appear in R0. Hence, we conclude that a does not impact the size of R0. However, as we have shown in Figure 11, a smaller value of a reduces the size of total infected humans and mosquitoes in the malaria-endemic equilibrium point. Therefore, it can be concluded that although the level of community readiness to carry out fumigation does not affect the final state of population (endemic or not), it is clear that the higher the community readiness (the smaller the value of a), then the smaller the total size of the infected population in malaria-endemic equilibrium.Figure 11
Simulations showing the effect of fumigation saturation parameter (a) on the total of infected human (left) and mosquitoes (right). We use same parameter values as in Table 2, except u1=0.1, u2=0.2, p=0.8, and a varying: a=20 (blue), a=50 (green), a=100 (red), a=200 (cyan), and a=300 (black). With this set of parameter, we have that R0=1.23.
(a)(b)


## 5.3. Effect of Different Fumigation Strategy

As we have mentioned in sensitivity analysis onR0, we find that fumigation does not affect the size of R0, but it can reduce the size of malaria-endemic equilibrium when fumigation intervention increases, as shown in Figure 12.Figure 12
Effect of fumigation on the endemic size of total infected human (left) and mosquitoes (right). Three different values ofu2 are given: 0 (red), 0.5 (blue), and 1 (green) do not change the value of R0, which is always 1.77.
(a)(b)Now, we conduct our simulation with three different scenarios, based on the measured fumigation control depending on the implementation time. In the 1st and the 2nd scenarios, we choose u2 to be changed depending on the time interval, using the following step function:
(33)u21stscenario=0.3,t≤500,0.9,500≤t≤2000,0.3,2000≤t≤10000,u22ndscenario=0.3,t≤10,0.9,10≤t≤1510,0.3,1510≤t≤10000,while the 3rd scenario when u2=0.3 for all time t∈0,10000. The result is given in Figure 13. We can see from Figure 13 that when there exist two stable malaria-endemic equilibrium points, then proper fumigation intervention may change the direction of stability of the system, which in our numerical experiment is from the large endemic size into small endemic size. When the improvement of fumigation is given too early (2nd scenario), then after the fumigation intervention loosened back into 0.3, then the dynamic of total infected human goes back to the large endemic size. On the other hand, when the intervention is given several times after the first implementation (1st scenario), then the dynamic of total infected humans is continuously going to the small endemic equilibrium. Based on this, it is necessary to consider the time for implementing an appropriate increase in fumigation intervention so that the solution dynamics can be directed to a smaller endemic point if the bistability phenomenon appears.Figure 13
Simulations on the total infected human showing the effect of early (green), late (red), and no change (blue) of fumigation intervention.


## 6. Conclusions

Malaria has long been a critical health problem in various parts of the world. Every year, hundreds of millions of people are at risk of becoming infected with malaria, with the majority of cases occurring in Africa. The disease is spread due to the bite of a female Anopheles mosquito and is caused by five different types of Plasmodium. Different types of Plasmodium that infect give different symptoms/serious illness that appears in patients with malaria. Various interventions have been and are being researched, such as vaccination, treatment, vector control with fumigation, and use of insecticide-treated bed net. The high number of cases in various parts of the world until now indicates that our understanding of malaria is still not sufficient to help us optimally control the spread of malaria.In this research, we introduce a new malaria model that considers two important factors: the use of a new treatment (tafenoquine) to prevent relapse and a saturated fumigation function. The fundamental properties, the existence and stability criteria of the equilibrium points, and how they relate to the basic reproduction number are analyzed in detail. We use Descartes’s rule of signs to show a possible number of malaria-endemic equilibrium points when the basic reproduction is less or larger than one. We find that it is possible to have more than one endemic equilibrium when the basic reproduction number is larger than one. Our bifurcation analysis shows how our model consistently exhibits a forward bifurcation at the basic reproduction number equal to one. However, our numerical simulations show forward bifurcation phenomena with hysteresis. This phenomenon results in the emergence of three malaria-endemic equilibrium for a basic reproduction number larger than one.Our sensitivity analysis shows that tafenoquine has a big potential to control the spread of malaria by preventing the possibility of exposed individuals from relapsing. Furthermore, we also find that although fumigation does not affect the basic reproduction number, it can reduce the number of infected individuals at malaria-endemic equilibrium. Furthermore, a numerical investigation on implementing a high intensity of fumigation in a short time intervention interval may lead to a final switching condition if the forward with hysteresis phenomena appears. We find that when fumigation is implemented in a proper time interval, the direction of endemic equilibrium can be “kicked down” into the smaller size of malaria-endemic equilibrium, which is easier to control with other intervention strategies. We hope that the results of our research in this article can provide another perspective in evaluating the possibility of implementing tafenoquine and fumigation in the field.



---

*Source: 1135452-2022-06-28.xml*