An Analysis of the Impact of Stay-at-Home Measures on the Occurrence of Vaccine Shortages

  Nurunnahar

Department of Mathematics, Feni University, Feni, Bangladesh.

DOI: 10.4236/aid.2024.142030                

Abstract

COVID-19, a contagious respiratory disease, presents immediate and unforeseen challenges to people worldwide. Moreover, its transmission rapidly extends globally due to its viral transmissibility, emergence of novel strains (variants), absence of immunity, and human unawareness. This framework introduces a revised epidemic model, drawing upon mathematical principles. This model incorporates a modified vaccination and lockdown approach to comprehensively depict an epidemic’s transmission, containment, and decision-making processes within a community. This study aims to provide policymakers with precise information on real-world situations to assist them in making informed decisions about the implementation of lockdown strategies, maintenance variables, and vaccine availability. The suggested model has conducted stability analysis, strength number analysis, and first and second-order derivative analysis of the Lyapunov function and has established the existence and uniqueness of solutions of the proposed models. We examine the combined effects of an effective vaccination campaign and non-pharmaceutical measures such as lockdowns and states of emergency. We rely on the results of this research to assist policymakers in various countries in eradicating the illness by developing more innovative measures to control the outbreak.

Keywords

Epidemic Dynamics, Non-Pharmaceuticals Interventions, Vaccine Shortage, Multi-Waving

1. Introduction

The initial documentation of COVID-19, a communicable respiratory illness, originated in Wuhan, Hubei Province, Republic of China, on 31 December 2019 [1] . The lockdown/shutdown policy has been widely adopted globally to restrict the rapid spread and transmission of the virus due to its effective containment measures and the absence of adequate treatment options. In addition, the scarcity and postponement of vaccines are significant factors that justify implementing the lockdown policy, initially introduced in Wuhan in January 2020 [2] . The lockdowns, implemented with different degrees of severity, help policymakers reduce the loss of human lives if the effectiveness of the lockout measures is as anticipated. Furthermore, there is a substantial global scarcity of vaccines, except for a handful of affluent nations, and the politicians of these countries are apprehensive about mitigating the vaccine deficit for their citizens. Furthermore, the vaccine effectiveness rate plays a crucial role in discussions on vaccination. Consequently, the disease proliferates globally and poses a serious menace to humanity. To incorporate scenarios involving vaccine shortages and lockdown policies, we have utilized a modified SEIR (susceptible-exposed-infected-recovered) model [3] to represent epidemic dynamics accurately.

The inception of mathematical modeling of infectious illnesses may be dated back to the early 1900s [4] . Given the increasing risk of population illnesses in the human population, disease modeling exhibition has emerged as a crucial component of epidemic management [3] [5] - [10] . The mathematical model will enable us to calculate the possible effects of facilitation and the dynamics of infectious illnesses, allowing us to predict the development of an epidemic or pandemic. Various lockdown, quarantine, awareness, and vaccine models have been created to clarify COVID-19 situations. To capture how people’s behavior changes during economic shutdowns and the concept of immunity in the COVID-19 pandemic, Kabir and Tanimoto [11] suggest using a modeling method based on evolutionary game theory. In their study, Kabir et al. [12] briefly discussed how individuals use masks and other protective measures during an epidemic. These activities provide immediate benefits to the wearer and positively affect others. Alam et al. [13] , Higazy et al. [14] , Ullah et al. [15] , and Higazy and Alyami [16] conducted a comparative analysis of quarantine and isolation policies, transmission models concerning the ABO blood group, various pandemic scenarios excluding vaccination rates, and a genetic algorithm-based control strategy model for epidemic transmission. Das et al. [17] devise and evaluate a mathematical methodology to examine the patterns of COVID-19 transmission when comorbidity is present.

Islam et al. [18] quantify the impact of lockout and isolation measures on covid-19 situations in Bangladesh. Chowdhury et al. [19] examine quarantine and social-distancing measures’ theoretical and statistical effectiveness in potentially curbing the spread of new coronavirus cases in emerging or low-income countries. Das et al. [20] present a mathematical model that considers comorbidity to analyze the spread of COVID-19. They also offer a technique based on optimal control to minimize virus transmission. However, the studies mentioned above failed to consider the cumulative effect of the lockdown and vaccination scarcity despite the potential to yield intriguing findings. To address these deficiencies, we utilize a mathematical epidemiological model to guide our research on the combined impact of lockdown measures and vaccination shortages.

The vaccination program is widely employed to manage and eliminate the spread of the dangerous illness. Vaccination is an efficient method to reduce the spread of contagious illnesses and is crucial for implementing public health strategies. In addition, it confers direct protection to individuals, mitigates illness transmission, and offers herd immunity to the community [21] . Vaccination programs are often voluntary, allowing individuals to choose whether or not to get vaccinations [22] [23] . Several studies [24] [25] [26] [27] [28] have examined mass vaccination tactics, including mandatory vaccine programs. In addition, conventional and enduring vaccination methods are also essential; refer to references [24] [25] . Furthermore, the scarcity of vaccines is a highly significant problem in the management of infectious illnesses and has a direct influence on individuals’ choices regarding immunization. A shortage is improbable if an ample supply of immunizations reaches a distinct balance.

In addition, if there is a lack of accessible vaccinations to meet the demand, there might be an imbalance, even though scarcity is inherently self-perpetuating. When this situation arises, a significant number of individuals delay being vaccinated. Moreover, it exacerbates societal apprehension and unease during an epidemic, impacting individuals’ choices on vaccination. Under these circumstances, the interest among those chosen to receive the vaccine would intensify due to the increasing number of persons opting for vaccination. Choi et al. [29] and MacIntyre et al. [30] provide a mathematical modeling method for developing immunization plans for COVID-19 in Korea and NSW, Australia. Chen [31] does a theoretical analysis of volunteer immunizations and the issue of vaccine shortages. Li et al. [32] , Kahwati et al. [33] , Fairbrother et al. [34] , and Allison et al. [35] provided a concise analysis of how vaccination shortages affect several epidemic illnesses.

Ultimately, this analysis offers a structure for elucidating the collective impact of lockdown measures and vaccination campaigns, considering the conventional availability of vaccines. Furthermore, it portrays the collective actions that swiftly eliminated the plague from society. The suggested model elucidates the impact of reduced vaccination effectiveness and scarcity on the ultimate magnitude of the pandemic. Extended lockdown periods can amplify an individual’s inclination to venture outside, accelerating disease transmission. Hence, the synergistic effect of the efficient immunization campaign and non-pharmaceutical measures such as lockdowns, shutdowns, and states of emergency are more successful ways of managing the pandemic.

2. Model Formulation

The COVID-19 model is based on the standard SEIR model. The whole population is partitioned into nine subgroups (Figure 1): Susceptible $S\left(t\right)$ , Lock-down $L\left(t\right)$ , Exposed $E\left(t\right)$ , Infected $I\left(t\right)$ , Recovered $R\left(t\right)$ , Vaccinated $V\left(t\right)$ , Vaccinated exposed $E_V\left(t\right)$ , Vaccinated infected $I_V\left(t\right)$ and Vaccinated Recovered $R_V\left(t\right)$ .

Therefore, the proposed vaccination model is represented by the subsequent system of nonlinear ordinary differential equations:

$\frac{\text{d}S}{\text{d}t}=-\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\delta S\left(t\right)-lS\left(t\right)+l_{d}L\left(t\right),$ (1.1)

$\frac{\text{d}L}{\text{d}t}=lS\left(t\right)-\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-l_{d}L\left(t\right),$ (1.2)

$\frac{\text{d}E}{\text{d}t}=\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)+\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E\left(t\right),$ (1.3)

$\frac{\text{d}I}{\text{d}t}=\alpha E\left(t\right)-\gamma I\left(t\right),$ (1.4)

$\frac{\text{d}R}{\text{d}t}=\gamma I\left(t\right),$ (1.5)

$\frac{\text{d}V}{\text{d}t}=\delta S\left(t\right)-\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right),$ (1.6)

$\frac{\text{d}{E}_V}{\text{d}t}=\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E_V\left(t\right),$ (1.7)

$\frac{\text{d}{I}_V}{\text{d}t}=\alpha E_V\left(t\right)-\gamma I_V\left(t\right),$ (1.8)

$\frac{\text{d}{R}_V}{\text{d}t}=\gamma I_V\left(t\right).$ (1.9)

The total population,

$N\left(t\right)=S\left(t\right)+L\left(t\right)+E\left(t\right)+I\left(t\right)+R\left(t\right)+V\left(t\right)+E_V\left(t\right)+I_V\left(t\right)+R_V\left(t\right).$ (1.10)

Let us assume that if $x$ is the state variable vector, then we can write,

$x={\left(S\left(t\right),L\left(t\right),E\left(t\right),I\left(t\right),R\left(t\right),V\left(t\right),E_V\left(t\right),I_V\left(t\right),R_V\left(t\right)\right)}^{\prime }$ and $f:R^9\to R^9$ .

Figure 1. Schematic diagram of the model in which the population is divided into nine states: Susceptible $S\left(t\right)$ , Lockdown $L\left(t\right)$ , Exposed $E\left(t\right)$ , Infected $I\left(t\right)$ , Recovered $R\left(t\right)$ , Vaccinated $V\left(t\right)$ , Vaccinated exposed $E_V\left(t\right)$ , Vaccinated infected $I_V\left(t\right)$ and Vaccinated Recovered $R_V\left(t\right)$ . All parameters value and their biological significance are listed in Table 1.

The suggested model right side (Equations (1.1)-(1.9)) is thus a continuously differentiable function on $R^9$ . A unique elucidation of (1.1)-(1.9) occurs at every initial condition and continues for the maximum existence interval [36] . As a result, the suggested model has a clear biological significance. The model’s solution is also positive $\forall t\ge 0$ and bounded by the entire population $N\left(t\right)$ according to [37] (Equation (1.10)). Therefore, each compartment is regarded as one of the nine potential states at any time.

3. Vaccine Shortage

If $V_o$ is the amount of available vaccine, then

$\delta =\{\begin{array}{ll}0\hfill & \text{if}V\left(t\right)>V_o,\hfill \\ \delta \hfill & \text{otherwise}.\hfill \end{array}$ (2.1)

i) Susceptible individuals, $S\left(t\right)$ : In the beginning, the susceptible part of the overall population is exposed to the infected individuals (Equation (1.1)). Those who lose immunity due to an earlier infection add to the vulnerable group and are reduced through vaccination (moving to class V with an efficacy rate $\left(1-\eta \right)$ at the rate, δ), following the lockdown level l, infection (moving to class E).

ii) Lockdown individuals, $L\left(t\right)$ : Lockdown state members are the individuals who have adhered to the lockdown-related guidelines. It refers to susceptible people who stay at home and are not infected by the virus. The lockdown maintenance or obedience factor $\left(1-q\right)$ has a crucial role in the prevention of infection. Here, $q=0$ this indicates that no people maintain lockdown rules, and $q=1$ , the lockdown policy works appropriately. Using the Heaviside function, we assess the lockdown open and shut mechanism.

$l\left[d_{start},d_{end}\right]=\{\begin{array}{ll}0,\hfill & t\notin \left[d_{start},d_{end}\right],\hfill \\ 1,\hfill & t\in \left[d_{start},d_{end}\right].\hfill \end{array}$ (2.2)

where, $d_{start}$ = lockdown starting time, and $d_{end}$ = lockdown ending time.

iii) Exposed individuals, $E\left(t\right)$ : Infected people who are still vulnerable after being vaccinated increase the number of the exposed and not vaccinated yet; also, the lockdown maintenance or obedience factor. The onset of infection reduces the exposed population (moving to class $I\left(t\right)$ ) at the rate α. During the exposed phase, humans maintain a low infectivity level. This compartment technically reflects the slightly infectious stage.

iv) Infected individuals, $I\left(t\right)$ : A fraction of exposed people evolve an infected class at the rate α increases the population of an infective class. Recovery from the disease at the rate γ reduces the number of infective persons.

v) Recovered individuals, $R\left(t\right)$ : The recovered population increased from infective class at the rate γ.

