Dynamics analysis and optimal control of SIVR epidemic model with incomplete immunity.

In this paper, we establish an SIVR model with diffusion, spatially heterogeneous, latent infection, and incomplete immunity in the Neumann boundary condition. Firstly, the threshold dynamic behavior of the model is proved by using the operator semigroup method, the well-posedness of the solution and the basic reproduction number ℜ 0 are given. When ℜ 0 < 1 , the disease-free equilibrium is globally asymptotically stable, the disease will be extinct; when ℜ 0 > 1 , the epidemic equilibrium is globally asymptotically stable, the disease will persist with probability one. Then, we introduce the patient's treatment into the system as the control parameter, and the optimal control of the system is discussed by applying the Hamiltonian function and the adjoint equation. Finally, the theoretical results are verified by numerical simulation.

Introduction

The SARS in 2003, the Zika virus (ZIKV) invasion in 2013, the H7N9 avian influenza pandemic, and the emergence of the Dengue virus in the world, these recurrent infectious diseases and various emerging infectious diseases have been challenging modern life and medical standards [1]. For example, COVID-19, which broke out in 2019, is still affecting the world. As of August 24, 2021, the cumulative number of COVID-19 cases and deaths has reached 212,357,898 and 4,439,843. Therefore, how to prevent and control the occurrence and spread of infectious diseases is one of the hot issues today.

From the perspective of mathematics, the study of infectious diseases usually starts according to the transmission mechanism of diseases, which is analyzed by establishing mathematical models. The earliest epidemic model was established by Kermack and Mckendrick. They established the plague susceptibility infection removal model (SIR) [2] in 1927 and the plague susceptible infected susceptible model (SIS) [3] in 1932, respectively. Since the establishment of SIR and SIS models, most of the subsequent research is based on the standard SIR model. Among the existing prevention and treatment methods for infectious diseases, vaccine injection is one of the fast and effective methods. For example, in the prevention and control of COVID-19, vaccine injection can reduce the infection rate of the Delta variant virus to a certain extent. Therefore, an increasing number of researchers take vaccine injection into account in the process of modeling infectious diseases to make the model close to the actual situation.

Among them, Chen et al. in [4] established the susceptibility vaccination–infection isolation–recovery (SVIQR) model and susceptibility–vaccination–infection isolation (SVIQS) model, respectively. The basic reproduction number of the two models was given. Furthermore, the global attractivity and the global asymptotic stability of the solutions were proved by the Lyapunov function method. And the existence of backward bifurcation was also proved. In [5], Kribs-Zaleta and Velasco-Hernández studied a simple SIV model with inoculation and demonstrated the backward bifurcation of solutions to some parameter values. At the same time, complete bifurcation analysis of the model was given under the condition that the vaccine reduces the basic reproduction number. Liu et al. in [6] studied the following SIVR system: Here, λ is the constant update rate of the susceptible host and α is the rate at which susceptible individuals are vaccinated. and are the transmission rates of infected people in contact with susceptible groups and vaccinated groups, respectively. Since vaccinated individuals may have partial immunity during vaccination, it is assumed that . μ is the host mortality per compartment. γ and δ are the recurrence rates of vaccinated people and infected people, respectively. All of these parameters are assumed to be positive. In [6], the authors gave the threshold dynamics of system (1.1) by using the basic reproduction number, showing that reducing the number of infected individuals by vaccination can control the disease. In addition, many researchers have studied infectious disease models with immunization from the perspective of age structure [7] and pulse vaccination [8].

The above models are all established in a homogeneous space environment. However, in practice, the transmission of some diseases is often related to spatial location. For example, the transmission rate of COVID-19 in Asia is different from that in North America. In [9], Wu et al. discussed a class of spatially heterogeneous host-pathogen models. The authors used the basic reproduction number to discuss the threshold dynamic behavior of the models and gave the asymptotic behavior of the models. In [10], a reaction-diffusion model of SVIR infection in a spatially heterogeneous environment was proposed. The authors gave the proof of the extinction and persistence of the disease by giving a basic reproduction number. In [11], the authors established an SIVS epidemic model with a degree-dependent transmission rate and incomplete vaccination on a scale-free network. The global asymptotic stability of the equilibrium and the global attractivity of the unique endemic equilibrium were proved. In addition, the effects of various immunization programs such as unified immunization, target immunization, and acquaintance immunization were studied and compared.

