An Optimal Control Model to Understand the Potential Impact of the New Vaccine and Transmission-Blocking Drugs for Malaria: A Case Study in Papua and West Papua, Indonesia

Malaria is one of the major causes of a high death rate due to infectious diseases every year. Despite attempts to eradicate the disease, results have not been very successful. New vaccines and other treatments are being constantly developed to seek optimal ways to prevent malaria outbreaks. In this article, we formulate and analyze an optimal control model of malaria incorporating the new pre-erythrocytic vaccine and transmission-blocking treatment. Sufficient conditions to guarantee local stability of the malaria-free equilibrium were derived based on the controlled reproduction number condition. Using the non-linear least square fitting method, we fitted the incidence data from the province of Papua and West Papua in Indonesia to estimate the model parameter values. The optimal control characterization and optimality conditions were derived by applying the Pontryagin Maximum Principle, and numerical simulations were also presented. Simulation results show that both the pre-erythrocytic vaccine and transmission-blocking treatment significantly reduce the spread of malaria. Accordingly, a high doses of pre-erythrocytic vaccine is needed if the number of infected individuals is relatively small, while transmission blocking is required if the number of infected individuals is relatively large. These results suggest that a large-scale implementation of both strategies is vital as the world continues with the effort to eradicate malaria, especially in endemic regions across the globe.


Introduction
The primary cause of any infectious disease is the spread of pathogens such as bacteria, viruses, parasites, or fungus. Some of these diseases can spread from person to person, either through direct or indirect contact. A disease that spreads in the human population with an intermediate vector is called a vector-borne disease. According to the world health organization (WHO) [1], vector-borne infection contributes to over 17% of total infectious disease cases worldwide, with an estimated 700,000 deaths annually. Malaria is a vectorborne disease that threatens millions of people every year. Most of the cases come from Africa, southeast Asia, the eastern Mediterranean, and West Pacific [2]. Globally, there were over 600,000 deaths reported out of over 240 million detected cases in 2020. In Indonesia, over 74% of all detected cases are from the provinces of Papua and West Papua [3].
Malaria spreads through the bite of female Anopheles mosquitoes from previously infected individuals with Plasmodium (P). The five types of Plasmodium that cause malaria are P. falciparum, P. vivax, P. malariae, P. ovale, and P. knowlesi [4]. Globally, among all these types of Plasmodium, P. falciparum is primarily commonly found in African countries, while P.vivax exists in other parts of the world (that is, outside Africa) [5]. Several forms of intervention to control the spread of malaria have been introduced for many years worldwide. Most forms of intervention aim to control the population of mosquitoes in the field, such as with the use of fumigation, larvacide, or insecticide-treated bed nets [6]. Other than vector control, malaria prevention also focuses on developing a new vaccine to prevent new infections in the human population. Based on the life cycle of Plasmodium in the human body, the malaria vaccine is divided into three types, namely pre-erythrocytic vaccine, blood-stage vaccine, and transmission-blocking vaccine [7]. The pre-erythrocytic vaccine was designed to eliminate sporozoites immediately after mosquitoes inject Plasmodium into the human body. It will also block the Plasmodium from going to the human liver. According to the WHO report [5], RTS, S/AS01 is one type of pre-erythrocytic vaccine that is recommended by WHO to kill P. falciparum in children population. The blood-stage vaccine is introduced to target the asexual stage of Plasmodium (merozoites) in human red blood cells. On the other hand, the transmission-blocking vaccine is given such that further infection can be avoided [8]. Other than the transmission-blocking vaccine, a type of transmission-blocking drug has been used to eliminate the sexual stages of Plasmodium in human or mosquito bodies.
Mathematical models have been used for many years to understand the mechanisms of malaria transmission, since such model approaches can be used to give a visual interpretation of any possible intervention that can be implemented in the field to control malaria transmission. These approaches provide a scientific background before a final decision from the government should be taken. From the early work by Ross [9] in 1911, many mathematical models were introduced by authors to help a better understanding of malaria transmission. Macdonald in 1957 [10] used his model to estimate the infection and recovery rates of malaria. He found that reducing the number of mosquitoes in an endemic area is an inefficient malaria control strategy. In the early 1980s-1990s, Anderson and May [11] and Aron and May [12] constructed their malaria model based on the assumption that immunity to malaria is independent of the duration of exposure. Okosun et al., in [13], concluded from their mathematical model that a combination between insecticides and transmission-blocking treatment is the most cost-effective interventions to control malaria. An optimal vaccination and bed net mathematical model have been introduced by Prosper et al. in [14]. Their analytical result reveals that increasing the case detection strategy may reduce the chance of backward bifurcation phenomena in their model. An analysis of the potential impact of pre-erythrocytic vaccine from clinical data was discussed by the author in [15]. Similar to Prosper et al. in [14], Woldegerima et al. in [16] also found a possible backward bifurcation from their model on the impact of transmission-blocking drugs. Their model projects an approximately 82% death rate is reduced by 2035 if 35% of the population in Sub-Saharan Africa receives a transmission-blocking drug with an efficacy of 95%. Kuddus and Rahman in [17] analyze the dynamics of malaria using incidence data in Bangladesh from 2001 to 2014. They found that as infection rate has the greatest impact on the basic reproduction number compared to other model parameters, it is important to reduce this infection rate, such as using insecticide-treated bed nets, spraying insecticides, clearing stagnant water, etc. In a more recent mathematical model, authors include more recent facts on malaria transmission such as the effect of vector-bias [18], asymptotic carriers [19], age-structured [20], competitive strains [21], seasonal factor [22,23], and coinfection of malaria with COVID-19 [24]. Furthermore, intervention models also have been widely introduced by authors, such as the use of fumigation [25], insecticide-treated bed nets [26], and vaccines with waning immunity [27], or transmission-blocking vaccines [28].
The authors in [29] concluded that a better understanding of Plasmodium development in the pre-erythrocytic stage could lead to the discovery of new antigenic targets for enhanced immunization strategies. These findings may aid in the development of new malaria immunization techniques that are more successful and useful. In 2015, Glaxo-SmithKline's (GSK's) RTS, S/AS01 pre-erythrocytic vaccine received a positive scientific response on the quality of this vaccine in combating malaria transmission [30]. Based on the above description, we consider it essential to see the potential impact of the preerythrocytic vaccine (RTS, S/AS01) and transmission-blocking drug as a combination of control to reduce the spread of malaria using a mathematical model approach.
Our model divided the susceptible population based on a condition, whether they use a pre-erythrocytic vaccine or not. Furthermore, we also consider the possibility that the transmission-blocking drug cannot kill all Plasmodium in either the sexual or asexual stage in the human body. We used incidence data from Papua and West Papua to calibrate our model by determining the best-fit parameter on our model. Sensitivity analysis on the basic reproduction number and the numerical simulation on the dynamics of each compartment are also conducted to see the impact of pre-erythrocytic vaccine and transmission-blocking drugs in the malaria control strategy. An optimal control model is then analyzed to determine the best possible scenario for combining pre-erythrocytic vaccine and transmission-blocking drugs. The novelty of our model lies in the originality of the model, where we discuss the potential impact of the pre-erythrocytic vaccine. Furthermore, we also use incidence data from two provinces in Indonesia, namely Papua and West Papua, to calibrate our model parameters. These incidence data have never been used in any malaria mathematical model.
The organization of this paper is as follows. In Section 2, we construct our model based on our assumptions, and estimate the parameter values on our model by fitting the output of our model with incidence data in Papua and West Papua. Dynamical analysis regarded the positiveness of the solution of our model, the existence and stability of equilibrium points, and the controlled reproduction number discussed in Section 3. Sensitivity analysis is given in Section 4. Existence of optimal solution and the characterization of the optimal control problem are given in Sections 5 and 6, respectively. Finally, some conclusions are given in Section 7.