vi) Vaccinated individuals, $V\left(t\right)$ : The number of people who have been vaccinated has grown as more people have been immunized. The vaccinated people were reduced by a factor $\left(1-\eta \right)\beta$ with $0\le \eta \le 1$ .

vii) Vaccinated exposed individuals, $E_v\left(t\right)$ : The vaccinated exposed population increases if the mass susceptible individuals take part under the vaccinated program and are reduced by the onset of infection (moving to class $I_V\left(t\right)$ ) at the rate α.

viii) Vaccinated infected individuals, $I_V\left(t\right)$ : The population of vaccinated infective class increased by a fraction of vaccinated exposed individuals becoming infective class at the rate α. The vaccinated infected individuals are reduced by recovery from the disease at the rate γ.

ix) Vaccinated recovered individuals, $R_V\left(t\right)$ : The recovered population increased from the vaccinated infective class at the rate γ.

Table 1 explains the biological significance of the parameters and the values preferred for them.

4. Mathematical Analysis

The present section studies the proposed model’s positivity and boundedness, basic (R₀) reproduction number, the threshold of points of equilibrium (disease-free and endemic), stability of disease-free (E₀) and endemic ( $E_{\ast }$ ) equilibrium point, effective (R_e) reproduction number, the relation between R₀ and crucial harmonies of vaccination, adequate strength number, the first and second-order Lyapunov functions (LF), existence, and uniqueness theorem for validating the model’s stability.

Table 1. List of parameters, variables, and their biological meanings.

4.1. Models Positivity and Boundedness [38] [39] [40]

Let the assumed state vector be represented as $x=\left(x_1,x_2,x_3,x_4,x_5,x_6,x_7,x_8,x_9\right)=\left(S,L,E,I,R,V,E_V,I_V,R_V\right)$ . Then Equations (1.1)-(1.9) are possible to signify in the form of some initial set of values

$x^{\prime }=f\left(x,t\right),x\left(0\right)=x_0$ (2.3)

It is evident that the above-defined function f is locally Lipschitz for the first argument and continuous for the second in ${ℝ}^9\times ℝ$ , which implies that the solution $x\left(t\right)$ holds in $t\in \left(0,T\right)$ for some T, according to the system of nonlinear ODEs existence and uniqueness theorem [38] [41] . For the invariance of sets under a flow, we used the Bony-Brezis theorem [42] . The theorem applies to smooth manifolds. There are some corner points at the boundary of the invariant set. Shortly, by making an argument that if we have an initial condition at the corner points, it cannot go outside; for instance, the direction of the flow is towards the interior at those points. Thus, T might be re-defined as the supremum

overall mentioned periods. Finally, for $T<\infty$ , ${\lim }_{t\to T^{-}}\Vert x\left(t\right)\Vert \to \infty$ .

We have to show that under the flow $x\left(t\right)$ for $t\in \left(0,T\right)$ , if

$x_0\in \Lambda =\left\{x\in {ℝ}^9:x_i\ge 0\text{for}i=1,2,3,\cdots ,9,{\displaystyle \underset{i=1}{\overset{9}{\sum }}{x}_i}\le N\left(0\right)\right\},$

the set Λ is positively invariant, which concludes that $T=\infty$ , i.e., for initial conditions $x_0\in \text{Λ}$ , the solution $x\left(t\right)\in \text{Λ}$ holds globally with time [38] [40] [41] .

Theorem 1: The closed set

$\Lambda :=\left\{x=\left(S,L,E,I,R,V,E_V,I_V,R_V\right)\in {ℝ}^9:x_i\ge 0\text{for}i=1,2,3,\cdots ,9,{\displaystyle \underset{i=1}{\overset{9}{\sum }}{x}_i}\le N\left(0\right)\right\},$

is entirely uniform under the flow generated by the Equations (1.1)-(1.9). As a result, the solution $x\left(t\right)\in \text{Λ}$ occurs globally in the time given initial conditions $x_0\in \text{Λ}$ .

Proof: Let the boundary segment is ${\Pi }_i,i=1,2,\cdots ,10$

${\Pi }_i=\left\{x\in \Lambda :x_i=0\right\},i=1,2,\cdots ,9$

${\Pi }_{10}=\left\{x\in \Lambda :{\displaystyle \underset{i=1}{\overset{10}{\sum }}{x}_i}=0\right\}.$

It is evident that $\partial \Lambda ={\displaystyle {\cup }_{i=1}^{10}{\Pi }_i}$ .

To complete the proof of the above invariance of the set Λ, it is enough to prove $\forall$ inward normal, $n\cdot x^{\prime }\left(t\right)\ge 0$ on $\partial \Lambda$ . The inward normal on ${\Pi }_i$ for $i=1,2,\cdots ,9$ is unambiguously provided by $n_i=v_i=\left[0,\cdots ,1,\cdots ,0\right]$ , where only the i-component is nonzero, but an inward normal on ${\Pi }_{10}$ is given by $n_{10}=\left(-1,-1,\cdots ,-1\right)$ .

Now, on ${\Pi }_i$ for $i=1,2,\cdots ,9$

$v_1\cdot x^{\prime }=l_d{x}_2\ge 0,\forall x\in {\Pi }_1,$

$v_2\cdot x^{\prime }=l{x}_1\ge 0,\forall x\in {\Pi }_2,$

$v_3\cdot x^{\prime }=\beta x_1\left(x_4+x_8\right)+\left(1-q\right)\beta x_2\left(x_4+x_8\right)\ge 0,\forall x\in {\Pi }_3,$

$v_4\cdot x^{\prime }=\alpha x_3\ge 0,\forall x\in {\Pi }_4,$

$v_5\cdot x^{\prime }=\gamma x_4\ge 0,\forall x\in {\Pi }_5,$

$v_6\cdot x^{\prime }=\delta x_1\ge 0,\forall x\in {\Pi }_6,$

$v_7\cdot x^{\prime }=\left(1-\eta \right)\beta x_6\left(x_4+x_8\right)\ge 0,\forall x\in {\Pi }_7,$

$v_8\cdot x^{\prime }=\alpha x_7\ge 0,\forall x\in {\Pi }_8,$

$v_9\cdot x^{\prime }=\gamma x_8\ge 0,\forall x\in {\Pi }_9,$

while $N={\displaystyle {\sum }_{i=1}^9{x}_i}$ is readily seen to satisfy $\frac{\text{d}N}{\text{d}t}=0$ , which follows that on ${\Pi }_{10}$ such that

$n_{10}\cdot x^{\prime }=0.$

Thus, for given initial conditions $x_0\in \Lambda$ , solution $x\left(t\right)\in \Lambda$ holds globally over time on the positively invariant domain Λ.

4.2. Derivation of the Basic Reproduction Number (R₀)

Before commencing the formal study, we will provide a comprehensive explanation of the methodology used to calculate the basic reproduction number (R₀), which is crucial in epidemiological modeling due to its significant contribution to understanding stability conditions. Research has shown that if this numerical value is less than one, it indicates stable circumstances, whereas values more than one indicate unstable situations. However, it has been stressed that the value may be determined using other methods, with the next-generation matrix techniques [43] being widely recognized. Thus, the detailed process

$F{V}^{-1}\left(\equiv F_1{V}_1^{-1},F_2{V}_2^{-1}\right)$ of finding the value of the basic reproduction number (R₀) is

$F_1=\left[\begin{array}{cc}0& \beta +\left(1-q\right)\beta \\ 0& 0\end{array}\right],V_1=\left[\begin{array}{cc}\alpha & 0\\ -\alpha & \gamma \end{array}\right].$

$\therefore F_1{V}_1^{-1}=\frac{1}{\alpha \gamma }\left[\begin{array}{cc}\alpha \beta & \alpha \left\{\beta +\left(1-q\right)\beta \right\}\\ 0& 0\end{array}\right].$

The eigenvalues of $F_1{V}_1^{-1}$ are ${\lambda }_1=\frac{\beta }{\gamma }$ and ${\lambda }_2=\frac{\beta +\left(1-q\right)\beta }{\gamma }$ .

Correspondingly, for the vaccination part,

$F_2=\left[\begin{array}{cc}0& \beta \left(1-\eta \right)\\ 0& 0\end{array}\right],V_2=\left[\begin{array}{cc}\alpha & 0\\ -\alpha & \gamma \end{array}\right].$

$\therefore F_2{V}_2^{-1}=\frac{1}{\alpha \gamma }\left[\begin{array}{cc}\alpha \beta & \alpha \left(1-\eta \right)\beta \\ 0& 0\end{array}\right].$

The eigenvalues of $F_2{V}_2^{-1}$ are ${\lambda }_3=\frac{\beta }{\gamma }$ and ${\lambda }_4=\frac{\left(1-\eta \right)\beta }{\gamma }$ .

Thus, the basic reproduction number (R₀) is

$R_0=\frac{\beta }{\gamma }+\frac{\left(1-q\right)\beta }{\gamma }+\frac{\left(1-\eta \right)\beta }{\gamma }.$ (2.4)

4.3. Analysis of Disease-Free and Endemic Equilibrium Point Threshold

For this model, the disease-free equilibrium point threshold is the solution of the nonlinear system of Equations (1.1)-(1.9). Thus,

$\frac{\text{d}S}{\text{d}t}=\frac{\text{d}L}{\text{d}t}=\frac{\text{d}E}{\text{d}t}=\frac{\text{d}I}{\text{d}t}=\frac{\text{d}R}{\text{d}t}=\frac{\text{d}V}{\text{d}t}=\frac{E_V}{\text{d}t}=\frac{I_V}{\text{d}t}=\frac{R_V}{\text{d}t}=0,$

gives $E_0=\left(N\left(=1\right),0,0,0,0,0,0,0,0\right)$, the disease-free equilibrium point threshold of the current model.

The proposed model endemic equilibrium point $E_{*}=\left(S^{*},L^{*},E^{*},I^{*},R^{*},V^{*},E_{v^{*}},I_{v^{*}},R_{v^{*}}\right)$ is a solution of the following system:

$0=-\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\delta S\left(t\right)-lS\left(t\right)+l_{d}L\left(t\right),$

$0=lS\left(t\right)-\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-l_{d}L\left(t\right),$

$0=\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)+\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E\left(t\right),$

$0=\alpha E\left(t\right)-\gamma I\left(t\right),$

$0=\gamma I\left(t\right),$

$0=\delta S\left(t\right)-\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right),$

$0=\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E_V\left(t\right),$

$0=\alpha E_V\left(t\right)-\gamma I_V\left(t\right),$

$0=\gamma I_V\left(t\right).$

Solving the above system of equations setting $I\ne 0,I_V\ne 0$ , the endemic equilibrium point is

$E_{*}=\left(S^{*},L^{*},E^{*},I^{*},R^{*},V^{*},E_{v^{*}},I_{v^{*}},R_{v^{*}}\right)$ ,

where

$\begin{array}{l}{S}^{*}=\frac{\gamma I^{*}}{\beta I^{*}+\delta +l},L^{*}=\frac{\gamma }{\left(1-q\right)\beta }-\frac{\gamma I^{*}}{\left(1-q\right)\left(\beta I^{*}+\delta +l\right)},E^{*}=\frac{l_d{L}^{*}}{\alpha },I^{*}=\frac{l_d{L}^{*}}{\gamma },\\ R^{*}=l_d{L}^{*},V^{*}=\frac{\gamma }{\left(1-\eta \right)\beta },E_{v^{*}}=\frac{\delta S^{*}}{\alpha },I_{v^{*}}=\frac{\delta S^{*}}{\gamma },R_{v^{*}}=\delta S^{*}.\end{array}$

4.4. Stability of Disease-Free (E₀) Equilibrium Point

In this section, we will show that for $R_0<1$ the disease-free $(E_0)$ equilibrium point isasymptoticallystable locally. The disease will be eradicated and continue biologically in society when the primary reproduction number is less than and more significant than unity.

Theorem 2: If $R_0<1$ , then the unique disease-free equilibrium E₀ is locally asymptotically stable. If $R_0>1$ , the unique disease-free equilibrium is unstable.

Proof: To validate the local stability, the Jacobian matrix of the proposed system (1.1)-(1.9) is

