Fractional SIAQR model with time dependent infection rate

Abstract

In this study we have proposed the SIAQR model with time-dependent infection rate. The model’s well-posedness is demonstrated, and the proposed model is then extended to the fractional mathematical model by employing the Caputo-Fabrizio fractional derivative operator. We have also used the Lipschitz condition and the linear growth condition to determine the conditions under which the model has a unique solution. The numerical solutions are presented. Furthermore, using graphical representations we will see how the number of infected, recovered individuals vary as the order of fractional derivative varies.

2010 Mathematics Subject Classification:

• Primary 92B05
• 92C60
• Secondary 26A33

Keywords:

• Fractional derivative
• numerical solution
• SIAQR model
• existence and uniqueness

1. Introduction

As we know that mathematical models are capable of decision-making, saving lives, assisting in policy, and many more. These are helpful in understanding the conditions needed to sustain lives and provide us with ways to study and predict the behavior of the spread. The concept of derivatives and integrals plays a lot in the formulation of these mathematical models. In this work, we will be concerned with the SIAQR mathematical model. A SIAQR model is an epidemiological model that determines the number of people infected with a transmissible infection in a closed population over a period of time. This class of model derives from the fact that they involve equations relating the number of susceptible people S(t), the number of people infected symptomatically I(t), the number of people who are asymptomatically infected A(t), the number of people who have been quarantined Q(t), and who have recovered R(t). Recently the concept of fractional derivatives (Kilbas et al., 2006; Miller and Ross, 1993; Oldham and Spanier, 1974) and integrals are proved to be very useful in making a better understanding of real-world problems which exhibit non-local behaviors. The first and most widely accepted definition of the singular kernel was proposed by Riemann-Liouville, later on, many new developments were made in the field of fractional calculus which involves the definition of the fractional derivative with the non-singular kernel, one of which was given by Caputo. This definition is based on the concept of the power law. The literature is full-fledged with many definitions of fractional derivatives and integrals. Applications can be seen in the literature (Ahmad et al., 2021; Ain & He, 2019; Ain et al., 2021, 2022a, 2022b, 2022c; Alkahtani & Alzaid, 2020, 2021; Alshomrani et al., 2021; Al-Smadi et al., 2021; Anjum et al., 2021a, 2021b; Atangana & Aguilar, 2018; Atangana & Araz, 2021; Atangana & Atangana, 2020; Djeddi et al., 2020; Hasan et al., 2022; Jarad et al., 2018; Karaji & Nyamoradi, 2020; Khan & Atangana, 2020; Modi et al., 2021; Momani et al., 2021; Naik et al., 2020, 2021, 2022; Ucar et al., 2019). In this article, we consider the SIAQR model with time-dependent infection rate. An analysis of this SIAQR model is presented in detail. Using the Caputo definition of fractional derivative, we will study the fractional SIAQR model. In Section 1, an introduction along with some definitions of the fractional derivatives are given. In Section 2, we have presented the SIAQR model with a time-dependent infection rate. In Section 3, we have presented the existence and the uniqueness of the system of solution of the SIAQR model. In Section 4, the SIAQR model with time-dependent infection rate in sense of the Caputo fractional derivative operator is given and numerical simulations are performed for a better understanding of the model.

Following are some definitions of fractional derivatives and integrals.

Definition 1.

For an integrable function f on $\mathbb{R},\gamma >0$ the Riemann-Liouville fractional integral of order γ is given as (1.1) $${}^{RL}{\mathcal{I}}_t^{\gamma }f(t)=\frac{1}{\Gamma (\gamma )}{\int }_0^{t}f(\tau ){(t-\tau )}^{\gamma -1}d\tau ,$$(1.1) where ${}^{RL}{\mathcal{I}}_t^{\gamma }$ denotes the Riemann-Liouville fractional integral of order γ (Miller and Ross, 1993).

Definition 2.

For an integrable function f on $\mathbb{R},0<\gamma <1$ the Riemann-Liouville fractional derivative of order γ is given as (1.2) $${}^{\mathcal{R}\mathcal{L}}{\chi }_t^{\gamma }f(t)=\frac{1}{\Gamma (1-\gamma )}\frac{d}{d\tau }{\int }_0^t\frac{f(\tau )}{{(t-\tau )}^{\gamma }}d\tau ,$$(1.2) where ${}^{RL}{\chi }_t^{\gamma }$ denotes the Riemann-Liouville fractional derivative of order γ (Miller and Ross, 1993).

Definition 3.

For an integrable function f on $\mathbb{R},0<\gamma <1$ the Caputo fractional derivative of order γ is given as (1.3) $${}^{\mathcal{C}}{\chi }_t^{\gamma }f(t)=\frac{1}{\Gamma (1-\gamma )}{\int }_0^t\frac{1}{{(t-\tau )}^{\gamma }}\frac{d}{d\tau }f(\tau )d\tau ,$$(1.3) where ${}^{\mathcal{C}}{\chi }_t^{\gamma }$ denotes the Caputo fractional derivative of order γ (Miller and Ross, 1993).

Theorem 1.

Assuming that there exists positive constants L and $\overline{L}$ such that the following holds

• Lipschitz condition: $\forall v_1,v_2\in \mathbb{R}$ and $\forall t\in [t_0,T]$ (1.4) $$|g(t,v_1(t))-g(t,v_2(t))|\le L|v_1-v_2|.$$(1.4)

