Dynamical analysis of an anthrax disease model in animals with nonlinear transmission rate

1.   Introduction

Anthrax is a highly lethal infectious disease caused by Bacillus anthracis, a Gram-positive, rod-shaped bacterium capable of forming spores ^[1]. It can lead to the sudden and unforeseen death of infected animals, including livestock, and occasionally poses a fatal threat to humans. Generally, there are two ways in which anthrax can be transmitted to animals: by ingesting a sufficient number of spores in soil (or on plants) and via stable flies, Stomoxys calcitrans, that feed on infected blood bleeding out of the carcasses of animals that died of anthrax infection ^[2,3]. Utilizing pre-exposure vaccinations 2–4 weeks before the expected endemic season, anthrax can be controlled, and, is the only effective measure ^[2]. Anthrax spores' exceptional resistance to pH extremes, heat and cold, desiccation, and a variety of chemical agents contributes to their long-term viability of up to two centuries ^[4]. Moreover, the anthrax infection often becomes more lethal with no symptoms in cattle ^[2]. In such a case, the period (from incubation to death) may be crucial for preventing anthrax epidemics.

Anthrax poses a threat not only to herbivorous animals but also to carnivores that consume the contaminated remains of anthrax-infected prey. However, the infection and mortality rates among carnivores are generally lower compared to herbivores. The World Anthrax Data Site, as reported by the World Health Organization (WHO), has documented anthrax outbreaks in animals across nearly 200 countries, particularly in regions with high animal and human populations. For instance, in Etosha National Park, a range of herbivores, from elephants to ostriches, succumbed to anthrax, and the precise cause remained unknown (see ^[5,6]). In November 2004, the deaths of three captive cheetahs (Acinonyx jubatus) were recorded in Jwana Game Reserve in Jwanen and Botswana ^[5,6]; they had been fed meat from anthrax-infected dead red hartebeest (Alcelaphus buselaphus) (revealed later). The sporadic outbreaks and epizootics of anthrax infection among livestock and wild animals are comparatively high in Africa, the Middle East, and Asia ^[7]. In Bangladesh and India, two countries in southern Asia, anthrax is extremely common, and both human and animal infections and outbreaks are frequently documented ^[8,9]. The most common cutaneous form of anthrax infection can have a mortality rate of up to 20% if no treatment measures are applied; it may reduce to $< \; 1\%$ with proper treatment measures ^[10]. According to the emergency prevention system (EMPRES-i) data from the food and agriculture organisation of the UN, seventeen anthrax outbreaks were reported during 2009–2019 only in India and Bangladesh (see Figure 1). During the period 2009–2023, several anthrax outbreaks were reported globally (see Figure 2).

The mathematical study of the prey-predator ecosystem was first introduced by Lotka and Volterra ^[11,12] (working separately), well known as the Lotka-Volterra model. After the popularity of the Lotka-Volterra model, Kermack and Kendrick ^[13] first introduced the term epidemiology in mathematical modelling and explored a SEIRS compartment model.

Following ^[13], researchers have drawn much attention to this domain in recent decades ^{[14,15,16,17,18]}. In ^[14], Asamoah et al. presented transmission dynamics of Q fever in cattle herds, a bacterial infection caused by Coxiella burnetii. They have investigated the model for optimal disease control by introducing time-dependent vaccination, environmental hygiene and culling. Extending the model proposed in ^[14] by incorporating direct transmission by both asymptomatic and symptomatic cattle, different shedding rates, treatment class, the rate of relapse and seasonal fluctuations, Asamoah et al. ^[15] determined the reproduction number of Q fever under different scenarios. They have also investigated the model for optimal control strategies under time-dependent treatment, disinfectant and separate facilities for animal birthing. Over the last few years, Dengue, a viral disease carried by mosquitoes and triggered by the dengue virus, has become a notable global health issue, especially in tropical and subtropical regions. A mathematical model on the transmission dynamics of dengue fever was presented in ^[16]. The parametric expressions associated with the spread of the disease were determined, and the transmission dynamics were explored under different time-dependent control variables.

