Optimal control of non-autonomous SEIRS models with vaccination and treatment

We study an optimal control problem for a non-autonomous SEIRS model with incidence given by a general function of the infective, the susceptible and the total population, and with vaccination and treatment as control variables. We prove existence and uniqueness results for our problem and, for the case of mass-action incidence, we present some simulation results designed to compare an autonomous and corresponding periodic model, as well as the controlled versus uncontrolled models.

1. Introduction. Frequently, decision makers must balance the effort put in treatment and vaccination in order to give the best response to outbreaks of infectious diseases. Optimal control is an important mathematical tool that can be used to find the best strategies [16,23,27]. Real-world problems, under recent investigation with optimal control techniques, include Ebola [2] and Zika [6].
In [10], Gaff and Schaefer studied three distinct autonomous epidemic models, having vaccination and treatment as control variables. They established existence and, in some small interval, uniqueness of solutions for the optimal control problems. In that paper, one of the main questions under study is to know if the underlying epidemic structure has a significant impact on the obtained optimal control strategy. One of the models discussed in [10] was the SEIR, which is one of the most studied models in epidemiology. In the family of SEIRS models, it is assumed that the population is divided in four compartments: additionally to the infected class I, the susceptible class S and the recovered class R, present in SIR models, an exposed class E is also considered, in order to divide the infected population in the group of individuals that are infected and can infect others (the infective class) and the individuals that are infected but are not yet able to infect other individuals (the exposed or latent class). Such division of the population is particularity suitable to include several infectious diseases, e.g. measles and, assuming vertical transmission, rubella [17]. When there is no recovery, the model can be also used to describe diseases such as Chagas' disease [29]. It is also appropriate to model hepatitis B and AIDS [17], and Ebola [23]. Although influenza can be modeled by a SEIRS model [5], due to the short latency period, it is sometimes better to use the simpler SIRS formulation [7].
In [10] the parameters of the considered models are assumed to be independent of time. This is not very realistic in many situations, in particular due to periodic seasonal fluctuations. A classical example of seasonal patterns of incidence, exhibited by some infectious diseases, is given by data on weekly measles notification in England and Wales during the period 1948-1968 [1]. Other examples occur in several childhood diseases such as mumps, chicken-pox, rubella and pertussis [19]. It is also worth to mention that some environmental and demographic effects can be non-periodic. For example, for some diseases like cholera and yellow fever, it is known that the size of the latency period may decrease with global warming [26]. In this work, we consider a controlled SEIRS model with vaccination and treatment as control variables but, unlike [10], we let the parameters of our model to be time dependent. One of our main objectives is to discuss the effect of seasonal behavior on the optimal strategy.
Another important aspect of our work is that we consider a model with a general incidence function given by some function ϕ of the susceptible, the infective and the total population. This allows us to prove simultaneously results on existence and uniqueness of optimal solution for models with several different incidence functions that have already been considered in the context of SEIR/SEIRS models. In particular, our setting includes not only Michaelis-Menten incidence functions, considered for instance in [4,11,15,21,30,33], but also incidence functions that are not bilinear, which are appropriate to include saturation effects as well as other non-linear behaviors [18,35]. Additionally, and in contrast with [10], we assume that immunity can be partial and thus a fraction of the recovered individuals return to the susceptible class.
The paper is organized in the following way. In Section 2, we describe our nonautonomous SEIRS model, including treatment and vaccination as control variables, and we formulate the optimal control problem under investigation. Then, in Section 3, we discuss the question of existence of optimal solutions. The optimal controls are characterized in Section 4 with the help of Pontryagin's minimum principle and, in Section 5, we present a result concerning uniqueness of the optimal control. We end with Section 6 of numerical simulations. 2. The non-autonomous SEIRS model and the optimal control problem. In practice, evolution on the number of susceptible, exposed, infective and recovered, depends on some factors that can be controlled. Two of the main factors are: treatment of infective and vaccination of susceptible. For this reason, we consider a non-autonomous SEIRS model including treatment, denoted by T, and vaccination, denoted by V, as control variables. Namely, we consider the optimal control problem where κ 1 , κ 2 , κ 3 and S 0 , E 0 , I 0 , R 0 are non-negative, the state variables are absolutely continuous functions, i.e., (S(·), E(·), I(·), R(·)) ∈ AC [0, t f ]; R 4 , and the controls are Lebesgue integrable, i.e., (T(·), V(·)) ∈ L 1 ([0, Moreover, the following conditions hold: C1) Λ, β, µ, η, ε, γ ∈ C 1 ([0, t f ]; R) are ω-periodic; C2) ϕ ∈ C 2 (R 3 ; R); C3) ϕ(0, y, x) = ϕ(x, y, 0) = 0 for all x ∈ R + 0 and y ∈ R + ; f (x, y, z) = ϕ(x, y, z)/z, if z > 0, and f (x, y, 0) = lim z→0 ϕ(x, y, z)/z, is continuous and bounded. The total population N (t) = S(t) + E(t) + I(t) + R(t), t ∈ [0, t f ], is not constant (see Remark 1). It can be also seen, from the equations that define our control system, that we assume treatment to be applied to infective individuals only, moving a fraction of them from the infected to the recovered compartment, and that vaccination is applied to susceptible individuals only, also moving a fraction of them to the recovered class. Moreover, note that in our problem (P), besides general nonautonomous parameters, we consider a general incidence function. As examples of such incidence functions, we mention the ones obtained by ϕ(S, N, I) = SI/(1+αI), considered for instance in [25,34], ϕ(S, N, I) = I p S q , considered in [12,14,18], and ϕ(S, N, I) = SI p /(1 + αI q ), considered in [13,24]. Naturally, one needs to specify the incidence function, the other parameters and functions of the model, in order to carry out simulations and investigate particular situations. This is done in Section 6, where we compare an autonomous problem with a corresponding periodic model, as well as uncontrolled and corresponding controlled models. Before that, we prove results on existence of optimal solution (Theorem 3.2), necessary optimality conditions (Theorem 4.1), and uniqueness of solution (Theorem 5.1) for (P).
3. Existence of optimal solutions. Problem (P) is an optimal control problem in Lagrange form: In the above context, we say that a pair (x, u) ∈ AC ([t 0 , t 1 ]; R n ) × L 1 ([t 0 , t 1 ]; U ) is feasible if it satisfies the Cauchy problem in (1). We denote the set of all feasible pairs by F. Next, we recall a theorem about existence of solution for optimal control problems (1), contained in Theorem III.4.1 and Corollary III.4.1 in [9].
Theorem 3.1 (See [9]). For problem (1), suppose that f and L are continuous and there exist positive constants C 1 and C 2 such that, for t ∈ R, x, x)u, and L(t, x, ·) is convex on U ; g) L(t, x, u) ≥ c 1 |u| β − c 2 , for some c 1 > 0 and β > 1. Then, there exist (x * , u * ) minimizing J on F. Before applying Theorem 3.1 to obtain an existence result to our problem (P), we make a useful remark. Remark 1. Adding the equations in (P), we obtain that N (t) = Λ(t) − µ(t)N (t), which describes the behavior of the total population N (t). Since N 0 = N (t 0 ) = S 0 + E 0 + I 0 + R 0 , we have Consider the problem obtained by replacing function ϕ in (P) by some bounded and twice continuously differentiable function ψ such that ψ(S, N, I) = ϕ(S, N, I) for (S, N, I) ∈ [0, K] 3 (and maintaining the initial condition, the cost functional and the set of admissible controls). By (2), this new problem has exactly the same solutions as problem (P).
Theorem 3.2 (Existence of solutions for the optimal control problem (P)). There exists an optimal control pair (T * , V * ) and a corresponding quadruple (S * , E * , I * , R * ) of the initial value problem in (P) that minimizes the cost functional J in (P) Proof. We show that the problem obtained from (P) by replacing the function ϕ by some of the functions ψ in Remark 1 satisfies the conditions in Theorem 3.1. By Remark 1, using conditions C1) and C3), we immediately obtain a) and b). Conditions c) and d) are immediate from the definition of F and since U = [0, τ max ] × [0, ν max ]. By (2), we conclude that all state variables are in the compact set {(x, y, z, w) ∈ (R + 0 ) 4 : 0 ≤ x + y + z + w ≤ K} and condition e) follows. Since the state equations are linearly dependent on the controls, we obtain f). Finally, L is convex in the controls since it is quadratic.
and we establish g) with c 1 = min{k 2 , k 3 }. Thus, taking into account Remark 1, the result follows from Theorem 3.1.

