Dynamics of Salmonella Infection

Abstract

In this chapter, we propose a mathematical epidemic model, with integer and fractional order to describe the dynamics of Salmonella infection in animal herds. We investigate the qualitative behaviors of such model and find the conditions that guarantee the asymptotic stability of disease‐free and endemic steady states. To assess the severity of the outbreak, as well as the strength of the medical and/or behavioral interventions necessary for control, we estimate basic reproduction number R0. This threshold parameter specifies the average number of secondary infections caused by one infected individual during his/her entire infectious period at the start of an outbreak. We also provide an unconditionally stable implicit scheme for the fractional‐order epidemic model. The theoretical and computational results give insight into the modelers and infectious disease specialists.

1. Introduction

Mathematical epidemic models, for Salmonella infections, provide a comprehensive framework for understanding the disease transmission behaviors and for evaluating the effectiveness of different intervention strategies [1, 2]. We recall here that the Salmonella infection, a major zoonotic disease, is transmitted between humans and other animals. Reports conducted by the National Center for Emerging and Zoonotic Infectious Diseases (NCEZID) revealed that the number of people infected by Salmonella, over the past few years, has remained increasing. The most commonly developed symptoms of Salmonella include diarrhea, fever, and abdominal cramps that appear 12–72 hours after infection. The infected people usually recover without medical aid within a period of 4–7 days [3, 4]. However, hospitalization may be needed for some infected people in the case of severe diarrhea. Salmonella is found living in the intestinal tracts of not only humans but also other creatures such as birds. The transmission of bacterium to humans occurs through the ingestion of food that has been contaminated with animal feces. These contaminated foods are commonly from an animal source, such as beef, poultry, milk, or eggs [5]. However, vegetables and other foods may also become contaminated. Additionally, foods that have been contaminated are almost impossible to detect while eating, due to their normal taste and smell. Therefore, Salmonella is considered as a serious problem for the public health throughout the world. There are no doubts that mathematical modeling of Salmonella infection plays an important role in gaining understanding of the transmission of the disease in a specific environment and to predict the behavior of any outbreak. Furthermore, mathematical analysis leads to determining the nature of equilibrium states and to suggest recommended actions to be taken by decision makers to control the spreading of the disease. The objective of this work is to adopt the fractional‐order epidemic model to describe the dynamics of Salmonella infections in animal herds.

Fractional‐order (or free‐order) differential models have been successfully applied to system biology, physics, chemistry, and biochemistry, hydrology, medicine, and finance (see, e.g., [6–12] and the references therein). In many cases, they are more contestant with the real phenomena than the integer‐order models, because the fractional derivatives and integrals enable the description of the memory and hereditary properties inherent in various materials and processes. Hence, there is a growing need to study and use the fractional‐order differential and integral equations in epidemiology and biological systems with memory [13]. However, analytical and closed solutions of these types of fractional equations cannot generally be obtained. As a consequence, approximate and numerical techniques are playing an important role in identifying the solution behavior of such fractional equations and exploring their applications (see, e.g., [14–16] and the references therein).

A large number of work done on modeling biological systems have been restricted to integer‐order ordinary (or delay) differential equations (see, e.g., [17–22]). In Ref. [23], the authors proposed the classical Susceptible‐Infected‐Recovered (SIR) model. The authors in Ref. [24] introduced a new compartment into the SIR model, which is called cross‐immune compartment to be called SIRC model. The added compartment cross‐immune $C(t)$ describes an intermediate state between the fully susceptible $S(t)$ and the fully protected $R(t)$ one. A fractional‐order SIRC model of influenza, a disease in human population, was discussed in Ref. [25]. In the present chapter, we consider the fractional‐order SIRC model associated with evolution of Salmonella infection in animal herds. However, we will take into account the disease‐induced mortality rate $m$ in the model. Qualitative behavior of the fractional‐order SRIC model is then investigated. Numerical simulations of the fractional‐order $\text{SRIC}$ model are provided to demonstrate the effectiveness of the proposed method by using implicit Euler's method.

Definitions of fractional‐order integration and fractional‐order differentiation/integration are given in Appendix.

2. Construction of the model

Assume that the Salmonella infection spreads in animal herds which are grouped as four compartments, according to their infection status: $S(t)$ is the proportion of susceptible at time $t$ (individuals that do not have the infection), $I(t)$ is the proportion of infected individuals (that have the infection), $R(t)$ is the proportion of recovered individuals (that recovered from the infection and have temporary immunity), and $C(t)$ is the proportion of cross‐immune individuals at time $t$. The total number of animals in the herd is given by $N=S+I+R+C$. We consider that initially all the animals are susceptible to the infection. Once infected, a susceptible individual leaves the susceptible compartment and enters the infectious compartment where it then becomes infectious. The infected animals pass into the recovered compartment. After recovery from an infection animals, the individuals enter a new class $C(t)$. Therefore, we consider the disease transmission model consists of nonnegative initial conditions together with system of equations.