The Model
Based on the malaria status in the human and mosquitoes population, we divided the total human population (denoted by N) into eight compartments as follows:

1.
Susceptible without vaccine (S 1 ) consists of a group of individuals susceptible to malaria who have not received a pre-erythrocytic vaccine yet.

2.
Susceptible with a vaccine (S 2 ) consists of a group of individuals who are also susceptible to malaria but have already received a pre-erythrocytic vaccine. 3.
Exposed without vaccine (E 1 ) consists of a group of newly infected individuals from S 1 who have not yet gotten the pre-erythrocytic vaccine. We assume that individuals in this compartment are in the leaver stage. Hence, although these individuals do not show any symptoms yet, we assume that they can transmit Plasmodium to susceptible mosquitoes. 4.
Exposed with vaccine (E 2 ) consists of a group of newly infected individuals from S 2 . Although these individuals have already gotten vaccinated, the description is still the same with E 1 .

5.
Infected (I) consists of a group of individuals who have already gotten infected by malaria and show their symptoms. 6.
Infected individuals undergo transmission-blocking treatment (T), defined as a group of individuals who already get infected, show symptoms, but get a transmissionblocking treatment. We assume that this treatment can kill the sexual Plasmodium (gametocytes). 7.
Recovered but carrier (R 1 ) consists of individuals who recovered from malaria (do not show symptoms anymore) and have a temporal immunity but still have asexual Plasmodium inside their bodies. Hence, this group of individuals can still transmit Plasmodium to the mosquito.

8.
Fully recovered (R 2 ) consists of a group of individuals who recovered from malaria and succeeded in the transmission-blocking treatment process. Hence, unlike in R 1 , individuals in R 2 lack sexual and asexual Plasmodium in their blood. Therefore, R 2 compartment does not transmit malaria anymore.
On the other hand, we only divide the mosquitoes population (denoted by M) into two compartments, namely the susceptible and infected mosquitoes, denoted by V and W, respectively.
We use the following assumptions to construct the model: 1.
The rate of new individuals only came from newborns with a constant rate of Π h . We ignore migration from our model.
The infected individuals in E 1 , E 2 , I, T, and R 1 are capable to transmit Plasmodium to mosquitoes with a constant rate β v . 4.
The pre-erythrocytic vaccine is given to S 1 individuals with a constant rate u 1 to give temporal protection from mosquito bites that can lead to malaria infection. Hence, we assume that the transmission rate of S 2 is less than S 1 (β h2 < β h1 ).

