Optimal control analysis of vector-host model with saturated treatment.

Vector-host infectious diseases remain a challenging issue and cause millions of deaths each year globally. In such outbreaks, many countries especially developing or underdevelopment faces a situation where the number of infected individuals is getting larger and the medical facilities are limited. In this paper, we construct an epidemic model to explore the transmission dynamics of vector-borne diseases with nonlinear saturated incidence rate and saturated treatment function. This type of incidence rate, as well as the saturated treatment function, is also known as the Holling type II form and describes the effect of delayed treatment. Initially, we formulate a mathematical model and then present the basic analysis of the model including the positivity and boundedness of the solution. The threshold quantity R 0 is presented and the stability analysis of the system is carried out for the model equilibria. The global stability results are shown using the Lyapunov function of Goh-Voltera type. The existence of backward bifurcation is discussed using the central manifold theory. Further, the global sensitivity analysis of the model is carried out using the Latin Hypercube sampling and the partial rank correlation coefficient techniques. Moreover, an optimal control problem is formulated and the necessary optimality conditions are investigated in order to eradicate the disease in a community. Four strategies are presented by choosing different set of controls combination for the disease minimization. Finally, the numerical simulations of each strategy are depicted to demonstrate the importance of suggesting control interventions on the disease dynamics and eradication. © Società Italiana di Fisica and Springer-Verlag GmbH Germany, part of Springer Nature 2020.

Introduction

Vector-host diseases are regarding as a major part of infectious diseases caused by pathogens and transmitted by insects, ticks, bacteria and protozoa. The disease transmitting biological agents are called vectors. The vectors carry the infection from one human to another or from animals to humans. Vector-host diseases are considered to be 17% of all infectious diseases which results about one million deaths each year. Mostly, the outbreaks of these infectious diseases are reported in the tropical and sub-tropical regions and areas where the sanitation systems are problematic and safe drinking-water is unavailable [1]. Currently, about 14 vector-borne diseases are reported in the USA that are of national public health concern. Dengue fever is one of the rapidly spreading worldwide mosquito-borne viral infections [2, 3]. It is a disease of huge complexity and becomes a major health problem in recent years. Due to the rapid growth in the dengue incidence in tropical and sub-tropical countries, especially in urban and semi-urban regions, this infectious disease is posing a serious threat to the world’s population. In some areas of Asia and Latin America, dengue is a major cause of illness and mortality among children. A recent study indicates that there are 390 million dengue infection cases each year [2]. This viral infection is reported in only nine countries before 1970. Now the dengue has become a global problem and cases have been reported in more than 125 countries with almost 400 million infections globally every year [1]. In 2015, the worldwide outbreak of dengue was increased dramatically. In the USA alone, 2.35 million infective cases were recorded. In Philippine, about 169,000 new cases were reported and Malaysia exceeded 111,000 suspected cases of dengue. In Brazil, over 1.5 million infected individuals were reported in 2015, almost 3 times higher than in 2014. Also in India, the recorded cases were over 15,000 in 2015. More details of the global burden of dengue can be found in [2]. Malaria is another most prevalent vector-borne infection transmitted by the Anopheles mosquitoes. In 2015, more than 3.2 billion people were at risk, most of them, children under the age of five years. Malaria transmission was found in 95 countries and territories of the world. Globally 100 million people are affected by diseases such as Chagas disease, leishmaniasis, and schistosomiasis [1]. Most of these diseases are controllable by adapting protective actions. The Zika virus infection is another deadly mosquito transmitted disease. This virus is in fact carried by the same vector that transmits dengue virus and some other viral infections. In October 2015, an outbreak of Zika epidemic has been recorded initiated from Brazil and has spread to most parts of South and Central America [4]. Currently, the mosquito-transmitted Zika infection is reported in a total of 86 countries and territories around the globe.

