# Dynamic Analysis of a VSEIR Model with Vaccination Efficacy and Immune Decline

Authors: Xiaotao Han; Hua Liu; Xiaofen Lin; Yumei Wei; Ma Ming
Journal: Advances in Mathematical Physics
 (2022)

---

## Abstract

Vaccination is an effective way to prevent the spread of infectious diseases. In this study, we formulate a VSEIR mathematical model to explore the effects of vaccination rate, vaccine efficacy, and immune decline on the COVID-19 transmission. The existence and stability criteria of equilibrium states were determined by analyzing the model. Model analysis was performed. One of the interesting phenomena involved in this issue is that diseases may or may not die out when the basic reproduction number falls below unity (i.e., a backward bifurcation may exist and cause multistability). The disease eventually becomes endemic in the population when the basic reproduction number exceeds one. By comparing different vaccination rates, vaccine efficacy, and infection rate factors, the diseases can be eliminated, not only by vaccines but also by strict protective measures. In addition, we used the COVID-19 number of reported cases in Xiamen in September 2021 to fit the model, and the model and the reported data were well matched.


---

## Body

## 1. Introduction

Infectious diseases have always posed a severe threat to human health. Many experts and scholars are exploring how to more accurately assess the spread of infectious diseases and assist medical workers in formulating scientific strategies for disease prevention and control. Mathematical models provide scientific tools that researchers can use to understand the causes and patterns of disease transmission and help to identify appropriate control interventions and control policies. Li et al. simulated the third wave of COVID-19 epidemics in Pakistan and proposed that controlling the relevant parameters could significantly reduce the number of infections [1]. Shen et al. established an optimal control model, and the results show that the control strategy can minimize the number of infected individuals [2]. More references will be given in a later discussion. In addition, the theory of fractional differential equations has been developed relatively maturely and applied to many fields [3–10]. He et al. proposed a new fractional-order discrete-time SIR model with vaccines, which showed that there is more complex dynamic system behaviour under the action of vaccines [11]. Altaf Khan et al. demonstrated the global properties of fractional-order Zika models using Lyapunov function theory in a fractional environment [12].Severe acute respiratory syndrome (SARS) in 2003 [13], the 2009 influenza A (H1N1) outbreak [14], Ebola virus disease in 2014 [15], and the COVID-19 pandemic in 2019 [16] each had a substantial impact on human health and on global economic and social behaviour. For emerging infectious diseases, owing to the lack of effective treatment due to insufficient understanding of the pathogenic virus, the world’s people are pinning their hopes on vaccines. Although vaccines have been and continue to be developed through the unremitting efforts of scientists, effective and equitable management poses a considerable challenge.In some developed and developing countries, there is sufficient financial and technical support to dominate the world in the development of new vaccines. In some of the least developed countries, health conditions are relatively poor, and they need vaccines to curb the spread of diseases. However, owing to vaccine research and development technology and public health expenditure constraints, there can be a shortage of vaccine stocks or even no vaccines available. So far, seven COVID-19 vaccines have been put into use, mainly the vaccines produced in China and some European and American countries. The difference is that the vaccines in Europe and American are mainly supplied to rich developed countries, while the Chinese vaccine flows mainly to relatively poor developing countries. High-income countries stocked vaccines while developing countries struggle to get them. Nearly 70 low-income countries worldwide have a vaccination rate of only 1 in 10 people [17]. The international community is also actively cooperating in vaccination programs to promote equitable distribution of vaccines. The problem of uneven distribution of vaccines is difficult to solve in the short term. Therefore, more realistic vaccination rates need to be explored. Some people hope to see more information about the effects of new vaccines, and the public has paid close attention to the data on safety and efficacy. Most people considering vaccination seek better understanding of the vaccine before making an appointment. As the COVID-19 pandemic has progressed, more people have become willing to vaccinate, the scope of vaccination has expanded, and the effects of vaccine use have become more visible. Hence, in modeling infectious diseases, the impact of different vaccination rates on disease transmission should be considered.Most known pathogenic viruses are RNA viruses. Owing to the specific molecular structure of the genetic material of RNA viruses, they readily experience replication errors when reproducing and so produce new viral subtypes. These new subtypes can quickly affect the immune effect of vaccines. According to the World Health Organization, only 2 of the more than 100 human papillomaviruses (HPV) that currently exist (types 16 and 18) cause 70% of cases of cervical cancer, although several other low-risk subtypes can cause other lesions. Three different vaccines have been used against these highly carcinogenic variants [18]. On November 26, 2021, the World Health Organization announced a particular variant of the new coronavirus called Omicron [19]. At that time, there was preliminary evidence of increased transmission and resistance to vaccines and an increased risk of reinfection. Since then, these traits have been confirmed. Therefore, the impact of vaccination on disease transmission cannot be assessed without considering viral mutations. After vaccination, the human body produces large amounts of antibodies, which can prevent the virus from multiplying after exposure. However, the body’s natural metabolism and apoptosis can cause the antibody levels to decline. The immune effect on the disease is therefore weakened.The vaccine’s effectiveness cannot reach 100% [20]; the infection is temporary. Studies have shown that six weeks after the COVID-19 vaccine, the antibody levels in the vaccinated people begin to decline, and within 10 weeks, they drop by more than 50%. Regardless of the vaccinated person’s age, chronic disease status, or sex, the decline in antibody levels remains consistent across all populations [21]. That is, the immunity gained from either vaccines or recovery from infection is temporary. Therefore, when simulating the disease transmission process, both vaccine-induced immunity loss and the disappearance of natural immunity must be taken into consideration.Vaccines are the most effective preventive measure to control the development of the disease, but they are affected and limited by many factors in the implementation process. Scholars in many countries have tried to find the best vaccination strategies for vaccines from different perspectives and have published many practical research results. For example, Zheng et al. established an extended SEAIRD partition model to describe the epidemic dynamics of single- and double-dose vaccination [22]. Deng et al. and Chen et al. produced a Filippov epidemic model to consider the saturation function of vaccination and treatment strategies. They explored media effects on the prevalence of vaccine-containing diseases [23, 24]. Acuña-Zegarra et al. considered secondary infections after vaccination, established a SIRSI model, and analyzed the dynamic behaviour of the disease at immunization [25]. Omame et al. established a gender-sensitive HPV model that studied and rigorously analyzed the effects of treatment and vaccination on their transmission dynamics [26]. Vaccination strategies cannot be developed without taking into account vaccine availability, vaccination rates, efficacy, and the actual situation of loss of natural immune response and viral mutation [27–29]. In addition, Parsamanesh et al. established a series of epidemiological models on vaccination and theoretically demonstrated the local and global stability of the disease-free equilibrium point and the endemic equilibrium point, as well as various possible bifurcations [30–35].While many researchers have studied the effects of vaccination on disease transmission, one of the vaccinations and vaccine decline was sometimes overlooked in previous models. So, we consider improving existing mathematical models to model the impact of vaccination and immunity decline (the decline of natural and vaccine-guided immunity due to viral mutations over time) on vaccine protection and disease transmission.The organization of the paper is as follows. Section2 presents a mathematical expression for an epidemiological model with linear vaccination rates considering both immune decline and viral variation. In the next section, the positive and bounded nature of the model solution is demonstrated. The presence and basic reproduction number of disease-free equilibrium points and the uniqueness of endemic equilibrium points are given in Section 4. The local stability of equilibrium points is discussed in Section 5. In Sections 6 and 7, mathematical proofs of the existence of forward and backward bifurcation and related numerical simulations are given, respectively. In the final section, we summarize the conclusions and discuss further work.


