Optimal control and analysis of a stochastic SEIR epidemic model with nonlinear incidence and treatment

1.   Introduction

Mathematical models are very effective tools in studying the dynamical behavior of infectious diseases ^[1,2,3,4,5]. When there is no way to eradicate diseases completely, researchers are always looking for and developing the best methods to control the spread of diseases. A lot of mathematical models have been presented for control effects of infectious diseases^[6,7,8]. Disease control is mainly considered from two aspects: vaccine and treatment. Some researchers have considered vaccine control, such as ^[9,10,11], while others considered treatment control. Among them, in order to measure the effect of delayed treatment of the infected, Zhang and Liu ^[12] proposed the form of saturated treatment function $T(I), \; T(I) = \frac{\gamma I}{1+\alpha I}$, where $\gamma > 0$, $\alpha > 0$. A treatment function containing both the control and the infected is proposed and is defined as $T(u, I) = \frac{\phi uI}{1+\alpha uI}$ in ^[13], which better reflects the characteristics of natural epidemics.

Many human epidemics, such as measles, smallpox, epidemics, and dengue fever, are represented by the SEIR model ^[14,15]. In particular, the literature ^[15] takes into account the Crowley-Martin-type incidence rate and Holling type Ⅱ treatment rate, and proposes the following model:

where the total population $N$ is divided into four parts: the susceptible$(S)$, exposed$(E)$, infected$(I)$, and recovered$(R)$, $\frac{\beta_{1}SI}{(1+\alpha_{1}S)(1+\gamma I)}$ represents the transmission population of disease from $S$ to $I$ by the Crowley-Martin incidence rate, $\frac{\beta_{2}I}{1+\alpha_{2}I}$ represents the treatment rate of the infected population, and $\frac{auI}{1+buI}$ is the saturated treatment function of infected population where $u$ is treatment control. Other parameters and their definition are shown in Table 1. $b$ is a nonnegative quantity, and other parameter are all positive. Neglecting the fourth equation, they considered an equivalent model where the basic reproduction number was described as

They have performed the stability and bifurcation analysis of the model system. If $R_{0} < 1,$ then symtem $(1.1)$ has a unique disease-free equilibrium $P^{0}(\frac{\Lambda}{\mu}, 0, 0)$, which is locally asymptotically stable. Conversely, if $R_{0} > 1$, then system $(1.1)$ has two equilibrium points: one disease-free equilibrium $P^{0}$ that is unstable, and another endemic equilibrium $P^{*}(S^{*}, E^{*}, I^{*})$ that is locally asymptotically stable, where $(S^{*}, I^{*})$ is presented numerically, $E^{*} = \frac{\beta_{1}S^{*}I^{*}}{(\mu+\varphi)(1+\alpha_{1}S^{*})(1+\gamma I^{*})}$.

In real life, the spread of diseases is inevitably affected by environmental white noise, as it is an integral part of nature, therefore, considering deterministic models no longer fits the actual needs. Some scholars have studied the dynamical behaviors of epidemic models affected by white noise, such as ^{[16,17,18,19]}. Hence, we incorporate white noise perturbations into model $(1.1)$. We propose the following stochastic SEIR epidemic model with nonlinear incidence and treatment.

where independent standard Brownian motions are expressed as $B_{i}(t) \; (i$$=$$1, 2, 3, 4)$, and $\sigma_{i} \; (i = 1, 2, 3, 4)$ are positive constants that represent the intensity of the environment white noise, respectively.

This paper is organized as follows: In Section 2, we prove that the proposed model has a unique positive solution and the solution is global. In Sections 3 and 4, we study the dynamical behaviors of the proposed model in terms of the existence of stationary distribution and the extinction of the disease, respectively. In Section 5, we discuss the optimal control problem of the proposed model. In Section 6, we give a series of numerical simulations. Finally, in Section 7, conclusions are given.

Let $(\Omega, \mathcal{F}, \{\mathcal{F}_{t}\}_{t\geq 0}, \mathbb{P})$ be a complete probability space with a filtration $\{\mathcal{F}_{t}\}_{t\geq 0}$ that satisfies the usual conditions (i.e., it is increasing and right continuous while $\mathcal{F}_{0}$ contains all $\mathbb{P}$-null sets).

Consider the stochastic differential equation (SDE) of $n$-dimensional of the form

where $F(t, X):\mathbb{R_{+}}\times\mathbb{R}^{n}\rightarrow \mathbb{R}^{n}$ and $G(t, X):\mathbb{R_{+}}\times\mathbb{R}^{n}\rightarrow \mathbb{R}^{n\times m}$ are measurable functions and $B(t)$ is $\mathbb{R}^{m}$-valued standard Brownian motion. Given $V(X, t)\in C^{2, 1}(\mathbb{R}^{n}\times\mathbb{R}_{+}, \mathbb{R}_{+})$, we define the operator $\mathcal{L}V$ corresponding to the SDE $(1.3)$ by

