Modeling, analysis and numerical solution to malaria fractional model with temporary immunity and relapse

The present paper deals with a fractional-order mathematical epidemic model of malaria transmission accompanied by temporary immunity and relapse. The model is revised by using Caputo fractional operator for the index of memory. We also recommend the utilization of temporary immunity and the possibility of relapse. The theory of locally bounded and Lipschitz is employed to inspect the existence and uniqueness of the solution of the malaria model. It is shown that temporary immunity has a great effect on the dynamical transmission of host and vector populations. The stability analysis of these equilibrium points for fractional-order derivative α and basic reproduction numberR0 is discussed. The model will exhibit a Hopf-type bifurcation. The two control variables are introduced in this model to decrease the number of populations. Mandatory conditions for the control problem are produced. Two types of numerical method via Laplace Adomian decomposition and Runge–Kutta of fourth order for simulating the proposed model with fractional-order derivative are presented. To validate the mathematical results, numerical simulations, sensitivity analysis, convergence analysis, and other important studies are given. The paper is finished with some conclusions and discussion.


Introduction
Malaria is a vector-borne disease that affects the developing countries of the world. It is a life-threatening illness and is caused by a female mosquito plasmodium parasite. As per World Health Organization (WHO) reports that malaria is deadly to host life and remains a dangerous infectious epidemic disease [46]. This endemic disease is an economic load on the countries. From time to time various awareness camps at the international level have been organized by WHO about these vital diseases. Besides, several national and international conferences, seminars, workshops, etc. have also been arranged by different universities to find out the possible results for the management of this endemic disease. Several researchers from epidemiological backgrounds have mathematically researched the propagation mechanisms of the disease to understand and capture the complex relationship between a susceptible host and populations of infected vectors. Mathematical modeling of infectious diseases has also played a crucial role in understanding the complex mechanisms of control strategies and the different transmission parameters of diseases. Ross gave the first epidemiological model of malaria [37]. To demonstrate malaria infection dynamics between vector and host populations, he used the deterministic compartmental epidemic models. A review, as well as a full survey on the modeling of malaria infection, is considered [23]. Since bifurcation is a key rule in computational models of viral diseases and is studied by a lot of researchers prepared a structure model to show the effects of relapse rate in malaria disease [2-4, 15, 19, 26, 39, 40]. The fractional-order derivative is considered as the memory index. The classical epidemic model does not give any information about the learning mechanism of the host population or memory which shows the disease transmission. A periodic mathematical epidemic model has been considered by Ross-MacDonald [16]. They focus on the effects of spatial and temporal heterogeneity on illness's dynamical transmission. Thus the given mathematical epidemic model incorporated periodic variation in vector seasonal and ecology host shifting to catch changes of malaria disease spread among different areas.
The additive compound matrices approach of the endemic equilibrium to show the global stability analysis. Their numerical simulations show that the endemic equilibrium approaches the disease-free equilibrium points. A deterministic mathematical compartmental model to investigate the effect of drug resistance in malaria infection and some authors formulated an optimal control model for malaria and cholera co-infection [29,31]. We have observed the effectiveness of drugs in a dynamical transmission of malaria model [1]. A model to study infectious disease co-infection transmission of malaria and HIV [7]. Furthermore one studied the spread of malaria infection with an optimal control model [28,43,44]. Thus it is formulated as a co-infection model for meningitis and malaria among children [20]. They observed that in the global sense when the threshold quantity R 0 is less than unity, the disease-free equilibrium might show some stability. Mathematical modeling of infectious disease's dynamical process and other main areas of studies such as finance, engineering, and economics has extensively been explored using the theory and applications of the classical differential equations.
But in modern times, the theory and applications of fractional calculus have become extremely important and beneficial in modeling biological processes and another field of studies due to the memory property of fractional derivatives. Many authors have contributed significantly to compartmental mathematical modeling of infectious disease dynamics using fractional differential equations. There has been formulated a vectortransmitted disease model with fractional differential equations. [38]. The fractionalorder mathematical models for malaria disease infection were studied [33]. Further, there has been analyzed the fractional-order co-infection models for TB and HIV [32]. Also, the fractional-order models for HIV infection have been studied and analyzed by many researchers who formulated a fractional-order competition model for the love triangle [5,12,16,21,22]. One also studied the backward bifurcation in a fractional-order model with vaccination parameter [13].
A fractional-order model for the cholera infection was formulated [18]. Moreover, the homotopy analysis method and Runge-Kutta method of fourth order to evaluate the numerical analytical and results for fractional-order childhood disease models were discussed [5]. Therefore the epidemic model with fractional order for influenza has been analyzed and formulated and one studied the fractional-order SIRC model with Salmonella bacteria infection [14,17,36]. The deadly Ebola disease which killed a lot of people in some parts of Africa has been modeled by a system of fractional-order derivatives [6]. An endemic model with a constant population has been numerically and analytically studied by using Caputo fractional derivatives [30].
The rest of the paper is organized as follows. Some basic preliminaries on fractional calculus are presented in Sect. 2. Section 3 is devoted to the model description. The basic properties of the proposed model are given in Sect. 4. The analysis is discussed in Sect. 5. Optimality is presented in Sect. 6. The numerical solutions are obtained in Sect. 7. The numerical simulations, sensitivity analysis, and convergence analysis are discussed in Sect. 8. A discussion is in Sect. 9 and finally, the conclusion is drawn in Sect. 10.