The mathematical modeling approach gains much interest in recent years and is helpful to explore the complex dynamics of any real world problem including infectious diseases. These models are used significantly to set an appropriate controlling strategy to eradicate the disease in a community [5-8]. Using numerical simulation from these models one can predict the disease transmission, i.e., a disease incidences under certain assumptions, and a threshold for epidemic outbreaks can also be obtained from the results of these simulations. An epidemic model based on non-integer order differential operators with singular or non-singular kernel provides a better understanding o the disease dynamics. A number of mathematical models formulated for the transmission dynamics of various infectious diseases are available in the literature for example see [9-15] and references therein. In the last few decades, numerous epidemic models were studied by researchers to analyze and predict the transmission patterns of various vector-borne disease outbreaks in different regions and countries. Anderson et al. [16] formulate an epidemic model considering the geographic spread of dengue fever. A new transmission model was formulated and studied in [3] to explore the dynamics of dengue fever in mainland China. Recently, new mathematical models of Zika virus infection in various countries were formulated in [4, 17]. A combined study of dengue and chikungunya dynamical models have been formulated and analyzed in [18]. There are many mathematical models available in the literature that characterized the dynamics of different vector-borne infections [19-24]. The implementation of mathematical modeling approach to the transmission patterns of the novel COVID-19 pandemic is recently studied in [25, 26].

In the literature, various mathematical models of infectious diseases with bilinear incidence rate have been suggested where S is suspectable and I are infected population [19, 20]. In [27], the Capasso and Serio for the first time introduced the incidence rate in the form , where measures the force of infection in a completely susceptible individuals and accounts the inhibition effects due to the behavioral change of susceptible people when their number increases or due to the crowding effect of the infected population. The saturated incidence rate is more reasonable than the bilinear incidence rate and is used in many models like [21, 22, 24, 28]. Further, the saturated treatment function is biologically more feasible to use in a situation of sudden outbreaks of a highly infectious disease. In such cases, the population lacks enough medical facilities to overcome the situation, some infected individuals are delayed in getting treatment because far distances to be covered and the availability of sufficient materials therein such as beds, isolation centers e.t.c.

Motivated by the above discussion, in the current study, we present a new transmission model in order to explore the dynamics and possible control of the vector-host epidemics. The saturated incidence and treatment rates are utilized to formulate the proposed model. Initially, we construct the proposed compartmental model with no control measures by considering four different classes for humans and three compartments for vector (mosquitoes) population. Then we carry out the global sensitivity analysis and reformulate the model with three appropriate control measures. The remaining sections are organized in the following pattern: The model description is discussed briefly in Sect. 2. The basic analysis of the model including stability results as well as the existence of backward bifurcation is investigated in Sect. 3. Section 4 contains the global sensitivity analysis. Characterization of the optimal control is presented in Sect. 5. The impact of various control strategies is depicted graphically in Sect. 6. Finally, a brief concluding remarks of the present investigations are written in Sect. 7.

Model formulation

In this section, we briefly present the construction of the proposed mathematical model that divides the host population denoted by and vector population denoted by in different mutually exclusive subclasses. The host population is divided into four distinct subclasses which are susceptible , exposed , infectious and denotes the recovered humans at any time instant t so that,The total population of vectors is categorized further into three subclasses, namely, the susceptible, , exposed , and infectious vectors, at any time instant t, soThe recovery in the vector population does not exist, because the infection period of mosquitoes ends with their death. The susceptible humans are generated by the constant recruitment rate . It is decreased by the effective contact rate , where denotes the disease transmission rate from infected vectors to susceptible humans and the parameter shows saturation level of the force of infection. The class of susceptible individuals is further reduced due to the natural mortality at the rate . The susceptible individuals acquired infection after effective interaction with infected vectors and they first enter to the exposed human class. This class is reduced by the natural mortality rate and becomes infectious at the rate after the completion of incubation period. The population of infected humans is generated by the infection of exposed individuals at the rate of and decreased by the natural death rate , the disease-related death rate and become recovered at the rate , where b denotes delay in treatment. The recovered human population is increased by joining the infected individuals after successful treatment while decreased by the natural death rate .

