On a Novel Dynamics of a SIVR Model Using a Laplace Adomian Decomposition Based on a Vaccination Strategy

Abstract

In this paper, we introduce a SIVR model using the Laplace Adomian decomposition. This model focuses on a new trend in mathematical epidemiology dedicated to studying the characteristics of vaccination of infected communities. We analyze the epidemiological parameters using equilibrium stability and numerical analysis techniques. New mathematical strategies are also applied to establish our epidemic model, which is a pandemic model as well. In addition, we mathematically establish the chance for the next wave of any pandemic disease and show that a consistent vaccination strategy could control it. Our proposal is the first model introducing a vaccination strategy to actively infected cases. We are sure this work will serve as the basis for future research on COVID-19 and pandemic diseases since our study also considers the vaccinated population.

1. Introduction, Motivation, and Objectives

1.1. Introductory Aspects

Mathematical models of epidemic describe the number of infected, death, or recovered cases. These models provide information to individuals in organizations and healthcare to understand how infectious diseases grow. Epidemic models are classified into two types: (i) stochastic and (ii) deterministic. In the stochastic case, one uses random variables and arrives at the solutions based on probability distributions [1,2]. In the determinist case, one employs differential equations to arrive at the solutions with the change rate of each compartment depending on the unit time, which is often measured in days [3,4]. Recently, many mathematicians have started to study fractional differential equations (FDE) [5] because the journey of fractional calculus is wonderful to read both in the aspects of history and the theory of calculus, beginning with the contributions of Riemann–Liouville, then of Caputo, and later of Atangana–Balenau in the Caputo (ABC) sense [6,7].

One can estimate the basic reproduction number of a disease and attain the global stability of the solutions of the system by proving the local stability of equilibrium points. Disease-free and disease-dependent equilibrium points, feasibility region analyses, as well as the existence of positive solutions can also be used. Numerical solutions may be brought up by applying the Laplace Adomian decomposition method (LADM). If the problem is of a linear and ordinary differential equation (ODE), we are not required to seek numerical results since analytical results are available. Nevertheless, for nonlinear problems, only very few cases possess analytical results and most of them do not possess such results. We utilize numerical approaches as the Runge–Kutta method (RKM) in our research.

Note that the RKM does not allow us to obtain solutions of an FDE. The LADM permits us to obtain semi-analytical solutions. In [8], the LADM was utilized to study a fractional model of the human immunodeficiency virus (HIV) and measles [9]. In the initial period of fractional calculus, some of the problems for FDEs involved Riemann–Liouville fractional derivatives only. Today, mathematicians are interested in this type of derivative [10]. The stability of such FDEs was analyzed in [11]. Nonetheless, if such nonlinear problems are of fractional order, it is hard to solve them step by step. Several methods are available in the literature, such as FDEs, to solve these problems. Some of them are the differential transform method, the homotopy perturbation method, and the LADM. We prefer the LADM because it uses Adomian polynomials, the basic idea of Laplace transformations, and their inverse to obtain the solutions. Previous studies that motivated us to conduct the present research are summarized below. When one wants to model any biological or medical problem using fractional calculus, one must be careful in estimating the corresponding parameters [12].

1.2. Motivation and State of the Art

The motivation for framing our proposal is as follows. Considering the susceptible, infectious, and recovered (SIR) model, one can assume a new compartment to analyze what does happen when someone is vaccinated during the period of infection; that is, by considering a susceptible, infectious, recovered, and vaccinated (SIRV) model. One can ask: Does such a disease model exist? The response is yes. We consider some people with type I or type II diabetes who are only in the category “I”. However, doctors prescribe to these diabetes patients several vaccines such as those for influenza, pneumococcal infection, Tdap, hepatitis B, and Zoster to prevent seasonal flu, pneumococcal diseases, tetanus, diphtheria, whooping cough, hepatitis B, and shingles, respectively. In some cases, vaccines may be used to treat active viral infections. The concept behind this is to improve immunity by providing vaccinations against the virus that causes the symptoms and the infectious disease.

Let us consider rabies, a dangerous neurological disorder, as an example. Often, rabies is a severe infection that spreads through the saliva of some animals, especially mammals infected with rabies. It takes up to two or three weeks for mammal or animal rabies to infect our central nervous system. Medically speaking, one is suspected of having rabies when bitten by such animals. Earlier, one may start with symptoms. By the symptoms and lab report, the infection can be confirmed. This period of two or three weeks of infectious time is enough to vaccinate an infected individual. The vaccines are given to such an actively infectious person to boost the immune response and prevent the virus from penetrating the nerve tissue. Thus, vaccinations during the actively infectious stage of rabies prevent the dangerous neurological effects of the disease. First, one recovers from the affected bite and then from the infections. The same approach is also applied to treat Ebola, which is known to be the fastest-spreading and deadliest virus on the globe.

The diseases spread by animals, such as bats and monkeys, are responsible for the death of 80% of infected individuals within two weeks. By applying newly found vaccines, victims can strengthen their immune systems and control the active virus to save their life. Observe that one can present an SIVR model where the susceptible individual becomes infected. Furthermore, for a situation where an infected person undergoes vaccination therapies, what happens mathematically in the number of infected cases and what happens with the recovered cases as time increases can be studied using a SIRV model.

Note that when somebody has diabetes, they undergo to many vaccines to prevent them from any new infections. When one is infected with rabies, vaccines prevent worsening severity of the disease, which results in damage to the central nervous system. However, both of these diseases are different. We are motivated by these diseases and wish to create a new SIVR model. Thus, our aim is to frame the model in which the vaccinations are given to actively infected persons to recover quickly and prevent them from acquiring severe infections [13].

Some SVIR models have been proposed [14,15]. Applications of FDEs were predominantly equipped with epidemic models, and their analysis [16,17,18,19]. Some epidemic models, such as those for severe acute respiratory syndrome [20], HIV [21,22,23], tumor/cancer [24,25], tuberculosis [26], hepatitis B [27], and COVID-19 [28,29,30,31,32] were studied using fractional derivatives. Moreover, some physical systems such as the visco-elastic damper model [33] were analyzed with FDEs. Other than this, the theory of periodic reactions [34] and the advanced analysis of SIVR models [35,36] have also paved a path on which our time-dependent disease model can proceed.

1.3. Objective and Description of Sections

According to our bibliographical review, to the best of our knowledge, there are no SVIR models that describe a pandemic based on a new strategy called the constant and consistent application of vaccination. After the pandemic gradually decreases, a consistent depletion vaccination strategy can be analyzed to mathematically state how it could help to show depletion in the number of infections in society. Therefore, the main objectives of our investigation are as follows: (i) to formulate a new pandemic model with SIVR cases; (ii) to study strategies related to constant/consistent application and consistent depletion vaccinations; and (iii) to analyze how our model could help to mathematically show depletion in the number of infections in society.

This paper is organized as follows. The presentation of the SIVR model is in Section 2. This section also introduces the necessary preliminary knowledge. In Section 3, we provide equilibrium values and stability, as well as a feasibility region analysis. The positive solutions are outlined in Section 4, with a sensitivity analysis also being presented in this section. Section 5 reports our numerical results. In Section 6, the findings obtained are discussed. Finally, conclusions of the present work are provided.

2. Background

2.1. SIVR Epidemic Model

Consider the following equations for building up the desired mathematical model:

$$\begin{array}{c}\hfill \begin{array}{ccc}{S}^{\prime }\left(t\right)\hfill & =\hfill & -aS\left(t\right)I\left(t\right),I^{\prime }\left(t\right)=aS\left(t\right)I\left(t\right)-(b+r)I\left(t\right),\hfill \\ V^{\prime }\left(t\right)\hfill & =\hfill & (b-c-r)I\left(t\right)-rV\left(t\right),R^{\prime }\left(t\right)=cI\left(t\right)-rR\left(t\right),\hfill \end{array}\end{array}$$

(1)

where $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$, and $V\left(t\right)$ in the SIRV model stated in (1) permit us to obtain the number of susceptible, infected, recovered, and vaccinated people at time t, respectively. Additionally, a, b, c, and r are the rates at which the susceptible become infected, the infected become vaccinated, the vaccinated recovered, and the infected are recovered, respectively. Note that

$$N\left(t\right)=S\left(t\right)+I\left(t\right)+V\left(t\right)+R\left(t\right),$$

(2)

where $N\left(t\right)$ stated in (2) is the total population at time t. To frame the model as seen in Figure 1, we assume several values given in

Table 1. Notations/symbols employed in the SIVR model and values for a, b, c, and r as in [3,37].

Notations/Symbols	Definition
S	Susceptible population
I	Infected population
V	Vaccinated population
R	Recovered population
t	Time instant
$S_0$	Initial susceptible population, $S_0=400$ (fixed)
$I_0$	Initial infected population, $I_0=300$ (fixed)
$V_0$	Initial vaccinated population, $V_0=200$ (fixed)
$R_0$	Initial recovered population, $R_0=100$ (fixed)
$S_{\mathrm{critical}}$	Threshold point or epidemic critical community size
${\mathcal{R}}_0$	Basic reproduction number
a	Rate at which those susceptible become infected, $a=0.0004$ (fixed)
b	Rate at which those infected become vaccinated, $b=0.03$ (fixed)
c	Rate at which those infected recover via vaccination, $c=0.01$ (fixed)
r	Rate at which those infected recover naturally, $r=0.02$ (fixed)