5.
The pre-erythrocytic vaccine is not for a lifetime. Hence, after a period of α −1 , individuals in S 2 will return to S 1 . 6.
The transmission-blocking treatment is given to individuals in I with a constant rate u 2 to cure malaria and wipe out sexual and asexual Plasmodium from their blood. 7.
The transmission-blocking treatment is not always successful in curing infected individuals. 8.
The description of all parameters is given in Table 1 and assumed to be nonnegative.  [34] µ v Natural death rate of mosquitos 30 21 30 21 [35] The construction of the model is based on the transmission diagram in Figure 1. All newborns in human and mosquito populations are assumed to be susceptible, with a rate of Π h and Π v , respectively. On the other hand, we assume that each compartment has a natural death rate of µ h and µ v for human and mosquito populations, respectively. According to [32,33], the total inhabitants in Papua and West papua are 3,453,430 and 981,822, respectively. Our main purpose is to understand the potential impact of the pre-erythrocytic vaccine on reducing the possible transmission of malaria in the human population. This vaccine was designed to clean the sporozoite from the human body right after the mosquito injects the Plasmodium into the human body and block the sporozoites' invasion of the human liver cell. Hence, we include the rate of this pre-erythrocytic with a constant rate of u 1 to the susceptible population (S 1 ). This vaccine is not for a lifetime [37]. Hence, after a certain period of time, individuals in S 2 will return to S 1 . We denote this dropout rate with α. The transmission of malaria comes from the bite of infected mosquitoes to susceptible humans. We assume that the successful transmission rates of S 1 and S 2 are denoted by β h1 and β h2 , respectively. Let b be the average number of mosquitoes bite per month, then bβ h1 and bβ h2 are the averages of possible success transmission in S 1 and S 2 each month due to the bite of one infected mosquitoes, respectively. We use a ratio-dependent function to model our infection term. Therefore, the number of total infections in S 1 is given by . We assume that the total human population is constant. Therefore, bβ h1 N is constant, and we denote it by β h1 . Therefore, we have that the total number of new infections in the human population from the S 1 compartment is given by β h1 S 1 W. Using the same approach, the number of new infections in S 2 is given by β h2 S 2 W. Exposed individuals without vaccine (E 1 ) and with vaccine (E 2 ) move to Infected class (I) due to progression rate δ 1 and δ 2 , respectively. Note that due to the effect of pre-erythrocytic vaccine that can suppress invasion of sporozoites in hepatocytes [8], we have δ 2 < δ 1 . Without treatment, we assume that infected individual I can recover with a constant rate γ 1 .
Another essential factor that is considered in this article is the use of transmissionblocking treatment for infected individuals I. We assume that individuals in I get a treatment with a constant rate u 2 , which will transfer them into compartment T. This treatment is given in a duration of γ −1 2 . Hence, after a period of γ −1 2 , the result of this treatment should be evaluated. If the treatment successfully kills all sexual and asexual parasites, then they are transferred to the compartment of R 2 with a rate of pγ 2 . If the treatment only partially succeeds where only sexual parasites could be eliminated but not the asexual parasites, these individuals will go to R 1 with a rate of qγ 2 . Finally, if the treatment fails, then individuals in T go back to I with a rate of (1 − p − q)γ 2 . We assume that individuals in R 1 still could infect mosquitoes, but they are immune to the symptoms of malaria. Further, we also assume that they have a temporal immunity of κ −1 . On the other hand, individuals in R 2 cannot transfer Plasmodium to mosquitoes and have a temporal immunity of ϑ −1 .
Unlike the human population, the modeling of the mosquito population is not as complex. It only involves newborn Π v , natural death rate µ v , and infection. We assume that susceptible mosquitoes will get infected if they bite infected individuals in compartments E 1 , E 2 , I, T, and R 1 . Since the status of parasites in each mentioned compartment is in different stages, we assume the transmission rate for the mosquitoes are different based on their parasites status. Gametocytes are sexual forms of parasites produced by blood-stage merozoites. However, for P. vivax, gametocytes can be produced by liver stage merozoites [39]. Therefore, humans still in the pre-erythrocyte stage can also transmit susceptible mosquitoes. Exposed humans who received the pre-erythrocytic vaccine (E 1 ) will have fewer parasites in the liver, so the likelihood of liver-stage merozoites producing gametocytes will be less than that of exposed unvaccinated humans (E 2 ). In addition, according to [39], gametocytes can still be found in peripheral blood after several weeks of asexual parasite infection being cleared, so recovered humans (R 1 ) can still transmit susceptible mosquitoes. Humans in treatment can transmit susceptible mosquitoes because humans in treatment (T) are people who are still infected or have not been fully treated. Hence, the transmission rate of E 1 , E 2 , R, T to susceptible mosquitoes are given by ζ e2 β v , ζ e1 β v , ζ r β v , and ζ t β v , respectively. Note that ζ e1 , ζ e2 , ζ r , and ζ t are the correction parameters due to parasites' stage in human body, and satisfies the following inequality: Therefore, the total of new infection in mosquitoes population is given by From all mentioned parameters above, we must understand that some parameters, especially in the mosquito population, might depend on time or other factors such as weather. For example, mosquitoes are more active when the temperature is low, which means that mosquitoes will be more active in biting humans when the temperature is relatively low [40]. Based on this explanation, all parameter values in our model as in Table 1 are average values.
In many mathematical models for disease control, an optimal control approach is commonly used by many authors to understand the behavior of each intervention as a time-dependent variable under a specific budget limitation problem [41][42][43]. Our proposed model has two different forms of intervention: the pre-erythrocytic vaccine and transmission-blocking treatment, denoted by u 1 and u 2 , respectively. In this article, these two forms of intervention will be treated as a time-dependent variable to pursue the best strategy to suppress the spread of malaria. Hence, we have u 1 = u 1 (t) and u 2 = u 2 (t).
Based on assumptions above, we describe the dynamics of malaria under the effect of pre-erythrocytic vaccine and transmission-blocking treatment by the following system of ordinary differential equations: with positive initial conditions, and time measured in months. Next, we deduce the related cost which needed to be minimized in the context of malaria control strategy.