The class of susceptible vectors is generated by the recruitment and is decreased by moving to exposed class after the effective contact with infected humans at the rate , where is the disease transmission rate from human to vector and is the level of saturation. This class is further decreased by the natural death rate . The population of exposed mosquitoes is increased by joining of susceptible vectors after effective contact with infected humans and decreased by the natural death rate , and become infectious at the rate . The exposed vectors join the infected class and become infectious after completion of the incubation period. This vector class is also decreased by the natural death rate . Thus, keeping the above discussion in mind, the proposed vector-host epidemic model can be formulated mathematically by the following system of differential equations:subject to the initial conditionsThe detail biological description with numerical values is given in Table 2 and the transmission diagram of both host and vector population is depicted in Fig. 1.

Table 2

Biological description and numerical values for the model (1) parameters

Parameter	Description	Value (day^-1)	References
Λ_H	Recruitment rate of hosts	0.002	Assumed
Λ_V	Recruitment rate of vectors	0.4	Assumed
β_1	Host contact rate	0.00375	[33]
β_2	Vectors contact rate	0.375	[33]
d_H	Natural mortality rate in host population	1/(365*67.7)	[34, 35]
δ_H	Disease-induced death in infected host class	0.001	Assumed
d_V	Natural mortality rate in vectors	1/15	[35]
α_1	Saturation in hosts	0.03	Assumed
α_2	Saturation in vectors	0.03	Assumed
γ	Treatment rate	0.3288330	[36, 37]
ξ_V	Transmission from E_H to I_H class	1/10	[35, 38]
ξ_H	Transmission from E_V to I_V class	0.001	Assumed
b	Delay in treatment	0.1	Assumed

Fig. 1

Flowchart describing the transmission of the vector-host disease model (1)

Model analysis

Here, we present the important and necessary results of the vector-host transmission model formulated in (1). First, we will show that the model posses a positivity as well as boundedness solution.

Positivity and boundedness

The total dynamics of only human’s population is demonstrated byThe total population will be bounded and remain positive for all finite time as the above initial conditions make sure that . Further, the total vector dynamics isWhen then the total dynamics of human and vectors approach as . Hence, the biologically feasible region for the proposed model is

Invariant region

Proposition 1

The closed set is bounded as well as positively invariant. Moreover, it attracts all the model solutions

Proof

To proceed the desire proof, let us considerTaking its time derivative that isFurther, we can prove thatFrom (6), it follows that , which means that is positively invariant set. Further, solving (6) using integrating factor technique, we haveThus as , , which concludes that is an attractor and biologically feasible set.

The basic reproduction number

The vector-host transmission model formulated in (1) possess two types of equilibria, i.e., the infection free or disease free equilibrium (DFE) and the endemic equilibrium (EE). The DFE denoted by and is given asFurther, the threshold quantity also called the basic reproduction number is an important quantity which provides information about the future spread of a disease. It is usually denoted by and accounts that how much an infectious person infect other healthy individuals. In epidemiological perspective, this threshold parameter is stated as: “the average number of secondary infections generated by a typical single infective individual in the entirely susceptible individuals”. We use a well known next-generation method [29] to find . The corresponding matrices for the next-generation method are

Since, , so

Existence of endemic equilibrium

The present section explores the existence of EE of the vector-host epidemic model described in (1). Let us denote the EE bywhere the following expressions are obtained from the solution of system (1) at the steady state:whereSubstituting the values in (7) into the first equation of (1) and after some manipulation, it is observed that the non-zero equilibria obeys the following quadratic equation in :whereClearly, and is positive whenever and negative when . Hence, the following consequences are drawn:

Theorem 1

The model (1) possess

a unique EE if and ;

a unique EE if , or ;

two EE if , and ;

and no EE otherwise.

Local stability analysis

The present section explores the local stability analysis of both the DFE and EE. For the local asymptotical stability (LAS) at the DFE, we state the following theorem.

Theorem 2

If , then the DFE of the system (1) is LAS on and it is unstable otherwise.

