Global stability and optimal vaccination control of SVIR models

1.   Introduction

Infectious diseases remain a significant contributor to overall human mortality. According to data reported by the WHO, a large number of people die from infectious diseases every year. In 2019, new measles cases worldwide increased to 869 770, and since 2016, the death rate from measles has increased by almost 50%. In 2019 alone, the number of measles deaths was about 207 500 ^[1]. In 2022, there were major outbreaks of monkeypox in the Americas and Western Pacific, with 87 929 confirmed cases and 146 deaths in 111 countries by June 2023 ^[2]. COVID-19 also remains a significant risk to people's lives, with at least 767 million cases confirmed and over 6.9 million fatalities worldwide by June 4, 2023 ^[3]. These data show that infectious diseases pose a severe threat to human lives.

Vaccination is one of the most important means of reducing both mortality and the number of infected persons. For most infectious diseases, vaccination could reduce the number of infected people, thus decreasing the mortality rate. New antibodies are produced when an unfamiliar etiologic agent or disease enters the body. Many vaccines are or contain a weakened or inactivated part of a specific pathogen that, when injected into the body, forms an antigen and induces an immune response ^[4]. For example, from 2000 to 2017, measles immunization has averted an estimated 21.1 million deaths, and vaccination efforts contributed to an 80% decline in measles deaths worldwide, making the measles vaccine one of the most effective products in public health ^[5]. After years of research, three vaccines (MVA-BN, LC16, and OrthopoxVac) have been approved for the prevention of monkeypox ^[6]. Coronavirus vaccination has been shown to increase immunity to COVID-19 and reduce the incidence of the disease. As of today, at least 70.3% of the global population has received at least a single dose of coronavirus vaccine ^[7].

Many scholars have explored epidemic dynamical models of infectious diseases. In 1927, Kermack and McKendrick first studied the deterministic compartment model and proposed the susceptible-infectious-recovered (SIR) model ^[8]. Other scholars have since developed improved epidemic dynamical models on this foundation, including the susceptible-infectious-susceptible (SIS) and susceptible-exposed-infectious-recovered (SEIR) models ^[9]. Considering the susceptible-vaccinated-infectious-recovered (SVIR) model with vaccination age, Duan et al. showed the influence of vaccine effectiveness and vaccination age on disease transmission ^[10]. By establishing a SVIR model with two different vaccination strategies, Liu et al. demonstrated that diseases can be eliminated under appropriate vaccination strategies ^[11].

Guzzi et al. investigated some graph-based epidemiological models, and their results showed that these models significantly improved disease transmission control ^[12]. Petrizzelli et al. considered a network-based model to assess the impact of different vaccination strategies and showed that topological awareness and age-based vaccination strategies were able to mitigate the spread of virus ^[13]. Han et al. formulated the vaccinated-susceptible-exposed-infectious-recovered (VSEIR) model with vaccination and immunity decline, and investigate the impact of vaccination on the transmission of COVID-19 ^[14]. Further, they considered the effects of the Omicron variant, and the susceptible-exposed-asymptomatic infected-infected-recovered (SEI_AIR) model was developed to explore prevention and control strategies to contain COVID-19 ^[15]. Song et al. applied an susceptible-exposed-asymptomatic infected-symptomatic infectious-isolated-recovered (SEAIQR) model with vaccination and isolation delays, and their results showed that decreasing the isolation delays and increasing the vaccination rate were effective in preventing the spread of COVID-19 ^[16]. Zhang et al. established a multi-patch model of heterosexual transmission and investigated the impact of population migration on HIV /AIDS spreading ^[17]. Considering different age groups, Chen et al. investigated the transmission dynamics model and optimal control strategy of COVID-19 in Hong Kong ^[18]. To study the monkeypox outbreak in Nigeria, Li et al. developed a new mathematical model to study the optimal control measures to contain the epidemic ^[19].

There is always some level of discrepancy between the actual data and the parameter values. However, sensitivity analysis can be used to measure the robustness of the predicted values with respect to the parameter values provided ^[20]. By analyzing the sensitivity of the basic reproduction number, Ndaïrou et al. found the two parameters with the largest impact on ${R_0}$ to be the contact rate and the mortality rate from disease ^[21]. To identify the vital parameters in the mechanism of the Zika virus, Biswas et al. performed a sensitivity analysis of important threshold numbers, and their results produced a scheme to control the transmission of Zika virus ^[22]. Optimal control can help decision-makers develop the most cost-effective solutions to minimize losses to people and society. Li et al. used optimal control theory to develop an optimal solution to this problem and to realize the aim of the lowest possible number of infections and the lowest possible vaccination rate for a given vaccination period ^[23]. Wang et al. studied the time-varying optimal control model, considering seasonal factors, and they addressed the issue numerically using the symplectic pseudospectral method, showing that seasonality changes may cause optimal vaccination strategies to also shift ^[24].

Our model takes into account the vaccination rate, vaccine effectiveness and the vaccine-induced decline in the immune rate, thus analyzing the stability of equilibria. We also performed a sensitivity analysis on these threshold parameters. We found solutions to the optimal vaccination control problem using the Pontryagin maximum principle, and we chose the proper parameters to prove that the theoretical results hold.

This paper is organized as follows. The formulation of the model and the basic reproduction number are presented in Section 2. The stability of equilibria is analyzed in Section 3. The bifurcation analysis of the model is presented in Section 4. The sensitivity analysis of the threshold parameters is described in Section 5. The optimal vaccination model is shown in Section 6. Numerical simulation is presented in Section 7. The conclusion is presented in Section 8.