• Linear growth condition: $\forall (v,t)\in \mathbb{R}\times [[t_0,T]]$ (1.5) $$|g(t,v){|}^2\le \overline{L}(1+|v{|}^2).$$(1.5)

Then the Cauchy problem with Caputo derivative admits a unique solution (Atangana, 2020).

2. SIAQR model with time dependent infection rate

In the SIAQR model, we have considered the five classes of individuals as S: Susceptible group of individuals I: Infected symptomatic group of individuals A: Asymptomatically infected group of individuals, Q: Quarantined group of individuals and R: Recovered group of individuals. The following assumptions are made in this model:

1. We assume that the total population N(t) does not change much with time i.e. $N(t)=S(t)+I(t)+A(t)+Q(t)+R(t).$

2. There is the interaction among every person.

3. The infected symptomatically and asymptomatically infected group of individuals are transmitting the infection to a susceptible group of individuals.

4. It is also assumed that after recovery the infected individuals develop immunity and eventually the recovered individuals are prone to lose their immunity and return to the susceptible class of individuals after a period of time.

5. Further, a time dependence on the parameter β₁ and β₂, the transmission rate of infection, is introduced in this model, so the effect of transmission will be given as ${\beta }_1 e^{zt}$ and ${\beta }_2 e^{zt}.$ This choice of decreasing exponential function is justified by step by step implementation of rules and by a rise in consciousness and awareness in individuals. Here, β₁ and β₂ is the infection rate at the beginning of the infection and z determines the change in the time of infection (Iannni & Rossi, 2020).

Keeping in mind the aforementioned assumptions the SIAQR model is formulated as: $$\begin{array}{c}\frac{dS}{dt}=\delta N(t)-{\beta }_1{e}^{-z\mathrm{ }t}\frac{S(t)I(t)}{N}-{\beta }_2{e}^{-z\mathrm{ }t}\frac{S(t)A(t)}{N}+\omega R(t)-\delta S(t),S(0)=S_0,\\ \frac{dI}{dt}={\beta }_1{e}^{-z\mathrm{ }t}\frac{S(t)I(t)}{N}-r_{1}I(t)-qI(t)-\delta I(t),I(0)=I_0,\\ \frac{dA}{dt}={\beta }_2{e}^{-z\mathrm{ }t}\frac{S(t)A(t)}{N}-r_{2}A(t)-qA(t)-\delta A(t),A(0)=A_0,\\ \frac{dQ}{dt}=qI(t)+qA(t)-r_{3}Q(t)-\delta Q(t),Q(0)=Q_0,\\ \frac{dR}{dt}=r_{1}I(t)+r_{2}A(t)+r_{3}Q(t)-\omega R(t)-\delta R(t),R(0)=R_0,\end{array}$$ where, δ defines the birth and natural death rate, $r_1,r_2,r_3,$ represents the recovery rates for infected symptomatically, asymptomatically infected, and quarantined groups of individuals respectively, q is the rate by which infected people are quarantined. β₁ and β₂ are the infection rate for the individuals infected symptomatic and asymptomatically. Further, a time dependence on the infection rate β₁ and β₂ is introduced with an exponential decay behaviour $e^{-z\mathrm{ }t},$ where, z is the rate of decrease.

3. SIAQR model with time-dependent infection rate in sense of Caputo fractional derivative operator

We now extend the above formulated SIAQR model to the fractional SIAQR model in sense of Caputo fractional derivative, which is given as follows: (3.1) $$\begin{array}{c}{}^{\mathcal{C}}{\chi }_t^{\gamma }S(t)=\delta N(t)-{\beta }_1{e}^{-z\mathrm{ }t}\frac{S(t)I(t)}{N}-{\beta }_2{e}^{-z\mathrm{ }t}\frac{S(t)A(t)}{N}+\omega R(t)-\delta S(t),S(0)=S_0,\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }I(t)={\beta }_1{e}^{-z\mathrm{ }t}\frac{S(t)I(t)}{N}-r_{1}I(t)-qI(t)-\delta I(t),I(0)=I_0,\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }A(t)={\beta }_2{e}^{-z\mathrm{ }t}\frac{S(t)A(t)}{N}-r_{2}A(t)-qA(t)-\delta A(t),A(0)=A_0,\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }Q(t)=qI(t)+qA(t)-r_{3}Q(t)-\delta Q(t),Q(0)=Q_0,\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }R(t)=r_{1}I(t)+r_{2}A(t)+r_{3}Q(t)-\omega R(t)-\delta R(t),R(0)=R_0,\end{array}$$(3.1)

Here ${}^{\mathcal{C}}{\chi }_t^{\gamma }$ denotes the Caputo fractional derivative of order $\gamma .$

Theorem 2.

Assuming that all the parameters defined in the above formulated model are positive, the solution of the model $S(t),I(t),A(t),Q(t),R(t)$ with non-negative initial conditions are positive $\forall t>0.$

Proof.

We consider the first equation (3.2) $$\begin{array}{c}\frac{dS}{dt}=\delta N(t) -{\beta }_1 e^{-z\mathrm{ }t} \frac{S I}{N}-{\beta }_2 e^{-z\mathrm{ }t} \frac{S A}{N} +\omega R(t) -\delta S(t) \ge -\delta S(t),\\ \Rightarrow \\ \frac{dS}{dt}=\ge -\delta S(t).\end{array}$$(3.2)

On solving, we get $S(t)=C_1 e^{-\delta t},$ where C₁ is constant of integration.