In another study ^[17], it has been reported that the combination of human prevention strategies and vector control measures in control interventions can potentially lower the prevalence of Dengue within a community. In ^[18], Asamoah et al. presented another dengue transmission model with four time-dependent control measures: treated bednets, treatment (prophylactics), insecticides, and vaccination. Their findings highlighted that vaccination and treatment had the lowest additional cost compared to the benefits they provided. Mathematical modelling is of two types: deterministic and stochastic. In the first type, the output is fully estimated by the parameter values and the initial conditions. In contrast, the second type involves randomness, and variable states are not described by unique values but rather by probability distributions. Deterministic models are generally studied using the theories of dynamical systems. Non-linearity plays a significant role in complex natural systems, which can be studied using the theories of dynamical systems. Due to the non-linearities, researchers pay much interest in exploring the inner dynamics of the models from the dynamical system approach.

Due to the long-time survivability of anthrax spores in the soil, no country can claim to be completely free of the disease. However, anthrax outbreaks are more common in underdeveloped countries. From a mathematical modelling perspective, the pioneering work of anthrax modelling was carried out by Furniss and Hahn, who investigated an anthrax epizootic in the Kruger National Park, South Africa ^[19]. In 1983, the same authors proposed another deterministic model with some threshold results ^[20]. Friedman and Yakubu modified the model proposed by Hahn and Furniss ^[20], including factors viz. migration, growth and the natural mortality of the host ^[21]. Following the model by Hahn and Furniss ^[20], Mushayabasa proposed a new model by incorporating a fixed time delay. Their model also includes the role of human effort (as a control measure) in decontaminating infected carcases, including the soil or area from the environment ^[22]. Mushayabasa et al. presented two different models exploring transmission and control of anthrax infection in two frameworks. In the first framework, they explored the impact of the incubation period and carcass disposal. In the second framework, they investigated the effect of vector population on the spread and control of anthrax ^[23]. Kimathi and Wainaina ^[24] investigated the transmission dynamics of anthrax in animal populations with vaccination as a control measure. Osman et al.^[25] explored the co-dynamics of Anthrax and Listeriosis infections and studied the qualitative and quantitative relationship between them under prevention and treatment measures. Their study revealed that Anthrax infection could contribute to an increased risk of Listeriosis, but Listeriosis infection is not associated with the risk of anthrax. According to their findings, anthrax infection is related to a higher chance of Listeriosis infection, while Listeriosis infection is not related to an increased risk of anthrax. Motivated by the work of Kimathi and Wainaina ^[24], Rezapour et al. proposed a new fractional ordered model in the sense of the Caputo-Fabrizio fractional derivative ^[26]. Determining the basic reproduction number $\mathcal{R}_0$ of the fractional system, they analyzed the asymptotic stability of the disease-free and endemic equilibrium states.

In recent years, there has been a heightened interest among researchers in disease modelling, including the modelling of conditions such as cancer and HIV. The COVID-19 pandemic, caused by the SARS-CoV-2 virus, has recently attracted significant attention, with over 21 million confirmed cases and 758,000 deaths reported. Mathematical models have played a crucial role in forecasting various waves during the COVID-19 pandemic, and this continues to be an active area of research ^[27]. A comparative analysis of the mathematical models of COVID-19 pandemic can be found in ^[28]. Our literature review reveals that, in contrast to diseases like HIV or cancer, researchers have given relatively minimal attention to modelling anthrax disease. In underdeveloped regions, anthrax remains a persistent threat, which serves as a compelling motivation for us to delve into the mathematical modelling of anthrax. Our goal is to construct a comprehensive model that incorporates all potential biotic and abiotic factors. One prominent strategy for controlling anthrax involves implementing a vaccination policy. However, there can be delays in this process due to factors such as asymptomatic carriers among animals or limitations in the effectiveness of manpower in vaccination efforts. Such scenarios further motivate our mathematical exploration of the disease's underlying dynamics.