$J=\left[\begin{array}{ccccccccc}-\beta \left(I+I_V\right)-\delta -l& l_d& 0& -\beta S& 0& 0& 0& -\beta S& 0\\ l& -\beta \left(1-q\right)\left(I+I_V\right)-l_d& 0& -\beta \left(1-q\right)L& 0& 0& 0& -\beta \left(1-q\right)L& 0\\ \beta \left(I+I_V\right)& \beta \left(1-q\right)\left(I+I_V\right)& -\alpha & \beta S+\beta \left(1-q\right)L& 0& 0& 0& \beta S+\beta \left(1-q\right)L& 0\\ 0& 0& \alpha & -\gamma & 0& 0& 0& 0& 0\\ 0& 0& 0& \gamma & 0& 0& 0& 0& 0\\ \delta & 0& 0& -\beta \left(1-\eta \right)V& 0& -\beta \left(1-\eta \right)\left(I+I_V\right)& 0& -\beta \left(1-\eta \right)V& 0\\ 0& 0& 0& \beta \left(1-\eta \right)V& 0& \beta \left(1-\eta \right)\left(I+I_V\right)& -\alpha & \beta \left(1-\eta \right)V& 0\\ 0& 0& 0& 0& 0& 0& \alpha & -\gamma & 0\\ 0& 0& 0& 0& 0& 0& 0& \gamma & 0\end{array}\right].$

At the DFE point $E_0$ , we have

$J\left(E_0\right)=\left[\begin{array}{ccccccccc}-\delta -l& l_d& 0& -\beta & 0& 0& 0& -\beta & 0\\ l& -l_d& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& -\alpha & \beta & 0& 0& 0& \beta & 0\\ 0& 0& \alpha & -\gamma & 0& 0& 0& 0& 0\\ 0& 0& 0& \gamma & 0& 0& 0& 0& 0\\ \delta & 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& -\alpha & 0& 0\\ 0& 0& 0& 0& 0& 0& \alpha & -\gamma & 0\\ 0& 0& 0& 0& 0& 0& 0& \gamma & 0\end{array}\right]$

The characteristic equation $\left|J\left(E_0\right)-\lambda I\right|=0$ , has nine roots, which are

${\lambda }_1={\lambda }_2={\lambda }_3=0,{\lambda }_4=-\alpha ,{\lambda }_5=-\gamma ,$

${\lambda }_6=\frac{1}{2}\left(-\alpha -\gamma -\sqrt{{\alpha }^2+4\alpha \beta -2\alpha \gamma +{\gamma }^2}\right),$

${\lambda }_7=\frac{1}{2}\left(-\alpha -\gamma +\sqrt{{\alpha }^2+4\alpha \beta -2\alpha \gamma +{\gamma }^2}\right),$

${\lambda }_8=\frac{1}{2}\left(-l_d-l-\delta -\sqrt{{\left(l_d+l+\delta \right)}^2+4l_d\delta }\right)\text{and}$

${\lambda }_9=\frac{1}{2}\left(-l_d-l-\delta +\sqrt{{\left(l_d+l+\delta \right)}^2+4l_d\delta }\right).$

i.e., all nine eigenvalues are less than or equal to zero. Therefore, according to the Routh-Hurwitz criterion, the model is locally asymptotically stable at the unique DFE point whenever $R_0<1$ and unstable whenever $R_0>1$ .

4.5. Derivation of the Effective Reproduction Number (R_e)

In general, a population in the real world can infrequently be entirely susceptible to infection. Some contacts will be immune because of previous infections, which have provided long-term immunity or are a consequence of prior vaccinations. Consequently, all communications will not infect, and the mean per secondary case of disease is less than the number of basic reproduction numbers. Thus, the effective reproductive number (R_e) is the mean per infectious case in a population of susceptible and non-susceptible people. When $R_e>1$ , the number of cases will increase and start an epidemic. When $R_e=1$ , the disease is endemic, and when $R_e<1$ , there will be a decline in the number of cases. The effective reproduction number is sensitive to the multiple of the basic reproduction number and the fraction of the host population. Therefore,

$R_e=\frac{\beta }{\gamma }S\left(t\right)+\frac{\left(1-q\right)\beta }{\gamma }L\left(t\right)+\frac{\left(1-\eta \right)\beta }{\gamma }V\left(t\right).$ (2.5)

4.6. Relation between R₀ and Crucial Harmonies of Vaccination

In the recommended model, $\frac{\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)}{N}$ is used to compute the sum of secondary infections of susceptible and vaccinated infected people per unit of time. On top of that, it is equal to $\frac{\beta S\left(t\right)}{N}$ for the single infected individual. Furthermore, $\frac{1}{\gamma }$ and $\frac{\beta }{\gamma }$ indicates the lifespan of an infectious individual and the total secondary infections of susceptible individuals that may produce one infected individual in a disease-free society, respectively, whereas, for vaccinated individuals, both of the values are $\frac{1}{\gamma }$ and $\frac{\left(1-\eta \right)\beta }{\gamma }$ . Analogously, secondary infections of vaccinated individuals equal $\frac{\beta S\left(t\right)I_V\left(t\right)}{N}$ and per unit of time, the sum is $\frac{\beta S\left(t\right)}{N}$ . Realistically, due to the mass vaccination rate $\frac{1}{\delta }$ , the basic reproduction number must be a decreasing function. Opposite scenarios for less vaccine efficacy rate $\left(1-\eta \right)$ .

4.7. Adequate Strength Number

Besides the reproduction number commonly employed in epidemiology to determine whether or not a disease will spread, we introduce the adequate strength number to measure wave tendency in the epidemiology model. We use the well-established next-generation matrix technique to calculate the adequate strength number by assuming the total population is constant (N). Then,

$\beta SI\equiv \frac{\beta SI}{N},E\equiv \frac{E}{N}.$

Now,

$\begin{array}{l}\frac{\partial }{\partial I}\left\{\beta S\frac{I}{N}+\left(1-q\right)\beta L\frac{I}{N}-\alpha \frac{E}{N}\right\}\\ =\left\{\beta S+\left(1-q\right)\beta L\right\}\left(\frac{N-\stackrel{\dot{}}{N}I}{N^2}\right)-\alpha \left(\frac{N\stackrel{\dot{}}{E}-\stackrel{\dot{}}{N}E}{N^2}\right)\\ =\left\{\beta S+\left(1-q\right)\beta L\right\}\left(\frac{N-\stackrel{\dot{}}{N}I}{N^2}\right)-\alpha \frac{\left\{\beta S+\left(1-q\right)\beta L\right\}I-\alpha E-\stackrel{\dot{}}{N}E}{N^2}\\ =\left\{\beta S+\left(1-q\right)\beta L\right\}\left(\frac{N-\stackrel{\dot{}}{N}I}{N^2}\right)-\alpha \left\{\beta S+\left(1-q\right)\beta L\right\}\frac{I}{N^2}+\frac{{\alpha }^{2}E}{N^2}+\alpha \frac{\stackrel{\dot{}}{N}E}{N^2}.\end{array}$ (2.6)

Similarly,

$\begin{array}{l}\frac{\partial }{\partial I_V}\left\{\beta S\frac{I_V}{N}+\left(1-q\right)\beta L\frac{I_V}{N}-\alpha \frac{E}{N}\right\}\\ =\left\{\beta S+\left(1-q\right)\beta L\right\}\left(\frac{N-\stackrel{\dot{}}{N}{I}_V}{N^2}\right)-\alpha \left\{\beta S+\left(1-q\right)\beta L\right\}\frac{I_V}{N^2}+\frac{{\alpha }^{2}E}{N^2}+\alpha \frac{\stackrel{\dot{}}{N}E}{N^2}.\end{array}$ (2.7)

$\begin{array}{l}\frac{\partial }{\partial I}\left\{\left(1-\eta \right)\beta V\frac{I}{N}-\alpha \frac{E_V}{N}\right\}\\ =\left(1-\eta \right)\beta V\left(\frac{N-\stackrel{\dot{}}{N}I}{N^2}\right)-\left(1-\eta \right)\alpha \beta V\frac{I}{N^2}+\frac{{\alpha }^2{E}_V}{N^2}+\alpha \frac{\stackrel{\dot{}}{N}{E}_V}{N^2}.\end{array}$ (2.8)

$\begin{array}{l}\frac{\partial }{\partial I_V}\left\{\left(1-\eta \right)\beta V\frac{I_V}{N}-\alpha \frac{E_V}{N}\right\}\\ =\left(1-\eta \right)\beta V\left(\frac{N-\stackrel{\dot{}}{N}{I}_V}{N^2}\right)-\left(1-\eta \right)\alpha \beta V\frac{I_V}{N^2}+\frac{{\alpha }^2{E}_V}{N^2}+\alpha \frac{\stackrel{\dot{}}{N}{E}_V}{N^2}.\end{array}$ (2.9)

In this case,

$F_{SN\left(1\right)}=\left[\begin{array}{cc}0& -\frac{\beta +\left(1-q\right)\beta }{N^2}-\frac{\alpha \beta +\left(1-q\right)\alpha \beta }{N^2}\\ 0& 0\end{array}\right],$ (2.10)

$F_{SN\left(2\right)}=\left[\begin{array}{cc}0& -\frac{\left(1-\eta \right)\beta }{N^2}-\frac{\left(1-\eta \right)\alpha \beta }{N^2}\\ 0& 0\end{array}\right],$ (2.11)

and

$V_{SN\left(1\right)}=\left[\begin{array}{cc}-\frac{{\alpha }^2}{N^2}& 0\\ -\alpha & \gamma \end{array}\right],$ (2.12)

$V_{SN\left(2\right)}=\left[\begin{array}{cc}-\frac{{\alpha }^2}{N^2}& 0\\ -\alpha & \gamma \end{array}\right].$ (2.13)

$\therefore F_{SN\left(1\right)}{V}_{SN\left(1\right)}^{-1}=-\frac{N^2}{\gamma {\alpha }^2}\left[\begin{array}{cc}-a_{11}& a_{12}\\ 0& 0\end{array}\right],$ (2.14)

where

$a_{11}=\frac{\alpha \beta +\left(1-q\right)\alpha \beta }{N^2}+\frac{{\alpha }^2\beta +\left(1-q\right){\alpha }^2\beta }{N^2},$

$a_{12}=\left\{\frac{\beta +\left(1-q\right)\beta }{N^2}+\frac{\alpha \beta +\left(1-q\right)\alpha \beta }{N^2}\right\}\frac{{\alpha }^2}{N^2}.$

Analogously,

$F_{SN\left(2\right)}{V}_{SN\left(2\right)}^{-1}=-\frac{N^2}{\gamma {\alpha }^2}\left[\begin{array}{cc}-\frac{\left(1-\eta \right)\alpha \beta }{N^2}-\frac{\left(1-\eta \right){\alpha }^2\beta }{N^2}& \left\{\frac{\left(1-\eta \right)\beta }{N^2}+\frac{\left(1-\eta \right)\alpha \beta }{N^2}\right\}\frac{{\alpha }^2}{N^2}\\ 0& 0\end{array}\right].$ (2.15)

Thus, the four Eigenvalues are

${\lambda }_{SN\left(1\right)}=\frac{\beta +\left(1-q\right)\beta }{\alpha \gamma }+\frac{\beta +\left(1-q\right)\beta }{\gamma },{\lambda }_{SN\left(2\right)}=-\frac{\beta +\left(1-q\right)\beta }{\gamma N^2}-\frac{\alpha \beta +\left(1-q\right)\alpha \beta }{\gamma N^2},$

${\lambda }_{SN\left(3\right)}=\frac{\left(1-\eta \right)\beta }{\alpha \gamma }+\frac{\left(1-\eta \right)\beta }{\gamma },{\lambda }_{SN\left(4\right)}=-\frac{\left(1-\eta \right)\beta }{\gamma N^2}-\frac{\left(1-\eta \right)\alpha \beta }{\gamma N^2}.$

At that juncture

