Assessing the Impact of Relapse, Reinfection and Recrudescence on Malaria Eradication Policy: A Bifurcation and Optimal Control Analysis

In the present study, we propose and analyze an epidemic mathematical model for malaria dynamics, considering multiple recurrent phenomena: relapse, reinfection, and recrudescence. A limitation in hospital bed capacity, which can affect the treatment rate, is modeled using a saturated treatment function. The qualitative behavior of the model, covering the existence and stability criteria of the endemic equilibrium, is investigated rigorously. The concept of the basic reproduction number of the proposed model is obtained using the concept of the next-generation matrix. We find that the malaria-free equilibrium point is locally asymptotically stable if the basic reproduction number is less than one and unstable if it is larger than one. Our observation on the malaria-endemic equilibrium of the proposed model shows possible multiple endemic equilibria when the basic reproduction number is larger or smaller than one. Hence, we conclude that a condition of a basic reproduction number less than one is not sufficient to guarantee the extinction of malaria from the population. To test our model in a real-life situation, we fit our model parameters using the monthly incidence data from districts in Central Sumba, Indonesia called Wee Luri, which were collected from the Wee Luri Health Center. Using the first twenty months’ data from Wee Luri district, we show that our model can fit the data with a confidence interval of 95%. Both analytical and numerical experiments show that a limitation in hospital bed capacity and reinfection can trigger a more substantial possibility of the appearance of backward bifurcation. On the other hand, we find that an increase in relapse can reduce the chance of the appearance of backward bifurcation. A non-trivial result appears in that a higher probability of recrudescence (treatment failure) does not always result in the appearance of backward bifurcation. From the global sensitivity analysis using a combination of Latin hypercube sampling and partial rank correlation coefficient, we found that the initial infection rate in humans and the mosquito infection rate are the most influential parameters in determining the increase in total new human infections. We expand our model as an optimal control problem by including three types of malaria interventions, namely the use of bed net, hospitalization, and fumigation as a time-dependent variable. Using the Pontryagin maximum principle, we characterize our optimal control problem. Results from our cost-effectiveness analysis suggest that hospitalization only is the most cost-effective strategy required to control malaria disease.


Introduction
In many parts of the world, diseases spread by intermediary vectors such as mosquitoes, are still a serious problem for policymakers. Some examples of mosquito-borne diseases seasonal factors [21,22], impact of vaccines and transmission-blocking drugs [23], social hierarchy [24], etc.
Although there are many mathematical models for malaria transmission, not much research considers recurrent phenomena in their model. Niger and Gumel [25] considered reinfection in their model and found that reinfection can trigger a backward bifurcation phenomenon. Chamcod and Beier [26] considered relapses, treatment, seasonality, and G6PD in their model. Their results indicate that increasing the deficiency of G6PD will increase the number of infected humans. Ghosh et al. [27] considered two types of recurrent on their proposed model, namely relapse and reinfection. Cost-effectiveness analysis was considered in their model, to find the most effective strategies for the malaria eradication program. Woldegerima et al. [5] considered relapse in their model with transmissionblocking drug (TBD) intervention. They used incidence data from Sub-Saharan Africa to validate their model and found that TBD combined with other interventions could suppress the spread of malaria in a few years. Recently, Wang et al. [28] introduced an age-structure malaria model considering vaccination and contact prevention as eradication efforts. To date, not many mathematical models include three types of infection in malaria: recrudescence, relapse, and reinfection in a single model. Recently, [29] introduced their malaria model involving the three types of infection mentioned before. A mathematical analysis regarding the equilibrium points, the basic reproduction number, the backward bifurcation phenomena, and the global stability of the endemic equilibrium was discussed in detail. However, the authors in [29] did not yet include limited hospital capacity in their model. This factor is essential since health services will be able to be provided more optimally if the hospital capacity is also adequate in number.
Motivated by the above discussions, this paper aims to fill the gap in the mentioned references, where we consider all recurrent phenomena (relapse, reinfection, and recrudescence) together in one model. Limitation of hospital capacity to treat infected individuals is accommodated into our model as a saturated hospitalization function. We analyze the existence and stability criteria of our equilibrium points thoroughly. Bifurcation analysis on our model shows three possible types of bifurcation phenomena at basic reproduction numbers equal to one, namely forward bifurcation, backward bifurcation, and forward bifurcation with hysteresis. To calibrate our mathematical model, we estimated our parameter using malaria incidence data from the Wee Luri Health Center in Indonesia. In addition, we consider our second model, namely the malaria intervention model, as an optimal control problem. Cost-effectiveness analyses were computed to determine the most cost-effective strategy necessary to control malaria leading to its eradication. The paper is organized as follows: Formulation of the malaria model as a system of ordinary differential equations is given in Section 2 along with its fundamental properties. In Section 3, we analyze the existence criteria of our model and how it relates to the basic reproduction number. We also conduct parameter estimation on the infection parameters in Section 3. Bifurcation analysis using Castillo-Song theorem [30] is conducted analytically in Section 4 and supported with some numerical experiments to give a visualization of our results. Optimal control characterization of our intervention model is given in Section 6 and followed with optimal control simulations. We wrap up our work with discussion and conclusion in Section 8.

