Theory and Modern Applications

# Crowding effects on the dynamics of COVID-19 mathematical model

## Abstract

A disastrous coronavirus, which infects a normal person through droplets of infected person, has a route that is usually by mouth, eyes, nose or hands. These contact routes make it very dangerous as no one can get rid of it. The significant factor of increasing trend in COVID19 cases is the crowding factor, which we named “crowding effects”. Modeling of this effect is highly necessary as it will help to predict the possible impact on the overall population. The nonlinear incidence rate is the best approach to modeling this effect. At the first step, the model is formulated by using a nonlinear incidence rate with inclusion of the crowding effect, then its positivity and proposed boundedness will be addressed leading to model dynamics using the reproductive number. Then to get the graphical results a nonstandard finite difference (NSFD) scheme and fourth order Runge–Kutta (RK4) method are applied.

## Introduction

Coronaviruses cause a common group of infections that results in common cold type symptoms being one of the old classes of viruses in human history, but COVID19 became most disastrous, resulting in the highest death tolls in its track records. In history it became deadly in the shape of Serious Intense Respiratory Conditions (SARS) and Middle East Respiratory Disorder (MERS). The Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2) opened new doors for researchers and on February 12, 2020, it was formally named Novel Coronavirus disease 2019 (COVID-19) by the World Health Organization (WHO). The spread mechanism is very simple: when somebody infected sneezes or coughs, the droplets enter nearby people and soon enter their bodies through contact routes like mouth, nose, and later their lungs, and start damaging their respiratory system as it is operated through the lungs. Later it was claimed by researchers that contacting infected surfaces also causes transmission of the virus. There are many researchers stating the possible lifetime of virus on different surfaces; WHO also gave many guidelines on that. The National Institute of Health in our country states that the virus that causes coronavirus illness 2019 (COVID-19) is stable for many hours to days in aerosols and on surfaces. So far researchers are giving different time frames for the life of the virus on different surfaces. The patients are categorized in three types in start by WHO. (i) Suspect case: These are the patients with acute respiratory illness and who have been in contact with a COVID19 infected person. (ii) Exposed case: These are the patients who are infected but so far showing no symptoms of COVID19. (iii) Confirmed case: These are infectious people with authentic laboratory assertion with all the symptoms, and they are put in isolation to protect other people and the community.

Measuring the spread of a disease is a significant factor in epidemic models, which is determined by incidence rates. Generally, it is thought of as the newly infected rate per unit of time. For this purpose, $$\beta SI$$ with an extension of $$\beta SI/(1+\alpha I)$$ are used. The saturated nonlinear incidence $$Sf(I)$$ approach was used by Capasso and Serio  in the case of a cholera epidemic model. It is commonly observed in most communicable diseases that the acquired immunity may disappear after some time such as with pertussis, influenza, malaria and cholera (time may vary from disease to disease) and this puts individuals at risk again; see [2, 3]. It is also reported that in some situations the recovered patients may get mixed with a susceptible group of people with the belief of having a transient antibody [4, 5]. The SIRS model designed by Chen, having a standard incident, disease related standard incidence and death, and we have transference from the infected group to the susceptible class. So far, many research domains are in process to address the current pandemic using different mathematical models like differential equations (in integer and fractional form) (see  and the references therein).

Getting motivation from the above work, in this work a new approach is used to address COVID19 by taking the “crowding impact” in considerations. The theme behind this work is how crowding of infected individuals affects the susceptible class or population. In mathematical models this effect is addressed by a nonlinear incidence rate.

The structure of paper is as follows: in Sect. 2, we will show the proposed model, Sect. 3 is the qualitative analysis whereas Sect. 4 is all about the numerical solution of the proposed model. In the last section our conclusion and future work will be presented.

## Model formulation

In order to illustrate the crowding effects in a COVID-19 mathematical model, we proposed a model with recruitment rate μN for the susceptible individual and nonlinear incidence rate $$\frac{\beta SI}{1+\alpha I}$$. The complete dynamics is given by

\begin{aligned} \textstyle\begin{cases} \frac{dS}{dt} = \mu N-\frac{\beta SI}{1+\alpha I}-\mu S, \\ \frac{dI}{dt} =\frac{\beta SI}{1+\alpha I}-(\gamma +\mu )I, \\ \frac{dR}{dt} = \gamma I-\mu R, \end{cases}\displaystyle \end{aligned}
(1)