Preliminaries on fractional calculus
Some basic definitions of fractional calculus areas are given in the following.

Definition 2.1
The Riemann-Liouville fractional integral of the function f : R + → R exists for order α > 0 in two forms, forward and backward. For the closed interval [u, v], the two integral are defined as [34] RL where is the gamma function.

Definition 2.2
The Riemann-Liouville fractional derivative of the function f : R + → R have also in two forms, forward and backward. This α is solved by the use of Lagrange's rule. To find the nth-order derivative over the integral of (nα), the derivative is obtained. If n > α, where n is the smallest integer, then the derivatives are

Definition 2.3
Necessitated by the drawback of the Riemann-Liouville fractional derivative a definition was given [8]: where α ∈ (n -1, n), and n ∈ N.
Obviously, C 0 D α t f (t) approaches towards D α t f (t) whenever the fractional order tends to one. C 0 D α t f (t) and C 0 D α t g(t) exist almost everywhere [11] and letting r a , r b ∈ R, then  [45], C 0 D α t e * (t) = f t, e * (t) ⇔ f t, e * t = 0, where 0 < α < 1.

Model description
In this section, we formulate a deterministic mathematical model of malaria epidemic model by dividing total population into two mutually exclusive classes namely, the vector (mosquito) and host (human). To make the dynamical transmission of malaria model with temporary immunity β and relapse H , two classes of mosquito population are considered, the mosquito susceptible class, at time t, is denoted by M S (t); mosquito infectious class M I (t). Hence the total mosquito population 1 is given by On the other hand, three classes of host(human) population; human susceptible H S (t), human infectious H I (t) and human recovered H R (t) are considered so that the total host population 2 is The notations and parametric values of the state variables used in the formulation of dynamical transmission of model are given in Tables 1 and 2. Let 1 be the recruitment rate of the vector population, 2 be the recruitment rate of the host population, δ M be the rate of infection in mosquitoes population, δ H be the rate of infection in humans population, d M be the natural death rate of mosquito population,  Induced death rate in humans population 0.1 [25] β The temporary immunity in hosts population 0.05-1 Estimated Figure 1 Schematic diagram of malaria with temporary immunity d H be the natural death rate of the human population, φ H be the vaccination rate in the host population, H be the relapse rate of malaria in the host population, and θ H be the induced death rate in the host population. The schematic diagram of malaria with temporary immunity is portrayed in Fig. 1.
The flows from the susceptible to infected classes of mosquitoes and humans populations depend on the transmission probabilities δ M , δ H and the number of infectious and susceptibles of each species. In M S class, the recruitment rate of mosquitoes population is 1 . δ M is the rate of infection between M S to H I . Similarly in H S class the recruitment rate of the human population is 2 . δ H is the rate of infection between H S to M I . β is the rate of flow between H R and H S , φ H is the rate of flow between H I and H R and θ H is the rate for H I class, respectively.