Motivated by the recent development of epidemic modeling, the optimal control problem is often discussed in some cases, optimal control theory is one of the important branches of mathematical optimization, which is often used to study how to find a control for a dynamic system in a period of time to optimize the objective function. Thus we consider two different models based on (1.1). The first is a direct extension of (1.1). A reaction-diffusion SIVR model is established based on spatial heterogeneity with incomplete immunity. The well-posedness of the system is discussed by using the operator semigroup method. At the same time, the global dynamic behavior of the system solution is discussed by analyzing the basic reproduction number. In the second model, as [12, 13], it is assumed that the spread of a disease can be influenced by decision-makers. That is, decision-makers can control the response rate to a certain extent by increasing the treatment ability or the efficiency of drug treatment. Therefore, by further expanding the model, we obtain the control system under the assumption of limited control resources. Considering the progress of medical technology, the targeted treatment for patients will be gradually developed, so we will consider the targeted treatment for patients as a control parameter in the system and discuss its optimal control problem. In addition, we analyze the optimal control of the system by using the Hamilton equation and the adjoint equation. In the proof, we may encounter the following problems:

How to determine the basic regeneration number of the system. In the presence of the diffusion term, it is necessary to select an appropriate method to represent the basic reproduction number , which is an important prerequisite for discussing the dynamic behavior of the system by using as the threshold value.

Can the existence of optimal control be obtained? Because of the existence of diffusion terms, it is difficult to define the adjoint equation and the Hamiltonian function of the control system. At the same time, there are some requirements for the selection of parameters in the numerical simulation.

In view of the above problems, this article is organized as follows. In Sect. 2, an SIVR model with incomplete immunity and spatial heterogeneity is established. Furthermore, the well-posedness of the model is derived, meanwhile, the global existence and global attractiveness of the solution are proved. Section 3 is devoted to studying the threshold dynamic behavior of the system. The extinction or persistence of diseases is analyzed by using the basic reproductive number as the threshold. In Sect. 4, the optimal control of the control system is analyzed by taking the treatment for the patient as the control parameter in the system. Meanwhile, the optimal control problem is discussed by using the Hamiltonian function and adjoint equation. Finally, in Sect. 5, the corresponding results are verified by numerical simulation.

Model formulation and well-posedness

In this paper, the spatial heterogeneity of the spread of infectious diseases and spatial diffusion is considered. In addition, for vaccines, we consider vaccination rates in susceptible individuals and the effectiveness of the vaccine. Based on model (1.1), we can establish the following epidemic model of SVIR with incomplete immunity. The meanings of parameters in the system (2.1) are shown in Table 1.

Table 1

Description of parameters of the model

Parameter	Biological implication
D_i (i = 1,2,3,4)	Diffusion coefficient in susceptibility, infection, vaccination, recovery path
Λ(x)	Recruitment rate of the susceptible host
r(x)	Vaccination coverage rates of susceptible persons
α(x)	Transmission between infected and vaccinated hosts
β(x)	Transmission between infected and susceptible hosts
d_i(x) (i = 1,2,3,4)	Mortality of susceptible, infected, vaccinated, and recovered hosts
γ(x)	Recovery rate of infected persons
K(x)	Half-saturation concentration
η(x)	Effectiveness of vaccine

Remark 1

We considered that there is a vaccine coverage rate for susceptible path S, and unvaccinated susceptible persons will be injected into the infected path with a transmission rate . The susceptible person who has been vaccinated enters the vaccinated compartment, suppose the effectiveness rate of the vaccine to be , and if the vaccine fails, the vaccinator will also be injected into the infected path since the inoculated host has some resistance to the virus after being vaccinated. Thus this propagation process is assumed to obey a half-saturation rate and .

In addition, because does not appear in the first three equations of (2.1), we denote system (2.1) as with the initial value and boundary conditions It is sufficient to determine the dynamics of (2.1). Here, Ω is a smooth bounded region in . Define a Banach space with the supremum norm and . Next, we mainly analyze the dynamic behavior of system (2.2).

Let () be the -semigroup associated with , , . For any , is given by the following formula: and where () is the Green function associated with the operator in Ω̅ subject to the boundary condition. With [14, Section 7], are compact and strongly positive. Set Then we can rewrite (2.2) as the following integral equation: for and .