where the infection force of the disease is expressed by βI, the crowding effect by $$\frac{1}{1+\alpha I}$$, the susceptible population by $$S(t)$$, the infectious population by $$I(t)$$ and the recovered population by $$R(t)$$ at time t. Furthermore, μ is the death rate, γ is the recovery rate and β is the transmission coefficient. $$N(t)$$ is the total constant population,

\begin{aligned} N(t)=S(t)+I(t)+R(t), \end{aligned}
(2)

with

\begin{aligned} N(0)=S(0)+I(0)+R(0). \end{aligned}
(3)

Furthermore, from

\begin{aligned} S(0)\geq 0,\qquad I(0)\geq 0,\qquad R(0) \geq 0, \end{aligned}

we have

\begin{aligned} S(t)\geq 0,\qquad I(t)\geq 0,\qquad R(t) \geq 0. \end{aligned}
(4)

So the solution possesses the property of positivity.

## Dynamics of the model

### Theorem 1

(Positivity)

Suppose that the related solution of the initial data $$(S(0), I(0), R(0))\in \mathbf{R}_{+}^{3}$$ is $$(S(t), I(t), R(t) )$$. Then for model (1) the positively invariant set is $$\mathbf{R}^{3}_{+}$$.

### Proof

By setting $$\Lambda =\mu N$$ and $$\lambda _{1}=\frac{\beta I}{1+\alpha I}$$, the first equation of system (1) implies that

\begin{aligned} \frac{dS}{dt}=\Lambda -\lambda _{1}S-\mu S. \end{aligned}