The Mathematical Model
To develop a mathematical model which describes the transmission process of malaria considering relapse, reinfection, and recrudescence, we divided our human population based on their health status into seven classes of the human population, namely the susceptible class of human at time t, denoted by S(t); the exposed class of human at time t, denoted by E(t); the dormant and latent classes of humans at time t, denoted by D(t) and L(t), respectively; the infected class at time t, denoted by I(t); the treated class of human at time t, denoted by T(t); and the recovered class of human, denoted by R(t). Hence, the total of the human population at time t, denoted by N(t), is given by N(t) = S(t) + E(t) + D(t) + L(t) + I(t) + T(t) + R(t). (1) In addition, the mosquito population is divided into susceptible and infected classes of mosquitoes, which is denoted by U(t) and V(t), respectively. Note that we do not consider a recovered class for the mosquito population due to the short life expectancy of female anopheles mosquitoes, which is only ±10 − 14 days [31]. Hence, the total mosquito population (denoted by M(t)) is given by We use the transmission diagram shown in Figure 1 to construct our model. The derivation of the model is given in detail as follows. The individuals of the susceptible human class are assumed to be recruited at a constant number ∆ h . We assume that all recruitment rates of humans from newborns are entering the susceptible class only. However, several reports [32,33] indicate a possibility of mother-to-children vertical transmission. We assume that no migration is involved in our model, both in the human and mosquito populations. In addition, the susceptible population may decrease due to acquiring malaria through direct contact with an infected mosquito V with a successful contact rate ofβ 1 . We consider the standard incidence rate for the interaction between humans and mosquitoes in the form ofβ 1 SV N . These individuals are then transferred to the exposed class. In the exposed class, Plasmodium is in the sporozoites form. Hence, individuals in E are not yet capable of transmitting malaria to the susceptible mosquito. We assume that the intrinsic incubation period of E is denoted byξ −1 and takes about 7-30 days [34]. In our model, it is possible that sporozoites in E may transformed into different type of Plasmodium, namely hypnozoites (transferred into D with a proportion of κ 1 ) and scizont (transferred into L with a proportion of 1 − κ 1 ).
Plasmodium in the form of sporozoites is assumed to be dormant in the liver of individuals in class D. The dormant period of Plasmodium in the human liver varies depending on the type of Plasmodium that infects the human (P. vivax or P. ovale). In our model, we assume that exposed individuals who do not experience a dormant stage are given in a proportion of κ 1 and transferred to L, while a 1 − κ 1 proportion of E that experience the dormant are stage transferred to D. As we mentioned earlier, relapse occurs due to reactivation of Plasmodium in the liver. Hence, after a dormancy period ofη 1 , individuals in D go to L due to relapse phenomena. In addition, infected individuals in the dormant class can also experience reinfection with a transmission rate ofβ 2 due to direct contact with the infected mosquito. Hence, the total number of recruitment on latent class L is given by κ 1ξ E +ηD +β 2 DV/N.
After passing through the schizont phase, Plasmodium will develop to attack the human body and enter the gametocyte phase. In this phase, humans will begin to show symptoms of malaria such as fever and flu-like illness, muscle aches, tiredness, depending on each individual's response to malaria. In our model, groups with characteristics such as this are in the infected group (I). The number of individuals in I increases due to transition from the latent group at a rate of¯ . Our model assumes that a portion of I will receive hospitalized care to get better treatment. The hospitalization action gives rise to a new problem in many countries that have a low hospital bed capacity ratio. Hospital bed capacity ratio is defined as the average number of beds available in hospitals per 1000 population in a country [35]. For example, Nigeria, as the country with the highest malaria cases in 2019 (27% cases worldwide) [2] had a hospital bed capacity ratio of only 0.5 [35]. This means that there are only five beds available in a hospital for 10,000 people in Nigeria on average. In addition, four other countries with the highest malaria cases in the world in 2020, namely the Democratic Republic of the Congo (12%), Uganda (5%), Mozambique (4%), and Niger (3%) [2], only have hospital bed capacity ratios of 0.8, 0.5, 0.7, and 0.4, respectively [35]. Therefore, when the number of cases increases significantly, the number of people who can be hospitalized will certainly not be optimal. To accommodate the limited hospitalization resources, we use the following assumption: Let the hospitalization intervention as a function depending on the number of infected individuals, f (I). We assume that increasing the number of infected individuals will decrease the hospitalization rate. Hence, we have that d f (I) dI < 0. In addition, we assume lim I→∞ f (I) = 0, which describes a situation where hospitalization is impossible to conduct when the number of infected individuals is out of control. Hence, we choose f (I) =γ 1+ωI I, whereγ andω present the hospitalization rate and half-saturation parameters, respectively. The next class is the treated class. This class increased due to hospitalization from I individuals and decreased after the treatment period, denoted byᾱ −1 . In our model, we assume that after following the treatment program, an individual in T will go to the R class due to treatment success with a proportion of 1 − κ 2 and return to I due to recrudescence with a proportion of κ 2 . As we mentioned before, recrudescence occurs when the number of Plasmodium in the treated individual can no longer be detected and then becomes detectable again. Our model assumes that this recrudescence leads to malaria treatment failure. The last class is the recovered class, which increases due to treatment success from I. It decreases due to the end of the temporary period of immunity against malaria, whose rate is denoted byφ. We assume that all classes of the human population decreased due to the natural death rate, denoted byμ h . In addition, we neglect the death rate induced by malaria.
Unlike the human population, the mosquito population is divided into two classes in our model, namely susceptible (U) and infected (V) mosquitoes. We assume that susceptible mosquitoes can get infected by malaria due to biting infected humans I and T, with rates of β 3 and β 4 , respectively. Both of these mosquito classes decreased due to the natural mosquito death rate, denoted by µ v . A recovered class in mosquitoes was not considered in our proposed model. This is because the infected female Anopheles mosquito will continue to be infected for the rest of her life [36]. Based on transmission diagram in Figure 1 and using all mentioned assumptions above, the mathematical model for malaria transmission, considering relapse, reinfection, recrudescence, and saturated hospitalization rate is given in Supplementary File (Equation (S1)). Please see Table 1 for a description of the parameters.