Using the fact that $S(0)=S_0>0,$ we get (3.3) $$S(t)=S_0 e^{-\delta \mathrm{ }t}\ge 0.$$(3.3)

In a similar way, we consider the second equation (3.4) $$\begin{array}{c}\frac{dI}{dt}=-{\beta }_1 e^{-z\mathrm{ }t} \frac{S I}{N}-r_1 I(t)-q I(t)-\delta I(t)\ge -(r_1+q+\delta )I(t),\\ \Rightarrow \\ \frac{dI}{dt}\ge -(r_1+q+\delta )I(t).\end{array}$$(3.4)

On solving, we get $I(t)=C_2 e^{-(r_1+q+\delta )t},$ where C₂ is constant of integration.

Using the fact that $I(0)=I_0>0,$ we get (3.5) $$I(t)=I_0 e^{-(r_1+q+\delta )t}\ge 0.$$(3.5)

On a similar account, we get (3.6) $$A(t)=A_0 e^{-(r_2+q+\delta )t} \ge 0,$$(3.6) (3.7) $$Q(t)=Q_0 e^{-(r_3+\delta )t} \ge 0,$$(3.7) (3.8) $$R(t)=R_0 e^{-(w+\delta )t} \ge 0.$$(3.8)

□

4. Existence and uniqueness of the fractional SIAQR model

In this section we will prove the existence and the uniqueness for the SIAQR model. For convenience, we write the SIAQR model as (4.1) $$\begin{array}{c}{}^{\mathcal{C}}{\chi }_t^{\gamma }S(t)=g_1(t,S),\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }I(t)=g_2(t,I),\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }A(t)=g_3(t,A),\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }Q(t)=g_4(t,Q),\\ {}^{\mathcal{C}}{\chi }_t^{\gamma }R(t)=g_5(t,R),\end{array}$$(4.1) where (4.2) $$\begin{array}{c}{g}_1(t,S)=\delta N(t) -{\beta }_1 e^{-z\mathrm{ }t} \frac{S(t) I(t)}{N}-{\beta }_2 e^{-z\mathrm{ }t} \frac{S(t) A(t)}{N} +\omega R(t)-\delta S(t),\\ g_2(t,I)={\beta }_1 e^{-z\mathrm{ }t} \frac{S(t) I(t)}{N}-r_1 I(t)-q I(t)-\delta I(t),\\ g_3(t,A)={\beta }_2 e^{-z\mathrm{ }t} \frac{S(t) A(t)}{N}-r_2 A(t)-q A(t)-\delta A(t),\\ g_4(t,Q)=q I(t)+q A(t)-r_3 Q(t)-\delta Q(t),\\ g_5(t,R)=r_1 I(t)+r_2 A(t)+r_3 Q(t)-\omega R(t)-\delta R(t).\end{array}$$(4.2)

To prove the existence and the uniqueness of the solution of the SIAQR model, we use the concept which was recently proposed by Atangana (2020) and hence prove the following theorem.

Theorem 3.

Assuming that there exists positive constants $L_1,L_2,L_3,L_4,L_5$ and $\overline{L_1},\overline{L_2},\overline{L_3},\overline{L_4},\overline{L_5}$ such that the following holds

1. (4.3) $$\begin{array}{c}|g_1(t,S(t))-g_1(t,S_1(t))|\le L_1|S-S_1|,\\ |g_2(t,I(t))-g_2(t,I_1(t))|\le L_2|I-I_1|,\\ |g_3(t,A(t))-g_3(t,A_1(t))|\le L_3|A-A_1|,\\ |g_4(t,Q(t))-g_4(t,Q_1(t))|\le L_4|Q-Q_1|,\\ |g_5(t,R(t))-g_5(t,R_1(t))|\le L_5|R-R_1|.\end{array}$$(4.3) (4.4) $$\begin{array}{c}|g_1(t,S(t)){|}^2\le \overline{L_1}(1+|S{|}^2),\\ |g_2(t,I(t)){|}^2\le \overline{L_2}(1+|I{|}^2),\\ |g_3(t,A(t)){|}^2\le \overline{L_3}(1+|A{|}^2),\\ |g_4(t,Q(t)){|}^2\le \overline{L_4}(1+|Q{|}^2),\\ |g_5(t,R(t)){|}^2\le \overline{L_5}(1+|R{|}^2).\end{array}$$(4.4)

Proof.

Consider (4.5) $$|g_1(t,S)-g_1(t,S_1)|=\left|\left(-{\beta }_1{e}^{-z\mathrm{ }t}\frac{I(t)}{N(t)}-{\beta }_2{e}^{-z\mathrm{ }t}\frac{A(t)}{N(t)}-\delta \right)(S(t)-S_1(t))\right|.$$(4.5)

Define the norm as $\left|\right|w|{|}_{\infty }=Su{p}_{t\in [0,T]}|w|,$ we get (4.6) $$|g_1(t,S)-g_1(t,S_1)|\le \left({\beta }_1\frac{\left|\right|I\left|\right|}{\left|\right|N\left|\right|}+{\beta }_2\frac{\left|\right|A\left|\right|}{\left|\right|N\left|\right|}\right)|S-S_1|.$$(4.6)

Taking $\left({\beta }_1\frac{\left|\right|I\left|\right|}{\left|\right|N\left|\right|}+{\beta }_2\frac{\left|\right|A\left|\right|}{\left|\right|N\left|\right|}\right)=L_1,$ we have (4.7) $$|g_1(t,S)-g_1(t,S_1)\le L_1|S-S_1|.$$(4.7)

