Modeling and Optimal Control on the Spread of Hantavirus Infection

Abstract

In this paper, optimal control theory is applied to a system of ordinary differential equations representing a hantavirus infection in rodent and alien populations. The effect of the optimal control in eliminating the rodent population that caused the hantavirus infection is investigated. In addition, Pontryagin’s maximum principle is used to obtain the necessary condition for the controls to be optimal. The Runge–Kutta method is then used to solve the proposed optimal control system. The findings from the optimal control problem suggest that the infection may be eradicated by implementing some controls for a certain period of time. This research concludes that the optimal control mathematical model is an effective method in reducing the number of infectious in a community and environment.

1. Introduction

Hantavirus infection is a zoonotic disease associated with hemorrhagic fever with renal syndrome (HFRS), hantavirus cardiopulmonary syndrome (HCPS), and hantavirus pulmonary syndrome (HPS) in humans [1]. Hantaviruses are transmitted between their vertebrate hosts through aggressive interactions or the inhalation of excreta [2]. Hantaviruses are hosted by a variety of rodents and disseminated but do not affect the hosts themselves. The emergence of hantavirus is often linked with ecological interactions among reservoirs, viruses, and humans that are heavily influenced by densities and structures of rodent population, environmental factors associated with survival of the virus outside rodent dwellings, and anthropogenic risk behavior factors. These factors allow prolonged contact to host habitats and areas supporting virus survival [3].

Pathogenic hantaviruses, in nature carried by a specific rodent host species, can cause severe diseases in humans with mortality rates ranging from 12% (HFRS) to 40% (HCPS) [4]. A simple mathematical model describing the spread of the hantavirus infection in rodents has been developed by [5], wherein the model takes into account the temporal and spatial characteristics of this infection. This model has been explored extensively by [4,6,7,8,9,10,11,12].

In a natural ecosystem, rodents share the resource and environment with another species (that shall be called an "alien" species). Consequently, the competition between rodents and alien species should be taken into account in a mathematical model of hantaviruses. This results in competition for resources [8]. As stated by [13] rodents undergo two types of competition: intraspecific, among their own species, and interspecific, with aliens.

Biodiversity has the effect of increasing competition pressure on a host species, which may result in the reduction or elimination of hantavirus [8]. Considerable work has been done on modeling the effects of biodiversity on the spread of disease using simple mathematical models based on systems of ordinary differential equations. These works include [14,15,16].

According to [17], harvesting is defined as reduction of the population due to hunting or capturing individuals, which in effect removes individuals from the population. Such a policy has been used to stabilize the population in an environment with limited resources or environmental parameters. For additional discussion on this matter, see [12,18].

Hantavirus and leptospirosis live in rodents. Both diseases are caused by thrombotic microangiopathy with hemorrhagic phenomena, as well as hepatic and pulmonary involvement [19]. Given that there are similarities between these two diseases, we use leptospirosis as a study guide in studying optimum control strategies to control transmission of hantavirus. According to [20], optimal control is the establishment of control and state trajectories for a dynamic scheme over a time period to minimize a performance index.

Mathematical models of leptospirosis disease based on optimum control strategies have been developed. Such a model was developed by [21] in order to study the prevention of a leptospirosis-infected vector and human population using multiple control variables. Okosun et al. [22] investigated the effectiveness of leptospirosis control measures, preventive vaccination, and treatment of infective humans in reducing disease transmission. In their paper, the authors showed that there is local and global stability of the proposed model. They also analyzed the constant control parameters, calculated the basic reproduction number, and investigated the existence and stability of equilibria. They applied the optimal control method using Pontryagin’s maximum principle to determine the necessary conditions for the control measures for leptospirosis disease.

In our present paper, the model of [8] of the effect of alien species on the hantavirus infection, which was extended by [12] to incorporate harvesting, is now extended to incorporate optimum control strategies. This is to control the transmission of the spread of the hantavirus. Optimal control theory is employed to diminish the proportion of the population of infected rodents and aliens using various controls. According to [23], optimal control theory is a tool which enables decisions involving complex dynamical systems to be made. The aim of this paper is to develop a deterministic mathematical model for hantavirus transmission and optimal control in rodent and alien populations.