$SN=\frac{\beta +\left(1-q\right)\beta }{\alpha \gamma }+\frac{\beta +\left(1-q\right)\beta }{\gamma }+\frac{\left(1-\eta \right)\beta }{\alpha \gamma }+\frac{\left(1-\eta \right)\beta }{\gamma }.$

Or,

$SN=\frac{\beta }{\alpha \gamma }\left(3-q-\eta \right)\left(1+\alpha \right).$ (2.16)

As stated above, the reproduction number (effective) must reflect the stability of the disease-free equilibrium state but failed to present wave tendency. According to equation (2.16), the strength number (SN) acts as a positive number (SN>0) for the parameters q, η, and α, all ranging from 0 to 1. Thus, SN > 0 determines that the spread has enough power to begin the revival stage, meaning more than one wave will occur.

4.8. The First Derivative of the Lyapunov Function (LF)

Let us assume models independent variables for the endemic LF, $\left\{S,L,E,I,R,V,E_V,I_V,R_V\right\}$ , $L_f<0$ is the detrimental equilibrium point $E_{*}$ .

Theorem 3: For the value of the basic reproductive number $R_0-1>0$ , the proposed $SLEIRV{E}_V{I}_V{R}_V$ model endemic equilibrium point $E_{*}$ is globally asymptotically stable.

Proof: The LF may be represented as follows to prove the theorem above:

$\begin{array}{l}{L}_f\left(S,L,E,I,R,V,E_V,I_V,R_V\right)\\ =\left(S-S^{*}-S^{*}\log \frac{S^{*}}{S}\right)+\left(L-L^{*}-L^{*}\log \frac{L^{*}}{L}\right)+\left(E-E^{*}-E^{*}\log \frac{E^{*}}{E}\right)\\ +\left(I-I^{*}-I^{*}\log \frac{I^{*}}{I}\right)+\left(R-R^{*}-R^{*}\log \frac{R^{*}}{R}\right)+\left(V-V^{*}-V\log \frac{V^{*}}{V}\right)\\ +\left(E_V-E_V^{*}-E_V^{*}\log \frac{E_V^{*}}{E_V}\right)+\left(I_V-I_V^{*}-I_V^{*}\log \frac{I_V^{*}}{I_V}\right)+\left(R_V-R_V^{*}-R_V^{*}\log \frac{R_V^{*}}{R_V}\right).\end{array}$ (2.17)

Differentiating both sides in terms of t, we get,

$\begin{array}{l}\frac{\text{d}{L}_f}{\text{d}t}={\stackrel{\dot{}}{L}}_f=\left(\frac{S-S^{*}}{S}\right)\stackrel{\dot{}}{S}+\left(\frac{L-L^{*}}{L}\right)\stackrel{\dot{}}{L}+\left(\frac{E-E^{*}}{E}\right)\stackrel{\dot{}}{E}+\left(\frac{I-I^{*}}{I}\right)\stackrel{\dot{}}{I}+\left(\frac{R-R^{*}}{R}\right)\stackrel{\dot{}}{R}\\ +\left(\frac{V-V^{*}}{V}\right)\stackrel{\dot{}}{V}+\left(\frac{E_V-E_V^{*}}{E_V}\right){\stackrel{\dot{}}{E}}_V+\left(\frac{I_V-I_V^{*}}{I_V}\right){\stackrel{\dot{}}{I}}_V+\left(\frac{R_V-R_V^{*}}{R_V}\right){\stackrel{\dot{}}{R}}_V.\end{array}$ (2.18)

From Equations (1.1)-(1.9), putting the value of $\stackrel{\dot{}}{S},\stackrel{\dot{}}{L},\stackrel{\dot{}}{E},\stackrel{\dot{}}{I},\stackrel{\dot{}}{V},{\stackrel{\dot{}}{E}}_V,{\stackrel{\dot{}}{I}}_V,{\stackrel{\dot{}}{R}}_V$ in equation (2.18), we have,

$\begin{array}{c}\frac{\text{d}{L}_f}{\text{d}t}=\left(\frac{S-S^{*}}{S}\right)\left\{-\beta S\left(I+I_V\right)-\delta S-lS+l_{d}L\right\}\\ +\left(\frac{L-L^{*}}{L}\right)\left\{lS-\left(1-q\right)\beta L\left(I+I_V\right)-l_{d}L\right\}\\ +\left(\frac{E-E^{*}}{E}\right)\left\{\beta S\left(I+I_V\right)+\left(1-q\right)\beta L\left(I+I_V\right)-\alpha E\right\}\end{array}$

$\begin{array}{c}+\left(\frac{I-I^{*}}{I}\right)\left(\alpha E-\gamma I\right)+\left(\frac{R-R^{*}}{R}\right)\gamma I\\ +\left(\frac{V-V^{*}}{V}\right)\left\{\delta S-\left(1-\eta \right)\beta V\left(I+I_V\right)\right\}\\ +\left(\frac{E_V-E_V^{*}}{E_V}\right)\left\{\left(1-\eta \right)\beta V\left(I+I_V\right)-\alpha E_V\right\}\\ +\left(\frac{I_V-I_V^{*}}{I_V}\right)\left(\alpha E_V-\gamma I_V\right)+\left(\frac{R_V-R_V^{*}}{R_V}\right)\gamma I_V.\end{array}$ (2.19)

In Equation (2.19), substitute $S-S^{*}$ , $L-L^{*}$ , $E-E^{*}$ , $I-I^{*}$ , $V-V^{*}$ , $E_V-E_V^{*}$ , $I_V-I_V^{*}$ , $R_V-R_V^{*}$ instead of $S,L,E,I,R,V,E_V,I_V,R_V$ . Then,

$\begin{array}{c}\frac{\text{d}{L}_f}{\text{d}t}=\left(\frac{S-S^{*}}{S}\right)[-\beta \left(S-S^{*}\right)\left\{\left(I-I^{*}\right)+\left(I_V-I_V^{*}\right)\right\}-\delta \left(S-S^{*}\right)-l\left(S-S^{*}\right)\\ +l_d\left(L-L^{*}\right)]+\left(\frac{L-L^{*}}{L}\right)[l\left(S-S^{*}\right)-\left(1-q\right)\beta \left(L-L^{*}\right)\{\left(I-I^{*}\right)\\ +\left(I_V-I_V^{*}\right)\}-l_d\left(L-L^{*}\right)]+\left(\frac{E-E^{*}}{E}\right)[\beta \left(S-S^{*}\right)\left\{\left(I-I^{*}\right)+\left(I_V-I_V^{*}\right)\right\}\\ +\left(1-q\right)\beta \left(L-L^{*}\right)\left\{\left(I-I^{*}\right)+\left(I_V-I_V^{*}\right)\right\}-\alpha \left(E-E^{*}\right)]\end{array}$

$\begin{array}{c}+\left(\frac{I-I^{*}}{I}\right)\left\{\alpha \left(E-E^{*}\right)-\gamma \left(I-I^{*}\right)\right\}+\left(\frac{R-R^{*}}{R}\right)\gamma \left(I-I^{*}\right)\\ +\left(\frac{V-V^{*}}{V}\right)\left[\delta \left(S-S^{*}\right)-\left(1-\eta \right)\beta \left(V-V^{*}\right)\left\{\left(I-I^{*}\right)+\left(I_V-I_V^{*}\right)\right\}\right]\\ +\left(\frac{E_V-E_V^{*}}{E_V}\right)\left[\left(1-\eta \right)\beta \left(V-V^{*}\right)\left\{\left(I-I^{*}\right)+\left(I_V-I_V^{*}\right)\right\}-\alpha \left(E_V-E_V^{*}\right)\right]\\ +\left(\frac{I_V-I_V^{*}}{I_V}\right)\left\{\alpha \left(E_V-E_V^{*}\right)-\gamma \left(I_V-I_V^{*}\right)\right\}+\left(\frac{R_V-R_V^{*}}{R_V}\right)\gamma \left(I_V-I_V^{*}\right).\end{array}$ (2.20)

After simplifying, we may write,

$\frac{\text{d}{L}_f}{\text{d}t}={\Pi }_1-{\Pi }_2,$ (2.21)

where

$\begin{array}{c}{\Pi }_1=\beta \frac{{\left(S-S^{*}\right)}^2}{S}{I}^{*}+\beta \frac{{\left(S-S^{*}\right)}^2}{S}{I}_V^{*}+l_{d}L+l_d\frac{S^{*}}{S}{L}^{*}+lS+l\frac{L^{*}}{L}{S}^{*}+\beta \left(1-q\right)L{I}^{*}\\ +\beta \left(1-q\right)L{I}_V^{*}+\beta \left(1-q\right)L^{*}I+\beta \left(1-q\right)L^{*}{I}_V+\beta \left(1-q\right)\frac{E^{*}}{E}LI\\ +\beta \left(1-q\right)\frac{E^{*}}{E}L{I}_V+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}^{*}+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}_V^{*}+\beta SI+\beta S{I}_V\\ +\beta S^{*}{I}^{*}+\beta S^{*}{I}_V^{*}+\beta \frac{E^{*}}{E}S{I}^{*}+\beta \frac{E^{*}}{E}S{I}_V^{*}+\beta \frac{E^{*}}{E}{S}^{*}I+\beta \frac{E^{*}}{E}{S}^{*}{I}_V\end{array}$

$\begin{array}{c}+\beta \left(1-q\right)LI+\beta \left(1-q\right)L{I}_V+\beta \left(1-q\right)L^{*}{I}^{*}+\beta \left(1-q\right)L^{*}{I}_V^{*}+\beta \left(1-q\right)\frac{E^{*}}{E}L{I}^{*}\\ +\beta \left(1-q\right)\frac{E^{*}}{E}L{I}_V^{*}+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}I+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}_V+\alpha E+\alpha \frac{I^{*}}{I}{E}^{*}+\gamma I\\ +\gamma \frac{R^{*}}{R}{I}^{*}+\delta S+\delta \frac{V^{*}}{V}{S}^{*}+\beta \left(1-\eta \right)\frac{{\left(V-V^{*}\right)}^2}{V}{I}^{*}+\beta \left(1-\eta \right)\frac{{\left(V-V^{*}\right)}^2}{V}{I}_V^{*}\\ +\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}V{I}^{*}+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}V{I}_V^{*}+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}{V}^{*}I+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}{V}^{*}{I}_V\\ +\gamma I_V+\gamma \frac{R_V^{*}}{R_V}{I}_V^{*}+\alpha E_V+\alpha \frac{I_V^{*}}{I_V}{E}_V^{*},\end{array}$

and

$\begin{array}{c}{\Pi }_2=\beta \frac{{\left(S-S^{*}\right)}^2}{S}I+\beta \frac{{\left(S-S^{*}\right)}^2}{S}{I}_V+\left(\delta +l\right)\frac{{\left(S-S^{*}\right)}^2}{S}+l_d{L}^{*}+l_d\frac{S^{*}}{S}L+l{S}^{*}\\ +l\frac{L^{*}}{L}S+\beta \left(1-q\right)LI+\beta \left(1-q\right)L{I}_V+\beta \left(1-q\right)L^{*}{I}^{*}+\beta \left(1-q\right)L^{*}{I}_V^{*}\\ +\beta \left(1-q\right)\frac{E^{*}}{E}L{I}^{*}+\beta \left(1-q\right)\frac{E^{*}}{E}L{I}_V^{*}+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}I+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}_V\\ +l_d\frac{{\left(L-L^{*}\right)}^2}{L}+\beta S{I}^{*}+\beta S{I}_V^{*}+\beta S^{*}I+\beta S^{*}{I}_V+\beta \frac{E^{*}}{E}SI+\beta \frac{E^{*}}{E}S{I}_V\\ +\beta \frac{E^{*}}{E}{S}^{*}{I}^{*}+\beta \frac{E^{*}}{E}{S}^{*}{I}_V^{*}+\beta \left(1-q\right)L{I}^{*}+\beta \left(1-q\right)L{I}_V^{*}+\beta \left(1-q\right)L^{*}I\end{array}$

