Stochastic Reaction Network Modeling

and Optimal Control for Covid-19

Sidy Ly, Lena Tendeng, Mouhamadou A. M. T. Balde, Diene Ngom,

and Diaraf Seck

Abstract In this paper, we aim to build and study a stochastic SEIRS model with
vaccination in which we consider two types of infected: undetected and confirmed.
The latter are hospitalized and receive treatment. We also derive the Hamilton Jacobi
Bellman (HJB) equation by using dynamic programming tools, and we study it from
theoretical and numerical points of view.

Keywords Stochastic differential equations · Stochastic optimization · Dynamic

programming principle · Hamilton-Jacobi-Bellman equations · Numerical
simulations · Covid-19

1 Introduction

The Covid-19 pandemic was first reported in Wuhan in December 2019, and more
than two years after, the disease has still swept across. Several measures have been
adopted around the world to limit the evolution of the pandemic. In this work, we
seek to reduce disease transmission through optimal control. Data used for models
are generally incomplete or tainted with errors; hence a stochastic model helps
to take into account random phenomena. We propose a stochastic model obtained
by the reactions network approach. We prove the existence and uniqueness of the
solution of the stochastic differential equation deriving from the model. Control
is applied to the detection rate v2 and the cure rate γ of the model, and we use the
dynamic programming principle tools to derive the Hamilton Jacobi-Bellman (HJB)

S. Ly () · L. Tendeng · M. A. M. T. Balde · D. Seck

Laboratory of Mathematics of Decision and Numerical Analysis, University of Cheikh Anta
Diop., Dakar, Senegal
e-mail: sidy1.ly@ucad.edu.sn; mouhamadouamt.balde@ucad.edu.sn; diaraf.seck@ucad.edu.sn
D. Ngom
Dpartement de Mathmatiques UFR des Sciences & Technologies, Universit Assane Seck de
Ziguinchor, Ziguinchor, Senegal
e-mail: dngom@univ-zig.sn

© The Author(s), under exclusive license to Springer Nature Switzerland AG 2022 121
D. Seck et al. (eds.), Nonlinear Analysis, Geometry and Applications, Trends
in Mathematics, https://doi.org/10.1007/978-3-031-04616-2_5
122 S. Ly et al.

equation. We prove the existence and the uniqueness of solution of such equation as
in [3–7] and [8]. HJB equation in high dimension is numerically difficult to resolve,
so we use in this paper the Crank-Nicholson scheme.
The paper is organized as follows: in Sect. 2, we build the stochastic model
and prove the existence and uniqueness of the SDE solution. The optimal control
problem is studied in Sect. 3 and numerical simulations given at least in Sect. 4.

2 Modeling