2. The Optimal Control Model

In this section, the basic model of hantavirus proposed and developed by Abramson and [5] is considered. The model proposed a single rodent species with no movement from one location to another. In this model, the total population of rodents is divided into two groups: one is susceptible and the other is infected. The model is given by:

$$\begin{array}{l}\frac{d{r}_s}{dt}=br-c{r}_s-\frac{r_{s}r}{k(t)}-a{r}_s{r}_i,\\ \frac{d{r}_i}{dt}=-c{r}_i-\frac{r_{i}r}{k(t)}+a{r}_s{r}_i.\end{array}$$

where the susceptible rodent is denoted by $r_s$, the infected rodent is denoted by $r_i$, and the total population of rodents is given by $r\left(t\right)=r_s\left(t\right)+r_i\left(t\right)$. The birth rate is denoted by b, the natural mortality rate is denoted by c, the transmission rate of hantavirus is denoted by a (also called the aggression parameter), and the environmental parameter is denoted by k. Based on the Abramson and Kenkre model, the infection dies away for $k<k_c$, where k_c = b/(a(b − c)) and the infection thrives since there is an increase in edible resources for $k>k_c$. For additional discussion on this matter, see [5].

As mentioned earlier, [8] developed a model which incorporated biodiversity. Yusof et al. [12] included the effects of harvesting in the model of Peixoto and Abramson. The model is of the form:

$$\begin{array}{l}\frac{d{r}_s}{dt}=br-c{r}_s-\frac{r_s}{k}\left(r+q{z}_a\right)-a{r}_s{r}_i-E{r}_s,\\ \frac{d{r}_i}{dt}=-c{r}_i-\frac{r_i}{k}\left(r+q{z}_a\right)+a{r}_s{r}_i-E{r}_i,\\ \frac{d{z}_a}{dt}=\left(\beta -\gamma \right)z_a-\frac{z_a}{k}\left(z_a+\epsilon r\right).\end{array}\}$$

(1)

where the total rodent population, $r\left(t\right)$$\left(\mathrm{where}r\left(t\right)=r_s\left(t\right)+r_i\left(t\right)\right)$ is divided into susceptible $\left(r_s\left(t\right)\right)$ and infected $\left(r_i\left(t\right)\right)$. Then, E represents the harvesting efforts, $q$ is the influence of the alien population, and $z_a\left(t\right)$ is the population of the alien. For the alien population, $\beta ,\gamma \mathrm{and}\epsilon$ are the corresponding parameters to obtain resources from other species, and k is an environmental parameter.

The model of the optimum control problem applied to a dynamical leptospirosis as described in Khan et al. [21] can be written as:

$$\begin{array}{l}\frac{d{S}_h}{dt}=b_1-{\mu }_h{S}_h-{\beta }_2{S}_h{I}_v\left(1-u_1\left(t\right)\right)-{\beta }_1{S}_h{I}_h\left(1-u_2\left(t\right)\right)+{\lambda }_h{R}_h,\\ \frac{d{I}_h}{dt}={\beta }_2{S}_h{I}_v\left(1-u_1\left(t\right)\right)+{\beta }_1{S}_h{I}_h\left(1-u_2\left(t\right)\right)-\left({\mu }_h+{\delta }_h+{\gamma }_h+u_2\left(t\right)\right)I_h,\\ \frac{d{R}_h}{dt}={\gamma }_h{I}_h-\left({\mu }_h+{\lambda }_h\right)R_h+u_1\left(t\right)S_h+u_2\left(t\right)I_h,\\ \frac{d{S}_v}{dt}=b_2-{\gamma }_v{S}_v-{\beta }_3{S}_v{I}_h-{\in }_1{u}_3\left(t\right)S_v,\\ \frac{d{I}_v}{dt}={\beta }_3{S}_v{I}_h-\left({\gamma }_v+{\delta }_v\right)I_v-{\in }_2{u}_3\left(t\right)I_v\end{array}\}$$

(2)