$\begin{array}{c}+\beta \left(1-q\right)L^{*}{I}_V+\beta \left(1-q\right)\frac{E^{*}}{E}LI+\beta \left(1-q\right)\frac{E^{*}}{E}L{I}_V+\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}^{*}\\ +\beta \left(1-q\right)\frac{E^{*}}{E}{L}^{*}{I}_V^{*}+\alpha \frac{{\left(E-E^{*}\right)}^2}{E}+\alpha E^{*}+\alpha \frac{I^{*}}{I}E+\gamma \frac{{\left(I-I^{*}\right)}^2}{I}+\gamma I^{*}\\ +\gamma \frac{R^{*}}{R}I+\delta S^{*}+\delta \frac{V^{*}}{V}S+\beta \left(1-\eta \right)\frac{{\left(V-V^{*}\right)}^2}{V}I+\beta \left(1-\eta \right)\frac{{\left(V-V^{*}\right)}^2}{V}{I}_V\\ +\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}VI+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}V{I}_V+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}{V}^{*}{I}^{*}+\beta \left(1-\eta \right)\frac{E_V^{*}}{E_V}{V}^{*}{I}_V^{*}\\ +\gamma I_V^{*}+\gamma \frac{R_V^{*}}{R_V}{I}_V+\alpha E_V^{*}+\alpha \frac{I_V^{*}}{I_V}{E}_V+\gamma \frac{{\left(I_V-I_V^{*}\right)}^2}{I_V}+\alpha \frac{{\left(E_V-E_V^{*}\right)}^2}{E_V}.\end{array}$

It is evident that, $\frac{\text{d}{L}_f}{\text{d}t}<0$ if ${\Pi }_1<{\Pi }_2$ . However, for $S=S^{*}$ , $L=L^{*}$ , $E=E^{*}$ , $I=I^{*}$ , $R=R^{*}$ , $V=V^{*}$ , $E_V=E_V^{*}$ , $I_V=I_V^{*}$ , $R_V=R_V^{*}$ we may write, $0={\Pi }_1-{\Pi }_2$

$\Rightarrow \frac{\text{d}{L}_f}{\text{d}t}=0.$ (2.22)

We conclude that for the recommended model, the leading compact invariant set in

$\left\{\left(S^{*},L^{*},E^{*},I^{*},R^{*},V^{*},E_V^{*},I_V^{*},R_V^{*}\right)\in \Gamma :\frac{\text{d}{L}_f}{\text{d}t}=0\right\}$ (2.23)

is the endemic equilibrium point $\left\{E_{*}\right\}$ . Finally, it is evident that according to Lasalle’s invariance, if ${\Pi }_1<{\Pi }_2$ , $E_{*}$ is globally asymptotically stable in Γ.

4.9. The Second Derivative of the Lyapunov Function (LF)

Generally, the first derivative of LF assists researchers in checking the global stability of the models. However, it helps in knowing crucial information like disease sequence but not well enough to comprehend the variabilities. As a result, second derivative analysis is essential for further information, for example, curvature and sign. The second derivative, we believe, will give more details.

$\begin{array}{c}\frac{\text{d}{\stackrel{\dot{}}{L}}_f}{\text{d}t}=\frac{\text{d}}{\text{d}t}\{\left(1-\frac{S^{*}}{S}\right)\stackrel{\dot{}}{S}+\left(1-\frac{L^{*}}{L}\right)\stackrel{\dot{}}{L}+\left(1-\frac{E^{*}}{E}\right)\stackrel{\dot{}}{E}+\left(1-\frac{I^{*}}{I}\right)\stackrel{\dot{}}{I}+\left(1-\frac{R^{*}}{R}\right)\stackrel{\dot{}}{R}\\ +\left(1-\frac{V^{*}}{V}\right)\stackrel{\dot{}}{V}+\left(1-\frac{E_V^{*}}{E_V}\right){\stackrel{\dot{}}{E}}_V+\left(1-\frac{I_V^{*}}{I_V}\right){\stackrel{\dot{}}{I}}_V+\left(1-\frac{R_V^{*}}{R_V}\right){\stackrel{\dot{}}{R}}_V\}\\ ={\left(\frac{\stackrel{\dot{}}{S}}{S}\right)}^2{S}^{*}+{\left(\frac{\stackrel{\dot{}}{L}}{L}\right)}^2{L}^{*}+{\left(\frac{\stackrel{\dot{}}{E}}{E}\right)}^2{E}^{*}+{\left(\frac{\stackrel{\dot{}}{I}}{I}\right)}^2{I}^{*}+{\left(\frac{\stackrel{\dot{}}{R}}{R}\right)}^2{R}^{*}+{\left(\frac{\stackrel{\dot{}}{V}}{V}\right)}^2{V}^{*}+{\left(\frac{{\stackrel{\dot{}}{E}}_V}{E_V}\right)}^2{E}_V^{*}\\ +{\left(\frac{{\stackrel{\dot{}}{I}}_V}{I_V}\right)}^2{I}_V^{*}+{\left(\frac{{\stackrel{\dot{}}{R}}_V}{R_V}\right)}^2{R}_V^{*}+\left(1-\frac{S^{*}}{S}\right)\ddot{S}+\left(1-\frac{L^{*}}{L}\right)\ddot{L}+\left(1-\frac{E^{*}}{E}\right)\ddot{E}+\left(1-\frac{I^{*}}{I}\right)\ddot{I}\\ +\left(1-\frac{R^{*}}{R}\right)\ddot{R}+\left(1-\frac{V^{*}}{V}\right)\ddot{V}+\left(1-\frac{E_V^{*}}{E_V}\right){\ddot{E}}_V+\left(1-\frac{I_V^{*}}{I_V}\right){\ddot{I}}_V+\left(1-\frac{R_V^{*}}{R_V}\right){\ddot{R}}_V.\end{array}$ (2.24)

Here, the second derivative of Equations (1.1)-(1.9) is

$\ddot{S}=-\beta \stackrel{\dot{}}{S}\left(I+I_V\right)-\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\delta \stackrel{\dot{}}{S}-l\stackrel{\dot{}}{S}+l_d\stackrel{\dot{}}{L}$

$\ddot{L}=l\stackrel{\dot{}}{S}-\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)-\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-l_d\stackrel{\dot{}}{L}$

$\ddot{E}=\beta \stackrel{\dot{}}{S}\left(I+I_V\right)+\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)+\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)+\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha \stackrel{\dot{}}{E}$

$\ddot{I}=\alpha \stackrel{\dot{}}{E}-\gamma \stackrel{\dot{}}{I}$

$\ddot{R}=\gamma \stackrel{\dot{}}{I}$

$\ddot{V}=\delta \stackrel{\dot{}}{S}-\left(1-\eta \right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)-\left(1-\eta \right)\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)$

${\ddot{E}}_V=\left(1-q\right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)+\left(1-q\right)\beta V\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha {\stackrel{\dot{}}{E}}_V$

${\ddot{I}}_V=\alpha {\stackrel{\dot{}}{E}}_V-\gamma {\stackrel{\dot{}}{I}}_V$

${\ddot{R}}_V=\gamma {\stackrel{\dot{}}{I}}_V.$

Hence,

$\begin{array}{c}\frac{\text{d}{\stackrel{\dot{}}{L}}_f}{\text{d}t}={\left(\frac{\stackrel{\dot{}}{S}}{S}\right)}^2{S}^{*}+{\left(\frac{\stackrel{\dot{}}{L}}{L}\right)}^2{L}^{*}+{\left(\frac{\stackrel{\dot{}}{E}}{E}\right)}^2{E}^{*}+{\left(\frac{\stackrel{\dot{}}{I}}{I}\right)}^2{I}^{*}\\ +{\left(\frac{\stackrel{\dot{}}{R}}{R}\right)}^2{R}^{*}+{\left(\frac{\stackrel{\dot{}}{V}}{V}\right)}^2{V}^{*}+{\left(\frac{{\stackrel{\dot{}}{E}}_V}{E_V}\right)}^2{E}_V^{*}+{\left(\frac{{\stackrel{\dot{}}{I}}_V}{I_V}\right)}^2{I}_V^{*}\\ +{\left(\frac{{\stackrel{\dot{}}{R}}_V}{R_V}\right)}^2{R}_V^{*}+\left(1-\frac{S^{*}}{S}\right)\left\{-\beta \stackrel{\dot{}}{S}\left(I+I_V\right)-\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\delta \stackrel{\dot{}}{S}-l\stackrel{\dot{}}{S}+l_d\stackrel{\dot{}}{L}\right\}\\ +\left(1-\frac{L^{*}}{L}\right)\left\{l\stackrel{\dot{}}{S}-\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)-\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-l_d\stackrel{\dot{}}{L}\right\}\end{array}$

$\begin{array}{c}+\left(1-\frac{E^{*}}{E}\right)\{\beta \stackrel{\dot{}}{S}\left(I+I_V\right)+\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)+\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)\\ +\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha \stackrel{\dot{}}{E}\}+\left(1-\frac{I^{*}}{I}\right)\left(\alpha \stackrel{\dot{}}{E}-\gamma \stackrel{\dot{}}{I}\right)+\left(1-\frac{R^{*}}{R}\right)\gamma \stackrel{\dot{}}{I}\\ +\left(1-\frac{V^{*}}{V}\right)\left\{\delta \stackrel{\dot{}}{S}-\left(1-\eta \right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)-\left(1-\eta \right)\beta V\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)\right\}\\ +\left(1-\frac{E_V^{*}}{E_V}\right)\left\{\left(1-\eta \right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)+\left(1-\eta \right)\beta V\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha {\stackrel{\dot{}}{E}}_V\right\}\\ +\left(1-\frac{I_V^{*}}{I_V}\right)\left(\alpha {\stackrel{\dot{}}{E}}_V-\gamma {\stackrel{\dot{}}{I}}_V\right)+\left(1-\frac{R_V^{*}}{R_V}\right)\gamma {\stackrel{\dot{}}{I}}_V.\end{array}$ (2.25)

and

$\begin{array}{c}\frac{{\text{d}}^2{L}_f}{\text{d}{t}^2}=\stackrel{\dot{}}{\Pi }\left(S,L,E,I,R,V,E_V,I_V,R_V\right)\\ +\left(1-\frac{S^{*}}{S}\right)\left\{-\beta \stackrel{\dot{}}{S}\left(I+I_V\right)-\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\delta \stackrel{\dot{}}{S}-l\stackrel{\dot{}}{S}+l_d\stackrel{\dot{}}{L}\right\}\\ +\left(1-\frac{L^{*}}{L}\right)\left\{l\stackrel{\dot{}}{S}-\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)-\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-l_d\stackrel{\dot{}}{L}\right\}\\ +\left(1-\frac{E^{*}}{E}\right)\{\beta \stackrel{\dot{}}{S}\left(I+I_V\right)+\beta S\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)+\left(1-q\right)\beta \stackrel{\dot{}}{L}\left(I+I_V\right)\end{array}$

$\begin{array}{c}+\left(1-q\right)\beta L\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha \stackrel{\dot{}}{E}\}+\left(1-\frac{I^{*}}{I}\right)\left(\alpha \stackrel{\dot{}}{E}-\gamma \stackrel{\dot{}}{I}\right)+\left(1-\frac{R^{*}}{R}\right)\gamma \stackrel{\dot{}}{I}\\ +\left(1-\frac{V^{*}}{V}\right)\left\{\delta \stackrel{\dot{}}{S}-\left(1-\eta \right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)-\left(1-\eta \right)\beta V\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)\right\}\\ +\left(1-\frac{E_V^{*}}{E_V}\right)\left\{\left(1-\eta \right)\beta \stackrel{\dot{}}{V}\left(I+I_V\right)+\left(1-\eta \right)\beta V\left(\stackrel{\dot{}}{I}+{\stackrel{\dot{}}{I}}_V\right)-\alpha {\stackrel{\dot{}}{E}}_V\right\}\\ +\left(1-\frac{I_V^{*}}{I_V}\right)\left(\alpha {\stackrel{\dot{}}{E}}_V-\gamma {\stackrel{\dot{}}{I}}_V\right)+\left(1-\frac{R_V^{*}}{R_V}\right)\gamma {\stackrel{\dot{}}{I}}_V.\end{array}$ (2.26)