4.
Characterization of the optimal controls. Now we address the question of how to identify the solutions predicted by Theorem 3.2. We do this with the help of the celebrated Pontryagin Maximum Principle [22]. This is possible because all control minimizers T * and V * of problem (P) are in fact essentially bounded. Indeed, Theorem 3.1 requires U to be a closed set, not necessarily bounded, and it may happen, in general, that the optimal controls predicted by Theorem 3.1 are in L 1 but not in L ∞ and do not satisfy the Pontryagin Maximum Principle [28]. However, in our case, U is compact and we can conclude that the L 1 optimal controls (T * , V * ) of problem (P), predicted by Theorem 3.2, are in fact in L ∞ , as required by the necessary optimality conditions [22]. Moreover, our optimal control problem (P) has only given initial conditions, with the state variables being free at the terminal time, that is, This implies that abnormal minimizers [3] are not possible in our context and we can fix the cost multiplier associated with the Lagrangian L to be one: for our optimal control problem (P), the Hamiltonian is given by In what follows, we use the operator ∂ i to denote the partial derivative with respect to the ith variable.
Theorem 4.1 (Necessary optimality conditions for the optimal control problem (P)). If ((S * , E * , I * , R * ), (T * , V * )) is a minimizer of problem (P), then there exist for almost all t ∈ [0, t f ], with transversality conditions Furthermore, the optimal control pair is given by Proof. Direct computations show that equations (3) are a consequence of the adjoint system of the Pontryagin Minimum Principle (PMP) [22]. Similarly, (4) are directly given by the transversality conditions of the PMP. It remains to characterize the controls using the minimality condition of the PMP [22]. The minimality condition and thus, on this set, If t ∈ int{t ∈ [t 0 , t 1 ] : V * (t) = 0}, then the minimality condition is Analogously, if t ∈ int{t ∈ [t 0 , t 1 ] : T * (t) = 0}, then the minimality condition is Therefore, we obtain (5) and (6).