Parameter Estimation
In this section, we fitted the model against monthly data of malaria incidence from Wee Luri Health Center, Central Sumba, Indonesia. We use the nonlinear-squared fitting by minimizing the model solutions and the data [45]. In our model, we fit the number of treated individuals (T(t)) using the best fit parameter β 1 , β 2 , β 3 , β 4 and γ. Let Θ = (θ 1 , θ 2 , . . . , θ 5 ) present β 1 , β 2 , β 3 , β 4 , and γ, respectively, and be the set of parameters that minimize the sum of squared error between the data T t i = T t 1 , T t 2 , . . . , T t n and the model solutions f (t i , Θ) of our malaria model, which iŝ The model's solution is obtained by solving the model using the ode45 function in MATLAB. The other parameter values are given in Table 1. The fitted parameter values are given in Figure 2 and a plot of the data and model's solutions is given in Figure 3. It can be seen that the data is still in the range of the confidence intervals given in Figure 3.

Mathematical Analysis
In this chapter, analysis of the mathematical model begins by non-dimensionalizing the model on both variables and parameters. After the non-dimensionalization process, an analysis of the existence of all equilibrium points is discussed in relation to the basic reproduction number (R 0 ). In addition, we analyzed the effects of three types of recurrences on the size of R 0 .

Non-Dimensional Model
To facilitate an easier analytical study, we utilized a non-dimensionalization process of the variables and parameters in the original model (Equation (S1) in Supplementary File). By assuming a constant total population and using some new transformations (see the Supplementary Section S1 for the details), the model can be simplified into less dimensions (please see Equation (S4) in the Supplementary File).
With this non-dimensionalization process, we reduce our system from 9 to 7 variables and from 17 to 14 parameters. Since our non-dimensionalized model presents the number of humans and mosquitoes, then the value should be non-negative. Hence, it is crucial to show that the non-dimensionalized model will always give a non-negative solution for all time t ≥ 0. These properties are stated in the following theorem.
and v ≥ 0, then any solution of the non-dimensionalized model will always be non-negative for all time t ≥ 0.
Proof. Please see Supplementary Section S1 for the proof.
Since the total population is always constant, and each variable in the non-dimensionalized malaria model (see Supplementary Section S1, Equation (S4)) always non-negative, then we have that each variable is bounded at [0, 1].

Malaria-Free Equilibrium and the Basic Reproduction Number
In this section, we aim to determine the malaria-free equilibrium point of the nondimensionalized model, and the respective basic reproduction number (R 0 ). The malariafree equilibrium point is given by Next, we calculate the basic reproduction number, which is defined as the expected number of secondary malaria cases caused by one primary malaria case in a single infection period in a completely susceptible population [4]. The basic reproduction number provides a dimensionless threshold, which becomes the endemic indicator of the respected epidemiological models. In many epidemiological models [46][47][48][49], it was discovered that there is a chance that the disease will disappear when the basic reproduction number is less than one, and always exist when it is larger than one. Hence, it is important to find the R 0 , so we can determine whether malaria will exist or disappear from the population. The basic reproduction number is taken from the spectral radius of the next-generation matrix of the respective model. Using next-generation matrix approach [50][51][52], the R 0 is given by It is easy to show that our model satisfies the five conditions described in [50]. Hence, using the result in [50], we understand that whenever the basic reproduction number is less than one, we will have a chance to eliminate malaria from the environment. On the other hand, if the basic reproduction number is larger than one, then we will always have malaria in the environment. From the form of R 0 in (5), we can see that the reinfection parameter β 2 does not appear. Hence, we can conclude that reinfection phenomena do not affect the size of the basic reproduction number. In addition, we can see that R 0 is linearly proportional to all infection parameters, except for β 2 (β 1 , β 3 , β 4 ). To give further analysis on the impact of recurrence phenomena on the malaria transmission on our proposed model, we will analyze the basic reproduction number of a simple case using our proposed model.