In this research, we propose an SIR model with an additional compartment V for the vaccinated population. The anthrax models discussed above involve bilinear disease transmission rates. However, data and evidence observed for many diseases show that disease transmission dynamics are not always as simple as presented in such models. In reality, transmission often saturates at high levels of either susceptible or infectious individuals. In the initial stages of an outbreak, when a substantial portion of the population is vulnerable, the disease can increase swiftly due to the abundance of potential hosts for the pathogen. Nevertheless, as more individuals contract the disease and develop immunity, the reservoir of susceptible hosts diminishes, making it increasingly challenging for the pathogen to identify new hosts. Consequently, this leads to a decrease in the transmission rate. Nonlinear transmission functions, such as the Holling type-II and Michaelis-Menten type functions, offer a way to illustrate the saturation effect. These functions allow the transmission rate to increase at low levels of susceptibles or infectious individuals but level off or decrease as those levels become high.

In recent years, researchers proposed many nonlinear incidence rates ^{[29,30,31,32,33]}, and these models have attracted a considerable readership. Based on our literature review, no anthrax transmission model has been proposed with a saturation effect. Our proposed model is novel, and it is an extension of the model proposed by Kimathi and Wainaina ^[24], in terms of a nonlinear transmission rate with saturation. In Section 2, we formulate the model with a nonlinear ratio-dependent disease transmission function. In Section 3, we investigate the non-negativity and boundedness of the solutions of the system. In Section 4, we list the possible equilibria, and the basic reproduction number of the epidemic is determined therein. Sections 5 and 6 deal with the stability analysis of the equilibria and analysis of transcritical bifurcations. In Section 7, we extend our model into a fractional-order model. In Section 8, we conduct numerical simulations using different plausible parameter sets to validate the analytical results. At the end of the paper, in Section 9, we conclude our study with a brief discussion of the biological implications of the mathematical results. Finally, the essential findings of this work and future scopes of the model are given in Section 10.

2.   Construction of the model

In this section, we propose an extended version of an anthrax transmission model initially presented by Kimathi and Wainaina ^[24]. In their model ^[24], the total population is broken down into four compartments, viz., susceptible, infective, recovered and vaccinated. The total animal population is given by

where $S (t)$ denotes the sub-population vulnerable to the infection, $I(t)$ denotes the sub-population with anthrax symptoms, $R (t)$ denotes the subpopulation that has developed temporary immunity after recovering from anthrax infection, and $V (t)$ represents the subpopulation vaccinated against anthrax infection. With these four compartments, the SIRV model of anthrax transmission in animals, explored by Kimathi and Wainaina ^[24], is represented by the following system:

The following presumptions constitute the basis of Model (2.1):

$(1)$ It is assumed that the animals are recruited at a constant rate $\lambda > 0$ and are susceptible to anthrax infection. The natural mortality rate of animals in each compartment is assumed to be $\mu > 0$.

$(2)$ The disease transmission is considered to be bilinear with a transmission rate $\beta > 0$. Due to high mortality rates associated with the disease, the disease-induced mortality rate is assumed to be $\theta > 0$.

$(3)$ Susceptible animals are vaccinated at the rate $\gamma > 0$, and animal recovery rate is assumed to be $\alpha > 0$.

$(4) \; \sigma > 0$ and $\tau > 0$ denote the waning recovery rate and waning immunity of vaccinated animals, respectively.

Note: Even if there is very little human involvement in the animal-to-animal transmission of anthrax sickness, discussing the disease transmission among animals primarily becomes crucial.

We consider a nonlinear ratio-dependent disease transmission function $\frac{\beta SI}{m_1 S+m_2 I}$, where $\beta$ is the disease transmission rate, and $m_1$, $m_2$ are the half saturation constants. With this assumption, the following system of ordinary differential equations (ODEs) describes the transmission dynamics of Anthrax (transmission flowchart is in Figure 3).

with initial conditions $S(0)\geq0$, $I(0)\geq0$, $R(0)\geq0$, and $V(0)\geq0$.

3.   Mathematical analysis

This section consists of the fundamental results, viz., non-negativity and boundedness of the solutions, to establish the physiological validity of the model.

3.1.   Non-negativity

Theorem 3.1. The system (2.2) is positively invariant.

