A SUPERINFECTION MODEL ON MALARIA TRANSMISSION: ANALYSIS ON THE INVASION BASIC REPRODUCTION NUMBER

unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. Abstract. Malaria is one of the many kinds of vector-borne diseases which threaten many developing countries around the world every year. Malaria is caused by more than one type of Plasmodium, which allows a superinfection between two kinds of Plasmodium in the human body. This article presents a mathematical model that describes the superinfection between Plasmodium Falciparum and Plasmodium Vivax. The model is developed as a system of a nonlinear ordinary differential equation which accommodates several essential factors, such as birth and death rate, infection process, superinfection phenomenon, recovery rate, etc. Mathematical analysis regarding the existence and stability of fixed points is discussed followed by the construction of the respective ”local” basic reproduction numbers and the ”invasion” basic reproduction numbers between Plasmodium. We found that the malaria-free equilibrium point will be locally stable if both local basic reproduction numbers are less than unity. Our results also indicate that although the ”local” basic reproduction number exhibits the existence of a single Plasmodium equilibrium, it is still possible that this equilibrium is not stable if the invasion basic reproduction number is not larger than unity. Some numerical experiments were conducted to obtain a visual interpretation of the analytical results.


INTRODUCTION
Malaria became an endemic many years ago and remains so in multiple countries all over the world. Several interventions have been conducted in many countries, such as the use of insecticides mosquito bed nets (ITNs) [1,2], use of antimalarial drugs [3,4], and fumigation strategies for larva stage [5,6] or adult mosquito stage [2].
In 2018, according to a WHO report, it was estimated that more than 400.000 death cases occurred from a total of more than 200 million human cases. Most of the cases were caused by P.
Falciparum and P. Vivax. This disease is a vector-borne disease. In this case, the malaria vector is the female Anopheles mosquito [7]. However, the details of the interaction between humans, Plasmodium, and Anopheles to effectively fight the infection are still not well understood until today.
Although many reports state that the number of cases in many countries has decreased year by year, there are still many specific phenomena in malaria that are required to be understood better, such as re-infection [8], malaria drug resistance [9], superinfection [10], and many more.
In regions of high malaria transmission, humans can be exposed to several hundred infected mosquitoes per year [11], which can potentially bring several types of Plasmodium, such as P.
Vivax and P. Falciparum. Many theories have been introduced to describe how superinfection in malaria occurs [12,13]. They stated that superinfection occurs when a single individual human hosts more than one Plasmodium species, which might occur from consecutive bites or a single bite of a mosquito with multiple Plasmodium [14]. Nevertheless, it is well established that superinfection occurs in areas that have a high endemicity of malaria [15], which suggests that superinfection originates from consecutive infectious bites.
Most of these studies construct models using a compartmental model with an ordinary differential equation approach. Mathematical analysis of the existence and stability of the equilibrium points was conducted alongside with the basic reproduction number on these references. On the other hand, several mathematical models have been established to understand how malaria spreads among the human population, which involves many aspects, such as reinfection and relapse [34,35] asymptomatic case and superinfection [36,37], optimal control approach for prevention and endemic reduction scenario [2], etc. Mathematical model analysis has been conducted by the references as mentioned earlier, such as the existence and local stability of the equilibrium points, bifurcation phenomena, and how the basic reproduction number (R 0 ) plays an important role in determining the long-term behavior of their model. Essentially, the model states that the disease would die out whenever R 0 < 1 and persist whenever R 0 > 1.
Different from previous references, in this article, we focus on the purpose of understanding how superinfection between P. Falciparum and P. Vivax occurs and how superinfection could determine the existence or persistence of malaria. In superinfection between these Plasmodium, P. Vivax can dominated P. Falciparum within the human body, leading to superinfection. Instead of "local" basic reproduction number, our research also intends to find the "invasion" basic reproduction number between these types of Plasmodium.
The paper is organized as follows: In section 2, the mathematical model construction will be explained. Some assumptions based on facts and many biological reports are stated to construct the transmission diagram and the system of the equation, which describes the malaria transmission model with superinfection phenomena. Mathematical analysis of the malaria-free and malaria-endemic equilibrium points will be analyzed in section 3. The basic reproduction number (R 0 ) will be constructed using the next-generation matrix, and the invasion basic reproduction number calculated from the local stability criteria of single-Plasmodium endemic equilibrium will also be discussed in this section. Some numerical experiments for the sensitivity analysis of the basic reproduction number and the simulation for the autonomous system will be conducted in section 4. Our work will finally be wrapped up in the final section 5.

