Next Article in Journal
Fractal Behavior of a Ternary 4-Point Rational Interpolation Subdivision Scheme
Next Article in Special Issue
Modeling and Simulation of a Hydraulic Network for Leak Diagnosis
Previous Article in Journal
A Computational Method with MAPLE for a Piecewise Polynomial Approximation to the Trigonometric Functions
Previous Article in Special Issue
Differential Evolution Algorithm for Multilevel Assignment Problem: A Case Study in Chicken Transportation
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Time Needed to Control an Epidemic with Restricted Resources in SIR Model with Short-Term Controlled Population: A Fixed Point Method for a Free Isoperimetric Optimal Control Problem

Department of Mathematics and Computer Sciences, Faculty of Sciences Ben M'Sik, Hassan II University of Casablanca, Casablanca 20000, Morocco
*
Author to whom correspondence should be addressed.
Math. Comput. Appl. 2018, 23(4), 64; https://doi.org/10.3390/mca23040064
Submission received: 10 October 2018 / Revised: 20 October 2018 / Accepted: 20 October 2018 / Published: 22 October 2018
(This article belongs to the Special Issue Optimization in Control Applications)

Abstract

:
In this paper, we attempt to determine the optimal duration of an anti-epidemic control strategy which targets susceptible people, under the isoperimetric condition that we could not control all individuals of this category due to restricted health resources. We state and prove the local and global stability conditions of free and endemic equilibria of a simple epidemic compartmental model devised in the form of four ordinary differential equations which describe the dynamics of susceptible-controlled-infected-removed populations and where it is taken into account that the controlled people cannot acquire long-lived immunity to move towards the removed compartment due to the temporary effect of the control parameter. Thereafter, we characterize the sought optimal control and we show the effectiveness of this limited control policy along with the research of the optimal duration that is needed to reduce the size of the infected population. The isoperimetric constraint is defined over a fixed horizon, while the objective function is defined over a free horizon present under a quadratic form in the payoff term. The complexity of this optimal control problem requires the execution of three numerical methods all combined together at the same time, namely, the forward–backward sweep method to generate the optimal state and control functions, the secant method adapted to the isoperimetric restriction, and, finally, the fixed point method to obtain the optimal final time.

1. Introduction

1.1. Background

Many epidemiological models have been interested in the study of the dynamics of susceptible, infected and removed individuals who belong to a sample of a population threatened by an infection. Many theoretical models in epidemiology have been devised to show the effect of different anti-epidemic control strategies when they are followed to prevent transmission of a particular type of infection to the susceptible population. As examples of these control approaches, we can cite Refs. [1,2] where the authors introduced an awareness control function in their models and which aimed to prevent the susceptible people from Human Immunodeficiency Virus infection and Acquired Immune Deficiency Syndrome (HIV/AIDS) epidemic. Roy et al. [3] treated the idea of awareness control programs in HIV/AIDS prevention after the addition of a new variable in their models and which defines the number of individuals in the aware class. Other examples of control models have used vaccination of the susceptible individuals as a control policy, while considering the number of vaccinated people as additional compartment in their systems (see [3,4,5,6,7]). Most optimal control strategies suggested for preventing an infection to spread consider a fixed final time. However, health policy-makers need very often to know when it is appropriate to stop their anti-epidemic measures as this information is important for managing their medical resources [8], and then, it becomes not reasonable in such situations to study dynamics of an epidemic under control without an estimation of the final time. The present paper tries to find the optimal value of this mentioned variable through a free horizon optimal control approach applied to an epidemic model with four compartments, namely classes of susceptible, controlled, infected and removed people.

1.2. Formulation of the Problem of Interest for this Investigation

Here, we devise a simple generalized model where the control utilized has only a temporary effect on the immunity of the targeted population so the susceptible people under control do not acquire long-lived immunity to move to the removed class, while taking into account the presence of an equation which describes the evolution of the number of the controlled people. First, we study the local and global stability of our epidemic model, which is devised in the form of four ordinary differential equations, and wherein a control is introduced as a constant parameter; and, second, we seek the optimal duration needed for reaching the goal of our strategy while determining the optimal value of this control when it is changed to a function of time under the hypothesis that this anti-epidemic preventive measure can reach only a specific fraction of susceptible people due to the limited health resources. We should note that authors of [8,9] tried to find the optimal final time in a first case through their plots and checked the values where this optimal function verified the obtained additional necessary condition, and in two other cases via discrete numerical schemes. We believe this work is more interesting, as it provides a more precise numerical method. In fact, as the necessary condition on final time also represents here a fixed point equation, the incorporation of the fixed point method better facilitates our task, as this technique seems more accurate and convincing since it meets the theoretical aspect of the found condition.

1.3. Literature Survey

Zhou and Fan [10] discussed different forms of functions introduced in epidemic models to explore the impact of limited medical resources in the transmission of infectious diseases. Abdelrazec et al. [11] introduced, in a mathematical model of dengue fever, a function in place of the recovery rate for similar purpose. More recently, Yu et al. proposed an optimal control approach to investigate the optimal distribution policy of the limited vaccination resources based on the research of a parameter introduced in their model and which minimizes the basic reproduction number [12].

1.4. Scope and Contribution of this Study

We present here an optimal control approach which treats the problem of limited resources differently to the three above-mentioned references. In fact, our method considers a constraint, so-called “isoperimetric”, which is used on the control function as done in [13] for the resolution of a dosage problem, on the fraction of controlled variables as done in [14] in the case of an epidemic model, or when adapted to a discrete-time SIRS epidemic model as in [15], and it supposes also that the final time or horizon of the objective function is free (non-fixed) as used in many applications (see, for example, Ref. [8], which discusses the problem of optimal duration needed for reaching the intravesical therapy goal, and Ref. [9], where it is explained why such considerations are very important to health-policy makers and managers in the health sector when there is an epidemic that is controlled through awareness of the susceptible class).

