A fractional SIS model for a random process about infectivity and recovery

1. Introduction

Some processes of diseases depend on the present situation of the system and its history.¹ The classical SIR model can’t accommodate this. Therefore, in recent years, people have paid much attention to the study of introducing fractional derivatives into infectious disease models. Fractional time derivatives, including the history of the function, can also be used to merge the history of the system.^1,2,3,4,5,6 Integral derivatives concerning fractional derivatives are not the only ones, average examples are the Caputo derivative and the Riemann–Liouville derivative.⁷ For the most part, fractional derivatives are incorporated into compartment model derivatives by simply substituting Caputo for integer derivatives.^{3,8,9,10,11,12,13} The model and its analysis in math is very interesting, there is some motivation for incorporating memory effects,¹⁴ and many method equations for solving fractional differentiation have been developed.¹⁵ However, as Angstmann et al.¹⁶ pointed out, an actual epidemic model is no a priori reason to provide a physical model of the system. While these models may be capable of fitting the data, the underlying assumption is that these models may be in violation of the conservation of mass due to the positivity and dimensional consistency of some parameters.¹⁷

To solve such problems, based on the research of Angstmann et al.,^5,6,14,15 by thinking about the governing equations of potential stochastic processes, we acquire a common criterion for building fractional models. In the SIS model, the population is divided into two parts: the susceptible population (S) and the infected population (I). In the random process view of the SIS model, an individual starts in S and transitions to I with a certain probability after touching the infected. It is then probable for the infected individual to recover and re-transition to S. In the SIS model, their probabilities are linked to exponential wait time density. This model could be built through SIS compartment has to the continuous time random walk (CTRW).^5,14,15 Parameters in occurring fractional equations are plainly defined, whereupon we agree with the dynamics of individuals in the population. Furthermore, the fractional equation on the dimension is the same.

In addition, Wang et al.¹⁸ investigated the fractional SEIR model’s dynamics. Wang et al.¹⁸ investigated the fractional epidemiological model of COVID-19. Though they only show the existence of fractional model equilibria, and conducted some studies on the local stability of the equilibrium point.⁵ However, to our knowledge, there seems to be little work on the global stability of model equilibrium points. Therefore, in this paper, based on Monteiro and Mazorche’s¹⁹ dynamic study on frSIS infectious disease model, we discussed the local and global stabilities of the equilibrium point about frSIS model, considering the function of fractional infectivity and fractional recovery rate.

This paper mainly studies the derivation of the stochastic process of the frSIS model considering fractional contagion and fractional recovery, as well as the local and global stability of the equilibrium model. In Sec. 2, we provide the process of deriving the frSIS model. In Sec. 3, we studied $R_0$ of the frSIS model and proved the presence of the equilibrium point in model (25). In Sec. 4, we demonstrated the stability of $E_0$. In Sec. 5, we supply a detailed proof of the stability of the equilibrium point for endemic diseases $E^{\ast }$. In Sec. 6, we brought some numerical simulations to assist our detections. In Sec. 7, we provided an unsophisticated discussion and synopsis of the main detections.

2. Model Derivation

One of the most simple SIS model basic assumptions is that the transition of one through each compartment has nothing to do with the amount of time since the individual into a compartment. The hypothesis in mathematics represents time spent in every compartment is exponentially distributed. In some diseases with chronic contagion potential, like HPV, there is proof that the distribution about infection time has a power law tail.⁷ When we derive the generalized SIS model using CTRW^20,21 below, we introduce any time in the infection chamber.

According to the standard methods, we divide the population into two parts: S and I.²² This group consists of individuals born in the S and receiving targeted CTRW in the S and I compartments until their death and exclusion from consideration. Like standard compartments, individuals can only move from S to I and then to S. When a person is infected, it switches to the cubicle I, and when a person recovers from an infection, a transition occurs to S. The fractional diffusion equation is derived in Refs. 23 and 24. Fractional reaction–diffusion equation,^25,26 fractional order Fokker–Planck equation^27,28,29 and fractional chemotaxis diffusion equation^30,31 received very good research, and for each of these systems to provide the clear physical motives. We derive evolution equations for the SIS model and frSIS model from CTRW, comparable to deriving fractional Fokker–Planck equation with reactions³² and the master equation of CTRW with reactions on the network.³³

First, we consider individuals ill since time ${\mathcal{𝓉}}^{\prime }$. Supposing probability of the infected individual infecting a specific S within the time $\mathcal{𝓉}$ interval to $\mathcal{𝓉}+\mathrm{\Delta }\mathcal{𝓉}$ is $\zeta (\mathcal{𝓉},{\mathcal{𝓉}}^{\prime })\mathrm{\Delta }\mathcal{𝓉}+o(\mathrm{\Delta }\mathcal{𝓉})$, where $\zeta (\mathcal{𝓉},{\mathcal{𝓉}}^{\prime })$ is the disease transmission rate of each infected person, which depends on the infected person’s time of infection $\mathcal{𝓉}-{\mathcal{𝓉}}^{\prime }$ and current time $\mathcal{𝓉}$. For the convenience of expression, we have made $\mathcal{𝓉}-{\mathcal{𝓉}}^{\prime }=\mathrm{\Delta }\mathcal{𝓉}$ and $\zeta (\mathcal{𝓉},{\mathcal{𝓉}}^{\prime })=\zeta (\mathrm{\Delta }T)$, then the abbreviations in the following text are the same. The number of S at time t is $S(\mathcal{𝓉})$. For convenience, we will abbreviate $S(\mathcal{𝓉})$ as S. So in a time interval $\mathcal{𝓉}$ to $\mathcal{𝓉}+\mathrm{\Delta }\mathcal{𝓉}$, the number of newly infected individuals is expected to be $\zeta (\mathrm{\Delta }T)S\mathrm{\Delta }\mathcal{𝓉}+o(\mathrm{\Delta }\mathcal{𝓉})$.