MODEL CONSTRUCTION
In this section, we propose a mathematical model for Malaria transmission dynamics that focuses on the superinfection phenomenon between P. Falciparum and P. Vivax. To construct the model, some assumptions need to be stated as follows: (1) Population of humans and mosquitoes are closed. With this assumption, no migration is allowed into the human or mosquito population.
(2) The total human (N) and mosquito (M) populations are constant. The death rate induced by malaria will be ignored in this article.
(3) Only two strains of Plasmodium is involved; (P. Falciparum and P. Vivax).   A constant recruitment rate A h is taken for the susceptible human compartment. Each compartment in the human population has an outflow rate caused by natural death, which is denoted byμ h . The recovery rate allows infected individuals in I 1 and I 2 into the recovered compartment, with different rates γ 1 and γ 2 . The transmission of malaria is assumed to follow mass action law between susceptible humans with infected mosquitoes. Letβ 1 andβ 2 represent the disease transmission rate of P. Falciparum and P. Vivax, respectively. The transmission rateβ 2 is reduced or enhanced by the parameter δ . If δ < 1, then the superinfection is reducedβ 2 . If δ > 1, the superinfection is enhancedβ 2 .
Similarly, a constant recruitment rate of A v appears in the susceptible mosquito population. Each compartment of mosquitoes will be reduced by natural death rateμ v . Susceptible mosquitoes will get infected after biting infected human who has P. Falciparum (I 1 ) or P. Vivax (I 2 ) with the probability of successful infection rate given byη 1 andη 2 , respectively. Since the life expectancy of a mosquito is relatively short (between 3 to 4 weeks), no recovered compartment and superinfection categories have been included in the mosquito population.
Using the flow chart of the model as shown in Figure 1, finally, the system of the non-linear ordinary differential equation which describes the transmission of malaria with superinfection is given as follows : The description of parameters and the unity of system (1) are provided in Table 1. System (1) completed with a non-negative initial conditions given by :

Par. Description Unit
A h Per-capita birth rate of human  (4) with their description and unit.
Before we further analyze system (1), the following basic properties of our model need to be stated to guarantee a proper definition of our model in the biological meanings.
. Given that the initial condition X(t = 0) ≥ 0, then there exists a solution of model 1 which is non-negative for all t > 0.
Proof. From system (1), we obtain the following: The rates mentioned above are all non-negative over their boundary planes of the non-negative cone R 7 + . Therefore, we have the direction of the vector fields tended inward from their boundary. Consequently, we have, starting from the non-negative initial conditions, all solutions of system (1) positive for all time t > 0. Hence, the theorem is proved.
where M is the total mosquito population. Since we also assume that M is constant, we have With the above results, the system of equations in (1) is mathematically well-posed and epidemiologically reasonable since all variables remain non-negative for all positive t. Furthermore, since we have all solutions non-negative and a total of mosquito and human population constant, we also have each variable bounded above by A h µ h for the human population and A v µ v for mosquitoes.