Finally, replacing the value of $\stackrel{\dot{}}{S},\stackrel{\dot{}}{L},\stackrel{\dot{}}{E},\stackrel{\dot{}}{I},\stackrel{\dot{}}{R},\stackrel{\dot{}}{V},{\stackrel{\dot{}}{E}}_V,{\stackrel{\dot{}}{I}}_V,{\stackrel{\dot{}}{R}}_V$ in equation (2.26), we have,

$\frac{{\text{d}}^2{L}_f}{\text{d}{t}^2}={\Sigma }_1-{\Sigma }_2,$ (2.27)

where ${\Sigma }_1$ and ${\Sigma }_2$ are the summation of all positive and negative terms.

Therefore,

$\begin{array}{l}\frac{{\text{d}}^2{L}_f}{\text{d}{t}^2}>0\text{if}{\Sigma }_1>{\Sigma }_2,\\ \frac{{\text{d}}^2{L}_f}{\text{d}{t}^2}<0\text{if}{\Sigma }_1<{\Sigma }_2,\\ \frac{{\text{d}}^2{L}_f}{\text{d}{t}^2}=0\text{if}{\Sigma }_1={\Sigma }_2.\end{array}$ (2.28)

5. Existence and Uniqueness

The present sub-section examines the existence and uniqueness of the proposed model’s solution through the concept of classical calculus.

Theorem 4: If ${\theta }_i$ and ${\bar{\theta }}_i$ are the positive constants, then

i) $\forall i\in \left\{1,2,3,\cdots ,9\right\}$

${\left|f_i\left(x_i,t\right)-f_i\left({x^{\prime }}_i,t\right)\right|}^2\le {\theta }_i{\left|x_i-{x^{\prime }}_i\right|}^2.$ (2.29)

ii) $\forall \left(x,t\right)\in {ℝ}^9\times \left(0,T\right)$

${\left|f_i\left(x_i,t\right)\right|}^2\le {\bar{\theta }}_i\left(1+{\left|x_i\right|}^2\right)\text{or}{\bar{\theta }}_i{\left|x_i\right|}^2.$ (2.30)

We may represent the current model is as follows

$\begin{array}{c}\frac{\text{d}S\left(t\right)}{\text{d}t}=-\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\delta S\left(t\right)-lS\left(t\right)+l_{d}L\left(t\right)\\ =f_1\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),\end{array}$

$\begin{array}{c}\frac{\text{d}L\left(t\right)}{\text{d}t}=lS\left(t\right)-\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-l_{d}L\left(t\right)\\ =f_2\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),\end{array}$

$\begin{array}{c}\frac{\text{d}E\left(t\right)}{\text{d}t}=\beta S\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)+\left(1-q\right)\beta L\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E\left(t\right)\\ =f_3\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),\end{array}$

$\frac{\text{d}I\left(t\right)}{\text{d}t}=\alpha E\left(t\right)-\gamma I\left(t\right)=f_4\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),$

$\frac{\text{d}R\left(t\right)}{\text{d}t}=\gamma I\left(t\right)=f_5\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),$

$\frac{\text{d}V\left(t\right)}{\text{d}t}=\delta S\left(t\right)-\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)=f_6\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),$

$\frac{\text{d}{E}_V\left(t\right)}{\text{d}t}=\left(1-\eta \right)\beta V\left(t\right)\left(I\left(t\right)+I_V\left(t\right)\right)-\alpha E_V\left(t\right)=f_7\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),$

$\frac{\text{d}{I}_V\left(t\right)}{\text{d}t}=\alpha E_V\left(t\right)-\gamma I_V\left(t\right)=f_8\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right),$

$\frac{\text{d}{R}_V\left(t\right)}{\text{d}t}=\gamma I_V\left(t\right)=f_9\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right).$

To begin, we will show that the given function $f_1\left(t,S,L,E,I,R,V,E_V,I_V,R_V\right)$ satisfies

${\left|f_1\left(S_1,t\right)-f_1\left(S_2,t\right)\right|}^2\le {\theta }_1{\left|S_1-S_2\right|}^2.$ (2.31)

Therefore,

$\begin{array}{c}{\left|f_1\left(S_1,t\right)-f_1\left(S_2,t\right)\right|}^2={\left|-\beta \left(I+I_V\right)\left(S_1-S_2\right)-\delta \left(S_1-S_2\right)-l\left(S_1-S_2\right)\right|}^2\\ ={\left|\left\{-\beta \left(I+I_V\right)-\delta -l\right\}\left(S_1-S_2\right)\right|}^2\\ \le \left\{2{\beta }^2\left({\left|I\right|}^2+{\left|I_V\right|}^2\right)+2{\delta }^2+2l^2\right\}{\left|S_1-S_2\right|}^2\\ \le \left\{2{\beta }^2\left(\underset{0\le t\le T}{\sup }{\left|I\right|}^2+\underset{0\le t\le T}{\sup }{\left|I_V\right|}^2\right)+2{\delta }^2+2l^2\right\}{\left|S_1-S_2\right|}^2\\ \le \left\{2{\beta }^2{\Vert I\left(t\right)\Vert }_{\infty }^2+2{\beta }^2{\Vert I_V\left(t\right)\Vert }_{\infty }^2+2{\delta }^2+2l^2\right\}{\left|S_1-S_2\right|}^2\\ \le {\theta }_1{\left|S_1-S_2\right|}^2,\end{array}$

where ${\theta }_1=2{\beta }^2{\Vert I\left(t\right)\Vert }_{\infty }^2+2{\beta }^2{\Vert I_V\left(t\right)\Vert }_{\infty }^2+2{\delta }^2+2l^2$ .

In the same way, the other compartments may be shown to meet the inequality mentioned above.

Secondly, we shall demonstrate that

${\left|f_1\left(S,t\right)\right|}^2\le {\bar{\theta }}_1\left(1+{\left|S\right|}^2\right).$ (2.32)

Then

$\begin{array}{c}{\left|f_1\left(S,t\right)\right|}^2={\left|-\beta S\left(I+I_V\right)-\delta S-lS+l_{d}L\right|}^2\\ ={\left|\left\{-\beta \left(I+I_V\right)-\delta -l\right\}S+l_{d}L\right|}^2\\ \le \left\{2{\beta }^2\left({\left|I\right|}^2+{\left|I_V\right|}^2\right)+2{\delta }^2+2l^2\right\}{\left|S\right|}^2+2l_d^2{\left|L\right|}^2\\ \le \left\{2{\beta }^2\left(\underset{0\le t\le T}{\sup }{\left|I\right|}^2+\underset{0\le t\le T}{\sup }{\left|I_V\right|}^2\right)+2{\delta }^2+2l^2\right\}{\left|S\right|}^2+2l_d^2\underset{0\le t\le T}{\sup }{\left|L\right|}^2\\ \le \left\{2{\beta }^2\left({\Vert I\left(t\right)\Vert }_{\infty }^2+{\Vert I_V\left(t\right)\Vert }_{\infty }^2\right)+2{\delta }^2+2l^2\right\}{\left|S\right|}^2+2l_d^2{\Vert L\left(t\right)\Vert }_{\infty }^2\\ \le {\bar{\theta }}_1\left(1+{\left|S\right|}^2\right),\end{array}$

implies that

$\frac{2{\beta }^2\left({\Vert I\left(t\right)\Vert }_{\infty }^2+{\Vert I_V\left(t\right)\Vert }_{\infty }^2\right)+2{\delta }^2+2l^2}{2l_d^2{\Vert L\left(t\right)\Vert }_{\infty }^2}<1,$

where ${\bar{\theta }}_1=2l_d^2{\Vert L\left(t\right)\Vert }_{\infty }^2$ .

In the same way, as mentioned earlier, we can also show that inequality holds for the other compartments. To summarize, our system’s solution exists and is unique, as described in [15] [36] .

6. Results and Discussions

According to [36] , the proposed model is biologically significant, and its solution is positive for all $t\ge 0$ . On top of that, specified by the entire population $N\left(t\right)$ [37] . In accordance, we prove that the proposed model is globally stable. On top of that, we deduce the point of equilibrium, which is ${\epsilon }_0=\left(0,0,0,0,0,0,0,0,0\right)$ . Also, discuss the basic and effective reproduction number and their relationship with critical vaccine proportion, which demonstrates the real scenarios of vaccine efficacy rate and basic reproduction number’s value. We know that a model’s fundamental reproduction number cannot identify whether or not it has created waves. To assist in detecting the waves, we introduced a new technique called strength number, which was produced using the next-generation matrix by calculating the second derivative of infectious classes. If the policymakers of different countries maintain the proposed model’s strategy, they can control the epidemic without making any waves. In addition, we also analyzed the LF’s first and second derivatives. The second derivative of the LF informs us of the curvature-based on its sign, whereas the first derivative tells us about the progression of the disease. Furthermore, the clarification of the proposed system of nonlinear equations is unique according to [15] and [44] .

Initially, we focused on the result of numerical simulation for the time-evolving curve about the endemic steadiness of the proposed model. The time series of susceptible, vaccinated, infected, lockdown, and recovered individuals have been portrayed in Figure 2, which presents the changing behavior of the controlling parameters δ, η, l, and q, respectively. The baseline system values are defined as the default case ( $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ ) presented in Figure 2(i). Due to lower vaccine effectiveness and vaccination rate (Figure 2(ii)), the disease incidence shows a similar tendency as in the default case; the vaccine does not work. However, Figure 2(iii) reveals that increasing the vaccination rate and vaccine effectiveness reduced the pick of infected individuals and the final epidemic size (recovered). Interestingly, the vaccinated emerges at a sporadic peak before stabilizing at equilibrium. As time passes, some people cannot maintain their health due to the vaccine’s ineffectiveness (50 percent are perfectly immune, and the remaining are non-immune).

Consequently, with the higher rate of the vaccination program and the high efficacy rate, the vaccine can control or eradicate the disease (see Figure 2(iv) and Figure 2(v)). Besides, suppose the vaccine is not available. In that case, the policymakers of different countries need to require alternate policies to control the diseases, such as lockdowns, shutdowns, states of emergency, and mask-wearing, that can help stop the spread of the COVID-19 virus. To represent the impact of lockdown irrespective of vaccination, we displayed Figure 2(vi) and Figure 2(vii) for the settings $q=0.3$ and $q=0.8$ (where $l=0.01$ and $\delta =\eta =0$ ). As shown in Figure 2(vi), the infected individuals remain unchanged with the lower lockdown maintain factor q (=0.3). However, the number of infected individuals decreased with increasing q (Figure 2(vii)). Thus, when the lockdown works properly (higher q), people are more compliant with maintaining the lockdown, which reduces the infected number of individuals. One interesting phenomenon of multi-wave characteristics observed in Figure 2(vii) may be directed by the lockdown strategy.

Figure 2. Final epidemic size (R(∞)) colored with black, suspected susceptible (S(∞)) colored with blue, lockdown (L(∞)) colored with violet, infected (I(∞)) colored with red and vaccinated (V(∞)) with green. Parameters used are (i) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.0$ , $q=0.0$ , $l=0.0$ and $\delta =0.0$ . (ii) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.5$ , $q=0.0$ , $l=0.0$ and $\delta =0.01$ . (iii) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.8$ , $q=0.0$ , $l=0.0$ and $\delta =0.05$ . (iv) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.95$ , $q=0.0$ , $l=0.0$ and $\delta =0.05$ . (v) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.95$ , $q=0.0$ , $l=0.0$ and $\delta =0.1$ . (vi) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.0$ , $q=0.3$ , $l=0.01$ and $\delta =0.0$ . and (vii) $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $\eta =0.0$ , $q=0.8$ , $l=0.01$ and $\delta =0.0$ .

Overall, Figure 3, Figure 4, and Figure 5 display the 2D heat maps of the vaccine efficacy rate (η) (x-axis) versus lockdown maintenance factor (q) (y-axis), which illustrate the final epidemic size (FES), vaccination coverage (VC) and Lockdown Individuals (LDI), respectively of the epidemic at the equilibrium point $t\to \infty$ . Moreover, panels (A-*), (B-*), and (C-*) show the results under the lockdown rate $l=0.1,0.5$ and 0.9, respectively, whereas panels (*-i), (*-ii) and (*-iii) present the results for vaccination rate $\delta =0.1,0.5$ and 0.9, respectively. Each panel in Figure 3 is partitioned into two equilibrium states: disease-free equilibrium (blue) and endemic equilibrium (deep red). As expected, reducing vaccine efficacy and lockdown maintenance factors increased the FES; the disease spread quickly. Nevertheless, the opposite tendency was found for higher η and q values.