with initial conditions $S_h\left(0\right)\ge 0$, $I_h\left(0\right)\ge 0$, $R_h\left(0\right)\ge 0$, $S_v\left(0\right)\ge 0$, $I_v\left(0\right)\ge 0$.

Here the total human population, $\left(N_h\left(t\right)=S_h\left(t\right)+I_h\left(t\right)+R_h\left(t\right)\right)$ is divided into susceptible, $S_h$, infected, $I_h$, and recovered, $R_h$. Then, the total vector population, $\left(N_v\left(t\right)=S_v\left(t\right)+I_v\left(t\right)\right)$ is divided into two groups: one is susceptible, $S_v$, and the other is infected, $I_v$. The parameters of Equation (2) are described in

Table 1. Description of the parameters used in Equation (2).

Symbol	Description
b1	The recruitment rate of human population
β1, β2	The mediate transmission coefficients
λh	The immune human once again is susceptible at constant rate
μh	The natural mortality rate for human
δh	The disease related death rate from infected class occurs at human population
γh	The recovery rate for human from the infections
b2	The recruitment rate for vector population
γv	The death rate of vector
β3	The disease carrying of susceptible vector per host per unit time
δv	The disease death rate of vector populations
∈1, ∈2	The rate of positive constants.
u1(t)	Human should cover all cuts, water dry, full-cover boots, shoes and long sleeve shirts when handling animals
u2(t)	Wash hands thoroughly on a regular basis and shower after work
u3(t)	Clean up both work place and home

.

For the human population, the associated force of infections has decreased by factors of $\left(1-u_1\left(t\right)\right)$ and $\left(1-u_2\left(t\right)\right)$, respectively. They assume that the natural death rate of vector population increases at a rate proportional to $u_3\left(t\right)$ where ${\in }_1>0$ and ${\in }_2>0$ are the positive constant rates.

Using [15] as a reference, optimal control theory was applied to the hantavirus infection with the presence of an alien and hantavirus. The new hantavirus infection with optimal control model is given by:

$$\begin{array}{l}\frac{d{r}_s}{dt}=br-c{r}_s-\frac{r_s}{k}\left(r+q{z}_a\right)-a{r}_s{r}_i-E\left(t\right)r_s,\\ \frac{d{r}_i}{dt}=-c{r}_i-\frac{r_i}{k}\left(r+q{z}_a\right)+a{r}_s{r}_i-E\left(t\right)r_i,\\ \frac{d{z}_a}{dt}=\left(\beta -\gamma \right)z_a-\frac{z_a}{k}\left(z_a+\left(1-u\left(t\right)\right)\epsilon r\right)\end{array}\}$$

(3)

where $E\left(t\right)$ and $u\left(t\right)$ are the fraction of the population of rodents removed for each time period, i.e., harvesting efforts and biodiversity refers to the distribution and abundance of rodent and alien populations within an ecosystem, respectively. The value $\left(1-u\left(t\right)\right)$ is used to diminish the force of the hantavirus infection. The control variables being engaged, $E\left(t\right)$ and $u\left(t\right)$, are bounded Lebesgue measurable functions on the interval. The main objective in this control variable is to minimize both populations of the rodents.

3. Analyzing the Model

In this section, we analyze the optimal control and the model proposed in Equation (3). In this analysis, we follow the work of [24]. Pontryagin’s maximum principle was used to govern the necessary conditions for the optimal control of the spread of the hantavirus infection. It is desired to minimize the number of rodents spreading the hantavirus infection using harvesting and biodiversity controls. Therefore, the goal is to find a pattern to minimize the number of rodents spreading this infection using the functional given by:

$$J\left(E\left(t\right),u\left(t\right)\right)={{\displaystyle \int }}_0^T\left[\left(A_1{r}_s+A_2{r}_i\right)+\frac{1}{2}\left(B_1{E}^2+B_2{u}^2\right)\right]dt.$$

We thus use the bounded Lebesgue measurable functions and define our objective functional as:

$$J\left(E\left(t\right),u\left(t\right)\right)=\min {\displaystyle {\int }_0^T\left[\left(A_1{r}_s+A_2{r}_i\right)+\frac{1}{2}\left(B_1{E}^2+B_2{u}^2\right)\right]}dt$$