Proof. Let $(S(t), I(t), R(t), V(t))$ be the solutions of the system (2.2) with positive initial conditions $(S(0), I(0), R(0), V(0))$. Integrating both sides of the susceptible prey isocline in (2.2) from $0$ to $t$, gives,

Similarly, one can easily verify that $I(t) > 0$, $R(t) > 0$ and $V(t) > 0$.

Thus, any solution starting with positive initial conditions $(S(0), I(0), R(0), V(0))$ in the interior of $\mathbb{R}_4^+$ remains there for all future time. □

Remark 3.1. From an ecological standpoint, a system with non-negative solutions signifies the persistence or survival of species within specific ecosystems. In such systems, species populations remain at levels greater than zero, indicating that these species continue to exist and contribute to the ecological balance of their respective environments.

3.2.   Boundedness

Theorem 3.2. The orbits of system (2.2) are uniformly bounded, i.e., there exists a bounded set $\mathcal{B}$ such that for every orbit $(S(t), I(t), R(t), V(t))$ of (2.2) there is a time $t_0$ such that $(S(t), I(t), R(t), V(t)) \in \mathcal{B}$ for all $t \geq t_0$.

Proof. Let us define a function $X(t) = S(t)+I(t)+R(t)+V(t)$. Then,

Now choose any $m$ with $0 < m < \min \left\{\mu, \left(\mu +\theta \right) \right\}$. Then,

The differential inequality mentioned above can be expressed in the following way

Now, by applying Lemma 2 on page 27 in Birkhoff and Rota^[34], we obtain

For any $\epsilon > 0$ define

Then, for every orbit of (2.2), there is a time $t_0$ such that $(S(t), I(t), R(t), V(t)) \in \mathcal{B}$ for all $t \geq t_0$. □

Remark 3.2. In ecology, the concept of bounded solutions is crucial for understanding the dynamics of ecosystems. Bounded solutions signify that within an ecosystem, the populations of different species are restrained and do not spiral into infinite growth. This ecological phenomenon results from various interconnected factors that collaborate harmoniously to maintain population sizes within reasonable limits, thereby mitigating the risk of overpopulation scenarios that could ultimately result in the collapse of the ecosystem.

4.   Equilibrium points and basic reproduction number

This section presents the existence conditions for the possible equilibria and their stability analysis. The basic reproduction number of the epidemic is determined, and with its help, the stability of the disease-free and endemic equilibrium is discussed. The number of secondary infections that are anticipated to be caused by an index case in a community that is entirely susceptible is known as the basic reproduction number ($\mathcal{R}_0$) and is a measure of the potential for infection spread within a population ^[35,36,37]. If $\mathcal{R}_0 < 1$, a few infected individuals entering an entirely susceptible population will fail to replace themselves, and the disease will not spread; otherwise, $\mathcal{R}_0 > 1$. $\mathcal{R}_0$ may no longer be a reliable indicator of disease transmission after the disease has started to spread within a susceptible population, as factors favoring transmission will change. In this study, we utilize the next-generation matrix method in the determination of the basic reproduction number. The ODE system (2.2) has the two following biologically feasible equilibrium points.

$(a)$ Disease-free equilibrium and basic reproduction number $\mathcal{R}_0$

At the disease-free equilibrium, there are no infections, and no recoveries, that is, $I = R = 0$, which is given by $E_0(S_0, I_0 = 0, R_0 = 0, V_0)$, where

To determine the basic reproduction number of the epidemic, we utilize the approach of the next-generation matrix method ^[38]. The next generation matrix is given by $\textbf{K} = -T\Sigma^{-1}$, where the matrix $T$ contains the transmission terms, and $\Sigma$ contains the transition terms. For our model (2.2),

Thus, the next generation matrix is calculated as,

We obtain the basic reproduction number as

where $\rho$ is the spectral radius. When there are no infections, $I = 0$ and

This implies

$(b)$ The endemic equilibrium

The endemic equilibrium state is where the anthrax disease persists and cannot be eradicated from the animal population: $E^*(S^*, I^*, R^*, V^*)$ where

The endemic equilibrium $E^*$ exists under the sufficient condition

