# Vaccination Control in a Stochastic SVIR Epidemic Model

Authors: Peter J. Witbooi; Grant E. Muller; Garth J. Van Schalkwyk
Journal: Computational and Mathematical Methods in Medicine
 (2015)

---

## Abstract

For a stochastic differential equation SVIR epidemic model with vaccination, we prove almost sure exponential stability of the disease-free equilibrium forR
0
<
1, where R
0 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.


---

## Body

## 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 theunderlying 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 Section2 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 R
0
<
1, R
0 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(
Ω
,
F
,
F
t
t
≥
t
0
,
P
) and let W
(
t
) 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 (S), under vaccination (V), infected with the disease (I), and removed (R).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 β
S
I, for a constant β which is called thecontact rate. A fraction α
(
t
) of the susceptible class is being vaccinated at time t. 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 0
≤
ρ
≤
1. If ρ
=
0 the vaccination is perfectly effective while ρ
=
1 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(1)
S
t
+
V
t
+
I
t
+
R
t
=
1
∀
t
≥
0
.The stochastic SVIR model is given as follows:(2)
d
S
t
=
μ
-
μ
S
t
-
β
I
t
S
t
-
α
t
S
t
d
t
-
σ
S
t
I
t
d
W
t
,
d
V
t
=
α
t
S
t
-
ρ
β
V
t
I
t
-
μ
V
t
d
t
-
ρ
σ
V
t
I
t
d
W
t
,
d
I
t
=
β
S
t
I
t
+
ρ
β
V
t
I
t
-
γ
I
t
-
μ
I
t
d
t
+
σ
S
t
+
ρ
V
t
I
t
d
W
t
,
d
R
t
=
γ
I
t
-
μ
R
t
d
t
.Of course the parameters μ, γ, and β are positive constants. The parameter σ determines the intensity of the stochastic perturbation. If σ
=
0, then the model reduces to the underlying deterministic model.Ifα is constant then the model above is identical to that in [8], and the point(3)
E
0
=
S
0
,
I
0
,
V
0
,
R
0
=
m
,
0,1
-
m
,
0
with
m
=
μ
μ
+
αis the disease-free equilibrium point of the underlying deterministic model and is the only equilibrium point of the stochastic model.For somen
∈
N, some x
0
∈
R
n, and an n-dimensional Wiener process B
(
t
), consider the general n-dimensional stochastic differential equation(4)
d
x
t
=
F
x
t
,
t
d
t
+
G
x
t
,
t
d
B
t
,
x
0
=
x
0
.A solution to the above equation is denoted by x
(
t
,
x
0
). We assume that F
(
t
,
0
)
=
G
(
t
,
0
)
=
0 for all t
≥
0, so that the origin point is an equilibrium of (4).Definition 1.
The equilibriumx
=
0 of the system equation (4) is said to bealmost surely exponentially stable if, for all x
0
∈
R
n,(5)
lim
t
→
∞
⁡
sup
⁡
1
t
ln
⁡
x
t
,
x
0
<
0
a.s.Let us denote by L the differential operator associated with the function displayed in (4), defined for a function U
(
t
,
x
)
∈
C
1,2
(
R
×
R
n
) by (6)
L
U
=
∂
U
∂
t
+
F
trp
∂
U
∂
x
+
1
2
Trc
G
trp
∂
2
U
∂
x
2
G
.Here T
r
c means trace and t
r
p denotes the transpose of a matrix.For the remainder of this section we shall study stability and for this purpose we regardα
(
t
) to be a positive constant, α
(
t
)
=
α for all t
>
0. 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 w
∈
Ω for which the coordinates are positive for all t
≥
0. The following observation is towards the proof of the stability theorem. Recall that we have introduced the number m
=
μ
(
μ
+
α
)
-
1 in expression (3). Let us define, for any constants a
>
0 and c
>
0, the stochastic processes: (7)
z
t
=
a
S
t
-
m
2
+
I
t
+
c
R
t
,
u
1
t
=
-
2
a
S
t
-
m
σ
S
t
I
t
z
t
,
u
2
t
=
σ
S
t
+
ρ
V
t
I
t
z
t
.Note that I
(
t
)
/
z
(
t
)
≤
1 and S
(
t
)
+
ρ
V
(
t
)
≤
S
(
t
)
+
V
(
t
)
≤
1. Now we calculate(8)
∫
0
t
u
1
2
s
d
s
≤
∫
0
t
4
a
2
S
s
-
m
2
σ
2
S
2
s
d
s
≤
4
a
2
σ
2
t
,
∫
0
t
u
2
2
s
d
s
≤
∫
0
t
σ
2
S
s
+
ρ
V
s
2
d
s
≤
σ
2
t
.Therefore by the strong law of large numbers for martingales (see [19]) it follows that (9)
lim
t
→
∞
⁡
sup
⁡
1
t
∫
0
t
u
1
s
+
u
2
s
d
W
s
=
0
.Proposition 2.
For constantsa
>
0 and c
>
0, let z
(
t
)
=
a
(
S
(
t
)
-
m
)
2
+
I
(
t
)
+
c
R
(
t
), and let U
(
t
)
=
ln
⁡
z
(
t
). Then z
(
t
) almost surely (a.s.) converges exponentially to 0 if (10)
lim
t
→
∞
⁡
sup
⁡
L
U
t
<
0
a
.
s
.
.Proof.
Using the Itô formula and withu
1
(
t
) and u
2
(
t
) defined as above, we can express U
(
t
) as (11)
U
t
=
U
0
+
∫
0
t
L
U
s
d
s
+
∫
0
t
u
1
s
+
u
2
s
d
W
s
.We have shown that (12)
lim
t
→
∞
⁡
sup
⁡
1
t
∫
0
t
u
1
s
+
u
2
s
d
W
s
=
0
a.s.
,and hence the claim of the proposition follows readily.The following observation from [9] is useful in exponential stability analysis.Lemma 3.
Fork
∈
N, let X
(
t
)
=
(
X
1
(
t
)
,
X
2
(
t
)
,
…
,
X
k
(
t
)
) be a bounded R
k-valued function. Let (
t
0
,
n
) be any increasing unbounded sequence of positive real numbers. Then there is a family of sequences (
t
l
,
n
) such that, for each l
∈
{
1,2
,
…
,
k
}, (
t
l
,
n
) is a subsequence of (
t
l
-
1
,
n
) and the sequence X
l
(
t
l
,
n
) converges to the largest limit point of the sequence X
l
(
t
l
-
1
,
n
).Now we present our stability theorem. Recall thatR
0 is the basic reproduction number of the underlying deterministic model, and R
0
=
β
/
(
γ
+
μ
). Also recall from expression (3) that m
=
μ
(
μ
+
α
)
-
1.Theorem 4.
IfR
0
<
1, then the disease-free equilibrium point (
S
0
,
I
0
,
V
0
,
R
0
)
=
(
m
,
0,1
-
m
,
0
) is almost surely exponentially stable.Proof.
The conditionR
0
<
1 is equivalent to β
-
(
γ
+
μ
)
<
0. Choose any number c
>
0 such that β
-
(
γ
+
μ
)
+
γ
c
<
0. Choose a number a
>
0 such that(13)
β
-
γ
+
μ
+
γ
c
+
a
2
β
+
3
σ
2
<
0
.Now let z
(
t
)
=
a
(
S
(
t
)
-
m
)
2
+
I
(
t
)
+
c
R
(
t
) and U
(
t
)
=
ln
⁡
z
(
t
). It suffices to prove that z
(
t
) converges to 0 exponentially (a.s.). To this end, by Proposition 2 it suffices to prove that (14)
lim
t
→
∞
⁡
sup
⁡
L
U
t
<
0
a.s.
.Now we calculate L
U
(
t
). The latter can be expressed in the form (15)
L
U
t
=
A
t
+
B
1
t
+
B
2
t
+
B
3
t
,with A
(
t
), B
1
(
t
), B
2
(
t
), and B
3
(
t
) as follows: (16)
A
t
=
2
a
m
-
S
t
z
t
μ
+
α
S
t
-
m
+
β
I
t
S
t
+
I
t
z
t
β
S
t
+
ρ
β
V
t
-
γ
-
μ
+
c
z
t
γ
I
t
-
μ
R
t
,
B
1
t
=
2
a
z
t
-
4
a
2
S
t
-
m
2
2
z
2
t
σ
S
t
I
t
2
,
B
2
t
=
-
2
a
S
t
-
m
z
2
t
σ
S
t
+
ρ
V
t
I
t
σ
S
t
I
t
,
B
3
t
=
-
1
2
z
2
t
σ
S
t
+
ρ
V
t
I
t
2
.With reference to A
(
t
) we note that since ρ
≤
1, we have (17)
β
S
t
+
ρ
β
V
t
≤
β
S
t
+
ρ
V
t
≤
β
S
t
+
V
t
≤
β
,and therefore we obtain the inequality: (18)
A
t
≤
2
a
m
-
S
t
z
t
μ
+
α
S
t
-
m
+
β
I
t
S
t
+
I
t
z
t
β
-
γ
-
μ
+
c
z
t
γ
I
t
-
μ
R
t
.In view of Lemma 3 we can define the following limits for a suitable increasing, unbounded sequence (
t
n
): (19)
q
=
lim
n
→
∞
⁡
S
t
n
-
m
2
z
t
n
,
i
=
lim
n
→
∞
⁡
I
t
n
z
t
n
,
r
=
lim
n
→
∞
⁡
R
t
n
z
t
n
,
s
=
lim
n
→
∞
⁡
S
t
n
,
v
=
lim
n
→
∞
⁡
V
t
n
,and with (20)
limsup
t
→
∞
⁡
L
U
t
=
lim
n
→
∞
⁡
L
U
t
n
.In particular then we have (21)
a
q
+
i
+
c
r
=
1
,
0
≤
s
≤
1
,
0
≤
v
≤
1
.Let (22)
Λ
=
lim
n
→
∞
⁡
L
U
t
n
.We find an upper bound for B
1
(
t
)
+
B
2
(
t
)
+
B
3
(
t
). Since S
2
I
≤
1 we have (23)
B
1
t
≤
1
2
2
a
z
t
z
2
t
σ
S
t
I
t
2
≤
a
σ
2
I
t
z
t
,
B
2
t
≤
2
a
σ
2
I
t
z
t
2since |
S
(
t
)
-
m
|
<
1. Noting that B
3
(
t
)
≤
0 and i
<
1, we obtain (24)
lim
n
→
∞
⁡
B
1
t
n
+
B
2
t
n
+
B
3
t
n
≤
3
σ
2
a
i
.Therefore Λ satisfies the following inequality: (25)
Λ
≤
-
2
μ
+
α
a
q
+
i
-
2
a
β
s
s
-
m
+
β
s
+
ρ
β
v
-
γ
+
μ
+
c
γ
i
-
c
μ
r
+
3
a
σ
2
i
.Now we note that 2
a
β
s
(
s
-
m
)
≤
2
a
β, and so we can write (26)
Λ
≤
-
2
μ
+
α
a
q
+
i
β
-
γ
+
μ
+
γ
c
+
2
a
β
+
3
a
σ
2
-
r
c
μ
.The coefficients of q, i, and r are negative (see (13)). Furthermore, q, i, and r cannot all be zero since (27)
a
q
+
i
+
c
r
=
1
.Therefore Λ
<
0, and the proof is complete.We run two sets of simulations of the deterministic and stochasticS and I-trajectories over T
=
300 years. In the simulations, we use small perturbation parameters and small values of V. 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μ
=
0.016
,
ρ
=
0.1
,
α
=
0.3 and σ
=
0.65while the initial conditions areS
0
=
0.7, V
0
=
0.05 and I
0
=
0.15.In the first run (see Figures 1 and 2) we use β
=
0.55 and γ
=
0.6, and we obtain R
0
=
0.893
<
1. In the second run (see Figures 3 and 4) we choose β
=
0.2 and γ
=
0.1, which yields R
0
=
1.724
>
1.Figure 1
A simulation ofS for the deterministic and stochastic models with R
0
=
0.893.Figure 2
A simulation ofI for the deterministic and stochastic models with R
0
=
0.893.Figure 3
A simulation ofS for the deterministic and stochastic models with R
0
=
1.724.Figure 4
A simulation ofI for the deterministic and stochastic models with R
0
=
1.724.ForR
0
<
1, Theorem 4 guarantees the disease-free equilibrium to be almost surely exponentially stable and indeed the I-curves in Figure 2 appear to converge to 0. For R
0
>
1 the theorem did not predict almost sure exponential stability and in fact the I-curves in Figure 4 (and many simulations with R
0
=
1.724 not shown) certainly do not show any clear intention of converging to 0.