## 2. Model

In this paper, we refine traditional vaccination models, taking into account the effectiveness of vaccines and the loss of immunity (vaccine-induced and natural). An epidemiological model with vaccination was proposed to reflect the epidemiological dynamics of epidemic transmission with preventive measures. We consider dividing the population into five segments: vaccinated susceptibleVt, unvaccinated susceptible St, exposed but incapacitated Et, symptomatic and capable of infection It, recovered Rt.In vaccinated populations, due to the immune response caused by vaccination, these individuals cannot be infected with the virus, let alone transmit it. Over time, vaccine-induced antibodies gradually decrease, and some individuals return to the susceptible population and may be infected by the virus. Susceptible people become exposed by coming into contact with infected people. Here, an exposed person indicates that they have been infected but cannot infect others. These individuals are in the incubation period, and when they pass the incubation period, they become symptomatic infected people. The infected population is transferred to the recovering population through treatment, but because natural immunity also disappears over time, a subset of individuals in the recovered population become susceptible again. A path map of disease transmission is given in Figure1.Figure 1
Flow diagram of disease transmission routes.The mathematical form of the corresponding VSEIR model is expressed by the following ODES nonlinear system:(1)dVtdt=γSt−1−θβVtIt−δ+dVt,dStdt=Λ+δ1Rt+δVt−βStIt−γ+dSt,dEtdt=βStIt+1−θβVtIt−ε+dEt,dItdt=εEt−d+α+γ1It,dRtdt=γ1It−δ1+dRt.The total populationNt is
(2)Nt=Vt+St+Et+It+Rt.The initial condition for model (1) takes the form
(3)V0>0,S0>0,E0>0,I0>0,R0>0.Here,γ represents the maximum vaccination rate of susceptible people. Parameter θ indicates the effective rate of vaccine; then, 1−θ suggests the probability of vaccination failure. We know that both vaccine-induced immunity and natural immunity will gradually disappear so that δ and δ1 are vaccine-induced immunity and natural immunity loss rate, respectively. ε represents the progression rate from incubation to infection. The model parameter γ1 indicates the treatment rate of infected individuals. Parameters d and α are natural mortality and disease-related mortality rate, respectively. β is the contact rate between susceptible people and infected people. Λ indicates the recruitment of susceptible population rate.


## 3. The Positive and Bounded of the Solution