1.5. Organization of the Paper

Based on the theory of mathematical epidemiology in [16], the spread of an epidemic can be described mathematically by SIR models which in turn have been developed later to extended forms such as SEIR [17], SIRS [18], or SIS [19], where each letter refers to a class of individuals. A class of controlled people can also be considered, as done in [20] where authors added a vaccination compartment in a model of pertussis and tuberculosis, and in [21] where they studied nonfatal diseases, and [22,23] in case of influenza. Sharomi and Malik [24] represented an other form of SIR model with an additional equation corresponding to the vaccinated category in the case the vaccination is not 100 % effective; such considerations can also be found in [4]. Based on similar assumptions as in the two last mentioned references, we devise our present model.
In the following parts of the paper, we start with the presentation of our mathematical model and study its stability in cases of free and endemic equilibria. Furthermore, we seek the optimal value of the free horizon considered in the objective function, along with the determination of the optimal value of the control function. Finally, we discuss our numerical results.

2. The Mathematical Model and Stability

In this section, we consider a mathematical model with the four following main compartments:
  • S is the number of susceptible people to infection or who are not yet infected.
  • C S is the number of susceptible people who are temporarily controlled so they cannot move to the removed class due to the limited effect of control.
  • I is the number of infected people who are capable of spreading the epidemic to those in the susceptible and temporary controlled categories.
  • R is the number of removed people from the epidemic.