Let the survival function $\mathrm{\Phi }(\mathrm{\Delta }T)$ be the probability that one is ill in time ${\mathcal{𝓉}}^{\prime }$ and sickened in time $\mathcal{𝓉}$ too. We stipulate that the prerequisite for being ill is to have contact with patients. Let $\tilde{f}(\mathcal{𝓉})$ represent the flow of people joining $I(\mathcal{𝓉})$ at time $\mathcal{𝓉}$, and for convenience, we will abbreviate $I(\mathcal{𝓉})$ as I. Then recursively for construction as

Initial condition $i(-{\mathcal{𝓉}}^{\prime },0)$ when time 0 said still infected individual number, and these individuals infected in earlier time. Therefore,

Then Eq. (1) could be redescribed as

The disease transmission rate $\zeta (\mathrm{\Delta }T)$ of an infected person is a function in time $\mathcal{𝓉}$ and time $\mathrm{\Delta }\mathcal{𝓉}$ of ill. This indicates that the disease transmission rate is related to both the external transmissibility of the time and the intrinsic transmissibility of the disease in the natural course. Therefore, the rate of disease transmission could be redescribed as

where $\varpi (\mathcal{𝓉})\ge 0$ and $\varrho (\mathcal{𝓉})\ge 0$ represent external and internal infectivities, respectively.

We stipulate that infected individuals are not allowed to leave I until they have died or recovered from the disease. These processes are assumed to be independent of each other, the survival function $\mathrm{\Phi }(\mathrm{\Delta }T)$ remaining in I because it is still infected at time $\mathcal{𝓉}$ could be shown as

where $\varphi (\mathrm{\Delta }\mathcal{𝓉})$ is probability an infected individual didn’t recover before time $\mathcal{𝓉}$, recovered at time ${\mathcal{𝓉}}^{\prime }$ and transitioned to R. Similarly $\vartheta (\mathrm{\Delta }T)$ is the probability that a patient will not die before time $\mathcal{𝓉}$. Supposing survival function $\vartheta (\mathrm{\Delta }T)$ of death process is signified in following form: and there are

An individual in I at time t must enter the compartment earlier, instead of leaving the compartment. Therefore, we said the number of individuals in I infection cubicle in terms of the flow into I and the survival function, i.e.

where function $I_0(\mathcal{𝓉})$ denotes the number of people are ill in time 0 and ill in time $\mathcal{𝓉}$ too. In combination with original case $i(-{\mathcal{𝓉}}^{\prime },0)$, the specific form of function $I_0$ can be signified as

By performing differential differentiation on Eq. (6), the main equation is obtained :

where $\psi (\mathcal{𝓉})=-\frac{d\varphi (\mathcal{𝓉})}{d\mathcal{𝓉}}$ is probability density function about $\varphi (\mathcal{𝓉})$. Equations (2)–(5) are substituted into Eq. (8) to obtain :

Let us eliminate a in Eq. (9) and obtain $\tilde{f}(\mathcal{𝓉})$ generalized equation. Now we combine Eq. (5) and rewrite Eq. (6) as

Since Eq. (10) is in kind of convolution, through Laplace transform $\mathcal{ℒ}\{\cdot \}$ from $\mathcal{𝓉}$ to s Eq. (10), we obtain Combining Eq. (5) again, the first integral in Eq. (9) is obtained through Laplace transform and we have Combining Eqs. (11) and (12) yields Define as the infectious memory kernel, where ${\mathcal{ℒ}}^{-1}\{\cdot \}$ is inverse Laplace transform from s to $\mathcal{𝓉}$. Similarly, combining Eq. (8), the third integral in Eq. (9) can be obtained by the same Laplace transform from $\mathcal{𝓉}$ to s, so we get Define as the healing rate memory kernel.

Substituting Eqs. (13) and (15) into Eq. (9) yields

Let initial condition be $i(-\mathcal{𝓉},0)=i_0\mathrm{\Delta }(\mathcal{𝓉})$, where $i_0$ is the normal number and $\mathrm{\Delta }(\mathcal{𝓉})$ is Dirac Delta function: Therefore, Eq. (17) is simplified as Similarly, the main equation for S can be obtained as