Suppose that the solution exists of system (1) for a certain interval $$J\in [0;+\infty [$$, then the above equation can be solved, for all $$t\in J$$, as

\begin{aligned} & \frac{dS}{dt}+(\lambda _{1}+\mu )S= \Lambda, \\ &\frac{d}{dt} \bigl(S(t)e^{\mu t+\int _{0}^{t}\lambda _{1}(s)\,ds} \bigr) \geq \Lambda e^{\mu t+\int _{0}^{t}\lambda _{1}(s)\,ds}, \end{aligned}

which implies that

\begin{aligned} & S(t)e^{\mu t+\int _{0}^{t}\lambda _{1}(s)\,ds}-S(0)\geq \int _{0}^{t} \Lambda e^{\mu t+\int _{0}^{w}\lambda _{1}(u)\,du}\,dw, \\ &S(t)e^{\mu t+\int _{0}^{t}\lambda _{1}(s)\,ds}\geq S(0)+ \int _{0}^{t} \Lambda e^{\mu w+\int _{0}^{w}\lambda _{1}(u)\,du}\,dw, \\ &S(t)\geq S(0)e^{-\mu t-\int _{0}^{t}\lambda _{1}(s)\,ds}+e^{-\mu t- \int _{0}^{t}\lambda _{1}(s)\,ds}\times \int _{0}^{t}\Lambda e^{\mu t+ \int _{0}^{w}\lambda _{1}(u)\,du}\,dw>0. \end{aligned}

Hence, $$\forall t\in J$$, S(t) is positive. Now the second equation of the model (1) implicitly shows that

\begin{aligned} \frac{dI}{dt} =\frac{\beta SI}{1+\alpha I}-(\gamma +\mu )I\geq -( \gamma +\mu )I, \end{aligned}

or

\begin{aligned} \frac{dI}{dt}\geq -(\gamma +\mu )I, \end{aligned}

which can be written as

\begin{aligned} \frac{dI}{I}\geq -(\gamma +\mu )\,dt\quad (I\neq 0). \end{aligned}

By integrating, we get

\begin{aligned} &\ln {I}\geq -(\gamma +\mu )t+C, \\ &I\geq e^{-(\gamma +\mu )t+C}, \\ &I\geq C_{1}e^{-(\gamma +\mu )t}, \end{aligned}

at $$t=0$$,

\begin{aligned} I\geq I(0)e^{-(\gamma +\mu )0}\geq 0. \end{aligned}

Therefor, for all values of t, $$I(t)$$ is positive. Similarly, the last equation of system (1) implies that

\begin{aligned} \frac{dR}{dt} =\gamma I-\mu R\geq -\mu R, \end{aligned}

or

\begin{aligned} \frac{dR}{dt}\geq -\mu R, \end{aligned}

which can be written as

\begin{aligned} \frac{dR}{R}\geq -\mu dt \quad(R\neq 0). \end{aligned}

By integrating, we get

\begin{aligned} &\ln {R}\geq -\mu t+K, \\ &R\geq e^{-\mu t+K}, \\ &R\geq K_{1}e^{-\mu t}, \end{aligned}

at $$t=0$$,

\begin{aligned} R\geq R(0)e^{-\mu 0}\geq 0. \end{aligned}

Hence, $$R(t)$$ also is positive in the given interval. □

### Existence and uniqueness of the solution

The first-order ODE in general form is

\begin{aligned} \acute{z}=g(t,z),\qquad z(t_{0})=z_{0}. \end{aligned}
(5)

With the help of the theorem below, we can establish the existence and uniqueness of the solution for the considered model.

### Theorem 2

(Uniqueness of solution)

Let us use D to denote the domain

\begin{aligned} \vert t-t_{0} \vert \leq a,\qquad \vert z-z_{0} \vert \leq b,\quad z=(z_{1},z_{2}, \ldots,z_{n}), z_{0}=(z_{10},z_{20}, \ldots,z_{n0}), \end{aligned}
(6)

and assume that the Lipschitz condition is satisfied by $$h(t,z)$$:

\begin{aligned} \bigl\Vert h (t, z_{1}) - h(t, z_{2}) \bigr\Vert \leq c \Vert z_{1}- z_{2} \Vert , \end{aligned}
(7)

and the two pairs $$(t, z_{1})$$ and $$(t, z_{2})$$ are in D, where c is a positive constant. Then there exists a constant $$\delta > 0$$ such that for the interval $$|t-t_{0}|\leq \delta$$ there exists a unique continuous vector solution $$z(t)$$ of the system (5). It should be noticed that condition (6) is satisfied with

\begin{aligned} \frac{\partial h_{i}}{\partial z_{j}},\quad i,j=1,2,\ldots,n, \end{aligned}
(8)

in the domain D, being continuous and bounded.

### Lemma 1

If the continuous partial derivative of $$h (t, z)$$ (i.e., $$\frac{\partial h_{i}}{\partial z_{j}}$$) exists for a bounded closed convex domain , then, for , it satisfies a Lipschitz condition. We are interested in the domain

\begin{aligned} 1\leq \epsilon \leq \Re. \end{aligned}
(9)

Hence, a solution in the form of condition (10) is searched:

\begin{aligned} 0< \Re < \infty. \end{aligned}
(10)

Now the existence theorem can be proved as follows.

### Theorem 3

Assume D represents the domain of (6) in such a manner that (7) and (8) hold. Then the bounded solution in domain D of (1) exists.

### Proof

Let

\begin{aligned} & h_{1}=\mu N-\frac{\beta SI}{1+\alpha I}-\mu S, \end{aligned}
(11)
\begin{aligned} & h_{2}=\frac{\beta SI}{1+\alpha I}-(\gamma +\mu )I, \end{aligned}
(12)
\begin{aligned} &h_{3}=\gamma I-\mu R. \end{aligned}
(13)

To show that $$\frac{\partial h_{i}}{\partial z_{j}}$$ $$j = 1, 2, 3$$ are continuous and bounded, the following partial derivatives for the proposed model are performed. By taking the partial derivative of Eq. (11) we have

\begin{aligned} & \frac{\partial h_{1}}{\partial S}=- \frac{\beta I}{1+\alpha I}-\mu, \\ &\biggl\vert \frac{\partial h_{1}}{\partial S} \biggr\vert = \biggl\vert - \frac{\beta I}{1+\alpha I}-\mu \biggr\vert < \infty, \\ & \frac{\partial h_{1}}{\partial I}=- \frac{\beta S}{(1+\alpha I)^{2}} \\ &\biggl\vert \frac{\partial h_{1}}{\partial I} \biggr\vert = \biggl\vert - \frac{\beta S}{(1+\alpha I)^{2}} \biggr\vert < \infty, \end{aligned}

and

\begin{aligned} &\frac{\partial h_{1}}{\partial R}=0, \\ &\biggl\vert \frac{\partial h_{1}}{\partial R} \biggr\vert = \vert 0 \vert < \infty. \end{aligned}

For class I, from Eq. (12)

\begin{aligned} & \frac{\partial h_{2}}{\partial S}=\frac{\beta I}{1+\alpha I} , \\ &\biggl\vert \frac{\partial h_{2}}{\partial S} \biggr\vert =\biggl\vert \frac{\beta I}{1+\alpha I} \biggr\vert < \infty, \\ & \frac{\partial h_{2}}{\partial I}=\frac{\beta S}{(1+\alpha I)^{2}}-( \gamma +\mu ), \\ &\biggl\vert \frac{\partial h_{2}}{\partial I} \biggr\vert = \biggl\vert \frac{\beta S}{(1+\alpha I)^{2}}-(\gamma +\mu ) \biggr\vert < \infty, \end{aligned}

and

\begin{aligned} & \frac{\partial h_{2}}{\partial R}=0, \\ &\biggl\vert \frac{\partial h_{2}}{\partial R} \biggr\vert = \vert 0 \vert < \infty. \end{aligned}

Similarly, for class R, from Eq. (13)

\begin{aligned} & \frac{\partial h_{3}}{\partial S}=0, \\ &\biggl\vert \frac{\partial h_{3}}{\partial S} \biggr\vert = \vert 0 \vert < \infty, \\ & \frac{\partial h_{3}}{\partial I}=\gamma, \\ &\biggl\vert \frac{\partial h_{3}}{\partial I} \biggr\vert = \vert \gamma \vert < \infty, \end{aligned}

and

\begin{aligned} & \frac{\partial h_{3}}{\partial R}=-\mu, \\ &\biggl\vert \frac{\partial h_{3}}{\partial R} \biggr\vert =\vert -\mu \vert < \infty, \end{aligned}

therefore, it can be concluded that all partial derivatives are bounded in the considered domain and are continuous. Hence, from Theorem 2, it is proved that there exists a unique solution of system (1) in D. □

### Equilibria and the reproduction number

For finding the equilibrium points of model (1), we take the algebraic system

\begin{aligned} \begin{aligned} & \mu N-\frac{\beta SI}{1+\alpha I}-\mu S = 0, \\ &\frac{\beta SI}{1+\alpha I}-(\gamma +\mu )I = 0, \\ &\gamma I-\mu R = 0.\end{aligned} \end{aligned}
(14)

Two solutions of the system (1) are obtained by some algebraic manipulations, one is $$D_{0} = (N, 0, 0)$$, a disease free equilibrium (DFE) point and the second will be discussed after computing the reproduction number of the model (1). The reproductive number is computed with the help of the next generation matrix approach presented by van den Driessche and Watmough .

Let $$x = (I, S)$$ and rewrite the model (1) for the susceptible and infected classes in the general form

\begin{aligned} \frac{dx}{dt} = \mathcal{F}(x) - \mathcal{V}(x), \end{aligned}
(15)

where

\begin{aligned} \mathcal{F}(x) = \begin{pmatrix} \frac{\beta SI}{1+\alpha I} \\ 0 \end{pmatrix}, \qquad\mathcal{V}(x) = \begin{pmatrix} (\gamma +\mu )I \\ -\mu N+\frac{\beta SI}{1+\alpha I}+\mu S \end{pmatrix}. \end{aligned}
(16)

Now the Jacobian of $$\mathcal{F}(x)$$ and $$\mathcal{V}(x)$$ of the disease free equilibrium point is

\begin{aligned} F = \begin{pmatrix} \beta N &0 \\ 0 &0 \end{pmatrix},\qquad V = \begin{pmatrix} \gamma +\mu & 0 \\ \beta N&\mu \end{pmatrix}, \end{aligned}
(17)

and further by using the idea of van den Driessche and Watmough , the reproduction number of the model (1) is follows:

\begin{aligned} \mathcal{R}_{0} = \rho \bigl(FV^{-1}\bigr) = \frac{\beta N}{\mu +\gamma }. \end{aligned}
(18)

### Theorem 4

For system (1), there exists a unique positive endemic equilibrium point $$D_{*}$$, if $$\mathcal{R}_{0}>1$$.

### Proof

By some algebraic manipulations, the second solution of the system (14) yields

\begin{aligned} &S_{*}=\frac{\mu \alpha N+\mu +\gamma }{\beta +\alpha \mu }, \\ &I_{*}=\frac{\mu (\mathcal{R}_{0} -1)}{\beta +\alpha \mu }, \\ &R_{*}=\frac{\gamma (\mathcal{R}_{0} -1)}{\beta +\alpha \mu }. \end{aligned}

It is clear from the values of $$I_{*}$$ and $$R_{*}$$ that is there exists a unique positive endemic equilibrium point $$D_{*}$$, if $$\mathcal{R}_{0}>1$$. □

### Theorem 5

The system (1) is locally stable related to the virus free equilibrium point $$E_{0}$$, if $$\mathcal{R}_{0}<1$$ and unstable if $$\mathcal{R}_{0}>1$$.

### Proof

For local stability the Jacobian of system (1) is

\begin{aligned} J= \begin{pmatrix} -\mu - \frac{\beta I}{1+\alpha I}& -\frac{\beta S}{(1+\alpha I)^{2}} & 0 \\ \frac{\beta I}{1+\alpha I} & \frac{\beta S}{(1+\alpha I)^{2}}-(\mu + \gamma )& 0 \\ 0 & \gamma & -\mu \end{pmatrix}. \end{aligned}
(19)

At $$E_{0}$$, the Jacobian becomes

\begin{aligned} J(E_{0})= \begin{pmatrix} -\mu & -\beta N & 0 \\ 0 & \beta N -(\mu +\gamma )& 0 \\ 0 & \gamma & -\mu \end{pmatrix}, \end{aligned}
(20)

from which it follows that the eigenvalues are $$\lambda _{1}= -\mu <0, \lambda _{3}= -\mu <0$$ and $$\lambda _{2}=\beta N -(\mu +\gamma )$$, implying that $$\lambda _{2}<0$$, if $$\mathcal{R}_{0}<1$$. So the system (1) is locally stable for $$\mathcal{R}_{0}<1$$ and unstable for $$\mathcal{R}_{0}>1$$. The proof is complete. □

### Theorem 6

If $$\mathcal{R}_{0}<1$$, then the DFE point of the system (1) is globally stable.

### Proof

For the proof of this theorem, first we construct the Lyapunov function L:

\begin{aligned} L(I)=\ln {\frac{I}{I_{0}}}. \end{aligned}
(21)

Differentiating Eq. (21) with respect to time, we have

\begin{aligned} \begin{aligned} \frac{d}{dt}\bigl(L(I)\bigr)={}&\frac{1}{I}\frac{dI}{dt}, \\ \frac{d}{dt}\bigl(L(I)\bigr)={}&\frac{1}{I} \biggl( \frac{\beta SI}{1+\alpha I}-(\gamma +\mu )I \biggr) \\ ={}&\frac{\beta S}{1+\alpha I}-(\gamma +\mu ) \\ ={}&(\gamma +\mu ) \biggl(\frac{\beta S}{(\gamma +\mu )(1+\alpha I)}-1 \biggr) \\ \leq {}&(\gamma +\mu ) \biggl(\frac{\beta S}{\gamma +\mu }-1 \biggr) \\ ={}&(\gamma +\mu ) (\mathcal{R}_{0}-1 ) \\ \leq{} &0 \quad\text{for } \mathcal{R}_{0}< 1. \end{aligned} \end{aligned}
(22)

Therefore, if $$\mathcal{R}_{0}<1$$, then $$\frac{d}{dt}(L(I))<0$$, which implies that, for $$\mathcal{R}_{0}<1$$, the DFE point of the system (1) is globally stable. □

### Theorem 7

For $$\mathcal{R}_{0}>1$$, the system (1) at the positive endemic equilibrium point $$E_{*}$$ is locally stable.

### Proof

The Jacobian matrix of system (1) is

\begin{aligned} J= \begin{pmatrix} -\mu - \frac{\beta I}{1+\alpha I}& -\frac{\beta S}{(1+\alpha I)^{2}} & 0 \\ \frac{\beta I}{1+\alpha I} & \frac{\beta S}{(1+\alpha I)^{2}}-(\mu + \gamma )& 0 \\ 0 & \gamma & -\mu \end{pmatrix}. \end{aligned}
(23)

At $$E_{*}$$, the Jacobian becomes

\begin{aligned} J(E_{*})= \begin{pmatrix} -\mu - \frac{\beta I_{*}}{1+\alpha I_{*}}& - \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}} & 0 \\ \frac{\beta I_{*}}{1+\alpha I_{*}} & \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}}-(\mu +\gamma )& 0 \\ 0 & \gamma & -\mu \end{pmatrix}, \end{aligned}
(24)