In our modeling approach, we choose to describe dynamics of variables S, C S , I and R at time t, based on the following differential system
S ˙ ( t ) = Π ( t ) β S ( t ) I ( t ) a θ S ( t ) μ S ( t ) C S ˙ ( t ) = a θ S ( t ) b β C S ( t ) I ( t ) μ C S ( t ) I ˙ ( t ) = β ( S ( t ) + b C S ( t ) ) I ( t ) γ I ( t ) μ I ( t ) R ˙ ( t ) = γ I ( t ) μ R ( t )
with initial conditions S ( 0 ) > 0 , C S ( 0 ) 0 , I ( 0 ) 0 and R ( 0 ) 0 , and where Π ( t ) = μ N ( t ) , with N ( t ) = S ( t ) + C S ( t ) + I ( t ) + R ( t ) as the total population size, gives the newborn people at time t; a θ ( 0 a 1 ) is the recruitment rate of susceptibles to the controlled class with θ defining the control parameter as a constant between 0 and 1 and “a” modeling the reduced chances of a susceptible individuals to be controlled; β = δ N ( t ) with δ the infection transmission rate, μ the natural death rate, b θ ( 0 b 1 ) the recruitment rate of controlled people to the infected class even in the presence of θ and “b” modeling the reduced chances of a temporarily controlled individual to be infected; and γ is the recovery rate. We note that the population size is constant because N ˙ ( t ) = S ˙ ( t ) + C ˙ S ( t ) + I ˙ ( t ) + R ˙ ( t ) = 0 . Hence, N ( t ) = N = a c o n s t a n t , and then, Π ( t ) = Π = a c o n s t a n t .
For the sake of readability, hereafter, we use S, C S , I and R as notations of the time functions S ( t ) , C S ( t ) , I ( t ) and R ( t ) .
Recalling that R 0 = β μ + γ is the basic reproduction number of the standard SIR model (see [25] where it is concluded that the disease free equilibrium E 0 is global asymptotically stable if R 0 1 , and there exists a global asymptotically stable and unique endemic equilibrium E + if R 0 > 1 ).
Since the two first equations are independent of the last equation, we only study the stability of the following differential system
S ˙ = Π β S I a θ S μ S C S ˙ = a θ S b β C S I μ C S I ˙ = β ( S + b C S ) I γ I μ I
A disease free equilibrium in our case can be defined as E 0 = ( S 0 , C S 0 , 0 ) where S 0 and C S 0 are obtained based on the assumptions S ˙ = 0 and I ˙ = 0 when there is no infection.
Explicitly, we have S ˙ = 0 when I = 0 , gives S 0 = Π μ + a θ . In addition, we have C S ˙ = 0 when I = 0 , gives C S 0 = a θ Π μ ( μ + a θ ) .
If lim t + C S ( t ) = 0 as a consequence of the case when θ = 0 , we define the basic reproduction number for our case by R 0 C which is the average new infections produced by one infected individual during his life cycle when the population is at E 0 .
Since I is the only infected compartment, then R 0 C = β ( S 0 + b C S 0 ) × 1 μ + γ . Thus, we have
R 0 C = β Π ( μ + a θ ) ( μ + γ ) + a b β θ Π μ ( μ + a θ ) ( μ + γ ) = β Π ( μ + a b θ ) μ ( μ + a θ ) ( μ + γ )
Now, we try to find the components of the endemic equilibrium E + = ( S + , C S + , I + ) where S + and C S + are obtained based on the assumptions S ˙ = 0 and C S ˙ = 0 when there is an infection.
Explicitly, we have S ˙ = 0 when I > 0 , which gives S + = Π μ + a θ + β I + . In addition, we have C S ˙ = 0 when I > 0 , which gives C S + = a θ S + μ + b β I + .
On the other part, we have I ˙ = 0 when I > 0 , which gives
β S + I + + b β C S + I + γ I + μ I + = 0 β Π μ + a θ + β I + + β a b θ S + μ + b β I + = γ + μ β Π ( μ + b β I + ) + a b β θ S + ( μ + a θ + β I + ) = ( γ + μ ) ( μ + a θ + β I + ) ( μ + b β I + ) β Π ( μ + b β I + ) + a b β θ Π = ( γ + μ ) ( μ + a θ + β I + ) ( μ + b β I + ) b β 2 Π I + + β Π ( μ + a b θ ) = ( γ + μ ) ( μ + a θ + β I + ) ( μ + b β I + ) ( γ + μ ) ( μ + a θ ) R 0 C μ ( γ + μ ) ( a θ μ + μ 2 ) = b β 2 ( γ + μ ) I + 2 + [ ( γ + μ ) ( β μ ( 1 + b ) + a b β θ ) b β 2 Π ] I +
Thus, we find that I + is the root of the function f ( I + ) = α 1 I + 2 + α 2 I + + ( 1 R 0 C ) α 3 where α 1 , α 2 and α 3 are constants.
In the following three theorems, we state and prove stability results on free and endemic equilibria.
Theorem 1.
E 0 always exists and is locally asymptotically stable if R 0 C < 1 (respectively, E 0 is unstable if R 0 C > 1 ).
Proof. 
The existence of E 0 is trivial.
For the stability of E 0 , we define the Jacobian Matrix associated to the system in Equation (2) by
β I a θ μ 0 β S a θ b β I μ b β C S β I b β I β ( S + b C S ) γ μ
At E 0 , (4) becomes
a θ μ 0 β S 0 a θ μ b β C S 0 0 0 β ( S 0 + b C S 0 ) γ μ
whose eigenvalues are
λ 1 = μ < 0 , λ 2 = ( μ + a θ ) < 0 λ 3 = β ( S 0 + b β C S 0 ) γ μ = ( γ + μ ) ( R 0 C 1 ) ,
which imply the local asymptotic stability of E 0 when R 0 C < 1 , and its instability when R 0 C > 1 . □
Theorem 2.
The differential system in Equation (2) admits E + = Π μ + a θ + β I + , a θ S + μ + b β I + , I + as the unique positive equilibrium and which is asymptotically stable when it exists, if and only if R 0 C > 1 .
Proof. 
First, we have
α 1 = b β 2 ( γ + μ ) > 0 , α 2 = [ ( γ + μ ) ( β μ ( 1 + b ) + a b β θ ) b β 2 Π ] , α 3 = μ ( γ + μ ) ( μ + a θ ) > 0
For the sufficiency of the existence and uniqueness of E + , so we have α 1 > 0 and since f ( 0 ) = ( 1 R 0 C ) α 3 < 0 if R 0 C > 1 , then f ( I + ) has two real roots, one is positive and the other is negative. For the necessity, let us assume that R 0 C 1 and prove that f ( I + ) has no positive roots. In this case, the first fraction in Equation (3) verifies
β Π ( μ + γ ) ( μ + a θ ) α 2 = ( b β ( μ + a θ ) + μ β ) ( γ + μ ) b β 2 Π ( b β ( μ + a θ ) + μ β ) ( γ + μ ) b β ( μ + a θ ) ( μ + γ ) = μ β ( γ + μ ) > 0 .
Thus, we have α 1 > 0 and since f ( 0 ) = ( 1 R 0 C ) α 3 0 , f ( I + ) is increasing and f ( I + ) > f ( 0 ) 0 , then we reach the non-positivity of the roots.
For the stability of E + , at E + , Equation (4) is defined as
β I + a θ μ 0 β S + a θ b β I + μ b β C S + β I + b β I + β ( S + + b C S + ) γ μ
= Π S + 0 β S + a θ a θ S + C S + b β C S + β I + b β I + 0
whose characteristic equation is λ 3 + σ 1 λ 2 + σ 2 λ + σ 3 and where
σ 1 = Π S + + a θ S + C S + σ 2 = a θ Π C S + + b 2 β 2 C S + I + + β 2 S + I + σ 3 = a b β 2 θ S + I + + a θ β 2 S + 2 I + C S + + b 2 β 2 Π C S + I + S + .
Hence, we have
σ 1 σ 2 σ 3 = Π S + a θ Π C S + + β 2 S + I + + a θ S + C S + a θ Π C S + + b 2 β 2 C S + I + a b β 2 θ S + I + = a θ Π 2 S + C S + + ( Π + a θ + β I + ) β 2 S + I + + a θ S + C S + a θ Π C S + + b 2 β 2 C S + I + a b β 2 θ S + I + = a θ Π 2 S + C S + + ( Π + β I + ) β 2 S + I + + a 2 θ 2 Π S + C S + 2 + a θ S + I + ( β b β ) 2 + a b β 2 θ S + I + > 0
Finally, based on the Routh–Hurwitz Criterion, we deduce the local asymptotic stability of E + . □
Theorem 3.
If R 0 C 1 , then E 0 is globally asymptotically stable. If R 0 C > 1 , then E + is globally asymptotically stable.
Proof. 
We suppose that R 0 C 1 and we prove that E 0 is globally asymptotically stable. Let us define the Lyapunov function by
L 0 = S S 0 S 0 ln S S 0 + C S C S 0 C S 0 ln C S C S 0 + I .
Its derivative is then defined by
L 0 ˙ = S ˙ + C S ˙ + I ˙ S 0 S ˙ S C S 0 C S ˙ C S = Π μ S μ C S μ I γ I S 0 Π S + μ S 0 + β S 0 I + a θ S 0 C S 0 a θ S C S + b β C S 0 I + μ C S 0 .
Since μ = a θ S 0 C S 0 and Π = μ S 0 + a θ S 0 , then this derivative becomes
L 0 ˙ = μ S + 2 μ S 0 + 3 a θ S 0 ( μ + γ β S 0 b β C S 0 ) I a θ C S 0 S C S S 0 S ( μ S 0 + a θ S 0 ) a θ S 0 C S C S 0 = μ S 0 S S 0 + S 0 S 2 a θ S 0 C S C S 0 + S 0 S + C S 0 S C S S 0 3 ( μ + γ ) ( 1 R 0 C ) I
Now, we have S S 0 + S 0 S 2 0 and C S C S 0 + S 0 S + C S 0 S C S S 0 3 0 due to the fact that arithmetic mean is larger than or equals to the geometric mean, and the equalities hold if S = S 0 and C S = C S 0 . Thus, L 0 ˙ 0 which implies the global asymptotic stability of E 0 based on Lyapunov–LaSalle’s invariance principle.
Similarly, we study the global asymptotic stability of E + by considering the following Lyapunov function
L + = S S + S + ln S S + + C S C S + C S + ln C S C S + + I I + I + ln I I + .
The derivative is then defined as
L + ˙ = S ˙ + C S ˙ + I ˙ S + S ˙ S C S + C S ˙ C S I + I ˙ I = Π μ S μ C S ( μ + γ ) I Π S + S + μ S + + β S + I + a θ S + a θ S C S + C S + b β C S + I + μ C S + β S I + b β C S I + + ( μ + γ ) I +
Since μ + γ = β S + + b β C S + , μ = a θ S + b β C S + I + C S + and Π = μ S + + β S + I + + a θ S + = μ S + + μ C S + + ( μ + γ ) I + , this derivative becomes
L + ˙ = 2 Π μ S a θ S + b β C S + I + C S + C S ( μ S + + β I + S + + a θ S + ) S + S + a θ S + a θ S C S + C S β S I + b β C S I + = 2 ( μ S + + β S + I + + a θ S + ) μ S a θ S + C S C S + + a θ S + μ S + 2 S β S + 2 I + S a θ S + 2 S a θ S C S + C S β S I + = ( μ S + + β S + I + ) S S + + S + S 2 a θ S + C S C S + + S + S + C S + S C S S + 3 0
which is the final result, sought to prove for deducing that E + is globally asymptotically stable. □