## 3. The Deterministic Optimal Control Problem

We now formulate and solve the deterministic version of the control problem. Recall from Section2 that in our SVIR model, α
(
t
) represents the fraction of the susceptible class being vaccinated at time t. We wish to design an optimal vaccination schedule, α
∗
(
t
), 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 [
0
,
T
]. The cost of the vaccination is assumed to be proportional to the square of α
(
t
).For the purposes of optimization we introduce the functionsf
1
(
t
), f
2
(
t
), and f
3
(
t
) appearing in the SDE system (2) as follows:(28)
f
1
t
=
μ
-
μ
S
t
-
β
I
t
S
t
-
α
t
S
t
,
f
2
t
=
α
t
S
t
-
ρ
β
V
t
I
t
-
μ
V
t
,
f
3
t
=
β
S
t
I
t
+
ρ
β
V
t
I
t
-
γ
I
t
-
μ
I
t
.Now we can formulate the optimization problem.Problem 5.
Minimize the objective function(29)
D
α
·
=
∫
0
T
α
2
t
+
c
I
t
d
tsubject to S
˙
(
t
)
=
f
1
(
t
), V
˙
(
t
)
=
f
2
(
t
), I
˙
(
t
)
=
f
3
(
t
),
S
(
0
)
=
S
0
≥
0, V
(
0
)
=
V
0
≥
0, I
(
0
)
=
I
0
≥
0.The control variable α
(
t
) is assumed to be a measurable function of time and bounded: 0
≤
α
(
t
)
≤
α
¯
≤
1.We solve Problem5 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 λ
1
(
t
), λ
2
(
t
), and λ
3
(
t
), 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 (30)
H
t
,
S
,
V
,
I
,
α
t
=
α
2
+
c
I
+
λ
1
f
1
t
,
S
,
V
,
I
,
α
+
λ
2
f
2
t
,
S
,
V
,
I
,
α
+
λ
3
f
3
t
,
S
,
V
,
I
,
α
.Theorem 6.
An optimal solution for Problem5 exists and satisfies the following system of differentiable equations:(31)
λ
1
˙
=
λ
1
μ
+
β
I
+
α
-
λ
2
α
-
λ
3
β
I
,
λ
2
˙
=
λ
2
ρ
β
I
+
μ
-
λ
3
ρ
β
I
,
λ
3
˙
=
-
c
+
λ
1
β
S
+
λ
2
ρ
β
V
-
λ
3
β
S
+
ρ
β
V
-
γ
-
μ
,with transversality conditions λ
1
(
T
)
=
0, λ
2
(
T
)
=
0, and λ
3
(
T
)
=
0. Furthermore, the optimal vaccination rate is given by(32)
α
∗
t
=
min
⁡
max
⁡
0
,
1
2
S
t
λ
1
t
-
λ
2
t
,
α
¯
.Proof.
In particular the Hamiltonian is convex with respect toα
(
t
) 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 λ
i
˙
(
t
) of the costate variables. The calculation of λ
1
˙
(
t
)
=
-
∂
H
/
∂
S, λ
2
˙
(
t
)
=
-
∂
H
/
∂
V, and λ
3
˙
(
t
)
=
-
∂
H
/
∂
Ileads 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, α
∗
(
t
). Since the function α
∗
(
t
) must optimize the Hamiltonian, we calculate(33)
∂
H
∂
α
=
2
α
-
S
λ
1
-
λ
2
.Consider now a fixed value of t. If 2
α
-
S
(
λ
1
-
λ
2
) is zero for some value α
(
t
) in the interval [
0
,
α
¯
], then the given value of α
(
t
) is optimal. If for every α
¯
¯
∈
[
0
,
α
¯
] we have 2
α
¯
¯
-
S
(
λ
1
-
λ
2
)
≥
0 (resp., 2
α
¯
¯
-
S
(
λ
1
-
λ
2
)
≤
0), then we must choose α
(
t
)
=
0 (resp., α
(
t
)
=
α
¯). Thus we must have (34)
α
∗
t
=
min
⁡
max
⁡
0
,
1
2
S
t
λ
1
t
-
λ
2
t
,
α
¯
.In Figure5 we plot for T
=
3 years the S, V, and I-trajectories of the deterministic model subject to the optimal control α
∗. We simulate α
∗ and α
∗
/
S
∗ for the deterministic model in Figure 6. In the simulations, we consider the parameters c
=
0.3, μ
=
0.016, β
=
0.55, ρ
=
0.1, γ
=
0.45, and α
¯
=
0.8 and initial conditions S
0
=
0.7, V
0
=
0.05, and I
0
=
0.2.Figure 5
A simulation ofS, V, and I for the deterministic model subject to the optimal control α
∗.Figure 6
A simulation ofα
∗ and α
∗
/
S
∗ for the deterministic model.In Figure6 the dashed curve represents α
∗
(
t
), 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 α
∗
(
t
) that minimizes the objective functional which for an initial state x
0 is defined by (35)
E
0
,
x
0
∫
0
T
α
2
s
+
c
I
s
d
s
.Here the expectation is obtained on the condition that the initial state (at time t
=
0) of the system is x
0. In step with the deterministic problem of earlier, we assume that there is a fixed constant α
¯
≤
1 with α
(
t
)
≤
α
¯ (a.s.). The class of admissible control laws is(36)
A
=
α
·
:
α
is adapted
,
and
0
≤
α
≤
α
¯
a.s.
.To solve this stochastic control problem, we define the performance criterion as follows:(37)
J
t
,
x
;
α
=
E
t
,
x
∫
t
T
α
2
s
+
c
I
s
d
s
,where the expectation is conditional on the state of the system being a fixed value x at time t. We define the value function as (38)
U
t
,
x
=
inf
α
·
∈
A
⁡
J
t
,
x
;
α
=
J
t
,
x
;
α
∗
.We determine a control law that minimizes the expected value J
:
A
→
R
+ 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 A as in (36) with J as in (37), find the value function(39)
U
t
,
x
=
inf
α
∈
A
⁡
J
t
,
x
;
α
,and an optimal control function (40)
α
∗
t
=
arg
⁡
inf
α
∈
A
⁡
J
x
;
α
t
∈
A
.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 (41)
α
∗
t
=
min
⁡
max
⁡
0
,
1
2
S
t
U
S
t
-
U
V
t
,
α
¯
.Proof.
We determine (41) via the dynamic programming approach. First we calculate L
U
(
t
): (42)
L
U
t
=
f
1
t
U
S
t
+
f
2
t
U
V
t
+
f
3
t
U
I
t
+
1
2
σ
S
t
I
t
2
U
S
S
t
+
1
2
ρ
σ
S
t
I
t
2
U
V
V
t
+
1
2
σ
S
t
+
ρ
V
t
I
t
2
U
I
I
t
+
ρ
σ
I
t
2
S
t
V
t
U
S
V
t
-
σ
I
t
2
S
t
S
t
+
ρ
V
t
U
S
I
t
-
ρ
σ
I
t
2
V
t
S
t
+
ρ
V
t
U
V
I
t
.Applying the Hamilton-Jacobi-Bellman theory (see, for instance, [24]) we must find the infimum: (43)
inf
α
∈
A
⁡
α
2
t
+
c
I
t
+
L
U
t
.For this purpose, we need to find the partial derivative of the expression (44)
α
2
t
+
c
I
t
+
L
U
twith respect to α, and this derivative should vanish. This leads to the equation:(45)
2
α
t
-
S
t
U
S
t
+
S
t
U
V
t
=
0
.We consider the bounds on α and by an argument similar to that in the proof of the deterministic case; the asserted expression for α
∗
(
t
) emerges.