which yields one eigenvalue $$\lambda =-\mu$$ and the characteristic equation

\begin{aligned} &\lambda ^{2}+ \biggl(\mu +\frac{\beta I_{*}}{1+\alpha I_{*}}- \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}}+(\mu +\gamma ) \biggr)\lambda \\ &\quad {}+ \biggl(\mu +\frac{\beta I_{*}}{1+\alpha I_{*}} \biggr) \biggl( \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}}+(\mu +\gamma ) \biggr)+ \biggl( \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}} \biggr) \biggl( \frac{\beta I_{*}}{1+\alpha I_{*}} \biggr)=0. \end{aligned}
(25)

It is clear, for $$\mathcal{R}_{0}>1$$, that

\begin{aligned} &\biggl(\mu +\frac{\beta I_{*}}{1+\alpha I_{*}}- \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}}+(\mu +\gamma ) \biggr) \\ &\quad= \biggl(\mu +\frac{\beta I_{*}}{1+\alpha I_{*}}- \frac{\gamma +\mu }{(1+\alpha I_{*})}+(\mu +\gamma ) \biggr)>0 \end{aligned}

and

\begin{aligned} \biggl(\mu +\frac{\beta I_{*}}{1+\alpha I_{*}} \biggr) \biggl( \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}}+(\mu +\gamma ) \biggr)+ \biggl( \frac{\beta S_{*}}{(1+\alpha I_{*})^{2}} \biggr) \biggl( \frac{\beta I_{*}}{1+\alpha I_{*}} \biggr)>0. \end{aligned}