$\begin{array}{l}\dot{S}(t)=\mu N+\eta C(t)-(\beta I(t)+\mu )S(t)\mathrm{,}\\ \dot{I}(t)=\beta S(t)I(t)+\sigma \beta C(t)I(t)-(\theta +m+\mu )I(t)\mathrm{,}\\ \dot{R}(t)=(1-\sigma )\beta C(t)I(t)+\theta I(t)-(\mu +\delta )R(t)\mathrm{,}\\ \dot{C}(t)=\delta R(t)-\beta C(t)I(t)-(\eta +\mu )C(t)\mathrm{.}\end{array}$E1

Here $\text{'}{.}^{\prime }=D=\frac{d}{dt}$. The parameter $\mu$ denotes the mortality rate in every compartment and is assumed to equal the rate of newborns in the population. $\beta$ is the contact rate and also called the transmission rate for susceptible to be infected. ${\eta }^{-1}$ is the cross‐immune period, while ${\theta }^{-1}$ is the infectious period and ${\delta }^{-1}$ is the total immune period. $\sigma$ represents the fraction of the exposed cross‐immune individuals who are recruited in a unit time into the infective subpopulation [24, 26]. The presented model (1) differs from existing model, we assume a disease induced mortality rate $m$; see the diagram of Figure 1.

Figure 1.

Schematic diagram of SIRC epidemic model for Salmonella infection.

2.1. Fractional‐order SIRC epidemic model

Most of biological systems have long‐range temporal memory. Modeling of such systems by fractional‐order (or arbitrary order) models provides the systems with long‐time memory and gains them extra degrees of freedom [27]. A large number of mathematical models, based on ordinary and delay differential equations with integer‐orders, have been proposed in modeling the dynamics of epidemiological diseases [18, 20, 28, 29]. In recent years, it has turned out that many phenomena in different fields can be described very successfully by models using fractional‐order differential equations (FODEs) [13, 6, 27]. This is due to the fact that fractional derivatives enable the description of the memory and hereditary properties inherent in various processes. Herein, we replace the integer‐order of the model (1) into a fractional‐order (or free‐order) and assume that $s(t)=S(t)/N\mathrm{,}i(t)=S(t)/N\mathrm{,}r(t)=R(t)/N\mathrm{,}c(t)=C(t)/N$, where $N$ is the total number of population. Then the model with a fractional‐order $\alpha$ ($0<\alpha \le 1$) takes the form

$\begin{array}{c}{D}^{\alpha }s(t)=\mu +\eta c(t)-(\beta i(t)+\mu )s(t)\mathrm{,}\\ D^{\alpha }i(t)=\beta s(t)i(t)+\sigma \beta c(t)i(t)-(\theta +m+\mu )i(t)\mathrm{,}\\ D^{\alpha }r(t)=(1-\sigma )\beta c(t)i(t)+\theta i(t)-(\mu +\delta )r(t)\mathrm{,}\\ D^{\alpha }c(t)=\delta r(t)-\beta c(t)i(t)-(\eta +\mu )c(t)\mathrm{.}\end{array}$E2

Here,

$D^{\alpha }f(t)=\frac{1}{\Gamma (n-\alpha )}{\left(\frac{d}{dt}\right)}^n{\displaystyle {\int }_0^t}{(t-s)}^{\alpha -n-1}f(s)ds\mathrm{.}$E3

When $0<\alpha \le 1$,

$D^{\alpha }f(t)=\frac{1}{\Gamma (1-\alpha )}{\displaystyle {\int }_0^t}\frac{f^{\prime }(s)}{{(t-s)}^{\alpha }}ds\mathrm{.}$E4

(The initial conditions $s(0)=s_0$, $i(0)=i_0$, $r(0)=r_0$ should be given.) We note that the fractional derivatives involve an integration and are nonlocal operators, which can be used for modeling systems with memory; see the Appendix.

2.2. Stability criteria for the epidemic SIRC model (2)

To find the equilibria of the model (2), we put $D^{\alpha }s(t)=D^{\alpha }i(t)=D^{\alpha }r(t)=D^{\alpha }c(t)=0$. We have disease‐free (infection‐free) equilibrium state ${\mathcal{E}}_0$ and endemic equilibrium state ${\mathcal{E}}_{+}$:

${\mathcal{E}}_0=(\mathrm{1,0,0,0})\text{and}{\mathcal{E}}_{+}=(s^{*}\mathrm{,}{i}^{*}\mathrm{,}{r}^{*}\mathrm{,}{c}^{*})\mathrm{,}$E5

where

$\begin{array}{c}{s}^{*}=\frac{\theta +m+\mu }{\beta }-\sigma (\frac{\delta \theta i^{*}}{(\mu +\delta \sigma )\beta i^{*}+(\mu +\delta )(\mu +\eta )})\mathrm{,}\\ r^{*}=\frac{\theta i^{*}(\beta i^{*}+\eta +\mu )}{(\mu +\delta \sigma )\beta i^{*}+(\mu +\delta )(\mu +\eta )}\mathrm{,}\\ c^{*}=\frac{\theta \delta i^{*}}{(\mu +\delta \sigma )\beta i^{*}+(\mu +\delta )(\mu +\eta )}\mathrm{.}\end{array}$E6

The positive endemic equilibrium ${\mathcal{E}}_{+}=(s^{*}\mathrm{,}{i}^{*}\mathrm{,}{r}^{*}\mathrm{,}{c}^{*})$ satisfies Eq. (2) and $i^{*}$ is the positive root of $A_1{i}^{*2}+A_2{i}^{*}+A_3$, where

$\begin{array}{c}{A}_1=-{\beta }^2[m(\mu +\delta \sigma )+\mu (\theta +\mu +\delta \sigma )]\mathrm{,}\\ A_2=\beta [\beta \mu (\mu +\delta \sigma )+\eta \theta \delta -(\theta +m+\mu )[(\mu +\delta )(\mu +\eta )+(\mu +\delta \sigma )]+\mu \delta \theta ]\mathrm{,}\\ A_3=\beta \mu (\mu +\delta )(\mu +\eta )\left[1-\left(\frac{\theta +m+\mu }{\beta }\right)\right]\mathrm{.}\end{array}$E7

The Jacobian matrix of the model (2) is

$J=\left(\begin{array}{cccc}-\beta i(t)-\mu & -\beta s(t)& 0& \eta \\ \beta i(t)& \beta s(t)+\sigma \beta c(t)-(\theta +m+\mu )& 0& \sigma \beta i(t)\\ 0& (1-\sigma )\beta c(t)+\theta & -(\mu +\delta )& (1-\sigma )\beta i(t)\\ 0& -\beta c(t)& \delta & -\beta i(t)-(\eta +\mu )\end{array}\right)\mathrm{.}$E8

2.3. The reproduction number ${\mathcal{R}}_0$

The basic reproduction number

The number of individuals infected by a single infected individual placed in a totally susceptible population.

${\mathcal{R}}_0$ that includes the indirect transmission may be obtained using next‐generation matrix method [30]. The spectral radius of the next generation matrix $(F{V}^{-1})$, which is the dominant eigenvalue of the same matrix, gives the value of ${\mathcal{R}}_0$. Then, the basic reproductive number ${\mathcal{R}}_0$ is obtained by the form

${\mathcal{R}}_0=\rho (F{V}^{-1})\mathrm{,}$E9

where the matrices $F={\left[\frac{\partial {ℱ}_i(x)}{\partial x_j}\right]}_{x=x_0}$ and $V={\left[\frac{\partial {\mathcal{V}}_i(x)}{\partial x_j}\right]}_{x=x_0}$. ${ℱ}_i(x)$, where $x$ is the set of all disease‐free states in the compartment $i$, is the rate of appearance of new infections in the compartment $i$, and $V_i(x)$ is the net transfer rate (other than infections) of the compartment $i$. The net transfer rate is given by ${\mathcal{V}}_i={\mathcal{V}}_i^{-}-{\mathcal{V}}_i^{+}\mathrm{,}$ where ${\mathcal{V}}_i^{-}$ is the rate of transfer of individuals out of the compartment $i$ and ${\mathcal{V}}_i^{+}$ is the rate of transfer of individuals into the compartment $i$ by all other means. Therefore, the disease transmission model consists of nonnegative initial conditions, $x_i(0)$, together with the following system of equations:

${x^{\prime }}_j=f_j(x)\equiv {ℱ}_j(x)-{\mathcal{V}}_j\mathrm{,}j\ge 1.$E10

From the model (2), we have

$\begin{array}{c}F=\left(\begin{array}{cc}\frac{\partial {ℱ}_1}{\partial i(t)}& \frac{\partial {ℱ}_1}{\partial r(t)}\\ \frac{\partial {ℱ}_2}{\partial i(t)}& \frac{\partial {ℱ}_2}{\partial r(t)}\end{array}\right)=\left(\begin{array}{cc}\beta s& 0\\ 0& 0\end{array}\right),\\ V=\left(\begin{array}{cc}\frac{\partial {\mathcal{V}}_1}{\partial i(t)}& \frac{\partial {\mathcal{V}}_1}{\partial r(t)}\\ \frac{\partial {\mathcal{V}}_2}{\partial i(t)}& \frac{\partial {\mathcal{V}}_2}{\partial r(t)}\end{array}\right)=\left(\begin{array}{cc}\theta +m+\mu & 0\\ -\theta & \mu +\delta \end{array}\right)\mathrm{.}\end{array}$E11

Since we have only two distinct stages namely $I(t)$ and $R(t)$; it follows that both $F$ and $V$ are $2\times 2$ square matrices. Furthermore, it can be noticed that $F$ is nonnegative and $V$ is nonsingular. The basic reproductive number ${\mathcal{R}}_0$ is the dominant eigenvalue of the matrix $F{V}^{-1}$, which is obtained by solving the characteristic equation $(F{V}^{-1})I-\Lambda I=0$ where $\Lambda$ is the eigenvalue and $I(t)$ is the identity matrix. At the disease‐free equilibrium, ${\mathcal{E}}_0=(\mathrm{1,0,0,0})$, we have

${\mathcal{R}}_0=\frac{\beta }{\theta +m+\mu }\mathrm{.}$E12

The following theorem states that ${\mathcal{R}}_0$ is a threshold parameter for the stability of the model (2).

Theorem 1 The disease‐free equilibrium is locally asymptotically stable and the infection will die out if ${\mathcal{R}}_0<1$ and is unstable if ${\mathcal{R}}_0>1$. Conversely, the endemic equilibrium ${\mathcal{E}}_{+}$ is stable when ${\mathcal{R}}_0>1$ and

$a_i>\mathrm{0,}i=\mathrm{1,2,3,4,}{a}_1{a}_2-a_3>0\text{and}{a}_1{a}_2{a}_3-a_1^2{a}_4-a_3^2>\mathrm{0,}$E13

where

$\begin{array}{c}{a}_1=(D_1+D_3+D_5)\mathrm{,}\\ a_2=(D_1{D}_3-D_4\delta +D_1{D}_5+D_3{D}_5+{\beta }^2{i}^{*}{s}^{*}+\sigma {\beta }^2{c}^{*}{i}^{*})\mathrm{,}\\ a_3=(D_1{D}_3{D}_5-D_1{D}_4\delta +D_3{\beta }^2{i}^{*}{s}^{*}+D_5{\beta }^2{i}^{*}{s}^{*}+{\beta }^2{c}^{*}\eta i^{*}-D_2\sigma \beta \delta i^{*}+\\ \sigma {\beta }^2{D}_1{c}^{*}{i}^{*}+\sigma D_3{\beta }^2{c}^{*}{i}^{*})\mathrm{,}\\ a_4=D_3{D}_5{\beta }^2{i}^{*}{s}^{*}-D_2\beta \delta \eta i^{*}+D_3{\beta }^2{c}^{*}\eta i^{*}-D_4{\beta }^2\delta i^{*}{s}^{*}-\sigma \beta \delta D_1{D}_2{i}^{*}+\sigma D_1{D}_3{\beta }^2{c}^{*}{i}^{*}\mathrm{,}\end{array}$E14

and

$\begin{array}{c}{D}_1=\beta i+\mu \mathrm{,}\\ D_2=(1-\sigma )\beta c^{*}+\theta \mathrm{,}\\ D_3=(\mu +\delta )\mathrm{,}\\ D_4=(1-\sigma )\beta i^{*}\mathrm{,}\\ D_5=\beta i^{*}+(\eta +\mu )\mathrm{,}\\ D_5=\beta i^{*}+\mu \mathrm{.}\end{array}$E15

Proof The disease‐free equilibrium is locally asymptotically stable if all the eigenvalues, ${\lambda }_i$ $i=\mathrm{1,2,3,4.}$ of the Jacobian matrix, $J({\mathcal{E}}_0)$ satisfy the following condition

$|\text{arg}({\lambda }_i)|>\frac{\alpha \pi }{2}\mathrm{.}$E16

where

$J({\mathcal{E}}_0)=\left(\begin{array}{cccc}-\mu & -\beta & 0& \eta \\ 0& \beta -(\theta +m+\mu )& 0& 0\\ 0& 0& -(\mu +\delta )& 0\\ 0& 0& \delta & -(\eta +\mu )\end{array}\right)\mathrm{.}$E17

The eigenvalues of the Jacobian matrix $J({\mathcal{E}}_0)$ are

${\lambda }_1=-\mu \mathrm{,}{\lambda }_2=\beta -(\theta +m+\mu )\mathrm{,}{\lambda }_3=-(\mu +\delta )\mathrm{,}{\lambda }_4=-(\eta +\mu )\mathrm{.}$E18

Hence ${\mathcal{E}}_0$ is locally asymptotically stable if ${\mathcal{R}}_0<1$ and is unstable if ${\mathcal{R}}_0>1$.

Now, we extend the analysis to endemic equilibrium ${\mathcal{E}}_{+}$. The Jacobian matrix $J({\mathcal{E}}_{+})$ evaluated at the endemic equilibrium is

$J({\mathcal{E}}_{+})=\left(\begin{array}{cccc}-\beta i^{*}-\mu & -\beta s^{*}& 0& \eta \\ \beta i^{*}& \beta s^{*}+\sigma \beta c^{*}-(\theta +m+\mu )& 0& \sigma \beta i^{*}\\ 0& (1-\sigma )\beta c^{*}+\theta & -(\mu +\delta )& (1-\sigma )\beta i^{*}\\ 0& -\beta c^{*}& \delta & -\beta i^{*}-(\eta +\mu )\end{array}\right)\mathrm{,}$E19

with characteristic equation

${\lambda }^4+a_1{\lambda }^3+a_2{\lambda }^2+a_3\lambda +a_4=0.$E20

Using Routh‐Hurwitz stability criteria [31], the endemic equilibrium ${\mathcal{E}}_{+}$ is locally asymptotically stable provided that

$a_i>\mathrm{0,}i=\mathrm{1,2,3,4,}{a}_1{a}_2-a_3>0\text{and}{a}_1{a}_2{a}_3-a_1^2{a}_4-a_3^2>0.$E21

This completes the proof.

3. Numerical method and simulations

Since most of the FODEs do not have exact analytic solutions, so approximation and numerical techniques must be used. In addition, most of resulting biological systems are stiff,

One definition of the stiffness is that the global accuracy of the numerical solution is determined by stability rather than local error and implicit methods are more appropriate for it.

therefore, efficient use of a reliable numerical method for dealing with such problems is necessary. In this section, we provide an implicit scheme to approximate the solutions of the fractional‐order epidemic model. We also verify that the approximate solution is stable and convergent.

Consider a biological system, with fractional‐order, of the form

$\begin{array}{cc}{D}^{\alpha }y(t)=f(t\mathrm{,}y(t))\mathrm{,}& t\in [0\mathrm{,}\text{T}]\mathrm{,}\\ y^{(k)}(0)=y^{(k)}(0)\mathrm{,}& k=\mathrm{0,1,2,\dots ,}m-1.\end{array}0<\alpha \le 1$E22

Here, $y(t)=[y_1(t)\mathrm{,}{y}_2(t)\mathrm{,}\dots \mathrm{,}{y}_n(t{)]}^T$ and $f(t\mathrm{,}y(t))$ satisfy the Lipschitz condition

$\parallel f(t\mathrm{,}y(t))-f(t\mathrm{,}x(t))\parallel \le K\parallel y(t)-x(t)\parallel \mathrm{,}K>\mathrm{0,}$E23

where $x(t)$ is the solution of the perturbed system.

Theorem 2 The FODE (22) has a unique solution if Lipschitz condition (23) is satisfied and

$M=\frac{K{T}^{\alpha }}{\Gamma (\alpha +1)}<1.$E24

Proof One can apply the fractional integral operator (given in the Appendix) to the differential Eq. (22) and incorporate the initial conditions. Thus, Eq. (22) can be expressed as

$y(t)={\displaystyle \sum _{k=0}^{m-1}}{y}_0^{(k)}\frac{t^k}{k!}+\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}{(t-s)}^{\alpha -1}f(s\mathrm{,}y(s))ds\mathrm{.}$E25