When a person suffers from a long-term persistent infection and has little chance of self-recovery, they are referred to as chronic infections. The assumption of the standard SIS model, that is, the waiting time of the exponential distribution, does not apply to this type of behavior. Therefore, as time goes on, we can reduce the probability of removing diseases to merge with chronic infection. In the power-law tail waiting time distribution, the expected waiting time diverges. In our SIS model, using such distribution will lead to falling into the infected individual room “rap”, until death. By using the power-law tail waiting time distribution, the SIS model is simplified to a group of infectious transition and recovery transition time of fractional derivative differential equations.

So that, we use equation Mittag-Leffler to represent $\psi (\mathcal{𝓉})$,²³

where $\tau$ is the scaling parameter. Due to the power law tail of this distribution,³⁴ there is $\psi (\mathcal{𝓉})\sim {\mathcal{𝓉}}^{-1-a}$. The dual parameter Mittag-Leffler function of $E_{a,b}(z)$ is defined as So survival function $\varphi (\mathcal{𝓉})$ is Let us substitute Eqs. (18) and (19) into the Laplace transform of cure rate memory kernel. Then Eq. (16) can be written as The Riemann–Liouville fractional derivative is defined by Ref. 35. To signify inverse Laplace transform in Eq. (20), we contemplate kind of the Laplace transform of the fractional derivative since $f(\mathcal{𝓉})$ is smooth around origin, there is then the convolution of the healing rate memory kernel can be expressed as

If the infectious memory kernel Eq. (14) also has a Laplace transform similar to Eq. (20), the fractional derivative can be combined with the infectious memory kernel, which can be satisfied by the following form of $\varrho (\mathcal{𝓉})$,

Using fractional derivatives as above, Eq. (22) can be expressed as The Laplace transform of the infectious memory core equation (14) can be expressed as so that

Finally, we obtain the equation system frSIS model :

Since model (24) is a nonautonomous dynamical system, we could use $\lambda (\mathcal{𝓉})$, $\varpi (\mathcal{𝓉})$, $\gamma (\mathcal{𝓉})$ as a normal number parameter to simplify them, i.e. $\lambda (\mathcal{𝓉})=\lambda$, $\varpi (\mathcal{𝓉})=\varpi$, $\gamma (\mathcal{𝓉})=\gamma$, $\vartheta (\mathcal{𝓉},0)=e^{-\gamma \mathcal{𝓉}}$, to obtain When $a=1$, $b=1$, model (25) is simplified to the classic SIS infectious disease model :

3. The Basic Regeneration Number $R_0$ and Random Ultimate Boundedness

Let $N(\mathcal{𝓉})=S+I$ be the total number of people and substitute it into model (25) to obtain

Therefore, there is $\frac{dN(\mathcal{𝓉})}{\lambda -\gamma N(\mathcal{𝓉})}=d\mathcal{𝓉}$. Integrating from 0 to t at both ends simultaneously yields the balance of all number of people is $N^{\ast }={lim}_{\mathcal{𝓉}\to \infty }N(\mathcal{𝓉})=\frac{\lambda }{\gamma }$.

Any nonnegative solution of model (24) will ultimately access set $𝔹$, where the specific form of $𝔹$ is given in Theorem 2. Therefore, if we merely center upon the solution of model (25) in $𝔹$, then boundary $\partial 𝔹$ of $𝔹$ is

Regeneration is the basic number of vulnerable groups and the expected number of infected individuals,^5,6 i.e.

where $N(\mathcal{𝓉})$ is the total number of people. From Eq. (7) and $i(-\mathcal{𝓉},0)=\mathrm{\Delta }(-\mathcal{𝓉})$, it can be concluded that Because of this integral is in standard form. Then its solution is Let $N^{\ast }=\frac{\lambda }{\gamma }$, we can get

Remark 3.1. From the expression of $R_0$, we acquire

(i)	item one: If $ln(\tau \gamma )>0$, then $\frac{\partial R_0}{\partial a}>0$. So, $R_0$ increases with the increase of a;
(ii)	item two: If $ln(\tau \gamma )=0$, then $\frac{\partial R_0}{\partial a}=0$. So, a does not affect the change of $R_0$;
(iii)	item three: If $ln(\tau \gamma )<0$, then $\frac{\partial R_0}{\partial a}<0$. So, $R_0$ decreases as a increases.

Similarly, we offer the relationship between $R_0$ and b.

Remark 3.2. From the expression of $R_0$, we acquire

(i)	item one: If $ln(\tau \gamma )>0$, then $\frac{\partial R_0}{\partial b}<0$. So $R_0$ decreases as b increases;
(ii)	item two: If $ln(\tau \gamma )=0$, then $\frac{\partial R_0}{\partial b}=0$. So b does not affect the change of $R_0$;
(iii)	item three: If $ln(\tau \gamma )<0$, then $\frac{\partial R_0}{\partial b}>0$. So $R_0$ increases with the increase of b.

For the existence of equilibrium state $(S^{\ast },I^{\ast })$, there must be the following limits^5,6 :

For the sake of certifying the existence of the equilibrium point, we adapt model (25) into Volterra integral equations. According to Eq. (21), we have