1.
Model without recurrence. When all recurrence phenomena are not involved in the model, i.e.,β 2 = 0,η = 0, and κ 2 = 0, then the transmission diagram in Figure 1 is reduced to the transmission diagram in Figure 4. Using the parameter transformation in Supplementary Section S1, the basic reproduction number of the malaria model using transmission diagram in Figure 4 is given by To give further interpretation of R standard 0 , Equation (6) can be rewritten as follows: It is clear to see that R standard 0 is a result of multiplication between the number of new infected humans, new infected mosquitoes, and the lifetime of the exposed and latent classes. It can be seen that the saturated parameter of the treatment term does not appear in R standard 0 .

2.
Model with reinfection only. When the malaria model in the transmission diagram in Figure 1 includes reinfection only, without relapse and recrudescence, then the transmission diagram becomes that depicted in Figure 5. Calculating the basic reproduction number of the non-dimensional form of the model from the transmission diagram in Figure 5, we have It can be seen that R reinfection 0 = R standard 0 , which means that reinfections do not change the standard basic reproduction number.

3.
Model with relapse only. With the same approach as before, when reinfection and recrudescence are not involved,β 2 = 0 and κ 2 = 0. Based on this, the transmission diagram in Figure 1 changes to Figure 6. The basic reproduction number of a non-dimensional model based on transmission diagram in Figure 6 is given by Since . This reduction was due to the dormant period experienced by infected individuals in the hypnozoite phase, which made them unable to directly infect healthy mosquitoes. As previously mentioned, malaria infection by Plasmodium Vivax and Ovale can result in a long dormant period of up to 2-3 years.

4.
Model with recrudescence only. When relapse and reinfection are not involved in the original model (Equation (S1) in the Supplementary File), then we haveβ 2 = 0,η = 0. Hence, the transmission diagram in Figure 1 is reduced to the one in Figure 7.
The basic reproduction number of the non-dimensional model based on transmission diagram in Figure 7 is given by It can be seen that since we may conclude that recrudescence will increase the standard basic reproduction number.
To summarize the biological meaning of the above special cases of the basic reproduction number, we have the following relations: 1.
The basic reproduction number when all recurrence phenomena are constructed as a multiplication between infection in humans, infection in mosquitoes, lifetime period of class e, and the ratio between the in-and outflow of class l. We call this basic reproduction number the standard basic reproduction number.

2.
The existence of reinfection phenomena does not change the size of the standard basic reproduction number. It means that increasing reinfection the rate will not affect the size of the standard basic reproduction number. However, it will increase the endemic size and the possibility for the existence of multiple endemic conditions in the environment. We discuss this in the next section of this article. 3.
The existence of relapse phenomena will reduce the size of the standard basic reproduction number. This is highly related to the duration of the dormant period of the hypnozoite inside the human body. 4.
The existence of recrudescence phenomena will increase the size of the standard basic reproduction number.

Existence of the Endemic Equilibrium
In this section, we analyze the existence of the non-trivial equilibrium points, namely the malaria-endemic equilibrium point. This equilibrium differs from the malaria-free equilibrium, wherein malaria-endemic equilibria coexist between all classes, which means that malaria may still exist among the human and mosquito populations. The malariaendemic equilibrium is given by where and i * is a positive solution of the following five-degree polynomial: , and b 5 can not be shown in this article due to the complexity of their forms. Since b 6 > 0 and b 1 < 0 ⇐⇒ R 0 > 1, we always have at least one malaria-endemic equilibrium point if R 0 > 1. Furthermore, since the malaria-endemic equilibrium point is taken from a five-degree polynomial, the possibility that the malaria-endemic equilibrium point is not unique must be considered. For this purpose, we find that our model may have an endemic equilibrium even though R 0 < 1. Furthermore, we also have shown that it is possible to find one or three malaria-endemic points when R 0 > 1, and none or two malaria-endemic equilibrium points when R 0 < 1. Please see Supplementary Section S2 for the mathematical details of these results.