where

Then, the Itô formula is obtained:

2.   Existence and uniqueness of the global positive solution

In this section, using the Lyapunov analysis method^[20], we first show that the system (1.2) has a unique local positive solution, then we show that this solution is global.

Theorem 1. If $(S(0), E(0), I(0), R(0))\in \mathbb{R}_{+}^{4}$ is any initial value to $(1.2)$, then $(S(t), E(t), \; I(t), R(t))$ is a unique existing positive solution to $(1.2)$ for $t\geq 0$ and the solution remains in $\mathbb{R}_{+}^{4}$ with probability 1.

Proof. Since the local Lipschitz condition is satisfied by system $(1.2)$, for any initial value $(S(0), E(0), I(0), R(0))\in \mathbb{R}_{+}^{4}$, there exists a unique local solution $(S(t), E(t), I(t), R(t))$ for $t\in[0, \tau_{e})$, where $\tau_{e}$ denotes the explosion time ^[21]. To prove that the solution is global, we only need to prove $\tau_{e} = + \infty$ a.s. To this end, let $k_{0}\geq1$ be a sufficiently large constant such that $S(0), E(0), I(0)$, and $R(0)$ lie within the interval $[\frac{1}{k_{0}}, k_{0}]$. For $k\geq k_{0}$, we define the stopping time as follows:

Clearly, $\tau_{\infty}\leq\tau_{k}$ a.s. If $\tau_{\infty} = +\infty$ a.s., then we have $\tau_{e} = +\infty$ a.s., and $(S(t), E(t), I(t), R(t))\in \mathbb{R}_{+}^{4}$ a.s. If this is false, then there exists a pair of constants $T > 0$ and $\varepsilon\in(0, 1)$ such that $\mathbb{P}\{\tau_{\infty \leq T}\} > \varepsilon$. Therefore, there is an integer $k_{1}\geq k_{0}$ satisfying

Define the $C^{2}$-function $V_{1}:\mathbb{R}_{+}^{4}\rightarrow \mathbb{R}_{+}^{4}$:

Applying the Itô formula, we obtain

where $K$ is a positive constant.

The following proofs are similar to references^[22]. □

3.   Stationary distribution

The unique stationary distribution of the stochastic SEIR model indicates that the persistence of the disease in the future under certain conditions, that is, the stochastic model fluctuates around the endemic equilibrium of the corresponding deterministic model.

Let $X(t)$ be a regular time-homogeneous Markov process described by the following stochastic differential equation in $\mathbb{R}^{d}$:

The diffusion matrix of the process $X(t)$ is defined as follows:

Lemma 1. ^[23] If there is a bounded open domain $D\subset E_{d}$ with regular boundary $\Gamma$, it has the following properties:

(i) The diffusion matrix $A(x)$ is strictly positive definite for all $x\in D$;

(ii) For any $x\in E_{d}\backslash D$, it has a nonnegative $C^{2}$ function $V$ such that $\mathcal{L}V$ is negative.

Then there exists a unique ergodic stationary distribution $\pi(\cdot)$ for the Markov process $X(t)$.

Theorem 2. For any initial value $(S(0), E(0), I(0), R(0))\in\mathbb{R}_{+}^{4}$, the system $(1.2)$ admits a unique ergodic stationary distribution $\pi(\cdot)$, if

Proof. In order to prove the theorem, we first verify that condition (i) in Lemma 1 holds. From $(1.2)$, we obtain that the diffusion matrix of system $(1.2)$ is

It is easy to see that the matrix $A$ is positive definite for any compact subset of $\mathbb{R}_{+}^{4}$. Therefore, condition (ⅰ) of Lemma 1 is satisfied.

Next, we verify that the condition (ⅱ) in Theorem 2 also holds. Define $C^{2}$ functions $V_{1}:\mathbb{R}_{+}^{4}\rightarrow \mathbb{R}$:

where

Making use of the Itô formula, we obtain

Set $V_{2}(S, E, R) = -\ln S-\ln E-\ln R$. Then, we have

Define

Then, we have

where

Define a $C^{2}$ function $Q:\mathbb{R}_{+}^{4}\rightarrow \mathbb{R}$ in the following form:

where $M > 0$ is sufficiently large and satisfies the condition

where

In addition, $Q(S, E, I, R)$ is continuous, and $(S_{0}, E_{0}, I_{0}, R_{0})$ is a minimum value point of $Q(S, E, I, R)$ in $\mathbb{R}_{+}^{4}$. Therefore, define a $C^{2}$ function

Clearly, $V$ is nonnegative. By the Itô formula and combining $(3.2)$–$(3.4)$, we get