This section is the fundamental property of model (1). When t>0, it is necessary to prove that the solution of model (1) under the condition that a nonnegative initial value is always nonnegative. The following theorem shows the positive and bounded nature of model (1) solution.Theorem 1 (positive).
For allt≥0, the solutions Vt, St, Et, It, and Rt of model (1) with initial condition (3) are positive.Proof of Theorem 1.
Define(4)Mt=minVt,St,Et,It,Rt,
whereVt, St, Et, It, and Rt are any positive solution of model (1) with initial condition (3). It is clear that M0>0. Assuming that there exists a t1>0 such that Mt1=0 and Mt>0 for all t∈0,t1. If Mt1=St1, then Vt≥0, Et≥0, It≥0, and Rt≥0 for all t∈0,t1.
It follows from the first equation of model (1) that
(5)dStdt≥Λ−βIt+γ+dSt,t∈0,t1.
Both sides of the equation are multiplied byexp∫0tβIt+γ+dds at the same time, which can be rewritten as
(6)ddtStexp∫0tβIt+γ+dds≥Λexp∫0tβIt+γ+dds.
Therefore, we get(7)St1≥S0exp−∫0t1βIt+γ+dds+exp−∫0t1βIt+γ+dds×Λ∫0t1exp∫0wβIs+γ+ddsdw>0,which leads to a contradiction. Therefore, we obtain St>0 for all t>0. Similar can get Vt>0,Et≥0,It≥0, and Rt≥0.Theorem 2 (bounded).
All solutions of model (1) with initial condition (3) is positively invariant in the region Ω.Proof of Theorem 2.
Add all equations in model (1) to get
(8)dNtdt≤Λ−dNt.We can get
(9)limsupt⟶∞Nt≤Λd.
Let(10)Ω=V0,S0,E0,I0,R00≤Vt+St+Et+It+Rt≤Λd.
Proof is completed.


## 4. The Existence of Equilibrium Points and the Basic Reproduction Number

### 4.1. Disease-Free Equilibrium Point and Basic Reproduction Number

Let the right side of the equation of model (1) be all zero, and we get the disease-free equilibrium point E0=V0,S0,0,0,0,where(11)V0=Λγdd+γ+δ,S0=Λδ+ddd+γ+δ.We give the basic reproduction numberR0 of model (1) using the method of the next-generation regeneration matrix [36].LetX=E,I,R,V,S; then, model (1) can be rewritten as
(12)dXdt=FX−VX,where
(13)FX=βSI+1−θβVI0000,VX=ε+dE−εE+d+α+γ1I−γ1I+δ1+dR−γS+1−θβVI+δ+dV−Λ−δ1R−δV+βSI+γ+dS.The derivatives ofFX and VX are estimated at the disease-free equilibrium point E0 as follows:
(14)DFE0=F000,DVX=V0J3J4.Here,(15)F=0βS0+1−θβV000,V=ε+d0−εd+α+γ1,(16)J3=0−γ101−θβV00βS0,J4=δ1+d000δ+d−γ−δ1−δγ+d.Therefore, the basic reproduction number of model (1) is determined by the spectral radius of FV−1, namely,
(17)R0=ρFV−1=βεS0+1−θV0ε+dd+α+γ1=Λβεδ+d+γ1−θdd+γ+δε+dd+α+γ1.Here,1/d+α+γ1 indicates the average duration of illness for an infected individual; ε/ε+d is the probability that an infected individual will survive the incubation period. The expression βεS0+1−θV0/ε+dd+α+γ1 actually represents the number of people an infected individual can infect a new patient during infection. So, when R0<1, the disease may go extinct, just maybe.


### 4.2. Uniqueness of the Existence of Endemic Equilibrium Points

By model (1) can get E∗=d+α+γ1/εI∗ and R∗=γ1/δ1+dI∗. The first and third forms of model (1) can be combined to obtain
(18)V∗=γε+dd+α+γ1βεδ+d+γ1−θ+β2ε1−θI∗,(19)S∗=ε+dd+α+γ1βε−γ1−θε+dd+α+γ1βεδ+d+γ1−θ+β2ε1−θI∗,whereI∗ is the positive root of the following equation:
(20)AI∗2+BI∗+C=0,where(21)A=δ1γ1β2ε1−θ−δ1+dε+dd+α+γ1β21−θ<0,(22)B=DR0R0−M,(23)C=dδ1+dε+dd+α+γ1γ+d+δR0−1,in which(24)D=d2d+γ+δ2ε+d2d+α+γ121−θδ1+dΛεδ+d+γ1−θ2,(25)M=δ+d+γ1−θd+γ+δ+δ+d+γ1−θ2ε+dd+α+γ1δ1+d−εγ1δ1dd+γ+δε+dd+α+γ11−θδ1+d>0.Theorem 3.
ForR0>1, model (1) has a unique endemic equilibrium point E1=V∗,S∗,E∗,I∗,R∗.Proof of Theorem 3.
Obviously,A<0 always holds, when R0>1 has C>0. From the relationship between the roots of the unary quadratic equations and the coefficients, we can know I1∗I2∗=C/A<0; therefore, when R0>1, equation (20) has a unique positive root E1=V∗,S∗,E∗,I∗,R∗. Proof is completed.Remark 4.
WhenR0>1, the endemic equilibrium point position of different situations is different, but there is always only one endemic equilibrium point.Equation (20) is a quadratic equation, in which the signs of the coefficients B and C can change the case of the root of the equation. The sign of the coefficient B is mainly determined by the relationship between R0 and M; obviously, M>0; we divide M into M>1, M=1, and M<1 three cases. We summarize the case of the positive root of equation (20) in Table 1, and the corresponding schematic diagram is shown in Figure 2 and M changes. Figure 2(a) is the case when M>1, Figure 2(b) is the case when M=1, and Figure 2(c) is the case when M<1.Table 1
Equation (4). The sign change of factors B and C and the number of possible equilibrium points of endemic diseases.
CasesCompareM and R0CoefficientsNo. of possible positive rootsBCM>1R0>M++11<R0<M−+1R0<1−−0M=1R0>1++1R0<1−−0M<1R0>1++1M<R0<1+−2R0<M−−0Figure 2
Schematic diagram of the symbolic change ofB when the relationship between R0 and M changes. (a) is the case when M>1, (b) is the case when M=1, and (c) is the case when M<1.
(a)(b)(c)