1.
Cost to implement pre-erythtocytic vaccine. The total cost to implement the preerythrocytic vaccine is given by where ω 1 is the weight parameter for pre-erythrocytic parameter and t f is the final time of simulation. The unit of ω 1 is 1 month 2 . In this article, we consider the nonlinearity term on this cost term as authors in [18,26].

2.
Cost to implement transmission-blocking treatment. Similar to u 1 , we also use a quadratic term for u 2 . Hence, the total cost of transmission blocking is given by where ω 2 is the weight parameter for transmission blocking. The unit of ω 2 is 1 month 2 .

3.
Cost related to the high number of infected individuals. Except the cost related to the implementation of pre-erythrocytic and transmission-blocking, the cost for malaria control strategy also comes from the cost related to the high number of infected individuals who were not treated (I and R 1 ). This cost is given by where ω 3 and ω 4 are the weight parameters for I and R 1 , respectively. ω 3 I and ω 4 R 1 respectively denote the cost associated with the high number of infected individuals, such as for maintaining health campaigns or any other related cost which comes as a consequence of the high number of infected individuals. The unit of ω 3 and ω 4 is 1 individual . Based on the above explanation, the cost function for our model is given by: There are no exact values of ω i for i = 1, 2, 3, 4. Choosing the values of ω i should balance the value of each components on J , since the range of u 1 and u 2 are very small compared to I and R 1 . Hence, we have to choose ω i so that no component dominates other components in the cost function. For example, in order to balance between ω 1 u 2 1 and ω 3 I, then ω 1 and ω 3 should satisfy ω 1 To simulate our optimal control problem, it is important that we use parameter values which can present the situation of malaria dynamics in Papua and West papua. Hence, it is important that we estimate our parameter values based on the incidence data in these areas.

Parameter Estimation
In this section, we discuss the incidence data that we use to estimate our parameters and the numerical scheme that has been used. The incidence data for malaria are taken from two provinces in Indonesia that have the highest number of malaria incidence every year, namely Papua and West Papua [44]. According to [32,33], the total number of inhabitants in Papua and West Papua was 3,435,430 and 981,822 in 2020, respectively. In both areas, the majority of malaria cases are caused by Plasmodium falciparum, and Plasmodium vivax, both of which cause routine health problems and morbidity in these areas [45].
The incidence data used in this article is the monthly new cases in both provinces, from January to December 2020. The data was collected by personal request from the Directorate of Prevention and Control of Vector and Zoonotic Infectious Diseases, Ministry of Health of the Republic of Indonesia. The data can be seen in Figure 2. In our model, we do not have a compartment that describes the monthly cases of malaria. Hence, we need to adapt our model to accommodate this. Hence, we transform the monthly cases in both provinces into the accumulated case. Next, we create a new compartment from our model, which describes the accumulated cases. From the transmission diagram given in Figure 1, the total detected cases are coming from the compartment T. Hence, the newly detected cases come from the intervention of transmission-blocking from I to T with a rate of u 2 . Hence, the dynamic of the monthly cases is given by Solving the above differential equation will give us the accumulated cases of malaria incidence. Our task is to minimize the Euclidean distance between c(t) from our numerical results with the incidence data, using the best-fit parameter of our model. This task can be formulated as follows.
where X is the set of best-fit parameters and initial conditions, while t f is the final time of available data. However, because the use of pre-erythrocytic vaccines had not been implemented in Papua and West Papua in 2020, parameter estimation is performed when the parameters u 1 , α, β h2 , δ 2 , and ζ e2 are zero and the S 2 and E 2 compartments are zero. The corresponding best-fit parameters are κ, ϑ, β h1 , b, δ 1 , p, q, u 2 , γ 1 , γ 2 , β v , ζ e1 , ζ t , ζ r , and all initial conditions of the system (1) (except S 2 (0) and E 2 (0)) and (3). We used a nonlinear optimization toolbox in MATLAB to conduct this parameter estimation, called fmincon function. To run the simulation, we need to give an estimation of an interval in which our parameter values must exist. Hence, we use the lower and upper bound of our parameter values, calculated as shown in Table 2. We estimate our parameter values for incidence data in Papua and West Papua, and the result is given in Table 1, while the fitted cumulated cases are in Figures 2 and 3. This finding indicates that malaria incidence will increase in Papua Province, which tends to the malaria-endemic equilibrium point. Meanwhile, malaria incidence in West Papua Province tends to the malaria-free equilibrium point. From our modeling analyses/results, we can deduce the following: 1.
Papua is more malaria-endemic than West Papua given the non-controlled basic reproduction number (R 0 | Papua = 1.75 and R 0 | West Papua = 1.53). See the formula of non-controlled basic reproduction number (R 0 ) in (12).

2.
Based on the value of p and q, we conclude that individuals in Papua have a slightly bigger chance of reaching the malaria elimination target if there is a continuous implementation of treatment efforts compared to West Papua.

3.
The infection rates from mosquitoes to humans (β h1 , and β h2 ) in West Papua is higher than in Papua. Hence, it is important to develop a media campaign that targets reducing the contact between humans and mosquitoes in West Papua than in Papua. Such campaigns may include the use of bed nets or mosquito repellent.   [34] µ v