3. Free Horizon Isoperimetric Optimal Control Approach

Now, we consider the mathematical model in Equation (1) with θ as a control function of time t.
Motivated by the desire to find the optimal time needed to reduce the number of infected people as much as possible while minimizing the value of the control θ ( t ) over a free (non-fixed) horizon t f , our objective is to seek a couple ( θ ( t ) , t f ) such that
J ( θ ( t ) , t f ) = min ( θ ( t ) , t f ) U × R + J ( θ ( t ) , t f )
where J is the functional defined by
J ( θ ( t ) , t f ) = t f 2 + 0 t f a I ( t ) + b 2 θ 2 ( t ) d t
and where the control space U is defined by the set
U = { θ ( t ) | 0 θ ( t ) 1 , θ ( t ) m e a s u r a b l e , t [ 0 , t f ] , t f f r e e }
where a and b represent constant severity weights associated to functions I and θ , respectively. Alkama et al. treated three cases of the form of the free horizon t f in the final gain function of their objective function when applying a free final time optimal control approach to a cancer model [9]. Here, we suppose that t f takes the quadratic form as formulated in Equation (6) to obtain a direct formula which characterizes t f . In fact, if t f is taken linear or the final gain function is zero, t f would just be approximated numerically due to the nature of necessary conditions in these two cases (see [9] for explanation).
Since managers of the anti-epidemic resources cannot well-predict whether their control strategy will reach the entire susceptible population over a fixed horizon T, we treat here an example where the number of targeted people in the susceptible class is equal for example to only a constant C = 3026 for T = 50 months. Hence, we try to find ( θ ( t ) , t f ) under the definition of the following isoperimetric restriction
0 T a θ ( t ) S ( t ) d t = C
In [13], the authors defined an isoperimetric constraint on the control variable only to model the total tolerable dosage amount of a therapy along the treatment period. In their conferences talks [26,27], Kornienko et al. and De Pinho et al. introduced state constraints in an optimal control problem that is subject to an S-Exposed-I-R differential system to model the situation of limited supply of vaccine based on work in [14] and where the isoperimetric constraint is defined on the product of the control and state variables.
In our case, to take into account the constraint in Equation (7) for the resolution of the optimal control problem in Equation (5), we consider a new variable Z defined as
Z ( t ) = 0 t a θ ( v ) S ( v ) d v
Then, we have Z ˙ ( t ) = a θ ( t ) S ( t ) . Using notations of the state variables in the previous section and keeping θ as a notation of θ ( t ) and Z in place of Z ( t ) , we study the differential system defined as follows
S ˙ = Π β S I a θ S μ S C S ˙ = a θ S b β C S I μ C S I ˙ = β ( S + b C S ) I γ I μ I R ˙ = γ I μ R Z ˙ = a θ S
instead of the model in Equation (1). We also note that, when the minimization problem in Equation (5) is under the constraint in Equation (7), the application of Pontryagin’s Maximum Principle would not be appropriate for this case, but the new variable Z has the advantage to convert Equations (5)–(7) to a classical optimal control problem under one restriction which is the system in Equation (9) only [28]. If we follow most optimal control approaches in the literature, the objective function in Equation (6) will be defined over a fixed time interval. However, t f is free here, and, to find the optimal duration needed to control an epidemic, it would be advantageous to managers of medical or health resources to control an epidemic before reaching the fixed time T for lesser costs. For this purpose, we need to assume that 0 t f T to guarantee the sufficient condition for an optimal θ in the case of a free horizon. This is because θ exists for the minimization problem in Equation (5) when Equation (6) is defined over T based on the verified properties of the sufficient conditions as stated in details in Theorem 4.1, pp. 68–69 of [29] and that can easily be verified for many examples as ours, and this implies in our case that the existence of an optimal control θ and associated optimal trajectories S , C S , I , R and Z comes directly from the convexity of the integrand term in Equation (6) with respect to the control θ and the Lipschitz properties of the state system with respect to state variables S , C S , I , R and Z. Then, it exists for any time in the interval [ 0 , T ] including t f . As regard the necessary conditions, we state and prove the following theorem.
Theorem 4.
If there exist optimal control u and optimal horizon t f which minimize Equation (6) along with the optimal solutions S , C S , I and R associated to the differential system in Equation (9), then there exist adjoint variables λ k , k = 1 , 2 , 3 , 4 , 5 as notations of λ k ( t ) and which satisfy the following adjoint differential system
λ 1 ˙ = λ 1 ( β I + μ + a θ ) a λ 2 θ β λ 3 I a θ λ 5 λ 2 ˙ = λ 2 ( b β I + μ ) b λ 3 β I λ 3 ˙ = a + λ 1 β S + b λ 2 β C S λ 3 ( β ( S + b C S ) μ γ ) λ 4 γ λ 4 ˙ = λ 4 μ λ 5 ˙ = 0
with the transversality conditions λ k ( t f ) = 0 , k = 1 , 2 , 3 , 4 and λ 5 ( t f ) = constant which should be determined.
Furthermore, the sought optimal control is characterized by
θ = min max 0 , a S ( λ 1 λ 2 λ 5 ) b , 1
while the sought optimal horizon is characterized by
t f = H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) 2
where H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) defines the Hamiltonian function as the sum of the integrand term of Equation (6) and the term λ 1 S ˙ + λ 2 C S ˙ + λ 3 I ˙ + λ 4 R ˙ + λ 5 Z ˙ at t f .
Moreover, t f is positive only when H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) is negative.
Proof. 
Let H be a notation of the Hamiltonian function H ( t , S ( t ) , C S ( t ) , I ( t ) , R ( t ) , Z ( t ) , λ 1 ( t ) , λ 2 ( t ) , λ 3 ( t ) , λ 4 ( t ) , λ 5 ( t ) , θ ( t ) ) in all time t. Then, we have
H = a I + b 2 θ 2 + λ 1 S ˙ + λ 2 C S ˙ + λ 3 I ˙ + λ 4 R ˙ + λ 5 Z ˙ = a I + b 2 θ 2 + λ 1 ( Π β S I a θ S μ S ) + λ 2 ( a θ S b β C S I μ C S ) + λ 3 ( β ( S + b C S ) I γ I μ I ) + λ 4 ( γ I μ R ) + a θ S λ 5
Using Pontryagin’s maximum principle [30], we have
λ 1 ˙ = H S = λ 1 ( β I + μ + a θ ) a λ 2 θ β λ 3 I a θ λ 5 λ 2 ˙ = H C S = λ 2 ( b β I + μ ) b λ 3 β I λ 3 ˙ = H I = a + λ 1 β S + b λ 2 β C S λ 3 ( β ( S + b C S ) μ γ ) λ 4 γ λ 4 ˙ = H R = λ 4 μ λ 5 ˙ = H R = 0
while the transversality conditions defined as minus the derivative of the final gain function with respect to the state variables S, C S , I and R. Since the final gain function in Equation (6) does not contain any term of these variables, then λ k ( t f ) = 0 , k = 1 , 2 , 3 , 4 and λ 5 ( t f ) is unknown but we are sure it is a constant since λ ˙ 5 ( t ) = 0 t [ 0 , t f ] . The solution of this problem is treated in the next section.
The optimality condition at θ = θ implies that H θ = 0 . Then, after setting S = S , we have
b θ a S λ 1 + a S λ 2 + a S λ 5 = 0 θ = a S ( λ 1 λ 2 λ 5 ) b
Taking into account the bounds of the control, we obtain,
θ = min max 0 , a S ( λ 1 λ 2 λ 5 ) b , 1
Now, let us prove the necessary conditions on t f . As J ( θ , t f ) reaches its minimum at θ and t f , we have
lim h 0 J θ , t f + h J θ , t f h = 0
with the consideration of the final gain function ϕ that we deduce it is defined in Equation (6) by ϕ ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) ) = t f 2 , while setting θ = θ and t f = t f , we obtain
lim h 0 1 h ϕ ( t f + h , S ( t f + h ) , C S ( t f + h ) , I ( t f + h ) , R ( t f + h ) , Z ( t f + h ) ) + 0 t f + h a I ( t ) + b 2 θ 2 ( t ) d t ϕ ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) ) 0 t f k I a I ( t ) + b 2 θ 2 ( t ) d t = 0 lim h 0 ϕ ( t f + h , S ( t f + h ) , C S ( t f + h ) , I ( t f + h ) , R ( t f + h ) , Z ( t f + h ) ) ϕ ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) ) h + 1 h t f t f + h a I ( t ) + b 2 θ 2 ( t ) d t = 0 ϕ t ( t f ) + ϕ S ( t f ) S ˙ ( t f ) + ϕ C S ( t f ) C S ˙ ( t f ) + ϕ I ( t f ) I ˙ ( t f ) + ϕ R ( t f ) R ˙ ( t f ) + ϕ Z ( t f ) Z ˙ ( t f ) + a I ( t ) + b 2 θ 2 ( t ) = 0 2 t f + H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) = 0 t f + H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) 2 = 0
Finally, we have
t f = H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) ) 2
Otherwise, the positivity of t f under the condition of negativity of
H ( t f , S ( t f ) , C S ( t f ) , I ( t f ) , R ( t f ) , Z ( t f ) , λ 1 ( t f ) , λ 2 ( t f ) , λ 3 ( t f ) , λ 4 ( t f ) , λ 5 ( t f ) , θ ( t f ) )
is trivial, but this is not a condition we should have necessarily for θ since the Hamiltonian could change signs any time along the interval of study. □