2.   Model formulation

Developing an epidemic dynamical model of vaccination, the total population is $N(\text{t})$, which consists of vaccinated $V\left(t \right)$, susceptible $S\left(t \right)$, infected $I\left(t \right)$, and recovered individuals $R\left(t \right)$, where

Suppose all the population is homogeneously mixed. After vaccination, susceptible people will carry antibodies to the vaccine and develop immunity. With the increase of time, the immunity of the vaccine will weaken and eventually disappear, and the vaccinated people become susceptible. Susceptible people become infected by coming into contact with infected people, and because immunity disappears over time, a subset of individuals in the vaccinated population become infected again. As the treatment progresses, some of the infected individuals may recover, while others may succumb to the disease. Assuming that the natural mortality rate is equal in each compartment, the schematic diagram of this model is in Figure 1, and the epidemic dynamical model is given by the following four equations:

where $\Lambda$ represents the susceptible individuals recruited by immigration or birth. $\beta$ represents the effective contact rate which results from contact between susceptible individuals and infected individuals, $d$ represents the natural death rate of each compartment. $\theta$ represents vaccine effectiveness rate, which is between 0 and 1. When $\theta = 0$, the vaccines are invalid, and when $\theta = 1$, it shows that the vaccines are completely active. $\gamma$ represents the vaccination rate of susceptible individuals. $b$ represents the cure rate of infected individuals. $\alpha$ represents the disease-induced mortality, since the immunity of the vaccine weakens over time. $\delta$ represents the loss rate of immunity. All parameters are non-negative.

2.1.   Invariant set

Theorem 2.1. The solutions of the system remain non-negative, when $S\left(0 \right) > 0$, $I\left(0 \right) \geqslant 0$, $V\left(0 \right) \geqslant 0$, $R\left(0 \right) \geqslant 0$, for all $t > 0$.

Proof. According to the expression of the system (2.2), we have

Simultaneously integrating from $0$ to $t$ for the above equation, then getting the analytic expression for $I\left(t \right)$,

Similarly,

Since $S\left(0 \right) > 0$, $I\left(0 \right) \geqslant 0$, $V\left(0 \right) \geqslant 0$, $R\left(0 \right) \geqslant 0$, that is, $S\left(t \right) > 0$, $I\left(t \right) \geqslant 0$, $V\left(t \right) \geqslant 0$, $R\left(t \right) \geqslant 0$.

Theorem 2.2. The closed set

is a bounded set.

Proof. According to ^[25], we have

by the comparison theorem ^[26],

In summary, the close set $\Gamma$ is a positively invariant set.

2.2.   The basic reproduction number

The basic reproduction number ${R_0}$ is an essential threshold parameter which can show the prevalence of the disease in the epidemic dynamical model. If ${R_0} < 1$, it will die out, and if ${R_0} > 1$, it will become prevalent in the population. In our paper, ${R_0}$ is determined via the next-generation matrix method ^[27].

When $\theta = \gamma = \delta = 0$, the system (2.2) reduced to the SIR model without the vaccine,

computing the regeneration matrix and transfer matrix for the SIR model,

differentiating ${\mathcal{F}_0}$ and ${\mathcal{V}_0}$, calculate the value of ${E_{00}} = \left({{S_{00}}, 0, 0} \right) = \left({\frac{\Lambda }{d}, 0, 0} \right)$ which is the disease-free equilibrium point of SIR model

Then, ${R_0} = \rho \left({{F_0}{V_0}^{ - 1}} \right) = \frac{{\beta {S_{00}}}}{{d + b + \alpha }} = \frac{{\beta \Lambda }}{{d\left({d + b + \alpha } \right)}}$, and $\frac{1}{{d + b + \alpha }}$ is the number of the secondary infections developed by an infected person throughout his or her infectious period ^[28].

To compute the system's disease-free equilibrium, let the left of system (2.2) equal to be zero, that is,

Let $x = {\left({I, R, V, S} \right)^T}$, and the system (2.2) can be written as $\dot x = \mathcal{F} - \mathcal{V}$, where the regeneration matrix $\mathcal{F}$ and the transfer matrix $\mathcal{V}$ is as follows:

differentiating $\mathcal{F}$ and $\mathcal{V}$, calculate the value of ${E_0}$,

where

since the basic reproduction number is the spectral radius of $F{V^{ - 1}}$, thus

So, the basic reproduction number of the system (2.2) is

The interpretation of Eq (2.6) is as follows: $\frac{\beta }{{\left({d + \alpha + b} \right)}}\frac{{\theta \Lambda \gamma }}{{d(\delta + d + \gamma)}}$ is the proportion of secondary infections attributed to the average vaccinated people with a vaccination rate $\gamma$, $\frac{\beta }{{\left({d + \alpha + b} \right)}}\frac{\Lambda }{d}$ is the number of secondary infections attributed to the susceptible individuals ^[29].

3.   Stability analysis

In the above, we defined the basic reproduction number of the system (2.2). In the following text, the local asymptotic stability (LAS) and globally asymptotically stable (GAS) of the equilibria is analyzed using the relationship between the magnitude of ${R_{vac}}$ and unity.

3.1.   Local stability of disease-free equilibrium

Since the equation for variable $R\left(t \right)$ is decoupled from the first three equations of system (2.2), it is just considering the following system (3.1).