which is a Volterra equation of the second kind. Define the operator $ℒ$, such that

$ℒy(t)={\displaystyle \sum _{k=0}^{m-1}}{y}_0^{(k)}\frac{t^k}{k!}+\frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}{(t-s)}^{\alpha -1}f(s\mathrm{,}y(s))ds\mathrm{..}$E26

Then, we have

$\begin{array}{l}\parallel ℒy(t)-ℒx(t)\parallel \le \frac{1}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}{(}^t\parallel f(s\mathrm{,}y(s))-f(s\mathrm{,}x(s))\parallel ds\\ \le \frac{K}{\Gamma (\alpha )}{\displaystyle {\int }_0^t}{(}^t\underset{s\in [\mathrm{0,}\text{T}]}{{\displaystyle \sup }}|y(s)-x(s)|ds\\ \le \frac{K}{\Gamma (\alpha )}\parallel y-x\parallel {\displaystyle {\int }_0^t}{s}^{\alpha -1}ds\\ \le \frac{K{\text{T}}^{\alpha }}{\Gamma (\alpha +1)}\parallel y-x\parallel {\text{T}}^{\alpha }\mathrm{.}\end{array}$E27

Then, we have

$\parallel ℒy(t)-ℒx(t)\parallel \le \text{M}\parallel y-x\parallel \mathrm{.}$E28

Using the Banach contraction principle, we can prove that that $ℒ$ has a unique fixed point which means that the problem has a unique solution. $\square$

Many efficient numerical methods have been proposed to solve the FODEs [14, 32]. Among them, the so‐called predictor‐corrector algorithm is a powerful technique for solving the FODEs, and considered as a generalization of the Adams‐Bashforth‐Moulton method. The modification of the Adams‐Bashfourth‐Moulton algorithm is proposed by Diethelm [14, 33–34] to approximate the fractional‐order derivative. However, the converted Volterra integral equation (25) is with a weakly singular kernel, such that regularization is not necessary anymore. In our case, the kernel may not be continuous, and therefore the classical numerical algorithms for the integral part of Eq. (25) are unable to handle the solution of Eq. (22). Therefore, we implement the implicit Euler's scheme to approximate the fractional‐order derivative.

Given fractional‐order model (Eq. (22)) and mesh points $\mathcal{T}=\{t_0\mathrm{,}{t}_1\mathrm{,}\dots \mathrm{,}{t}_N\}$, such that $t_0=0$ and $t_N=\text{T}$. Then a discrete approximation to the fractional derivative can be obtained by a simple quadrature formula, using the Caputo fractional derivative (42) of order $\alpha$, $0<\alpha \le 1$, and using implicit Euler's approximation as follows (see [15]):

$\begin{array}{c}{D}_{*}^{\alpha }{x}_i(t_n)=\frac{1}{\Gamma (1-\alpha )}{\displaystyle {\int }_0^t}\frac{d{x}_i(s)}{ds}{(}^{t_n}ds\\ \approx \frac{1}{\Gamma (1-\alpha )}{\displaystyle \sum _{j=1}^n}{\displaystyle {\int }_{(j-1)h}^{jh}}\left[\frac{x_i^j-x_i^{j-1}}{h}+O(h)\right]{(}^{n}ds\\ =\frac{1}{(1-\alpha )\Gamma (1-\alpha )}{\displaystyle \sum _{j=1}^n}\left\{\left[\frac{x_i^j-x_i^{j-1}}{h}+O(h)\right]\left[{(n-j+1)}^{1-\alpha }-{(n-j)}^{1-\alpha }\right]\right\}{h}^{1-\alpha }\\ =\frac{1}{(1-\alpha )\Gamma (1-\alpha )}\frac{1}{h^{\alpha }}{\displaystyle \sum _{j=1}^n}\left[x_i^j-x_i^{j-1}\right]\left[{(n-j+1)}^{1-\alpha }-{(n-j)}^{1-\alpha }\right]+\\ \frac{1}{(1-\alpha )\Gamma (1-\alpha )}{\displaystyle \sum _{j=1}^n}\left[x_i^j-x_i^{j-1}\right]\left[{(n-j+1)}^{1-\alpha }-{(n-j)}^{1-\alpha }\right]O(h^{2-\alpha })\mathrm{.}\end{array}$E29

Setting

$\mathcal{G}(\alpha \mathrm{,}h)=\frac{1}{(1-\alpha )\Gamma (1-\alpha )}\frac{1}{h^{\alpha }}\mathrm{,}\text{and}{\omega }_j^{\alpha }=j^{1-\alpha }-{(j-1)}^{1-\alpha }\mathrm{,}(\text{where}{\omega }_1^{\alpha }=1)\mathrm{,}$E30

then the first‐order approximation method for the computation of Caputo's fractional derivative is then given by the expression

$D_{*}^{\alpha }{x}_i(t_n)=\mathcal{G}(\alpha \mathrm{,}h){\displaystyle \sum _{j=1}^n}{\omega }_j^{\alpha }\left(x_i^{n-j+1}-x_i^{n-j}\right)+O(h)\mathrm{.}$E31

From the above analysis and numerical approximation, one arrives at the following Remark.

Remark 1 The presence of a fractional differential order in a differential equation can lead to a notable increase in the complexity of the observed behavior, and the solution continuously depends on all the previous states.

3.1. Stability and convergence

Here, we prove that the suggested numerical scheme of implicit difference approximation (Eq. (31)) is unconditionally stable. It follows then that the numerical solution converges to the exact solution as $h\to 0$.

In order to study the stability of the numerical method, let us consider a test problem of linear scaler fractional differential equation

$D_{*}^{\alpha }u(t)={\rho }_{0}u(t)+{\rho }_1\mathrm{,}u(0)=u_0\mathrm{.}$E32

such that $0<\alpha \le 1$, and ${\rho }_0<0$, ${\rho }_1>0$ are constants.

Theorem 3 The fully implicit numerical approximation (31), to test problem (32) for all $t\ge 0$, is consistent and unconditionally stable.

Proof We assume that the approximate solution of Eq. (32) is of the form $u(t_n)\approx U^n\equiv {\zeta }_n$, then Eq. (32) can be reduced to

$\left(1-\frac{{\rho }_0}{G_{\alpha \mathrm{,}h}}\right){\zeta }_n={\zeta }_{n-1}+{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}({\zeta }_{n-j}-{\zeta }_{n-j+1})+{\rho }_1/G_{\alpha \mathrm{,}h}\mathrm{,}n\ge 2.$E33

Or

${\zeta }_n=\frac{{\zeta }_{n-1}+{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}({\zeta }_{n-j}-{\zeta }_{n-j+1})+{\rho }_1/G_{\alpha \mathrm{,}h}}{(1-\frac{{\rho }_0}{G_{\alpha \mathrm{,}h}})}\mathrm{,}n\ge 2.$E34

Since $\left(1-\frac{{\rho }_0}{G_{\alpha \mathrm{,}h}}\right)\ge 1$ for all $G_{\alpha \mathrm{,}h}$, then

${\zeta }_1\le {\zeta }_0\mathrm{,}$E35

${\zeta }_n\le {\zeta }_{n-1}+{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}({\zeta }_{n-j}-{\zeta }_{n-j+1})\mathrm{,}n\ge 2.$E36

Thus, for $n=2$, the above inequality implies

${\zeta }_2\le {\zeta }_1+{\omega }_2^{(\alpha )}({\zeta }_0-{\zeta }_1)\mathrm{.}$E37

Using the inequality (35) and the positivity of the coefficients ${\omega }_2$, one gets

${\zeta }_2\le {\zeta }_1\mathrm{.}$E38

Repeating the process, we have from Eq. (36)

${\zeta }_n\le {\zeta }_{n-1}+{\displaystyle \sum _{j=2}^n}{\omega }_j^{(\alpha )}({\zeta }_{n-j}-{\zeta }_{n-j+1})\le {\zeta }_{n-1}\mathrm{.}$E39

Since each term in the summation is negative. Thus ${\zeta }_n\le {\zeta }_{n-1}\le {\zeta }_{n-2}\le \dots \le {\zeta }_0$. With the assumption that ${\zeta }_n=\left|U^n\right|\le {\zeta }_0=\left|U^0\right|\mathrm{,}$ which entails $\parallel U^n\parallel \le \parallel U_0\parallel$ and we have stability.

The above numerical technique can then be used both for both linear and nonlinear problems, and it may be extended to multiterm FODEs.

3.2. Numerical simulations

The approximate solutions of epidemic model (2) are displayed in Figures 2–4, and sensitivity of ${\mathcal{R}}_0$ to transmission coefficients is displayed in Figure 5. The numerical simulations are performed by Euler's implicit scheme discussed in Section 3. We choose different fractional‐order values ($0.5<\alpha <1$), and parameter values given in Table 1. The displayed solutions in Figure 4 confirm that the fractional order of the derivative plays the role of time‐delay (or memory) in the system.

Parameter	Description	Value	Reference
$\mu$	Replacement and exit rate (day${}^{-1}$)	0.011	[35]
$\beta$	Transmission rate of susceptible to be infected (animal${}^{-1}$ day${}^{-1}$)	0.15	[35]
$\theta$	Recovery rate of infected animals day${}^{-1}$	0.16	Assumed
$m$	Disease‐induced mortality rate (day${}^{-1}$)	0.041	Assumed
$\eta$	Cross‐immune period	0.5	[36]
$\sigma$	The average reinfection probability of $C(t)$	0.06	Assumed
$\delta$	The average time of appearance of new dominant clusters	1	Assumed
$N$	The total number of population	345	Assumed

Table 1.

List of parameters.

Figure 2.

Numerical simulation of the fractional‐order epidemic model (2), when $\alpha =0.8$, and ${\mathcal{R}}_0>1$ (Each infected individual infects more than one other member of the population and a self‐sustaining group of infectious individuals will propagate), with parameter values of Table 1.

Figure 3.

Phase plane portrait for the fractional‐order endemic model (2), in absence of $C(t)$ and $R(t)$ components, when $\alpha =0.7$ (left) and $\alpha =0.9$ (right) with ${\mathcal{R}}_0=0.5<1$. We note that solution paths approach the disease‐free equilibrium ${\mathcal{E}}_0=(\mathrm{1,0,0})$.

Figure 4.

Phase plane portrait for the classic fractional‐order endemic model (2) when $\alpha =1$ (left) and $\alpha =0.9$ (right) with ${\mathcal{R}}_0=1.2>1$. We note that solution paths approach the endemic equilibrium ${\mathcal{E}}_{+}$ given by Eq. (5).

Figure 5.

Sensitivity of ${\mathcal{R}}_0$ with respect to the transmission coefficients $\beta$ and $\theta$.

4. Discussion and conclusion

In this chapter, we provided a fractional‐order SIRC epidemic model with Salmonella infection. The model provides a comprehensive framework for understanding the disease transmission behaviors, as well as for evaluating the effectiveness of different intervention strategies. We derived the sufficient conditions to preserve the asymptotical stability of disease‐free and endemic steady states. The threshold parameter (reproduction number) ${\mathcal{R}}_0$ has been evaluated in terms of contact rate, recovery rate, and other parameters in the model. The threshold parameter ${\mathcal{R}}_0$ is very sensitive to transmission coefficients $\beta$ and $\theta$ that reflects that these parameters play an important role to assess the strength of the medical and behavioral interventions necessary for control. We provided an unconditionally stable method, using Euler's implicit method for the fractional‐order differential system. The solution of a fractional‐order model at any time $t^{*}$ continuously depends on all the previous states at $t\le t^{*}$.

It has been found that fractional‐order dynamical models are more suitable to model biological systems with memory than their integer‐orders. The presence of a fractional differential order into a corresponding differential equation leads to a notable increase in the complexity of the observed behavior. However, fractional‐order differential models are as stable as their integer‐order counterpart.

Acknowledgments

This work was supported by UAE University (NRF Project UAEU‐NRF‐7‐20886).

Appendix

Let $L^1=L^1[a\mathrm{,}b]$ be the class of Lebesgue integrable functions on $[a\mathrm{,}b]$, $a<b<\infty$.

Definition 1 The fractional integral of order $\beta \in {ℝ}^{+}$ of the function $f(t)$, $t>0$ ($f:{ℝ}^{+}\to ℝ$) is defined by

$I_a^{\nu }f(t)={\displaystyle {\int }_a^t}\frac{{(t-s)}^{\nu -1}}{\Gamma (\nu )}f(s)ds\mathrm{,}t>0.$E40

The fractional derivative of order $\alpha \in (n-\mathrm{1,}n)$ of $f(t)$ is defined by two ways:

• Riemann‐Liouville fractional derivative: Take fractional integral of order $(n-\alpha )$ and then take $n^{th}$ derivative,

• Caputo fractional derivative: Take $n^{th}$ derivative and then take a fractional integral of order $(n-\alpha )$

$D_a^{\alpha }f(t)=D_a^n{I}_a^{n-\alpha }f(t)\mathrm{,}{D}_{*}^n=\frac{d^n}{d{t}^n}\mathrm{,}n=\mathrm{1,2,}\dots$E41

$D_a^{\alpha }f(t)=I_a^{n-\alpha }{D}_a^{n}f(t)\mathrm{,}n=\mathrm{1,2,}\dots \mathrm{.}$E42

We notice that the definition of time‐fractional derivative of a function $f(t)$ at $t=t_n$ involves an integration and calculating time‐fractional derivative that requires all the past history, $\text{i}\text{.e}\mathrm{.}$_, all the values of $f(t)$ from $t=0$ to $t=t_n$. Caputo's definition, which is a modification of the Riemann‐Liouville definition, has the advantage of dealing properly with initial value problems. The following Remark addresses some of the main properties of the fractional derivatives and integrals (see [12, 36–39]).

Remark 2 Let $\nu \mathrm{,}\gamma \in {ℝ}^{+}$ and $\alpha \in (\mathrm{0,1})$. Then

1. If $I_a^{\nu }:L^1\to L^1$ and $f(t)\in L^1$, then $I_a^{\nu }{I}_a^{\gamma }f(t)=I_a^{\nu +\gamma }f(t)$;

2. ${{\displaystyle \lim }}_{\nu \to n}{I}_a^{\nu }f(x)=I_a^{n}f(t)$ uniformly on $[a\mathrm{,}b]$, $n=\mathrm{1,2,3,}\dots$, where $I_a^{1}f(t)={\displaystyle {\int }_0^t}f(s)ds$;

3. ${{\displaystyle \lim }}_{\nu \to 0}{I}_a^{\nu }f(t)=f(t)$ weakly;

4. If $f(t)$ is absolutely continuous on $[a\mathrm{,}b]$, then ${{\displaystyle \lim }}_{\alpha \to 1}{D}_{*}^{\alpha }f(t)=\frac{df(t)}{dt}$;

5. Thus $D_{*}^{\alpha }f(t)=\frac{d}{dt}{I}_{*}^{1-\alpha }f(t)$ (Riemann‐Liouville sense) and $D_{*}^{\alpha }f(t)=I_{*}^{1-\alpha }\frac{d}{dt}f(t)$ (Caputo sense).

The generalized mean value theorem and another property are defined in the following Remark [40].

Remark 3

1. Suppose $f(t)\in C[a\mathrm{,}b]$ and $D_{*}^{\alpha }f(t)\in C(a\mathrm{,}b]$ for $0<\alpha \le 1$, then we have

2. If (i) holds, and $D_{*}^{\alpha }f(t)\ge 0$ $\forall t\in [a\mathrm{,}b]$, then $f(t)$ is nondecreasing for each $t\in [a\mathrm{,}b]$. If $D_{*}^{\alpha }f(t)\le 0$ $\forall t\in [a\mathrm{,}b]$, then $f(t)$ is nonincreasing for each $t\in [a\mathrm{,}b]$.