4. Numerical Simulations

Based on the formulation of Equation (8), we have Z ( 0 ) = 0 and Z ( t f ) = C . Since the optimal control problem consists to resolve the two-point boundary value problem defined by the two systems in Equations (2) and (10), the differential system in Equation (2) will be numerically resolved forward in time because of its initial conditions and the value of Z(0) does not change, while the differential system in Equation (10) will be numerically resolved backward in time because of its final or transversality conditions but with the condition that Z ( t f ) varies depending on the value of k. Based on the numerical approach in [13], we propose also here to define a real function g such that k g ( k ) = Z ˜ f Z f and where Z ˜ f is the value of Z at t f for various values of k and Z f is the value fixed by C. This leads to the combination of the Forward–Backward–Sweep Method (FBSM) which resolves the two-point boundary value problem in Equations (2)–(10), with the secant-method to find the value of the root k of the function g [31]. The necessary condition on t f defined by the characterization in Equation (12), which leads to seek a fixed point of a real function F such that F ( t f ) = t f . We choose to solve this numerical problem differently to the method used in [8,9] using the fixed point method. In brief, the four steps of numerical calculus associated to the resolution of our free optimal control problem (5) under isoperimetric constraint (7), are described in Algorithm 1.
Algorithm 1: Resolution steps of the two-point boundary value optimal control problem (9) and (10).
Step 0:
     Guess an initial estimation to θ and t f i n a l .