The Jacobian of the vector-host transmission model (1) evaluated at is describe as follows:The characteristic equation of is as follows:whereClearly, if , then , for also it can be easily verified through some algebraic manipulations that . Thus, implementing Routh-Hurtwiz criteria the DFE of the system (1) is LAS about the equilibrium point .

In the next result, we analyze the LAS of the epidemic model (1) at EE .

Theorem 3

For 1, the EE is LAS in the region and unstable otherwise.

To proceed the proof, the Jacobian matrix of (1) at the EE is evaluated as follows:One of the eigenvalues is which has clearly negative real part. For the remaining eigenvalues, we define the following characteristic equation of :The coefficients involved in the (10) are as follows:where

Clearly, for Therefore, if and the necessary condition of Routh–Hurwitz criteria for the 6th-order characteristic polynomial in (10) is satisfied, then the EE of the model (1) will be GAS.

Global stability analysis

The present section analyzes the global asymptotic stability (GAS) results of the model (1). For DFE case, we establish the following result.

Theorem 4

If 1, then the DFE given by of the system (1) is LAS on the region and unstable for 1.

In order to prove the desire result, we utilized the Lyapunov function given below:where, , with are unknown constants to be selected according to the demand of the proof. After computing the time derivative of L(t) along the solutions of the system (1), we haveSimplifying we haveNow we choosePutting in (12) and simplifying we getHence, is negative if . Also note that, if and only if , , . Hence utilizing the LaSalle’s invariant principle [30], is GAS in the region .

Next, we present the GAS of the system (1) for the special case when . For this purpose, we derive the following relations at of the model (1):whereThe above equations will be used in the following stability result.

Theorem 5

For the case when , then the EE is GAS if , and condition (25) is hold.

In order to establish the desired proof, we define an appropriate Lyapunov function as follows:Differentiating along with the model (1) solution, we proceed as below:By direct calculations, we have:Using (14) and after some simplification and arrangement, it follows from (17-22)In the expressions (23), it is obvious thatand ifthen . Hence by the asymptotic stability theorem [30], the EE is GAS on .

Analysis of backward bifurcation

The present part investigates the existence of backward bifurcation phenomena in the vector-host model (1). For this purpose, we utilize the concept of the center manifold theory presented in [31]. Taking as the bifurcation parameter and setting , we obtain the following expressionFurther, by replacing , , , , , and and using the vector notation then the model (1) can be written equivalently as , where are as shown belowThe relative Jacobian of (1) evaluated at DFE with isThe matrix described by J as above has a simple zero eigenvalue. Moreover, we denote the corresponding right and left eigenvectors of the zero eigenvalue by and , respectively, where the components are evaluated as follows:andComputation of a: To evaluate a, using formulawe obtained a asComputation of b: To compute b make use of the formulawe obtained b asSince, therefore, the existence of the backward bifurcation is totally dependent on the sign of coefficient a.

Global sensitivity analysis

The aim of sensitivity results of is to study the effect of changes in the most influential parameters on the disease dynamics that greatly affects the spread of the infection. To examine the sensitivity of parameters of the model, we utilize the statistical techniques namely, the LHS coupled with PRCC as developed in [32]. This technique is one of the most effective techniques and can be applied successfully to acquire monotonic nonlinear relationships among the inputs and outcomes of the system. Utilizing this approach, we calculate the PRCC values of the parameters and its corresponding p-values, which can check the level of uncertainty in the dynamical model. The parameters of the model (1) with absolute maximum PRCC values and lowest p-values are known as dominant parameters. Figure 2 shows the graphical results of the PRCC values of the parameters used in the vector-host model. Table 1 has the numerical results of the PRCC and corresponding p-values of the vector-host model parameters. From Table 1 and Fig. 2, we notice that the contact rates and have positive maximum PRCC values and lowest p-values followed by and , respectively. On the other hand, has negative high PRCC and lowest p-value followed by , and , respectively. So, clearly one can observe that , , , and are the most dominant factors of the newly developed vector-host model. Hence, the spread of vector-host infectious disease can be controlled by reducing the contact rates , and enhancing the rates of , parameters.

Fig. 2

Graphical PRCC results for versus various model parameters