Bifurcation Analysis
The backward bifurcation phenomenon is believed to be one of the reasons why several types of diseases with relapse or reinfection periods are difficult to control. Many authors have studied a backward bifurcation phenomenon in the malaria model. In the non-dimensional version of our proposed malaria model, our previous analysis shows the possibility of the existence of backward bifurcation, which is indicated by the existence of an endemic equilibrium when R 0 < 1 under a specific condition (please see Supplementary Section S2). Epidemiologically, this means that the basic reproduction number can not be the only threshold in the determination of whether malaria exists or is extinct from the population.
In this section, our aim is to investigate the possible bifurcation phenomena in our model and determine another threshold that can be used to determine what type of bifurcation can arise depending on this new threshold. In order to do so, we use the well-known Castillo-Song bifurcation theorem [30]. For further details on the analysis of the bifurcation type from the non-dimensionalized malaria model, please see Supplementary Section S3. From the analysis that has been performed, we find that our model may undergo a backward bifurcation phenomenon at R 0 = 1.
The appearance of backward bifurcation phenomena (see Supplementary Section S3) indicates another threshold that can determine the existence or extinction of malaria from the field other than R 0 , namely the ω * . The emergence of the backward bifurcation phe-nomenon indicates that concluding the basic reproduction number as the main threshold in malaria control is no longer sufficient. This is because even though the basic reproduction number is already less than one, it is still possible that malaria is endemic in the community. In the following paragraph, we discuss this theorem's mathematical and biological consequences, backward and forward bifurcation, and hysteresis phenomena.
From the Castillo-Song bifurcation coefficient in Supplementary Section S3, we can see that the sign of A will determine the type of bifurcation that will appear from our model at R 0 = 1. When A < 0, we will have a forward bifurcation phenomenon. On the other hand, when A > 0, backward bifurcation is the phenomenon that we will obtain. To visualize this result, we use the parameter in Table 1, we calculate ω * , and we obtain ω * = 65.64. We choose four different values of ω, and the bifurcation diagrams are shown in Figures 8 and 9. From Figure 8, we show that the malaria-free equilibrium point is always stable when R 0 < 1 and unstable when R 0 > 1. When the saturation parameter ω, which presents the capacity of the hospital, is more significant than ω * , then backward bifurcation arises (Figure 8b). This indicates that malaria does not go extinct, even when R 0 is already smaller than one. In this case, we have a large size of the malaria-endemic equilibrium that is also stable for some values of R 0 < 1. When the number of the bed capacity is significant in the hospital (presented by a much smaller value of ω compared to ω * ), then forward bifurcation arises at R 0 = 1. Interesting results are shown in Figure 9. Even though we set the value of ω smaller than ω * , which will give us a forward bifurcation phenomenon, as we can see from Figure 9a that we may still have a stable endemic equilibrium for some values of R 0 < 1. The reason for this is that the forward bifurcation phenomena arise coupled with a hysteresis. As a result, we can have multiple stable malaria-endemic equilibria for some values of R 0 > 1, and bistability between the malaria-free equilibrium and the malaria-endemic equilibrium for some values of R 0 < 1.
In comparison to the simulation in Figure 9a where the bistability between the malariafree and malaria-endemic equilibrium points could arise when R 0 < 1, the result in Figure 9b is different. We can see from Figure 9b that we never had bistability between malaria-free and malaria-endemic equilibria. We only obtain bistability between two malaria-endemic equilibrium points for some value of R 0 > 1. The phenomena above give an indication of the difficulty of creating malaria-controlling policy if only relying on the magnitude of the basic reproduction numbers. The emergence of a stable endemic point is still possible when backward bifurcation occurs because the basic reproduction number can lead to policies that are not optimal. Furthermore, the emergence of forward bifurcation with hysteresis allows the emergence of a stable endemic point of a large size even though the basic reproduction number is larger but close to one.
From the illustration of the impact of hospital capacity on the type of bifurcation of our model in Figures 8 and 9, we can see that a smaller bed capacity in the hospital (larger ω) will trigger the appearance of backward bifurcation. Hence, more bed capacity in the hospital is better for malaria eradication policy. Furthermore, our aim in this article is also to find the impact of reinfection (β 2 ), relapse (η), and recrudescence (κ 2 ) on the type of bifurcation on our proposed model. To achieve this aim, we simulate the impact of β 2 , η, ω, and κ 2 on the sign of A. If A > 0, then backward bifurcation appears, and if A < 0, then forward bifurcation appears. Using the same parameter as shown in Table 1, the results are given in Figure 10. From Figure 10a, the larger the proportion of treated individuals that failed the treatment (larger κ 2 ), then the higher the chance that our model exhibits a backward bifurcation. On the other hand, there is an interval for β 2 depending on the value of κ 2 such that we can have a forward or backward bifurcation in our model. This means that a larger reinfection rate does not always trigger a backward bifurcation phenomenon. Figure 10b shows the impact of the treatment saturation parameter (ω) and relapse rate (η) on the bifurcation type of the proposed model. A smaller capacity of beds in the hospital (larger ω) will increase the chance of the appearance of backward bifurcation. On the other hand, with a larger relapse rate value, there exists a larger possibility that we can avoid a backward bifurcation. This illustration shows that relapse holds an important role in determining whether malaria will exist or become extinct when the basic reproduction number is already less than one, since it will trigger a condition where a basic reproduction number smaller than one will not guarantee the extinction of malaria.
(a) (b) Figure 10. Plot of A as a function between rate of reinfection in human (β 2 ) and recrudescence proportion (κ 2 ) in (a) and between treatment saturation rate and rate of relapse (η) (b).