For a positive and continuous function on Ω̅, define Thus, for the local solution of (2.2), we have the following.

Lemma 2.1

System (2.2) with any initial value ϕ for (where ) has a unique solution with . Moreover, is a classical solution.

The proof is shown in Appendix A.

In the remainder of this section we will prove the global existence and boundedness of the solution. Consider the following equation: where and , are positive and continuous functions on Ω. Thus we have the following.

Lemma 2.2

([9, Lemma 1])

System (2.4) admits a positive steady state which is unique and asymptotically stable. Furthermore, if , are constants, thus .

The following theorem proves the boundedness of the model.

Theorem 2.1

For , , the solution of system (2.2) satisfies where is the semiflow associated with the solution. Moreover, is ultimately bounded.

The proof is shown in Appendix B.

From what has been discussed above, we can get the following results.

Lemma 2.3

The semiflow admits a compact and global attractor.

Proof

With Theorem 2.1 we ensure the ultimate boundedness of system (2.2). Notice that the equation of (2.2) has the diffusion term, which ensures that is compact. With the direct consequence in [15, Theorem 2.4.6], we can complete the proof. □

In Sect. 3, we analyze the threshold dynamics of system (2.2). The global stability of disease-free equilibrium (DFE) and endemic equilibrium (EE) is analyzed by establishing the relationship between the basic reproduction number and the principal eigenvalue.

Threshold dynamics

Basic reproduction number

In this section, applying the methods in [16, Sect. 3], we give the basic reproduction number of (2.2). It is easy to see that (2.2) admits a disease-free equilibrium (DFE) , where , satisfy

Linearizing (2.2) at , we can get the following system: In order to discuss the basic reproduction number , we will focus on the linearized equation of infected person I.

Substituting , we consider the following subsystem:

Thus, define then the next generation operator is defined as . Set as the -semigroup associated with , then

We can use the variational formula to give the solution of the eigenvalue problem as By [16, Theorem 3.2], we have and is the spectral radius of .

Remark 2

If all parameters are all constants, , , we have

From Remark 2, we can see the relationship between and parameters. We have the following lemma on the impact of with .

Lemma 3.1

For the basic reproduction number , we have:

;

For , is a positive and strictly monotonic decline function;

for ;

for ;

For , there exists such that, for , and , ; For , for all .

Define the principal eigenvalue of (3.3) as . Thus, we have the following result.

Lemma 3.2

has the same sign as .

The proof is shown in Appendix C.

Extinction of disease

In this subsection we give the proof of , the stability of DFE.

Theorem 3.1

If , the disease-free equilibrium is globally asymptotically stable. Thus, for ,

Proof

Applying Lemma 3.1, we can infer that for the principal eigenvalue . By the equation of S, V in (2.2), with continuity, there exist and Thus, for , , the eigenvalue problem has a principal eigenvalue . Thus, (3.2) implies By the comparison principle, we can find a positive constant which satisfies where is a strong positive eigenfunction associated with . Since , we directly have Thus the equation of S, V in (2.2) is asymptotic to (3.1), with asymptotic autonomous semiflow theory in [17, Corollary 4.3], such that This completes the proof. □

Disease persistence

In this section, we prove the global asymptotic stability of the endemic equilibrium in the case of . First, using [10, Lemma 4.1], we get the following conclusion, which ensures that (2.2) has a positive epidemic equilibrium.

Lemma 3.3

Let , , be the solution of (2.2) with the initial value ϕ. If there exists a positive such that , thus for any , . Moreover, for S, V, and where , .

The proof is shown in Appendix D.

Next, we conclude this section by proving the global stability of the endemic equilibrium.

Theorem 3.2

For , (2.2) admits at least one positive steady state, and we can find positive ε for any with such that

Define two sets as and With Lemma 3.3, for , we can find that , , which implies and . Set here is an omega limit set.

Claim 1. .

With , we know that , , thus we have . Thus we can find that (2.2) is asymptotic to Then, with Lemma 2.2, on , S, I satisfy and , uniformly. Hence , .

Claim 2. satisfies where is a positive constant, here , will be defined in what follows.

First, by using [16, Lemma 2.2], there exists positive sufficiently small, the following eigenvalue problem admits a principal eigenvalue . To deal with nonlinear terms in the system, we can choose with Then, to the contrary, assume that there is positive such that, for , It follows that there exists for which satisfies Define for , which satisfies Thus, for , is the unique solution of (3.7). ( is the strong positive eigenfunction corresponding to ). For , with Lemma 3.3, it follows that . From the definition of the upper solution, for the solution of (3.8) and , we can find that