Classical integer model
The classical integer model of epidemic malaria between mosquito-to-human and vice versa is as follows:

Fractional-order mathematical model
In this subsection, we have modified the classical integer system (3.1) into fractional-order in sense of Caputo. The system of the non-linear fractional differential equations is given as

Humans population
The system (3.3) gives the dynamics of host populations, and all state variables and parameters are assumed to be non-negative.

Basic properties of the model
In this section for the well-posedness of the system (3.3), its existence and uniqueness criteria are given below. Let x = M S , M I , H S , H I , H R , and D is either a Riemann-Liouville or a Caputo fractional operator, then with initial conditions are x(t 0 ). Let α ∈ (0, 1), and g is piecewise continuous in time t and locally Lipschitz in x on [t 0 , ∞] × V to R n , where V ∈ R n is a domain.

Existence and uniqueness
The fact that g(t, x) is locally bounded and Lipschitz in x implies the existence and uniqueness of the solution of Eq. (4.1). Thus we have the following.
x) be the real-valued continuous function, α ≥ 0, and · be an arbitrary norm, Proof From Definition 2.1, and the properties of the norm function we have Proof Applying the fractional integral operator t 0 I α t on both sides of Eq. (4.1), thus it follows from Lemma 4.1 and the Lipschitz condition that There exists M(t) > 0, we have Using Laplace transform on both sides of Eq. (4.2), we have Using the inverse Laplace transform on both sides of Eq. (4.3), we have Hence the result is proved.

Invariant region and attractivity
The dynamical transmission of the system (3.3) will be analyzed in the biologically feasible region as,

Lemma 4.3 The biological feasible region
Proof Adding the first two equations of the system (3.3), the total vectors population, is The solution of the fractional system (3.3) is where E α,β is the Mittag-Leffler function. The behavior of this function is asymptotic, so it gives Therefore, from Eq. (4.6) 1 tends towards 1 /d M as time(t) tends to infinity. Further, the proof for the case of the host individuals is completely similar to the vector and hence it is omitted. Thus, for all positive values of t, all the values of the Caputo fractional derivative the system (3.3) with i.c.s of Z remain in Z. Therefore, the region Z is positively invariant for the system (3.3) and attracts all solutions in R 5 + .
To show the system (3.3) has non-negative solution, we have R 5

Positivity and boundedness
Proposition 4.5 The solution of the system (3.3) is positive, bounded for all x(0) ∈ R 5 + , for t > 0.
Proof To demonstrate the positivity, it is required that on every hyperplane bounding of R 5 + , from the system (3.3), we have Thus, by Corollary 4.4, the above target set has been achieved, that is, the solution will stay in R 5 + and thus we have the following biologically feasible region: Therefore all the terms of the sum are positive, then the solution of the system (3.3) is bounded.

Lemma 4.6
The feasible region Z is positively invariant for the system (3.3) with an initial condition in R 5 + . It is sufficient to show the dynamics of the aforesaid system in the region given in Z. So this region can be supposed epidemiologically and biologically.

The analysis
It is hard to find the stability of the equilibrium from the equations of the system (3.3) because the fractional derivatives do not follow the Leibniz rule [41]. So, we use the following transformation: (5.1) The modified system (5.2) is given as Thus the malaria system (3.3) is equivalent to the system (5.2), so it is enough to know more about the stability properties of the system (5.2).