Preliminary Results on the Positiveness and Boundedness of the Solution
The malaria model in system (1) involves human and mosquito populations. Hence, it is necessary to guarantee that all associated variables are positive with nonnegative parameters and initial conditions. Theorem 1. Let the initial conditions of all variables in system (1) be nonnegative as follows: Then the solution set is nonnegative for all t ≥ 0.
Since the initial conditions of R 1 , R 2 , and S 2 are nonnegative, then the first equation in system (1) satisfies: Hence, we have: Choosing the integrating factor as exp t 0 G(τ)dτ. gives us: Integrating both sides yields Therefore, Therefore, we can see that if S 1 (0) > 0, then we have that S 1 (t) ≥ 0 for all t > 0. A similar approach can be used to show that S 2 (t), E 1 1(t), E 2 (t), I(t), R 1 (t), R 2 (t), V(t), and W(t) are nonnegative for all t ≥ 0. This completes the proof.
It can be shown that our model in system (1) is well-defined biologically in the feasible region below. where

The Malaria-Free Equilibrium and the Controlled Reproduction Number
The first equilibrium point of malaria model in (1) is the malaria-free equilibrium, which is denoted by E * , and given by It can be seen that E * presents a condition where all infections have disappeared from both populations. In this condition, we also can see that the ratio of susceptible human beings with and without pre-erythrocytic vaccine is given by The purpose is, of course, to reduce this ratio, which is equivalent to increasing the number of individuals in S 2 . It can be seen that this ratio is only affected by three parameters, namely µ h , α, and u 1 . Of these three parameters, only the latter two are controllable. Therefore, we can see that reducing this ratio can be done by increasing the rate of vaccination, and the vaccine's duration could be prolonged. Hence, the higher the quality of the vaccine (i.e., its duration), the better.
We will now calculate the basic and controlled reproduction number of our model, which is denoted by R 0 and R C , respectively. The basic reproduction number defines the expected average number of secondary malaria cases caused by a single malaria-infected case during its infection period in a completely susceptible population. This threshold holds a vital role in many malaria models [48]. From its definition, R 0 helps us to understand the qualitative behavior of our model and whether the disease may persist or exist. Mostly, they found that malaria will disappear from the population if R 0 < 1, and exist if R 0 > 1. In several cases, mainly when backward bifurcation phenomena occur in their model [18,[49][50][51], it is still possible to find malaria even though R C is already less than one.
When no control intervention is given to our model, then system (1) reduced into: To find the basic reproduction number of the basic model (11), we use the well known next-generation matrix approach [52]. The respected basic reproduction number for the basic model (11) is given as follows.
With a similar approach, we can calculate the controlled reproduction number of our malaria problem. Implementing the next-generation matrix method into system (1), the controlled reproduction number is given by: where and As we can see from the expression on (13), R 0 is the combination of many paths of infection on our model. The explanation is given as follows.

1.
R E 1 0v shows the path of transmission for a new case in mosquito due to a bite from susceptible mosquitoes to an exposed human in compartment E 1 .

2.
R E 1 0h1 shows the transmission path which gives a new infection in humans without pre-erythrocytic vaccine (E 1 ) due to a bite from infected mosquitoes.

3.
R E 2 0v presents the transmission path for a new infection in the mosquito population after they bite exposed humans in E 2 .

4.
R E 2 0h2 presents the transmission path for a new infection in vaccinated human E 2 due to a bite from infected mosquitoes.

5.
R I 0v shows a transmission for a new infection in the mosquitoes population after biting an infected individual in I.

6.
R I 0h presents a transmission path for a new cases in I due to progression of E 1 and E 2 . 7.
R T 0v presents a transmission path for a new case in the mosquitoes population after biting an infected individual who is undergoing treatment (T). 8. R T 0h presents a transmission path for a new case in treated human population (T) due to the treatment rate from I.

9.
R R 1 0v presents a transmission path for a new case in the mosquito population after biting individuals in R 1 . 10. R R 1 0h presents a transmission path for new cases in the compartment of humans who partially succeed in conducting transmission-blocking treatment.

The Malaria-Endemic Equilibrium
The malaria-endemic equilibrium point of system (1) indicates the condition when malaria persists in the human and mosquito populations. The malaria-endemic equilibrium of system (1) will be determined when the number of infected individuals is not equal to zero. However, due to the complexity of the model, we can not explicitly show the form of the endemic equilibrium point as a function of other parameters. Hence, we write our endemic equilibrium point as a function of dependent variables I and W. The idea is as follows. We set the right hand side of system (1) equal to zero, and solve them backwardly with respect to each variables until we only have I and W left. This process will give us two different polynomial as a function of I and W. Hence, the malaria-endemic equilibrium of system (1) is given by: where with a = µ h 2 + (γ 1 + γ 2 + u 2 )µ h + γ 2 ((p + q)u 2 + γ 1 ) while I * and W * are taken from the positive intersection of the following polynomials: where a i (I), b i (W) for i = 0, 1, 2 has a complex form to be shown. We leave the existence analysis on this malaria-endemic equilibrium as an open problem to an interested reader. We show the existence of this malaria-endemic equilibrium numerically with the following scenario. To conduct this numerical experiment, we use the following parameter values: Using the above parameter values, we have the basic reproduction number of system (1) is 1.002169498, which is larger than one. Substituting the parameter values in Table 3 into G 1 and G 2 yield:  Next, we plot G 1 and G 2 in I − W plane. The result is given in Figure 4. Since we have an intersection between G 1 and G 2 in the first quadrant of I − W plane, then we have the malaria-endemic equilibrium point, which is given by: 20241, 200, 2, 1874, 862, 1960, 9898, 6870710, 151).