Step 1:
     Use the initial condition S ( 0 ) , C S ( 0 ) , I ( 0 ) , R ( 0 ) and Z ( 0 ) and the stocked values by θ and t f .
     Find the optimal states S , C S , I , R and Z which iterate forward in the two-point boundary value problem (2)–(10).
Step 2:
     Use the stocked values by θ and the transversality conditions λ k ( t f ) for k = 1 , 2 , 3 , 4 while searching the constant λ 5 ( t f ) using the secant-method.
     Find the adjoint variables λ k for k = 1 , 2 , 3 , 4 , 5 which iterate backward in the two-point boundary value problem (2)–(10).
Step 3:
     Update the control utilizing new S, C S , I, R, Z and λ k for k = 1 , 2 , 3 , 4 , 5 in the characterization of θ as presented in (11) while searching the optimal time t f characterized by (12) using the fixed point method.
Step 4:
     Test the convergence. If the values of the sought variables in this iteration and the final iteration are sufficiently small, check out the recent values as solutions. If the values are not small, go back to Step 1.
Figure 1 depicts the S C S I R Z dynamics in the absence and presence of the control and we can see that the number of susceptible people has increased linearly from its initial condition to a number higher than 92.5 individuals when we choose θ = 0 , while the optimal state S increases during the first months of the optimal control strategy and it decreases when we work with the characterization of Equation (11). Simultaneously, the number of removed people increases to only a value close to eight people while it reaches a value higher than this number with a maximal peak equaling to 17 when θ 0 . As regards to the number of infected people, it decreases from its initial condition to a value close to an important value of 30 individuals because of the natural death and recovery only, while it decreases towards a value very close to zero after the introduction of the control θ . We can see the relationship between the number of controlled people and the optimal values taken by θ so when this is increasing, the optimal state C S is also increasing. In fact, we can deduce that, with only small values of θ , we reach our goal by minimizing I function, and maximizing R function while the total number of the susceptible who received the control along T and which is represented by the function Z has not exceeded the imposed constant C. The dashed lines introduced in this figure show the highest fixed point value of the sought final time, and we can understand that, at this point, we have already reached our goal which concerns the minimization of the number of infected people and maximization of the number of removed people. The next figure gives more information about the obtained value of t f .
In Figure 2, we present dynamics of the functions S, I and R, and we can see the fixed points t f in the first plot above. The solution of the equation F ( t f ) = t f starts from an initial guess which equals zero, and increases to values that are very close or sometimes equal to 26 months (we note that, even if they appeared taking the value 26, this is not the case at all iterations but just because all values are very close to 26 with a small precision of about 10 4 ). As noted in this figure, for instance, the highest value of t f = 26.4081 found at iteration 292 among 1000. In the same figure, in the plot below, we observe that, at t f indicated by the dashed purple line, the number of infected people has already taken the direction towards zero values, while the number of removed people has already reached its positive peak and started to decrease because of the decrease of the optimal control function θ , as shown in the previous figure. This means that there is no need in this case to extend the optimal control approach for other months since, at t f , Equation (5) has been almost realized.
In Figure 3, the fixed points t f for different values of the control severity weight b suggest that, as the value of b increases, t f increases. In fact, the bigger is b , the lesser is the optimal control θ , which is important, as we can deduce from the formulation in Equation (11), and this is reasonable since, when θ is small, we need more time to control the epidemic. The obtained results in this figure can be summarized as follows:
  • When b = 60 : I ( t f ) = 6.734 with θ ( t f ) = 0.8658 (iteration 278), which implies that 83.165 % of infected people have left the I compartment.
  • When b = 70 : I ( t f ) = 6.0627 with θ ( t f ) = 0.8624 (iteration 295), which implies that 84.8425 % of infected people have left the I compartment.
  • When b = 80 : I ( t f ) = 4.7619 with θ ( t f ) = 0.8573 (iteration 332), which implies that 88.09525 % of infected people have left the I compartment.
In Figure 4, we show the impact of the initial condition of I function, namely I 0 , on fixed points t f , and we can deduce from the obtained optimal horizons that, as I 0 increases, t f increases, and this is reasonable since, when the number of infected people is important, the anti-epidemic measures need longer time for controlling the situation. The obtained results in this figure can be summarized as follows:
  • When I ( 0 ) = 50 : I ( t f ) = 8.4416 with θ ( t f ) = 0.8758 (iteration 278), which implies that 83.1168 % of infected people have left the I compartment.
  • When I ( 0 ) = 60 : I ( t f ) = 5.7018 with θ ( t f ) = 0.8662 (iteration 338), which implies that 90.497 % of infected people have left the I compartment.
  • When I ( 0 ) = 70 : I ( t f ) = 8.1909 with θ ( t f ) = 0.8802 (iteration 310), which implies that 88.2987 % of infected people have left the I compartment.

5. Conclusions