## 4.1. Disease-Free Equilibrium Point and Basic Reproduction Number

Let the right side of the equation of model (1) be all zero, and we get the disease-free equilibrium point E0=V0,S0,0,0,0,where(11)V0=Λγdd+γ+δ,S0=Λδ+ddd+γ+δ.We give the basic reproduction numberR0 of model (1) using the method of the next-generation regeneration matrix [36].LetX=E,I,R,V,S; then, model (1) can be rewritten as
(12)dXdt=FX−VX,where
(13)FX=βSI+1−θβVI0000,VX=ε+dE−εE+d+α+γ1I−γ1I+δ1+dR−γS+1−θβVI+δ+dV−Λ−δ1R−δV+βSI+γ+dS.The derivatives ofFX and VX are estimated at the disease-free equilibrium point E0 as follows:
(14)DFE0=F000,DVX=V0J3J4.Here,(15)F=0βS0+1−θβV000,V=ε+d0−εd+α+γ1,(16)J3=0−γ101−θβV00βS0,J4=δ1+d000δ+d−γ−δ1−δγ+d.Therefore, the basic reproduction number of model (1) is determined by the spectral radius of FV−1, namely,
(17)R0=ρFV−1=βεS0+1−θV0ε+dd+α+γ1=Λβεδ+d+γ1−θdd+γ+δε+dd+α+γ1.Here,1/d+α+γ1 indicates the average duration of illness for an infected individual; ε/ε+d is the probability that an infected individual will survive the incubation period. The expression βεS0+1−θV0/ε+dd+α+γ1 actually represents the number of people an infected individual can infect a new patient during infection. So, when R0<1, the disease may go extinct, just maybe.


## 4.2. Uniqueness of the Existence of Endemic Equilibrium Points

By model (1) can get E∗=d+α+γ1/εI∗ and R∗=γ1/δ1+dI∗. The first and third forms of model (1) can be combined to obtain
(18)V∗=γε+dd+α+γ1βεδ+d+γ1−θ+β2ε1−θI∗,(19)S∗=ε+dd+α+γ1βε−γ1−θε+dd+α+γ1βεδ+d+γ1−θ+β2ε1−θI∗,whereI∗ is the positive root of the following equation:
(20)AI∗2+BI∗+C=0,where(21)A=δ1γ1β2ε1−θ−δ1+dε+dd+α+γ1β21−θ<0,(22)B=DR0R0−M,(23)C=dδ1+dε+dd+α+γ1γ+d+δR0−1,in which(24)D=d2d+γ+δ2ε+d2d+α+γ121−θδ1+dΛεδ+d+γ1−θ2,(25)M=δ+d+γ1−θd+γ+δ+δ+d+γ1−θ2ε+dd+α+γ1δ1+d−εγ1δ1dd+γ+δε+dd+α+γ11−θδ1+d>0.Theorem 3.
ForR0>1, model (1) has a unique endemic equilibrium point E1=V∗,S∗,E∗,I∗,R∗.Proof of Theorem 3.
Obviously,A<0 always holds, when R0>1 has C>0. From the relationship between the roots of the unary quadratic equations and the coefficients, we can know I1∗I2∗=C/A<0; therefore, when R0>1, equation (20) has a unique positive root E1=V∗,S∗,E∗,I∗,R∗. Proof is completed.Remark 4.
WhenR0>1, the endemic equilibrium point position of different situations is different, but there is always only one endemic equilibrium point.Equation (20) is a quadratic equation, in which the signs of the coefficients B and C can change the case of the root of the equation. The sign of the coefficient B is mainly determined by the relationship between R0 and M; obviously, M>0; we divide M into M>1, M=1, and M<1 three cases. We summarize the case of the positive root of equation (20) in Table 1, and the corresponding schematic diagram is shown in Figure 2 and M changes. Figure 2(a) is the case when M>1, Figure 2(b) is the case when M=1, and Figure 2(c) is the case when M<1.Table 1
Equation (4). The sign change of factors B and C and the number of possible equilibrium points of endemic diseases.
CasesCompareM and R0CoefficientsNo. of possible positive rootsBCM>1R0>M++11<R0<M−+1R0<1−−0M=1R0>1++1R0<1−−0M<1R0>1++1M<R0<1+−2R0<M−−0Figure 2
Schematic diagram of the symbolic change ofB when the relationship between R0 and M changes. (a) is the case when M>1, (b) is the case when M=1, and (c) is the case when M<1.
(a)(b)(c)