Table 1

The numerical PRCC values and the corresponding p values of the parameters of vector-host model (1) relative to

Parameter	PRCC	p values
Λ_H	+0.3172	0.000
Λ_V	+0.3333	0.000
β_1	+0.8023	0.000
β_2	+0.7254	0.000
d_H	−0.1459	0.4152
δ_H	−0.6768	0.000
d_V	−0.9801	0.000
γ	−0.3613	0.000
ξ_H	0.3530	0.000
ξ_V	0.3251	0.000

Optimal control problem

The present section is devoted to formulating a control problem in order to mitigate the infection from a community. Under consideration of sensitivity results carried in the previous section, we introduce three time-dependent control variables in the newly proposed vector-host model to construct an optimal control problem (30). These time-depended control measures are denoted by , and , where represents the prevention control by using the bed-nets or mosquitoes repellent. The term (1-) reduces the contact rates , of the susceptible humans(vectors) with infected vectors(humans), respectively. The second control function is used to enhance the treatment of infected humans (increase ) and the third control function is used for insecticide spray against vectors population. These control measures play a major role to minimize or eradicate the disease from the human population. After incorporating the above-mentioned controls in the vector-host model, we obtain the following control problem:Subject to non-negative initial conditions stated with model (1). The corresponding objective function for the purposed model to minimize the infection is given as follows:The constants , for and in (31) represent the balancing cost factors and denotes the final time level. The objective functional considered in this study is of quadratic type since the intervention is nonlinear. Such type of approach is successfully utilized in the literature and references therein [39-42]. Further, the main goal is to investigate optimal controls , and such thatwith the associated control set given as follows:Furthermore, the Lagrangian of the system (30) describing the control problem is defined in (32) and the Hamiltonian is defined in (33) as follows:andwhere the adjoint variables are described by , for .

Optimal control solution

The present subsection explores the optimal solution of the control model (30) via a well-known Pontryagin’s Maximum Principle discussed in [43]. To proceed further, let , m=1,2,3 denote the corresponding optimal solution of the system (30), then in order to find the adjoint variables the necessary conditions stated below are implemented: The solution of optimal system (30) with the help of conditions demonstrated in (34) is derived in the following theorem.

Theorem 6

For the optimal controls denoted with , , and the respective solutions , , , , and of the control problem (30) that minimize over , there exists adjoint variables , where along with transversality conditions such thatFurthermore, the associated optimal controls for are given by

To obtain the desired expressions along with the transversality conditions stated in (35), we make use of the conditions demonstrated in (34) for the Hamiltonian function described in (33). Moreover, the settings , , , , , , are taken to obtain the desired results. Finally, the condition given in (34) for and are utilized to derive the expressions in equations (36).

Numerical simulations

This section presents the graphical interpretation of different control strategies to explore the impact of each control intervention on infection dynamics and minimization. The control and without control variables systems are solved numerically by utilizing the 4th-order Range-Kutta technique. The values of parameters mostly taken from the literature used in the simulation results are tabulated in Table 2. In Figs. 3, 4, 5 and 6, the solutions of the state variables in the control problem are denoted by blue dashed lines and the red bold lines are used for the solution curves of state variables when there is no control. The weight and balancing constants involved in the objective functional are chosen as , , , , and in order to perform the simulations. The time unit taken is per day. We develop four different strategies by selecting different combinations of the control variables. The graphical impact and discussion of each strategy on the disease dynamics and eradication are as follows:

Fig. 3

Graphical results of the model (30) with , and

Fig. 4

Graphical results of the model (30) with , and

Fig. 5

Graphical results of the model (30) with , and

Fig. 6

Graphical results of the model (30) with , ,

Strategy 1: treatment and insecticides spray controls ()