MODEL ANALYSIS
Before we analyze the existence and local stability of equilibrium in (1), we will non-dimensionalize the model first by using the following: s = S N , (1) can be simplified into a non-dimensional system as follows : Since r never appeared in the other equation except in dr/dτ and u = 1 − v 1 − v 2 , we may reduce model 3 into a more simplified form as follows.
Please note that with this non-dimensionalization, we can reduce the dimension of our system from seven to five dimensions, and the number of parameters can be reduced from eleven to seven parameters. Note also that all parameters and variables in system 4 are dimensionless.
Further in this section, we will analyze system 4 rather than the original one in system 1.
3.1. Disease-free equilibrium and the basic reproduction number R 0 . The first equilibrium of system 4 is the disease-free equilibrium, which presents a situation when all infected individuals do not exist in the population. The disease-free equilibrium of system 4 is calculated by setting the infectious classes i 1 , i 2 , v 1 , v 2 to zero. This gives us the following: 0, 0, 0, 0).
With the disease-free equilibrium (Ω 1 ) in hand, we are ready to construct the basic reproduction number (R 0 ) of system 4. Here, we use the next-generation matrix in the [38] approach to construct the R 0 . The Jacobian matrix of infectious classes in system 4 is given by the following: Deconstructing J into J = F +V where F is the transmission matrix and V is the transition matrix, gives us Therefore, the next-generation matrix is given by , and the basic reproduction number is given by The threshold R 1 and R 2 is the basic reproduction number of system 4 for the P. Falciparum and P. Vivax, respectively. The local stability of Ω 1 is summarized in the following theorem.
It is evident that all eigenvalues λ will be negative if Furthermore, we find that if R 0 > 1, then at least we have one positive eigenvalue of J Ω 1 .

Endemic equilibrium (single plasmodium).
The first endemic equilibrium point where only susceptible human, infected human, and mosquito by P. Falciparum exist, which is located in the boundary s − i 1 − v 1 plane is given by the following: It is evident that Ω 2 exists uniquely in the boundary of s − i 1 − v 1 plane if R 1 > 1. Using the Jacobian matrix of system (4), the criteria for the local stability of Ω 2 can be presented by the following theorems. and is locally stable if : R 1 2 is defined as the invasion reproduction number of P. Falciparum to P. Vivax.
Proof. The characteristic equation of the Jacobian matrix of system (4) evaluated at Ω 2 can be written as follows : given bys,ī 1 andv 1 , respectively. Therefore, all roots of equation (8) have a negative real part if and only if a 3 > 0 and a 6 > 0. To guarantee a positive value of a 3 , it should be that The positive value of a 6 is a direct implication of the positiveness of a 3 . This completes the proof.
The next equilibrium is the endemic equilibrium where only susceptible humans, infected humans, and mosquito infected by P. Vivax exist, which is located in s − i 2 − v 2 plane. This equilibrium is given by the following: It is apparent that Ω 3 exists if R 2 > 1. Similarly, using the Jacobian matrix of system (4), the local stability results near Ω 3 can be presented by the following theorems.
R 1 2 is defined as the invasion reproduction number of P. Vivax to P. Falciparum.
Proof. The characteristic equation of the Jacobian matrix of system (4) evaluated at Ω 3 can be written as follows : are given by s * , i * 2 , and v * 2 , respectively. Therefore, all roots of equation (10) have a negative real part if and only if b 3 > 0 and b 6 > 0, To guarantee positive value of b 3 , it should be that while the positive value of b 6 is a direct implication of the positiveness of b 3 . This completes the proof.

Coexistence endemic equilibrium.
The last endemic equilibrium is where all compartments are positive, which is given by the following: , where v * 1 = , and K 1 > 0, K 2 > 0 has a long expression to be written in this article. The existence criteria of Ω 4 and the local stability type (numerically) is given by the following theorem.
(1) Let R 4 > 1. Then a unique coexistence equilibrium exists if and only if R 2 1 < 1 and R 1 2 < 1, that is each strain of Plasmodium can not invade the equilibrium of the other (stable).
(2) Let R 4 < 1. Then a unique coexistence equilibrium exists if and only if R 2 1 > 1 and R 1 2 > 1, that is neither strain of Plasmodium can invade the equilibrium of the other (unstable).
Proof. The existence criteria of Ω 4 is a direct consequence of the positiveness for each component in Ω 4 , by dividing R 4 into two cases, > 1 and < 1.
To summarize the results of this section, please see Table 2 which presents the potential existence and stable equilibrium region of system 4 depending on R 1 , R 2 , R 1 2 , R 2 1 . Figure 2 illustrates Table 2 with a coexistence region in the R 1 − R 2 plane.