Disease-free equilibrium (DFE)
The DFE point of the system (3.3) is To know more about the stability of the equilibrium point E 0 , we build the Jacobian matrix J at E 0 of the system (5.2) as Since the system (5.2) is a fractional-order compartmental model with multi-orders except for α = 2 3 . So, the multi-orders can be mathematically among any numbers like rational, irrational, real, or complex. But in the literature, there is no existing method for checking the stability of the fractional derivatives compartmental model with multi-orders among irrational, real, or complex numbers. Therefore, we assume that α is a rational number [9]. The stability of E 0 is defined by the solutions of the following characteristic equation of J E 0 : where D and denotes the determinant and the diagonal matrix. But the stability of E 0 can be investigated by the use of the following result [42].
Remark 1 For all E 0 of the system (5.2), the λ are found from the following characteristic equation: where J(E 0 ) is the jacobian matrix evaluated at the DFE points and F is the LCM of the denominators of the rational numbers 5(p 1 , p 2 ), and p 1 , p 2 ∈ Q.
Thus, the stability in system (5.2) of the remaining factor depends on the nature of the values of the characteristic equation (5.5), and we have Hence, if all the values λ e of Eq. (5.6) yield | arg(λ e )| > π 2S , where, i = 0, 1, 2, . . . , 2S, then E 0 of the system (5.2) is locally asymptotically stable.
Therefore, the stability of E 0 of the system (5.2) of the remaining term depends on the nature of the solutions in the characteristic equation (5.5).
(iii) Let us discuss the mandatory condition for the existence of Hopf-type bifurcation for the system (5.2). The condition for the Hopf-type bifurcation of the system (5.2) is that if for some roots λ q and α * = R S ∈ (0, 1) of Eq. (5.6) then | arg(λ q )| = π 2S . At α * the system (5.2) shows a Hopf-type bifurcation [24]. Now, we will study about the sufficient conditions of a Hopf-type bifurcation for the existence of the system (5.2), let us suppose ∃ a value λ of Eq. (5.6) in the form λ = re iπ 2S , and r ∈ R + and i = √ -1. Substitute λ in Eq. (5.6), we have the two equations in the variable r: The above equations can be changed into the variable r as below: By Descartes's rules, the polynomial equation (5.8) has either zero or two positive root. If Eq. (5.8) having two positive solution, which are conjugate root of the form λ = re ± iπ 2S of Eq. (5.6) for some α * = R S ∈ (0, 1). Hence the system (5.2) shows a Hopf-type bifurcation at the order α * .

Remark 2 The variation of α gives
It is important to note here that we can get back the result of the system (3.1) after putting the value of α = 1. So, on the base of this fact at the fractionalorder derivative α = 1 of the system, Eq. (5.6) is of the form is known as a threshold quantity, and this is dimensionless and this quantity is epidemiologically known as the basic reproduction number. The basic reproduction number is simply called the number of secondary infections caused by the primary substance or the potential which measures the spread of the disease in a population.
It is concluded that if α = 1 we can go backward towards the classical result of the epidemiological view for our system (3.1). If R 0 < 1, the E 0 is locally asymptotically stable and if R 0 > 1, the E 0 is unstable and the disease persists in the population [10].
Remark 3 If α = 2 3 , then the system (5.2) converts into a one order fractional compartmental model. So from this fact we can conclude that our main stability results in Theorem 5.2 in terms of the threshold quantity (R 0 ).
We have the following result.

Proposition 5.3 When the order α = 2 3 , then the equilibrium value of threshold quantity
2) is unstable and the system (5.2) shows a Hopf-type bifurcation; and (iii) if R 0 = R c 0 , then the DFE E 0 is stable.

Positive equilibrium
The system (3.3) has a positive equilibrium point We address the locally asymptotically stable of the positive (endemic) equilibrium point of the system (3.3). The Jacobian matrix J(E 1 ) is given as The characteristic equation of the Jacobian matrix at E 1 is where  We have the following proposition.