## 5. Local Stability of the Equilibrium Point

In this section, we give the stability proof of model (1) about the disease-free equilibrium point E0 and the endemic equilibrium point E1.Theorem 5.
IfR0<1, the disease-free equilibrium point E0 of model (1) is locally asymptotically stable, and when R0>1, it is unstable.Proof of Theorem 5.
The Jacobian matrix atE0 is given as
(26)JE0=−δ+dγ0−1−θβV00δ−γ+d0−βS0δ100−ε+dβS0+1−θβV0000ε−d+α+γ10000γ1−δ1+d.
The characteristic equation is(27)λ+δ1+dλ+d+δ+γλ+dλ2+l1λ+l2=0,and the eigenvalue is λ1=−δ1+d, λ2=−d+δ+γ, λ3=−d, and λ4,5=1/2−l1±l12−4l2, where l1=ε+α+γ1+2d and l2=1−R0ε+dd+α+γ1. Therefore, when R0<1, all roots of the characteristic equation are with negative real part, and the disease-free equilibrium point E0 is locally asymptotically stable; when R0>1, λ5=1/2−l1+l12−4l2>0. At this time, the disease-free equilibrium point E0 is unstable. Proof is completed.
Next, forR0>1, we analyze the local stability of the unique endemic equilibrium point of model (1), and the Jacobian matrix corresponding to model (1) at E1=V∗,S∗,E∗,I∗,R∗ is
(28)JE1=−1−θβI∗−δ+dγ0−1−θβV∗0δ−βI∗−γ+d0−βS∗δ11−θβI∗βI∗−ε+dβS∗+1−θβV∗000ε−d+α+γ10000γ1−δ1+d.
The characteristic equation is(29)λ5+h1λ4+h2λ3+h3λ2+h4λ1+h5=0.
Here,(30)h1=ε+α+γ+δ+γ+δ+5d+2−θβI∗,(31)h2=2d+α+γ1+εδ1+d+γ+δ+2d+2−θβI∗ε+α+γ1+δ1+3d+1−θβI∗+δ+dβI∗+γ+d−δγ,(32)h3=γ+δ+2d+2−θβI∗2d+α+γ1+εδ1+d+1−θβI∗+δ+dβI∗+γ+d−δγε+α+γ1+δ1+3d+βI∗βεS∗+1−θ2V∗,(33)h4=1−θβI∗+δ+dβI∗+γ+d−δγ2d+α+γ1+εδ1+d+β2εS∗I∗1−θβI∗+δ+d+γ1−θ+β2ε1−θ2V∗I∗δ1+γ+2d+βI∗+β2εδ1−θV∗I∗+βεI∗βS∗δ1+d−δ1γ1,(34)h5=βεI∗βS∗δ1+d−δ1γ11−θβI∗+δ+d+γ1−θ+β2εδ1+d1−θV∗I∗δ+1−θγ+d.
It is easy to know thath1, h2, and h3 are all positive, and if S∗≥δ1γ1/βδ1+d, then h4,h5>0.
If there isΔ2>0 and Δ3>0, then the Routh–Hurwitz criterion shows that all roots of equation (29) have negative real part. That is, the unique endemic equilibrium point E1 is locally asymptotically stable, where
(35)Δ2=h11h3h2,Δ3=h110h3h2h1h5h4h3.
The mathematical expression ofΔ3 here is very complex and is not easy to derive directly. Nevertheless, we give a numerical example that shows that Δ3>0 does hold. Further discussion is summarized in the following theorem.Theorem 6.
IfΔ3>0 and S∗≥δ1γ1/βδ1+d, then the unique endemic equilibrium point E1 of model (1) is locally asymptotically stable.Proof of Theorem 6.
AssumingΔ3>0, now let us prove Δ2>0, and before that, for the convenience of proof, the following two conditions are already valid.
(1)
From equations (18) and (19), we can obtain βεS∗+1−θV∗=ε+dd+α+γ1(2)
By0<θ<1 gets βS∗+1−θβV∗>βS∗+1−θ2βV∗(36)Δ2=h1h2−h3=ε+α+γ1+δ1+3d2d+α+γ1+εδ1+d+ε+α+γ1+δ1+3d2γ+δ+2d+2−θβI∗+1−θβI∗+δ+dβI∗+γ+dε+α+γ1+δ1+3d−δγε+α+γ1+δ1+3d+γ+δ+2d+2−θβI∗2d+α+γ1+εδ1+d+γ+δ+2d+2−θβI∗2ε+α+γ1+δ1+3d+1−θβI∗+δ+dβI∗+γ+dγ+δ+2d+2−θβI∗−δγγ+δ+2d+2−θβI∗−γ+δ+2d+2−θβI∗2d+α+γ1+εδ1+d−1−θβI∗+δ+dβI∗+γ+d−δγε+α+γ1+δ1+3d−βI∗βεS∗+1−θ2V∗>ε+α+γ1+δ1+3d2d+α+γ1+εδ1+d+ε+α+γ1+δ1+3dγ+δ+2d+2−θβI∗h1+1−θβI∗+δ+dβI∗+γ+d−δγγ+δ+2d+2−θβI∗−βI∗βεS∗+1−θV∗=ε+α+γ1+δ1+3d2d+α+γ1+εδ1+d+ε+α+γ1+δ1+3dγ+δ+2d+2−θβI∗h1+1−θβI∗+δ+dβI∗+γ+d−δγγ+δ+2d+2−θβI∗−βI∗ε+dd+α+γ1>0
Thus, all roots of the eigen equation (29) have negative real part, and the unique endemic equilibrium point E1 of model (1) is locally asymptotically stable. Proof is completed.