By the comparison principle, we can find a small positive constant ς which satisfies Thus, for , with Lemma 3.2, which implies , then when . It means that is unbounded, which contradicts the previous proof. This completes the proof.

Here we give the definition of a distance function in the semiflow as follows: where . Hence, applying [10, Theorem 4.1], for any , we have and By Lemma 3.3, we can set thus . Set , we can get the result that the endemic equilibrium is uniformly persistent. Therefore, by [18, Theorem 4.7], for (2.2), there is at least one positive steady state of on . □

Optimal control

In the previous chapters, we have focused on the disease-free equilibrium and the infectious equilibrium of infectious diseases. But if there is a sudden outbreak, we need to control the impact of the disease at a lower level as far as possible, that is, to control the number of infected people. In addition to calling for vaccinations, governments often spend more money on treatment. The mathematical language to describe this method is the optimal control problem. The main aim of this section is to develop effective strategies for controlling the spread of infectious diseases. We hope that the number of infected people does not exceed the number of susceptible and effective vaccinators.

In this section, we introduce the control strategy to (2.2) and analyze its properties. For convenience, we rewrite as Λ, the same for other parameters. To complete our research, we analyze the control variables of the model (2.2). Therefore, the control variables are given as follows.

With the development of medical technology, infected patients can be treated better. Therefore, define represents the medical intervention for infected patients. Considering that medical resources are limited, we use for specific. Here, c is the cure rate and ω denotes the saturation constant.

From this, we give the control system of (2.2) as follows: with the boundary condition here u is measurable, other parameters are the same as in (2.2).

Define an objective function where and . Assume that the control set is convex, , are weight of each item. This objective function describes our goal to control the problem: to reduce the number of susceptible and infected people with minimal intervention costs. The value function is defined as

Define a Hilbert space and as the initial value of (4.1), which satisfies (IOC: Initial value of Optimal Control) Let , where With the method in [19, 20], we have the following lemma.

Lemma 4.1

For the initial value , , , system (4.1) admits a unique strong solution such that and . Furthermore, there exists a positive constant C independent of u, for any , where .

Lemma 4.1 ensures the existence and boundedness of the global solution of system (4.1). Next, according to [21], we can analyze the existence of optimal control of system (4.1).

Theorem 4.1

Let the initial value be defined in (IOC). Then there exists an optimal solution of the control system (4.1) corresponding to optimal control .

From the boundedness we have proved in Sect. 3, we can infer that is finite. Thus, for and ,we can find a sequence that is the solution to the following subsystem: with the initial condition and the boundary condition such that

With Lemma 4.1, for system (4.2), we can infer that

Since is compactly imbedded in , we can also get the compactness of , , and , , . Here, by using the Arzela–Ascoli theorem [22], for the compactness we proved in Sect. 2, uniformly in with respect to a subsequence denoted by . In addition, with the weak convergence of , , (with the boundedness in system (4.2)), we have weakly in . With (4.3), we have for and

Next, we focus on the second equation of (4.2). By direct calculation we have Thus Similarly, we can discuss the first and third equation of (4.3). For the subsequence , weakly in . With the convexity and closeness of , hence , we have

With the analysis above, we can give the conclusion that for , is the optimal solution associated with the optimal function of system (4.1), which completes the proof. □

The Hamiltonian function of the control system is given as follows: Next we give the adjoint equation for the control system (4.1)

Using the method in [23], give the partial derivative of the Hamiltonian function to u, substitute it into the optimal control solution , Let , the optimal control pair satisfying the minimum value of the objective function can be expressed as

Numerical simulation

In this section, we use numerical simulation to verify the stability of the system and the impact of controls on the development of the disease. The values of each parameter are shown in Table 2.

Table 2

Values for the parameters in numerical simulation