Multiplying $e^{\gamma \mathcal{𝓉}}$ on both sides of model (25) yields Then, integrating both sides simultaneously from 0 to $\mathcal{𝓉}$, we obtain We can organize and obtain Through Integration by parts, obtain Volterra integral equation of corresponding model (25) Let us take limit $\mathcal{𝓉}\to \infty$ on both sides of model (27) simultaneously, and obtain Because there is so we have

Let us calculate the disease-free equilibrium point $E_0=(\frac{\lambda }{\gamma },0)$ and the endemic equilibrium point $E^{\ast }=(\frac{\lambda }{\gamma R_0},\frac{\lambda (R_0-1)}{\gamma R_0})$ based on model (28). It’s easy to see that a and b have no effect on the disease-free balance $E_0$, but they have an impact on the endemic balance $E^{\ast }$.

4. Stability of the Disease-Free Equilibrium Point $E_0$

We provide relevant theorems on the stability of $E_0$ and then prove the conclusion.

Theorem 1. If $R_0<1$, then $E_0$ of model (25) is globally asymptotically stable in $𝔹$, where

While if $R_0>1$, then it is unstable.

We first prove the local asymptotic stability of $E_0$ and then certify the global asymptotic stability of $E_0$. Finally, we demonstrate an unstable situation.

4.1. Local asymptotic stability of $E_0$

Lemma 1. If $R_0<1$, then $E_0$ is locally asymptotically stable.

Proof. On account of

substituting $S=\frac{\lambda }{\gamma }$ into the above equation yields Following Oldham and Spanier,³⁶ we assume $I\propto e^{z\mathcal{𝓉}}$, it can be concluded that Dividing $e^{z\mathcal{𝓉}}$ at both ends yields Let $\mathrm{\Delta }{\mathcal{𝓉}}^{\prime }={\mathcal{𝓉}}^{\prime }-\mathcal{𝓉}$, we have

Assuming that Eq. (29) has a complex root $z=\alpha +i\beta ,\alpha \ge 0$, then

Due to Take ${\mathcal{𝓉}}_1>0$ as large as possible, so that for $\forall \epsilon >0(\epsilon <min\{{\gamma }^{-b},{\gamma }^{-a}\})$, we have $\epsilon$ as above, take ${\mathcal{𝓉}}_2>{\mathcal{𝓉}}_1$ as large enough to have then, we have

Combining Eqs. (29)–(31), for $\mathcal{𝓉}>{\mathcal{𝓉}}_2$ we get

this is a contradiction. Therefore, the root $z=\alpha +i\beta$ of Eq. (29) has a negative real root, which is $\alpha <0$, making $E_0$ locally asymptotically stable. The proof is completed. □

4.2. Global asymptotic stability of $E_0$

Lemma 2. If $R_0<1$, $E_0$ is globally asymptotically stable.

Proof. First, it is proven that if $R_0<1$, we obtain ${lim}_{\mathcal{𝓉}\to \infty }(S(\mathcal{𝓉}),I(\mathcal{𝓉}))=(\frac{\lambda }{\gamma },0)$ by inverting the parameters. Let $I^{\infty }:={limsup}_{\mathcal{𝓉}\to \infty }I(\mathcal{𝓉})$. Assuming there is $I^{\infty }>0$, then for $\forall \epsilon >0$, let’s take ${\mathcal{𝓉}}_3>0$ as large enough to

We take ${\mathcal{𝓉}}_4>{\mathcal{𝓉}}_3$ large enough to make $I<I^{\infty }+\epsilon$ and have $S\le \frac{\lambda }{\gamma }$. For Eq. (27), we can obtain Substituting $\mathrm{\Delta }{\mathcal{𝓉}}^{\prime }={\mathcal{𝓉}}^{\prime }-\mathcal{𝓉}$ into above equation yields where

If there is $\frac{\lambda \varpi {(\mathrm{\Delta }\mathcal{𝓉})}^{b-1}}{{\tau }^b\mathrm{\Gamma }(b)}>\frac{{(\mathrm{\Delta }\mathcal{𝓉})}^{a-1}}{{\tau }^a\mathrm{\Gamma }(a)},\mathrm{\Delta }\mathcal{𝓉}\in (0,\infty )$, then there is

Therefore, we take $\epsilon$ small enough that $\epsilon (1+R_0+I^{\infty }+\epsilon )<I^{\infty }(1-R_0)$, we have

If there is $\frac{\lambda \varpi {(\mathrm{\Delta }\mathcal{𝓉})}^{b-1}}{{\tau }^b\mathrm{\Gamma }(b)}<\frac{{(\mathrm{\Delta }\mathcal{𝓉})}^{a-1}}{{\tau }^a\mathrm{\Gamma }(a)},\mathrm{\Delta }\mathcal{𝓉}\in (0,\infty )$, then there is $Ĩ(\mathcal{𝓉})<0$. Similarly, if $\epsilon$ is small enough to satisfy $\epsilon (1+I^{\infty }+\epsilon )<I^{\infty }$, we have

Equations (32) and (33) contradict the definition of $I^{\infty }$. So, when $\mathcal{𝓉}\to \infty$, there is $I^{\infty }=0$ and $I\to 0$. Furthermore, based on $S=N(\mathcal{𝓉})-I$, ${lim}_{\mathcal{𝓉}\to \infty }S(\mathcal{𝓉})=\frac{\lambda }{\gamma }$ can be obtained. It is combined with the local asymptotic stability of $E_0$, and $E_0$ is globally asymptotic stable in $𝔹$. Proof completed. □