The tectonic compact set is

For the sake of discussion, let's divide $\mathbb{R}_{+}^{4}\backslash D$ into eight regions:

where $\varepsilon_{i}(0 < \varepsilon_{i} < 1, i = 1, 2, 3, 4)$ are positive numbers small enough to satisfy that the following conditions hold

where

In the following, we prove that the eight regions have $\mathcal{L}V(S, E, I, R)\leq-1$ for any $(S, E, I, R)\in D^{c}$.

Case 1. For any $(S, E, I, R)\in D_{1}$, by $(3.8)$, we have

Case 2. On $D_{2}$, by $(3.7)$ and $(3.9)$, we have

Case 3. When $(S, E, I, R)\in D_{3}$, by $(3.7)$ and $(3.10)$, we obtain

Case 4. On $D_{4}$, by $(3.7)$ and $(3.11)$, we get

Case 5. For any $(S, E, I, R)\in D_{5}$, by $(3.7)$ and $(3.12)$, we have

Case 6. On $D_{6}$, by $(3.7)$ and $(3.13)$, we have

Case 7. When $(S, E, I, R)\in D_{7}$, by $(3.7)$ and $(3.14)$, we obtain

Case 8. On $D_{8}$, by $(3.7)$ and $(3.15)$, we get

Thus, for sufficiently small positive numbers $\varepsilon_{i}(i = 1, 2, 3, 4)$, we obtain

Consequently, Theorem 2 holds. □

4.   Extinction of disease

In this section, we will demonstrate that under certain assumptions, the disease will become extinct.

Define a parameter

Theorem 3. Let $(S(t), E(t), I(t), R(t))$ be the solution of system $(1.1)$ with any given initial value $(S(0), E(0), I(0), R(0))\in\mathbb{R}_{+}^{4}$. If ${R}_{0}^{e} < 1,$ then

that is to say, $(E(t), I(t), R(t))$ exponentially converges to $(0, 0, 0)$ a.s.

Proof. Let $p(t) = \varphi E(t)+(\mu+\varphi)I(t)$. By the Itô formula, we obtain

where

By $R_{0}^{e} < 1$, we get

There is $0 < \eta < \min\{\frac{\sigma_{2}^{2}}{2}, \frac{\sigma_{3}^{2}}{2}\}$. Setting $\frac{\bar{\sigma}_{2}^{2}}{2} = \frac{\sigma_{2}^{2}}{2}-\eta$, $\frac{\bar{\sigma}_{3}^{2}}{2} = \frac{\sigma_{3}^{2}}{2}-\eta$, we can get

Combining $(4.3)$ and $(4.4)$, we obtain

By $(4.1)$, $(4.2)$, and $(4.5)$, we get

Integrating both sides of $(4.6)$ from $0$ to $t$ and dividing by $t$, we get

According to $(4.7)$, we have

The upper formula indicates that

According to $(1.2)$, we get ${\lim\limits_{t\rightarrow \infty}R(t) = 0 }\ a.s$. That shows that $(E(t), I(t), R(t))$ exponentially converges to $(0, 0, 0)$ a.s. We complete the proof of Theorem 3. □

5.   Stochastic optimal control

If sustained control is implemented, the processing level will remain at a relatively high level over time. From the previous sections, we conclude that the cost eradicating the disease may be too high. In order to eliminate the disease within a limited time, time-dependent control should be considered.

As in previous studies^[24], using the stochastic maximum principle as in ^[25], we find the characteristics of optimal control problem of model (1.3). Our objective is to minimize both the number of infectious individuals and the cost of treatment control; thus, we establish the following objective function.

where $A, B,$ and $C$, respectively, represent the weights of the relationship between the state variables $E, I,$ and $u$. The control set is given by $\Gamma = \{u$ is measurable and $0\leq u(t)\leq1$, for $t\in[0, t_{1}]\}$. According to the stochastic control theory in the book^[26] of $\emptyset$ksendal, we need to find an optimal control variable $u^{*}(t)$ that minimizes the objective functional when the initial state is $x_{0}$. We define the expectation of the initial state $x_{0}$ as

Let's assume that there is a fixed constant $\bar{u}(t)$ in the deterministic problem that $\bar{u}(t)\leq1$ with $u(t)\leq\bar{u}(t)\ a.s.$ The class of admissible control laws is

In order to obtain a solution of stochastic control, we define the expectation of the system at time $t$ and a fixed value of $x$ as follows:

Now, let's define the value function to be

We now define the control law of minimizing the expected value of $\mathcal{J}_{s}:\Pi\rightarrow \mathbb{R}_{+}$ given by $(5.3)$. The present solution formulated is the solution of the stochastic analogue we now describe for the optimal control problem.