## 6. Bifurcation Analysis atR0=1

Now, we focus on the dynamic behaviour of model (1) at R0=1. To do this, we choose β as the bifurcation parameter, and R0=1 gets the corresponding threshold of
(37)β∗=dd+γ+δε+dd+α+γ1Λεδ+d+γ1−θ.We use the central manifold theorem [36, 37] to explore the behaviour of model (1) in the vicinity. Firstly, let V=x1, S=x2, E=x3, I=x4, R=x5, dVt/dt=f1,dSt/dt=f2,dEt/dt=f3, dIt/dt=f4, and dRt/dt=f5, so X=x1,x2,x3,x4,x5T. Then, model (1) can be rewritten as
(38)dXdt=f1,f2,f3,f4,f5T.Obviously, whenR0=1, the Jacobian matrix JE0 corresponding to the disease-free equilibrium point E0 has a zero eigenvalue at E0,β∗ and the remaining eigenvalues have negative real parts. The right eigenvector and left eigenvectors of model (1) at E0,β∗ are W=w1,w2,w3,w4,w5T and V=v1,v2,v3,v4,v5, respectively. Over here,
(39)w1=γ2d+δ+γδ+d−d+δ+γ1−θγdd+δ+γ1−θ−γδ+d+εδ1γ1γd+δ+γ1−θdd+δ+γδ1+d1−θε+dd+α+γ1,(40)w2=γd1−θ−δ+dd+δ+γ1−θdd+δ+γ1−θ+εδ1γ1δ+dd+δ+γ1−θdd+δ+γδ1+d1−θε+dd+α+γ1,(41)w3=d+δ+γ1−θ1−θε+d,(42)w4=εd+δ+γ1−θ1−θε+dd+α+γ1,(43)w5=εγ1d+δ+γ1−θ1−θε+dd+α+γ1δ1+d,(44)v1=v2=v5=0,v4=1,v3=εε+d.Now, we calculate the second-order partial derivatives offkk=1,2,⋯,5 with respect to variables V,S,E,I,R and β at E0,β∗. We have
(45)∂2f1∂V∂I=−1−θβ∗,∂2f1∂I∂β=−1−θV0,∂2f2∂S∂I=−β∗,∂2f2∂I∂β=−S0,∂2f3∂V∂I=1−θβ∗,∂2f3∂S∂I=β∗,∂2f3∂I∂β=S0+1−θV0, and others are zero.Apply the central manifold theorem to determine the dynamic behaviour of the system. According to the literature [36–38], there are the following two formulas:
(46)a=∑k,i,j=15vkwiwj∂2fk∂xi∂xj0,β∗,b=∑k,i=15vkwi∂2fk∂xi∂β∗0,0.In this case, we have(47)a=v3w2w4∂2f3∂S∂I+v3w1w4∂2f3∂V∂I=v3w4β∗w2+1−θw1=v3w4β∗1−M,(48)b=v3w4∂2f3∂I∂β∗=v3w4S0+1−θV0>0.According to Equation (23), Equation (24), and Lemma 3 in reference [30], it can be known that the local dynamics of model (1) around E0 are determined by a and b. So, for b>0, if there is M>1, then a<0; if there is M<1, then a>0. The above discussion is summarized as follows.Theorem 7.
WhenM>1 (define M as Equation (25)), model (1) goes through forward (transcritical) bifurcation, and when M<1, it produces backward bifurcation.We chooseα=0.04, γ=0.08, γ1=0.094, δ=0.0001, δ1=0.7, ε=0.00015, θ=0.95, d=0.0071, Λ=10, and β∈0.0255,0.06. The intuitive results of the bifurcation situation are programmed by the MATLAB 2019b software. In Figure 3(a), when R0 exceeds unity, the stable disease-free equilibrium point becomes unstable and a stable endemic equilibrium point appears; that is, a forward (transcritical) bifurcation appears. We choose α=0.0001, γ=0.08, γ1=0.094, δ=0.005, δ1=0.7, ε=0.9, θ=0.95, d=0.00071, Λ=2, and β∈0.00025499,0.00035. In Figure 3(b), as R0 changes, when R0<1, a stable disease-free equilibrium point, an unstable endemic equilibrium point, and a stable endemic equilibrium point appear. When R0>1, there is an unstable disease-free equilibrium point and a stable endemic equilibrium point; that is, the backward bifurcation appears, at which R0<1 cannot be used as a threshold for judging the disappearance of the disease.Figure 3
(a) Forward (transcritical) bifurcation and (b) backward bifurcation.
(a)(b)