. The values in this table are used for illustrative purposes only and are not intended to represent true values for any specific disease or pandemic. The values of a, b, c, and r will depend on the specific characteristics of the disease or pandemic, as well as on the population being studied. These values are subject to change over time as more information becomes available about the disease or pandemic. Choice of the population size will affect the numerical results of the model. In this study, a population size of 100,000 was chosen for convenience, but the model can be applied to any population size with appropriate adjustments to the parameter values. Specifically, we have a population size to 100,000 and values of 0.004 × 100,000 = 400, 0.003 × 100,000 = 300, 0.002 × 100,000 = 200, and 0.001 × 100,000 = 100 for S, I, V, and R, respectively; that is, we take a very small part of $0.01\%$ of the total population and reduce it still further for transmission of S to I, I to V, recovery by vaccination (V to R), and natural recovery (I to R) to obtain a, b, c, and r, whose values are given in Table 1.

In general, it is recommended to use data-driven parameter values for epidemiological models whenever possible. These values can be obtained from historical data or from ongoing surveillance efforts. Using data-driven parameter values can improve the accuracy of the model and its predictions.

The theory of fractional calculus has recently influenced more researchers and academics. Diverse real-life phenomena can be easily modeled with the help of fractional calculus. Fractional calculus studies and establishes the derivatives formed by various orders between 0 and 1, often denoted as $\alpha$; that is, $0\le \alpha \le 1$. At $\alpha =0$, we obtain the algebraic nondifferential equation. At $\alpha =1$, we arrive at the first-order differential equations. For example, the half derivative is given by fractional-order $\alpha =0.5$.

In fractional calculus, we can also reach first-order derivatives. Because this calculus generalizes ordinary differential calculus, it is impossible to obtain half derivatives or any noninteger order derivatives in ordinary differential calculus. Where the epidemic modeling is concerned, the changes observed in every compartment are low for a unit time. Thus, instead of an ODE, we prefer an FDE. Furthermore, we must draw attention to one important aspect: the ordinary differential derivatives are local.

In contrast, the fractional-order derivatives are global since such derivatives are used for hiking the region of stability of any physical system. When considering the fractional derivative without the singular kernel, one must remember Caputo and Fabrizio [6]. The special feature of the ABC fractional derivative [7] is that it is the new fractional derivative with a nonlocal and nonsingular kernel. Thus, some necessary results in FDEs are presented here.

2.2. Preliminaries on Fractional Differential Equations

Next, we introduce some necessary definitions and explain our motivation for utilizing fractional derivatives. Fractional order operators are now widely used to model numerous mathematical and physical problems in different areas, such as physics, chemistry, and biology. Mathematical models from these areas have been studied and solved by employing different fractional order operators. Singular kernels seem to be the problem in the integrands of Riemann–Liouville and Caputo type fractional operators.

Such kernels are not smooth at modeling problems and coming up with finer solutions. To overcome this, Atangana and Baleanu introduced a new advanced operator based upon the Mittag-Leffler type of kernel [7]. The kernel given was nonlocal and nonsingular. The important advantage of this operator is that it is very helpful in a modeling biological dynamical systems since the Atangana–Baleanu fractional order derivative eliminates the difficulty of modeling any biological or physical problems with singularity. Various biological systems have been modeled and analyzed using the ABC operator.

Definition 1

([7]). The Mitag–Leffler function is the solution of the fractional ODE $D^{\alpha }y=ay$, for $0<\alpha <1$, where a generalized Mittag-Leffler function ${\mathcal{E}}_{\alpha }(-t^{\alpha })={\sum }_{k=0}^{\infty }{(-t)}^{\alpha k}/\Gamma (\alpha k+1)$ is considered as a nonlocal function, with Γ representing the standard gamma function.

Definition 2

([7]). The ABC fractional derivative of function f supported on $[0,t]$ is stated as

$${\displaystyle {}^{\mathrm{ABC}}{D}_{[0,t]}^{\alpha }\left\{f\left(t\right)\right\}=\frac{F\left(\alpha \right)}{(1-\alpha )}{\int }_0^t{f}^{\prime }\left(\mu \right){\mathcal{E}}_{\alpha }(-\frac{\alpha {(t-\mu )}^{\alpha }}{(1-\alpha )})\mathrm{d}\mu ,}$$

where $0<\alpha \le 1$ is the fractional order, $F\left(\alpha \right)=(1-\alpha )/\Gamma (1-\alpha )$ is a normalization function that satisfies $F\left(0\right)=F\left(1\right)=1$, as in the Caputo–Fabrizio case [6], and ${\mathcal{E}}_{\alpha }$ is the Mittag-Leffler function of order α defined as

$${\mathcal{E}}_{\alpha }\left(z\right)=\sum _{k=0}^{+\infty }\frac{z^k}{\Gamma (k\alpha +1)}.$$

The ABC derivative can be transformed to Caputo by using

$${\mathcal{E}}_{\alpha }(-\frac{\alpha {(t-\mu )}^{\alpha }}{(1-\alpha )})$$

to ${(t-\mu )}^{-\alpha }$, and

$$\frac{F\left(\alpha \right)}{(1-\alpha )}$$

to $1/\Gamma (1-\alpha )$.

Definition 3

([7]). The Laplace transform of the ABC derivative of Definition 2 is given by

$$\mathcal{L}\left\{{}^{\mathrm{ABC}}{D}_{[0,t]}^{\alpha }\left\{f\left(t\right)\right\}\right\}\left(u\right)=\frac{F\left(\alpha \right)}{(1-\alpha )}\left({\displaystyle \frac{u^{\alpha }\mathcal{L}\left\{f\left(t\right)\right\}\left(u\right)-u^{\alpha -1}}{u^{\alpha }+\alpha /\left(1-\alpha \right)}}\right),$$

where $\mathcal{L}\left\{f\left(t\right)\right\}$ is the Laplace transform of $f\left(t\right)$.

Definition 4

([38]). The Atangana–Baleanu integral of f on $[0,t]$ with $\alpha >0$ is expressed as

$${\displaystyle {}^{\mathrm{ABC}}{J}_{[0,t]}^{\alpha }\left\{f\left(t\right)\right\}=\frac{(1-\alpha )}{F\left(\alpha \right)}f\left(t\right)+\frac{\alpha }{F\left(\alpha \right)\Gamma \left(\alpha \right)}{\int }_0^{t}f\left(\mu \right){(t-\mu )}^{\alpha -1}\mathrm{d}\mu .}$$

Works such as [39,40,41,42] serve as a motivation and guidance for us to construct a fractional ordered system of differential equations. We recommend readers to refer to [3,43] for a better understanding of the various types of epidemic models that can be estimated using FDEs.

The system described in (1) considers the variables $S\left(t\right)$, $I\left(t\right)$, $R\left(t\right)$, and $V\left(t\right)$. This system can be reformulated by means of an FDE system stated as

$$\begin{array}{c}\hfill \begin{array}{ccc}{}^{\mathrm{ABC}}{D}^{{\alpha }_1}S\left(t\right)\hfill & =\hfill & -aS\left(t\right)I\left(t\right),\hfill \\ {}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t\right)\hfill & =\hfill & aS\left(t\right)I\left(t\right)-(b+r)I\left(t\right),\hfill \\ {}^{\mathrm{ABC}}{D}^{{\alpha }_3}V\left(t\right)\hfill & =\hfill & (b-c-r)I\left(t\right)-rV\left(t\right),\hfill \\ {}^{\mathrm{ABC}}{D}^{{\alpha }_4}R\left(t\right)\hfill & =\hfill & cI\left(t\right)-rR\left(t\right).\hfill \end{array}\end{array}$$

(3)

The total population size, at $t=0$, is given by $N\left(0\right)=N_0=S_0+I_0+V_0+R_0$. Now, we have to establish a few things. How do the models formulated in (1) and (3) relate real-life phenomena? Why is this model more relevant to any existing diseases? For establishing this, we can refer to the model stated in (1), where we consider an ODE for the derivatives; that is, the by-products that depend on the rate of changes concerning time are fully changed from 0-th to 1-st order, which frequently requires at least a unit time with a magnitude of days. However, similar to that presented in (3), the model is converted to the fractional ordered system, allowing for the by-products to indicate that derivatives can exist not only in zero and first orders but also in fractional orders; that is, very minute changes between zero and first orders. For an example, consider $D^{{\alpha }_1}$ used in the above equation to provide the half-order derivative when ${\alpha }_1=0.5$, which is similar for other orders ${\alpha }_2,{\alpha }_3,{\alpha }_4$. In contrast, one can adopt different orders to each compartment $S,I,V,R$ by taking different values of ${\alpha }_i$ for $i\in \{1,2,3,4\}$. Similarly, one can obtain fourth order changes in the model by assigning ${\alpha }_i=0.25$, for $i\in \{1,2,3,4\}$.