4.3. Instability of $E_0$

Subsequently, we will demonstrate the instability of $E_0$.

Lemma 3. If $R_0>1$, $E_0$ is unstable.

Proof. Substituting $z=\alpha >0$ into the characteristic equation (29) yields

Let Because of for $\forall \epsilon >0(\epsilon <min\{{(\gamma +\alpha )}^{-b},{(\gamma +\alpha )}^{-a}\})$, we take ${\mathcal{𝓉}}_5>0$ as large enough to obtain so we have ${lim}_{\alpha \to \infty }\phi (\alpha )>1,\mathcal{𝓉}>{\mathcal{𝓉}}_5$.

Moreover,

Therefore, for $\forall 0<\epsilon <min\{R_0-1,{(\gamma +\alpha )}^{-b},{(\gamma +\alpha )}^{-a}\}$, we take ${\mathcal{𝓉}}_6>{\mathcal{𝓉}}_5$ as large enough to make so there is $\phi (0)<0,\mathcal{𝓉}>{\mathcal{𝓉}}_6$. Based on the zero point theorem, in $(0,+\infty )$, as Eq. (34) has at least one positive root. Therefore, if $R_0>1$, then Eq. (29) has a positive root and $E_0$ is unstable anyway. End of proof. □

5. Stability of the Endemic Disease Equilibrium Point $E^{\ast }$

First, we provide a theorem on the stability of $E^{\ast }$ and then provide its proof process.

Theorem 2. If $R_0>1$, then $E^{\ast }$ is globally asymptotically stable at $\stackrel{\circ }{𝔹}$, where

The stability steps for proving $E_0$ are similar. First, we prove the local stability of $E^{\ast }$. Then prove the global stability of $E^{\ast }$.

5.1. Local asymptotic stability of $E^{\ast }$

So as to contemplate the local stability of local equilibrium $E^{\ast }$, following Hethcote and Driessche,³⁷ we consider Volterra integral equation of equation group equation (27) again, we have

Since ${lim}_{\mathcal{𝓉}\to \infty }S(\mathcal{𝓉})=S^{\ast }$, ${lim}_{\mathcal{𝓉}\to \infty }I(\mathcal{𝓉})=I^{\ast }$, we have Let $S=G+S^{\ast }$, $I=W+I^{\ast }$, then calculate them and substitute Eq. (35) into Eq. (27) to obtain where

Assuming

Therefore, Eq. (36) can be rewritten as The characteristic equation linearized by Eq. (37) is

Lemma 4 (Refs. 37 and 38). If the solution of Eq. (35) exists in $[0,\infty )$ and is bounded, $F(\mathcal{𝓉})\in C[0,\infty )$, $F(\mathcal{𝓉})\to 0$ as $\mathcal{𝓉}\to \infty$, $A(\mathcal{𝓉})\in L^1[0,\infty )$, $M(X)\in C^1(R^2)$, $M(0)=(0,0)$ and $\mathrm{\text{O}}=(0,0)$ is original point. Jacobian matrix $J=𝔻M(O)$ of M around O is nonsingular, and the eigenvalues of Eq. (38) have negative real parts, and for Eq. (37) $\mathrm{\text{O}}=(0,0)$, it is locally asymptotically stable.

Lemma 5. If $R_0>1$, $E^{\ast }$ is locally asymptotically stable.

Proof. It is easy to verify $F(\mathcal{𝓉})\in C[0,\infty )$ and ${lim}_{\mathcal{𝓉}\to \infty }F(\mathcal{𝓉})=0$, because

so $|J|=I^{\ast }>0$. Therefore, J is nonsingular.

On account of

by substituting the above equation into Eq. (38), it can be obtained that where

Next, we prove that if $R_0>1$, all roots of Eq. (39) have strictly negative real parts. To facilitate, let $\mu =\frac{\mathrm{\Lambda }}{\gamma }$ rewrite Eq. (39) to obtain

Let $\mu =\alpha +i\beta$ be the root of Eq. (40), where $\alpha \ge 0$. Due to the conjugation of complex roots, assuming $\beta \ge 0$, the real part $Re({(1+\alpha +i\beta )}^{b-a})\ge 1$. This indicates that $\mu$ cannot be the root of Eq. (39), that is, all roots of characteristic Eq. (39) have strict passive real parts, $E^{\ast }$ is locally asymptotically stable. Proof completed. □

5.2. Global asymptotic stability of $E^{\ast }$

Due to the in-depth work of Li et al.^39,40,41,42 on the geometric methods of global stability problems, we have determined that when $R_0>1$, all solutions of model (25) within $\stackrel{\circ }{𝔹}$ converge to $E^{\ast }$, which means $E^{\ast }$ is globally asymptotically stable in $\stackrel{\circ }{𝔹}$. Current applications can be determined in Refs. 19 and 43.