Hence, the system (1) is locally stable at $$E_{*}$$ for $$\mathcal{R}_{0}>1$$. The proof is complete. □

### Remark 1

Although the stability analysis of $$E_{*}$$ is an interesting and separate mathematical problem, while for prevention of the disease one is to find an effective strategy, the main focus of this work is on the specific condition $$\mathcal{R}_{0} < 1$$.

### $$\mathcal{R}_{0}$$ sensitivity analysis

From Theorem 6 it follows that we can control the parameters such that $$\mathcal{R}_{0} < 1$$. This leads to the best strategy to prevent and restrain the disease. In detail, when $$\mathcal{R}_{0} < 1$$, then

\begin{aligned} \lim_{t \to \infty } S(t) = N, \qquad \lim_{t \to \infty } I(t) = \lim _{t \to \infty } R(t) = 0, \end{aligned}

which shows that the spreading speed of the coronavirus can be reduced and prevented in the future. Hence, a sensitivity analysis of $$\mathcal{R}_{0}$$ is carried out to select the influential parameters to control the rapidly spreading current pandemic.

It is easy to verify that

\begin{aligned} \begin{aligned} &\frac{\partial \mathcal{R}_{0}}{\partial \beta } = \frac{N}{\mu +\gamma } > 0, \\ &\frac{\partial \mathcal{R}_{0}}{\partial \mu } = - \frac{\beta N}{(\mu +\gamma )^{2}}< 0, \\ &\frac{\partial \mathcal{R}_{0}}{\partial \gamma } = - \frac{\beta N}{(\mu +\gamma )^{2}} < 0. \end{aligned} \end{aligned}
(26)