3. Mathematical Analysis

3.1. Justification

In this paper, we present two methods for deriving the basic reproduction number, ${\mathcal{R}}_0$, for FDEs. One method is the next generation matrix, and the other is the traditional method. By using both methods, we establish that they are equally effective. Although the basic reproduction number is calculated by considering FDEs as normal algebraic equations, where the order of derivatives is equal to zero, these methods are applicable to both ODEs and FDEs. However, FDEs offer more flexibility in finding solutions at various orders, where $0<{\alpha }_i\le 1$, for $i\in \{1,2,3,4\}$. It is important to note that for any fractional order less than one, any values of $\mathcal{R}0$ and $S_{\mathrm{critical}}$ that do not exceed the unit order will be bound within the values of an integer order. For this reason, we have defined the analysis for an ODE. The fractional order is better because the epidemic will not grow suddenly but will rather grow fractionally over time. Even for orders less than one, there is a small change in the growth rate. This small change over a period can only be predicted using FDEs.

3.2. Estimation of the Basic Reproduction Number with a First Computation Method

Assume that the class of secondary susceptible individuals can be formed from both infected and vaccinated classes. Since we are considering a SIVR model and not a SIVRS model, we are not assuming that the recovered class will be susceptible again. Thus, we employ next-generation matrix techniques for large domains.

Let us consider

$$\begin{array}{ccc}\hfill {}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t\right)& =& aS\left(t\right)I\left(t\right)-(b+r)I\left(t\right),\hfill \\ \\ \hfill {}^{\mathrm{ABC}}{D}^{{\alpha }_3}V\left(t\right)& =& (b-c-r)I\left(t\right)-rV\left(t\right).\hfill \end{array}$$

(4)

The Jacobian matrix of the expressions formulated in (4) is given by

$$\begin{array}{c}\hfill G=\left(\begin{array}{cc}aS\left(0\right)-(b+r)& 0\\ b-c-r& -r\end{array}\right).\end{array}$$

(5)

Now, the matrix G stated in (5) can be decomposed as $G=G_1+G_2$, where

$$G_1=\left(\begin{array}{cc}aS\left(0\right)& 0\\ 0& 0\end{array}\right)$$

and

$$G_2=\left(\begin{array}{cc}-(b+r)& 0\\ b-c-r& -r\end{array}\right).$$

Let us calculate Y from $Y=-G_2$. Then, we have

$$\begin{array}{c}\hfill Y=\left(\begin{array}{cc}b+r& 0\\ r+c-b& r\end{array}\right).\end{array}$$

Hence, we find that

$$\begin{array}{c}\hfill \mathrm{adj}\left(Y\right)=\left(\begin{array}{cc}r& 0\\ b-c-r& b+r\end{array}\right),\left|Y\right|=(b+r)r,Y^{-1}=\left(\begin{array}{cc}\frac{1}{b+r}& 0\\ \frac{b-c-r}{r(b+r)}& \frac{1}{r}\end{array}\right).\end{array}$$

Therefore, we obtain

$$\begin{array}{c}\hfill G_1{Y}^{-1}=\left(\begin{array}{cc}\frac{aS\left(0\right)}{b+r}& 0\\ 0& 0\end{array}\right).\end{array}$$

Note that the basic reproduction number is calculated from ${\mathcal{R}}_0=\rho \left(G_1{Y}^{-1}\right)$, where $\rho \left(G_1{Y}^{-1}\right)$ is the spectral radius of the matrix $G_1{Y}^{-1}$ and is given by $max\left(\right|{\lambda }_i\left|\right)$, with ${\lambda }_i$ representing the eigenvalues of the matrix $G_1{Y}^{-1}$. Thus, ${\mathcal{R}}_0=S_{0}a/\left(b+r\right)$ since $S\left(0\right)=S_0$.

3.3. Estimation of the Basic Reproduction Number with a Second Computation Method

The advantage of the usual method is that it requires only one infected case to obtain $\mathcal{R}0$ and $S_{\mathrm{critical}}$, whereas the next-generation matrix method requires an additional case that directly depends on an infected case ($V\left(t\right)$ in our problem) to obtain a matrix format and calculate ${\mathcal{R}}_0$. However, by using the second method, we can consider Proposition 1.

Proposition 1.

Any epidemic model can be described as:

(i) 

Does not survive, if and only if, ${\mathcal{R}}_0<<S_0<<S_{\mathrm{critical}}$.

(ii) 

Endemic, if and only if, $S_0<<{\mathcal{R}}_0<<S_{\mathrm{critical}}$.

(iii) 

Pandemic, if and only if, $S_0>>S_{\mathrm{critical}}>>{\mathcal{R}}_0$.

The basic reproduction number is a predictive calculation for how many new cases arise from an infected individual. This is the rate of infection which leads to new susceptible individuals at time $t=0$. When ${}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t_0\right)=0$, we have $aS\left(t_0\right)I\left(t_0\right)-(b+r)I\left(t_0\right)=0$, and hence ${\mathcal{R}}_0=S_{0}a/\left(b+r\right)$ is the basic reproduction number and $S_{\mathrm{critical}}=\left(b+r\right)/a$. Using two different methods, we have obtained the same basic reproduction number ${\mathcal{R}}_0$. When $S_0<S_{\mathrm{critical}}$, the disease does not survive, and when $S_0>S_{\mathrm{critical}}$, the disease recurs; that is, there is an epidemic. We consider that $S_0=400$, $S_{\mathrm{critical}}=100$, and ${\mathcal{R}}_0=4$. As the basic reproduction number ${\mathcal{R}}_0$ is much less than $S_{\mathrm{critical}}$, but $S_0>>S_{\mathrm{critical}}>>{\mathcal{R}}_0$, or $400>>100>>4$, the disease does survive and increases rapidly. Therefore, our model replicates purely a pandemic model. Note that the initial population values are assumed to be similar to [3,37], which is very similar to most diseases with the same duration of time.

3.4. Equilibrium Points and Feasibility Region Analysis

In the system established in (3), we consider ${}^{\mathrm{ABC}}{D}^{{\alpha }_1}S\left(t\right)=0$, ${}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t\right)=0$, ${}^{\mathrm{ABC}}{D}^{{\alpha }_3}V\left(t\right)=0$, and ${}^{\mathrm{ABC}}{D}^{{\alpha }_4}R\left(t\right)=0$. Then, the disease free equilibrium points are either $F^{\left(0\right)}=(S^{\left(0\right)},I^{\left(0\right)},V^{\left(0\right)},R^{\left(0\right)})=(\left(b+r\right)/a,0,0,0)$ or $F^{\left(0\right)}=(S^{\left(0\right)},I^{\left(0\right)},V^{\left(0\right)},R^{\left(0\right)})=(0,0,0,0)$; that is, $F^{\left(0\right)}=(125,0,0,0)$ or $F^{\left(0\right)}=(0,0,0,0)$.

The disease-dependent equilibrium points, or the endemic equilibrium points denoted by $F^{*}=(S^{*},I^{*},V^{*},R^{*})$, are given by

$$F^{*}=(\frac{b+r}{a},I\left(0\right),\frac{rV\left(0\right)}{b-c-r},\frac{rR\left(0\right)}{c}).$$

The disease-dependent equilibrium points are computed as

$$F^{*}=(125,300,-1.15292\times {10}^{18},200).$$

The disease-dependent equilibrium points denoted by $F^{*}=(S^{*},I^{*},V^{*},R^{*})$ are established as

$$F^{*}=(\frac{b+r}{a},I\left(0\right),\frac{rV\left(0\right)}{b-c-r},\frac{rR\left(0\right)}{c}).$$

The disease-dependent equilibrium points are calculated as

$$F^{*}=(125,300,-1.15292\times {10}^{18},200).$$

Lemma 1 restricts the SIVR model solution to a feasible region, and then it detects when the outbreak occurs.

Lemma 1.

The solution of the model under consideration is restricted to the feasible region given by $E=\{(S,I,V,R)\in {\mathbb{R}}_{+}^4,0\le N\left(t\right)\le N_0\}$. Then, an epidemic outbreak occurs when $S_0>S_{\mathrm{critical}}$, where $S_{\mathrm{critical}}$ is the threshold phenomenon or an epidemic-critical community size.

Proof.

Let $N\left(t\right)=S\left(t\right)+I\left(t\right)+V\left(t\right)+R\left(t\right)$. Then, we have

$${}^{\mathrm{ABC}}{D}^{\alpha }N\left(t\right){=}^{\mathrm{ABC}}{D}^{{\alpha }_1}S\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_3}V\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_4}R\left(t\right).$$

Now, by summing all the terms presented in (3), we obtain

$$\begin{array}{ccc}\hfill {}^{\mathrm{ABC}}{D}^{\alpha }N\left(t\right)& =& -r\left(2I\right(t)+V(t)+R(t\left)\right),\hfill \\ \hfill {}^{\mathrm{ABC}}{D}^{\alpha }N\left(t\right)& \le & 0,\hfill \end{array}$$