In Figure 3, panel (A-*) for the fixed lockdown level $l=0.1$ and increasing vaccination program $\delta =0.1,0.5,0.9$ gradually decreases the FES, as expected. The shape of every heat map changed from oblique to almost parallel. However, the indisputable fact is that as the lockdown level is low, people do not stay at home; as a result, the disease spreads to the whole society even if the vaccination program is increased, a real-world phenomenon. Furthermore, in comparison to

Figure 3. Presented is the final epidemic size (FES) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show for the lockdown level (l), whereas panels (*-i), (*-ii), (*-iii) vaccination rate (δ), respectively. Parameters used are $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $l=0.1,0.5,0.9$ and $\delta =0.1,0.5,0.9$ .

Figure 4. Presented is the vaccination coverage (VC) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show the lockdown level (l), whereas panels (*-i), (*-ii), (*-iii) vaccination rate (δ), respectively. Parameters used are $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $l=0.1,0.5,0.9$ and $\delta =0.1,0.5,0.9$ .

Figure 3 with Figure 4 and Figure 5, panel (A-*), we see that in Figure 4, most of the people are covered with vaccination (high green), whereas, in Figure 5, no people are in the lockdown provision (higher gray). Thus, if the lockdown level is low, but the vaccination program is increased (no vaccine shortage) with a higher efficacy rate, the disease is eradicated from society. In this context, lockdown did not work enough. Again, for the fixed vaccination program $\delta =0.1$ and improving the lockdown level $l=0.1,0.5,0.9$ , Figure 3, in panels (A(i), B(i), C(i)), we see that FES slowly decreases, which is also practical. For example, countries like Japan and Bangladesh enhance lockdown levels from lower to higher and control the COVID-19 pandemic [45] [46] .

In comparison to Figure 3 with Figure 4 and Figure 5, panels (A(i), B(i), C(i)), we see that in Figure 4, vaccination did not work; a smaller number of the total population participates in the vaccination program. In contrast, most people are covered with the lockdown level (high violet), illustrated in Figure 5. Therefore, if a vaccination program is not high enough, i.e., a vaccine shortage exists in a society with a high efficacy rate, higher lockdown levels with high maintenance factors help policymakers control the COVID-19 pandemic. For the case in point, around 60% of the total population of the USA did not agree to

Figure 5. Presented are the lockdown Individuals (LDI) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show the lockdown level (l), whereas panels (*-i), (*-ii), (*-iii) vaccination rate (δ), respectively. Parameters used are $\beta =1.0$ , $\gamma =0.1$ , $\alpha =1/5$ , $l=0.1,0.5,0.9$ and $\delta =0.1,0.5,0.9$ .

take any vaccine and decided to maintain the lockdown policy [46] . However, according to Figure 2, a non-pharmaceutical intervention lockdown is not a permanent solution; it has economic issues. It is a one-seasonal solution because those under lockdown after a certain period become susceptible again, confirming that a highly accurate vaccination program is a permanent solution for controlling any epidemic.

Figure 3, panel (B-*) illustrates that the increasing value of lockdown level $l=0.5$ and vaccination program $\delta =0.1,0.5,0.9$ decreases the FES quicker than in Figure 3, panel (A-*). Every heat map’s size of the red region area unceasingly becomes smaller. When the lockdown level is medium with good maintenance and enhancing high efficacy vaccination programs, most people stay at home and are progressively vaccinated; the disease does not spread in society more quickly. Realistically, many countries (like Bangladesh) policymakers follow this strategy step by step to open essential offices/sectors, which also helps reduce economic loss [47] . Furthermore, panels (B-*) of Figure 4 and Figure 5 justify the illustration of Figure 3, panel (B-*), that many people participate in the vaccination program and come out from the lockdown provision day by day to fulfill their daily needs. On the other hand, we see that FES reduces significantly for the vaccination program $\delta =0.5$ and the lockdown level $l=0.1,0.5,0.9$ , Figure 3, in panels (A(ii), B(ii), C(ii)) as expected. In contrast, Figure 4 and Figure 5 elucidate that vaccination coverage (VC) reduces and the number of lockdown individuals (LDI) increases. Thus, the medium-level lockdown policy with a gradual vaccination program aids in controlling an epidemic.

Furthermore, when the lockdown level is high (shutdown, state of emergency) $l=0.9$ and maintenance, people who have locked their residences do not move elsewhere. In that case, increasing the rate of mass vaccination programs $\delta =0.1,0.5,0.9$ with a high efficacy rate considerably diminishes the FES of an epidemic (presented in Figure 3, panel C-*) compared to panel B-*, C-*, which assists policymakers in managing worst situations of the country. The government of India has overcome such a situation by applying this policy [48] . Moreover, Figure 4 revealed that in comparison with Figure 3, the VC portion reduces more than panel B-*, as expected. When any region has a higher level of lockdown with maintenance and an increasing high-efficacy vaccination program, people assume that the epidemic has died out. As stated earlier, the multi-waving phenomenon arises in Figure 2(vii). On the other hand, Figure 5 demonstrates that the number under the lockdown provision was reduced due to the increasing rate of the high-efficacy vaccination program. Again, more elevated level vaccination program $\delta =0.9$ with efficacy, and for the lockdown level $l=0.1,0.5,0.9$ , Figure 3, in panels (A(iii), B(iii), C(iii)), we see that FES turns into an endemic compared to panel (A-C(i-ii)). On top of that, the combined higher effect of the high efficacy vaccination program and increasing lockdown level eradication of the disease more quickly presented in Figure 4 and Figure 5, panels (A(iii), B(iii), C(iii)), validated the scenarios of Figure 3, panel (A-C(iii)). Realistically, if 70% of the people participate in the vaccination program, the disease automatically becomes controlled, even dying out. However, it is not possible for poor, developing, and under-developing countries, but it is likely for rich countries because of economic facts and the availability of vaccines. Therefore, the combined effect of the high efficacy vaccination program and lockdown level very shortly assists policymakers in eradicating the disease from society.

Finally, suppose we concentrate our attention diagonally. In that case, the combined effect of the high efficacy of the available vaccination program and lockdown level significantly eradicated the FES quickly, as portrayed in Figure 3. On top of that, Figure 4 and Figure 5 reveal that the combined effects work favorably. Moreover, it gives policymakers great hope in controlling the transmissible disease covid-19 from society.

Further, Figure 6 (FES), Figure 7 (VC), and Figure 8 (LDI) display the impact of vaccine shortage, vaccination rate, and lockdown effect as a form of 2D heat maps along with the vaccine efficacy rate (η) (x-axis) versus lockdown maintenance factor (q) (y-axis). Here, panels (A-*), (B-*), and (C-*) illustrate the outcomes under the settings $\left(l,\delta \right)=\left(0.0,0.9\right),\left(0.5,0.5\right)$ and $\left(0.9,0.9\right)$ correspondingly, whereas, panels (*-i), (*-ii), (*-iii), (*-iv) and (*-v) present the results

Figure 6. Presented is the final epidemic size (FES) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show the results under the lockdown rate $l=0.0,0.5,0.9$ and vaccination rate $\delta =0.9,0.5,0.9$ respectively, whereas, panels (*-i), (*-ii), (*-iii), (*-iv) and (*-v) present the results for vaccine availability rate $V_0=0.001,0.1,0.5,0.7,1.0$ , $\beta =1.0$ , $\gamma =0.1$ , and $\alpha =1/5$ respectively.

Figure 7. Presented is the vaccination Coverage (VC) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show the results under the lockdown rate $l=0.0,0.5,0.9$ and vaccination rate $\delta =0.9,0.5,0.9$ respectively, whereas, panels (*-i), (*-ii), (*-iii), (*-iv) and (*-v) present the results for vaccine availability rate $V_0=0.001,0.1,0.5,0.7,1.0$ , $\beta =1.0$ , $\gamma =0.1$ , and $\alpha =1/5$ respectively.

Figure 8. Presented are the lockdown Individuals (LDI) for vaccine effectiveness (η) and lockdown maintenance factor (q). Subpanels (A-*), (B-*), and (C-*) show the results under the lockdown rate $l=0.0,0.5,0.9$ and vaccination rate $\delta =0.9,0.5,0.9$ respectively, whereas, panels (*-i), (*-ii), (*-iii), (*-iv) and (*-v) present the results for vaccine availability rate $V_0=0.001,0.1,0.5,0.7,1.0$ , $\beta =1.0$ , $\gamma =0.1$ , and $\alpha =1/5$ respectively.

for vaccine availability rate $V_0=0.001,0.1,0.5,0.7$ and 1.0, respectively. Figure 6(A-i) and Figure 6(A-ii) represent higher FES as the value of $V_0=0.001$ , and $V_0=0.1$ when $l=0.0,\delta =0.9$ , which resulted in a significant vaccine shortage. However, the FES reduced gradually as the rate of vaccine availability increased, as displayed in Figures 6(A-iii)-(A-v) ( $V_0=0.5,0.7$ and 1.0). Furthermore, if the lockdown strategy is imposed with the vaccine program (Figure 6(B-*)), it exhibits less FES than the no-lockdown policy. Consequently, Figure 6(C-*) presents a reduced FES for total lockdown with a complete vaccine program. The country’s policymakers impose a 50% lockdown policy and continue 50% vaccination programs to control any epidemic disease when there is a vaccination shortage. Also, this strategy helps them to minimize economic loss. If the lockdown level is more stringent, $l=0.9$ , and the vaccination rate is mass $\delta =0.9$, there is no shortage of vaccines, and soon, the pandemic may be eradicated. Finally, Figure 7 and Figure 8 show that vaccinated and lockdown individuals reduced the disease when the effectiveness/efficacy and acceptance rate were high. In practice, if 70% of the people participate in the vaccination program, the disease automatically becomes controlled, reflecting our current results.

7. Conclusion

In response to the dangerous circumstances provided by the COVID-19 pandemic, we have introduced the $SLEIRV{E}_V{I}_V{R}_V$ epidemic model, which takes into account the scarcity of vaccinations. Furthermore, the suggested model holds great biological significance, and its solution is positive for all values of $t\ge 0$ . Furthermore, the total population $N\left(t\right)$ is explicitly indicated [37] . Furthermore, analyzes the concept of equilibrium, denoted as R₀ and R_e, and explore the correlation between critical vaccination proportions, strength number, and the first and second derivatives of LF. Ultimately, we have demonstrated that the suggested model exhibits global stability. The finite difference approach has been utilized for numerical simulation. Furthermore, when illustrating real-world situations, the heat maps are depicted even though they are entirely based on assumptions. Pragmatically, this examination is crucial for enhancing global real-life circumstances. Therefore, our comprehensive analysis indicates that the synergistic impact of implementing non-pharmaceutical interventions, including lockout measures and a widespread and efficient vaccination program, is very successful in eliminating the illness from both society and minimizing economic losses. As mentioned above, the findings illustrate the implementation of strict lockdown measures and a very efficient vaccination campaign that has not experienced any supply shortages. Following the supplementary expansion, we will utilize this model to gather insights on effectively managing the lockdown measures and vaccination program as governments strive to control the COVID-19 pandemic.

Acknowledgments

Anonymous reviewers gave substantially meaningful comments that improved our manuscript and made it more robust. We want to express our gratitude to them.

Conflict of Interest

The authors declare that they have no conflict of interest.

Fundings

This study has no funding.

Data

This article has no additional data.

Conflicts of Interest

The authors declare that they have no conflict of interest.

References

