NUMERICAL SIMULATION OF MALARIA TRANSMISSION MODEL CONSIDERING SECONDARY INFECTION

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. Malaria is an infectious disease caused by Plasmodium and transmitted through the bite of female Anopheles mosquitoes. This article constructs a mathematical model to understand the spread of malaria by considering the vector–bias phenomenon in the infection process, secondary infection, and fumigation as a means of malaria control. The model is constructed as a SIRI-UV model based on six-dimensional non-linear ordinary differential equations. Analysis of the equilibrium points with their local stability and sensitivity analysis of the basic reproduction number R0 is shown analytically and numerically. Based on the analytical studies, two types of equilibrium points were obtained, namely the disease-free equilibrium points and the endemic equilibrium points. We find that the disease-free equilibrium is locally stable if R0 < 1. Our proposed model shows the possibility of a forward bifurcation, backward bifurcation, or forward bifurcation with hysteresis. To support the interpretation of the model, a numerical simulation for the sensitivity of R0 and some autonomous simulations conducted to see how the change of parameter will affect the dynamics of our model. Simulation results show that the increasing of mortality rate on mosquitoes due to fumigation will increase the probability that malaria is eliminated.


INTRODUCTION
Malaria is a mosquito-borne disease caused by protozoan parasites in the genus Plasmodium which infect humans through the bites of infected female genus Anopheles mosquitoes. According to the World Health Organization's report, there were an estimated 219 million cases of malaria in 2017, compared to 217 million cases the year before. The total funding for malaria control and elimination reached an estimated of US$ 3.1 billion in 2017 [1]. The rate of malaria incidents is highest in the following five countries: Nigeria (25%), Congo Democratic Republic (11%), Mozambique (5%), India (4%) and Uganda (4%).
Five species of Plasmodium parasites can infect humans: (1) P. falciparum, (2) P. vivax, (3) P. ovale, (4) P. malariae, and (5) P. knowlesi. Clinically, P. falciparum and P. vivax are the biggest threat to humans compared to the other species. The common symptoms include fever, headache, shivering, cold sweat, aches, nausea and vomiting. Infection of P. falciparum can lead to other complications such as severe anemia, acute respiratory distress syndrome, kidney failure, and cerebral malaria [2].
Various efforts to prevent and control the spread of malaria have been carried out to reduce the prevalence of the disease. Some of these efforts include insecticide-treated mosquito nets and indoor spraying with residual insecticides [1]. In Indonesia, the government has implemented several programs to control mosquito vectors in the form of a week of anti-mosquito nets ("pekan kelambu anti nyamuk massal") and monitoring its use, training workers against malaria, and provision of diagnostic tools and anti-malarial drugs [3].
A mathematical model for malaria transmission was first introduced in 1916 [4] and further extended by Macdonald [5,6]. Models with acquired immunity were studied later [7,8,9].
Since the resistance to malaria is not fully acquired and declines through time, without new exposures, individuals may lose immune memory and become infected again [10]. Here, we use the SIS model to describe the dynamics of malaria in the human population. For a mosquito population, we use the SI model under the assumption that mosquitoes do not recover from malaria parasites, and also, the malaria parasites do not harm the host. This type of model is based on the Ross-Macdonald model [6]. Several blood-seeking mosquitoes search for their meals by the use of host odors, breath, and sweat [11,12,13].
A vector-bias model in malaria was first introduced in [14]. This model is an extension of the Ross-Macdonald model by considering the greater attractiveness of infectious humans to mosquitoes [5,6]. Following Kingsolver's work [15] incorporates an extrinsic incubation time in mosquitoes to study the dynamics of the disease in terms of a reproduction number. Chamchod and Britton combined a vector-bias term into the malaria transmission model and showed that the greater vector-bias affects greater attractiveness of infectious humans to mosquitoes affected malaria transmission [16].
After recovery, humans with primary infection of malaria can be infected malaria for the second time as a secondary infection. Generally, with a secondary infection, humans do not show clinical symptoms (asymptomatic). However, this infection may be fatal. Therefore, we include the compartment of secondary infection in our malaria transmission model by also considering vector-bias, which is the tendency of mosquitoes to bite humans. In detail, we discretized humans in our model into four compartments: Humans with primary infection (I 1 ) represents humans who received first-time infections and show clinical symptoms; humans with secondary infection (I 2 ) describes humans with re-infection without clinical symptoms; susceptible individuals (S), and individuals who have just recovered from malaria (R) [17].
The next section of this paper discusses the construction of malaria transmission model, which is followed by mathematical analysis of the model. The analysis is the existence and local stability of DF, basic reproduction number (R 0 ), endemic equilibrium EE, and bifurcation. Section 4 explains the numerical simulations that are carried out on the model. The numerical simulations consist of sensitivity of R 0 relative to some parameters, and an autonomous simulations. Finally, the last section is the conclusion.