Global Sensitivity Analysis on the Basic Reproduction Number
In this subsection, we carry out a global uncertainty analysis of basic reproduction number for the malaria model without control. As we have already shown in the previous section, the basic reproduction number could determine whether malaria will persist or become extinct in the population. Hence, it is essential to analyze the sensitivity of the basic reproduction number with respect to the model parameters. Global sensitivity analysis is used to study the relative changes in epidemic model parameters output, giving input/altered model parameters for the dynamical systems under investigation. Thus, we enabled modelers to identify key model variables that require controlling for the particular disease. To carry out the simulation, we used global sensitivity analysis, which combines the Latin-hypercube sampling and Partial Rank Cross Correlation (PRCC) technique as given in [53] with similar analysis in [54]. Using R-software, we performed 1000 simulations per run, and examined and evaluated the PRCC of the model parameters concerning in R 0 . The PRCC indicates the degree of monotonicity between the model's parameters in the R 0 . Thus, juxtaposing the values of PRCC gives an apparent effect and contribution of each model parameter on R 0 in our malaria model. The results from the numerical simulations are given in Figure 5 and Table 4. The sensitivity of the parameters to the basic reproduction number (12) is proportional to the absolute value of PRCC. As we can see in Table 4, we observe Π h , β h1 , µ h , µ v , β v , Π v , γ 1 , and ξ r . From these parameters, we can see that β h1 , β v , and µ v are the most significant parameters that can be controlled. Increasing µ v and reducing β h1 and β v will reduce R 0 . It can also been seen in Table 4 that all the parameters have p-values < 0.05 and thus are said to be significant besides ξ e1 , κ, and δ 1 . Therefore, it is essential to introduce intervention to reduce the infection rate with vaccine in susceptible populations, or increasing the mosquito death rate with fumigation.

Local Sensitivity Analysis on the Model Variables
The sensitivity methods can be used on infectious disease models to determine which variable or parameter is sensitive to a specific condition. In our case, identifying the key critical parameters of system (1) is an effective way to study the qualitative behaviour of the parameters which govern the model. In our proposed model (1), we have 10 compartments x i for i = 1, 2, ..., 10 and 22 parameters k j for j = 1, 2, ..., 22. Then, the local sensitivities can be calculated using three different techniques: non-normalisations, half-normalizations, and full-normalizations. Firstly, the equation of non-normalization sensitivities is given by where S x i k j is measured as sensitivity coefficients of each x i with respect to each parameter k j . Then, the formula of half-normalization sensitivities is defined below: Finally, the equation of full-normalization sensitivities is defined by In this article, we only perform a full-normalization sensitivity analysis on our proposed model in (1) for best-fit parameter of Papua and West-Papua incidence data. The results are given in Figures 6 and 7. We can see from Figures 6 and 7 that the progression rate from E 1 to I (δ 1 ) is the most influential parameter for all compartments, except S 2 and E 2 . The second most influential parameter is the transmission-blocking parameter (u 2 ). In contrast to u 2 , we can see that the impact of pre-erythrocytic (u 1 ) is not as sensitive as u 2 . This result is similar to our previous sensitivity analysis on R 0 .

Existence of the Optimal Solution
In this subsection we state and prove the existence of solutions for the optimal control ODE system given in Equation (1), before proving the existence of solutions for the optimal control problem. Firstly, it is necessary to establish the boundedness of our malaria model over the modeling time horizon. Now, letŜ 1 ,Ŝ 2 ,Ê 1 ,Ê 2 ,Î,T,R 1 ,R 2 ,V, andŴ denote the super-solutions with respect to each state variable in Equation (1). Therefore, we obtain Writing system (24) in a vector notation format, we get From the above matrices, it can be deduced that our model is linear and satisfies the properties of being semi-positive definite on Ω. Then, it follows that the super-solutions of the system (1) are bounded.
Applying Theorem 4.1 on pages 68-69 of Fleming and Rishel [56], we then show the existence of the optimal control model. Supposed

Theorem 2. Given the objective functional
which has associated admissible control while subject to the model initial conditions given in Theorem 1. There exists an optimal control pair u * 1 , u * 2 which minimizes the functional Therefore, the conditions enumerated below should be satisfied by our model.

1.
The admissible control parameters and each model state variable are non-empty.

2.
Control H is convex and bounded.

3.
Right hand side of our system is continuous and bounded above by a linear function in the state variables and the control parameters. 4.
The integrand of J (u 1 (t), u 2 (t)) is concave on H.

5.
There exists positive constants b 1 , b 2 > 0 and τ > 1 such that J (u k (t)) satisfies Proof. The proof follows for the stated conditions as follows: 1.
Applying the methodology of Theorem 9.2.1 on page 182 of [57], the first criteria is fulfilled as the solutions of our model system exist and are bounded in ω as also shown by the existence of the super-solutions.

2.
Using the definition of H, we have that set as bounded and closed. (2) is linear in u 1 and u 2 . with the aid of the result of the theorem as in [57], we can then deduce that the right-hand side of our system is continuous and bounded.

5.
Following the fact that the model state variable in the objective function, that is, I and R 1 are bounded as shown in 3. Above, there exits positive constants ϕ 1 , ϕ 2 > 0 such that the sum of (I + R 1 ) ≤ ϕ 2 . If we set ϕ 2 = max n = 3, 4 · · · , we have that in which τ = 2. This implies that there is a positive constant ϕ 1 , ϕ 2 > 0 and a τ > 1 such that the integrand component of J satisfies From the above proof, we have that our system satisfies the necessary conditions in Theorem 2, which is completed here.