In this section, we want to model the situation related to the COVID-19 pandemic
by the stochastic reaction network approach as in [1]. For this, we design by
X(t) = [S(t); E(t); U (t); V (t); H (t); R(t)]T the size of population in time t. In
this model, we study the variation of population 'X during an interval of time
't. This interval 't is assumed sufficiently small that it is impossible to have the
passage of several individuals from one compartment to .another at the same time.
Here, we designate by compartment every component of X(t), which means that we
have six of them. We also assume that natural mortality is negligible compared to
that caused by disease. In this model, newborns are not susceptible to contract the
disease, and those who are recovered are not immunized.
The coefficients used to switch from a comportiment to another as well as
the meaning of each comportiment are described in Table 1. The graph (Fig. 1)
summarizes in a concise way the model.
Under these assumptions, there are 12 possibilities for a population change 'Xi
for i = 1, . . . , 12 if we neglect multiple births, deaths or transformations in time 't
which have probabilities of order ('t)2 . These possibilities are listed in the Table 2
along with their corresponding probabilities. 'X1 = [−1; 1; 0; 0; 0; 0]T means that
one susceptible individual become exposed to the disease during time interval 't
and the probability of this event is given by p1 = (βu SU + βv SV )'t.

⎡ ⎤
δR − βu SU − βv SV − ηS
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ βu SU + βv SV − v1 E − v2 E ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
# ⎢ v E − d U − ξ U ⎥
12
⎢ 1 1 ⎥
E('X) = ⎢
pi 'X = ⎢
i ⎥ 't
⎥
⎢ ⎥
i=1 ⎢ v2 E − αV − d2 V ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ αV − d3 H − γ H ⎥
⎢ ⎥
⎣ ⎦
ξ U + γ H + ηS − δR
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 123

Table 1 Variables and parameters used in the model

Variable Description
βu Transmission rate of undetected
βv Transmission rate of confirmed
v1 Rate of incubation or latence of undetected
v2 Rate of incubation or latence of confirmed
d1 Death rate of undetected
d2 Death rate of confirmed
d3 Death rate of hospitalized patients
α Rate of hospitalized from confirmed
γ Cure rate of hospitalized
δ Re-susceptibility rate
η Vaccination rate
ξ Rate of recovered from undetected
S(t) Number of individuals susceptible to contracting the disease at time t
E(t) Number of individuals exposed to the disease at time t
U (t) Number of unknown infected at time t
V (t) Number of known infected at time t
H (t) Number of infected individuals hospitalized at time t
R(t) Number of individuals recovered from the disease at time t

Fig. 1 Model 4
124 S. Ly et al.

Table 2 Possible changes in Change Probability

the population with
corresponding probabilities 'X1 = [−1; 1; 0; 0; 0; 0]T p1 = (βu SU + βv SV )'t
'X2 = [0; 0; −1; 0; 0; 1]T p2 = ξ U 't
'X3 = [0; −1; 1; 0; 0; 0]T p3 = v1 E't
'X4 = [0; −1; 0; 1; 0; 0]T p4 = v2 E't
'X5 = [0; 0; −1; 0; 0; 0]T p5 = d1 U 't
'X6 = [0; 0; 0; −1; 0; 0]T p6 = d2 V 't
'X7 = [0; 0; 0; −1; 1; 0]T p7 = αV 't
'X8 = [0; 0; 0; 0; −1; 0]T p8 = d3 H 't
'X9 = [0; 0; 0; 0; −1; 1]T p9 = γ H 't
'X10 = [1; 0; 0; 0; 0; −1]T p10 = δR't
'X11 = [−1; 0; 0; 0; 0; 1]T p11 = ηS't
#
12
'X12 = [0; 0; 0; 0; 0]T p12 = 1 − pi
i=1

and
E('X('X)T ) =
⎡ ⎤
a11 −βu SU − βv SV 0 0 0 −δR − ηS
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ −βu SU − βv SV a22 −v1 E −v2 E 0 0 ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢0 −v1 E a33 0 0 −ξ U ⎥
⎢ ⎥
⎢ ⎥ 't.
⎢ ⎥
⎢0 −v2 E a44 −αV ⎥
⎢ 0 0 ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢0 0 0 −αV a55 −γ H ⎥
⎢ ⎥
⎣ ⎦
−δR − ηS 0 −ξ U 0 −γ H a66

where a11 = δR + βu SU + βv SV + ηS, a22 = βu SU + βv SV + v1 E + v2 E,

a33 = v1 E + d1 U + ξ U , a44 = v2 E + αV + d2 V , a55 = αV + d3 H + γ H , and
a66 = γ H + δR + ηS + ξ U .
We now define the expectation vector f and the 6 × 6 symmetric positive definite
covariance matrix G.

f(t, X) = E('X)/'t and G(t, X) = E('X('X)T )/'t.

Stochastic Reaction Network Modeling and Optimal Control for COVID-19 125

Refering to [1] that leads us to the stochastic differential equation system:

⎧
⎨ dX(t) = f(t, X)dt + G1/2 (t, X)dW(t)
(1)
⎩
X(0) = X0

with W(t) = [W1 (t); W2 (t); W3 (t); W4 (t); W5 (t); W6 (t)]T is the six-dimensional
Brownian motions.

3 Optimal Control Problem

3.1 Study of Stochastic Differential Equation

The aim of this section is to solve stochastic optimization problem by putting a

control on the detection and cure rate of the disease. It is thus about to replace
respectively v2 et γ by a1 (t)v2 et a2 (t)γ where a1 (t) is a parameter that measures
the effectiveness of detection campaigns and a2 (t) a parameter that measures the
effectiveness of treatment for patients 0 ≤ ai (t) ≤ 1 for i = 1; 2. It means that
fonctions f and g = G1/2 depend on the control a = (a1 ; a2 ). So system (1) can be
rewritten as the following form:
⎧
⎨ dX(t) = f(t, X, a)dt + G1/2 (t, X, a)dW(t)
(2)
⎩
X(0) = X0

where W = (W )t ≥0 is s the six-dimensional brownian motion starting at 0 and

adapted to a fixed filtration (F )t ≥0,
f : [0; T ] × [0; Pmax ]6 × [0; 1]2 −→ R6
G : [0; T ] × [0; Pmax ]6 × [0; 1]2 −→ R6×6
a = (a1 ; a2 ) is progressively measurable with values in [0; 1]2. We name by A the
set of all control processes and Pmax the total size of population.
Before solving our stochastic optimization problem, let us show that system (2)
has a unique solution.
Theorem 3.1 Assume that X0 is independent of the future of the Brownian motion
beyond time t = 0 and for any t ∈ [0, T [, a ∈ A and for X ∈ [0; Pmax ]6 a
progressively measurable process such that for any T > 0
 T
E( |X0 |2 ds) < ∞.
t

Then the system (2) has a unique solution.

126 S. Ly et al.

Proof For the proof we just need to show the K-lipschitz continuity and the linear
growth properties of f and g.
To simplify let’s set ϕ = f, G.
K-Lipschitz Continuity
ϕ is K-lipschitz continuity since it is continuously derivable in a compact set. Then
there is K > 0 such that:

||ϕ(t, x, a) − ϕ(t, x̃, a)|| ≤ K||x − x̃||, ∀ t ∈ [0, T ], x, x̃ ∈ [0, Pmax ]6 .

Linear Growth Property

We can see that the functions f, G can be written in the form f = (fk )k=1,··· ,6 ,
#6 #
6
with fk = k
αi,j xi xj + αik xi , and G = (Gl,m )l,m=1,··· ,6 , with Gl,m =
i,j =1 i=1
#
6
l,m
#
6
αi,j xi xj + αil,m xi . Where X = (xi )i=1,··· ,6 ∈ [0, Pmax ]6 .
i,j =1 i=1

#
6 #
6
||f (t, x, a)|| = max |fk | = max | k
αi,j xi xj + αik xi |
k k
i,j =1 i=1

#
6 #
6
≤ max( |αi,j
k
||xi xj | + |αik ||xi |)
k
i,j =1 i=1

≤ K1 (1 + ||x||),

#
6 #
6
with K1 = max(Pmax |αi,j
k
|+ |αik |).
k
i,j =1 i=1

#
6 #
6 #
6
l,m
#
6
||G(t, x, a)|| = max |Gl,m | = max | αi,j xi xj + αil,m xi |
l l
m=1 m=1 i,j =1 i=1

#
6 #
6
l,m
#
6
≤ max ( |αi,j ||xi xj | + |αil,m ||xi |)
l
m=1 i,j =1 i=1

≤ K2 (1 + ||x||),

#
6
l,m
#
6
with K2 = max(Pmax |αi,j |+ |αil,m |).
l
m,i,j =1 m,i=1
Taking K = max {K1 ; K2 } we find the result. 
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 127

3.2 Position of the Problem

The aim is to act on the system through detection campaigns and treatment offered to
hospitalized patients. The effectiveness of detection campaign and that of treatment
at time t depend respectively on two parameters a1 (t) and a2 (t) (0 ≤ ai (t) ≤ 1 for
i = 1; 2) which are the components of the control a ∈ A.
The idea is to reduce the disease’s spread and provide an effective treatment. This
is why we propose the following objective function
 T 
J (t, X, a) = E (L(s, X(s), a)) ds , (3)
t

where L : [0; T ] × [0; Pmax ]6 × [0; 1]2 −→ R is a continuous map.

Then our stochastic optimization problem is:

minJ (t, X, a) (4)

a∈A

under the constraints

⎧
⎨ dX(t) = f(t, X, a)dt + G1/2 (t, X, a)dW(t)
(5)
⎩
X(0) = X0

where f : [0; T ] × [0; Pmax ]6 × [0; 1]2 −→ R6 and G : [0; T ] × [0; Pmax ]6 ×
[0; 1]2 −→ R6×6 are given by:
⎡ ⎤
δR − βu SU − βv SV − ηS
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ βu SU + βv SV − v1 E − a1 (t)v2 E ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ v1 E − d1 U − ξ U ⎥
⎢ ⎥
f(t, X, a) = ⎢
⎢
⎥
⎥
⎢ a (t)v E − αV − d V ⎥
⎢ 1 2 2 ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ αV − d3 H − a2 (t)γ H ⎥
⎢ ⎥
⎣ ⎦
ξ U + a2 (t)γ H + ηS − δR
128 S. Ly et al.

G(t, X, a) =
⎡ ⎤
a11 −βu SU − βv SV 0 0 0 −δR − ηS
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢ −βu SU − βv SV a22 −v1 E −a1 (t)v2 E 0 0 ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢0 −v1 E a23 0 0 −ξ U ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢0 −a1 (t)v2 E −αV ⎥
⎢ 0 a44 0 ⎥
⎢ ⎥
⎢ ⎥
⎢ ⎥
⎢0 0 0 −αV a55 −a2 (t)γ H ⎥
⎢ ⎥
⎣ ⎦
−δR − ηS 0 −ξ U 0 −a2 (t)γ H a66
a11 = δR + βu SU + βv SV , a22 = βu SU + βv SV + v1 E + a1 (t)v2 E, a33 =
v1 E + d1 U + ξ U , a44 = a1 (t)v2 E + αV + d2 V , a55 = αV + d3 H + a2 (t)γ H ,
a66 = a2 (t)γ H + δR + ηS + ξ U
and W(t) = [W1 (t); W2 (t); W3 (t); W4 (t); W5 (t); W6 (t)]T is the six-dimensional
brownian motion.

3.3 Dynamic Programming Principle

To solve our stochastic optimization problem, let Z(t, X), known as the value
function, be the excepted value of the objective function J from t to T when an
optimal policy is followed from t to T , given Xt = x:

Z(t, x) = min J (t, x, a), . (6)

a∈A
*" +
T
with J (t, x, a) = E t L(s, Xs , a) ds .
The functional (6) gives the remaining optimal cost by assuming that we arrive
at x in time t < T and satisfies the following final condition:

Z(T , x) = 0 (7)

Let’s define the C 1,2 (!) the space of v : ! −→ R such that ∂t v, ∂xi v, ∂x2i v are all
continuous in t, xi , i = 1, · · · , 6. if ! = [0; T ] × R6 × [0; 1]2 then we have the
following proposition
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 129

Proposition 3.2 If the value function Z ∈ C 1,2 (!) then it is a solution of the
following Hamilton Jacobi Bellman equation:
⎧
⎪
⎪ ∂Z
⎨− + H (t, x, DZ, D 2 Z) = 0
∂t
(8)
⎪
⎪
⎩ Z(x, T ) = 0

Before proving proposition 3.2, let us recall It formula.

Theorem 3.3 (It Formula) Suppose that X(.) has a stochastic differential

dX = F dt + GdW,

for f ∈ L1 (0, T ), G ∈ L2 (0, T ). Assume u : R × [0; T ] −→ R is continuous and

∂u ∂ 2 u
that ∂u
∂t , ∂x , ∂x 2 exist and are continuous.
Set

Y (t) = u(X(t), t)

Then Y has the stochastic differential


∂u ∂u 1 ∂ 2u 2 ∂u
dY = + F+ G dt + GdW (9)
∂t ∂x 2 ∂x 2 ∂x

Remark 3.4 The expression (9) means for all 0 ≤ s ≤ r ≤ T

 
r ∂u ∂u 1 ∂ 2u 2 r ∂u
Y (r) − Y (s) = + F+ G dt + GdW almost surely.
s ∂t ∂x 2 ∂x 2 s ∂x

For more details see [2]

Proof Considering h > 0 and an arbitrary control (a1 ; a2 ) ∈ [0; 1]2 Then, by the
principle of optimality,
 t +h 
Z(t, x) ≤ E L(s, Xs , a)ds + Z(t + h, Xt +h ) . (10)
t

Assuming that Z is enough smooth, we have by It formula between t and t + h

 t +h 
∂Z 1
Z(t + h, Xt +h ) = Z(t, x) + + f.DZ + T rac G.D 2 Z (s, Xs )ds
 t +h t ∂t 2
+ DZ.G1/2 (s, X)dW(s)
t
(11)
130 S. Ly et al.

Substituting (11) into (10) and canceling term in Z, we have:

 t +h    
∂Z 1
0≤ E L(s, Xs , a) + + f.DZ + T rac G.D 2 Z (s, Xs ) ds .
t ∂t 2
(12)

Dividing by h and letting h −→ 0, we obtain the following inequality:


∂Z 1
0 ≤ L(t, Xt , a) + + f.DZ + T rac G.D 2 Z (t, Xt ). (13)
∂t 2

This is for any control (a1 ; a2 ) ∈ [0; 1]2 , so

∂Z 1
− (t, Xt ) − L(t, Xt , a) − f(t, Xt ).DZ(t, Xt ) − T rac G(t, Xt ).D 2 Z(t, Xt ) ≤ 0.
∂t 2
(14)

This mean that:

∂Z
− (t, Xt )+
∂t
 
1
sup −f(t, Xt ).DZ(t, Xt ) − T rac G(t, Xt ).D Z(t, Xt ) − L(t, Xt , a) ≤ 0.
2
a∈A 2
(15)

Other hand, we assume that a ∗ = (a1∗ ; a2∗ ) is an optimal control then by verification
theorem, we have:
 t +h 
Z(t, x) = E L(s, X∗s , a ∗ )ds + Z(t + h, X∗t +h ) . (16)
t

where X∗ is the state associated to a ∗ . By using similar arguments and regularity

condition for Z as in (10), we obtain:

∂Z 1
− (t, Xt )−L(t, Xt , a ∗ )−f(t, Xt ).DZ(t, Xt )− T rac G(t, Xt ).D 2 Z(t, Xt ) = 0.
∂t 2
(17)
Equation (17) combined with (14) shows that Z must satisfy the following equation:

∂Z
− (t, Xt )+
∂t
 
1
sup −f(t, Xt ).DZ(t, Xt ) − T rac G(t, Xt ).D Z(t, Xt ) − L(t, Xt , a) = 0.
2
a∈A 2
(18)
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 131

Let us introduce the Hamiltonian H of our problem: for p ∈ Rd and M ∈ Rd×d

 
1
H (t, x, p, M) = sup −L(t, Xt , a) − p.f − T rac(M.G) (19)
a∈A 2

Then the Hamilton-Jacobi equation can be rewritten as a terminal value

problem (8). 

3.4 Existence and Uniqueness of the HJB Equation

We denote by A the set of admissible controls defined on [0, T ].

Let us recall the value function:
⎧
⎪ inf J (t, x, a) ∀(t, x) ∈ [0, T [×[0; Pmax ]6
⎨ a∈ A
Z(t, x) = (20)
⎪
⎩
0 ∀(t, x) ∈ {T } × [0; Pmax ]6

Remark 3.5
1. U = [0, 1]2 is a Polish space.
2. The map ϕ = f, G1/2 or, L :
• is uniformly continuous
• satisfies the K-Lipschitz continuity and the linear growth property.
Theorem 3.6 The Hamilton Jacobi Bellman equation (8) has a unique solution
V ∈ C 1,2 (!) with ! = [0; T ] × [0; Pmax ]6 × [0; 1]2.
Proof Under the properties given recalled in the above remark, the existence and
uniqueness are ensured by some references like Theorem 1 p. 301 in [5] and
Theorem 6.2 p. 169 in [3]. 

4 Numerical Simulation

4.1 Numerical Scheme

We apply finite difference scheme of Crank-Nicholson. The Crank-Nicholson

scheme is universally stable in one dimension.
Let’s set the time step h = dt = T /N, with T the maximal time and N − 1 the
number of intervals for the time. We set H = dXi = maxi Xi /M, ∀i = 1, · · · , 6,
with X = (Xi ) = (S, E, U, V , H, R) the states. With M −1 the number of intervals
for the space.
132 S. Ly et al.

Let’s set a * the optimal control process such that (19) is satisfied. Then

∂Z 1
+ L(t, Xt , a * ) + f(t, Xt , a * ).DZ(t, Xt ) + T rac G(t, Xt , a * ).D 2 Z(t, Xt ) = 0
∂t 2
# 
1 #
6 6
∂Z * ∂Z ∂2Z
+ L(t, Xt , a ) +
*
fi (t, Xt , a ). + *
Gik (t, Xt , a ). =0
∂t ∂Xi 2 ∂Xi Xk
i=1 i,k=1

The space of states is 6 dimensions. So the discretization is complex and the

computation needs high memory and it is slow in time. Then for simplify let’s set
#6
∂Z ∂Z
the variable X̄ = αi Xi . Then = αi . Replacing Z(t, Xt ) by Z̄(t, X̄t )
∂Xi ∂ X̄
i=1
and dropping the “ ¯ ” for simplify, we obtain:

# 
1 #
6 6
∂Z ∂Z ∂ 2Z
+ L(t, X̄t , a * ) + αi fi (t, X̄t , a * ). + αi αk Gik (t, X̄t , a * ). =0
∂t
i=1
∂ X̄ 2
i,k=1
∂ X̄2

Using the above discretization and the Crank-Nicholson scheme we obtain:

! !
Zjl +1 − Zjl #
6
jl Zjl+1 − Zjl−1 Zjl+1 l−1
+1 − Zj +1
+ Lj l (a ) +
*
αi fi (a * ) . +
h 4H 4H
i=1
⎛ ⎞ !
1 # Zjl+1 − 2Zjl + Zjl−1 Zjl+1
+1 − 2Zj +1 + Zj +1
6 l l−1
+ ⎝ αi αk Gik (a * )⎠ .
jl
+ =0
2 2H 2 2H 2
i,k=1

Setting
)6 jl *
i=1 αi fi (a )h aj l h
Aj l = = (21)
4H 4H
)6 jl
1 i,k=1 αi αk Gik (a * )h 1 bj l h
Bj l = = , (22)
2 2H 2 2 2H 2
we obtain:

Zjl +1 − Zjl + Lj l (a * )h + Aj l Zjl+1 − Zjl−1 + Zjl+1 l−1

+1 − Zj +1

+Bj l Zjl+1 − 2Zjl + Zjl−1 + Zjl+1 l−1

+1 − 2Zj +1 + Zj +1 = 0
l
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 133

(−Aj l + Bj l )Zjl−1 l+1

+1 + (1 − 2Bj l )Zj +1 + (Aj l + Bj l )Zj +1
l

+(−Aj l + Bj l )Zjl−1 + (−1 − 2Bj l )Zjl + (Aj l + Bj l )Zjl+1 + Lj l h = 0

(23)

Let’s set

O h,H Zjl +1 = (−Aj l + Bj l )Zjl−1

+1 + (1 − 2Bj l )Zj +1 + (Aj l + Bj l )Zj +1
l l+1

P h,H Zjl = (−Aj l + Bj l )Zjl−1 + (−1 − 2Bj l )Zjl + (Aj l + Bj l )Zjl+1

Then the scheme (23) can be written as follows:

O h,H Zjl +1 + P h,H Zjl + Lj l h = 0.

Remark 4.1 In the next subsection showing the numerical result, we use αi =
#
6
1, i = 1, · · · , 6. Then X̄ = Xi = S + E + U + V + H + R.
i=1

Theorem 4.2 Let Z ∈ C 1,2 (!) be the unique solution of the HJB (8) and let Zjl
be the solution of the scheme (23). There are constants independent from h, H, V
such that

max |O h,H V (x l , tj +1 ) + P h,H V (x l , tj ) + L(x l , tj )| ≤ C max{|Vxt |, |Vxxt | }('t)2 + o(('x)2 )

l=1,··· ,M
j =1,··· ,N
(24)

Proof For simplicity, we set:

O h,H V (x l , tj +1 ) = (−Aj l + Bj l )V (x l−1 , tj +1 ) + (1 − 2Bj l )V (x l , tj +1 )

+ (Aj l + Bj l )V (x l+1 , tj +1 )

and

P h,H V (x l , tj ) = (−Aj l + Bj l )V (x l−1 , tj ) + (−1 − 2Bj l )V (x l , tj ) + (Aj l + Bj l )V (x l+1 , tj ).

Using Taylor expansion in space on V (x l−1 , tj +1 ) and V (x l+1 , tj +1 ), we obtain:

O h,H V (x l , tj +1 ) = V (x l , tj +1 ) + 2Aj l 'xVx (x l , tj +1 )

+ 2Bj l ('x)2 Vxx (x l , tj +1 ) + o(('x)2 ). (25)
134 S. Ly et al.

Again applying a Taylor expansion in time yields:

O h,H V (x l , tj +1 ) = V (x l , tj ) + 'tVt (x l , tj ) + 2Aj l 'xVx (x l , tj ) + 2Aj l 'x'tVxt (x l , tj )+

2Bj l ('x)2 Vxx (x l , tj ) + 2Bj l ('x)2 'tVxxt (x l , tj +1 ) + o(('x)2 ).