## 7. Numerical Simulation

In this section, we mainly use the MATLAB 2019b to draw numerical results. We present some different parameter combinations to give detailed numerical simulations of the various situations in Table1 and the stability of the disease-free equilibrium point.Case 1.R0<1.
We chooseα=0.0001, β=0.004, γ=0.1, γ1=0.094, δ=0.005, δ1=0.07, ε=0.000055, θ=0.8, d=0.0071, and Λ=5. A stable disease-free equilibrium point E0 is obtained as shown in Figure 4, at which point the disease will eventually perish over time, regardless of the initial size of the disease.Case 2.R0>1 and M>1.
We chooseα=0.0001, β=0.00125, γ=0.2, γ1=0.094, δ=0.002, δ1=0.08, ε=0.0055, θ=0.9, d=0.0071, and Λ=20. At this time, M=3.7189, R0=2.1136, Δ2=0.0402>0, and Δ3=1.692e−04>0. We have stable endemic equilibrium point E1=0.9258,0.0929,1.6145,0.0877,0.0947×1000 when 1<R0<M, as shown in Figure 5.Figure 4
Stable disease-free equilibrium pointE0.Figure 5
Stable disease-free equilibrium pointE1.We chooseα=0.0001, β=0.00525, γ=0.2, γ1=0.094, δ=0.002, δ1=0.08, ε=0.0055, θ=0.9, d=0.0071, and Λ=20. At this time, M=3.7189, R0=8.8773, Δ2=0.2992>0, and Δ3=0.0063>0. We have stable endemic equilibrium points E1′=0.0901,0.0352,2.4167,0.1313,0.1417×1000 when R0>M; it is shown in Figure 6.Case 3.R0>1 and M=1.
We chooseα=0.00011, β=0.00004, γ=0.008, γ1=0.0094, δ=0.0001, δ1=0.007364, ε=0.025, θ=0.95, d=0.00071, and Λ=2.5. At this time, R0=1.8405, Δ2=2.1585e−4>0, and Δ3=5.9043e−9>0. We have stable endemic equilibrium points E1′′=0.6713,0.2292,0.3927,0.9607,1.1184×1000 when R0>M, as shown in Figure 7.Case 4.R0>1 and M<1.
We chooseα=0.0001, β=0.00088, γ=0.08, γ1=0.094, δ=0.005, δ1=0.7, ε=0.3, θ=0.95, d=0.00071, and Λ=1.5. At this time, M=0.4816, R0=2.2163, Δ2=3.2675>0, and Δ3=1.2991>0. We have stable endemic equilibrium points E1′′′=0.1406,0.1010,0.3717,1.1761,0.1578×1000 when R0>M; it is shown in Figure 8.Case 5.R0<1 and M<1.
We chooseα=0.0001, β=0.00028, γ=0.08, γ1=0.094, δ=0.005, δ1=0.7, ε=0.9, θ=0.95, d=0.00071, and Λ=2. At this time, M=0.4333, R0=0.9417, Δ2=1.5240>0, Δ2′=2.5229>0, Δ3=0.1227>0, and Δ3′=0.7508>0. We have a stable endemic equilibrium points E1′′′′=1.1331,0.2822,0.1070,1.0154,0.1362×1000 and a unstable endemic equilibrium points E1′′′′′=2.4831,0.2147,1.1331,0.0863,0.0116×1000, when M<R0<1, as shown in Figure 9. At this time, R0<1 alone does not allow the disease to eventually die.Figure 6
Stable disease-free equilibrium pointE1′.Figure 7
Stable disease-free equilibrium pointE1″.Figure 8
Stable disease-free equilibrium pointE′1″.Figure 9
Stable endemic equilibrium pointE1′′′′ and unstable endemic equilibrium point E1′′′′′.IfR0>1, that is, a patient is infected with an average of more than one new disease during the period of infection, regardless of the initial size of the disease, the disease will eventually form endemic over time. We choose α=0.0001, β=0.005,0.0025,0.001, γ=0,0.2,0.5,0.95, γ1=0.094, δ=0.002, δ1=0.08, ε=0.0055, θ=0.7,0.9, d=0.0071, and Λ=20, to discuss the impact of exposure, vaccination rates, and vaccine efficacy on disease transmission when the disease progresses to endemic.When there is a low vaccine efficacyθ=0.7, as shown in Figure 10(a), if the protective measures (social distancing, wearing masks, etc.) are not strictly implemented, vaccination alone will not keep the number of patients low. As the protective measures become more stringent, the proportion of vaccinations will gradually increase, as shown in Figures 10(b) and 10(c). Then, there will be a significant decline in the number of people who will eventually be sick. When there is a high vaccine efficacy θ=0.9, as shown in Figure 11(a), if the protective measures (social distancing, wearing masks, etc.) are not strictly implemented, the rate of infection will still remain considerable. Only vaccinations are given at this time; although the number of patients will decline, it will remain high; as the protective measures become more stringent, the proportion of vaccinations will gradually increase, as shown in Figures 11(b) and 11(c). Then, there will be a significant decrease in the number of people who will eventually become sick, especially when the number of people affected decreases significantly. In the case of vaccination, it is also necessary to enact strict protective measures to drive the disease to extinction.Figure 10
Effects of vaccination rates and exposure rates when vaccines are less effective.
(a)(b)(c)Figure 11
Effects of vaccination rates and exposure rates when vaccines are more effective.
(a)(b)(c)To verify the plausibility of the model, we utilize the adaptive combination of Delayed Rejection and the Adaptive Metropolis (DRAM) algorithm to perform the Markov Chain Monte Carlo (MCMC) process, which is superior to the original MCMC method without good prior distribution [39, 40]. Reference [41] provides the MCMC toolbox. We ran the algorithm 10,000 times using MATLAB 2019b to fit the entire parametric simulation process based on case data reported in Xiamen, Fujian Province, from September 10, 2021, to September 29, 2021, which are derived from daily government reports. Before fitting, we first used the statistics of the Seventh National Census Bulletin of Xiamen Municipal Bureau of Statistics [42] to obtain the recruitment rate of susceptible populations Λ=31.7808/day, natural mortality d=0.0146 [43], vaccine-induced immune decline δ=0.0012396 [25], and natural immune decline δ1=0.00273973 [25]. Since data related to vaccination are difficult to assess, we assume that the vaccination rate γ=0.5 and the vaccine effective rate θ=0.5. Now, fit the model with the reporting data.First of all, six days of data were used to predict the development trend of the disease. As shown in Figure12, the development trend of the disease is roughly the same as the actual situation, but the time of no new cases lags behind the data reported by at least five days. If the fitted data is increased by two days and four days, as shown in Figures 13 and 14, the disease development trend can be better matched with the reported data, and the time of no new cases is advanced or lagged by two days. The prediction of the epidemic is acceptable. In addition, the fitting curve did not show a significant decline on days 5 and 6 because the model did not consider how changes in human behaviour would affect the spread of the disease; i.e., people would consciously and autonomously take precautions to avoid infection after the epidemic.Figure 12
6 days of data fit and prediction for the model.Figure 13
8 days of data fit and prediction for the model.Figure 14
10 days of data fit and prediction for the model.