Equation (26) can be used to obtain different parameters in such a way that $$\mathcal{R}_{0}$$ remains less than one. Hence, necessary actions can be taken on the basis of Eq. (26) to reduce the speed of the coronavirus.

## Numerical method and results

In order to get the numerical output from model (1), the NSFD method is utilized. The solution via the NSFD method is obtained via an iteration process [27, 28]. Assume the nonstandard ODE

\begin{aligned} w'_{k}=f[t,w_{1},w_{2}, \ldots,w_{n}], \end{aligned}

where $$k=1,2,,\ldots,n$$, then by the NSFD method

\begin{aligned} & w'_{1}=\frac{w_{1,k+1}-w_{1,k}}{h}, \\ &w'_{2}=\frac{w_{2,k+1}-w_{2,k}}{h}, \\ &\cdots \\ &w'_{n}=\frac{w_{n,k+1}-w_{n,k}}{h}. \end{aligned}

Now, for a numerical solution of system (1) using the NSFD method, it gives the following results:

\begin{aligned} &S_{k+1}=\frac{h\mu N+S_{k}}{1+h\mu +h\beta I_{k}/(1+\alpha I_{k})}, \end{aligned}
(27)
\begin{aligned} &I_{k+1}= \frac{h\beta I_{k}S_{k+1}/(1+\alpha I_{k})+I_{k}}{1+h(\mu +\gamma )}, \end{aligned}
(28)

and

\begin{aligned} R_{k+1}=\frac{h\gamma I_{k+1}+R_{k}}{1+h\mu }. \end{aligned}
(29)