Proposition 5.4
If Dis(f ) > 0 and Routh-Hurwitz criteria are satisfied, that is, , then E 1 is locally asymptotically stable.

Optimal control
We can put the control parameter on this serious problem to prevent its spreading. The following are the control parameters: h 1 : Infectious mosquitoes, which should be killed. h 2 : Malaria patients, who should be treated. Now, the objective function is where is the set of all compartmental variables, W 1 , W 2 , W 3 , W 4 , W 5 of the positive weight constants for the variables M S , M I , H S , H I , H R . Now, we will find every value of control variables from t = 0 to T f s.t., where M is the smooth function for the interval [0, 1]. Therefore, the Langrangian function related to the objective function F is given by The adjoint variables, ζ i = (ζ 1 , ζ 2 , ζ 3 , ζ 4 , ζ 5 ), for the system are calculated by taking the partial derivatives of ℵ with respect to each variable, Hence, this calculation gives H I (ζ 4ζ 5 ) 2u 2 . (6.4)

Numerical methods
In this section, we discuss the approximate solution via numerical methods such as Runge-Kutta method of fourth order (RK4) and Laplace Adomian decomposition method (LADM) for the system (3.3). For the sake of easiness we omit the RK4. But graphically we will observe it. We apply the Laplace transform on both sides of system (3.3), which yields After some simplification, using the initial conditions n 1 , n 2 , n 3 , n 4 , n 5 and taking the inverse Laplace to transform to system (7.1), we have 2) Suppose that the solutions x are in the form of infinite series that are given by  (7.6) . . . The Laplace Adomian decomposition method and Runge-Kutta method of fourth order both are used to find the approximate solution of the non-linear fractional-order system. We use both these methods to check the variations for biological state variables for fractional-order derivative α and parameters. The Laplace Adomian decomposition method is more efficient and accurate than the Runge-Kutta method. From Figs. 2(a&b) and 3(a&b), we have seen that the system (3.3) has more degree of freedom than system (3.1). The dynamical behavior of various compartments has been shown. From the graphics of the Laplace Adomian decomposition method, it is clear that the result obtained by using the Laplace domain decomposition method is very efficient and dramatically increased by the increase of terms. One more interesting point to be looked at is that for the use of the small interval of time, we have assumed comparatively small initial values. For the large interval of time, the initial value should be taken large so that the concerned population of vectors and hosts may not be negative and vice versa. Now we will make the graphs through the Runge-Kutta method from Figs. 4(a&b), 5(a&b), 6(a&b), 7(a&b) and 8(a&b). It is observed from Fig. 4(a) that when δ M increases, M S decreases. It means as the rate of infection in vector increases the number of susceptible mosquitoes naturally decreases. Next, it is observed from Fig. 4(b) that when δ M increases, M I increases. It means as the rate of infection in vectors increases the number of infectious mosquitoes naturally goes on increases. This is a natural phenomenon and we observe it correctly through graphs. Similarly, it is observed from Fig. 5(a) that when δ H decreases, H S increases. It means as the rate of infection in the host decreases the number of susceptible humans naturally goes on the increase. Next, it is observed from Fig. 5(b) that when δ H decreases, H I decreases. It means as the rate of infection in the host decreases the number of infectious humans naturally goes on decreases. Also, it is observed from Fig. 6(a) that when the rate of temporary immunity in human population β increases the H S class goes on increases, and from Fig. 6(b) it is observed that as β decrease H R decreases. Next from Fig. 7(a) it is observed that when the rate of vaccination φ H increases in the human population the number of infectious humans H I decreases and H R increases and this thing is seen in Fig. 7(b). Further, it is observed from Fig. 8(a&b) that as H increases H I increases and H R decreases. All these graphs can be observed naturally. Hence by the Runge-Kutta method, we have checked all the possible variations.   Fig. 9(a&b) that the biological parameter φ H which is the vaccination rate in the host population has a higher influence on R 0 than the parameters proportional to the transmission probabilities δ M , δ H . From Fig. 10(a&b) it is seen that the biological parameter φ H which is the vaccination rate in the host population has higher influence on R 0 than the parameter proportional to the transmission probabilities H and the biological parameter δ M , that is, the rate of infection in vectors population has higher influence on R 0 than the parameter proportional to the transmission probabilities H . Moreover from Fig. 11(a&b) one sees that the biological parameter δ H , that is, the rate of infection in the host population has higher influence on R 0 than the parameter proportional to the transmission probabilities H and the biological parameter φ H which is the vaccination rate in hosts population has higher influence on R 0 than the parameter proportional to the transmission probabilities δ M . Thus from the sensitivity analysis, we can observe that controlling the sexual route of dynamical transmission of this malaria disease is very important if we plan to control the transmission of this illness. The numerical simulation is performed to hold up our analytical discoveries. The values of the biological parameters are depicted in Table 2. It is noted that recruitment rate parameters are in per year, and all other biological parameters are in per day. To show the stability of disease-free equilibrium point, we consider the biological parametric values of different parameters that are involved in the basic reproduction number are listed as 1 Table 2. For these parameters R 0 = 0.679 and the disease-free equilibrium point E 0 = (1.66, 0, 200, 0, 0). The behavior of these points can be observed.