Let $\mathcal{X}\mapsto f(\mathcal{𝒳})$ be a $C^1$ function of $\mathcal{𝒳}$ in the open set $𝔻\subset R^n$. Considering autonomous systems in $R^n$

Supposing $\mathcal{𝒳}(\mathcal{𝓉},{\mathcal{𝒳}}_0)$ be the solution of system (41), where $\mathcal{𝒳}(\mathcal{𝓉},0)={\mathcal{𝒳}}_0$.

If for a certain Diagonal matrix $H=\mathrm{\text{diag}}({\epsilon }_1,\dots ,{\epsilon }_n)$, where ${\epsilon }_i$ is equal to 1 or $-1$, and HJH has nondiagonal lines for all Jacobian matrices $J=\frac{\partial f}{\partial \mathcal{𝒳}}$ of $\mathcal{𝒳}\in 𝔻$ and f, then the system (41) is considered competitive in $𝔻$. If there is $\mathcal{𝒳}(\mathcal{𝓉},{𝕂}^{\ast })\subset 𝕂$ for every compact set ${𝕂}^{\ast }\subset 𝔻$ and t is big enough, then the set $𝕂$ is called absorption in $𝔻$ for (41). Let us assuming that system (41) has an equilibrium point ${\mathcal{𝒳}}^{\ast }$, which is not general.

Additionally, let us assume that system (41) has a periodic solution $\mathcal{𝒳}(\mathcal{𝓉})=p(\mathcal{𝓉})$ with least period to $\varpi >0$ and orbit $\mathcal{𝒪}=\{p(\mathcal{𝓉}):0\le \mathcal{𝓉}\le \varpi \}$. Let us define linear system

where $\frac{\partial f^{[2]}}{\partial \mathcal{𝒳}}$ is the second additive composite matrix of the Jacobian matrix J.⁴²

An essential feature about a competitive system shows Poincaré–Bendixson property.

Lemma 6 (Poincaré–Bendixson property⁴⁴). Let us assume $\mathrm{\Upsilon }$ is the nonempty, closed and bounded limit set of a system (41) without equilibrium points, then $\mathrm{\Upsilon }$ is a closed orbit.

In Refs. 40 and 42, we provide the general principles of global stability.

Lemma 7. Suppose that

(i)	item one: $𝔻$ is simply connected;
(ii)	item two: there is a compact absorbing set $𝕂\subset 𝔻$;
(iii)	item three: system (41) satisfies Poincaré–Bendixson property;
(iv)	item four: ${\mathcal{𝒳}}^{\ast }$ is a unique equilibrium of system (41) in $𝔻$ provided it is stable.

Then ${\mathcal{𝒳}}^{\ast }$ of system (41) is globally asymptotically stable in $𝔻$.

To prove ${\mathcal{𝒳}}^{\ast }$ is a single equilibrium point in $𝔻$, it is essential to exclude the existence of periodic solutions. This is arrived at proving that every periodic solution of system (41) is asymptotically stable in orbit. The pursuing criteria for asymptotic stability of the periodic orbit of system (41) are given in Ref. 39.

Lemma 8. The sufficient condition for $\mathcal{O}=\{p(\mathcal{𝓉}):0\le \mathcal{𝓉}\le \varpi \}$ of system (41) to have asymptotic phase asymptotic orbital stability is that Eq. (42) is asymptotic stable.

Define a-exponential function⁴⁵

in view of the relationship between a-exponential function and Riemann–Liouville derivative,⁴⁵ we can get therefore we have

Let

Assuming and $\eta$ are positive numbers. By substituting Eq. (44) into $S\le \frac{\lambda }{\gamma }$, it can be calculated that Substituting Eq. (43) into the Jacobian matrix of equation model (25) yields Select $H=\mathrm{\text{diag}}(-1,1)$ and calculate based on the above equation According to inequality Eq. (45), it can be concluded that model (25) is competitive within $\stackrel{\circ }{𝔹}$.

Lemma 9. If $R_0>1$, $E^{\ast }$ is globally asymptotically stable.

Proof. As can be seen, $\stackrel{\circ }{𝔹}$ is simply concatenated and (i) holds.

Absorption in the $\stackrel{\circ }{𝔹}$ of the existence of the compact set is equivalent to proving model (25) is uniformly persistent. As a matter of fact, we can observe $\partial 𝔹$ is a positive invariant set of model (25). There is always a balance of only $E_0=(\frac{\lambda }{\gamma },0)$, it corresponds to disease free status. The utmost invariant set in $\partial 𝔹$ is a single instance $\{E_0\}$ and is isolated. According to Ref. 15, we can recognize consistent persistence of model (25) is equivalent to the instability of $E_0$. Combined with the instability of $E_0$ in the case of $R_0>1$,¹⁵ we can conclude that if $R_0>1$, model (25) is uniformly persistent. Therefore, $\partial 𝔹$ is the compact absorbing set in $𝔹$. (ii) is established.