Equivalently, condition (c) implies $\mathcal{R}_0 > 1$ and $\Lambda < 0$, which is the existence condition for the endemic equilibrium state.

5.   Stability analysis of the equilibrium points

In this section, we investigate the stability of the system at the equilibrium points. For local stability analysis, we utilize the linearization technique around each equilibrium point, giving a small perturbation. We investigate whether the model systems return to their original equilibrium point or converge to another equilibrium point or attractor.

5.1.   Stability at $E_0$

The characteristic equation of the Jacobian matrix at the equilibrium point $E_0$ has the roots

The disease-free equilibrium is locally asymptotically stable if $\Psi_4 < 0$, which implies

equivalently $\mathcal{R}_0 < 1$. Hence, the disease-free equilibrium state is locally asymptotically stable if and only if $\mathcal{R}_0 < 1$. Hence, we get the following theorem.

Theorem 5.1. The disease-free equilibrium state is locally asymptotically stable if $\mathcal{R}_0 < 1$ and is unstable if $\mathcal{R}_0 > 1$.

5.2.   Stability at $E^*$

The stability of the coexistence equilibrium point is one of the most important parts of mathematical modelling. A stable coexistence equilibrium indicates that infection is there in the ecosystem. Such an equilibrium point is necessary for the study of the dynamics of ecosystems. By varying different factors associated with the coexistence equilibrium point, one can estimate the threshold parameters for an infection-free equilibrium state. In this section, we analyze the local stability of the coexistence equilibrium by utilizing the Routh-Hurwitz theorem. The characteristic equation of the Jacobian matrix at $E^*$ has the form

According to the Routh-Hurwitz theorem, $E^*$ is locally assymptotically stable if $A_1, A_2, A_3, A_4 > 0$,

and

6.   Transcritical bifurcation

Analysis of local bifurcations, e.g., transcritical and Hopf bifurcations, is a necessary part of the dynamical analysis of nonlinear systems. Using the theories of local and global bifurcations, one can estimate the threshold parameters associated with a system above or below which the system is stable or unstable. Transcritical bifurcation analysis is a valuable tool for understanding the dynamics of disease eradication. It aids in establishing the crucial eradication thresholds, designing efficient strategies for control, optimizing resource allocation, and keeping a vigilant eye on the advancement toward the goals of disease elimination. In this section, we investigate the threshold parameters analytically for which the system undergoes a transcritical bifurcation. Transcritical bifurcation is a local bifurcation where one equilibrium point interchanges its stability with another equilibrium point. This type of bifurcation has a significant role in ecological systems, e.g., an interior equilibrium point interchanges its stability with a disease-free when a particular feasible parameter is varied. At the equilibrium $E_1$, the Jacobian of the system (2.2) has eigenvalues

and

This equilibrium will be non-hyperbolic if one of the eigenvalues becomes zero. $\lambda_4 = 0$ implies

or

or

which are the points of transcritical bifurcation. At the point $\beta^{BP}$, the endemic equilibrium $E^*$ collides with disease-free equilibrium $E_0$, resulting in the emergence of a stable disease-free equilibrium as well as an unstable disease-free equilibrium.

Theorem 6.1. At the critical parameter

the system (2.2) around the disease-free equilibrium $E_0$ undergoes a transcritical bifurcation.

Proof. The Jacobian of the system (2.2) at the disease-free equilibrium state $E_0$ has an eigenvalue equal to zero at the critical parameter

Due to the zero eigenvalue, the eigen analysis technique fails to predict the nature of the equilibrium state at the critical value $\beta^{BP}$. Therefore we use Sotomayor's theorem ^[39] to investigate the nature of the equilibrium $E_0$ at the $\beta^{TC}$. Rewrite the system (2.2) as

where

and

Let

and

be, respectively, the eigenvectors of $J_{E_0}$ and $[J_{E_0}]^T$ corresponding to the zero eigenvalue at $\beta = \beta^{TC}$, where $J$ represents the Jacobian matrix of the system (2.2).

Here,

Then, we have,