Parameter	Data 1	Data 2	Source
D_1	1.25 × 10−4	1.25 × 10−4	[10]
D_2	1.25 × 10−4	1.25 × 10−4	[10]
D_3	1.25 × 10−4	1.25 × 10−4	[10]
Λ	0.4	0.8	Assume
β	0.6	0.75	[24]
η	0.72	0.62	[25]
α	0.5	0.7	[24], Assume
r	0.4	0.2	Assume
d_1	0.1595	0.1595	[26]
d_2	0.1815	0.2145	Assume
d_3	0.1595	0.1595	[26]
γ	0.9	0.75	[27], Assume
c	0.75	0.75	Assume
ω	0.5	0.5	Assume

Stability of equilibrium

In this section, we discuss the stability of the solution of system (2.2). First of all, use the method in [28] to make difference and solve system (2.2). Then, for and , take the data in Data 1 and Data 2, respectively. We get the simulation results obtained as follows.

Firstly, take the value in Data 1. As we can see from Fig. 1, when , from Theorem 3.1, the disease-free equilibrium is asymptotically stable. In fact, as shown in Fig. 1(b), for , converges to zero; in addition, , . This is the same conclusion given by Theorem 3.1.

Figure 1

, the density of S, I, V

Then, using Data 2, we get the numerical simulation when . With Theorem 3.2, we can get the uniform persistence of disease. Actually, as we can see in Fig. 2(b), with , approaches a constant positive value. This is consistent with our proof in Theorem 3.2.

Figure 2

, the density of S, I, V

Influence of control parameters on disease progression

In this section, we discuss the impact of control on a disease. The main idea of this section is that we solve the optimal control problem by applying the iterative method. Then the optimal system is obtained by using the state equation and adjoint equation given in Sect. 4. And by solving the optimal system, the optimal control strategy is obtained. Furthermore, the method in [28] is used to make difference and solve the control system and the adjoint equation. In order to control the susceptible population, the infected population, and the vaccinated population, targeted treatment of patients is taken as control, and the impact of targeted treatment on the susceptible population, the infected population, and the density of the vaccinated population in a long-term state is considered. Finally, through numerical simulation, the actual situation of each path in the original system (2.2) and the control system (4.1) is compared.

First, we define the objective function as follows. Let the objective function corresponding to the control system (4.1) be as follows: where , [29]. The values of other parameters are the same as in Sect. 5.1. The numerical simulation results are as follows.

In Fig. 3, for , under controlled conditions, the duration of the disease is shorter and the time of extinction is earlier than without control. At the same time, the overall density of susceptible and vaccinated people before reaching the stabilization point was also higher in the controlled condition than without control. In addition, according to Fig. 4, the control intensity reached the maximum in the early stage and gradually decreased with the weakening of the disease scale, reaching zero value when the disease disappeared.

Figure 3

Infectious path S, I, V without and with control ()

Figure 4

The optimal control when

As shown in Fig. 5, for , when control exists, the scale of the disease reaches a minimum earlier than without control, and when the disease eventually becomes endemic, the total scale of the disease is lower than without control. In addition, when the disease reaches a stable state, the density of susceptible and vaccinated persons is higher in the controlled condition than without control. Furthermore, by Fig. 6, when the disease is in its initial state of development, control rises, and when the disease reaches equilibrium and becomes endemic, control is maintained at a stable value along with the duration of the disease.

Figure 5

Infectious path S, I, V without and with control ()

Figure 6

The optimal control when

Conclusion and discussion

In this paper, a kind of SIVR infectious disease model including vaccine immunity and vaccine effectiveness is considered. The optimal control theory is applied to the study of the model, and the threshold dynamics and optimal control of the model are discussed. Firstly, we prove the well-posedness of the model, which provides a theoretical basis for the following discussion. Secondly, we give the basic reproduction number to analyze the dynamic behavior of the disease threshold. In addition, the Hamiltonian function and adjoint equation of the optimal control problem is given. Finally, the stability of the system solution is verified by numerical simulation and the number of infections can be reduced as much as possible, while the cost is reduced under the treatment control. In this paper, the parameters are assumed to be accurate; in fact, due to various uncertainties, each parameter may be inaccurate or random. At the same time, according to the changes in the parameters, it can be seen that the vaccination rate and the effective rate of the vaccine also have a certain impact on the control (see Fig. 7(a), (b)). In addition, because the near-optimal control is more flexible, it can adapt to different degrees of model uncertainty. Therefore, in future work, the near-optimal control problem of the epidemic model can be further studied by considering the influence of random parameters, noise, the vaccination rate, and the efficiency rate of the vaccine as the control parameter.

Figure 7

Optimal control under different parameter values