Abstract
For a stochastic differential equation SVIR epidemic model with vaccination, we prove almost sure exponential stability of the disease-free equilibrium for , where denotes the basic reproduction number of the underlying deterministic model. We study an optimal control problem for the stochastic model as well as for the underlying deterministic model. In order to solve the stochastic problem numerically, we use an approximation based on the solution of the deterministic model.
1. Introduction
The discovery of the first vaccine marked a major breakthrough in the battle against infectious diseases. The subsequent development of vaccines for various diseases has brought about remarkable results. Vaccination gained increasing popularity and success after it eradicated the smallpox outbreak of 1976 [1].
Vaccination is a commonly used method to control diseases such as measles, polio, and tuberculosis. Usually there are different schedules of dosage for different diseases and vaccines. For some diseases, doses should be taken by vaccinees several times and there must be some fixed time interval between two doses (see for instance Gabbuti et al. [2]). In a given population, the proportion of susceptibles who goes on to vaccination depends on different factors, one of which is the availability of the necessary resources.
There are numerous examples of vaccination models in the literature; see, for instance, the book of Brauer and Castillo-Chávez [3] or the journal papers [4–10]. With vaccination models one would be interested in the extent to which a vaccination program would reduce the basic reproduction number or how one can optimally roll out a vaccination program over time, in order to reach a certain target. In the latter case, optimal control theory is the obvious candidate to employ in the analysis [11–14]. It is usually the proportion of susceptibles who are admitted into vaccination, which is used as the control variable.
As a way of accommodating randomness in a compartmental epidemic model, several authors have proposed models with stochastic perturbation. This means that one modifies a system of ordinary differential equations (ODEs) by adding stochastic noise or stochastic perturbations, giving rise to a system of stochastic differential equations (SDEs) [7, 8, 15–18]. We shall refer to the original system of ODEs as the underlying deterministic model.
It is known (see the book of Mao [19], for instance) that the stability of a system can be improved by adding stochastic perturbations. SDE epidemic models have been studied by various authors until it was proved in research articles [9, 15, 20, 21] that stochastic perturbation improves stability of the disease-free equilibrium in the given models. SDE models with vaccination have been studied by Tornatore et al. [7, 8]. In particular in [8] it is shown that the model permits solutions that are almost surely positive, and a theorem on exponential stability in mean square stability is proved. The model of Tornatore et al. [8] will be further analyzed in this paper.
For SDE models in epidemiology, optimal control has not been studied (or at least not published) extensively. One of the reasons for this could very well be the difficulty with high dimensionality of the resulting partial differential equation (PDE) for the value function; see the paper of Sulem and Tapiero [22], for instance. A four-compartmental SIVR model such as in [7, 8] could easily lead to a PDE having the time variable together with three state variables. In this contribution we study the exponential stability of the stochastic model of Tornatore et al. [8] and control of vaccination both in the underlying deterministic model and in the stochastic model. In control problems, our goal is to characterize the vaccination rate (as control variable) which, on a finite time interval, minimizes the number of infected individuals balanced against the cost of vaccination. In the stochastic case we minimize an expected value. We use standard methods for the deterministic case and the Hamilton-Jacobi-Bellman equation in the stochastic case.
In Section 2 we study exponential stability of the disease-free equilibrium of the stochastic SIVR model that was introduced in [8]. We prove that almost sure exponential stability of the disease-free equilibrium prevails whenever , being the basic reproduction number of the underlying deterministic model. The deterministic control problem is treated in Section 3, and we find a numerical solution for the optimal control problem by using the fourth-order Runge-Kutta method. The optimal control of the SDE model is presented in Section 4. We observe a similarity in the form of the control in the two cases, stochastic versus deterministic. On this basis, in Section 5 we propose that the numerical solution of the control problem of the underlying deterministic model can be used to compute an approximate solution to the stochastic control problem, assuming the perturbation parameter to be small. We present computational examples to illustrate our findings.
2. The Stochastic Model and Stability
In this section we introduce the SVIR model and we analyse the stability of the disease-free equilibrium. Firstly, we formulate the necessary assumptions for modeling with SDEs. Let us assume having a filtered complete probability space and let be a one-dimensional Wiener process defined on this probability space.
We consider a stochastic SVIR model similar to the stochastic SIVR model of [8]. As in [8], the population is subdivided into four compartments/classes. These classes consist of all the individuals who are susceptible to the disease (), under vaccination (), infected with the disease (), and removed ().
It is assumed that births and deaths occur at the same constant rate and that all newborns enter the susceptible class. New infections occur at a rate , for a constant which is called the contact rate. A fraction of the susceptible class is being vaccinated at time . The vaccination may reduce but not completely eliminate susceptibility to the disease. Therefore the model includes a factor in the contact rate of vaccinated members, such that . If the vaccination is perfectly effective while means the vaccination has no effect at all. Furthermore, immunity to the disease is assumed to be permanent so that a fraction of infectives goes into the removed class.
The total population is constant and the variables are normalized so that The stochastic SVIR model is given as follows:Of course the parameters , , and are positive constants. The parameter determines the intensity of the stochastic perturbation. If , then the model reduces to the underlying deterministic model.
If is constant then the model above is identical to that in [8], and the pointis the disease-free equilibrium point of the underlying deterministic model and is the only equilibrium point of the stochastic model.
For some , some , and an -dimensional Wiener process , consider the general -dimensional stochastic differential equationA solution to the above equation is denoted by . We assume that for all , so that the origin point is an equilibrium of (4).
Definition 1. The equilibrium of the system equation (4) is said to be almost surely exponentially stable if, for all ,Let us denote by the differential operator associated with the function displayed in (4), defined for a function by Here means trace and denotes the transpose of a matrix.
For the remainder of this section we shall study stability and for this purpose we regard to be a positive constant, for all . In [8] it was shown that the system of SDEs has unique solutions which are almost surely positive. Thus we shall restrict ourselves to sample paths for which the coordinates are positive for all . The following observation is towards the proof of the stability theorem. Recall that we have introduced the number in expression (3). Let us define, for any constants and , the stochastic processes: Note that and . Now we calculateTherefore by the strong law of large numbers for martingales (see [19]) it follows that
Proposition 2. For constants and , let , and let . Then almost surely (a.s.) converges exponentially to if
Proof. Using the Itô formula and with and defined as above, we can express as We have shown that and hence the claim of the proposition follows readily.
The following observation from [9] is useful in exponential stability analysis.
Lemma 3. For , let be a bounded -valued function. Let be any increasing unbounded sequence of positive real numbers. Then there is a family of sequences such that, for each , is a subsequence of and the sequence converges to the largest limit point of the sequence .
Now we present our stability theorem. Recall that is the basic reproduction number of the underlying deterministic model, and . Also recall from expression (3) that .
Theorem 4. If , then the disease-free equilibrium point is almost surely exponentially stable.
Proof. The condition is equivalent to . Choose any number such that . Choose a number such thatNow let and . It suffices to prove that converges to exponentially (a.s.). To this end, by Proposition 2 it suffices to prove that Now we calculate . The latter can be expressed in the form with , , , and as follows: With reference to we note that since , we have and therefore we obtain the inequality: In view of Lemma 3 we can define the following limits for a suitable increasing, unbounded sequence : and with In particular then we have Let We find an upper bound for . Since we have since . Noting that and , we obtain Therefore satisfies the following inequality: Now we note that , and so we can write The coefficients of , , and are negative (see (13)). Furthermore, , , and cannot all be zero since Therefore , and the proof is complete.
We run two sets of simulations of the deterministic and stochastic and -trajectories over years. In the simulations, we use small perturbation parameters and small values of . All parameter values in the computations except and are the same in both scenarios and we also consider the same initial conditions. The common parameter values used in the computations are and while the initial conditions are , and .In the first run (see Figures 1 and 2) we use and , and we obtain . In the second run (see Figures 3 and 4) we choose and , which yields .
For , Theorem 4 guarantees the disease-free equilibrium to be almost surely exponentially stable and indeed the -curves in Figure 2 appear to converge to . For the theorem did not predict almost sure exponential stability and in fact the -curves in Figure 4 (and many simulations with not shown) certainly do not show any clear intention of converging to .
3. The Deterministic Optimal Control Problem
We now formulate and solve the deterministic version of the control problem. Recall from Section 2 that in our SVIR model, represents the fraction of the susceptible class being vaccinated at time . We wish to design an optimal vaccination schedule, , which minimizes a combination of the number of infectives on the one hand and the overall cost of the vaccination on the other hand, over a certain time horizon . The cost of the vaccination is assumed to be proportional to the square of .
For the purposes of optimization we introduce the functions , , and appearing in the SDE system (2) as follows:Now we can formulate the optimization problem.
Problem 5. Minimize the objective functionsubject to , , , , , .The control variable is assumed to be a measurable function of time and bounded: .
We solve Problem 5 using well-established control theory such as in the book [23] of Lenhart and Workman. We construct the Hamiltonian function, and to this end we introduce Lagrange multipliers , , and , also referred to as the costate variables. In the theorem below, the control variable, the state variables, and the costate variables are functions of time, but this dependence is suppressed except where required explicitly. The Hamiltonian of Problem 5 has the form
Theorem 6. An optimal solution for Problem 5 exists and satisfies the following system of differentiable equations:with transversality conditions , , and . Furthermore, the optimal vaccination rate is given by
Proof. In particular the Hamiltonian is convex with respect to and the existence of a solution follows. We check the first-order conditions for this optimization problem. We calculate the partial derivatives of the Hamiltonian with respect to the different state variables to obtain the time derivatives of the costate variables. The calculation of , , and leads to the equations asserted in the theorem. We now turn to the final part of the proof, which is about the form of the optimal control, . Since the function must optimize the Hamiltonian, we calculateConsider now a fixed value of . If is zero for some value in the interval , then the given value of is optimal. If for every we have (resp., ), then we must choose (resp., ). Thus we must have
In Figure 5 we plot for years the , , and -trajectories of the deterministic model subject to the optimal control . We simulate and for the deterministic model in Figure 6. In the simulations, we consider the parameters , , , , , and and initial conditions , , and .
In Figure 6 the dashed curve represents , that is, the optimal vaccination rate. In this illustration, the optimal vaccination roll-out starts with initial vaccination of approximately 0.11 and decreases gradually over the next 3 years.
4. The Stochastic Optimal Control Problem
In this section we formulate the stochastic version of the optimization problem and describe its solution. For stochastic control theory we refer to the book [24] of Øksendal. Our objective is to find an optimal vaccination rate that minimizes the objective functional which for an initial state is defined by Here the expectation is obtained on the condition that the initial state (at time ) of the system is . In step with the deterministic problem of earlier, we assume that there is a fixed constant with (a.s.). The class of admissible control laws isTo solve this stochastic control problem, we define the performance criterion as follows:where the expectation is conditional on the state of the system being a fixed value at time . We define the value function as We determine a control law that minimizes the expected value given by (37). We now formulate the stochastic analogue of the optimal control problem, subsequent to which we present the solution formulae.
Problem 7. Given the system (4) and given as in (36) with as in (37), find the value functionand an optimal control function
We can find an expression for the optimal vaccination strategy through the following theorem.
Theorem 8. A solution to the optimal vaccination problem stated in Problem (36) is of the form
Proof. We determine (41) via the dynamic programming approach. First we calculate : Applying the Hamilton-Jacobi-Bellman theory (see, for instance, [24]) we must find the infimum: For this purpose, we need to find the partial derivative of the expression with respect to , and this derivative should vanish. This leads to the equation:We consider the bounds on and by an argument similar to that in the proof of the deterministic case; the asserted expression for emerges.
5. Numerical Example
In the discussion that follows, the computations are done for years and we consider three different values of the perturbation parameter, . In particular we consider the -values: together with the following parameter values and initial conditions: , , , , , and , , , and .For each value of above, we use the results of the deterministic control problem to find an approximate numerical solution for the stochastic control problem. In particular, we use as a proxy for in the calculation of in this case. We note that the presence of makes into a stochastic variable even with the said proxy (in the stochastic case).
For each value of , we show in Table 1 the -values obtained as the average over 3000 runs made for different candidates for as we take . In these runs we take and of the underlying deterministic model.
Using a contact rate instead of , we choose the first candidate The adjustment on the contact rate is motivated by the general stabilizing effect of this type of stochastic perturbation and the perturbation being associated with the parameter .
We also take corresponding to the rate and with .
Furthermore, we consider corresponding to the contact rate .
Just for further comparison we choose a linear curve candidate , say, with and , also corresponding to a contact rate of .
Let us use the notation for the value of corresponding to the control . We notice that and are bigger than . Even more interesting is the pattern and . This creates the impression that is close to being a minimum. At least, with all the difficulties of a more precise solution, the choice of seems like a viable option in a real application. The inequalities and has been tested for other parameter choices too, with small. The simulations tested all revealed the same behaviour.
In Figures 7–9 we simulate the optimal solutions , , and of the deterministic and stochastic () models. We use the same values for the parameters , , , , , , , and initial conditions , , and as in the discussion preceding these simulations.
An important point to note about our approximation is that it fully accommodates the stochasticity (embodied and concentrated in the factor of the expression for ). Therefore it gives a very good approximation, at least in the sense that in Table 1, the minimum value of consistently corresponds to .
6. Conclusion
This paper presents some further insights into a stochastic vaccination model introduced in the paper of Tornatore et al. [8]. Our investigations covered two important aspects: exponential stability of the disease-free equilibrium and optimal control of vaccination.
Regarding stability, the main result of our paper has a particularly simple formulation. Essentially it says that we have almost sure exponential stability whenever the basic reproduction number of the underlying deterministic model is below unity. It will be good to know how an increase in vaccination rate would perhaps lead to better stability of the stochastic model. Nevertheless, the result as it stands is a good assurance.
As for the control problem, on the basis of a popular type of objective functional, we designed an efficient strategy for the roll-out of vaccination. It roughly amounts to minimizing the infections, balanced against the cost of vaccination. We exploit a similarity between the forms of the optimal controls for the stochastic model and the underlying deterministic model; then we use the relative simplicity of the latter to find approximate numerical solutions for the stochastic optimal control. Numerical simulation enables us to assess the feasibility of the option we followed, for a specific example. A more formal approach to the numerical solution of the optimal control problem is far more intricate and labour-intensive, and our method is a workable alternative. This could be the starting point for more sophisticated approximation methods.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The authors acknowledge financial support by the National Research Foundation (NRF) of South Africa.