Given the system $(1.2)$ and $\Pi$ as in $(5.2)$ with $\mathcal{J}_{s}$ as in $(5.3),$ find the value of the function

and an objective function

By the following theorem, the optimal control $u^{*}(t)$ is obtained.

Theorem 4. A solution to the optimal control problem presented in problem $(5.2)$ is of the form

Proof. We calculate $\mathcal{L}\mathcal{U}(t)$:

Applying the Hamilton-Jacobi-Bellman theory^[26], the minimum of $(5.4)$ can be obtained as

In order to obtain the optimal control solution, consider the following expression:

Take the partial derivative of $(5.7)$ with respect to $u$ and set it equal to 0. Thus, the equation is obtained,

Considering the bounds of $u$, we can get an expression for $u^{*}(t)$. □

6.   Numerical simulations

In this section, we illustrate the theoretical results with example. By using Milstein's method^[27], the discrete equations of system $(1.3)$ are described by

where $\xi_{1k}, \xi_{2k}, \xi_{3k},$ and $\xi_{4k}(k = 1, 2, \cdots)$ are independent Gaussian random variables subject to $N(0, 1)$, and $\sigma_{i}(i = 1, 2, 3, 4)$ is the intensity of white noise.

We choose $\Lambda = 1.2, \mu = 0.004, \beta_{1} = 0.0134, \beta_{2} = 0.025, \alpha_{1} = 0.09, \alpha_{2} = 0.02, \gamma = 0.015$, $\varphi = 0.019, \nu = 0.02, a = 0.052, b = 0.01$, the initial value $S(0) = 58, E(0) = 15, I(0) = 20, R(0) = 20,$ and the step size $\Delta t = 0.01.$

In Figure 1, we choose $u = 0.66, \sigma_{1} = \sigma_{2} = 0.05, \sigma_{3} = 0.04, \sigma_{4} = 0.1$ to get that $R_{0}^{s} = 1.0029 > 1$, satisfying the condition of Theorem 2. The result of the graph is consistent with our conclusion in Theorem 2.

Figure 2 shows the stochastic epidemic system $(1.2)$ with $u = 0.66, \sigma_{1} = 0.3, \sigma_{2} = 0.49$, $\sigma_{3} = 0.49$, $\sigma_{4} = 0.3$, and we get that ${R}_{0}^{e} = 0.9945 < 1$, which satisfies the conditions of Theorem 3; this is consistent with our conclusion in Theorem 3. When the intensities of white noises are sufficiently large, the disease of the stochastic epidemic system $(1.2)$ is extinct.

Figure 3(a) shows the extinction image when $u$ takes the variable and the other parameters take the same values as in Figure 2, while Figure 3(b) shows the corresponding trend of $u$ varies with time $t$.

7.   Conclusions

The dynamic behaviors of a class of SEIR epidemic models are studied. First, we prove the existence and uniqueness of a global positive solution for the stochastic model. Second, we explore that the positive solution of the model has a stationary distribution, and we investigate the sufficient conditions for the extinction of the stochastic SEIR epidemic system. Furthermore, we aimed to minimize the total cost of infection and treatment expenses by studying optimal control strategies. The existence of optimal solutions is proved by using the stochastic maximum principle, and the dynamic behavior of the model affected by $u$ is studied. It can be found through experiments that the disease becomes extinct faster when $u$ takes variable values than when it takes constant values, and the number of infections has significantly decreased. In addition, we present the trend of $u$ over time $t$ when the disease is extinct. Based on the changing trends, the public health system can dynamically adjust treatment strategies.

It is shown by detailed theoretical analysis that environmental white noise can control the spread of diseases to some extent, and different proportions of control therapies can be used at different times to achieve the purpose of controlling infectious diseases with the least cost. This provides a theoretical basis for the actual control of infectious diseases.

Although we have studied the optimal treatment control of infectious diseases and can implement different proportions of treatment control according to different time periods. It is very difficult to find the optimal control measures because the dynamics of disease transmission are very complex and influenced by many factors. In addition, in practice, measures to control epidemics cannot be singular, and multiple measures should be considered to jointly control the spread of diseases. Therefore, the control effect of implementing multiple measures together should be studied in combination with reality.

Author contributions

Jinji Du: conceptualization, investigation, methodology, formal analysis, original draft preparation; Chuangliang Qin: validation, formal analysis, writing-review and editing; Yuanxian Hui: resources, writing-review and editing. All authors have read and approved the final version of the manuscript for publication.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (No. 11971127), the Key Research Project of Higher School in Henan Province (Nos. 23B110014 and 25B110019) and the Science and Technology Innovation Youth Project of Zhumadian (No. QNZX202310).

Conflict of interest

The authors declare that they have no competing of interests regarding the publication of this paper.