Figures 13 show the solutions for $$S(t)$$, $$I(t)$$ and $$R(t)$$ obtained by NSFD, RK4, and ode45 for $$\mathcal{R}_{0}<1$$, when the contact rate is chosen in a small range, then the spread of the current coronavirus disease may be controlled. From Fig. 1, we see how the disease controls which results are leading to an increase of the susceptible class. On the other hand as the susceptibility is increasing the infection goes to extinction reaching stability in Fig. 2. The decrease in the infection class yields increase an in the recovered class as shown in Fig. 3. For the numerical solutions, we consider the initial values $$S=40, I=20, R=10$$ and for the parameters values from , for the remaining features of the model one can use the real-life data of some specific country or of the whole world.

## Conclusion

This work presented the crowding effects of infective individuals over the susceptible population specially for the current pandemic. The crowding effect is described by a nonlinear incidence rate in the mathematical model. In this work, the formulation of the model is presented keeping in mind the crowding effect, which is in a large range of interaction of the infected population with the susceptible population. The dynamics of the model is presented based on the reproductive number and one showed the local and global stability of the proposed model. For a numerical solution we used the nonstandard finite difference (NSFD) scheme and the fourth-order Runge–Kutta (RK4) method and the obtained results are shown graphically.

## Availability of data and materials

The authors confirm that the data supporting the findings of this study are available in the article.

## References

1. Capasso, V., Serio, G.: A generalization of the Kermack–Mckendrick deterministic epidemic model. Math. Biosci. 42, 43–61 (1978)

2. Li, J., Yang, Y., Xiao, Y., Liu, S.: A class of Lyapunov functions and the global stability of some epidemic models with nonlinear incidence. J. Appl. Anal. Comput. 6, 38–46 (2016)

3. Chen, J.: An SIRS epidemic model. Appl. Math. J. Chin. Univ. Ser. B 19, 101–108 (2004)

4. Mena-Lorca, J., Hethcote, H.W.: Dynamic models of infectious diseases as regulators of population size. J. Math. Biol. 30, 693–716 (1992)

5. Hethcote, H.W.: Qualitative analysis of communicable disease models. Math. Biosci. 28, 334–356 (1976)

6. Abdo, M.S., Hanan, K.S., Satish, A.W., Pancha, K.: On a comprehensive model of the novel coronavirus (COVID-19) under Mittag-Leffler derivative. Chaos Solitons Fractals 135, 109867 (2020)

7. Bekiros, S., Kouloumpou, D.: SBDiEM: a new mathematical model of infectious disease dynamics. Chaos Solitons Fractals 136, 109828 (2020)

8. Bocharov, G., Volpert, V., Ludewig, B., Meyerhans, A.: Mathematical Immunology of Virus Infections. Springer, Berlin (2018)

9. Brauer, F., van den Driessche, P., Wu, J.: Mathematical Epidemiology. Springer, Berlin (2008)

10. Brauer, F.: Mathematical epidemiology: past, present, and future. Infect. Dis. Model. 2, 113–127 (2017)

11. Cakan, S.: Dynamic analysis of a mathematical model with health care capacity for COVID-19 pandemic. Chaos Solitons Fractals 139, 110033 (2020)

12. Din, R., Shah, K., Ahmad, I., Abdeljawad, T.: Study of transmission dynamics of novel COVID-19 by using mathematical model. Adv. Differ. Equ. 2020, Article ID 323 (2020)

13. Higazy, M.: Novel fractional order SIDARTHE mathematical model of COVID-19 pandemic. Chaos Solitons Fractals 138, 110007 (2020)

14. Kumar, S., Cao, J., Abdel-Aty, M.: A novel mathematical approach of COVID-19 with non-singular fractional derivative. Chaos Solitons Fractals 139, 110048 (2020)

15. Okuonghae, D., Omame, A.: Analysis of a mathematical model for COVID-19 population dynamics in Lagos, Nigeria. Chaos Solitons Fractals 139, 110032 (2020)

16. Postnikov, E.B.: Estimation of COVID-19 dynamics “on a back-of-envelope”: Does the simplest SIR model provide quantitative parameters and predictions? Chaos Solitons Fractals 135, 109841 (2020)

17. Shah, K., Abdeljawad, T., Mahariq, I., Jarad, F.: Qualitative analysis of a mathematical model in the time of COVID-19 BioMed Res. Int. 2020, Article ID 5098598 (2020). https://doi.org/10.1155/2020/5098598

18. Torrealba-Rodriguez, O., Conde-Gutiérrez, R.A., Hernández-Javiera, A.L.: Modeling and prediction of COVID-19 in Mexico applying mathematical and computational models. Chaos Solitons Fractals 138, 109946 (2020)