which obviously means that $N\left(t\right)\le N_0$. In a similar fashion, when ${}^{\mathrm{ABC}}{D}^{\alpha }S\left(t\right)\le 0$, we have $S\left(t\right)\le S_0$. In the same way, we claim that $I\left(t\right)\le I_0$ and $R\left(t\right)\le R_0.$ Now, if $S_0<\left(b+r\right)/a$, then $I\left(t\right)\le 0$. Otherwise, if $S\left(t\right)>\left(b+r\right)/a$, then there arises $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$. The term ${}^{\mathrm{ABC}}{D}^{\alpha }N\left(t\right)=-r(2I\left(t\right)+V\left(t\right)+R\left(t\right))$ indicates that the natural recovery of the infection has a substantial effect in reducing the total population cases $N\left(t\right)$. Additionally, we arrive at $-2rI\left(t\right)<-r\left(V\right(t)+R(t\left)\right)$, implying that $2I\left(t\right)>V\left(t\right)+R\left(t\right)$. This is the hypothesis when the number of infected is twice as high as that of the number of vaccinated and recovered, which indicates that the infection spread is high. Nonetheless, when the natural recovery is supported by vaccination, the situation is good in overcoming a pandemic outbreak.    □

Let us gather some of the existing knowledge on the mathematical analysis of pandemic outbreaks. A system is globally asymptotically stable when local asymptotic stability is ensured at all equilibrium points. If a system is marginally stable at some equilibrium point, the controlled pandemic may lead to the next wave. If a system is completely unstable, it establishes a pandemic situation that is impossible to control using existing strategies. Remark 1 provides conditions for the existence of the next wave and whether it can be controlled or not.

Remark 1.

Not all pandemics can be controlled in a similar fashion to an epidemic, as they often exhibit waves of peaks after initial control measures. Currently, vaccines are administered before infection to protect susceptible individuals. Vaccines are not typically given during an infection, as they may initially reduce immunity before ultimately boosting it. However, in our model, we have implemented a vaccination strategy for the infected individuals and analyzed the effect of vaccination and natural recovery on the system. By considering these two recovery strategies, we can effectively control the existence of future waves, as demonstrated in our analysis of the model presented in (3). Although the first wave may be initially controlled, there is still a chance for subsequent waves to occur. In addition, we have to keep in mind the following:

(i) 

The present situation should be under control.

(ii) 

The next pandemic wave can be expected at any time.

(iii) 

If the next wave occurs, this should also be controlled.

If we linearize the expressions stated in (3) formed by the initial populations, then they seem to be the Jacobian matrix defined by

$$\begin{array}{c}\hfill J=\left(\begin{array}{cccc}-aI\left(t\right)& -aS\left(t\right)\hfill & \hfill 0& \hfill 0\\ aI\left(t\right)& aS\left(t\right)-(b+r)\hfill & \hfill 0& \hfill 0\\ 0& b-c-r\hfill & \hfill -r& \hfill 0\\ 0& c\hfill & \hfill 0& \hfill -r\end{array}\right).\end{array}$$

The characteristic polynomial of the above matrix is found to be: $p\left(\lambda \right)={\lambda }^4+0.05{\lambda }^3+0.0068{\lambda }^2+0.000244\lambda +2.4\times {10}^{-6}.$ To compute the corresponding eigenvalues, we solve the equation $p\left(\lambda \right)=0$, obtaining: $\overrightarrow{\lambda }=(-0.005+0.0772981\mathrm{i},-0.005-0.0772981\mathrm{i},-0.02,-0.02),$ where “i” is the imaginary number. Figure 2 and Figure 3 show graphical representations of the real parts of the eigenvalues of the matrix J. For the system formulated in (3), note that the four roots of the eigenvalues have only the negative real parts. This implies that the system we are examining is locally asymptotically stable for the initial populations. Now, the Jacobian matrix formed with the two disease-free equilibrium points can be stated as

$$\begin{array}{c}\hfill J\left(F^0\right)=\left(\begin{array}{cccc}-a{I}^0& -a{S}^0\hfill & \hfill 0& \hfill 0\\ a{I}^0& a{S}^0-(b+r)\hfill & \hfill 0& \hfill 0\\ 0& b-c-r\hfill & \hfill -r& \hfill 0\\ 0& c\hfill & \hfill 0& \hfill -r\end{array}\right).\end{array}$$

The characteristic polynomial when $F^0=(0,0,0,0)$ is: $p\left(\lambda \right)=0.00002\lambda +0.0024{\lambda }^2+0.09{\lambda }^3+{\lambda }^4$. With $p\left(\lambda \right)=0$, the real parts of the eigenvalues are: $\overrightarrow{\lambda }=(-0.05,-0.02$, $-0.02,0)$. Figure 4 and Figure 5 show the plots with the real and imaginary parts of the eigenvalues of the matrix $J\left(F^0\right)$, for $F^0=(0,0,0,0)$. Thus, the system presented in (3), formed with above disease-free equilibrium points, is not locally asymptotically stable but marginally stable, since all the eigenvalues are nonpositive other than the one which is zero.

The characteristic polynomial when $F^0=(125,0,0,0)$ is stated as: $p\left(\lambda \right)=0.0004{\lambda }^2+0.04{\lambda }^3+{\lambda }^4$. From computing the zeros of $p\left(\lambda \right)$, the real parts of the eigenvalues are given as: $\overrightarrow{\lambda }=(-0.02,-0.02,0,0)$. Figure 6 and Figure 7 show the graph of the real and imaginary parts of the eigenvalues of the matrix $J\left(F^0\right)$, for $F^0=(125,0,0,0)$. Thus, the system established in (3), formed with disease-free equilibrium points, is not locally asymptotically stable or unstable. However, note that it is marginally stable since all the eigenvalues are only zero and negatives without imaginary terms. The marginal stability infers that the system for the disease-free equilibrium is marginally stable but, at any time, may become unstable; that is, the possibility for the next wave is there, overcoming disease-free equilibrium points on pandemic situations.

Now, the Jacobian matrix formed with the two disease free equilibrium points is stated as

$$\begin{array}{c}\hfill J\left(F^{*}\right)=\left(\begin{array}{cccc}-a{I}^{*}& -a{S}^{*}\hfill & \hfill 0& \hfill 0\\ a{I}^{*}& a{S}^{*}-(b+r)\hfill & \hfill 0& \hfill 0\\ 0& b-c-r\hfill & \hfill -r& \hfill 0\\ 0& c\hfill & \hfill 0& \hfill -r\end{array}\right).\end{array}$$

The characteristic polynomial when $F^{*}=(125,300,-1.15292\times {10}^{18},200)$ is: $p\left(\lambda \right)=2.4\times {10}^{-6}+0.000288\lambda +0.0112{\lambda }^2+0.16{\lambda }^3+{\lambda }^4$. From the solution of $p\left(\lambda \right)=0$, the eigenvalues are: $\overrightarrow{\lambda }=(-0.06+0.0489898\mathrm{i},-0.06-0.0489898\mathrm{i},-0.02,-0.02)$.

Figure 8 and Figure 9 show the plots with the real and imaginary parts of the eigenvalues of the matrix $J\left(F^{*}\right)$, for $F^{*}=(125,300,-1.15292\times {10}^{18},200)$. Thus, the system given in (3), formed with disease-dependent equilibrium points, is locally asymptotically stable since all the eigenvalues have negative real parts.

Figure 10a,b depicts the phase portraits of the three-dimensional SIV (susceptible–infected–vaccinated) and SIR (susceptible–infected–recovered) models, respectively. The fact that both models have a stable equilibrium at the center of the box and exhibit similar patterns indicates that they are interdependent.

The system is stable for both the initial population and the disease, depending on the equilibrium points. However, the disease-free equilibrium points are only marginally stable, not asymptotically stable. Therefore, we can infer that the outbreak can be controlled. Nevertheless, there is always a possibility of the next wave of transition occurring at any time. Even in such as situation, it can be controlled with the present strategy that we are following.

4. Positive Solutions and Sensitivity Analysis

4.1. Existence of Positive Solutions

Next, we justify the positive solutions obtained. Three important results satisfying for all time $t\in [0,\infty )$ are given as follows:

(i)

The solution path traced by our SIVR model stated in (3) is always positive (bounded from below by zero) and is also bounded from above by the population size for all positive initial conditions at any time t.

(ii)

The total population remains constant even for a long time period for any arbitrary positive initial conditions.

(iii)

The solutions are all positive, bounded from below by zero, and also bounded from above by the population size.

From (3), we get the solutions of $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$ and $R\left(t\right)$ expressed as