In the first strategy, we simulate the control model (30) except the prevention control , which means that people are not aware to use the bed-nets or other repellents to reduce the effective contacts. The combined effect of and control measures in the prevention of the disease are illustrated graphically in Fig. 3a–g. The sub-plots 3(a), 3(d) and 3(e) show that the susceptible individuals of both the populations and recovered of the host population are increasing very well in a comparison of without control case. The individuals in the infected humans and vectors are vanishing after 30 and 50 days as shown in 3(c) and 3(g), respectively. The population of exposed vectors has vanished after 40 days see 3(f), whereas the exposed humans are although decreased but not vanished up to 100 days. Figure 3h describes the control profile. The graphical interpretation of this case reveals that this strategy is useful and can be used for disease eradication.

Strategy 2: prevention and insecticides spraying controls ()

Here, in this case, we analyze the combined influence of prevention through bed-nets and use of insecticides spraying control measures. To get the desired interpretation, we ignore the treatment control, i.e., in the above control model (30). The resulting numerical simulation results of this strategy are depicted in Fig. 4a–g. From the graphical results, we observe that the susceptible and exposed humans population show almost the same behavior as in the first strategy see 4(a) and (b). Moreover, from the remaining graphical results, we can see that this strategy is not effective and is not suitable for a disease eradication in a community.

Strategy 3: prevention and treatment controls ()

In this case, we consider the simultaneous impact of the prevention control and treatment control on the disease dynamics while ignoring the insecticide spray control . The simulation results of this case are illustrated in Fig. 5 with sub-graphs (a)–(g). The graphical interpretation of this strategy demonstrates that this combination of control measures is better than the strategies 1 and 2 as the individuals in the susceptible and recovered humans increased significantly with the activation of control measures see Fig. 5a, d. Similarly, the susceptible vectors are also increased very well as shown in Fig. 5e. On the other hand, a slight faster decrease is observed in the population of exposed and infected humans than the previous two strategies see Fig. 5b and (c). Unlike strategy 1, the decrease in infected vectors is comparatively slower than strategy 1 as depicted in Fig. 5f. The corresponding control profile of this case is depicted in 5(h). Overall, this strategy is useful and can be used for the disease eradication.

Strategy 4: prevention, treatment and insecticides spraying controls ()

From the above discussion, we observe that although the first and third strategies are useful against the infected population both humans and vectors but these strategies are not effective collectively for all compartments. So, the above strategies are not significantly helpful for the eradication of the infection. Therefore, now we activate all the possible three controls at a time, i.e., , , . The simulation results of this case are depicted in Fig. 6a–g and the control profile is shown in Fig. 6h. From the comparison of this final strategy with all the previous three strategies, we notice that this strategy is more effective for the elimination of the infection from the humans as well as vectors population.

Conclusion

Vector-host or vector-borne diseases remain a major public health issue and cause thousands deaths each year. More than 129 countries are at risk of these infections. In this paper, we investigated the dynamics of vector-host infectious disease using a new deterministic mathematical model. Unlike the existing vector-host models, we used nonlinear incidence rate and saturated treatment rate to formulate the proposed model. This type of treatment rate function is also known as Holling type II which accounts for the delay in treatment. Initially, the model is constructed with no control measures. The basic mathematical analysis such as the basic reproduction number, equilibria, boundedness, and positivity of the vector-host model is explored. Further, the local and global stability analysis of the DFE and EE are proven with different stability approaches. Using central manifold theory, the existence of the bifurcation phenomena in the model is discussed. Besides the theoretical analysis of the model, we carried out the global sensitivity of the basic reproduction number for different relevant important parameters. Keeping the sensitivity analysis in mind, we reformulated the control model by incorporating three-time dependent control variables namely, prevention via bed-nets , treatment , and the insecticide spraying . Moreover, we presented four strategies using different combinations of the aforementioned control variables. The impact of each strategy on the disease dynamics is analyzed graphically and discussed in detail. From the graphical interpretations, it is found that although all strategies except strategy 2 are useful and can be used for disease minimization. But, from the simulation of strategy 4, we conclude that the implementation of all control interventions at the same time is the most useful strategy against the disease eradication. We believe that the analysis conducted in this paper will be beneficial for decision-making and health authorities to combat the disease prevalence. In the future, we will extend the present model to fractional order to explore more insights about the disease dynamics.