Global Sensitivity Analysis
We perform a global sensitivity analysis using the combination of Latin hypercube sampling (LHS) and partial rank correlation coefficient (PRCC) [53,54] to determine the influential parameters of the model as in [55][56][57]. We measure against the increasing number of infected individuals which is Figure 11 shows that the parameters β 1 , β 4 , and µ v are the most influential parameters, where the first two have positive relationships and the latter has a negative relationship. The results are similar to the other work, where the transmission-related parameters and the death rate of mosquitoes are the most influential parameters [55]. However, our results provide detailed information regarding which transmission parameters govern the dynamics of infected individuals: the transmission rate from mosquito to susceptible humans (β 1 ) and the transmission rate from human to mosquitoes (β 4 ). On the other hand, the influence of the transmission rate from mosquitoes to dormant individuals is not as strong as the previously mentioned transmission parameters. The results indicate that the transmission from mosquitoes to susceptible human plays an essential role in determining the dynamics of malaria transmission. Furthermore, the death rate of mosquitoes, which are generally influenced by climatic factors, can reduce the number of infected individuals. This suggests that further work to analyze the effects of climatic factors is important and is the subject of future work.  Figure 11. PRCC values when measured against an increasing number of infected individuals over time.

Characterization of the Optimal Control Problem
In this section, we modify our proposed model in Equation (S1) of the Supplementary File to include three different interventions which depend on time, namely the use of bed net (u 1 (t)), hospitalization (u 2 (t)), and fumigation (u 3 (t)).

1.
Use of a bed net. Use of bed net is reportedly successful in reducing malaria incidence worldwide [58]. Bed nets provide protection to humans from the bite of a mosquito. Let us assume u 1 represents the proportion of the human population who use a bed net. Hence, u 1 N and (1 − u 1 )N represent the total human population who use and do not use a bed net, respectively. The successful transmission rate for humans who use bed nets now read as δβ i , where 1 − δ presents the efficacy of bed nets in reducing the number of successful bites. Note that in δ ∈ (0, 1), a smaller δ represents a bed net with better quality. Based on this assumption, the total number of new infections for non-users of bed nets in susceptible populations is given by and the total number of new infections for the bed net users is given by Therefore, total of new infections in the susceptible population is given by Note that if the entire human population used a bed net and the quality of the bed net could provide 100% protection against mosquitoes' bites (δ = 0), then no infections would occur in the field (in this case, we have f (u 1 , β 1 ) = 0). On the other hand, if all humans used a bed net but the efficacy level (protection from mosquito bites) is not 100%, then there is still a possibility that new infections occur, which given by Furthermore, when not all humans use a bed net, but the efficacy level of the bed net is 100%, then the total number of new infections is given by A similar approach is applied to the reinfection term and the new infection in mosquito population term, which involveβ 2 ,β 3 , andβ 4 . Note that when u 1 = 1, the infection term is reduced into the standard model in Equation (S1) of the Supplementary File. For the sake of simplification, we use another interpretation of bed net use, the term ((1 − u 1 ) + u 1 δ), as follows. The term ((1 − u 1 ) + u 1 δ) can be rewritten as where ζ = 1 − δ represents the effective bed net utilization rate. If ζ = 0 or (equivalently) δ = 1, then the bed net is useless, and regardless of the number of people who use the bed net, there will not be any impact on malaria prevention. In contrast, if ζ = 1, which is equivalent to the condition δ = 0, then the bed net can always provide protection to humans from mosquito bites. The larger the utilization rate u 1 ζ, the stronger the impact of bed net usage in the malaria prevention program. Therefore, instead of using the expression as in (13), we use the following expression to show the impact of bed net usage For another modeling approach on the application of bed nets in malaria intervention, please see [17,59].

2.
Hospitalization. In endemic malaria areas, hospitalization is the most frequently used outbreak control effort. However, this effort is difficult to execute continuously at a high intensity. Therefore, instead of using the constant hospitalization rate of γ, we use the new term u 2 (t), which represents the time-dependent treatment rate.

3.
Fumigation. For many types of vector-borne disease, including malaria, vector control programs are the most common intervention to control the spread of the disease. Hence, we introduce u 3 (t) as the additional death rate of the mosquito population due to fumigation, where the intervention depends on time.
Based on the above descriptions of three time-dependent interventions, the original model in Equation (S1) of the Supplementary File now reads as shown in Equation (S19) of the Supplementary Section S4. Our aim is to minimize the number of infected humans with an optimal intervention (as low as possible). Mathematically, the task reads as minimizing the following objective functional: where ω 1 , ω 2 , ω 3 , and ω 3 are the weights of the objective function for E, D, L, and I, respectively. Furthermore, ϕ j for i = 1, 2, 3 are the weight parameters for the control variables. ω i and ϕ j are the weight parameters that can balance each term in J . We use the quadratic cost function for J to describe the cost of control efforts. quadratic form is a common way to express cost functions in many mathematical epidemiology models with optimal control [4,17]. This quadratic function can describe a nonlinear cost increase related to the implementation of control efforts in the field.

Optimal Control Characterization
Using Pontryagin's maximum principle [60], the optimal solutions of u 1 , u 2 , and u 3 are given by u * 1 = min 1, max 0, u * 3 = min 1, max 0, Please see the Supplementary Section S4 for details on the derivation of the optimal control characterization.