$$\begin{array}{ccc}\hfill S\left(t\right)& =& S\left(t_0\right)+{\mathrm{e}}^{-a{\int }_{t_0}^{t}I\left(k\right)\mathrm{d}k},0<S\left(t\right)<+\infty ,\hfill \\ \hfill I\left(t\right)& =& I\left(t_0\right)+{\mathrm{e}}^{{\int }_{t_0}^t(aS\left(k\right)-(b+r))\mathrm{d}k},0<I\left(t\right)<+\infty ,\hfill \\ \hfill V\left(t\right)& =& {\mathrm{e}}^{-r(t-t_0)}(V\left(t_0\right)+{\int }_{t_0}^t(b-c-r)I\left(k\right){\mathrm{e}}^{r(k-t_0)}\mathrm{d}k),0<V\left(t\right)<+\infty ,\hfill \\ \hfill R\left(t\right)& =& {\mathrm{e}}^{-r(t-t_0)}(R\left(t_0\right)+{\int }_{t_0}^t\left(c\right)I\left(k\right){\mathrm{e}}^{r(k-t_0)}\mathrm{d}k),0<R\left(t\right)<+\infty ,\hfill \end{array}$$

where “e” is the exponential or Euler function.

We know that the total population is the sum of all cases; that is,

$$N\left(t\right)=S\left(t\right)+I\left(t\right)+V\left(t\right)+R\left(t\right).$$

Additionally, from (3), we have

$$\begin{array}{ccc}\hfill {}^{\mathrm{ABC}}{D}^{{\alpha }_i}N\left(t\right)& =& {}^{\mathrm{ABC}}{D}^{{\alpha }_1}S\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_2}I\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_3}V\left(t\right){+}^{\mathrm{ABC}}{D}^{{\alpha }_4}R\left(t\right),\hfill \\ \hfill {}^{\mathrm{ABC}}{D}^{{\alpha }_i}N\left(t\right)& =& -r\left(2I\right(t)+V(t)+R(t\left)\right),\hfill \\ \hfill {}^{\mathrm{ABC}}{D}^{{\alpha }_i}N\left(t\right)& =& -r\left(N\right(t)+I(t)-S(t\left)\right).\hfill \end{array}$$

Observe that the solution of $N\left(t\right)$ is given by

$$N\left(t\right)={\mathrm{e}}^{-r(t-t_0)}(N\left(t_0\right)-{\int }_{t_0}^{t}r(I\left(k\right)-S\left(K\right)){\mathrm{e}}^{r(k-t_0)}\mathrm{d}k).$$

Then, we have $0<N\left(t\right)<+\infty$. Thus, the total population is always positive at any time $t\in [0,\infty )$.

Notice that, as t approaches infinity, $N(\infty )$ approaches $N\left(t_0\right)$. Hence, $N\left(t\right)$ approaches zero as t tends to infinity, which implies that the total number of cases vanishes for infinitely large time t. Furthermore, from (i) and (ii), it is clear that $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$ are positively bounded and that this boundedness is independent of t.

The relationship between $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$ over time indicates that the system is strictly positive, meaning that it is bounded from below by zero and from above by the population size, as shown in Figure 11. This figure presents a parametric plot of the SIVR model with respect to t to illustrate the interdependence between the variables in the system. Since the parametric plot is not based on true values of t, the graphical plot shows the limits of all the variables as t approaches infinity. The parametric plot is represented in black to differentiate it from the solution curves, which would be shown in various colors. A solution would have four curves, each representing one of the four variables and showing the common path shared by all of them over time. The curve displayed in the figure traces a path that is common to all four variables.

4.2. Sensitivity Analysis

Next, we compare how the state variables influence the basic reproduction number ${\mathcal{R}}_0$ of the infection. The expression used to calculate the sensitivity index of a parameter $\theta$ is defined as

$$\begin{array}{ccc}\hfill {\mathcal{I}}_{\theta }& =& r_{\theta }{\mathcal{R}}_0\hfill \\ & =& \frac{\theta }{{\mathcal{R}}_0}\frac{\partial {\mathcal{R}}_0}{\partial \theta }.\hfill \end{array}$$

(6)

When applying (6) to the model proposed in the present study, we obtain that

$$\begin{array}{ccc}\hfill {\mathcal{I}}_a& =& \frac{(b+r)}{S\left(0\right)}\frac{S\left(0\right)}{(b+r)}=1,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill {\mathcal{I}}_b& =& -\frac{b}{(b+r)}=-0.6,\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill {\mathcal{I}}_r& =& -\frac{r}{(b+r)}=-0.4.\hfill \end{array}$$

(9)

From the values of ${\mathcal{I}}_a$, ${\mathcal{I}}_b$, and ${\mathcal{I}}_r$ stated in (7)–(9), respectively, we can say that the basic reproduction of secondary infection is primarily influenced by a, more than by b or r.

On the one hand, when sensitivity index ${\mathcal{I}}_a$ increases or decreases by $1\%$, ${\mathcal{R}}_0$ also increases or decreases by $1\%$. On the other hand, the rates b and c are less dominant in causing secondary infections. In Figure 12, we can see a bar plot with the values of the three indexes.

5. Numerical Simulations

5.1. Graphical Representation of $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$

Next, we show the results that were obtained from the computational experiments carried out. The software we used was MathWorks SIMULINK and Wolfram Mathematica. First-order derivatives and fractional derivatives were computed.

The system presented in (1) was solved numerically using SIMULINK, which is a Matlab add-on. The numerical method utilized for the approximations was the classical fourth-order RKM. In Figure 13, we can see, through the block diagram shown, the graphical representation of (1). As a result of the simulations carried out, the values for $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$, with ordinary derivatives, are presented in

Table 2. Simulations with the SIVR model for $t\in \{0,10,20,\dots ,60\}$ at ${\alpha }_i=1$ for $i\in \{1,2,3,4\}$.

t	$\mathit{S}\left(\mathit{t}\right)$	$\mathit{I}\left(\mathit{t}\right)$	$\mathit{V}\left(\mathit{t}\right)$	$\mathit{R}\left(\mathit{t}\right)$
0	400	300	200	100
10	81.833	419.819	163.746	118.002
20	19.0757	300.543	134.064	129.423
30	7.22778	191.081	109.762	127.787
40	3.93508	118.373	89.8658	118.29
50	2.70459	72.7313	73.5759	105.276
60	2.14891	44.5381	60.2388	91.3616

for $t\in [0,60]$ and in

Table 3. Simulations with the SIVR model for $t\in \{0,100,200,\dots ,1000\}$ at ${\alpha }_i=1$, for $i\in \{1,2,3,4\}$.

t	$\mathit{S}\left(\mathit{t}\right)$	$\mathit{I}\left(\mathit{t}\right)$	$\mathit{V}\left(\mathit{t}\right)$	$\mathit{R}\left(\mathit{t}\right)$
0	400	300	200	100
100	1.57405	6.19969	27.0671	45.7696
200	1.4975	0.0443779	3.66313	6.46463
300	1.49697	0.000317467	0.49575	0.876828
400	1.49697	$2.27097\times {10}^{-6}$	0.0670924	0.118679
500	1.49697	$1.62182\times {10}^{-8}$	0.00907998	0.0160616
600	1.49697	$1.10376\times {10}^{-10}$	0.00122884	0.00217371
700	1.49697	$2.4074\times {10}^{-12}$	0.000166305	0.000294179
800	1.49697	$2.56041\times {10}^{-13}$	0.00002250720	0.00003981320
1000	1.49697	$8.23794\times {10}^{-11}$	$4.12206\times {10}^{-7}$	$7.29126\times {10}^{-7}$

for $t\in [0,1000]$. Figure 14 shows the corresponding plots.

5.2. Laplace Adomian Decomposition Method

To solve the proposed SIVR model with fractional derivatives, we use the LADM method [9]. The functions $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$ are calculated as

$$\begin{array}{c}\hfill \begin{array}{ccc}{S}_{k+1}\hfill & =\hfill & {\mathcal{L}}^{-1}\left\{{\displaystyle -\frac{a}{u^{{\alpha }_1}}\mathcal{L}\left\{A_k\right\}\left(u\right)}\right\},\hfill \\ I_{k+1}\hfill & =\hfill & {\mathcal{L}}^{-1}\left\{{\displaystyle \frac{a}{u^{{\alpha }_2}}\mathcal{L}\left\{A_k\right\}\left(u\right)-\frac{(b+r)}{u^{{\alpha }_2}}\mathcal{L}\left\{I_k\right\}\left(u\right)}\right\},\hfill \\ V_{k+1}\hfill & =\hfill & {\mathcal{L}}^{-1}\left\{{\displaystyle \frac{(b-c-r)}{u^{{\alpha }_3}}\mathcal{L}\left\{I_k\right\}\left(u\right)-\frac{r}{u^{{\alpha }_3}}\mathcal{L}\left\{V_k\right\}\left(u\right)}\right\},\hfill \\ R_{k+1}\hfill & =\hfill & {\mathcal{L}}^{-1}\left\{{\displaystyle \frac{c}{u^{{\alpha }_4}}\mathcal{L}\left\{I_k\right\}\left(u\right)-\frac{r}{u^{{\alpha }_4}}\mathcal{L}\left\{R_k\right\}\left(u\right)}\right\},\hfill \end{array}\end{array}$$

(10)

where $A_k$ stated in (10) is called the Adomian polynomial and defined as