In this paper, we have determined the optimal duration needed for controlling an epidemic based on a free horizon optimal control approach with an isoperimetric constraint and which has been applied to a four-compartmental epidemic model where it is supposed that the controlled population does not reach the removed class due to the temporary effect of the control. The isoperimetric restriction which has been proposed to define the number of susceptible people who receive the control along the anti-epidemic measures period, allowed us to find the optimal horizon of the optimal control strategy when there are limited resources devised to fight against a disease. In the numerical simulations, we used the fixed point method since the necessary condition on the free horizon led to a fixed point equation. Our results prove their usefulness, since, at the obtained optimal horizons for different values of parameters and initial conditions on infection, the infected population size has been reduced and this presents an advantage of the followed control approach to managers of the health resources even when these are limited.

Author Contributions

All authors contributed equally to this work. All authors read and approved the final manuscript.

Acknowledgments

The authors would like to thank all Managing Editors and members of the Editorial Board who were responsible for dealing with this paper, and the anonymous referees for their valuable comments and suggestions, improving the content of this paper.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zakary, O.; Rachik, M.; Elmouki, I. On the impact of awareness programs in HIV/AIDS prevention: An SIR model with optimal control. Int. J. Comput. Appl. 2016, 133, 1–6. [Google Scholar] [CrossRef]
  2. Zakary, O.; Larrache, A.; Rachik, M.; Elmouki, I. Effect of awareness programs and travel-blocking operations in the control of HIV/AIDS outbreaks: A multi-domains SIR model. Adv. Differ Equ. 2016, 2016, 169. [Google Scholar] [CrossRef]
  3. Roy, P.K.; Saha, S.; Al Basir, F. Effect of awareness programs in controlling the disease HIV/AIDS: An optimal control theoretic approach. Adv. Differ. Equ. 2015, 2015, 217. [Google Scholar] [CrossRef]
  4. Rodrigues, H.S.; Monteiro, M.T.T.; Torres, D.F. Vaccination models and optimal control strategies to dengue. Math. Biosci. 2014, 247, 1–12. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Kumar, A.; Srivastava, P.K. Vaccination and treatment as control interventions in an infectious disease model with their cost optimization. Commun. Nonlinear Sci. Numer. Simul. 2017, 44, 334–343. [Google Scholar] [CrossRef]
  6. Liu, X.; Takeuchi, Y.; Iwami, S. SVIR epidemic models with vaccination strategies. J. Theor. Biol. 2008, 253, 1–11. [Google Scholar] [CrossRef] [PubMed]
  7. Nainggolan, J.; Supian, S.; Supriatna, A.K.; Anggriani, N. Mathematical model of tuberculosis transmission with reccurent infection and vaccination. J. Phys. Conf. Ser. 2013, 423, 012059. [Google Scholar] [CrossRef]
  8. Zakary, O.; Rachik, M.; Elmouki, I. How much time is sufficient for benefiting of awareness programs in epidemics prevention? A free final time optimal control approach. Int. J. Adv. Appl. Math. Mech. 2017, 4, 26–40. [Google Scholar]
  9. Alkama, M.; Larrache, A.; Rachik, M.; Elmouki, I. Optimal duration and dosage of BCG intravesical immunotherapy: A free final time optimal control approach. Math. Methods Appl. Sci. 2018, 41, 2209–2219. [Google Scholar] [CrossRef]
  10. Zhou, L.; Fan, M. Dynamics of an SIR epidemic model with limited medical resources revisited. Nonlinear Anal. Real World Appl. 2012, 13, 312–324. [Google Scholar] [CrossRef]
  11. Abdelrazec, A.; Bélair, J.; Shan, C.; Zhu, H. Modeling the spread and control of dengue with limited public health resources. Math. Biosci. 2016, 271, 136–145. [Google Scholar] [CrossRef] [PubMed]
  12. Yu, T.; Cao, D.; Liu, S. Epidemic model with group mixing: Stability and optimal control based on limited vaccination resources. Commun. Nonlinear Sci. Numer. Simul. 2018, 61, 54–70. [Google Scholar] [CrossRef]
  13. Elmouki, I.; Saadi, S. BCG immunotherapy optimization on an isoperimetric optimal control problem for the treatment of superficial bladder cancer. Int. J. Dyn. Control 2016, 4, 339–345. [Google Scholar] [CrossRef]
  14. Neilan, R.M.; Lenhart, S. An Introduction to Optimal Control with an Application in Disease Modeling. In Modeling Paradigms and Analysis of Disease Trasmission Models; Gumel, A.B., Lenhart, S., Eds.; American Mathematical Society: Providence, RI, USA, 2010; pp. 67–82. [Google Scholar]
  15. El Kihal, F.; Abouelkheir, I.; Rachik, M.; Elmouki, I. Optimal Control and Computational Method for the Resolution of Isoperimetric Problem in a Discrete-Time SIRS System. Math. Comput. Appl. 2018, 23, 52. [Google Scholar] [CrossRef]
  16. Kermack, W.O.; McKendrick, A.G. A Contribution to the Mathematical Theory of Epidemics. Proc. R. Soc. A 1927, 115, 700–721. [Google Scholar] [CrossRef] [Green Version]
  17. Korobeinikov, A. Lyapunov functions and global properties for SEIR and SEIS epidemic models. Math. Med. Biol. J. IMA 2004, 21, 75–83. [Google Scholar] [CrossRef] [Green Version]
  18. Gonçalves, S.; Abramson, G.; Gomes, M.F. Oscillations in SIRS model with distributed delays. Eur. Phys. J. B 2011, 81, 363. [Google Scholar] [CrossRef]
  19. Gray, A.; Greenhalgh, D.; Mao, X.; Pan, J. The SIS epidemic model with Markovian switching. J. Math. Anal. Appl. 2012, 394, 496–516. [Google Scholar] [CrossRef] [Green Version]
  20. Kribs-Zaleta, C.M.; Velasco-Hernández, J.X. A simple vaccination model with multiple endemic states. Math. Biosci. 2000, 164, 183–201. [Google Scholar] [CrossRef] [Green Version]
  21. Kribs-Zaleta, C.M.; Martcheva, M. Vaccination strategies and backward bifurcation in an age-since-infection structured model. Math. Biosci. 2002, 177, 317–332. [Google Scholar] [CrossRef]
  22. Alexander, M.E.; Bowman, C.; Moghadas, S.M.; Summers, R.; Gumel, A.B.; Sahai, B.M. A vaccination model for transmission dynamics of influenza. SIAM J. Appl. Dyn. Syst. 2004, 3, 503–524. [Google Scholar] [CrossRef]
  23. Shim, E. A note on epidemic models with infective immigrants and vaccination. Math. Biosci. Eng. 2006, 3, 557. [Google Scholar] [CrossRef] [PubMed]
  24. Sharomi, O.; Malik, T. Optimal control in epidemiology. Ann. Oper. Res. 2017, 251, 55–71. [Google Scholar] [CrossRef]
  25. Hethcote, H.W. The mathematics of infectious diseases. SIAM Rev. 2000, 42, 599–653. [Google Scholar] [CrossRef]
  26. Kornienko, I.; Paiva, L.T.; De Pinho, M.D.R. Introducing state constraints in optimal control for health problems. Procedia Technol. 2014, 17, 415–422. [Google Scholar] [CrossRef]
  27. De Pinho, M.D.R.; Kornienko, I.; Maurer, H. Optimal Control of a SEIR Model with Mixed Constraints and L1 Cost. In Proceedings of the 11th Portuguese Conference on Automatic Control, Porto, Portugal, 21–23 July 2014; pp. 135–145. [Google Scholar]
  28. Lenhart, S.; Workman, J.T. Optimal Control Applied to Biological Models; CRC Press: New York, NY, USA, 2007. [Google Scholar]
  29. Fleming, W.H.; Rishel, R.W. Deterministic and Stochastic Optimal Control; Springer Science & Business Media: New York, NY, USA, 2012. [Google Scholar]
  30. Pontryagin, L.S. Mathematical Theory of Optimal Processes; CRC Press: New York, NY, USA, 1987. [Google Scholar]
  31. Gumel, A.B.; Lenhart, S. Modeling Paradigms and Analysis of Disease Transmission Models; American Mathematical Society: Providence, RI, USA, 2010. [Google Scholar]