Optimal Control Simulation
This section presents the numerical optimal control simulation to numerically solve the model with control and without control. Employing the backward and forward sweep as described in [61], we obtained the optimality of the model, which is comprised of the model without control, the adjoint system (see Supplementary Section S4), and the optimality conditions (17). To ascertain which strategy or combination gives the efficient methods of controlling malaria disease spread, we considered the following strategies/scenarios for our simulation results as enumerated below:

Strategy 1. Single Intervention: Use of Bed Net Only
The intervention of bed net usage only was simulated against all the model state variables, and the results are given in Figure 12. We observe in Figure 12 that the use of bed net as a control strategy only impacts the dormant population and the total population consisting of the exposed, dormant, latent, and infected individuals, while Figure 12a-d shows that the control u 1 does not have an impact on the control of the exposed, dormant, latent, and infected population, respectively. Figure 12f shows the control profile for strategy 1.

Strategy 2. Single Intervention: Use of Hospitalization Only
The use of hospitalization as the only control measure when simulating the model with and without optimal control is depicted in Figure 13. We see that this control strategy has no impact on the exposed (Figure 13a) and latent (Figure 13c) population as well as the control profile, as indicated in Figure 13f. It can be seen that hospitalization has a negligible impact on the dormant individuals and significantly affects both infected humans and the total population (that is, the total population of E, D, L, and I) as depicted in Figure 13d,e. The biological implication of hospitalization, as seen in Figure 13c, shows a decrease in the number of infected individuals with an increase in time. In other words, there is a reduction in the number of humans infected with malaria disease.  Simulation results for strategy 1 shows the dynamics of exposed compartment, dormant compartment, latent compartment, infected compartment, total of exposed, dormant, latent, and infected population, and the dynamic of control u 1 for subfigure (a)-(f), respectively. Simulation results for strategy 2 shows the dynamics of exposed compartment, dormant compartment, latent compartment, infected compartment, total of exposed, dormant, latent, and infected population, and the dynamic of control u 2 for subfigure (a)-(f), respectively. 6.3.3. Strategy 3. Single Intervention: Use of Fumigation Only Figure 14 shows the implementation of fumigation as the only control method with its control profile. The optimal control time series trajectories for the dormant population and the total combined population show a great decline in the number of individuals with an increase in time (See Figure 14b,e). Additionally, we observe that strategy 3 has a small impact on controlling exposed, infected, and latent individuals, as seen in Figure 14a,d,c, respectively. Simulation results for strategy 3 shows the dynamics of exposed compartment, dormant compartment, latent compartment, infected compartment, total of exposed, dormant, latent, and infected population, and the dynamic of control u 3 for subfigure (a)-(f), respectively.

Strategy 4. Double Intervention: Combination of Hospitalization and Fumigation
While Figure 15 shows the implementation of the combined strategy consisting of hospitalization and fumigation, Figure 15f shows its control profile. Clearly observable in Figure 15b,e, the implementation of a combined strategy reduced the number of the dormant, infectious, and the total combined populations while having a lesser impact on the exposed/latent individuals, as seen in Figure 15a,c. We can deduce that using the combined strategy is the best control measure to eliminate malaria infections. The shaded area in the control profile shows from our modeling that malaria can be controlled when u 2 = 0.2 and u 3 = 0.4, rather than when the values are equivalent. However, this signifies that the spread of malaria infection in any community with the disease can be prevented or eradicated by implementing both fumigation to control/decrease mosquito birth rate and hospitalization of infected individuals. Additionally, we observe that strategy 4 has a small impact in controlling exposed, infected, and latent individuals, as seen in Figure 15a,d,c respectively. Simulation results for strategy 4 shows the dynamics of exposed compartment, dormant compartment, latent compartment, infected compartment, total of exposed, dormant, latent, and infected population, and the dynamic of control u 1 and u 2 for subfigure (a)-(f), respectively.

Cost-Effectiveness Analysis
Here, we present the cost-effectiveness analysis for the model with and without control, intending to determine the best control strategies among all the scenarios. The total averted infection and total cost for each scenario are given in Table 2. The costs of the scenarios are obtained from the numerical output for the objective function Equation (15), arranged in incremental order. ICER = Difference in total costs by strategies A and B Difference in the total number of averted infection by strategies A and B . (17) Using the values in Table 2, we calculate the incremental cost-effectiveness ratio (ICER) by comparing each successive strategy with each other by employing (17) as defined in [62]. We arrange the strategies in increasing order from the lowest to the highest ICER values (see Table 2). Note that in order to make decisions, higher ICER values are eliminated in each step [62]. First, we calculate and compare the ICER values of strategies 1 and 3, where the ICER for strategy 1 is calculated by comparing the baseline total averted infections and total cost as given below Considering the results from the calculations of ICER(1) and ICER(3), we reject ICER(1), as it has a higher cost than ICER(3). Next, we compute and compare the ICER values for strategies 3 and 2. Using the total averted infections and total cost for strategy 2, we find the new ICER value using strategy 2 as a baseline while comparing it to the value of ICER(3) as follows: Similarly, comparing ICER(2) and ICER(4) using Equation (19), we eliminate strategy 4 since strategy 2 is cheaper than strategy 4. We conclude that to reach elimination/control of malaria, it is most cost-effective to first implement strategy 2. We repeat the above iteration to help establish the next most cost-effective strategy. Using the simulation results from the re-computation, we find that strategy 3 is the subsequent most cost-saving intervention after strategy 2 followed by strategies 1 and 4. Our results suggest that strategy 4 is the least cost-effective intervention required for malaria control as projected by our modeling framework.