Characterization of the Optimal Control Problem
The Pontryagin's Maximum principle [58] allows us to utilize costate functions to transform the optimization problem to the problem of determining the point-wise minimum of u * 1 and u * 2 of the Hamiltonian. This Hamiltonian function is built using the cost function in (2) and the malaria model in (1). With this, we derive: where X = (S 1 , S 2 , E 1 , E 2 , I, T, R 1 , R 2 , V, W) and λ j for i = 1, 2, . . . 10 are the associated adjoints for the model variables S 1 , S 2 , E 1 , E 2 , I, T, R 1 , R 2 , V, and W, respectively. Next, we will calculate the optimality system of our optimal control problem. The following results is a direct consequences of application of the Pontriagin's Maximum Principle for bounded controls [59]. The adjoint system can be written as follows: and with transversality conditions λ j (T) = 0 for j = 1, 2, . . . 10. The optimality conditions requires that ∂H ∂u 1 = ∂H ∂u 2 = 0. Hence, according to our model and cost function, and following the lower and upper bounds, we have: u * 2 = min max 0,

Numerical Simulation of the Optimal Control Problem
From the previous analysis, our optimal control problem consists of the state system, which is the malaria model in (1) with initial condition X(0), the costate system in (29) with transversality condition λ j (T) = 0, and optimal conditions in (30). To determine the optimal trajectory of u * 1 and u * 2 , we have to determine the solution of the state and costate system. Since the state system has an initial condition, while the costate has the final condition, we cannot solve these systems directly forward in time. Thus, we will solve it using the "forward-backward sweep method" [59]. The algorithm is as follows. An initial guess for the control variables should be made. Then, with this initial guess, we solve the state system in (1) numerically forward in time using a Runge-Kutta method. Having the initial guess of control variables and the solution of the state system for t ∈ [0, T], we solve our costate system in (29) backward in time, also with the Runge-Kutta method. Having the solution of state and costate for t ∈ [0, T], we update the control trajectories using the equation in (30). This process is repeated until the convergence criteria are satisfied. In our numerical implementation, the convergence criteria are when the Euclidean distance between i th and (i + 1) th iteration for costate, state, and control parameters are smaller than our chosen tolerance.
Unless stated otherwise, to run the simulations in this section, we used parameter values as shown in Table 1, with the following initial conditions: and W(0) = 20%N, with N is the total of human population in Papua (3,435,430). The values assigned for the weight constant are ω 1 = 1, ω 2 = 1, ω 3 = 5 × 10 11 , and ω 4 = 10 12 . In this section, the impact of control trajectories with various scenarios is studied. We studied the impact on the change of the dynamic of each sub-population (with and without controls), the total of susceptible and infected populations, and the cost related to the scenario.

Different Combination of Interventions
For the first scenario of optimal control simulations, we conducted our experiment based on the combination of controls that have been used. Let Scenarios 1, 2, and 3 be scenarios when only pre-erythrocytic drugs are used (u 1 > 0, u 2 = 0), transmissionblocking drugs are used (u 1 = 0, u 2 > 0), and both drugs are used (u 1 > 0, u 2 > 0), respectively. The illustration on the dynamics of system (1) and the control profiles for Scenarios 1, 2, and 3 are given in Figures 8, 9, and 10, respectively. From Figure 8, we can see that with a profile of u 1 of Scenario 1 in Figure 8k, the number of infected individuals is always decreasing except E 2 , as a consequence of the existence of individuals who are already vaccinated (S 2 ). It is clearly observed that the pre-erythrocytic vaccine should be maintained at one from the start of the simulation for a long time period to increase the number of vaccinated individuals. With a large number of vaccinated individuals, fewer people will be infected by mosquitoes. Hence, the total number of infected individuals is always smaller than when no intervention is used. The total cost for Scenario 1 is 4.31889 × 10 12 , with an average of total avoided infected individuals (E 1 + E 2 + I + T) compared to without intervention scenario is 34,679 individuals. Please see Figure 11 for the dynamics of total susceptible and infected individuals, with and without control.
Simulation results of Scenario 2, when only transmission-blocking is used as an intervention strategy, are given in Figure 9. Similar to Scenario 1, the intervention of transmission-blocking successfully reduces the number of infected individuals significantly. The trajectories of u 2 for Scenario 2 can be seen in Figure 9l, where intervention should be given in its maximum effort from the beginning of simulation for 10 months and decreasing slowly as time passes. The total cost for Scenario 2 is 1.33452 × 10 12 , which is almost 70% smaller than Scenario 1. The average number of total infected inverted in Scenario 2 is 160,114 individuals. Compared to Scenario 1, the average number of averted cases for Scenario 2 is almost 361% larger.   Figure 10 shows the dynamics of each subpopulation in system (1), and it controls trajectories for Scenario 3. It reflects the success of the combination between the preerythrocytic vaccine and transmission-blocking treatment. This scenario suggests the intervention of transmission-blocking should be given in its maximum effort, much longer than the pre-erythrocytic vaccine. With the total cost for Scenario 3 being 1.334401 × 10 12 , the average number of avoided infections is 170,680 individuals. The comparison of all scenarios in one figure can be seen in Figure 11. We can see that reducing the number of infected individuals between Scenarios 2 and 3 is almost similar, which is more successful than Scenario 1. Hence, with the set of parameters and initial conditions in this article, we can conclude that transmission-blocking treatment is more significant in reducing the number of infected individuals than the pre-erythrocytic vaccine. We have to state that the estimation of parameters on our model does not include the number of individuals who already use both forms of intervention. Hence, it is important to do another parameter estimation when the data for the number of individuals who use both interventions are already available and re-simulate our optimal control simulation.  Figure 11. Dynamics of total susceptible human (a) total of infected human (b) for four different strategies (no intervention (red), u 1 and u 2 (blue), u 1 only (black), and u 2 only (magenta)).