19. Yousaf, M., Muhammad, S.Z., Muhammad, R.S., Shah, H.K.: Statistical analysis of forecasting COVID-19 for upcoming month in Pakistan. Chaos Solitons Fractals 138, 109926 (2020)

20. Zeb, A., Alzahrani, E., Erturk, V.S., Zaman, G.: Mathematical model for coronavirus disease 2019 (COVID-19) containing Isolation class. BioMed Res. Int.. https://doi.org/10.1155/2020/3452402

21. Zhang, Z.: A novel Covid-19 mathematical model with fractional derivatives: singular and nonsingular kernels. Chaos Solitons Fractals 139, 110060 (2020). https://doi.org/10.1016/j.chaos.2020.110060

22. Zhang, Z., Zeb, A., Hussain, S., Alzahrani, E.: Dynamics of COVID-19 mathematical model with stochastic perturbation. Adv. Differ. Equ. 2020, 451 (2020). https://doi.org/10.1186/s13662-020-02909-1

23. Abdon, A.: Blind in a commutative world: simple illustrations with functions and chaotic attractors. Chaos Solitons Fractals 114, 347–363 (2018). https://doi.org/10.1016/j.chaos.2018.07.022

24. Abdon, A.: Fractional discretization: the African’s tortoise walk. Chaos Solitons Fractals 130, 109399 (2020). https://doi.org/10.1016/j.chaos.2019.109399

25. Abdon, A.: Fractal-fractional differentiation and integration: connecting fractal calculus and fractional calculus to predict complex system. Chaos Solitons Fractals 102, 396–406 (2017)

26. van den Driessche, P., Watmough, J.: Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission. Math. Biosci. 180, 29–48 (2002)

27. Mickens, R.E.: A SIR-model with square-root dynamics: an NSFD scheme. Differ. Equ. Appl. 16(2–3), 209–216 (2010)

28. Mickens, R.E.: Nonstandard Finite Difference Models of Differential Equations. World Scientific, Singapore (1994)

## Acknowledgements

This research was supported by Higher Education Institutions of Anhui Province (No. KJ2020A0002).

## Funding

Supported by Higher Education Institutions of Anhui Province under (KJ2020A0002).

## Author information

Authors

### Contributions

The authors equally contributed in preparing this manuscript. All authors read and approved the final manuscript.

### Corresponding author

Correspondence to Anwar Zeb.

## Ethics declarations

### Competing interests

The authors declare that they have no competing interests.

## Rights and permissions
{"ft_lang_label":"__label__en","ft_lang_prob":0.85907835,"math_prob":0.9988852,"size":20000,"snap":"2022-27-2022-33","text_gpt3_token_len":5165,"char_repetition_ratio":0.12012403,"word_repetition_ratio":0.032608695,"special_character_ratio":0.27115,"punctuation_ratio":0.15947181,"nsfw_num_words":0,"has_unicode_error":false,"math_prob_llama3":0.99989223,"pos_list":[0,1,2],"im_url_duplicate_count":[null,2,null],"WARC_HEADER":"{\"WARC-Type\":\"response\",\"WARC-Date\":\"2022-08-14T03:32:55Z\",\"WARC-Record-ID\":\"<urn:uuid:f2988130-79b4-4050-b1cf-8a96b91cd3e7>\",\"Content-Length\":\"275041\",\"Content-Type\":\"application/http; msgtype=response\",\"WARC-Warcinfo-ID\":\"<urn:uuid:ab51628a-5034-4ff5-b842-0d5fd0465f35>\",\"WARC-Concurrent-To\":\"<urn:uuid:49d9712b-1b01-4865-9dbb-eaee3c399400>\",\"WARC-IP-Address\":\"146.75.32.95\",\"WARC-Target-URI\":\"https://advancesindifferenceequations.springeropen.com/articles/10.1186/s13662-020-03137-3\",\"WARC-Payload-Digest\":\"sha1:PF2TXNS6F4MEQQMEHM4KJIQU3KNZY426\",\"WARC-Block-Digest\":\"sha1:BQZ3YUFGAOQO5KXWHYLBDHJWNER73DR4\",\"WARC-Identified-Payload-Type\":\"text/html\",\"warc_filename\":\"/cc_download/warc_2022/CC-MAIN-2022-33/CC-MAIN-2022-33_segments_1659882571993.68_warc_CC-MAIN-20220814022847-20220814052847-00623.warc.gz\"}"}