$$A_k=\frac{1}{k!}\left({\mathrm{d}}^k\right(\sum _{l=0}^k\left({\lambda }^l{S}_l{\lambda }^l{I}_l\right))/\mathrm{d}{\lambda }^k){|}_{\lambda =0},k\in \{0,1,\dots ,n\},$$

with $A_0=S_0{I}_0,A_1=S_0{I}_1+S_1{I}_0$, $A_2=S_0{I}_2+S_1{I}_1+S_2{I}_0$, and so on, with ${\mathcal{L}}^{-1}\left\{f\left(t\right)\right\}$ representing the inverse Laplace transform of $f\left(t\right)$. Now, we set

$$S\left(t\right)=\sum _{k=0}^{\infty }{S}_k,I\left(t\right)=\sum _{k=0}^{\infty }{I}_k,V\left(t\right)=\sum _{k=0}^{\infty }{V}_k,R\left(t\right)=\sum _{k=0}^{\infty }{R}_k$$

where

$$\begin{array}{c}\hfill \begin{array}{ccc}S\left(t\right)\hfill & =\hfill & -\frac{5.28t^{{\alpha }_1+{\alpha }_2}}{\Gamma ({\alpha }_1+{\alpha }_2+1)}+\frac{1.5552t^{2{\alpha }_1+{\alpha }_2}}{\Gamma (2{\alpha }_1+{\alpha }_2+1)}+\frac{0.6336t^{2{\alpha }_1+{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (2{\alpha }_1+{\alpha }_2+1)}-\frac{0.297216t^{3{\alpha }_1+{\alpha }_2}}{\Gamma (3{\alpha }_1+{\alpha }_2+1)}\hfill \\ & & -\frac{0.076032t^{3{\alpha }_1+{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (3{\alpha }_1+{\alpha }_2+1)}-\frac{0.076032t^{3{\alpha }_1+{\alpha }_2}\Gamma (2{\alpha }_1+{\alpha }_2+1)}{\Gamma (2{\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (3{\alpha }_1+{\alpha }_2+1)}-\hfill \\ & & \frac{0.110592t^{3{\alpha }_1+{\alpha }_2}\Gamma (2{\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_1+{\alpha }_2+1)\Gamma (3{\alpha }_1+{\alpha }_2+1)}-\frac{0.5808t^{{\alpha }_1+2{\alpha }_2}}{\Gamma ({\alpha }_1+2{\alpha }_2+1)}+\frac{0.272448t^{2{\alpha }_1+2{\alpha }_2}}{\Gamma (2{\alpha }_1+2{\alpha }_2+1)}\hfill \\ & & +\frac{0.101376t^{2{\alpha }_1+2{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}+\frac{0.069696t^{2{\alpha }_1+2{\alpha }_2}\Gamma ({\alpha }_1+2{\alpha }_2+1)}{\Gamma ({\alpha }_2+1)\Gamma ({\alpha }_1+{\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}\hfill \\ & & +\frac{0.069696t^{2{\alpha }_1+2{\alpha }_2}\Gamma ({\alpha }_1+2{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma (2{\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}-\frac{0.063888t^{{\alpha }_1+3{\alpha }_2}}{\Gamma ({\alpha }_1+3{\alpha }_2+1)}+\frac{5.76t^{2{\alpha }_1}}{\Gamma (2{\alpha }_1+1)}\hfill \\ & & -\frac{0.6912t^{3{\alpha }_1}}{\Gamma (3{\alpha }_1+1)}+\frac{0.082944t^{4{\alpha }_1}}{\Gamma (4{\alpha }_1+1)}-\frac{48.t^{{\alpha }_1}}{\Gamma ({\alpha }_1+1)}+400,\hfill \\ I\left(t\right)\hfill & =\hfill & -\frac{5.76t^{{\alpha }_1+{\alpha }_2}}{\Gamma ({\alpha }_1+{\alpha }_2+1)}+\frac{0.6912t^{2{\alpha }_1+{\alpha }_2}}{\Gamma (2{\alpha }_1+{\alpha }_2+1)}-\frac{0.082944t^{3{\alpha }_1+{\alpha }_2}}{\Gamma (3{\alpha }_1+{\alpha }_2+1)}-\frac{1.2672t^{{\alpha }_1+2{\alpha }_2}}{\Gamma ({\alpha }_1+2{\alpha }_2+1)}\hfill \\ & & -\frac{0.6336t^{{\alpha }_1+2{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma ({\alpha }_1+2{\alpha }_2+1)}+\frac{0.262656t^{2{\alpha }_1+2{\alpha }_2}}{\Gamma (2{\alpha }_1+2{\alpha }_2+1)}+\frac{0.076032t^{2{\alpha }_1+2{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}\hfill \\ & & +\frac{0.076032t^{2{\alpha }_1+2{\alpha }_2}\Gamma (2{\alpha }_1+{\alpha }_2+1)}{\Gamma (2{\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}+\frac{0.110592t^{2{\alpha }_1+2{\alpha }_2}\Gamma (2{\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_1+{\alpha }_2+1)\Gamma (2{\alpha }_1+2{\alpha }_2+1)}-\frac{0.209088t^{{\alpha }_1+3{\alpha }_2}}{\Gamma ({\alpha }_1+3{\alpha }_2+1)}\hfill \\ & & -\frac{0.069696t^{{\alpha }_1+3{\alpha }_2}\Gamma ({\alpha }_1+{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma ({\alpha }_2+1)\Gamma ({\alpha }_1+3{\alpha }_2+1)}-\frac{0.069696t^{{\alpha }_1+3{\alpha }_2}\Gamma ({\alpha }_1+2{\alpha }_2+1)}{\Gamma ({\alpha }_2+1)\Gamma ({\alpha }_1+{\alpha }_2+1)\Gamma ({\alpha }_1+3{\alpha }_2+1)}\hfill \\ & & -\frac{0.069696t^{{\alpha }_1+3{\alpha }_2}\Gamma ({\alpha }_1+2{\alpha }_2+1)}{\Gamma ({\alpha }_1+1)\Gamma (2{\alpha }_2+1)\Gamma ({\alpha }_1+3{\alpha }_2+1)}+\frac{3.63t^{2{\alpha }_2}}{\Gamma (2{\alpha }_2+1)}+\frac{0.3993t^{3{\alpha }_2}}{\Gamma (3{\alpha }_2+1)}+\frac{0.043923t^{4{\alpha }_2}}{\Gamma (4{\alpha }_2+1)}+\frac{33.t^{{\alpha }_2}}{\Gamma ({\alpha }_2+1)}+300,\hfill \\ V\left(t\right)\hfill & =\hfill & \frac{1.9984\times {10}^{-17}{t}^{{\alpha }_1+{\alpha }_2+{\alpha }_3}}{\Gamma [{\alpha }_1+{\alpha }_2+{\alpha }_3+1]}-\frac{2.39808\times {10}^{-18}{t}^{2{\alpha }_1+{\alpha }_2+{\alpha }_3}}{\Gamma [2{\alpha }_1+{\alpha }_2+{\alpha }_3+1]}+\frac{4.39648\times {10}^{-18}{t}^{{\alpha }_1+2{\alpha }_2+{\alpha }_3}}{\Gamma [{\alpha }_1+2{\alpha }_2+{\alpha }_3+1]}+\hfill \\ & & \frac{2.19824\times {10}^{-18}/\Gamma [{\alpha }_1+{\alpha }_2+1]t^{{\alpha }_1+2{\alpha }_2+{\alpha }_3}}{\Gamma [{\alpha }_1+1]\Gamma [{\alpha }_2+1]\Gamma [{\alpha }_1+2{\alpha }_2+{\alpha }_3+1]}-\frac{3.9968\times {10}^{-19}{t}^{{\alpha }_1+{\alpha }_2+2{\alpha }_3}}{\Gamma [{\alpha }_1+{\alpha }_2+2{\alpha }_3+1]}-\hfill \\ & & \frac{1.14492\times {10}^{-16}{t}^{{\alpha }_2+{\alpha }_3}}{\Gamma [{\alpha }_2+{\alpha }_3+1]}-\frac{1.25941\times {10}^{-17}{t}^{2{\alpha }_2+{\alpha }_3}}{\Gamma [2{\alpha }_2+{\alpha }_3+1]}-\frac{1.38535\times {10}^{-18}{t}^{3{\alpha }_2+{\alpha }_3}}{\Gamma [3{\alpha }_2+{\alpha }_3+1]}+\hfill \\ & & \frac{2.28983\times {10}^{-18}{t}^{{\alpha }_2+2{\alpha }_3}}{\Gamma [{\alpha }_2+2{\alpha }_3+1]}+\frac{2.51882\times {10}^{-19}{t}^{2{\alpha }_2+2{\alpha }_3}}{\Gamma [2{\alpha }_2+2{\alpha }_3+1]}-\frac{4.57967\times {10}^{-20}{t}^{{\alpha }_2+3{\alpha }_3}}{\Gamma [{\alpha }_2+3{\alpha }_3+1]}+\frac{0.08t^{2{\alpha }_3}}{\Gamma [2{\alpha }_3+1]}\hfill \\ & & -\frac{0.0016t^{3{\alpha }_3}}{\Gamma [3{\alpha }_3+1]}+\frac{0.000032t^{4{\alpha }_3}}{\Gamma [4{\alpha }_3+1]}-\frac{4.t_3^{\alpha }}{\Gamma [{\alpha }_3+1]}+200,\hfill \\ R\left(t\right)\hfill & =\hfill & \frac{1.t^{{\alpha }_4}}{\Gamma [1+{\alpha }_4]}+\frac{0.33t^{{\alpha }_2+{\alpha }_4}}{\Gamma [1+{\alpha }_2+{\alpha }_4]}-\frac{0.0576t^{{\alpha }_1+{\alpha }_2+{\alpha }_4}}{\Gamma [1+{\alpha }_1+{\alpha }_2+{\alpha }_4]}+\frac{0.006912t^{2{\alpha }_1+{\alpha }_2+{\alpha }_4}}{\Gamma [1+2{\alpha }_1+{\alpha }_2+{\alpha }_4]}+\frac{0.0363t^{2{\alpha }_2+{\alpha }_4}}{\Gamma [1+2{\alpha }_2+{\alpha }_4]}-\hfill \\ & & \frac{0.012672t^{{\alpha }_1+2{\alpha }_2+{\alpha }_4}}{\Gamma [1+{\alpha }_1+2{\alpha }_2+{\alpha }_4]}-\frac{0.006336t^{{\alpha }_1+2{\alpha }_2+{\alpha }_4}}{\Gamma [1+{\alpha }_1+{\alpha }_2]\Gamma [1+{\alpha }_1]\Gamma [1+{\alpha }_2]\Gamma [1+{\alpha }_1+2{\alpha }_2+{\alpha }_4]}+\frac{0.003993t^{3{\alpha }_2+{\alpha }_4}}{\Gamma [1+3{\alpha }_2+{\alpha }_4]}-\hfill \\ & & \frac{0.02t^{2{\alpha }_4}}{\Gamma [1+2{\alpha }_4]}-\frac{0.0066t^{{\alpha }_2+2{\alpha }_4}}{\Gamma [1+{\alpha }_2+2{\alpha }_4]}+\frac{0.001152t^{{\alpha }_1+{\alpha }_2+2{\alpha }_4}}{\Gamma [1+{\alpha }_1+{\alpha }_2+2{\alpha }_4]}-\frac{0.000726t^{2{\alpha }_2+2{\alpha }_4}}{\Gamma [1+2{\alpha }_2+2{\alpha }_4]}\hfill \\ & & +\frac{0.0004t^{3{\alpha }_4}}{\Gamma [1+3{\alpha }_4]}+\frac{0.000132t^{{\alpha }_2+3{\alpha }_4}}{\Gamma [1+{\alpha }_2+3{\alpha }_4]}-\frac{8.\times {10}^{-6}{t}^{4{\alpha }_4}}{\Gamma [1+4{\alpha }_4]}+100.\hfill \end{array}\end{array}$$

(11)

To model the formulation stated in (3), we apply the LADM method to obtain the solution. The equations produced for $\alpha \in [0,1]$ at $t\in [0,1]$ are huge, involving eleven values for $\alpha$ with eleven values for t, producing 121 terms for the four cases. Thus, we have shown them in Figure 15, Figure 16, Figure 17, Figure 18 and Figure 19.

Using the solutions to the ODEs for ${\alpha }_1={\alpha }_2={\alpha }_3={\alpha }_4=1$ and by utilizing the fourth order Laplace Adomian decomposition method (note that increasing the order of the method, we can gain a better accuracy of the solution), as given by

$$\begin{array}{c}\hfill \begin{array}{ccc}S\left(t\right)\hfill & =\hfill & 400-48t+0.24t^2+0.2584t^3-0.00403t^4,\hfill \\ I\left(t\right)\hfill & =\hfill & 300+33t-1.065t^2-0.24065t^3+0.00703812t^4,\hfill \\ V\left(t\right)\hfill & =\hfill & 200-4t+0.04t^2-0.000266667t^3+1.33333\times {10}^{-}6t^4,\hfill \\ R\left(t\right)\hfill & =\hfill & 100+t+0.155t^2-0.00458333t^3-0.000578708t^4.\hfill \end{array}\end{array}$$

(12)

Table 4. $S\left(t\right)$ at ${\alpha }_1\in [0,1]$ and $t\in [0,1]$.

	${}^{\mathbf{ABC}}{\mathit{D}}^{{\mathbf{\alpha }}_1}\mathit{S}\left(\mathit{t}\right)$
${\mathbf{\alpha }}_{\mathbf{1}}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{0}.\mathbf{1}$	$\mathit{t}=\mathbf{0}.\mathbf{2}$	$\mathit{t}=\mathbf{0}.\mathbf{3}$	$\mathit{t}=\mathbf{0}.\mathbf{4}$	$\mathit{t}=\mathbf{0}.\mathbf{5}$	$\mathit{t}=\mathbf{0}.\mathbf{6}$	$\mathit{t}=\mathbf{0}.\mathbf{7}$	$\mathit{t}=\mathbf{0}.\mathbf{8}$	$\mathit{t}=\mathbf{0}.\mathbf{9}$	$\mathit{t}=\mathbf{1}$
0	353.4	353.4	353.4	353.4	353.4	353.4	353.4	353.4	353.4	353.4	353.4
0.1	400.0	360.8	358.0	356.4	355.2	354.2	353.4	352.7	352.1	351.6	351.1
0.2	400.0	367.5	362.8	359.7	357.4	355.6	354.0	352.6	351.4	350.3	349.3
0.3	400.0	373.5	367.4	363.3	360.1	357.4	355.1	353.1	351.3	349.6	348.0
0.4	400.0	378.6	371.9	367.0	363.1	359.7	356.7	354.1	351.7	349.4	347.3
0.5	400.0	382.9	375.9	370.6	366.1	362.2	358.7	355.5	352.5	349.7	347.1
0.6	400.0	386.5	379.6	374.1	369.3	364.9	361.0	357.3	353.8	350.5	347.4
0.7	400.0	389.5	382.9	377.4	372.4	367.8	363.4	359.3	355.4	351.7	348.1
0.8	400.0	391.8	385.8	380.4	375.4	370.6	366.0	361.6	357.4	353.3	349.2
0.9	400.0	393.7	388.3	383.2	378.2	373.4	368.7	364.1	359.5	355.1	350.7
1	400.0	395.2	390.4	385.6	380.9	376.1	371.3	366.6	361.9	357.2	352.5

,

Table 5. $I\left(t\right)$ at ${\alpha }_2\in [0,1]$ and $t\in [0,1]$.

	${}^{\mathbf{ABC}}{\mathit{D}}^{{\mathbf{\alpha }}_2}\mathit{I}\left(\mathit{t}\right)$
${\mathbf{\alpha }}_{\mathbf{2}}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{0}.\mathbf{1}$	$\mathit{t}=\mathbf{0}.\mathbf{2}$	$\mathit{t}=\mathbf{0}.\mathbf{3}$	$\mathit{t}=\mathbf{0}.\mathbf{4}$	$\mathit{t}=\mathbf{0}.\mathbf{5}$	$\mathit{t}=\mathbf{0}.\mathbf{6}$	$\mathit{t}=\mathbf{0}.\mathbf{7}$	$\mathit{t}=\mathbf{0}.\mathbf{8}$	$\mathit{t}=\mathbf{0}.\mathbf{9}$	$\mathit{t}=\mathbf{1}$
0	330.1	330.1	330.1	330.1	330.1	330.1	330.1	330.1	330.1	330.1	330.1
0.1	300.0	325.7	327.3	328.3	329.1	329.7	330.1	330.6	330.9	331.2	331.5
0.2	300.0	321.5	324.5	326.3	327.8	328.9	329.9	330.7	331.4	332.1	332.7
0.3	300.0	317.7	321.6	324.2	326.2	327.8	329.3	330.5	331.6	332.6	333.5
0.4	300.0	314.4	318.8	321.9	324.4	326.5	328.4	330.0	331.5	332.9	334.1
0.5	300.0	311.5	316.2	319.6	322.5	325.0	327.2	329.3	331.1	332.8	334.4
0.6	300.0	309.1	313.7	317.4	320.6	323.4	325.9	328.2	330.4	332.5	334.4
0.7	300.0	307.2	311.6	315.3	318.6	321.6	324.4	327.0	329.5	331.9	334.1
0.8	300.0	305.6	309.7	313.3	316.6	319.8	322.8	325.6	328.4	331.0	333.5
0.9	300.0	304.3	308.0	311.5	314.8	318.0	321.1	324.1	327.1	329.9	332.7
1	300.0	303.3	306.6	309.8	313.0	316.2	319.4	322.5	325.6	328.7	331.7

,

Table 6. $V\left(t\right)$ at ${\alpha }_3\in [0,1]$ and $t\in [0,1]$.

	${}^{\mathbf{ABC}}{\mathit{D}}^{{\mathbf{\alpha }}_3}\mathit{V}\left(\mathit{t}\right)$
${\mathbf{\alpha }}_{\mathbf{3}}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{0}.\mathbf{1}$	$\mathit{t}=\mathbf{0}.\mathbf{2}$	$\mathit{t}=\mathbf{0}.\mathbf{3}$	$\mathit{t}=\mathbf{0}.\mathbf{4}$	$\mathit{t}=\mathbf{0}.\mathbf{5}$	$\mathit{t}=\mathbf{0}.\mathbf{6}$	$\mathit{t}=\mathbf{0}.\mathbf{7}$	$\mathit{t}=\mathbf{0}.\mathbf{8}$	$\mathit{t}=\mathbf{0}.\mathbf{9}$	$\mathit{t}=\mathbf{1}$
0	196.1	196.1	196.1	196.1	196.1	196.1	196.1	196.1	196.1	196.1	196.1
0.1	200.0	196.7	196.5	196.3	196.2	196.2	196.1	196.0	196.0	195.9	195.9
0.2	200.0	197.3	196.9	196.6	196.4	196.3	196.1	196.0	195.9	195.8	195.7
0.3	200.0	197.8	197.3	196.9	196.7	196.4	196.2	196.1	195.9	195.8	195.6
0.4	200.0	198.2	197.7	197.2	196.9	196.6	196.4	196.2	195.9	195.8	195.6
0.5	200.0	198.6	198.0	197.6	197.2	196.8	196.6	196.3	196.0	195.8	195.6
0.6	200.0	198.9	198.3	197.8	197.4	197.1	196.7	196.4	196.1	195.9	195.6
0.7	200.0	199.1	198.6	198.1	197.7	197.3	197.0	196.6	196.3	196.0	195.7
0.8	200.0	199.3	198.8	198.4	197.9	197.6	197.2	196.8	196.4	196.1	195.8
0.9	200.0	199.5	199.0	198.6	198.2	197.8	197.4	197.0	196.6	196.3	195.9
1	200.0	199.6	199.2	198.8	198.4	198.0	197.6	197.2	196.8	196.4	196.0

and

Table 7. $R\left(t\right)$ at ${\alpha }_4\in [0,1]$ and $t\in [0,1]$.

	${}^{\mathbf{ABC}}{\mathit{D}}^{{\mathbf{\alpha }}_4}\mathit{R}\left(\mathit{t}\right)$
${\mathbf{\alpha }}_{\mathbf{4}}$	$\mathit{t}=\mathbf{0}$	$\mathit{t}=\mathbf{0}.\mathbf{1}$	$\mathit{t}=\mathbf{0}.\mathbf{2}$	$\mathit{t}=\mathbf{0}.\mathbf{3}$	$\mathit{t}=\mathbf{0}.\mathbf{4}$	$\mathit{t}=\mathbf{0}.\mathbf{5}$	$\mathit{t}=\mathbf{0}.\mathbf{6}$	$\mathit{t}=\mathbf{0}.\mathbf{7}$	$\mathit{t}=\mathbf{0}.\mathbf{8}$	$\mathit{t}=\mathbf{0}.\mathbf{9}$	$\mathit{t}=\mathbf{1}$
0	101.3	101.3	101.3	101.3	101.3	101.3	101.3	101.3	101.3	101.3	101.3
0.1	100.0	101.0	101.1	101.2	101.2	101.2	101.3	101.3	101.3	101.3	101.4
0.2	100.0	100.8	101.0	101.1	101.1	101.2	101.2	101.3	101.3	101.4	101.4
0.3	100.0	100.6	100.8	100.9	101.0	101.1	101.2	101.3	101.3	101.4	101.4
0.4	100.0	100.5	100.7	100.8	100.9	101.0	101.1	101.2	101.3	101.4	101.4
0.5	100.0	100.4	100.6	100.7	100.8	100.9	101.0	101.1	101.2	101.3	101.4
0.6	100.0	100.3	100.5	100.6	100.7	100.9	101.0	101.1	101.2	101.3	101.4
0.7	100.0	100.2	100.4	100.5	100.6	100.8	100.9	101.0	101.1	101.2	101.3
0.8	100.0	100.2	100.3	100.4	100.6	100.7	100.8	100.9	101.0	101.2	101.3
0.9	100.0	100.1	100.3	100.4	100.5	100.6	100.7	100.8	101.0	101.1	101.2
1	100.0	100.1	100.2	100.3	100.4	100.5	100.7	100.8	100.9	101.0	101.2

report the values obtained from (11) for the functions $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$, respectively, with ${\alpha }_i\in \{0.0,0.1,0.2,\dots ,0.9,1.0\}$ and $t\in [0,1]$. Figure 15 shows the plots of the functions $S\left(t\right)$, $I\left(t\right)$, $V\left(t\right)$, and $R\left(t\right)$, with ${\alpha }_i\in \{0.0,0.1,0.2,\dots ,0.9,1.0\}$ and $t\in [0,1]$. Figure 16, Figure 17 and Figure 18 show, respectively, the plots of the same functions with $t\in [1,2]$, $t\in [2,5]$ and $t\in [0,30]$. Again, the values of ${\alpha }_i$ are taken from the set $\{0.0,0.1,0.2,\dots ,0.9,1.0\}$. Figure 19 provides the corresponding 3D plots. Additionally, we give the ODEs stated in (12).

6. Results, Discussion, and Conclusions

In the present investigation, the solutions obtained with the LADM were used to plot fractional illustrations for various values of $\alpha$ as seen in Figure 15, Figure 16, Figure 17 and Figure 18. Other than stability analysis, the feasibility region analysis, the existence of positive solutions (see Figure 11), and the sensitivity analysis (see Figure 12) were also performed. From the sensitivity analysis, we showed that the coefficient of S, that is, the constant a, is the most influential parameter of the basic reproduction number ${\mathcal{R}}_0.$ New propositions and theorems were stated to claim that the epidemic model is replicating the pandemic spread. From Figure 14a, we observed that the number of recovered cases increased more than did the number of infected cases at around day #40 due to a consistent decrease in vaccination. At approximately day #50, all infected cases were fully vaccinated. We presented the depletion vaccination strategy model because it is not necessary to increase daily vaccination cases when the number of infected and susceptible cases is reduced. Additionally, when the spread is severe, more people can be vaccinated, but when the spread is less, people may not be interested in vaccinating unnecessarily. Moreover, we mathematically established the chance of the next wave of any pandemic disease and showed that it could be controlled by a consistent vaccination strategy.

Typically, vaccines have been given to uninfected individuals to prevent them from contracting COVID-19. Individuals who have recovered from the virus may not need to be vaccinated for a certain period of time since they have natural immunity. Vaccines can temporarily reduce immunity before boosting it again. When a person is already infected, they are generally not recommended to receive vaccines for other diseases that require a strong immune system. However, for diseases like rabies, individuals may be vaccinated even while they have an active infection. Currently, vaccines for actively infected COVID-19 cases are not available in most countries. The usual recommendation is self-quarantine or medical quarantine under doctor’s supervision. We are interested in expanding our current research to specifically investigate the implications of vaccines for actively infected COVID-19 cases.

We utilized the LADM to obtain numerical simulations in our investigation, but other methods such as homotopy perturbation and differential transform can also be used to visualize the rate of change in the spread of epidemics over time.

Figure 15, Figure 16, Figure 17 and Figure 18 depicted the fractional changes (fixed fractional order) in daily SIVR cases over a period of time in days, while Figure 19 showed the daily changes in SIVR cases (varying fractional order). The ordinary plots of these SIVR cases were provided in Figure 14.

From Figure 15, Figure 16, Figure 17 and Figure 18, it was evident that from day 0 until day 30, the susceptible count decreased and increased for ${\alpha }_1$ over time t, while the infected (${\alpha }_2$ over time t) and recovered population (${\alpha }_3$ over time t) increased and decreases, respectively. In contrast, the vaccinated cases (${\alpha }_4$ over time t) decreased. These trends reflect the biological dynamics of the SIVR model.

We presented a novel SIVR model that employs consistent depletion vaccination strategies to control infections. The model was analyzed using necessary mathematical theories and validated with numerical simulations. Our motivation for developing the SIVR model stemmed from diseases such as rabies, Ebola, and diabetes, for which compartmental models, including susceptible, infected, vaccinated, and recovered cases, are commonly used.

Initially, the model was formulated as an ordinary differential system. However, to capture the fractional changes that occur in our SIVR pandemic cases, we reformulated it as a fractional ordered system. Our analysis showed that consistent depletion vaccination strategies were highly effective in reducing infections, and after a certain period of time, the diseases decreased steadily and eventually vanished completely.

Based on our findings, we recommend using fractional modeling to capture the very small changes (fractional models) in infectious disease dynamics that lead to larger changes observed in ordinary differential models.

In the future, we plan to explore the immunity developed after COVID-19 infection, as discussed in [44], by considering the SIR model without the susceptible population. Additionally, we will investigate models that incorporate the dead population, as discussed in [37,45], as well as other related methods [46].