## 5. Numerical Example

In the discussion that follows, the computations are done forT
=
3 years and we consider three different values of the perturbation parameter, σ. In particular we considertheσ-values: σ
∈
{
0.15,0.30,0.45
}together with the following parameter values and initial conditions:c
=
0.3, μ
=
0.016, β
=
0.55, ρ
=
0.1, γ
=
0.45, and α
¯
=
0.8,
S
0
=
0.7, V
0
=
0.05, and I
0
=
0.2.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 λ
1
-
λ
2 as a proxy for U
S
-
U
V in the calculation of α
∗ in this case. We note that the presence of S
(
t
) makes U into a stochastic variable even with the said proxy (in the stochastic case).For each value ofσ, we show in Table 1 the U
(
0
)-values obtained as the average over 3000 runs made for different candidates for α
^ as we take α
=
S
α
^. In these runs we take λ
1 and λ
2 of the underlying deterministic model.Table 1
U
(
0
) for different candidates of α
^ and σ
α
^
i
σ
=
0.15
σ
=
0.3
σ
=
0.45
α
^
-
3
29.84672
34.22981
36.21729
α
^
-
2
29.26219
34.09944
36.15703
α
^
-
1
28.96307
34.07355
36.14670
α
^
0
28.91028
33.93111
35.95781
α
^
1
29.52337
34.34593
36.26951
α
^
2
30.15975
34.56046
36.36480
α
^
3
31.01890
34.80816
36.47420
α
^
4
29.03129
34.07210
36.05649
α
^
5
30.71396
35.24172
36.82271Using a contact rateβ
0
=
β
-
(
σ
2
/
2
) instead of β, we choose the first candidate (46)
α
^
0
t
=
1
2
λ
1
t
-
λ
2
t
.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(47)
α
^
i
t
=
α
^
0
t
+
ϵ
i
-
ϵ
i
3
t
,
for
i
=
-
3
,
-
2
,
-
1,1
,
2,3corresponding to the rate β and with ϵ
=
0.001.Furthermore, we consider(48)
α
^
4
t
=
1
2
λ
1
t
-
λ
2
t
,corresponding to the contact rate β.Just for further comparison we choose a linear curve candidateα
^
5
(
t
), say, (49)
α
^
5
t
=
η
+
φ
twith η
=
0.0001 and φ
=
0.00015, also corresponding to a contact rate of β.Let us use the notationU
i
(
0
) for the value of U
(
0
) corresponding to the control S
α
^
i. We notice that U
4
(
0
) and U
5
(
0
) are bigger than U
0
(
0
). Even more interesting is the pattern U
-
3
(
0
)
>
U
-
2
(
0
)
>
U
-
1
(
0
)
>
U
0
(
0
) and U
0
(
0
)
<
U
1
(
0
)
<
U
2
(
0
)
<
U
3
(
0
). This creates the impression that U
0
(
0
) is close to being a minimum. At least, with all the difficulties of a more precise solution, the choice of α
^
0 seems like a viable option in a real application. The inequalities  U
-
3
(
0
)
>
U
-
2
(
0
)
>
U
-
1
(
0
)
>
U
0
(
0
) and U
0
(
0
)
<
U
1
(
0
)
<
U
2
(
0
)
<
U
3
(
0
) has been tested for other parameter choices too, with σ small. The simulations tested all revealed the same behaviour.In Figures7–9 we simulate the optimal solutions S
∗, V
∗, and I
∗ of the deterministic and stochastic (σ
=
0.3) models. We use the same values for the parameters c, μ, β, ρ, γ, σ, α
¯, and initial conditions S
0, V
0, and I
0 as in the discussion preceding these simulations.Figure 7
A simulation ofS
∗
(
t
) for the deterministic and stochastic models.Figure 8
A simulation ofV
∗
(
t
) for the deterministic and stochastic models.Figure 9
A simulation ofI
∗
(
t
) for the deterministic and stochastic models.An important point to note about our approximation is that it fully accommodates the stochasticity (embodied and concentrated in the factorS of the expression for α
∗
(
t
)). Therefore it gives a very good approximation, at least in the sense that in Table 1, the minimum value of U
(
0
) consistently corresponds to α
^
0.


## 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.



---

*Source: 271654-2015-05-24.xml*