Thus, the system satisfies all the conditions of Sotomayor's theorem for transcritical bifurcation. Therefore, the disease-free equilibrium $E_0$ of the system (2.2) undergoes a transcritical bifurcation at the parameter $\beta = \beta^{BP}$. □

Theorem 6.2. At the critical parameter

the system (2.2) around the disease-free equilibrium $E_0$ undergoes a transcritical bifurcation.

Proof. Similar to Theorem 6.1. □

Theorem 6.3. At the critical parameter

the system (2.2) around the disease-free equilibrium $E_0$ undergoes a transcritical bifurcation.

Proof. Similar to Theorem 6.1. □

7.   Fractional-order model

In recent decades, mathematical modelling incorporating fractional-order differential equations (FDEs) received much attention compared to integer-order differential equations. The main reason behind the popularity of fractional-order differential equations is their memory effect or non-local effects.

The subsequent state of any function depends on both its current state and all earlier states in FDEs ^[40,41]. Biological systems that exhibit the realistic biphasic decline behavior of infection or diseases are related to the memory effect, which can be explored using FDEs. It can be used to simulate a variety of universal phenomena more accurately ^[42,43]. According to research by Heymans et al., initial conditions described in terms of Riemann-Liouville fractional derivatives can have physical significance ^[42]. Moreover, conventional differential equations cannot estimate data between two different integer values. Several forms of fractional-order operators were employed in the recent literature to circumvent these limitations ^[44,45].

Definition 7.1. ^[43] Let $g$ be any function such that $g\in \mathcal{C}^n([t_0, +\infty), \mathbb{R})$. Then, the Caputo fractional derivative of $g$ having order $\epsilon$ is defined by

where $\Gamma(\cdot)$ is the Gamma function, $n$ is a non-negative integer such that $n-1 < \epsilon < n$, and $t\geq t_0$. In particular, when $0 < \epsilon < 1$,

In model (2.2), as the internal memory effects of the biological system of the anthrax infection are not included, we extend the proposed model (2.2) into a fractional order model utilizing Caputo fractional derivatives. To modify the existing model (2.2), we convert the first-order ordinary derivative into the Caputo derivative of fractional order $\epsilon\in (0, 1]$ as follows:

The fractional order system (7.1) has similar equilibrium points as the ODE system (2.2). A sufficient condition for the local asymptotic stability of the equilibrium point $E_0$ is

(see ^[46,47]). Similarly, the sufficient condition for the local asymptotic stability of the equilibrium point $E^*$ is

(see ^[46,47]), where $\varrho$ are the roots of the polynomial (5.1).

8.   Computational results

8.1.   Computational results for the ODE system (2.2)

In this section, we conduct numerical simulations to compare the analytical findings using Matlab and MatCont in ^[48]. Matlab's inbuilt "ODE45" function is utilized to simulate the solutions of the system. For an initial population $(204,189,130,340)$ and parameters provided in Table 1, the proposed system (2.2) has an endemic equilibrium

(see Figure 4). Considering $\mathcal{R}_0$ as a function of $\beta$, $\alpha$ and $\theta$, we plotted a 3D parametric region for which $\mathcal{R}_0 < 1$ (see Figure 5). For the parameters within this region and other parameters as mentioned in Table 1, the proposed system (2.2) always has a disease-free equilibrium.

8.1.1.   Effect of contact rate $\beta$

Initiating from the endemic equilibrium point $E^*$ with varying parameter $\beta$, we observe that $E^*$ undergoes a transcritical bifurcation at $\beta = \beta^{BP}\approx 0.035$. To visualize the transcritical bifurcation, we utilize the Matlab based continuation software MatCont 7.3. Figure 6 illustrates the transcritical bifurcation where the endemic equilibrium $E^*$ exchanges its stability with a disease-free equilibrium $E_0$.

8.1.2.   Effect of disease-induced death rate $\theta$

Initiating from the endemic equilibrium point $E^*$ with varying parameter $\theta$, we observe that $E^*$ undergoes a transcritical bifurcation at $\theta = \theta^{BP}\approx 0.08$. To visualize the transcritical bifurcation, we utilize the Matlab based continuation software MatCont 7.3. Figure 7 illustrates the transcritical bifurcation where the endemic equilibrium $E^*$ exchanges its stability with a disease-free equilibrium $E_0$.