(4)

where $B_1$ and $B_2$ are weight constants for harvesting efforts and biodiversity, respectively. For the objective functional, the weight constant of the susceptible rodent is denoted by $A_1$ and the weight constant of infected rodents is denoted by $A_2$. Here, we seek the optimal control problem $\left(E*\left(t\right),u*\left(t\right)\right)$ such that:

$$J\left(E*\left(t\right),u*\left(t\right)\right)=\min \left\{J\left(E\left(t\right),u\left(t\right)\right)|E\left(t\right),u\left(t\right)\in U\right\}$$

where the control set is defined as $U=\left\{\left(E\left(t\right),u\left(t\right)\right)\mathrm{are}\mathrm{Lebesque}\mathrm{measureable},0\le \left(E\left(t\right),u\left(t\right)\right)\le 1,t\in \left[0,T\right]\right\}$.

The Lagrangian for the optimal control problem is defined by

$$L=A_1{r}_s+A_2{r}_i+\frac{1}{2}\left(B_1{E}^2+B_2{u}^2\right)$$

The necessary conditions that an optimal solution must satisfy the maximum principle of [25]. This principle converts the model of Equations (3) and (4) into a problem of minimizing point-wise a Hamiltonian, H, with respect to $E\left(t\right)$ and $u\left(t\right)$ as:

$$\begin{gathered} H=L\left(r_s,r_i,E\left(t\right),u\left(t\right)\right)+{\lambda }_1\frac{d{r}_s}{dt}+{\lambda }_2\frac{d{r}_i}{dt}+{\lambda }_3\frac{d{z}_a}{dt} \\ =A_1{r}_s+A_2{r}_i+\frac{1}{2}{B}_{1}E{\left(t\right)}^2+\frac{1}{2}{B}_{2}u{\left(t\right)}^2+{\lambda }_1\frac{d{r}_s}{dt}+{\lambda }_2\frac{d{r}_i}{dt}+{\lambda }_3\frac{d{z}_a}{dt} \\ =A_1{r}_s+A_2{r}_i+\frac{1}{2}{B}_{1}E{\left(t\right)}^2+\frac{1}{2}{B}_{2}u{\left(t\right)}^2+{\lambda }_1\left(br-c{r}_s-\frac{r_s}{k}\left(r+q{z}_a\right)-a{r}_s{r}_i-E\left(t\right)r_s\right)+ \\ {\lambda }_2\left(-c{r}_i-\frac{r_i}{k}\left(r+q{z}_a\right)+a{r}_s{r}_i-E\left(t\right)r_i\right)+{\lambda }_3\left(\left(\beta -\gamma \right)z_a-\frac{z_a}{k}\left(z_a+\left(1-u\left(t\right)\right)\epsilon r\right)\right) \end{gathered}$$

(5)

where ${\lambda }_1$, ${\lambda }_2$, and ${\lambda }_3$ are the adjoint variables or co-state variables.

From Equation (5), the following adjoint equations can be obtained:

$$\begin{gathered} \frac{d{\lambda }_1}{dt}=-\frac{\delta H}{\delta r_s}=-A_1-{\lambda }_1\left(b-c-\frac{1}{k}\left(2r_s+r_i+q{z}_a\right)-a{r}_i-E\left(t\right)\right)-{\lambda }_2\left(-\frac{r_i}{k}+a{r}_i\right)- \\ {\lambda }_3\left(-\frac{z_a}{k}\left(\left(1-u\left(t\right)\right)\epsilon \right)\right) \\ \frac{d{\lambda }_2}{dt}=-\frac{\delta H}{\delta r_i}=-A_2-{\lambda }_1\left(b-\frac{r_s}{k}-a{r}_s\right)-{\lambda }_2\left(-c-\frac{1}{k}\left(r_s+2r_i+q{z}_a\right)+a{r}_s-E\left(t\right)\right)- \\ {\lambda }_3\left(-\frac{z_a}{k}\left(\left(1-u\left(t\right)\right)\epsilon \right)\right) \\ \frac{d{\lambda }_3}{dt}=-\frac{\delta H}{\delta z_a}={\lambda }_1\left(\frac{q}{k}{r}_s\right)+{\lambda }_2\left(\frac{q}{k}{r}_i\right)-{\lambda }_3\left(\left(\beta -\gamma \right)-\frac{1}{k}\left(2z_a+\left(1-u\left(t\right)\right)\epsilon r\right)\right) \end{gathered}$$