5.
Uniqueness of the optimal control. In this section we show that the optimal solution of (P) is unique. The result is new even in the particular case where all the parameters of the model are time-invariant (autonomous case), extending the local uniqueness results in [8,10], which are valid only in a small time interval, to uniqueness in a global sense, that is, uniqueness of solution of (P) along all the time interval [0, t f ] where the optimal control problem is defined. The proof of Theorem 5.1 is nontrivial, lengthy and technical. It can be summarized as follows: 1. By contradiction, we assume that there are two distinct optimal pairs of state and co-state variables ξ = (S, E, I, R, p 1 , p 2 , p 3 , p 4 ) and ξ * = (S * , E * , I * , R * , contains the time range of the optimal control problem (P), then we are done. Otherwise, taking for initial conditions at time T the values of the state trajectories at the right-end of the interval [0, T ], we obtain uniqueness for the interval [T, 2T ], because the estimates that allow us to obtain T are only related with the maximum value of the parameters and the bounds for the state and co-state variables on the invariant region Γ of the new control problem, and are therefore the same as the ones already obtained for [0, T ] 5. Iterating the procedure, we obtain uniqueness throughout the interval [0, t f ] after a finite number of steps.
Theorem 5.1 (Uniqueness of solution for the optimal control problem (P)). The solution of the optimal control problem (P) is unique.
Analogously, we obtain where Finally, we have all the bounds needed to prove our result. Adding equations (13), (14), (15), (16), (20), (21), (22) and (23), we obtain, for the sum of the left-hand sides, which is equivalent to We now choose α so that α > C + C and note that α− C C > 1. Subsequently, we choose T such that Then, It follows that α − C − Ce 3αT > 0, so that inequality (26) can hold if and only if we have and p 4 (t) =p 4 (t). With this, the uniqueness of the optimal control is established on the interval [0, T ]. If T ≥ t f , then we have uniqueness on the whole interval. Otherwise, we can obtain uniqueness on [T, 2T ] for the optimal control problem whose initial conditions on time T coincide with the values of S, E, I and R on the end-time of the interval [0, T ] (note that, by the forward invariance of the set Γ, and since the constants C and C in (27) depend only on the values of the several state and co-state variables on Γ, we still have the same T ). Proceeding in the same way, we conclude, after a finite number of steps, that we have uniqueness on the interval [0, t f ]. The proof is complete.
6. Numerical simulations. In what follows, the incidence into the exposed class of susceptible individuals and the birth function Λ (t) are chosen, for illustrative purposes, to be β(t)ϕ(S, N, I) = 0.56(1 − per · cos(2πt + 0.26))SI (28) and Λ (t) = 0.05 + 0.05 per · cos(2πt) (29) with per ∈ [0, 1[. The remaining parameter functions, µ (t), η (t), ε (t) and γ (t), are assumed constant. Their values, as well as several other used in this section, were taken from [31] and [32] and are presented in Table 1. Note that we are in the autonomous case when per = 0 and in a periodic situation for 0 < per < 1. The solutions to the optimal control problems (P) here considered exist (Theorem 3.2), are unique (Theorem 5.1) and can be found using the optimality conditions given in Theorem 4.1. Precisely, our method to solve the problem consists to use the state equations (the four ordinary differential equations of problem (P)), the initial conditions, the adjoint equations (3) and the transversality conditions (4) with the controls given by (5) and (6), which are substituted into the state and adjoint equations. The state equations are solved with the initial conditions of Table 1, while the adjoint system is solved backward in time, with the change of variable Then, the two systems of equations constitute an initial value problem, which is solved numerically with an explicit 4th and 5th order Runge-Kutta method through the ode45 solver of MATLAB. The procedure is briefly described in Algorithm 1. 5. If the relative error is smaller than a given tolerance for all the variables (< 1% in our case), then stop. Otherwise, go to step 2.
In each plot of Figures 1 and 2, we present two sets of trajectories (distinguished by using dashed and continuous lines), in order to easily compare the uncontrolled and optimally controlled situations, the former mentioned by the suffix "-u" in the figures' legend. The behavior of our SEIRS model with both per = 0 (autonomous case, Figure 1) and per = 0.8 (periodic case, Figure 2), show the effectiveness of optimal control theory. Indeed, in Figures 1 and 2 we observe that, if we apply treatment and vaccination as given by Theorem 4.1 (optimally controlled case), then the number of exposed and infected individuals is significantly lower, as well as the number of susceptible individuals, while the number of recovered is significantly higher. It can be also seen that the susceptible and recovered classes have a very different behavior in the controlled and uncontrolled situations.
In Figures 3 and 4, we have the same trajectories as in Figures 1 and 2. They illustrate the effect of the periodicity of Λ(t) (29) and β(t) (28) in the classes S, E, I and R of individuals. The effect is perceptible in susceptible and exposed classes, since the periodic functions are present in these classes. From our results, we claim that the periodicity effect is "softened" in the transition between classes.  (28) and (29)): uncontrolled (dashed lines) versus optimally controlled (continuous lines).    In Figure 5, we show the optimal controls: treatment T * (t) of infective individuals (left side) and vaccination V * (t) of susceptible (right side). According to the minimality conditions (5) and (6)  takes place in the susceptible class. Treatment occurs in the infective class and, as we have seen, in this class the periodicity is not so perceptible. As a consequence, periodicity is only slightly perceptible in the treatment control variable T * (t).  In Figures 6 to 9, we present the behavior of infected individuals I * (t), treatment T * (t) and vaccination V * (t) optimal controls, when we vary the parameters µ, γ, ε and η, respectively, maintaining, in each case, the initial values and the remaining parameters as in Table 1. In all Figures 6-9, we varied the respective parameter (µ, γ, ε and η) from 0 to 0.1 in steps of length 0.01.
Referring to Figure 6, where the variation of µ is analyzed, we can see that the effect of periodicity is more perceptible in vaccination than in treatment. In the infected class, for low values of µ (low mortality), we observe that the number of infected increases. This is justified by the difference between birth and death.
Concerning Figure 7, where we vary the rate of recovery γ, the effect of periodicity is analogous to the previous situation: the bigger the value of γ, more infected individuals recover and thus the faster the infected class decreases.
In Figure 8, one can see the effect of the variation of the infectivity rate ε. The effect of periodicity is also in this case analogous to the previous situations. Indeed, when we have a higher value of ε, we have a faster transition of exposed individuals  into the infected class and this is the reason why we observe in Figure 8 that an increase on the value of ε leads to an increase in the infected class.
Finally, in Figure 9, the variation of the loss of immunity rate η is highlighted. We conclude that the periodicity effect is similar to the previous considered scenarios. However, the variation of η is the one that less influences the behavior of the three variables I * (t), T * (t) and V * (t).  It is worth mentioning that, in the simulations done and in the range of parameters considered, the variation of the periodic parameter per has a very small effect on the obtained cost functional J . Namely, we saw numerically that J (I, T, V) per=v1 − J (I, T, V) per=v2 ≤ 0.000329537, for v 1 , v 2 ∈ {0, 0.8}.