Next consider (4.8) $$|g_2(t,I)-g_2(t,I_1)|=\left|\left({\beta }_1{e}^{-z*t}\frac{S}{N}-r_1-q-\delta \right)(I(t)-I_1(t))\right|.$$(4.8)

Again defining the norm as $\left|\right|w|{|}_{\infty }=Su{p}_{t\in [0,T]}|w|,$ we get (4.9) $$|g_2(t,I)-g_2(t,I_1)|=\left({\beta }_1\frac{\left|\right|S\left|\right|}{\left|\right|N\left|\right|}+r_1+q+\delta \right)|I-I_1|.$$(4.9)

Take $\left({\beta }_1\frac{\left|\right|S\left|\right|}{\left|\right|N\left|\right|}+r_1+q+\delta \right)=L_2,$ finally we get (4.10) $$|g_2(t,I)-g_2(t,I_1)|=L_2|I-I_1|.$$(4.10)

Next consider, (4.11) $$|g_3(t,A)-g_3(t,A_1)|=\left|\left({\beta }_2{e}^{-z*t}\frac{S}{N}-r_2-q-\delta \right)(A(t)-A_1(t))\right|.$$(4.11)

On a similar account, we get (4.12) $$|g_3(t,A)-g_3(t,A_1)|=L_3|A-A_1|,$$(4.12) where $L_3=\left({\beta }_2\frac{\left|\right|S\left|\right|}{\left|\right|N\left|\right|}+r_2+q+\delta \right).$

Similarly, we can show (4.13) $$|g_4(t,Q)-g_4(t,Q_1)|=L_4|Q-Q_1|,$$(4.13) and (4.14) $$|g_5(t,R)-g_5(t,R_1)|=L_5|R-R_1|.$$(4.14)