Convergence analysis
The solution in Eq. (7.8) is a series and is rapidly convergent as well as converges uniformly to the exact value. To check the convergence of the series (7.8), we will use the classical technique [35]. For the sufficient conditions of convergence of this method, we have the following theorem. Theorem 8.1 Let K be a Banach space and F : K → K be a non-linear contractive operator such that for all x, x ∈ K , F(x) -F(x ) ≤ u xx , 0 < u < 1. Then F has a unique point x such that Fx = x. By the use of Adomian decomposition method the series given in Eq. (7.8) can be written as , t = 1, 2, 3, . . . , and suppose that x 0 ∈ B r (x), where B r (x) = {x ∈ K : xx < r}, then we have (a) x t ∈ B r (x); (b) lim t→∞ x t = x.
Proof For (a), we use mathematical induction for t = 1, to get Suppose that this result is true for s -1, then that is, x tx ≤ u t x 0x ≤ u t r < r therefore, it implies that x t ∈ B r (x).

Discussion
In this paper, we analyzed numerically the mathematical malaria model. The existence and uniqueness of the malaria model are discussed. We have observed that the threshold quantity is similar to the basic reproduction number and is denoted by R 0 . We also have studied the possible equilibrium points and basic reproduction number of this disease. Next, we have studied the stability analysis of these equilibrium points to the fractionalorder derivative α and R 0 . We have determined the fractional-order derivative dependent threshold values for R 0 , below which E 0 is always stable, above which E 0 is unstable, and also at which the model shows a Hopf-type bifurcation. We analytically verify our results when the fractional-order derivative α = 2 3 and this is the special result in our paper. The numerical analysis work has been done by the use of the Laplace Adomian decomposition and Runge-Kutta method to compute an approximate solution of the malaria epidemic model. We also discuss the optimality of two biological control parameters, and present a sensitivity analysis of four parameters behavior on R 0 . We also have discussed the solutions of fractional differential equations. in the form of infinite series and their convergence analysis.

Conclusion
The fractional-order mathematical model of malaria transmission with temporary immunity and relapse via the Laplace Adomian decomposition method is presented. This numerical method is a fruitful tool to solve the non-linear models that are widely used in applied mathematics and engineering. Further, we have given a convergence result for this numerical method. Also, the aforesaid methods give good numerical results for the nonlinear fractional-order model as compared to the other numerical methods like homotopy analysis and homotopy perturbation. These two methods have an extra parameter namely h on which the solutions depend. But our proposed numerical methods do not need any kind of parameter, so it is easy to understand and then implement. We observe that the solution found through the LADM closely agrees with those obtained by other methods, like RK4.