Discussion And Conclusions
There are many mathematical models discussing malaria transmission, such as [4,14,15,17,21,23] and many more, but none of the mentioned references discuss the impact of relapse, recrudescence, and reinfection all together in one model. Recently, the authors of [29] introduced a complex nine-dimensional system of ordinary differential equations model to describe malaria transmission under the impact of relapse, recrudescence, and reinfection. They found that backward bifurcation phenomena can be avoided in the absence of reinfection. Different than [29], here we introduced a mathematical model for malaria transmission considering the impact of relapse, reinfection, recrudescence, and limited hospital bed capacity which affect the recovery rate. Malaria incidence data from the Wee Luri district in Central Sumba, Kupang, Indonesia, were used to estimate the model parameters. We investigated the existence and the local stability of the malaria-free equilibrium point utilizing the concept of the basic reproduction number [50]. We found that malaria can be driven to extinction if the basic reproduction number is less than one and always exists if the basic reproduction number is larger than one. The malaria-endemic equilibrium can not be shown explicitly. A gradient analysis and Descartes's rules of sign change were then used to show the malaria-endemic equilibrium when the basic reproduction number is larger than one. The concept of center manifold theory was used to establish the local stability criteria of the malaria-endemic equilibrium. Using the Castillo-Song bifurcation theorem [30], we determined another threshold other than basic reproduction number which can determine the existence or extinction of malaria from the system. Our analysis shows that our model may exhibit backward or forward bifurcation with hysteresis, triggering multiple malaria-endemic equilibria around the basic reproduction number equal to one. The existence of backward bifurcation makes malaria eradication more difficult since a basic reproduction number of less than one is no longer a sufficient indicator of the endemicity level.
From sensitivity analysis on the second threshold (ω * ) which may determine the type of bifurcation on our model, our findings suggest that backward bifurcation may likely occur given the following: (1) the hospital bed capacity is sufficiently small, (2) the relapse rate from dormant to latent compartment decreases, and (3) the proportion of treated individuals who failed the treatment is sufficiently large. Furthermore, we find that a more significant value of the reinfection rate cannot always guarantee the existence of a backward bifurcation phenomenon. We also studied the global sensitivity analysis of the total number of infected individuals (I). It was observed that the primary infection of humans and mosquitoes and the death rate of mosquitoes are the most influential parameters for determining the increase in the number of newly infected individuals. Hence, reducing the infection rate and increasing the mosquito's death rate could be an option for reducing the number of infected individuals. Additionally, an increase in reinfection, recrudescence, and relapse increases the number of newly infected individuals.
Further, an optimal control problem was formulated by modifying the malaria model without control, incorporating three control variables: the use of a bed net, hospitalization, and fumigation. The model consists of these interventions as a time-dependent variable to minimize the number of infectious individuals and the cost of implementing the intervention strategy. The characterization of the optimal control problem was done using the Pontryagin maximum principle [60]. We solved our optimal control problem using the forward-backward iterative method until the convergence criteria are achieved. Numerical simulations for each control variable depict the fact that the use of bed net, hospitalization, and fumigation can decrease the spread of malaria (see . However, we can see from Figure 12 that the use of a bed net is insignificant in reducing the number of infectious individuals when compared to hospitalization and fumigation intervention. Hence, we combine hospitalization and fumigation as the fourth strategy, and the result can be seen in Figure 15. Our sensitivity analysis indicates that the combination of hospitalization and fumigation is the most cost-effective strategy. Implementing hospitalization is the most cost-effective strategy if a single intervention is preferable to a combination of two or more interventions. Another interesting observation is that if the number of infected individuals already decreases, then the intensity of the intervention could be reduced to obtain a lower implementation cost. Therefore, we may conclude that intervention with hospitalization, fumigation, and bed nets can be used as a critical intervention strategy for controlling malaria, although bed net use is not as effective as the other two interventions. Conclusively, we predict that implementing the interventions which comprise the simultaneous use of hospitalization, fumigation, and bed nets will play a crucial role in the control of malaria disease. However, our model has some limitations. As mentioned by several authors [4,63], vector bias has an essential role in determining the success of malaria eradication programs. Another important factor that needs to be discussed is the seasonal impact on malaria transmission, as mentioned by many authors [64,65]. Including these factors in our model will be the future quest in order to obtain a better understanding of malaria transmission.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.