8.1.3.   Effect of animal recovery rate $\alpha$

Initiating from the endemic equilibrium point $E^*$ with varying parameter $\alpha$, we observe that $E^*$ undergoes a transcritical bifurcation at $\alpha = \alpha^{BP}\approx 0.075$. To visualize the transcritical bifurcation, we utilize the Matlab based continuation software MatCont 7.3. Figure 8 illustrates the transcritical bifurcation where the endemic equilibrium $E^*$ exchanges its stability with a disease-free equilibrium $E_0$.

8.2.   Computational results for the FDE system (7.1)

Throughout this section, we compare the analytical findings using a biologically plausible parameter set. Approximate solutions for our fractional-order system are determined using the generalized Adams-Bashforth-Moulton type predictor-corrector scheme ^[49]. For the same initial population $(204,189,130,340)$ and parameters provided in Table 1 and for $\epsilon = 0.9$ the proposed system (7.1) has an endemic equilibrium $E^* = (2265.22, 4154.55, 1370.13, 5628.09)$. With different values of $\epsilon$, time series solutions of the system (7.1) are drawn with the same parameters provided in Table 1 (see Figure 9).

On increasing the order $\epsilon$, initially the recovered population decreases and then increases after certain time (Figure 10). For $\beta = 0.035$, other parameters mentioned in Table 1 and $\epsilon = 0.9$, for the system (7.1), at the endemic equilibrium $E^*$,

Therefore, in the fractional case $0 < \alpha < 1$, the endemic equilibrium point is locally asymptotically unstable. This means that the endemic equilibrium undergoes a transcritical bifurcation at $\beta = 0.035$. For $\beta = 0.005$ with other parameters as mentioned in Table 1, time series solutions of the system (7.1) are drawn in Figure 11. It is observed that the fractional order $\epsilon$ has a significant role in the infection-free state. An increase in the order $\epsilon$ leads to an anthrax infection free ecosystem.

9.   Discussion

In this work, we have extended the anthrax transmission model proposed by Kimathi and Wainaina ^[24], by incorporating a nonlinear ratio-dependent type incidence function. Our study aims to study mathematical models of transmission dynamics of anthrax disease in the sense of ordinary differential equations and fractional order differential equations in the Caputo derivative sense. Positiveness and boundedness of the solutions of the system are discussed. The biologically feasible equilibrium points of the systems are determined, and their existence criteria are determined. The stability of the disease-free equilibrium point is analysed using basic reproduction number $\mathcal{R}_0$.

It is obtained that when $\mathcal{R}_0 < 1$, the disease-free equilibrium $E_0$ is locally asymptotically stable, and, unstable when $\mathcal{R}_0 > 1$. The stability of the endemic equilibrium is explored using the Routh-Hurwitz criterion. Transcritical bifurcations experienced by the endemic equilibrium of the ODE system are discussed using Sotomayor's theorem (Theorems 6.1–6.3). Again, we extended the proposed model (2.2) using Caputo fractional order derivatives in model (7.1). Both the models have the same equilibrium points, so we only discussed the stability criteria of the equilibrium points for the model (7.1). Using some biologically feasible parameter values, we conducted numerical simulations in Matlab software. For the menitioned parameter values the ODE system (2.2) has an endemic equilibrium point

Considering $\mathcal{R}_0$ as a function of contact rate ($\beta$), animal recovery rate ($\alpha$) and disease induced death rate ($\theta$), a 3D region is plotted for which $\mathcal{R}_0 < 1$ (Figure 5). From Figures 6–8, it is observed that with a decrease in the contact rate $\beta$, the endemic equilibrium $E^*$ undergoes transcritical bifurcations at

respectively. In an epidemiological sense, below the threshold contact rate ($\beta^{BP}$), i.e., when $\beta < \beta^{BP}$, the system will have a stable disease-free equilibrium state. To eliminate anthrax from the ecosystem, farmers, ecologists, and policymakers must work towards reducing the contact rate to a level below this crucial threshold $\beta^{BP}$. Also, above the threshold value $\alpha^{BP}$, i.e., $\alpha > \alpha^{BP}$, the system becomes anthrax infection free. To ensure an anthrax-free ecosystem, it is imperative for these stakeholders to attain a recovery rate that surpasses this minimum threshold $\alpha^{BP}$.