THE CONSTRUCTION OF MATHEMATICAL MODEL
To begin our construction, let us divide the human population into four classes based on their health status: susceptible (S); clinical symptomatic malaria individuals resulting in a susceptible individual after the first infection (I 1 ); recovered (R); asymptomatic infection resulting from recovered individuals (I 2 ); while the mosquito population is divided into two classes: susceptible (U) and infected (V ). I 2 compartment is only for non-primary infections, while the primary infection individuals are included in I 1 . In this paper, we assume that individuals who have Hence, based on transmission diagram in Fig.1, the dynamics of malaria with secondary infection is written as follows: where all parameters are positive and described as in Table 1. Denote η ≥ 1, and η is called a bias parameter. When η = 1, there is no vector-bias effect in the model.
Proof. Consider the first equation in model (1), i.e. dS dt + S( is a non decreasing function. Therefore S(t) ≥ S(0) for t ≥ 0. In the similar way, it can be proven that Next, the rate of change in the total of human population over time is By the variable separation method, Similarly, the rate of change in the mosquito population is By the variable separation method, This result can be written as Hence, N is bounded by A µ and M is bounded by B κ+θ . (1) can be stated in a biologically-feasible region as follows. Consider the region

Invariant regions. The malaria model
where Based on subsection 2.1, the region Ω in model (1) is positively-invariant with non-negative initial conditions in R 6 + .

MODEL ANALYSIS
This section discusses the analysis of the deterministic model (1), that consists of: the existence and local stability of disease-free equilibrium (DF), basic reproduction number (R 0 ) and the endemic equilibrium points (EE).
3.1. Existence and local stability of disease-free equilibrium. The malaria model (1) has a disease-free equilibrium point given by This steady state presents the extinction of malaria in both populations, human and mosquito.
One should note that the total of human population in (2) is absorbed in the susceptible population, likewise as in mosquito population. Next, we analyze the local stability of DF by linearizing the system (1) using the Jacobian matrix. The Jacobian matrix of system (1) is given by The corresponding eigenvalues are other two eigenvalues are determined from the second degree polynomial as follows: Based on the Routh-Hurwitz criterion, the following result holds: The malaria free steady state (DF) of system (1) always exists and is locally asymp- 3.2. Basic reproduction number(R 0 ). In this section, the basic reproduction number of the system (1) is discussed. In many mathematical models, the disease will always die out whenever the basic reproduction number is less than one [20,21,22,23,24]. The basic reproduction number (R 0 ) is determined from the system (1) by using the next generation matrix approach [25]. For that purpose, we define the transmission matrix (F ) and the transition matrix (V ) of system (1), respectively, as follows: By using both matrices above, the basic reproduction number of system (1) is given by: where ρ is the spectral radius operator. Please note that according to Theorem (2) , and the form of basic reproduction number in (4), we have the following corollary. Corollary 1. Malaria will be eradicated and disappear from population whenever the basic reproduction number is less than unity, and unstable otherwise. 3.3. Endemic equilibrium and bifurcation. The endemic equilibrium of system (1), that is EE, is difficult to determine explicitly. Therefore, we show this equilibrium as a function of I 1 , which is given by where , , µ + τ 1 , and b 2 = δ k 1 + k 2 τ 2 , with I * 1 is taken from the positive root of the following polynomial with coefficients a 0 , a 1 , a 2 , a 3 , a 4 can be seen in Appendix 1.
To determine the number of positive roots in (6), first we show that a 0 as a function of R 0 .
From Appendix (1), we have that Hence, it is easy to see that a 0 > 0 ↔ R 0 > 1.
Next, we can determine the number of possible positive real roots of the polynomial in equation (6) using the Descartes Rule of Signs in Table 2.
(2) If a 0 < 0 ⇔ R 0 > 1, and a 4 > 0 whereas for a 1 , a 2 and a 3 can be positive or negative then the eq.(6) has one or three positive roots.
Note that all variables in eq.(5), depend on I 1 . When the solution of (6), that is I * 1 , and I * 1 > 0, then all other variables will also be positive. So, based on the result of point 2, we have the   Table 1, except β h which we choose as the bifurcation parameter. The results can be seen in Fig. 2   FIGURE 2. Bifurcation diagram of system (1) using polynomial (6).
The relation between β h and I 1 can be distinguished based on the value of β h . Note that Fig. 2 has two fold points, namely at β h1 and β h3 . Furthermore, we have that R 0 = 1 when β h = β h2 .
From numerical calculation, we have the following results: If β h < β h1 , then system (1) has no endemic equilibrium. If β h ∈ [β 1 , β 2 ], (see Table 3  Using a similar approach with the previous description on the yellow region, the green region has an unstable DF, three endemic equilibria: one (in the curve of ω 3 ) is unstable, and the other two (small (in the curve ω 2 ) and large one (in the curve ω 4 )) are stable. Finally, the red region, which R 0 > 1, has one endemic equilibrium which is stable. The dashed line shows that the endemic equilibrium point is unstable, and the solid line shows that the endemic equilibrium point is stable. We summarized our numerical experiments in Table 3. It indicates how sensitive is β h is to the stability of the equilibrium point. Domain We close our analytical section with the following remark.
Remark 1. Our proposed malaria model shows the usual behavior regarding the relationship between the disease-free equilibrium and R 0 , i.e., it will always be stable whene R 0 < 1. An interesting feature occurs when our model shows the possibility of the existence of endemic equilibrium even though the disease-free stable. Hence, bistability phenomena appear, which makes the long time behavior of our model highly dependent on the initial condition. On the other hand, it is also possible that our model shows the existence of multiple stable endemic equilibrium, especially when R 0 is larger but close to one.

NUMERICAL SIMULATIONS
In this section, a sensitivity analysis is conducted to find how sensitive R 0 is to the robustness of parameters in system (1). Also, some numerical simulations for the autonomous system will be carried out from the mathematical model of malaria transmission in system (1) for several possible scenarios that might appear in the field.
4.1. Sensitivity of R 0 . In this section, numerical simulations are performed to investigate the sensitivity of R 0 related to controllable parameters such as recovery rate (τ 1 ) fumigation rate (θ ) and infection rate (β h ). Taking the derivative of R 0 with respect to these parameters will yield the following: .
These equations tell us that increasing θ and τ 1 will reduce R 0 . For the implementation in the field, increasing τ 1 is related to increased quality of medical treatment, drugs, health services, and other forms of intervention in a purpose to accelerate the human recovery rate. Increasing θ is not only related to large-scale fumigation but should also consider the quality of the chemical composition of insecticides. The chosen insecticide should effectively kill the adult mosquito and must not be harmful to the ecosystem. Furthermore, the possibility of resistance to insecticide must be considered carefully, since many reports show the phenomenon of mosquito resistance. On the other hand, reducing infection rate β h is related to intervention that is not only that intended to protect people from infection but also the genetic modifications of mosquitoes related to their capability to inject Plasmodium should be contemplated. Many means of intervention have massively campaigned to reduce β h such as the distribution of mosquito bed nets made from insecticide, the use of mosquito repellent, and using Wolbachia.  can be seen that increasing τ 1 and θ will significantly reduce R 0 . From our analytical result, to achieve a stable disease-free equilibrium, R 0 should be as small as possible less than one.
Therefore, implementation fumigation needs to be carried out partially with a good quality of medical treatment of malaria. Fig.3 (b) shows that carelessness in choosing a combination between fumigation and medical treatment cannot provide malaria eradication. To be specific, for a constant τ 1 , θ should be larger than its lower bound θ * (where θ * fulfill R 0 (θ * 1 ) = 1), which depend on τ 1 and is given by : On the other hand, for a specific value of θ , τ 1 should be larger from its lower bound τ * 1 (where τ * 1 fulfill R 0 (τ * 1 ) = 1) which is given by : Then, for every parameter p ∈ V , the NFSI of R 0 is defined as: Using a Definition 1, the sensitivity index of R 0 with respect to each of the parameters R 0 is as follows: Substitute all parameters in Table 1, we show the sensitivity diagram in Fig. 5. A positive index in Fig. 5 presents the proportional relation to the change of parameters with R 0 . For example, increasing β v , β h , A and β will increase R 0 . Furthermore, the value of sensitivity present how 1% change of the parameter will change R 0 . For example, increasing τ 1 for 10% will reduce R 0 for 4.99%. According to this result, we have the following remark.  Table 1 with the following initial conditions:  Table 4. From the above results, malaria does not cause an epidemic if the success rate for infection from mosquitoes to humans (β h ) is very small at 0.001, this can be seen in Fig.6 (Symptomatic Infected Human) and Table 4. Furthermore, the dynamics will continue to remain endemic when the value of β h = 0.042, i.e. I 1 = 3. Using similar parameters and initial conditions, the dynamics of mosquitoes is shown in Fig.7 with a table of equilibrium given in Table 5.  Similar interpretation given for Fig. 7 that β h must be reduced as much as possible to find a condition when all infected mosquitoes disappear. Please see Table 5 for the equilibrium point of mosquito with various values of β h in Fig. 7.

4.2.2.
Example 2 : Variation of θ . In this section, the analysis of the autonomous simulation by varying θ is discussed. In this simulation, it is shown that a proper value of the fumigation rate can eradicate the mosquito population from the population. We conduct the simulation using four values of θ , i.e. θ = 0.9 (high intensity); θ = 0.2 (medium intensity); θ = 0.08 (low intensity) and with a value of β h = 0.08. The simulation results result are given in Fig. 8 with a table of equilibrium is in Table 6. Based on our simulation and analytical result in the previous section, there is a minimum rate of fumigation that can make our system tend to the small size of endemic equilibrium.
To illustrate this, please see the green and blue curves in Fig. 8, which represents that the fumigation does not cross the minimum threshold of fumigation. The green curve show how a minimum effort of fumigation makes the infected mosquito tend to its endemic equilibrium point. On the other hand, when the fumigation rate is given in a larger value (0.2), the number of infected mosquitoes decreased significantly even without reaching the disease-free equilibrium.
Hence, since R 0 (θ = 0.2) = 0.8221 < 1 still gives a stable endemic equilibrium, we conclude that our system exhibits a backward bifurcation in this scenario. Furthermore, with significant intervention of fumigation given in the field (please see the red curve), the number of infected humans and mosquitoes decreased significantly.

CONCLUSIONS
In this paper, the transmission model for malaria by considering secondary infection has been discussed analytically and numerically. The basic reproduction number is calculated using the next-generation matrix approach. The stability of equilibrium points has been shown to be related to the basic reproduction number. We found that malaria can be eradicated if the basic reproduction number is less than unity. However, it is observed that the model might still have a stable endemic equilibrium of malaria, even though the basic reproduction number is less than unity. Hence, the policymakers should not rely solely on the threshold of the basic reproduction number being equal to one. This is because secondary infections could be fatal, leading to backward bifurcation phenomena. Our local sensitivity analysis and numerical experiments indicate that increasing the death rate of mosquitoes is the most effective way to reduce the basic reproduction number.
However, the existence of the backward bifurcation and/or forward bifurcation with hysteresis in our model suggests that malaria may still exist even though the basic reproduction number is already less than one. Further analysis needs to investigate to understand how to avoid the existence of this type of bifurcation. As another possible improvement for the model in this article, application of optimal control could be used to find the optimal strategy of malaria eradication policy. The optimal control in many epidemiological models intended to minimize number of infected individuals with minimal cost for intervention. This task is usually read as minimizing the following cost function: where x and u are vector values for infected populations, respectively, ω as the weight parameter, and T is the final time simulation. Minimizing this cost function is subject to its respective disease transmission modelẋ = f(x) with initial condition x(0) > 0. Pontriagin's Maximum principle used to characterize the optimal control problem, and can be solved using the well known forward-backward sweep method [27]. Please see [28,29,30,31] for examples of applications of optimal control problems in a ODE based disease transmission model, PDE in [32], time delay in [33], discrete model in [34], and many more.