Figure 1. S C S I R Z dynamics in the absence and presence of the control in the two cases θ = 0 and θ 0 . Parameters values: Π = 6.45 , a = 0.06 , b = 0.001 , β = 0.0003 , and μ = 0.05 , γ = 0.1 . Initial conditions: S ( 0 ) = 89 , C S ( 0 ) = 0 , I ( 0 ) = 40 , and R ( 0 ) = 0 . Severity weights constants: a = 1 and b = 50 .
Figure 1. S C S I R Z dynamics in the absence and presence of the control in the two cases θ = 0 and θ 0 . Parameters values: Π = 6.45 , a = 0.06 , b = 0.001 , β = 0.0003 , and μ = 0.05 , γ = 0.1 . Initial conditions: S ( 0 ) = 89 , C S ( 0 ) = 0 , I ( 0 ) = 40 , and R ( 0 ) = 0 . Severity weights constants: a = 1 and b = 50 .
Mca 23 00064 g001
Figure 2. SIR dynamics with the precision of the optimal horizon t f with the same parameters, initial conditions and severity weights constants as in Figure 1.
Figure 2. SIR dynamics with the precision of the optimal horizon t f with the same parameters, initial conditions and severity weights constants as in Figure 1.
Mca 23 00064 g002
Figure 3. t f for different values of b with the same parameters, initial conditions and severity weights constants as in Figure 1.
Figure 3. t f for different values of b with the same parameters, initial conditions and severity weights constants as in Figure 1.
Mca 23 00064 g003
Figure 4. t f for different values of the initial condition I ( 0 ) with the same parameters, severity weights constants and, initial conditions for S ( 0 ) , C S ( 0 ) , and R ( 0 ) as in Figure 1.
Figure 4. t f for different values of the initial condition I ( 0 ) with the same parameters, severity weights constants and, initial conditions for S ( 0 ) , C S ( 0 ) , and R ( 0 ) as in Figure 1.
Mca 23 00064 g004

Share and Cite

MDPI and ACS Style

Abouelkheir, I.; El Kihal, F.; Rachik, M.; Elmouki, I. Time Needed to Control an Epidemic with Restricted Resources in SIR Model with Short-Term Controlled Population: A Fixed Point Method for a Free Isoperimetric Optimal Control Problem. Math. Comput. Appl. 2018, 23, 64. https://doi.org/10.3390/mca23040064

AMA Style

Abouelkheir I, El Kihal F, Rachik M, Elmouki I. Time Needed to Control an Epidemic with Restricted Resources in SIR Model with Short-Term Controlled Population: A Fixed Point Method for a Free Isoperimetric Optimal Control Problem. Mathematical and Computational Applications. 2018; 23(4):64. https://doi.org/10.3390/mca23040064

Chicago/Turabian Style

Abouelkheir, Imane, Fadwa El Kihal, Mostafa Rachik, and Ilias Elmouki. 2018. "Time Needed to Control an Epidemic with Restricted Resources in SIR Model with Short-Term Controlled Population: A Fixed Point Method for a Free Isoperimetric Optimal Control Problem" Mathematical and Computational Applications 23, no. 4: 64. https://doi.org/10.3390/mca23040064

APA Style

Abouelkheir, I., El Kihal, F., Rachik, M., & Elmouki, I. (2018). Time Needed to Control an Epidemic with Restricted Resources in SIR Model with Short-Term Controlled Population: A Fixed Point Method for a Free Isoperimetric Optimal Control Problem. Mathematical and Computational Applications, 23(4), 64. https://doi.org/10.3390/mca23040064

Article Metrics

Back to TopTop