[1]	(2020) COVID-19 Coronavirus Pandemic. https://www.worldometers.info/coronavirus/#repro
[2]	Li, R.Y.M., Yue, X.G. and Crabbe, M.J.C. (2021) COVID-19 in Wuhan, China: Pressing Realities and City Management. Frontiers in Public Health, 8, Article 596913. https://doi.org/10.3389/fpubh.2020.596913
[3]	Anderson, R.M. and May, R.M. (1979) Population Biology of Infectious Disease: Part I. Nature, 280, 361-367. https://doi.org/10.1038/280361a0
[4]	Ross, R. (1916) An Application of the Theory of Probabilities to the Study of a Priori Pathometry—Part I. Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 92, 204-230. https://doi.org/10.1098/rspa.1916.0007
[5]	Kermack, W.O. and McKendrick, A.G. (1927) A Contribution to the Mathematical Theory of Epidemics, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character, 115, 700-721. https://doi.org/10.1098/rspa.1927.0118
[6]	Bailey Norman, T.J. (1975) The Mathematical Theory of Infectious Diseases and Its Applications. 2nd Edition, Griffin, London.
[7]	Anderson, R.M. (1982) Population Dynamics of Infectious Diseases: Theory and Applications. Chapman and Hall, London.
[8]	Brauer, F., Castillo-Chavez, C. and Castillo-Chavez, C. (2012) Mathematical Models in Population Biology and Epidemiology. Springer, New York. https://doi.org/10.1007/978-1-4614-1686-9
[9]	Martcheva, M. (2015) An Introduction to Mathematical Epidemiology. Springer, New York. https://doi.org/10.1007/978-1-4899-7612-3
[10]	Michael, Y.L. (2017) An Introduction to Mathematical Modeling of Infectious Diseases. Springer, New York.
[11]	Kabir, K.M.A. and Tanimoto, J. (2020) Evolutionary Game Theory Modeling to Represent the Behavioral Dynamics of Economic Shutdowns and Shield Immunity in the COVID-19 Pandemic. Royal Society Open Science, 7, Article ID: 201095. https://doi.org/10.1098/rsos.201095
[12]	Kabir, K.M.A., Risa, T. and Tanimoto, J. (2021) Prosocial Behavior of Wearing a Mask during an Epidemic: An Evolutionary Explanation. Scientific Reports, 11, Article No. 12621. https://doi.org/10.1038/s41598-021-92094-2
[13]	Alam, M., Kabir, K.M.A. and Tanimoto, J. (2020) Based on Mathematical Epidemiology and Evolutionary Game Theory, Which Is More Effective: Quarantine or Isolation Policy? Journal of Statistical Mechanics: Theory and Experiment, 2020, Article ID: 033502. https://doi.org/10.1088/1742-5468/ab75ea
[14]	Higazy, M., Allehiany, F.M. and Mahmoud, E.E. (2021) Numerical Study of Fractional-Order COVID-19 Pandemic Transmission Model in the Context of ABO Blood Group. Results in Physics, 22, Article ID: 103852. https://doi.org/10.1016/j.rinp.2021.103852
[15]	Ullah, M.S., Higazy, M. and Kabir, K.M.A. (2021) Modeling the Epidemic Control Measures in Overcoming COVID-19 Outbreaks: A Fractional-Order Derivative Approach. Chaos, Solitons & Fractals, 155, Article ID: 111636. https://doi.org/10.1016/j.chaos.2021.111636
[16]	Higazy, M. and Alyami, M.A. (2020) New Caputo-Fabrizio Fractional-Order SEIASqEqHR Model for COVID-19 Epidemic Transmission with Genetic Algorithm-Based Control Strategy. Alexandria Engineering Journal, 59, 4719-4736.
[17]	Das, P., Nadim, S.S., Das, S., et al. (2021) Dynamics of COVID-19 Transmission with Comorbidity: A Data Driven Modelling Based Approach. Nonlinear Dynamics, 106, 1197-1211. https://doi.org/10.1007/s11071-021-06324-3
[18]	Islam, M.S., Ira, J.I., Kabir, K.M.A. and Kamrujjaman, M. (2020) Effect of Lockdown and Isolation to Suppress the COVID-19 in Bangladesh: An Epidemic Compartments Model. Journal of Mathematical Analysis and Applications, 4, 83-93. https://doi.org/10.26855/jamc.2020.09.004
[19]	Chowdhury, A., Kabir, K.M.A. and Tanimoto, J. (2020) How to Quarantine and Social-Distancing Policy Can Suppress the Outbreak of Novel Coronavirus in Developing or under Poverty Level Countries: A Mathematical and Statistical Analysis. https://doi.org/10.21203/rs.3.rs-20294/v1
[20]	Das, P., Upadhyay, R.K., Misra, A.K., et al. (2021) Mathematical Model of COVID-19 with Comorbidity and Controlling Using Non-Pharmaceutical Interventions and Vaccination. Nonlinear Dynamics, 106, 1213-1227. https://doi.org/10.1007/s11071-021-06517-w
[21]	Anderson, R.M. and May, R.M. (1985) Vaccination and Herd Immunity to Infectious Diseases. Nature, 318, 323-329. https://doi.org/10.1038/318323a0
[22]	Struchiner, C.J., Halloran, M.E. and Spielman, A. (1989) Modeling Malaria Vaccines I: New Uses for Old Ideas. Mathematical Biosciences, 94, 837-113. https://doi.org/10.1016/0025-5564(89)90073-4
[23]	Galvani, A.P., Reluga, T.C. and Chapman, G.B. (2007) Long-Standing Influenza Vaccination Policy Is in Accord with Individual Self-Interest But Not with the Utilitarian Optimum. Proceedings of the National Academy of Sciences of the United States of America, 104, 5692-5697. https://doi.org/10.1073/pnas.0606774104
[24]	Anderson, R.M. and May, R.M. (1992) Infectious Diseases of Humans. Oxford Academic, Oxford. https://doi.org/10.1093/oso/9780198545996.001.0001
[25]	Coburn, B.J., Wagner, B.G. and Blower, S. (2009) Modeling Influenza Epidemics and Pandemics: Insights into the Future of Swine Flu (H1N1). BMC Medicine, 7, Article No. 30. https://doi.org/10.1186/1741-7015-7-30
[26]	Longini, I.M., Halloran, J.M.E., Nizam, A. and Yang, Y. (2004) Containing Pandemic Influenza with Antiviral Agents. American Journal of Epidemiology, 159, 623-633. https://doi.org/10.1093/aje/kwh092
[27]	Viboud, C., Boelle, P.Y., Carrat, F., Valleron, A.J. and Flahault, A. (2000) Prediction of the Spread of Influenza Epidemics by the Method of Analogues. American Journal of Epidemiology, 158, 996-1006. https://doi.org/10.1093/aje/kwg239
[28]	Viboud, C., Boelle, P.Y., Pakdaman, K., Carrat, F., Valleron, A.J. and Flahault, A. (2004) Influenza Epidemics in the United States, France, and Australia, 1972-1997. Emerging Infectious Diseases, 10, 32-39. https://doi.org/10.3201/eid1001.020705
[29]	Choi, Y., Kim, J.S., Kim, J.E., Choi, H. and Lee, C.H. (2021) Vaccination Prioritization Strategies for COVID-19 in Korea: A Mathematical Modeling Approach. International Journal of Environmental Research and Public Health, 18, Article 4240. https://doi.org/10.3390/ijerph18084240
[30]	MacIntyre, C.R., Costantino, V. and Trent, M. (2020) Modelling of COVID-19 Vaccination Strategies and Herd Immunity, in Scenarios of Limited and Full Vaccine Supply in NSW, Australia (Under Review). https://doi.org/10.1101/2020.12.15.20248278
[31]	Chen, F. (2018) Voluntary Vaccinations and Vaccine Shortages: A Theoretical Analysis. Journal of Theoretical Biology, 446, 19-32. https://doi.org/10.1016/j.jtbi.2018.02.032
[32]	Li, Z., Xu, J., Xu, J., Tan, T. and Zhang, C. (2020) Current Situation, Causes, and Countermeasures to NIP Vaccine Shortages in Guangzhou, China. Human Vaccines & Immunotherapeutics, 16, 76-79. https://doi.org/10.1080/21645515.2019.1644883
[33]	Kahwati, L.C., Elter, J.R., Straits-Tröster, K.A. Kinsinger, L.S. and Davey, V.J. (2007) The Impact of the 2004-2005 Influenza Vaccine Shortage in the Veterans Health Administration. Journal of General Internal Medicine, 22, 1132-1138. https://doi.org/10.1007/s11606-007-0249-6
[34]	Fairbrother, G., Broder, K., Staat, M.A., Schwartz, B., Heubi, C., Hiratzka, S., Walker, F.J. and Morrow, A.L. (2007) Pediatricians’ Adherence to Pneumococcal Conjugate Vaccine Shortage Recommendations in 2 National Shortages. Pediatrics, 120, e401-e409. https://doi.org/10.1542/peds.2007-0359
[35]	Allison, T.C., Wells, M.S.K., Seib, K., Kudis, A., Hannan, C., Orenstein, W.A., Whitney, E.A.S., Hinman, A.R., Buehler, J.W., Omer, S.B. and Berkelman, R.L. (2012) Lessons Learned from the 2007 to 2009 Haemophilus influenzae Type B Vaccine Shortage. Journal of Public Health Management and Practice, 18, E9-E16. https://doi.org/10.1097/PHH.0b013e31821dce27
[36]	Walter, W. (1998) Ordinary Differential Equations. Springer, New York. https://doi.org/10.1007/978-1-4612-0601-9
[37]	Nabi, K.N., Abboubakar, H. and Kumar, P. (2020) Forecasting of Covid-19 Pandemic: From Integer Derivatives to Fractional Derivatives. Chaos, Solitons & Fractals, 141, Article ID: 110283. https://doi.org/10.1016/j.chaos.2020.110283
[38]	Kamrujjaman, M., Mahmud, M.S., Ahmed, S., Qayum, M.O., Alam, M.M., Hassan, M.N., Islam, M.R., Nipa, K.F. and Bulut, U. (2021) SARS-CoV-2 and Rohingya Refugee Camp, Bangladesh: Uncertainty and How the Government Took over the Situation. Biology, 10, Article 124. https://doi.org/10.3390/biology
[39]	Mahmud, M.S., Kamrujjaman, M., Adan, M.M.I.Y., et al. (2022) Vaccine Efficacy and SARS-CoV-2 Control in California and US during the Session 2020-e2026: A Modeling Study. Infectious Disease Modelling, 7, 62-81. https://doi.org/10.1016/j.idm.2021.11.002
[40]	Kamrujjaman, M., Saha, P., Islam, M.S., et al. (2022) Dynamics of SEIR Model: A Case Study of COVID-19 in Italy. Results in Control and Optimization, 7, Article ID: 100119. https://doi.org/10.1016/j.rico.2022.100119
[41]	Islam, M.R., Gray, M.J. and Peace, A. (2021) Identifying the Dominant Transmission Pathway in a Multi-Stage Infection Model of the Emerging Fungal Pathogen Batrachochytrium salamandrivorans on the Eastern Newt. In: Teboh-Ewungkem, M.I. and Ngwa, G.A., Eds., Infectious Diseases and Our Planet, Springer, Cham, 193-216. https://doi.org/10.1007/978-3-030-50826-5_7
[42]	Yorke, J.A. (1967) Invariance for Ordinary Differential Equations. Mathematical Systems Theory, 1, 353-372. https://doi.org/10.1007/BF01695169
[43]	Diekmann, O. and Heesterbeek, J.A.P. (2000) Mathematical Models in Population Biology and Epidemiology. Springer, New York.
[44]	Atangana, A. (2021) Mathematical Model of Survival of Fractional Calculus, Critics and Their Impact: How Singular Is Our World? Advances in Difference Equations, 2021, Article No. 403. https://doi.org/10.1186/s13662-021-03494-7
[45]	Wikipedia: COVID-19 Pandemic in Bangladesh.
[46]	AAAS (2020) Just 50% of Americans Plan to Get a COVID-19 Vaccine. Here’s How to Win over the Rest.
[47]	Pew Research Center (2020) US Public Now Divided over Whether to Get COVID-19 Vaccine.
[48]	BBC News (2021) COVID: India Sees the World’s Highest Daily Cases Amid Oxygen Shortage.