(26)

Applying the Taylor expansion in space to P h,H V (x l , tj ), we obtain:

P h,H V (x l , tj ) = −V (x l , tj ) + 2Aj l 'xVx (x l , tj ) + 2Bj l ('x)2 Vxx (x l , tj ) + o(('x)2 ).

(27)

By summing, we get:

O h,H V (x l , tj +1 ) + P h,H V (x l , tj ) = 'tVt (x l , tj ) + 4Aj l 'xVx (x l , tj ) + 4Bj l ('x)2 Vxx (x l , tj )

+2Aj l 'x'tVxt (x l , tj ) + 2Bj l ('x)2 'tVxxt (x l , tj +1 ) + o(('x)2 )

(28)

Replacing Vt by using the (8), we obtain:

O h,H V (x l , tj +1 ) + P h,H V (x l , tj ) + L(x l , tj ) = 2Aj l 'x'tVxt (x l , tj )+

2Bj l ('x)2 'tVxxt (x l , tj +1 ) + o(('x)2) (29)

Using (21) and (22), we obtain:

|O h,H V (x l , tj+1 ) + P h,H V (x l , tj ) + L(x l , tj )| = 2ajl ('t)2 Vxt (x l , tj )+

bjl ('t)2 Vxxt (x l , tj+1 ) + o(('x)2 )

|O h,H V (x l , tj+1 ) + P h,H V (x l , tj ) + L(x l , tj )| ≤ max{|ajl ||Vxt |, |bjl ||Vxxt |}('t)2 + o(('x)2 ).
j,l
(30)

Finally, we can deduce (24). 

Remark 4.3 The Crank-Nicholson scheme (23) is consistent and is second order in
time and space, then by the Lax theorem it is convergent.

4.2 Covid19 Case Study: SDE

Some tests are performed using Covid19 data. In this subsection, we show numerical
results of the SDE (5).
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 135

Fig. 2 Stochastic vs deterministic differential equation: we show multiple SDE curves and
deterministic curves. Ydet stands for the deterministic curve of the states Y and Ystoc stands for
the stochastic curve of the state Y. (a) Test 1. (b) Test 2. (c) Test 3

The initial population is set to 8·106 and the initial states to: S0 = 785.7353·103,
E0 = 150.42 · 103 , U0 = 1.2043 · 103 , V0 = 7.4412 · 103 , H0 = 4.8786 · 103 ,
R0 = 1.9154 · 103 (Fig. 2).
The following table gives the parameters for the three tests we performed:
• Test 1: low vaccination rate η = 0.00127269, low re-susceptibility rate δ =
6.36347 10−6 and high hospitalization rate α = 0.7.
• Test 2: high vaccination rate η = 0.0127269, high re-susceptibility rate δ =
6.36347 10−4 and high hospitalization rate α = 0.07.
• Test 3: high vaccination rate η = 0.0127269, high re-susceptibility rate δ =
6.36347 10−2 and low hospitalization rate α = 0.07.
136 S. Ly et al.

Table 3 Table of parameters: values of parameters for different tests

Parameters Values of test 1 Values of test 2 Values of test 3
δ 6.36347 · 10−6 6.36347 · 10−4 6.36347 · 10−2
βu 2.65884 · 10−8 2.65884 · 10−8 2.65884 · 10−8
βv 2.53223 · 10−8 2.53223 · 10−8 2.53223 · 10−8
η 0.00127269 0.0127269 0.0127269
v1 0.141481 0.141481 0.141481
v2 0.141481 0.141481 0.141481
d1 0.001 0.001 0.001
d2 0.001 0.001 0.001
d3 0.001 0.001 0.001
ξ 0.05 0.05 0.05
α 0.7 0.07 0.07
γ 0.0667 0.0667 0.0667

4.3 Covid19 Case Study: The Optimal Control Problem

In this subsection, we show numerical results of the Eq. (8) associated with the
SDE (5). The values of parameters are given in the Table 3.
We consider two cases with the three tests. The first case consists on using the
functional L(t, X(t), u(t)) = E(t)/R(t) (R(t) = 0 for any t ∈ [0; T ]) and the
second is about the functional L(t, X(t), u(t)) = AE(t)+BH (t)+Cu21 (t)+Du22 (t)
where u = (u1 , u2 ) is the 2-dimensional control. In the numerical simulation we use
A = B = C = D = 1.
We use bang-bang method that is we search for optimal controls within the finite
set of possible controls.
The SDE is solved by using Euler-Maruyama approximation. Then we get X̄ =
#6
αi Xi .
i=1
Finally, we find Z(t, X̄) by solving HJB.

5 Discussion

Simulations performed in this paper are based on actual data from covid in Senegal.
To achieve them, we used two functionals in the optimal control problem. The first
one being L(t, X(t), u(t)) = E(t)/R(t) (R(t) = 0 for any t ∈ [0; T ]) represents
ratio between exposed individuals and recovered ones at time t.
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 137

Second expression is L(t, X(t), u(t)) = AE(t) + BH (t) + Cu21 (t) + Du22 (t) and
it represents the weighted sum between exposed individuals, hospitalized ones and
square of Euclidian norm of the control. These simulations were made with three
tests of each functional.
First test: Low vaccination rate, low re-susceptibility rate, and low hospitalization
rate.
We notice that with the first function, the control’s bang-bang effect is
very intense (Fig. 3c). That means that many efforts must be put into the
detection campaign of the disease and the effectiveness of the treatment
so that there should be fewer exposed individuals and more recovered
ones (Fig. 3a).The value function curve (Fig. 3e) shows that our goal
will be reached if the forces on the control are provided.
On the other hand, the controls bang-bang effect is less intense with
the second function. We remark that after 50 days, there is no detection
campaign and that the proposed treatment is maximum (Fig. 3d). This
situation is justified by the fact that after 50 days, there are no longer any
exposed individuals and that the number of recovered ones increases.
The curve of the value function reflects exactly this phenomenon
(Fig. 3f).
Second test: high vaccination rate, high re-susceptibility rate, and high hospital-
ization rate.
We observe the same situation as in the first test, i.e., controls bang-
bang effect is very intense (Fig. 4c). Thus translating that many efforts
are needed in the detection campaigns and in the treatments to reach
our goal (Fig. 4e). With the second function, we notice that, to achieve
our objective (Fig. 4f), i.e., to reduce the progression of the disease,
it is necessary that after 30 days to focus on the effectiveness of the
treatment than on the detection campaign because after 30 days it is
equal to zero (Fig. 4d).
Third test: high vaccination rate, high re-susceptibility rate, and low hospitaliza-
tion rate.
We remark that with the first functional, even with a few hospitalized
patients, many efforts are required in the detection and the treatment
(Fig. 5a) and (Fig. 5c). While with the second functional, we must
provide many efforts but not as intense as those of the first test to
achieve our goal (Fig. 5d).
138 S. Ly et al.

Fig. 3 Optimal states, controls and value function for test 1 with the first case on the left and
the second case on the right. (a) The states vs time in days. (b) The states vs time in days. (c)
The controls in respect to time in days. (d) The controls in respect to time in days. (e) The value
function at boundary vs time in days. (f) The value function at boundary vs time in days. (g) Value
function in 3D. (h) Value function in 3D
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 139

Fig. 4 Optimal states, controls and value function for test 2 with the first case on the left and
the second case on the right. (a) The states vs time in days. (b) The states vs time in days. (c)
The controls in respect to time in days. (d) The controls in respect to time in days. (e) The value
function at boundary vs time in days. (f) The value function at boundary vs time in days. (g) Value
function in 3D. (h) Value function in 3D
140 S. Ly et al.

Fig. 5 Optimal states, controls and value function for test 3 with the first case on the left and
the second case on the right. (a) The states vs time in days. (b) The states vs time in days. (c)
The controls in respect to time in days. (d) The controls in respect to time in days. (e) The value
function at boundary vs time in days. (f) The value function at boundary vs time in days. (g) Value
function in 3D. (h) Value function in 3D
Stochastic Reaction Network Modeling and Optimal Control for COVID-19 141

6 Conclusion

In this work, we built a stochastic model for COVID-19 on which we applied the
optimal control theory intending to reduce disease transmission. We prove the exis-
tence and uniqueness of solutions of both SDE governing the model and Hamilton
Jacobi equation. Numerical simulations of the SDE show that simultaneous and
adequate choice of vaccination and detection rates decreases disease transmission.
The numerical resolution of the Hamilton Jacobi equation shows that the effort
required for treating patients is far greater than that needed for detecting cases of
COVID-19. As perspectives to this work, we will seek to simulate the Hamilton
Jacobi equation by another method where a variable change will not be performed
to reduce the dimension. Likewise, we aim to investigate viscosity solutions and
machine learning tools to solve the Hamilton Jacobi PDE.

Acknowledgments This work was completed with the support of the NLAGA project.

References

1. E. Allen, Modeling with Ito Stochastic Differential Equations (Springer, Dordrecht, 2007)
2. L.C. Evans, Partial Differential Equations. Graduate Studies in Mathematics, vol. 19 (Ameri-
can Mathematical Society, Providence, 2000)
3. W.H. Fleming, R.W. Rishel, Deterministic and Stochastic Optimal Control (Springer, Berlin,
1975)
4. N.V. Krylov, Controlled Diffusion Processes (Springer, New York, 1980)
5. N.V. Krylov, Nonlinear Elliptic and Parabolic Equations of the Second Order (Reidel,
Dordrecht, 1987)
6. O.A. Ladyzhenskaya, V.A. Solonikov, N.N. Uralceva, Linear and Quasi linear Equations of
Parabolic Type (AMS, Providence, 1968)
7. S. Peng, Generalized dynamic programming principle and Hamilton-Jacobi-Bellman equation.
Stochastics Stochastics Rep. 38, 119–134 (1992)
8. J. Yong, X.Y. Zhou, Stochastic Controls Hamiltonian Systems and HJB Equations (Springer,
Berlin, 1999)