NUMERICAL EXPERIMENTS
In this section, we will analyze our malaria model (1). The first subsection will analyze the sensitivity and elasticity of the local basic reproduction numbers (R 1 , R 2 ) and the invasion basic reproduction number (R 1 2 , R 2 1 ) followed by the numerical simulation of the original model in (1) for some different scenario which presents a possibility in the field.

Region
Long-term Behaviour Competitive Outcome R 1 > 1, R 2 > 1 i 1 and v 1 persist P. Falciparum dominates Depend on R 1 and R 2 i 2 and v 2 persist

Sensitivity and elasticity analysis of the basic reproduction number. Elasticity anal-
ysis of the basic reproduction numbers is essential to determine how they will change with respect to the change in some parameters. To assess the sensitivity of the basic reproduction number, instead of using the basic reproduction number from the non-dimensionalized model, we will transform it into the original form of the basic reproduction number and analyze its sensitivity and the elasticity.
2 and R 2 1 , we have the transformed basic reproduction numbers in to the following form (for the sake of written simplification, we ignored the "bar" symbol in each parameter): , The elasticity of R 0 with respect to parameter Γ presents the percentage change for R 0 with respect to the percentage change in the parameter Γ. The elasticity of R 0 respect to Γ is defined as follows [39] : If E Γ R 0 is positive, then R 0 increases with respect to Γ and decreases when E Γ R 0 is negative.
It is apparent from Table 3 that the natural death rate of mosquito (µ v ) is the most influential parameter in R 1 , R 2 , R 2 1 and R 1 2 . It is also evident that µ v is inversely proportional to R 1 , R 2 , R 1 2 but directly proportional to R 2 1 . These results indicate that controlling the death rate of mosquitoes is the most effective method to control the spread of malaria. Several things can be done to increase the death rate of mosquitoes, such as with genetic change in the mosquito, using fumigation to kill the adult mosquitoes, and many more. Furthermore, we find that all parameters may change the invasion reproduction number R 1 2 , where the mosquito death rate is the most influential parameter on it. However, it is interesting that not all parameters could change the value of R 2 1 , and the elasticity of µ v to this reproduction number is not as high compared to R 1 2 .

Autonomous-system simulation.
In this section, some numerical simulations will be conducted to determine how the change in parameters in the model (which presents the effort by the government to eradicate malaria from the field) affects the level of malaria endemicity.

CONCLUSIONS
A mathematical model of malaria transmission considering the superinfection phenomenon has been presented in this work. A mathematical analysis was performed to find the existence and local stability criteria of all equilibrium points. Two local basic reproduction numbers and two invasion basic reproduction numbers between Plasmoidum were found to crucial in determining the existence and local stability of equilibrium points. We found that even though the local basic reproduction number of each Plasmodium was larger than one, the existence and the stability of the equilibrium points still depended on the invasion basic reproduction number.
We also observed that the presence of superinfection in malaria dynamics might trigger the possibility of the coexistence of two Plasmodium in the field, which can increase the difficulties of malaria prevention and endemic reduction strategies.
In this article, we focused on the superinfection phenomenon. We did not consider several important factors in malaria transmission, such as the vector-bias phenomena [40], relapse and reinfection [34,35], treatment failure, and many more. Therefore, for further research, the model proposed in this article can be modified by adding the aforementioned factors. A timedependent intervention also can be used for future malaria prevention models by constructing the model as an optimal control problem.

ACKNOWLEDGMENTS
This research is financially supported by the Indonesian RistekBRIN 2021 with PDUPT research grant scheme.