This can be achieved either by adopting a proper vaccination policy or by increasing awareness among the farmers. Policymakers need to formulate a well-structured strategy to accomplish this. Moreover, by achieving a minimum threshold disease-induced death rate $\theta^{BP}$, i.e., $\theta > \theta^{BP}$, a stable anthrax infection-free state can be obtained. Thus, the disease-induced death rate $\theta$ has a positive impact on the ecosystem. The higher the rate $\theta$ is the lower the infection will be. The extended fractional order model, (7.1), is studied for different values of the fractional order $\epsilon$ (see Figure 9). Figure 11 shows the impact of the fractional order $\epsilon$ on the recovered population. Figure 11 revealed that an increase in the order $\epsilon$ leads to an anthrax infection-free ecosystem.

10.   Conclusions and future scopes

The present investigation shows that the implementation of a proper vaccination policy can effectively control anthrax transmission in the animal population. We determined the expression for the basic reproduction number $\mathcal{R}_0$ of the infection and obtained that the disease is eradicated from the system when $\mathcal{R}_0 < 1$. We also determined additional threshold values for contact rates, disease-induced mortality rates, and rates of animal recovery, which could influence the feasibility of disease eradication. It has been observed that higher vaccination rates can boost the rate of animal recovery and effectively manage disease-induced mortality, potentially leading to the elimination of the disease from the ecosystem.

In the work of Kimathi et al. ^[24], they employed a bilinear function for disease transmission, although a more realistic scenario would benefit from a nonlinear transmission function. Their observations indicated that maintaining a vaccination rate below a critical threshold allows the anthrax disease to persist. Among the various parameters, recruitment and contact rates were found to be the most influential in determining the basic reproductive number. Rezapour et al. ^[26] expanded upon the model presented in ^[24] by incorporating Caputo-Fabrizio fractional derivatives but maintained the same bilinear transmission function. They have obtained the basic reproduction number to be dependent on several parameters, such as recruitment, recovery, disease-induced death rate, waning immunity rate of vaccinated animals and vaccination rate.

In contrast, in the present model, the basic reproduction number is dependent on contact rate, disease-induced mortality rate, and rate of animal recovery. Additionally, several anthrax transmission models have been introduced in various studies ^{[19,20,21,22,23,25,50,51]}, predominantly centered on the assumption of bilinear disease transmission. Our proposed model differs from the others in terms of a nonlinear ratio-dependent disease transmission function, and we extended its feasibility by introducing Caputo fractional derivatives. The dynamical behavior of the model is studied for different values of the fractional order. Computational results illustrated that for an increase in the fractional order $\epsilon$, the period required for disease eradication in the ecosystem decreases.

The conclusions drawn from the proposed model are completely based on theoretical analysis. Experimental validation will help identify any necessary adjustments to the foundational assumptions. The work in this paper can be extended to study the transmission dynamics under various nonlinear transmission functions with saturation effect. The extensions can be made using both ordinary and fractional-ordered differential equations. Fractional-order systems provide a more accurate and versatile framework for modelling and simulating real-world phenomena compared to traditional integer-order systems. For example, there are the Caputo-Fabrizio fractional derivatives ^[26] and the most recent Atangana-Baleanu-Caputo derivatives ^[44]. The proposed model also can be analyzed under time-delay. Time delay in vaccination modelling accounts for the period between vaccine eligibility and actual vaccination, which can vary due to logistical, behavioral, and immune response factors.

Researchers interested in this area can expand upon our model to investigate how delays in vaccination rates impact the outcomes. Moreover, with more accurate real-world data on anthrax transmission, all the available mathematical models can be explored.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to express their sincere appreciation to Royal Global University, Assam, India for providing financial support through the Institutional Research Project 2022.

Conflict of interest

The authors does not have any conflicts of interest.