Different Initial Condition of Population
In this subsection, we analyze the impact of the initial condition of the population on the dynamics of control variables. We use parameter values and initial condition of Scenario 3 in Section 6.1 as the base case, and change the initial conditions as follows: 150%N, and W(0) = 50%N, with N being the total human population in Papua (3,435,430). Let us call this Scenario 4. In Scenario 4, we can see that the number of infected individuals in the human and mosquito population is larger than in Scenario 3. The result of this scenario is shown in Figure 12. We can see in Figure 12k that the dynamics of u 1 are only given in its maximum effort for just a few days of simulation, and start to decrease rapidly compared to in Scenario 3 (Please see Figure 10k). Since the rate of vaccine is the maximum given in a concise time period, then the increased number of S 2 in Scenario 4 is not as significant as in Scenario 3. The dynamics of u 2 , which present the transmission-blocking treatment, should be given its maximum value over a longer period than u 1 to reduce the number of infected individuals. Consequently, the number of treated individuals increases very rapidly and decreases as they begin to recover. The cost for Scenario 4 is 4.6189 × 10 12 , which is much larger than Scenario 3. Although the number of infected cases avoided in Scenario 4 is 1,177,307, which is more significant than in Scenario 3, the cost to achieve this result is almost 330% larger than in Scenario 4. Hence, we can conclude that when the number of infected individuals is relatively large at the beginning of the intervention period, most of the intervention should be focused on reducing cases, which in our model is the implementation of transmission-blocking treatment. With more infected individuals at the beginning of the intervention period, more costs are needed to control the spread of malaria.

Different Initial Basic Reproduction Number
In this subsection, the optimal control simulation aims to see the effect of the initial value of the basic reproduction number when pre-erythrocytic intervention and transmissionblocking treatment is not yet implemented. In this case, the basic reproduction number of system (1) is given by (12).
For this simulation, we used parameter values as from Table 1, except that u 1 = 0, u 2 = 0, and also other parameters that related to the control parameters are zero as the base case (case 3 in Table 5). We run our simulation in seven different cases, and the differences are based on the value of β h1 , β h2 , and β v . Hence, we also have a different initial basic reproduction number R 0 , as shown in Table 5. It can be seen that the greater the value of R 0 , the higher the costs to reduce the high number of cases in the field. If we pay attention to Figure 13k,l, the higher the R 0 value, the longer the transmission-blocking treatment intervention should be given to its maximum value, and at the same time, the pre-erythrocytic intervention should be reduced in intensity to avoid high intervention costs. On the other hand, when R 0 is smaller, the pre-erythrocytic vaccine should be given for a longer period at its maximum level to avoid an increase in the number of malaria-infected individuals.  Figure 13. Time series dynamics for; humans (a-h), mosquitoes (i,j), and controls (k,l), but different initial R 0 when no control applied. The color in the legend are explained in Table 5.

Conclusions
In this paper, we formulated and analyzed a deterministic compartmental model for malaria transmission. The model consists of eight compartments for humans and two compartments for mosquitoes. The malaria-free equilibrium is locally asymptotically stable if the basic reproduction number is less than unity. The basic reproduction number of our model presents the combination of infection paths in our transmission model, such as the bite of susceptible mosquitoes to the exposed human (with/without the vaccine), infected, treated, or carrier recovered individuals, or the bite of infected mosquitoes to susceptible humans (with/without vaccine).
Our model parameter values were estimated using incidence data from the Papua and West Papua provinces in Indonesia. We found that without implementing any type of intervention, the basic reproduction number of malaria in Papua and West Papua is greater than one, which suggests that a massive intervention should be made to reduce the spread of malaria. Our sensitivity analysis of the basic reproduction number shows that transmission-blocking is more sensitive to reducing the basic reproduction number compared to the pre-erythrocytic vaccine.
Several scenarios on the optimal control simulation were carried out based on the combination of interventions (first scenario), different initial conditions (second scenario), and the value of the basic reproduction number (third scenario). From the first simulation, in order to understand the most effective combination strategies, we found whether using both vaccine and treatment as a single or combined form of intervention will be more effective compared to other possible combinations. Our analysis shows that using both interventions is the most successful strategy in reducing the number of new infections. However, the number of averted infections with this strategy is only slightly different compared to that of implementing transmission-blocking treatment as a single form of intervention, but with a cheaper cost of implementation. On the other hand, the preerythrocytic vaccine is not recommended as a single form of intervention, as the result in reducing the number of infected individuals is not as significant as other strategies. From the second scenario, we concluded that more cost is needed to control malaria when the number of infected individuals is already high. For the last scenario, we run our simulation for various values of the basic reproduction number to indicate that the level of endemicity of malaria when vaccines and treatment are not yet implemented. Our results suggest the larger the basic reproduction number of an area, the larger the scale of the implementation of transmission-blocking treatment is needed to reduce the spread of malaria. Additionally, the pre-erythrocytic vaccine should be implemented at its maximum rate for a short period to minimize the intervention cost.

Data Availability Statement:
The data that been used in this research is from personal request to Papua and West Papua Provincial Health Office.