Subsequently, we want to certify (iii). Let us assume $\mathrm{\Upsilon }$ is the omega limit set of model (25) in $\stackrel{\circ }{𝔹}$. If $\mathrm{\Upsilon }$ doesn’t contain $E^{\ast }$, then it does not contain equilibrium because $E^{\ast }$ is the only internal equilibrium. Lemma 6 will point to $\mathrm{\Upsilon }$ as a closed orbit. Let us assume $\mathrm{\Upsilon }$ contains $E^{\ast }$, according to Lemma 5, we understand whenever $E^{\ast }$ exists in $\stackrel{\circ }{𝔹}$, it is asymptotically stable, and any orbit approaching $E^{\ast }$ must converge to $E^{\ast }$. Therefore, $\mathrm{\Upsilon }=\{E^{\ast }\}$. That is to say, model (25) satisfies the Poincaré–Bendixson property.

Final verification (iv), due to Lemma 8, the second additive composite matrix of J is

The second composite system of model (25) with periodic solution $(S,I)$ is

To prove asymptotic stability of Eq. (46), consider Lyapunov function

The orbital $\mathcal{𝒪}$ of periodic solutions $(S,I)$ is consistent with the $\partial 𝔹$ is a positive distance apart. Therefore, there exists a constant $c>0$. Then we have for all $P\in R$ and $(S,I)\in \mathcal{𝒪}$. The right derivative of the solution of the V function along $(S,I)$ to inequality Eq. (47) and $(S,I)$ can be estimated as where Thus, $(S,I)$ is periodic of minimal period $\varpi$, Together with inequality equation (49), it means that when $\mathcal{𝓉}\to \infty$, there is $V(\mathcal{𝓉})\to 0$. According to inequality equation (48), when $\mathcal{𝓉}\to \infty$, there is $p(\mathcal{𝓉})\to 0$. Therefore Eq. (46) is asymptotically stable, if minimal period $\varpi >0$. When $\varpi \to 0$, inference is also established.

Therefore, combining Lemma 7, it could be concluded that $E^{\ast }$ is globally asymptotically stable in $\stackrel{\circ }{𝔹}$. □

6. Numerical Results

We provide numerical results to display the complex disease dynamics of the frSIS model (25). First, we will provide a numerical solution for model (25).

6.1. Numerical methods

First, we hand discrete calculation process of fractional operators ${}_0{D}_{\mathcal{𝓉}}^{1-a}(\frac{I}{\vartheta (0,\mathcal{𝓉})})$ and ${}_0{D}_{\mathcal{𝓉}}^{1-b}(\frac{I}{\vartheta (0,\mathcal{𝓉})})$. We assume a is $[0,1)$ and b is $(0,1]$. The time interval $[0,\mathcal{𝓉}]$ is discretized into $0={\mathcal{𝓉}}_0<{\mathcal{𝓉}}_1<\cdots <{\mathcal{𝓉}}_n=\mathcal{𝓉}$, where $\mathrm{\Delta }T$ is a constant. If $\mathrm{\Delta }T=1$, there is ${\mathcal{𝓉}}_i=i\mathrm{\Delta }T=i$ for all $i\in \{0,1,\dots ,n\}$. Therefore, we can perform discretization on Riemann–Liouville fractional derivative in model (25), as follows¹⁸:

where $\vartheta (0,{\mathcal{𝓉}}_m)=e^{-\gamma {\mathcal{𝓉}}_m},\vartheta (0,0)=1$.

In model (25), denote

following the Euler iterative method, frSIS model (25) could be discretized into the following form: where $i=0,1,2,\dots ,n-1$, and

6.2. Numerical simulation of disease dynamics

6.2.1. Disease extinction caused by fractional order a

We think about the frSIS model with constant parameters (25)

Example 6.1. We take a by $a=0.54$, 0.56 and 0.58, b keep unchanged as $b=1$. The results of equilibrium between the basic reproductive number $R_0$ and a model (25) with different a are shown in Table 1. When $a=0.54$, 0.56 and 0.58, $R_0<0$ is 0.9931, 0.9445 and 0.8978, respectively. This means that model (25) has a unique $E_0=(4,0)$. In Fig. 1, we show a plot of S and I with time for $a=0.54$, $0.56$ and $0.58$ under the initial conditions $S(0)=9$ and $I(0)=1$. It is simple to perceive that S as a whole drops sharply in a short period and then slowly rises to a stable state. I as a whole rises to a maximum for a short time and then slowly declines back to a stable state. It is easy to see that with $b=1$ held constant, the smaller a, the faster S decreases in the descending phase and the faster I rises in the ascending phase. We can conclude that in this case, keeping b unchanged and controlling the size of a can more effectively control the dynamics of disease extinction.

6.2.2. Disease extinction caused by fractional order b

Example 6.2. We take b as $b=0.2$, 0.4 and 0.6, and a remains unchanged as $a=0$. The results of equilibrium between the basic reproductive number $R_0$ and a model (25) with different b are shown in Table 2. When $b=0.2$, 0.4 and 0.6, $R_0<0$ is 0.2731, 0.4972 and 0.9051, respectively. This means that model (25) has a unique $E_0=(4,0)$. In Fig. 2, we show S and I in relation to time for $b=0.2$, 0.4 and 0.6 under the initial conditions $S(0)=9$ and $I(0)=1$. It is simple to perceive that S as a whole drops sharply in a short period and then slowly rises to a stable state. I as a whole rises to a maximum for a short time and then slowly declines back to a stable state. It is easy to see that if $a=0$ remains unchanged, the greater b, the faster S will fall in the descending stage and the faster I will rise in the ascending stage. We can conclude that in this case, keeping a unchanged and controlling the size of b can more effectively control the dynamics of disease extinction.