## 8. Discussion

We here present a mathematical model of VSEIR showing declines in vaccine efficacy and immunity. The model considers differences in vaccination rates due to people’s willingness to vaccinate and the constraints of national economic strength and scientific research conditions.When measuring the impact of vaccines on disease transmission, vaccination rates must not be ignored. No vaccine’s effectiveness can reach 100%. Viral mutations can also affect the protective effect of vaccines. Therefore, the rate of effectiveness of vaccines must not be disregarded in the study of the spread of infectious diseases.The antibodies produced after vaccination gradually decrease in concentration with the person’s metabolic activity, and when the antibody level reaches the resistance threshold against the virus, immunity has been lost. That is, the individual may be reinfected.Based on this, considering vaccination and immunity decline, we analyzed the VSEIR model. Usually, in infectious disease models, the disease goes extinct whenR0<1. But our result is that when R0<1, the disease may not become extinct. In dynamic systems, this result manifests itself as a backward bifurcation when R0 passes through 1 under certain conditions. The main reason for this result is that people who are vaccinated have less awareness of protection, which increases the risk of infection. In the actual epidemic prevention process, even if individuals are vaccinated, active protective measures should be taken to reduce the risk of infection. Furthermore, when R0>1, no matter how large the initial size of the disease is, it eventually becomes endemic. Under the same parameters, through the comparison of different exposure rates, vaccination rates, and vaccine efficiency, results have shown that vaccination alone cannot control the spread of the disease, and strict protective measures are needed to keep the number of patients low and then eliminate the disease.We also fitted the model with the MCMC method based on the reported data and obtained 8 or 10 days of data to predict the disease, which can have a prediction result similar to the reported data.Our study also has limitations. Firstly, if more complete vaccination data are available, the effect of multiple doses of vaccinations or boosted vaccination on disease transmission can be considered. Secondly, when studying the impact of vaccination on disease transmission, it is also crucial to assess the effects of human behavioural changes and migration on disease. Moreover, as the complexity of the model increases, more complex dynamic behaviour may result.



---

*Source: 7596164-2022-09-16.xml*