Theorem 3.1. When ${R_{vac}} < 1$, ${E'_0}$ is LAS, and when ${R_{vac}} > 1$, it is unstable.

Proof. The disease-free equilibrium of the system (3.1) is as follows:

Calculating the Jacobian matrix at ${E'_0}$,

where $\beta {S_0} + \left({1 - \theta } \right)\beta {V_0} - (d + \alpha + b) = \left({{R_{vac}} - 1} \right)(d + \alpha + b)$.

The characteristic equation is

then solve Eq (3.3), and the eigenvalues are as follows:

that is, when ${R_{vac}} < 1$, all eigenvalues are negative, which shows that ${E'_0}$ is LAS, and otherwise, it is unstable.

3.2.   Global stability of disease-free equilibrium

In this context, we analyze the global stability of ${E'_0}$. ^[30] proposed a method to prove the global asymptotically stable of the disease-free equilibrium, by dividing system (3.1) into two subsystems, that is,

where $X$ denotes uninfected individuals, and $Z$ denotes infected individuals. If system (3.1) satisfies (H1) and (H2), the disease free equilibrium is globally asymptotically stable.

(H1) If $\frac{{{\text{d}}X}}{{{\text{d}}t}} = F\left({X, 0} \right)$, ${X^*}$ is the disease free equilibrium;

(H2) $G\left({X, Z} \right) = FZ - \hat G\left({X, Z} \right)$, where $\hat G\left({X, Z} \right) \geqslant 0$ for $G\left({X, Z} \right) \in \Omega$, and $F = {D_Z}G\left({{X_0}, 0} \right)$ is a M-matrix (i.e., the M-matrix is the non-diagonal element is non-negative), $\Omega$ is positive invariant set. In this paper, it can obtain the following Theorem 3.2.

Theorem 3.2. If ${R_{vac}} < 1$, ${E'_0}$ is GAS.

Proof. Let $X$ denotes the uninfected individuals, $Z$ denotes the infected individuals. So, we can get $X{\text{ = }}{\left({V, S} \right)^T}$, and $Z = I$.

From the Eq (3.2), we have computed ${E'_0}$. Here, it is noted as ${U_0} = \left({{X_0}, 0} \right)$, where

For $\frac{{{\text{d}}X}}{{{\text{d}}t}} = F\left({X, 0} \right)$, assumption (H1) is satisfied, and when $t \to + \infty$, $X \to {X_0}$, i.e., ${X_0}$ is GAS.

Computing matrix $F$ at ${X_0}$,

where $F$ is a M-matrix.

from assumption (H2),

thus

Since $\Gamma$ is a positive invariant set, i.e., $S \leqslant \frac{{\Lambda (\delta + d)}}{{d(\delta + d + \gamma)}}$, $V \leqslant \frac{{\Lambda \gamma }}{{d(\delta + d + \gamma)}}$, then $\hat G\left({X, Y} \right) \geqslant 0$. That is, when $t \to + \infty$, then $I\left(t \right) \to {I_0}$.

Consequently, if ${R_{vac}} < 1$, ${E'_0}$ is GAS.

3.3.   Existence of endemic equilibrium

From the system (3.1), the expressions for calculating the endemic equilibrium are as follows. ${I^*} \ne 0$, so by the third expression of the system (3.1), where ${S^*} = \frac{{\left({d + b + \alpha } \right) - \left({1 - \theta } \right)\beta {V^*}}}{\beta }$, bringing ${S^*}$ into the first equation

so

bringing ${S^*}$ and ${V^*}$ into the second equation, set ${I^*}$ as the positive root of Eq (3.4).

where

obviously, $A < 0$. Clearly, it shows that ${R_{vac}} > 1$$\Leftrightarrow$$C > 0$; ${R_{vac}} = 1$$\Leftrightarrow$$C = 0$; ${R_{vac}} < 1$$\Leftrightarrow$$C < 0$.

Since $A < 0$ combined with the above equivalence relationship, the results show that: when ${R_{vac}} = 1$, i.e., $C = 0$, Eq (3.4) has one positive root provided $B > 0$, otherwise Eq (3.4) has no positive root; when ${R_{vac}} > 1$ i.e., $C > 0$, Eq (3.4) has a unique positive root ${I_1} = \frac{{ - B - \sqrt {{B^2} - 4AC} }}{{2A}}$. Next, we will consider the case of the roots of the equation, when ${R_{vac}} < 1$, i.e., $C < 0$.

Let $\Delta = \sqrt {{B^2} - 4AC} = 0$, that is,

Set

And then get the following equivalence:

Thus, when it satisfies ${R_{vac}} < 1$ and $B \leqslant 0$, Eq (3.4) does not exist positive root, and when it satisfies ${R_{vac}} < 1$ and $B > 0$, the following equivalence relationship is established:

(1) if ${R^*} < {R_{vac}} < 1$, there are two positive roots ${I_2} = \frac{{ - B - \sqrt {{B^2} - 4AC} }}{{2A}}$ and ${I_3} = \frac{{ - B + \sqrt {{B^2} - 4AC} }}{{2A}}$;

(2) if ${R_{vac}} = {R^*} < 1$, there is a unique positive root ${I_4} = - \frac{B}{{2A}}$;

(3) if ${R_{vac}} < {R^*}$, there is no positive root.

In summary, the following theorem can be derived.

Theorem 3.3. In the system (3.1), when ${R_{vac}} > 1$, there exists a unique endemic equilibrium ${E_1} = \left({{I_1}, {S_1}, {V_1}} \right)$; when ${R_{vac}} = 1$, there is an endemic equilibrium with $B > 0$; when ${R_{vac}} < 1$ and $B \leqslant 0$, there exists no endemic equilibrium; when ${R_{vac}} < 1$ and $B > 0$, there are two endemic equilibria ${E_2} = \left({{I_2}, {S_2}, {V_2}} \right)$ and ${E_3} = \left({{I_3}, {S_3}, {V_3}} \right)$ for ${R^*} < {R_{vac}}$, the unique endemic equilibrium ${E_4} = \left({{I_4}, {S_4}, {V_4}} \right)$ for ${R^*} = {R_{vac}}$, no endemic equilibrium for ${R_{vac}} < {R^*}$.

3.4.   Local stability of endemic equilibrium

Lemma 1. ^[31] Suppose ${A^*}$ is $3 \times 3$, a real matrix, if $tr\left({{A^*}} \right) < 0$, $\det \left({{A^*}} \right) < 0$, $\det \left({{A^*}^{[2]}} \right) < 0$. Then, all eigenvalues of ${A^*}$ have negative real parts.

In ^[32], for $n \times n$ matrix ${A^*} = \left({{a_{ij}}} \right)$, when $n = 2$, the second additive compound matrix of ${A^*}$ is ${A^*}^{[2]} = {a_{11}} + {a_{22}}$, when $n = 3$, ${A^*}^{[2]} = \left[{\begin{array}{*{20}{c}} {{a_{11}} + {a_{22}}} & {{a_{23}}} & { - {a_{13}}} \\ {{a_{32}}} & {{a_{11}} + {a_{33}}} & {{a_{12}}} \\ { - {a_{31}}} & {{a_{21}}} & {{a_{22}} + {a_{33}}} \end{array}} \right]$.

Theorem 3.4. When ${R_{vac}} > 1$ and $\left\{ {\delta - \left({1 - \theta } \right)\left[{\left({1 - \theta } \right)\beta {I_1} + \left({\delta + d} \right)} \right]} \right\} < 0$, ${E_1}$ is LAS.

Proof. Calculating the Jacobian matrix of the system (3.1),

we can compute $\beta S + \left({1 - \theta } \right)\beta V - \left({d + \alpha + b} \right){\text{ = 0}}$, then

Computing ${J^{\left[2 \right]}}$

then, calculate the value of the determinant of ${J^{\left[2 \right]}}$

When ${R_{vac}} > 1$, bringing ${E_1}$ into matrices (3.8) and (3.10), it can be calculated as $tr\left({{J_1}} \right) < 0$, $\det \left({{J_1}} \right) < 0$, if $\left\{ {\delta - \left({1 - \theta } \right)\left[{\left({1 - \theta } \right)\beta {I_1} + \left({\delta + d} \right)} \right]} \right\} < 0$ is satisfied, then $\det \left({J_1^{[2]}} \right) < 0$. According to Lemma 1, every eigenvalue of ${J_1}$ has a negative real part, it shows that ${E_1}$ is LAS.

3.5.   Global stability of endemic equilibrium

In this context, the geometric method is used to demonstrate GAS of endemic equilibriums. In Li and Muldowney ^[33], it shows the theoretical knowledge about how to prove GAS of endemic equilibriums, which considering the following differential system,

where $f\left(x \right):D \to {R^n}$ be ${C^1}$ continuous, $D \subset {R^n}$ be an open set, $x\left({t, {x_0}} \right)$ is uniquely determined by the initial value condition $x\left(0 \right) = {x_0}$. Suppose system $\dot x = f\left(x \right)$ satisfies the following three hypotheses:

(1) The open set $D$ is a simply connected,

(2) There exists a compact absorbing set $K \subset D$,

(3) There exists a unique endemic equilibrium $\bar x$ in open set $D$.

Suppose $|·|$ represents the vector norm induced by the $n \times n$ order matrix. The Lozinskii measure of matrix $W$ with respect to the induced matrix norm is defined by

Lemma 2. ^[34] The Lozinskii measure of a real $n \times n$ matrix ${M^*}$ for the matrix norm induced by the ${l_1}$ vector norm is given by

Theorem 3.5. If $\bar q \leqslant 0$, $\bar x$ is GAS.

In the system (2.2), we have the following lemma:

Lemma 3. ^[35,36] The system (2.2) is uniformly persistent, i.e., there exists constant $\varphi > 0$ such that

Proof. In the above section, ${E_0}$ exists on the boundary of region $\Gamma$, and Theorem 3.1 provides the stability condition for ${E_0}$, i.e., ${E_0}$ is stable when ${R_{vac}} < 1$, and otherwise ${E_0}$, is unstable, which leads to the instability of ${E_0}$. Since the instability of ${E_0}$ is equivalent to the uniform persistence of the system (2.2), it can be proved that the system (2.2) uniformly persists in the interior of $\Gamma$, that is, there exists a constant $\varphi > 0$ such that

with original values $\left({V\left(0 \right), S\left(0 \right), I\left(0 \right), R\left(0 \right)} \right) \in D$, which shows the hypothesis (2) holds.

In the system (2.2), when ${R_{vac}} > 1$, there exists ${E_1}$ which is a unique equilibrium in the region of $\Gamma$, and then hypotheses (1) and (3) are satisfied.

Theorem 3.6. When ${R_{vac}} > 1$, ${E_1}$ is GAS.

Proof. Consider the following subsystem (3.13):

the Jacobian matrix is as follows:

then, the second compound matrix is

Set matrix-value functions

thus ${P_f} = \left[{\begin{array}{*{20}{c}} {\frac{{V'I - I'V}}{{{I^2}}}} & 0 \\ 0 & {\frac{{V'I - I'V}}{{{I^2}}}} \end{array}} \right]$, ${P^{ - 1}} = \left[{\begin{array}{*{20}{c}} {\frac{I}{V}} & 0 \\ 0 & {\frac{I}{V}} \end{array}} \right]$, ${P_f}{P^{ - 1}} = \left[{\begin{array}{*{20}{c}} {\frac{{V'}}{V} - \frac{{I'}}{I}} & 0 \\ 0 & {\frac{{V'}}{V} - \frac{{I'}}{I}} \end{array}} \right]$,

then, $P{J_1}^{\left[2 \right]}{P^{ - 1}} = \left[{\begin{array}{*{20}{c}} Q & 0 \\ 0 & Q \end{array}} \right]$.

Clearly, we can obtain

where

The Lozinskii measure is determined as follows:

where ${g_1} = {\mu _1}\left({{W_{11}}} \right) + \left| {{W_{12}}} \right|$, ${g_2} = \left| {{W_{21}}} \right| + {\mu _1}\left({{W_{22}}} \right)$, obviously

Thus, we have

According to the system (2.2), computing that,

substitute the above equation into ${\mu _1}\left(W \right)$, that is

then the measure $\bar q$ is as follows

The above inequality indicates $\bar q \leqslant 0$, so $I \to {I_1}$, $V \to {V_1}$ with $t \to + \infty$.

Under expression of the system (2.2), we have

when $t \to + \infty$, we have

that is $\frac{{{\text{d}}S(t)}}{{{\text{d}}t}} + \left({\beta {I_1} + (\gamma + d)} \right)S = \Lambda + \delta {V_1}$, apparently a linear differential equation about $S$, thus $S \to \frac{{\Lambda + \delta {V_1}}}{{\beta {I_1} + \left({\gamma + d} \right)}} = {S_1}$, $R \to \frac{{b{I_1}}}{d} = {R_1}$ as $t \to + \infty$ ^[37].

Comprehensively stated, when $t \to + \infty$, we have $\left({V, S, I, R} \right) \to \left({{V_1}, {S_1}, {I_1}, {R_1}} \right)$, thus ${E_1}$ is GAS.

4.   Bifurcation analysis

Bifurcation analysis plays a vital role in controlling and eradicating epidemic diseases. In modeling epidemic diseases, when the basic reproduction number is less than unity, the disease will be eradicated. However, when bifurcation occurs, and the basic reproduction number exceeds unity, the endemic equilibrium will persist, i.e., susceptible individuals and infected individuals coexist. Thus, it is significant to analyze the bifurcation of the model in epidemiology.

4.1.   Forward bifurcation

Theorem 4.1. If ${R_{vac}} = 1$, the system (2.2) will undergo forward bifurcation at $\beta = {\beta ^*}$.

Proof. Compute the Jacobian matrix at ${E'_0}$,

when ${R_{vac}} = 1$, calculate the characteristic polynomial as follows:

that is

The three eigenvalues that can be calculated for this characteristic equation are

There is a zero eigenvalue and the rest of the roots are negative, in the light of the center manifold theorem ^[38], the system (3.1) exists the forward bifurcation.

When ${R_{vac}} = 1$, let bifurcation parameter

Then compute the eigenvector whose eigenvalue is zero, set the right eigenvector $w = {\left({{w_1}, {w_2}, {w_3}} \right)^T}$, and then we can obtain

getting the following relation:

Similarly, set the left eigenvector $v = \left({{v_1}, {v_2}, {v_3}} \right)$, and then we can obtain

it also satisfies $v \cdot w = 1$, that is, ${v_1}{w_1} + {v_2}{w_2} + {v_3}{w_3} = 1$. So we can get the left eigenvector as follows,

To determine the type of bifurcation, set $V = {x_1}$, $S = {x_2}$, $I = {x_3}$, $\frac{{{\text{d}}V}}{{{\text{d}}t}} = {f_1}\left(x \right)$, $\frac{{{\text{d}}S}}{{{\text{d}}t}} = {f_2}\left(x \right)$, $\frac{{{\text{d}}I}}{{{\text{d}}t}} = {f_3}\left(x \right)$, according to the system (3.1), it can obtain

According to $v$ and $w$ computing ${e_1}$ and ${e_2}$, where

Then, calculate the partial derivative at $\left({{{E'}_0}, {\beta ^*}} \right)$,

and others are zero. Then ${e_1}$ and ${e_2}$ are computed,

According to ^[39], the system (3.1) will undergo forward bifurcation. Since the fourth equation of system (2.2) is decoupled from the previous three, system (3.1) has the same stability as system (2.2), and there also exists a forward bifurcation for the system (2.2).

4.2.   Hopf bifurcation

Theorem 4.2. ^[40] Suppose ${R_{vac}} > 1$. When vaccination rate $\gamma$ exceeds thresholds ${\gamma ^*}$, there will be a Hopf bifurcation around ${E_1}$.

Proof. When ${R_{vac}} > 1$, calculate the Jacobian matrix of ${E_1}$, and then matrix (3.8) is as follows:

Calculate its characteristic equation,

where

Let ${\gamma ^*} = \gamma$, and ${a_1}\left({{\gamma ^*}} \right){a_2}\left({{\gamma ^*}} \right) - {a_3}\left({{\gamma ^*}} \right) = 0$, and the roots of characteristic Eq (4.6) are

${\lambda _1}$ and ${\lambda _2}$ are a pair of purely imaginary roots.

In general, the form of characteristic roots $\lambda$ is

thus, when ${\lambda _1} = {\kappa _1} + {\kappa _2}i$, characteristic Eq (4.6) is replaced by the following equation:

separate the real and imaginary parts, let $R\left({{\kappa _1}, \gamma } \right) = {\kappa _1}^3 - 3{\kappa _1}{\kappa _2}^2 + {a_1}{\kappa _1}^2 - {a_1}{\kappa _2}^2 + {\kappa _1}{a_2} + {a_3}$, then, compute the transversality conditions, by the implicit function theorem,

and ${\kappa _1} = 0$, ${\kappa _2} = \sqrt {{a_2}}$, then

since $\text{sgn} (\frac{{{\text{d}}{\kappa _1}}}{{{\text{d}}\gamma }}) = - 1$, and thus it generates Hopf bifurcation at $\gamma = {\gamma ^*}$.

5.   Sensitivity analysis

In determining how vaccines contribute to the reduction of mortality from disease in populations and control disease transmission of disease, we must explore the importance of vaccines in relation to different factors involved in disease transmission and epidemiology. Sensitivity analysis of ${R_{vac}}$ shows how each parameter contributes to the propagation of the disease. In Section 2.2, we calculate ${R_{vac}}$, which is correlated with the transmission and prevalence. In this paper, the normalized forward sensitivity index analysis is used to analyze the parametric variable's effect on ${R_{vac}}$.

Definition 1. ^[41] The normalized forward sensitivity index of ${R_{vac}}$ that depends differentiably on a parameter $\omega$, which is defined as

shows that ${R_{vac}}$ rises with $\omega$ increasing when ${O_\omega } > 0$, but ${R_{vac}}$ declines with increasing $\omega$ when ${O_\omega } < 0$.

The normalized forward sensitivity index of ${R_{vac}}$ with regard to $\beta, \Lambda, b, \alpha, \gamma, \theta$ as follows:

these equations show that $\beta$, $\Lambda$, $\delta$ has a negative effect on ${R_{vac}}$, and when $\beta$ and $\Lambda$ increase, ${R_{vac}}$ also increases. Conversely $b, d, \alpha, \gamma, \theta$ have a positive effect on ${R_{vac}}$, which shows that the number of individuals with secondary infections is lower, when $b, d, \alpha, \gamma, \theta$ increase, which shows that vaccines are very successful in halting the propagation of diseases.

6.   Optimize control

Optimal control theory is employed to identify method of to achieve the maximum performance at the minimum cost under different assumptions. In Section 5, we analyzed the sensitivity of ${R_{vac}}$. Therefore, in conjunction with the actual situation of infectious diseases transmission, we formulate the optimal vaccination control strategy ^[23].

In this section, to get as many people vaccinated as possible, the vaccination rate $\gamma$ is selected as the control variable $u\left(t \right)$. Since the vaccines produce are limited, thus the vaccinated people also are limited, it is possible to get $uS \leqslant \Omega$, where $\Omega$ is the maximum vaccination number.

Thus, the optimal vaccination control problem with inequality constraints is described as follows,

where ${t_f}$ represents the final time and ${t_f} \in {R^ + }$, ${P_1}$ represents the weight coefficient, and ${P_1} \in {R^ + }$.

In Problem $\left(P\right)$, there exists the inequality constraint ($0 \leqslant u \leqslant {u_{\max }}, 0 \leqslant S \leqslant {S_{\max }}, uS \leqslant \Omega$), and it is translated into equality constraints by non-negative parameters ${\eta _i}\left({i = 1, 2, 3, 4} \right)$, that is,

Thus, the Hamiltonian function of $Problem \left(P\right)$ is expressed as follows ^[42,43]:

where $\lambda = {\left[{{\lambda _V}, {\lambda _S}, {\lambda _I}, {\lambda _R}} \right]^T}$ is the cost variable, and $\mu = {\left[{{\mu _1}, {\mu _2}, {\mu _3}, {\mu _4}} \right]^T}$ is non-negative.

In order to prove the existence of the optimal vaccination strategy ${u_{^*}}$, we apply Theorem 2.2 in ^[44] and Theorem 5.1 in ^[45]. We need to check the following assumptions:

(1) The set of controls and corresponding state variables is nonempty.

(2) The admissible set is convex and closed.

(3) The right hand side of (6.1) is bounded by a linear function in the state and control variables.

(4) The integrand of the objective functional is convex.

(5) There exist constants ${C_1} > 0$, ${C_2} > 0$, and $\tau > 1$ such that the integrand of the objective functional satisfies $Y \geqslant {C_1}{\left| u \right|^\tau } - {C_2}$.

It is obvious that assumption 1 holds. Note that the solutions are bounded, so assumption 2 holds. To simplify the notation, we define $L\left({t, \phi, u} \right)$:

For assumption 3, there exist ${\phi _1}(t)$, ${\phi _2}(t)$, it follows that

where ${M_1} = 2u + d + 2\beta {I_2}$, ${M_2} = 2\left({1 - \theta } \right)\beta {V_1} + \beta {S_1} + \left({2 - \theta } \right)\beta {I_2} + 2b + d + \alpha$, ${M_4} = d$, ${M_3} = 2\left({1 - \theta } \right)\beta {I_2} + 2\delta + d$, $M = \max \left\{ {{M_1}, {M_2}, {M_3}, {M_4}} \right\}$. The assumption 3 holds.

For assumption 4, there exists the constant $k \in \left({0, 1} \right)$ and ${u_1}$, ${u_2}$, then it obtains

Hence, $Y\left({\left({1 - k} \right){u_1} + k{u_2}} \right) < \left({1 - k} \right)Y\left({{u_1}} \right) + kY\left({{u_2}} \right)$, it proves that assumption 4 holds.

For assumption 5, it is obviously that

and ${C_1} = 1$, $\tau = 2$, ${C_2} > 0$. The assumption 5 is satisfied. The existence of ${u_{^*}}$ is proved. In the next context, denote the optimal control of $Problem \left(P\right)$ as ${u_*}$ and the trajectory of the corresponding optimal control state as $x = \left[{{V^*}, {S^*}, {I^*}, {R^*}} \right]$. In ^[44], it shows that base on a set of necessary conditions that the optimal control and state must satisfy, which is the optimality conditions, the adjoint equations and the transversality conditions, solving the optimal control problems. Here is the specific solution step to find the necessary conditions.

The first-order necessary conditions can derive by utilizing the Pontryagin maximum principle ^[44]. Based on the conditions of optimality, co-state, and parametric variables, the two-point boundary value problems and nonlinear complementarity problems can formulated as follows ^[24]:

where $\eta \geqslant 0$, $\mu \geqslant 0$, ${\eta ^T}\mu = 0$.

The optimal conditions for the control variables $u\left(t \right)$ are given by:

Combining Eqs (6.3) and (6.4), bringing the optimal solution, and ${\dot{\lambda }}_{V}, {\dot{\lambda }}_{S}, {\dot{\lambda }}_{I}, {\dot{\lambda }}_{R}$ are as follows:

where the transversality conditions are determined as

For a particular optimal control problem, Eqs (6.4), (6.5) and (6.7) form a set of necessary conditions that an optimal control and state must satisfy.

By the Section 3.1 in ^[44], if our goals depend on the state at the terminal time, the objective function is

where $\phi \left({x\left({{t_1}} \right)} \right)$ represents the payoff term, which is a goal with respect to the final position or population level, $x\left({{t_1}} \right)$. Then, compute the necessary conditions, where the transversality conditions have changed. The transversality conditions are as

In the optimal vaccination control problem of this paper, the objective functional did not explicitly depend on the state at the terminal time. Thus, let the transversality conditions be stated as: ${\lambda _V}\left({{t_f}} \right) = {\lambda _S}\left({{t_f}} \right) = {\lambda _I}\left({{t_f}} \right) = {\lambda _R}\left({{t_f}} \right) = 0.$

Let the control variable $u$ meet $0 \leqslant u \leqslant {u_{\max }}$ and $u \leqslant \frac{\Omega }{S}$, and then $0 \leqslant u \leqslant \min \left\{ {{u_{\max }}, \frac{\Omega }{S}} \right\}$ is received. Set ${u_{\max }}$ is equal to unity and $\Omega < {S_{\max }}$, and thus, it can be obtained that $\frac{\Omega }{S} < {u_{\max }}$.

By solving Eq (6.5), the optimal vaccination control ${u_*}$ can be solved

In the following discussion, consider $\frac{\Omega }{S} > {u_{\max }}$, under the case with ${u_*} < \frac{\Omega }{{{S^*}}}$, where ${\mu _4} \equiv 0$. In order to determine the explicit expression without ${\mu _1}$, ${\mu _2}$, consider the following three cases:

(1) when $0 < {u_*} < {u_{\max }}$, ${\mu _1} = {\mu _2} = 0$, hence ${u_*} = \frac{{\left({{\lambda _S} - {\lambda _V}} \right){S^*}}}{2}$;

(2) when ${u_*} = 0$, ${\mu _2} = 0$, hence $0 = {u_*} = \frac{{\left({{\lambda _S} - {\lambda _V}} \right){S^*} + {\mu _1}}}{2}$, and when ${\mu _1} \geqslant 0$, there is $\left({{\lambda _S} - {\lambda _V}} \right){S^*} \leqslant 0$;

(3) when ${u_*} = {u_{\max }}$, ${\mu _1} = 0$, hence ${u_{\max }} = {u_*} = \frac{{\left({\lambda {}_S - {\lambda _V}} \right){S^*} - {\mu _2}}}{2}$, and when ${\mu _2} \geqslant 0$, there is $\left({{\lambda _S} - {\lambda _V}} \right){S^*} \geqslant {u_{\max }}$.

Summarize the three conclusions when $\frac{\Omega }{S} > {u_{\max }}$, the optimal vaccination control solution ${u_*}$ is as follows,

When $\frac{\Omega }{S} \leqslant {u_{\max }}$, the optimal vaccination control solution ${u_*}$ is as follows

Thus, the analytical formula for optimal control is:

7.   Numerical simulation

In this section, numerical simulations will be performed to demonstrate the results of the above theory using MATLAB.

For the parameter values: $b = 0.008$, $d = 0.0518$, $\theta = 0.8125$, $\beta = 0.0563$, $\delta = 0.111$, $\Lambda = 1$, $\alpha = 0.35$, $\gamma = 0.99$, and ${R_{vac}} = {\text{0}}{\text{.8016 < 1}}$ which satisfied the condition of the Theorem 3.2. It can be seen from Figure 2, which proved the validity of the Theorem 3.2.

The parameter values are as follows:$b = 0.0085$, $d = 0.005$, $\theta = 0.6$, $\beta = 0.29$, $\delta = 0.005$, $\Lambda = 10$, $\alpha = 0.0099$, $\gamma = 0.2$, we can obtain ${R_{vac}}{\text{ = 2}}{\text{.4716}} \times {\text{1}}{{\text{0}}^4}$, in Figure 3, it can observe that for any given different initial values, these initial points can eventually tend to be endemic equilibrium, that is, which proves the Theorem 3.6 numerically.

The bifurcation diagram is shown in Figure 4, with bifurcation parameter $\beta \in \left(0.15, 0.85\right)$, and other parameter values are $\alpha {\text{ = }}0.11$, $\gamma {\text{ = }}0.85$, $\delta {\text{ = }}0.005$, $\theta {\text{ = }}0.95$, $b = 0.943$, $d = 0.0718$, $\Lambda = 2$. In Figure 4, it can be seen that the system generates a forward bifurcation at ${R_{vac}} = 1$. When ${R_{vac}} < 1$, ${E_0}$ is stability, when ${R_{vac}} > 1$, there exists a stable ${E_1}$ and an unstable ${E_0}$.

In Section 4, we have shown that the system generates Hopf bifurcation when the parameter $\gamma$ crosses the threshold ${\gamma ^*}$. By choosing the appropriate parameters, the result of the Theorem 4.1 is proved numerically. In Figure 5(a), (b) we show the stability of the system (2). In Figure 5(a) $b = 0.85$, $d = 0.001$, $\theta = 0.3$, $\beta = 0.0006$, $\delta = 0.009$, $\Lambda = 100$, $\alpha = 0.95$, $\gamma = 0.065$, and ${R_0} = 23.4129$, which shows there exists the unique endemic equilibrium. In Figure 5(b) $\gamma = 0.45$ and all other parameters are the same.

With given parameter values, the normalized forward sensitivity index is obtained as shown in Table 1.

In Table 1, we choose the proper parameters, obtain the value of the sensitivity index. The index of $\beta$, $\Lambda$, $\delta$ is positive, where when $\beta$, $\Lambda$ increase 10%, ${R_{vac}}$ will increase 10%. Similarly, when $\delta$ increases 10%, ${R_{vac}}$ will increase 0.338%. Inversely, the index of $b, d, \alpha, \gamma, \theta$ is negative, which means that $b, d, \alpha, \gamma, \theta$ are all negative effect on the basic regeneration number, that is, ${R_{vac}}$ decreases with the increase of $b, d, \alpha, \gamma, \theta$. For example, if $\gamma$ increases 10%, ${R_{vac}}$ will decrease 5.198%. These results show that when changing the vaccination rate, vaccine effective rate values and the rate of vaccine-induced immunity, ${R_{vac}}$ also changes dramatically. Thus, to control the transmission of disease, people are encouraged to be active in vaccinations, and biologists should increase the vaccine effectiveness rate.

In Section 6, the vaccination rate $\gamma$ was selected as the control parameter, and the expression for the optimal solution was obtained. Set the parameter values as $A = 0.002$, ${u_{\max }} = 1$, ${R_{\max }} = 1100$, and the final time is ${t_f} = 100$. From Figure 6, it can be identified that the infected is apparently lower after adding control measures (red curve) than when no control measures were added, which indicates that the addition of control measures is beneficial to restrain the transmission of the disease, and the prevention and control strategy is effective.

8.   Conclusions

In this paper, an epidemic dynamical model of vaccination is developed. It includes vaccine effectiveness, vaccination rate, and the vaccine-induced immune decline rate. The qualitative theoretical analysis of the model leads us to the following conclusions.

Through the next-generation method, we calculated ${R_{vac}}$ and confirm the stability of ${E_0}$. When ${R_{vac}} < 1$, ${E_0}$ is GAS. A high vaccine effectiveness and vaccination rate and a low vaccine-induced immune decline were found to drive the elimination of the disease. When ${R_{vac}} > 1$, ${E_1}$ is GAS, that is, when certain threshold conditions for ${R_{vac}}$ are achieved, the propagation of the disease can be contained by employing vaccination, increased vaccine effectiveness, and similar factors.

Using the center manifold theory, we showed that the system will generate the forward bifurcation at ${R_{vac}} = 1$. As the bifurcation parameter ${\beta ^*}$ increases, the system will undergo the forward bifurcation, and with it, reflecting that the disease will not become extinct. Therefore, when controlling the transmission of the disease, it is possible to reduce the patient's contact with the external community by taking preventive measures such as wearing masks, lengthening social distance, and isolation. We chose the vaccination rate $\gamma$ as the bifurcation parameter, providing theoretical demonstrations of the existence of the Hopf bifurcation under certain parameter conditions.

Next, sensitivity analysis of ${R_{vac}}$ shows that the vaccination rates and vaccine effectiveness are fundamental to ${R_{vac}}$.Vaccine-induced immune decline rate also plays an essential role in disease containment. Thus, we can show that vaccinations are a vital component in containing the transmission of disease. Community policymakers use the outcomes of sensitivity analyses to develop strategies to limit the transmission of diseases, which encourages people to be vaccinated while investing more in vaccine research to improve vaccine effectiveness rate.

To find the most cost-effective strategy for controlling infectious diseases, we selected the vaccination rate $\gamma$ as the control parameter, and the optimal control function with respect to vaccination is obtained, which was solved using the Pontryagin maximum principle. Results showed that vaccination is effective in decreasing the number of infected individuals, and vaccination is an essential way to slow down the number of infected people.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

This work was supported by the Gansu Science and Technology Fund (20JR5RA512), the Research Fund for Humanities and Social Sciences of the Ministry of Education (20XJAZH006), the Fundamental Research Funds for the Central Universities (31920220041), the Leading Talents Project of State Ethnic Affairs Commission of China, the Gansu Provincial Education Department's Outstanding Graduate Student "Innovation Star" Project (2023CXZX-197).

Conflict of interest

The authors declare that they have no conflict of interest.