We now prove the second part of the above stated theorem. We first show $|g_1(t,S(t))|\le \overline{L_1}(1+|S{|}^2).$ Consider (4.15) $$\begin{array}{c}|g_1(t,S){|}^2={\left|\delta N(t)+\omega R(t)+\left(-{\beta }_1{e}^{-z\mathrm{ }t}\frac{I}{N}-{\beta }_2{e}^{-z\mathrm{ }t}\frac{A}{N}-\delta \right)S\right|}^2\\ \le 3{\delta }^2|N{|}^2+3{\omega }^2|R{|}^2+3\left({\beta }_1^2{e}^{-2zt}\frac{|I{|}^2}{|N{|}^2}+{\beta }_2^2{e}^{-2zt}\frac{|A{|}^2}{|N{|}^2}+{\delta }^2\right)|S{|}^2\\ \le 3{\delta }^2\left|\right|N|{|}_{\infty }^2+3{\omega }^2|\left|R\right|{|}_{\infty }^2+\left(3{\beta }_1^2\frac{\left|\right|I|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+{\beta }_2^2\frac{\left|\right|A|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+{\delta }^2\right)|S{|}^2\\ =3{\delta }^2\left|\right|N|{|}_{\infty }^2+3{\omega }^2|\left|R\right|{|}_{\infty }^2\left[1+\frac{\left(3{\beta }_1^2\frac{\left|\right|I|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+3{\beta }_2^2\frac{\left|\right|A|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+3{\delta }^2\right)|S(t){|}^2}{3{\delta }^2\left|\right|N|{|}_{\infty }^2+3{\omega }^2|\left|R\right|{|}_{\infty }^2}\right].\end{array}$$(4.15)

On taking $\overline{L_1}=3({\delta }^2|\left|N\right|{|}_{\infty }^2+{\omega }^2\left|\right|R|{|}_{\infty }^2)$ and under the condition $$\frac{\left(3{\beta }_1^2\frac{\left|\right|I|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+3{\beta }_2^2\frac{\left|\right|A|{|}_{\infty }^2}{\left|\right|N|{|}_{\infty }^2}+3{\delta }^2\right)}{3{\delta }^2\left|\right|N|{|}_{\infty }^2+3{\omega }^2|\left|R\right|{|}_{\infty }^2}\le 1,$$ we get (4.16) $$|g_1(t,S){|}^2\le \overline{L_1}(1+|S{|}^2).$$(4.16)

We now show (4.17) $$|g_2(t,I(t)){|}^2\le \overline{L_2}(1+|I(t){|}^2).$$(4.17)

Consider (4.18) $$\begin{array}{c}|g_2(t,I){|}^2={\left|\left({\beta }_1{e}^{-zt}\frac{S I}{N}-r_{1}I-qI-\delta I\right)\right|}^2\\ \le 2{\beta }_1^2{e}^{-2zt}\frac{\left|S{|}^2\right|I{|}^2}{|N{|}^2}+3(r_1^2+q^2+{\delta }^2)|I{|}^2\\ =\left(2{\beta }_1^2{e}^{-2zt}\frac{|S{|}^2}{|N{|}^2}+3(r_1^2+q^2+{\delta }^2)\right)|I{|}^2\\ \le (2{\beta }_1^2+3(r_1^2+q^2+{\delta }^2))I^2\\ \le |N{|}^2+(2{\beta }_1^2+3r_1^2+3q^2+3{\delta }^2)|I{|}^2\\ \le \left|\right|N|{|}_{\infty }^2+((2{\beta }_1^2+3r_1^2+3q^2+3{\delta }^2)|I{|}^2)\\ =\left|\right|N|{|}_{\infty }^2\left[1+\left(\frac{2{\beta }_1^2+3r_1^2+3q^2+3{\delta }^2}{\left|\right|N|{|}_{\infty }^2}\right)|I{|}^2\right].\end{array}$$(4.18)

On taking $\overline{L_2}=\left|\right|N|{|}_{\infty }^2,$ and under the condition $$\left(\frac{2{\beta }_1^2+3r_1^2+3q^2+3{\delta }^2}{\left|\right|N|{|}_{\infty }^2}\right)\le 1.$$

We get $|g_2(t,I(t)){|}^2\le \overline{L_2}(1+|I(t){|}^2).$

Next to show $|g_3(t,A(t)){|}^2\le \overline{L_3}(1+|A(t){|}^2).$

Consider (4.19) $$\begin{array}{c}|g_3(t,A){|}^2={\left|\left({\beta }_2{e}^{-zt}\frac{S A}{N}-r_{2}A-q A-\delta A\right)\right|}^2\\ \le 2{\beta }_2^2{e}^{-2zt}\frac{\left|S{|}^2\right|A{|}^2}{|N{|}^2}+3(r_2^2+q^2+{\delta }^2)|A{|}^2\\ =\left(2{\beta }_2^2{e}^{-2zt}\frac{|S{|}^2}{|N{|}^2}+3(r_2^2+q^2+{\delta }^2)\right)|A{|}^2\\ \le (2{\beta }_2^2+3(r_2^2+q^2+{\delta }^2))|A{|}^2\\ \le |N{|}^2+(2{\beta }_2^2+3r_2^2+3q^2+3{\delta }^2)|A{|}^2\\ \le \left|\right|N|{|}_{\infty }^2+((2{\beta }_2^2+3r_2^2+3q^2+3{\delta }^2)|A{|}^2)\\ =\left|\right|N|{|}_{\infty }^2\left[1+\left(\frac{2{\beta }_2^2+3r_2^2+3q^2+3{\delta }^2}{\left|\right|N|{|}_{\infty }^2}\right)|A{|}^2\right].\end{array}$$(4.19)

On taking $\overline{L_3}=\left|\right|N|{|}_{\infty }^2,$ and under the condition $$\left(\frac{2{\beta }_2^2+3r_2^2+3q^2+3{\delta }^2}{\left|\right|N|{|}_{\infty }^2}\right)\le 1.$$

We get $|g_3(t,A(t)){|}^2\le \overline{L_3}(1+|A(t){|}^2).$

Next to show $|g_4(t,Q(t)){|}^2\le \overline{L_4}(1+|Q(t){|}^2).$

Consider (4.20) $$\begin{array}{c}|g_4(t,Q){|}^2=|q I+q A-(r_3+\delta )Q{|}^2\\ \le 3q^2|I{|}^2+3q^2|A{|}^2+(3r_3^2+3{\delta }^2)|Q{|}^2\\ \le 3q^2\left|\right|I|{|}_{\infty }^2+3q^2|\left|A\right|{|}_{\infty }^2+3(r_3^2+3{\delta }^2)|Q{|}^2\\ =3q^2\left|\right|I|{|}_{\infty }^2+3q^2|\left|A\right|{|}_{\infty }^2\left[1+\frac{(3r_3^2+3{\delta }^2)|Q{|}^2}{3q^2\left|\right|I|{|}_{\infty }^2+3q^2|\left|A\right|{|}_{\infty }^2}\right].\end{array}$$(4.20)

On taking $\overline{L_4}=3q^2\left|\right|I|{|}_{\infty }^2+3q^2|\left|A\right|{|}_{\infty }^2$ and under the condition $$\frac{(3r_3^2+3{\delta }^2)}{3q^2\left|\right|I|{|}_{\infty }^2+3q^2|\left|A\right|{|}_{\infty }^2}\le 1.$$

We get $|g_4(t,Q(t)){|}^2\le \overline{L_4}(1+|Q(t){|}^2).$

Lastly we show $|g_5(t,R(t)){|}^2\le \overline{L_5}(1+|R(t){|}^2).$

Consider (4.21) $$\begin{array}{c}|g_5(t,R){|}^2=|r_{1}I+r_{2}A+r_{3}Q-(\omega +\delta )R{|}^2\\ \le 4r_1^2|I{|}^2+4r_2^2|A{|}^2+4r_3^2|Q{|}^2+4({\omega }^2+{\delta }^2)|R{|}^2\\ \le 4r_1^2\left|\right|I|{|}_{\infty }^2+4r_2^2|\left|A\right|{|}_{\infty }^2+4r_3^2\left|\right|Q|{|}_{\infty }^2+4({\omega }^2+{\delta }^2)|R{|}^2\\ =4r_1^2\left|\right|I|{|}_{\infty }^2+4r_2^2|\left|A\right|{|}_{\infty }^2+4r_3^2\left|\right|Q|{|}_{\infty }^2\left[1+\frac{(4{\omega }^2+4{\delta }^2)|R{|}^2}{4r_1^2\left|\right|I|{|}_{\infty }^2+4r_2^2|\left|A\right|{|}_{\infty }^2+4r_3^2\left|\right|Q|{|}_{\infty }^2}\right].\end{array}$$(4.21)

On taking $\overline{L_5}=4r_1^2\left|\right|I|{|}_{\infty }^2+4r_2^2|\left|A\right|{|}_{\infty }^2+4r_3^2\left|\right|Q|{|}_{\infty }^2$ and under the condition $$\frac{(4{\omega }^2+4{\delta }^2)}{4r_1^2\left|\right|I|{|}_{\infty }^2+4r_2^2|\left|A\right|{|}_{\infty }^2+4r_3^2\left|\right|Q|{|}_{\infty }^2}\le 1.$$

We get $$|g_5(t,R(t){|}^2\le \overline{L_5}(1+|R(t){|}^2).$$

Hence, by using Theorem 1, the existence and uniqueness of the SIAQR model is proved.

5. Derivation of numerical scheme for Caputo fractional SIAQR model using Atangana-Seda method

5.1. Atangana-Seda method with classical derivative

Consider the following Cauchy problem (Atangana & Araz, 2020). (5.1) $$\frac{dy(t)}{dt}=g(t,y(t)),$$(5.1) where g is a non-linear function.

To provide a numerical scheme to solve the above mentioned equation, we convert the above equation as (5.2) $$y(t)-y(0)={\int }_0^{t}g(\tau ,y(\tau ))d\tau .$$(5.2)

At point $t_{m+1}\Delta t,$ we have (5.3) $$y{(t)}_{m+1}-y(0)={\int }_0^{t_{m+1}}g(\tau ,y(\tau ))d\tau .$$(5.3)

At point $t_m\Delta t,$ we have (5.4) $$y{(t)}_m-y(0)={\int }_0^{t_m}g(\tau ,y(\tau ))d\tau .$$(5.4)

On considering the difference of above two equations, we have (5.5) $$y(t_{m+1})-y(t_m)={\int }_{t_m}^{t_{m+1}}g(\tau ,y(\tau )).$$(5.5)

The approximation of $g(t,y(t))$ as the Newton polynomial is given as (5.6) $$\begin{array}{c}{Q}_m(\tau )=g(t_{m-2,y(t_{m-2})})+\frac{g(t_{m-1},y(t_{m-1}))-g(t_{m-2},y(t_{m-2}))}{\Delta t}\left(\tau -t_{m-2}\right)\\ +\frac{g(t_m,y(t_m))-2g(t_{m-1},y(t_{m-1}))+g(t_{m-2},y(t_{m-2}))}{2{(\Delta t)}^2}\\ \times (\tau -t_{m-2})(\tau -t_{m-1}).\end{array}$$(5.6)

Substituting the approximated polynomial in equation (5.5)(5.5) $$y(t_{m+1})-y(t_m)={\int }_{t_m}^{t_{m+1}}g(\tau ,y(\tau )).$$(5.5) , we get (5.7) $$\begin{array}{c}{y}^{m+1}-y^m=g(t_{m-2},y(t_{m-2}))\Delta t+{\int }_{t_m}^{t_{m+1}}\frac{g(t_{m-1},y^{m-1})-g(t_{m-2},y^{m-2})}{\Delta t}\\ \times \left(\tau -t_{m-2}\right)d\tau \\ {\int }_{t_m}^{t_{m+1}}\frac{g(t_m,y^m)-2g(t_{m-1},y^{m-1})+g(t_{m-2},y^{m-2})}{2{(\Delta t)}^2}\\ \times (\tau -t_{m-2})(\tau -t_{m-1})d\tau .\end{array}$$(5.7)

On solving, we get (5.8) $$\begin{array}{c}{y}^{m+1}-y^m=g(t_{m-2},y(t_{m-2}))\Delta t+\frac{g(t_{m-1},y^{m-1})-g(t_{m-2},y^{m-2})}{\Delta t}\\ \times {\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})d\tau \\ \frac{g(t_m,y^m)-2g(t_{m-1},y^{m-1})+g(t_{m-2},y^{m-2})}{2{(\Delta t)}^2}\\ \times {\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})(\tau -t_{m-1})d\tau .\end{array}$$(5.8)

The values of the integrals mentioned in above equation are as follows (5.9) $$\begin{array}{c}{\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})d\tau =\frac{5}{2}(\Delta t).\\ {\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})(\tau -t_{m-1})d\tau =\frac{23}{6}(\Delta t).\end{array}$$(5.9)

Substituting values of integrals in equation (5.8)(5.8) $$\begin{array}{c}{y}^{m+1}-y^m=g(t_{m-2},y(t_{m-2}))\Delta t+\frac{g(t_{m-1},y^{m-1})-g(t_{m-2},y^{m-2})}{\Delta t}\\ \times {\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})d\tau \\ \frac{g(t_m,y^m)-2g(t_{m-1},y^{m-1})+g(t_{m-2},y^{m-2})}{2{(\Delta t)}^2}\\ \times {\int }_{t_m}^{t_{m+1}}(\tau -t_{m-2})(\tau -t_{m-1})d\tau .\end{array}$$(5.8) (5.10) $$\begin{array}{c}{y}^{m+1}=y^m+g(t_{m-2},y(t_{m-2}))\Delta t+\frac{g(t_{m-1},y^{m-1})-g(t_{m-2},y^{m-2})}{\Delta t}\frac{5}{2}(\Delta t)\\ \frac{g(t_m,y^m)-2g(t_{m-1},y^{m-1})+g(t_{m-2},y^{m-2})}{2{(\Delta t)}^2}\frac{23}{6}(\Delta t).\end{array}$$(5.10)

5.2. Atangana-Seda method with Caputo derivative

(5.11) $${}^{\mathcal{C}}{\chi }_t^{\gamma }H(t)=g(H,t)$$(5.11)

Applying fractional integral on the above system using definition 1, we get (5.12) $$H(t)=H(0)+\frac{1}{\Gamma (\gamma )}{\int }_0^{t}g(H,\tau ){(t-\tau )}^{\gamma -1}d\tau .$$(5.12)

At point $t=t_{m+1}$ (5.13) $$H(t_{m+1})=H(0)+\frac{1}{\Gamma (\gamma )}{\int }_0^{t_{m+1}}g(H,\tau ){(t_{m+1}-\tau )}^{\gamma -1}d\tau .$$(5.13)

The above can be written as (5.14) $$H(t_{m+1})=H(0)+\frac{1}{\Gamma (\gamma )}\sum _{r=2}^m{\int }_{t_r}^{t_{r+1}}g(H,\tau ){(t_{m+1}-\tau )}^{\gamma -1}d\tau .$$(5.14)

Using Newton’s polynomial approximation we write $g(H,\tau )$ as (5.15) $$\begin{array}{c}g(H,\tau )\approx g(H^{m-2},t_{m-2})+\frac{g(H^{m-1},t^{m-1})-g(H^{m-2},t^{m-2})}{\Delta t}(\tau -t_{m-2})\\ +\frac{g(H^m,t_m)-2g(H^{m-1},t_{m-1})+g(H^{m-2},t_{m-2})}{2{(\Delta t)}^2}(\tau -t_{m-2})(\tau -t_{m-1}).\end{array}$$(5.15)

Substituting $g(H,\tau )$ approximated through Newton’s polynomial in equation (5.14)(5.14) $$H(t_{m+1})=H(0)+\frac{1}{\Gamma (\gamma )}\sum _{r=2}^m{\int }_{t_r}^{t_{r+1}}g(H,\tau ){(t_{m+1}-\tau )}^{\gamma -1}d\tau .$$(5.14) , we get (5.16) $$\begin{array}{c}{H}^{m+1}=H^0+\frac{1}{\Gamma (\gamma )}\sum _{r=2}^{m}g(H^{r-2},t_{r-2}){\int }_{t_r}^{t_{r+1}}{(t_{m+1}-\tau )}^{\gamma -1}d\tau \\ +\frac{1}{\Gamma (\gamma )}\sum _{r=2}^m\frac{g(H^{r-1},t^{r-1})-g(H^{r-2},t^{r-2})}{\Delta t}{\int }_{t_r}^{t_{r+1}}(\tau -t_{r-2}){(t_{m+1}-\tau )}^{\gamma -1}d\tau \\ +\frac{1}{\Gamma (\gamma )}\sum _{r=2}^m\frac{g(H^m,t_m)-2g(H^{m-1},t_{m-1})+g(H^{m-2},t_{m-2})}{2{(\Delta t)}^2}\\ \times {\int }_{t_r}^{t_{r+1}}(\tau -t_{r-2})(\tau -t_{r-1}){(t_{m+1}-\tau )}^{\gamma -1}d\tau .\end{array}$$(5.16)

Lastly, on solving the integrals in the above expression, we finally get (5.17) $$\begin{array}{c}{H}^{m+1}=H^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^{m}g(H^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g(H^{r-1},t^{r-1})-g(H^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g(H^r,t_r)-2g(H^{r-1},t_{r-1})+g(H^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.17)

To find the numerical solution, we make use of the Atangana-Seda (Atangana & Araz, 2018, 2020; Gnitchogna & Atangana, 2018; Toufik & Atangana, 2017).

Now, using the aforementioned numerical scheme, we finally get the numerical solution for the Caputo fractional SIAQR model as (5.18) $$\begin{array}{c}{S}^{m+1}=S^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^m{g}_1(S^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g_1(S^{r-1},t^{r-1})-g_1(S^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g_1(S^r,t_r)-2g_1(S^{r-1},t_{r-1})+g_1(S^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.18) where (5.19) $$\begin{array}{c}{g}_1(t,S)=\delta N(t) -{\beta }_1 e^{-z\mathrm{ }t} \frac{S(t) I(t)}{N}-{\beta }_2 e^{-z\mathrm{ }t} \frac{S(t) A(t)}{N} +\omega R(t) -\delta S(t).\\ I^{m+1}=I^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^m{g}_2(I^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g_2(I^{r-1},t^{r-1})-g_2(I^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g_2(I^r,t_r)-2g_2(I^{r-1},t_{r-1})+g_2(I^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.19) where (5.20) $$\begin{array}{c}{g}_2(t,I)={\beta }_1 e^{-z\mathrm{ }t} \frac{S(t) I(t)}{N}-r_1 I(t)-q I(t)-\delta I(t).\\ A^{m+1}=A^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^m{g}_3(A^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g_3(A^{r-1},t^{r-1})-g_3(A^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g_3(A^r,t_r)-2g_3(A^{r-1},t_{r-1})+g_3(A^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.20) where (5.21) $$\begin{array}{c}{g}_3(t,A)={\beta }_2 e^{-z\mathrm{ }t} \frac{S(t) A(t)}{N}-r_2 A(t)-q A(t)-\delta A(t).\\ Q^{m+1}=Q^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^m{g}_4(Q^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g_4(Q^{r-1},t^{r-1})-g_4(Q^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g_4(Q^r,t_r)-2g_4(Q^{r-1},t_{r-1})+g_4(Q^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.21) where (5.22) $$\begin{array}{c}{g}_4(t,Q)=q I(t)+q A(t)-r_3 Q(t)-\delta Q(t).\\ R^{m+1}=R^0+\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +1)}\sum _{r=2}^m{g}_5(R^{r-2},t_{r-2})[{(m-r+1)}^{\gamma }-{(m-r)}^{\gamma }]\\ +\frac{{(\Delta t)}^{\gamma }}{\Gamma (\gamma +2)}\sum _{r=2}^m[g_5(Q^{r-1},t^{r-1})-g_5(Q^{r-2},t^{r-2})]\\ \times [{(m-r+1)}^{\gamma }(m-r+3+2\gamma )-{(m-r)}^{\gamma }(m-r+3+3\gamma )]\\ +\frac{{(\Delta t)}^{\gamma }}{2\Gamma (\gamma +3)}\sum _{r=2}^m[g_5(Q^r,t_r)-2g_5(Q^{r-1},t_{r-1})+g_5(Q^{r-2},t_{r-2})]\\ \times [{(m-r+1)}^{\gamma }(2{(m-r)}^2+(3\gamma +10)(m-r)+2{\gamma }^2+9\gamma +12)\\ -{(m-r)}^{\gamma }(2{(m-r)}^2+(5\gamma +10)(m-r)+6{\gamma }^2+18\gamma +12)].\end{array}$$(5.22) where $$g_5(t,R)=r_1 I(t)+r_2 A(t)+r_3 Q(t)-\omega R(t)-\delta R(t).$$

5.3. Numerical simulations

Numerical simulations for SIAQR model with ${\beta }_1=0.5$, $r_2=0.1$, $r_3=0.175$, $q\text{\hspace{0.17em}=\hspace{0.17em}0.01}$, $\delta =0.001$, $S(0)=1000$, $I(0)=20$, $A(0)=5$, $Q(0)=0$, $R(0)=0$.

• In Figure 1, we see the numerical simulations for z = 0.001 and in fractional order $\gamma =0.5.$ We see the effect of adding the time-dependent parameter z. The figure shows that susceptible, infected symptomatically, and infected asymptomatically decrease, and recovered individuals increase.

• In Figure 2, we have presented the numerical simulation at $\gamma =0.7.$ We see that as the order of fractional derivative increases from $\gamma =0.5$ to $\gamma =0.7,$ the decrease in the number of individuals who are infected symptomatically decreases well and the number of quarantined individuals increases.

• In Figure 3, we have presented the numerical simulation at $\gamma =0.899$ and z = 0.001. We see that as the order of fractional derivative increases from $\gamma =0.7$ to $\gamma =0.899,$ the decrease in the decrease in the number of individuals who are symptomatically infected is much greater compared to $\gamma =0.5$ and $\gamma =\mathrm{0.7.}.$ Furthermore, we can see that the number of asymptomatically infected individuals is also decreased.

• In Figure 4, we have presented the numerical simulation at $\gamma =0.5$ and z = 0.01. We see that as the value of z increases from 0.001 to 0.01, the decrease in the decrease in the number of individuals who are infected symptomatically and asymptomatically is much greater. This shows the effect of time dependency on the infection rate due to the parameter z, which determines the change in the time of infection.

• In Figure 5, we have presented the numerical simulation at $\gamma =0.7$ and z = 0.01. We see that as the value of z increases from 0.001 to 0.01. The number of recovered individuals increases well and the number of asymptomatically infected individuals is quite small.

• In Figure 6, we have presented the numerical simulation at $\gamma =0.899$ and z = 0.01. We see that as the value of z increases from 0.001 to 0.01. The number of recovered individuals is much more than the number of individuals who are symptomatically and asymptomatically infected.

• In Figure 7, we have presented the numerical simulation at $\gamma =0.899$ and z = 0.1. We see that due to the effect of the parameter z (change in infection time), there is a sharp increase in the number of quarantined individuals, and hence the recovery rate is also much higher.

Figure 1. Numerical simulations for SIAQR model for z = 0.001 and $\gamma =0.5.$

Figure 2. Numerical simulations for SIAQR model for z = 0.001 and $\gamma =0.7.$

Figure 3. Numerical simulations for SIAQR model for z = 0.001 and $\gamma =0.899.$

Figure 4. Numerical simulations for SIAQR model for z = 0.01 and $\gamma =0.5.$

Figure 5. Numerical simulations for SIAQR model for z = 0.01 and $\gamma =0.7.$

Figure 6. Numerical simulations for SIAQR model for z = 0.01 and $\gamma =0.899.$

Figure 7. Numerical simulations for SIAQR model for z = 0.1 and $\gamma =0.899.$

6. Conclusion

In this work, we have formulated the SIAQR model for Covid pestilence with a time-dependent infection rate. The SIAQR model is then extended to the fractional mathematical model within the framework of the Caputo derivative operator. The existence and uniqueness of the solution of the SIAQR system is carried out using the Lipschitz condition and the linear growth condition. The numerical scheme for the Caputo fractional SIAQR model is derived using the newly proposed Atangana-Seda method. Numerical simulations are performed for the SIAQR model with a time-dependent infection rate for different fractional order derivatives. Furthermore, we see that as the value of z increases, the number of infected individuals decreases. This is because the rate of infection is lower as the value of z increases. We also see that as the order of fractional derivatives increases, the number of infected populations decreases, and the number of recovered populations increases. Hence, the fractional derivatives provide us with a better understanding of the model.

Acknowledgements

Authors are thankful to reviewers for their constructive comments to improve this article. The researchers would like to acknowledge the Deanship of Scientific Research, Taif University for funding this work.