(6)

with transversality condition (or boundary condition):

$${\lambda }_1\left(t_{end}\right)={\lambda }_2\left(t_{end}\right)={\lambda }_3\left(t_{end}\right)=0$$

Then, the Hamiltonian, H, is differentiated, with respect to $E\left(t\right)$ and $u\left(t\right)$, and setting $\delta H/\delta E\left(t\right)=\delta H/\delta u\left(t\right)=0$.

$$\begin{gathered} \frac{\delta H}{\delta E\left(t\right)}=B_{1}E-{\lambda }_1{r}_s-{\lambda }_2{r}_i=0 \\ \frac{\delta H}{\delta u\left(t\right)}=B_{2}u\left(t\right)+\frac{{\lambda }_{\epsilon }\epsilon }{k}r{z}_a=0 \end{gathered}$$

Then, solving for optimal control $\left(\mathrm{i}.\mathrm{e}.\frac{\delta H}{\delta E\left(t\right)}=\frac{\delta H}{\delta u\left(t\right)}=0\right)$, yields:

$$\begin{gathered} E{\left(t\right)}^C=\frac{1}{B_1}\left({\lambda }_1{r}_s+{\lambda }_2{r}_i\right) \\ u{\left(t\right)}^C=-\frac{{\lambda }_{\epsilon }\epsilon }{k{B}_2}r{z}_a \end{gathered}$$

Therefore, the characterizations of the optimal controls $E{\left(t\right)}^{*}$ and $u{\left(t\right)}^{*}$ are given by:

$$\begin{array}{l}E\left(t\right)*=\max \left\{\min \left\{\frac{1}{B_1}\left({\lambda }_1{r}_s+{\lambda }_2{r}_i\right),1\right\},0\right\}\\ u\left(t\right)*=\max \left\{\min \left\{-\frac{{\lambda }_{\epsilon }\epsilon }{k{B}_2}r{z}_a,1\right\},0\right\}.\end{array}\}$$

(7)

By standard arguments which involve the bounds on the controls, it can be concluded that:

$$\begin{gathered} E\left(t\right)*=\{\begin{array}{cc}0& \mathrm{if}E{\left(t\right)}^C\le 0,\hfill \\ E{\left(t\right)}^C& \mathrm{if}0<E{\left(t\right)}^C<1,\hfill \\ 1& \mathrm{if}E{\left(t\right)}^C\ge 1,\hfill \end{array} \\ u\left(t\right)*=\{\begin{array}{cc}0& \mathrm{if}u{\left(t\right)}^C\le 0,\hfill \\ u{\left(t\right)}^C& \mathrm{if}0<u{\left(t\right)}^C<1,\hfill \\ 1& \mathrm{if}u{\left(t\right)}^C\ge 1.\hfill \end{array} \end{gathered}$$

4. Numerical Results and Discussion

In this paper, numerical results from the proposed model are discussed and compared with an existing model that is without optimum control (Equation (1)). First, the model without optimum control and with the optimum control model are solved using a fourth-order Runge–Kutta scheme with MATLAB software for 50 months. The optimal control set is obtained by solving the optimality system, which consists of the state and adjoint systems with the transversality conditions or boundary conditions [21,26,27,28]. The weight constant values in the objective functional are chosen to be $A_1=100$, $A_2=0.03$, $B_1=100$, and $B_2=140$. These values are similar to [19]. The complete list of the parameter values used for the numerical simulations of the optimum control model is tabulated in

Table 2. Parameter values in the numerical simulations.

Symbol	Description	Parameter value
$b$	The birth rate of rodent species	1
$c$	The death rate of rodent species	0.6
$a$	The transmission rate of hantavirus	0.1
$k$	The environmental parameter	Varied
$q$	Competitive effect of rodents on alien species	0.2
$\beta$	The birth rate of alien species	1
$\gamma$	The death rate of alien species	0.5
$\epsilon$	Competitive effect of rodents on alien species	0.1
$E\left(t\right)$	The harvesting efforts	Varied
$u\left(t\right)$	The biodiversity	Varied

, which is based on empirical studies by [8].

The dynamics of the models without optimum control (Equation (1)) and with optimum control (Equation (3)) are shown. Figure 1, Figure 2, and Figure 3 represent the susceptible rodent, infected rodent, and alien populations, respectively, in the two models, i.e., without control and with control. It is apparent in all of the figures that there are considerable differences between the number of rodents and alien populations under control as compared to cases without control.

In Figure 1, it can be observed that the susceptible rodent population $r_s$ will disappear and stabilize at a certain steady value when the two controls, i.e., harvesting efforts $E\left(t\right)$ and biodiversity $u\left(t\right)$, are applied. The susceptible rodent population behaves in the opposite way in the model without optimum control, wherein it increases sharply initially and reaches a certain maximum, before plunging and approaching an equilibrium quantity.

For the optimum control model (Equation (3)), Figure 2 shows the number of infected rodents $r_i$ decreases strongly and drops to extinction thereafter, while the population of infected rodents $r_i$ thrives and approaches an equilibrium stable value for the model without optimum control (Equation (1)). The results illustrated in Figure 2 clearly show that the activated control strategy is effective in eliminating the infected rodent population with the presence of harvesting efforts $E\left(t\right)$ and biodiversity $u\left(t\right)$ as optimum control, which is not the case without control. The use of harvesting efforts $E\left(t\right)$ and biodiversity $u\left(t\right)$ can eliminate the number of infected rodents, and this indicates that the hantavirus infection will die off with the presence of two controls.

The numerical simulation results in Figure 3 display a similar pattern to Figure 2 for the models without control (Equation (1)) and with control (Equation (3)). When the optimum control model is considered, we observe that the population of alien $z_a$ decreases with a slower rate and eventually drops to extinction. In the model without optimum control, the alien population can survive and stabilize at a certain steady value. Figure 4 shows harvesting efforts $E\left(t\right)$ and biodiversity $u\left(t\right)$ as control variables. The optimum control increases in month 50, and a small population of susceptible rodents remains, while both the infected rodent and alien populations die off.

Next, we consider the simulation result as the environmental parameter, k, is decreased to 20. The results for k = 20 are as shown in Figure 5, Figure 6, Figure 7 and Figure 8. Note that initial conditions for case k = 20 are given by $r_s=5$, $r_i=5$, and $z_a=5$. Figure 5 exhibits a similar transmission as Figure 1 as for the cases with and without control. In Figure 6, the illustrated result shows the number of infected rodents decreases strongly and drops to extinction. In this case, as the environmental parameter is decreased, the number of infected rodents in the system with no control will also drop to extinction, but at a slower rate compared to a system with control parameters. Figure 7 also shows similar trends to the result in Figure 2 but with a faster rate. In Figure 8, we show harvesting efforts and biodiversity as control variables. The optimum control for harvesting increases in month 50, while the biodiversity factor is reducing. A small population of susceptible rodents remains, while both the infected rodent and alien populations die off. Finally, the results from the two cases presented above indicate that effects on environmental capacity, such as availability of food sources, might also impact the spread of the infection.

5. Conclusions

In this paper, we studied the effects of the harvesting efforts $E\left(t\right)$ and biodiversity $u\left(t\right)$ as an optimum control on the spread of the hantavirus infection. First, we proposed a new mathematical model that takes into account the presence of an alien population and the hantavirus using optimum control theory. The theoretical studies of optimal control models and their numerical simulations are discussed in this paper. Then, the optimum controls are applied to minimize both populations of rodents. The results from this model clearly show that the optimum control is able to eliminate the infected rodent populations with the presence of harvesting efforts and biodiversity as the controls, as well as in the case without control.