6.2.3. Persistent disease caused by the combined action of fractional order a and b

We consider the frSIS model with constant parameters model (25)

Example 6.3. We take a as $a=0.1$, 0.2, 0.3, 0.4 and 0.5, and the corresponding b as $b=0.6$, 0.7, 0.8, 0.9 and 1. The results of equilibrium between the basic reproductive number $R_0$ and a model (25) with different a and b are shown in Table 3. When $a=0.1$, 0.2, 0.3, 0.4 and 0.5 and the corresponding $b=0.6$, 0.7, 0.8, 0.9 and 1, $R_0>0$ is 1.6487, 1.8307, 2.0012, 2.1561 and 2.2932, respectively. This means that model (25) has a unique $E^{\ast }=(2.4261,1.5739)$, $(2.1850,1.8150)$, $(1.9988,2.0012)$, $(1.8552,2.1448)$ and $(1.7443,2.2557)$, respectively. In Fig. 3, we show a plot of S and I in relation to time under the initial conditions $S(0)=9$ and $I(0)=1$. It is simple to perceive that S as a whole drops sharply in a short period and then slowly rises to a stable state. I as a whole rises to a maximum for a short time and then slowly declines back to a stable state. It is easy to see that with the slow increase of a and b, $R_0$ also slowly increases, the S in $E^{\ast }$gradually decreases, and the I gradually increases. We can conclude that, in this case, controlling the size of a and b simultaneously can more effectively control the dynamics of the disease persistence.

7. Conclusion and Discussion

We derived an frSIS model that opposes the classical SIS model by contemplating the power law tail distribution of disease infection rate and cure rate. During the derivation process, we demonstrated how fractional derivatives can be incorporated into the model without violating the physical properties of the model parameters. The biological motivation of this modeling method is to observe that the longer a person is infectious during the transmission of certain diseases, the more likely they are to remain infectious. We introduced the basic reproduction number $R_0$, which shows a and b The extinction and sustained impact on the disease. In short, we have summarized our main findings and their related biological significance.

(1)	The basic reproductive number $R_0$: the formula of $R_0$ and the threshold condition of whether the disease is extinct are obtained from mathematical analysis of frSIS model (25). It is worth noting that when proving $E_0$ (see Lemma 4) and $E^{\ast }$ (see Lemma 5), we first convert the frSIS model (25) into a nonlinear Volterra integral equation. Then perform an equivalent asymptotic behavior through the asymptotic behavior of linearization of the corresponding linear Volterra integral equation. In proving global stability of $E_0$ (see Lemma 5), we foremost confirm $E_0$ is a global attractor, then obtain global stability through local stability (see Lemma 4). In proving the global stability of local equilibrium $E^{\ast }$ (see Lemma 9), following Li et al. and using the Poincaré–Bendixson property (see Lemma 6) and the principle of global stability (see Lemma 7), we provide a geometric method for proving global stability. To be brief, we give a framework to demonstrate the global stability of the frSIS model (25).
(2)	The effects of a and b: For one thing, a and b will affect existence of the equilibrium point in model (25). For another, a and b will affect the value of $R_0$ (see Remarks 3.1 and 3.2). The numerical results of the relationship between S and I and time are shown in Figs. 1–3. Therefore, it can be concluded that the frSIS model (25) is caused by a and b.
(3)	Control strategy: Based on the results of numerical analysis, on the one hand, from Fig. 1, it can be seen that the frSIS model (25) tends toward $E_0$, which is the extinction of the disease. Therefore, while keeping b constant, appropriately increasing the value of a can quickly control the spread of the disease and make it disappear. Similarly, from Fig. 2, it can be seen that the frSIS model (25) tends toward $E_0$. Therefore, while keeping a constant, appropriately reducing the value of b can quickly control the spread of the disease and make it disappear. On the other hand, from Fig. 3, it can be seen that the frSIS model (25) tends toward $E^{\ast }$, meaning that the disease persists. Therefore, in the case of changes in a and b, appropriately reducing a and b can more effectively control the spread of infectious diseases.

For our frSIS model, we can extend the fractional order of infection rate and cure rate to more aspects such as mortality rate and permanent immunity. This is worth considering and researching in the future investigative.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Nos. 11361104, 12261104), the Youth Talent Program of Xingdian Talent Support Plan (No. XDYC-QNRC-2022-0514), the Yunnan Provincial Basic Research Program Project (No. 202301AT070016, No. 202401AT070036) and the Graduate Research and Innovation Fund of Yunnan University for Nationalities (2023SKY080). Hui Chen and Jia Li contributed equally to this paper.

ORCID

Hui Chen  https://orcid.org/0009-0003-7679-1734

Jia Li  https://orcid.org/0009-0003-6734-2398

Xuewen Tan  https://orcid.org/0000-0002-4710-2325