Deterministic, stochastic and fractional mathematical approaches applied to AMR

1.   Introduction

Antimicrobial resistance (AMR) is a critical global health challenge that compromises the effectiveness of antibiotics, leading to prolonged infections, increased mortality, and significant economic burdens ^[1]. AMR arises when bacteria develop mechanisms to evade the action of antibiotics, and reversing this resistance is essential to restore the efficacy of treatments ^[2]. The phenomenon of AMR reversal occurs when resistant bacteria lose their resistance mechanisms, often due to fitness costs associated with resistance. In antibiotic-free environments, sensitive bacteria may outcompete resistant ones, leading to the dominance of sensitive populations. Mechanisms such as reversible mutations, the loss of resistant plasmids, and the impact of fitness costs have been well-documented as drivers of AMR reversal ^[3,4,5]. Various mathematical models have described the dynamics of AMR, particularly in bacterial populations, both in vitro and in vivo ^[6,7]. These models provide important insights into how resistance evolves and under what conditions it can be reversed.

To enhance the understanding of AMR reversal dynamics, it is useful to consider analogies with epidemiological processes. In epidemiology, the propagation of an infectious disease is modeled by categorizing individuals either susceptible, infected, or recovered (SIR). Analogously, in AMR dynamics, sensitive bacteria can be conceptualized as susceptible, resistant bacteria as infected, and the reversal of resistance as a recovery process. This analogy provides a framework to examine AMR as an epidemiological process, particularly directly modeling bacterial populations, rather than in infected patients. However, it is crucial to note that while AMR reversal can occasionally align with this analogy, it does not invariably hold true because of the complexity of the resistance mechanisms and bacterial ecology ^[8,9].

Epidemiological models have long been an important tool for mathematical biology, starting with the pioneering work of Kermack and McKendrick in 1927, who proposed the SIR model to describe the spread of infectious diseases within a closed population ^[10]. Subsequent developments have introduced stochastic and fractional approaches to better capture the complexities of real-world dynamics. For instance, stochastic models incorporate randomness to account for population variability and inherent uncertainties in disease transmission ^{[11,12,13,14,15]}. Fractional models extend traditional models by incorporating fractional derivatives, which allow for the inclusion of memory effects and non-local interactions ^[16]. These approaches are particularly valuable for analyzing the stability of antimicrobial resistance (AMR) dynamics. Stability analysis is essential to understand whether resistance will persist, decline, or be eradicated, which is analogous to predicting the endemicity or eradication of infectious diseases. By identifying critical thresholds, such as the basic reproduction number in epidemiological models or resistance thresholds in AMR, these models enhance the capability of predicting the system behavior.

This study addresses AMR reversal using three modeling approaches: deterministic, stochastic, and fractional. To examine the qualitative characteristics of the three models and compare their outcomes, we introduce a simple mathematical model that has been adapted to explore the resistance reversal. To the best of our knowledge, no existing study has analyzed the qualitative properties of a model for AMR reversal from the perspectives of deterministic, stochastic, and fractional differential equations. Rather than proposing a single mathematical model that incorporates all three frameworks, our study focuses on comparing how the stability and long-term behavior of AMR reversal differ across these approaches. This comparative analysis reveals that the outcomes regarding AMR extinction and persistence can vary significantly, depending on the modeling approach employed.

We have divided this work into two main parts. First, we analyze a mathematical model of AMR reversal deterministically, stochastically, and fractionally. Then, the qualitative properties of the model are analyzed. More specifically, we ensured the existence and uniqueness of solutions to these models. In addition, we ensured the invariance of the solutions within a biologically relevant set. Moreover, we demonstrated the stability of the trivial solutions in these models. Finally, we provide the conditions for disease extinction, disease persistence, and stationary measures in some models. In the second part, we used assumed data from E. coli exposed to colistin to numerically validate the theoretical results.

2.   Mathematical model formulation

Let $S(t)$ denote the number of sensitive bacteria at time $t$, and $R(t)$ represent the number of resistant bacteria at time $t$. To formulate our mathematical model, we make the following basic assumptions, keeping in mind that the primary focus is on exploring the deterministic, stochastic, and fractional modeling approaches, thereby using AMR dynamics as an illustrative example, rather than going into biological details.

- We assume a turnover rate $\mu$ which represents that both sensitive and resistant bacteria die and are replaced at a proportional rate, while ignoring fitness costs of resistance.

- Let $\beta$ denotes the rate at which sensitive bacteria acquire plasmids through conjugation. Resistant bacteria lose their resistance genes and revert to a sensitive state at a rate $\gamma$.

- Sensitive bacteria reduce the rate of acquiring resistance by selective pressure as the population of the resistant bacteria increases, which is modeled by a functional response $g(R) = \beta\dfrac{R}{1+\epsilon R}$, where $\epsilon$ is a known constant called a saturation rate.

Consequently, our mathematical model is described by the following deterministic equations:

In this work, the term $\mathbb{R}_{0}^{+}$ represents non-negative real numbers. The equations in (2.1) are coupled by the assumption $N = S(t)+R(t)$ for all $t \in \mathbb{R}_{0}^{+}$. Thus, we can study one equation in terms of the resistant bacteria $R(t)$. More precisely, we will study the following deterministc equation:

In the stochastic model, we perturb the rate $\beta$ as follows:

where $\sigma$ is the rate of the stochastic perturbation and $B$ is a Brownian motion in a probability space $(\Omega, \mathcal{F}, \mathbb{P})$. Therefore, the stochastic differential equation takes the following form:

This type of perturbation appears in ^[12,13].

Finally, we introduce the fractional model, which we represent with the following fractional differential equation:

where $\dfrac{d^{\alpha}}{dt^{\alpha}}$ is the Caputo derivative of order $\alpha$, with $\alpha \in (0, 1)$. This type of fractional model equation can be found in ^[16,17].

In this work, we define the threshold parameters that characterize the qualitative properties of (2.2)–(2.4). The ideas for proving the qualitative properties of (2.3) were adapted from ^[18] and are presented here in terms of the stochastic threshold parameter.

3.   Theoretical results

3.1.   Defining thresholds

In order to facilitate the mathematical calculations of the three models (2.2)–(2.4) and make the results easier to understand, we define certain thresholds, similar to the basic reproduction numbers in epidemiological models, and introduce Lyapunov operators. These thresholds and operators guide the stability and persistence of plasmid-carrying bacteria within a population in different modeling contexts.

3.2.   The deterministic approach

This section explores the qualitative properties of the deterministic model (2.2). Let us note that if

then the equilibrium points of (2.2) are as follows:

We start by showing that $(0, N)$ is an invariant set. More precisely, if $R_0 \in (0, N)$, then there exists a unique global solution $R(t; R_0)$ of (2.2) such that $R(t; R_0) \in (0, N)$ for all $t \in \mathbb{R}_{0}^{+}$. For this end, we use Lyapunov operator techniques, which are typically applied in the case of SDEs, however these techniques can also be used in the deterministic case, where the stochastic noise is zero (see ^[19]).

Theorem 3.1 (Invariance). For any $R_0 \in (0, N)$, there exists a unique global solution of (2.2) invariant in $(0, N).$

Proof. Let $R_0 \in (0, N)$. Consider $D : = (0, N)$ and let $b(R) : = \beta \dfrac{R}{1+\epsilon R}(N-R) - (\gamma + \mu)R$, for $R \in \mathbb{R}^{+}$. Note that $b$ is locally Lipschitz because it is continuously differentiable. Now, consider $V : D \to \mathbb{R}$ defined by: $V(R) : = \dfrac{1}{R} + \dfrac{1}{N-R}$ and $D_n : = \Big(\dfrac{1}{n}, N - \dfrac{1}{n}\Big)$, for $n \in \mathbb{N}$. It is clear that $V \in C^1(D; \mathbb{R}_{0}^{+})$ and that $\{D_n\}_{n \in \mathbb{N}} \subset D$, $D_n$ are domains, $\overline{D_n} \subset D$, $D_n \uparrow D$ such that $V_n : = \inf_{R \notin D_n} {V(R)} \to \infty, \quad \text{as } n \to \infty.$ On the other hand, given $R \in D$, we have that:

Let us take $c : = \max \{\gamma + \mu, \beta N \}$. Thus,

Then, there exists a unique solution $R(\cdot\; R_0)$ of (2.2) in $\mathbb{R}_{0}^{+}$ such that $R(t; R_0) \in (0, N), \text{for all}\ t \in \mathbb{R}_{0}^{+}$. □

Now, we prove the asymptotic stability of the equilibrium points of the deterministic model (2.2) defined in (3.1) (see Definition.2). Intuitively, the parameter $\mathcal{K}_{0}^{d}$ acts as an indicator of the system's dynamics, as we will see in the following theorem.

Theorem 3.2 (Asymptotic stability). If $\mathcal{K}_{0}^{d} < 1$, then $R = 0$ is an asymptotically stable equilibrium point of (2.2). Moreover, if $\mathcal{K}_{0}^{d} > 1$, then $R = \xi_{d}$ is an asymptotically stable equilibrium point of (2.2).

Proof:

First, let's see that $R = 0$ is an asymptotically stable equilibrium point of (2.2). Consider $\delta > 0$ small enough and let $D = (-\delta, \delta)$. Note that

for $R \in D$. It is clear that $b$ is differentiable and

for $R \in D$. Thus, $b'(0) = \beta N - (\gamma + \mu) < 0$. By the Lyapunov linearization theorem for the deterministic case, we have that $R = 0$ is an asymptotically stable equilibrium point of (2.2).

Now, let's see that $R = \xi_{d}$ is an asymptotically stable equilibrium point of (2.2). Let $D : = (0, N)$ and consider $V: D \to \mathbb{R}$ defined by the following:

Clearly, $V \in C^1(D; \mathbb{R})$ and $V$ is positive definite at $R = \xi_{d}$. In addition, given $R \in D$, we have the following:

Thus, $L_d V(R) < 0$, for $R \neq \xi_{d}$. Then, by the Lyapunov stability theorem (Theorem 4.7 in (^[20]) for equilibrium points in the deterministic case, we have that $R = \xi_{d}$ is an asymptotically stable equilibrium point of (2.2). □

In the following theorem, we show that if the deterministic threshold condition $\mathcal{K}_{0}^{d} < 1$ is fulfilled, then the population of resistant bacteria will ultimately die out. Of course, this result is stronger than the previous one when $\mathcal{K}_{0}^{d} < 1$

Theorem 3.3 (Extinction). If $\mathcal{K}_{0}^{d} < 1$, then for any $R_0 \in (0, N)$, we have the following:

Proof:

Let's show that

We already know that there exists a unique solution $R(\cdot; R_0)$ in $\mathbb{R}_{0}^{+}$ of (2.2) such that $R(t; R_0) \in (0, N)$ for all $t \in \mathbb{R}_{0}^{+}$. Let us denote $R(\cdot; R_0)$ by $R$ and $D : = (0, N)$.

Now, consider the function $V : D \to \mathbb{R}$ defined by $V(R) : = \ln(R)$. It is clear that $V \in C^{1}(D; \mathbb{R})$.

Let $t > 0$. By the chain rule and integrating from 0 to $t$, we have the following:

In other words,

That is,

Then, dividing by $t$, we obtain the following:

Therefore,

Thus,

□

In the following theorem, we state that when the deterministic threshold $\mathcal{K}_{0}^{d} > 1$, the population of resistant bacteria will persist over time. In other words, the system described by (2.2) converges to the non-trivial equilibrium point $\xi_{d}$ defined in (3.1) as $t \to \infty$.

Theorem 3.4 (Persistence). If $\mathcal{K}_{0}^{d} > 1$, then there exists $\xi_{d} \in (0, N)$ such that for any $R_0 \in (0, N)$,

where $\xi_{d}$ is the equilibrium point defined in (3.1).

Proof:

Let $D : = (0, N)$ and consider $V : D \to \mathbb{R}$ defined by the following:

Clearly, $V \in C^{1}(D; \mathbb{R})$, $V$ is positive definite at $R = \xi_{d}$, and $L_{d} V(R) < 0$ for $R \neq \xi_{d}$. Additionally, we have: $V(R) \to \infty$ as $R \to 0^{+}$ and as $R \to N^{-}$. By the Barbashin-Krasovskii theorem (Theorem 4.2 in ^[20]), $R = \xi_{d}$ is a globally asymptotically stable equilibrium point of (2.2), meaning that $\lim\limits_{t \to \infty} R(t; R_0) = \xi_{d}$ for any $R_{0} \in (0, N)$. □

3.3.   The stochastic approach

In this section, we present a series of results that involve the numerical value $\mathcal{K}_0^s$. The proofs are direct modifications of the proofs in ^[18]. In ^[18], a family of stochastic differential equations is considered, depending on a general function $g(R)$ under certain conditions. In this article, we work with the function $g(R) = \beta \frac{R}{1 + \epsilon R}$ as a particular case of ^[18]. We have chosen to provide intuitive commentary on the results and to reference the work in ^[18] for the full proofs.

Now, in this section explores the qualitative properties of the stochastic model (2.3). Let us note that

and

where the only equilibrium point of (2.3) is $R = 0.$

The following result is fundamental in the stochastic modeling of population systems, as it guarantees that, the stochastic system has a unique global solution that remains within the interval $(0, N)$ over time regardless of the initial condition $R_0 \in (0, N)$. This implies that the model is robust to random fluctuations and ensures that the stochastic dynamics of the resistant bacterial population neither becomes extinct in finite time nor exceeds the maximum capacity $N$. In other words, the system is able to realistically model a sustained population range in the presence of randomness.

Theorem 3.5 (Invariance). For any $R_0 \in (0, N)$, there exists a unique global solution of the stochastic system (2.3) that remains invariant in $(0, N)$.

Proof:

See Theorem 2.1 in ^[18]. □

The following result states that if $\mathcal{K}_0^s < 1-\dfrac{\sigma^2 N^2}{(\gamma+\mu)}$, then, starting near extinction ($R = 0$), the resistant bacterial population converges to extinction with high probability (see Definition.4).

Theorem 3.6 (Asymptotic stability). If $\mathcal{K}_{0}^{s} < 1-\dfrac{\sigma^2 N^2}{(\gamma+\mu)}$, then $R = 0$ is an asymptotically stable equilibrium point in the probability of (2.3).

Proof:

See Theorem 2.2 in ^[18]. □

In the following theorem, we present sufficient conditions for the population of resistant bacteria to converge to extinction as $t \to \infty$. One condition that seems natural after the previous results is that $\mathcal{K}_{0}^{s} < 1$; this ensures that if we start close to the equilibrium point $R = 0$, then the system will converge to $0$ with a high probability. Another condition is needed to guarantee convergence to extinction starting from any point. This condition arises from the utilized Lyapunov function, so it may be subject to refinement. The condition is that $\sigma^2 N \leq \beta$; this suggests a certain restriction on $\beta$, in the sense that this parameter cannot be taken as arbitrarily small.

Theorem 3.7 (Extinction). If $\mathcal{K}_{0}^{s} < 1$ and $\sigma^2 N \leq \beta$, then for any $R_0 \in (0, N)$, we have the following:

Proof:

We already know that there exists a unique solution $R(\cdot, \cdot; R_0)$ in $\mathbb{R}_{0}^{+}$ of (2.3) such that $R(\cdot, t; R_0) \in (0, N),$ for all $t \in \mathbb{R}_{0}^{+},$ $\mathbb{P}$-a.s. Let us denote $R(\cdot, \cdot; R_0)$ by $R$ and define $D : = (0, N).$

Now, consider the functions $V: D \to \mathbb{R}$ defined by $V(R) : = \ln(R)$ and $\varphi: D \to \mathbb{R}$ defined by $\varphi(R) : = \dfrac{N - R}{1 + R}.$ It is clear that $V \in C^{2}(D; \mathbb{R})$ and $\varphi(R) \le N,$ for all $R \in (0, N).$

Let $t > 0.$ By Itô's formula (Theorem 6.2 in ^[21]), we have the following:

In other words,

Now, since $N \le \dfrac{\beta}{\sigma^2}$ and the function $f(u) : = \beta u - (\gamma + \mu) - \dfrac{1}{2} \sigma^{2} u^2$ is increasing for $u \le \dfrac{\beta}{\sigma^2},$ then $f(\varphi(R(s))) \le f(N),$ for all $s \in [0, t].$ In other words,

for all $s \in [0, t].$ The rest of the proof is similar to the proof of Theorem 2.3 in ^[18].

□

Before stating the theorem on the persistence of the bacterial population, we would like to provide a brief intuition of the result. In this context, the quantity $R(\cdot, t; R_0)$ represents the size of the resistant bacterial population over time $t$, starting from an initial condition $R_0$. A key objective is to understand whether this population tends to persist over time or, on the contrary, tends to vanish.

The theorem below indicates that, under certain conditions on the system parameters ($\mathcal{K}_{0}^{s} > 1$), there exists a persistence level $\xi_{s}$ such that the resistant bacterial population will repeatedly reach this level over time. In other words, despite the stochastic variability in the model, we can expect the population not only to remain above extinction thresholds but also to oscillate around $\xi_s > 0$. In summary, this implies a tendency towards persistence rather than extinction, thus providing a measure of stability for the resistant bacterial population under the effect of random factors.

The proof of the following theorem has subtle differences compared to the proof of Theorem 2.4 in ^[18]. We will highlight these differences in the demonstration, and when the tools are fully analogous, we will refer to ^[18].

Theorem 3.8. If $\mathcal{K}_{0}^{s} > 1$, then there exists $\xi_{s} \in (0, N)$ such that for any $R_0 \in (0, N)$, the following holds:

and

where $\xi_{s}$ is given by

In other words, the solution process $R(\cdot, \cdot\; R_0)$ visits the level $\xi_{s}$ infinitely often for any $R_0 \in (0, N)$.

Proof:

Let $D: = (0, N)$, and consider $V: D\to \mathbb{R}$ defined by $V(R): = \ln(R).$ It is clear that $V\in C^{2}(D; \mathbb{R}).$ Note that the stochastic Lyapunov operator applied to $V$ is given by the following:

Rewriting, we have the following:

where $\varphi(R): = \dfrac{N-R}{1+\epsilon R}.$ Consider the following quadratic function:

Note that $f(0) = -(\gamma+\mu)$, and since $\mathcal{K}_{0}^{s} > 1$, then $f(N) = \beta N -(\gamma+\mu)-\dfrac{1}{2}{\sigma^{2}}N^2 > 0.$ Therefore, there exists $\eta\in (0, N)$ such that $f(\eta) = 0.$ Now, there exists a unique $\xi_{s}\in(0, N)$ such that $\varphi(\xi_{s}) = \eta.$ Thus, $0 = f(\eta) = f(\varphi(\xi_{s})) = L_{s}V(\xi_{s}).$

Case 1. If $\dfrac{\beta }{\sigma^{2}} < N$, since $\varphi'(R) < 0,$ for all $R\in(0, N),$ there exists $m\in (0, N)$ such that:

- $L_{s}V$ is increasing in $(0, m)$ and $L_{s}V(R) > 0,$ for all $R\in(0, m).$

- $L_{s}V$ is decreasing in $(m, \xi_{s})$ and $L_{s}V(R) > 0,$ for all $R\in(m, \xi_{s}).$

- $L_{s}V$ is decreasing in $(\xi_{s}, N)$ and $L_{s}V(R) < 0,$ for all $R\in(\xi_{s}, N),$

where $L_{s}V(m)$ is the maximum value of $L_{s}V$ in $(0, N).$

Case 2. If $\dfrac{\beta }{\sigma^{2}}\ge N,$ since $\varphi'(R) < 0,$ then for all $R\in(0, N),$ $L_{s}V(0^{+}): = \lim\limits_{R\to 0^{+}}L_{s}V(R) = \beta N -(\gamma+\mu)-\dfrac{1}{2}{\sigma^{2}}N^2 > 0$ and $L_{s}V(N^{-}): = \lim\limits_{R\to N^{-}} L_{s}V(R) = -(\gamma +\mu) < 0$ :

- $L_{s}V$ is decreasing in $(0, \xi_{s})$ and $L_{s}V(R) > 0,$ for all $R\in(0, \xi_{s}).$

- $L_{s}V$ is decreasing in $(\xi_{s}, N)$ and $L_{s}V(R) < 0,$ for all $R\in(\xi_{s}, N).$

The rest of the proof is similar to proof of Theorem 2.4 of ^[18]. □

The following theorem establishes the stationary distribution of the solution $R$ of the stochastic system (2.3). This allows us to ensure that the stochastic model exhibits a stable and predictable behavior over the long term under certain conditions. Therefore, as the system evolves, the distribution of resistant bacteria will stabilize to a unique stationary distribution.

Theorem 3.9 (Stationary distribution). If $\mathcal{K}_{0}^{s} > 1$, then $\mu_{\infty}(\cdot)$ is the unique stationary distribution associated with the solution $R$ of (2.3).

Proof:

See Theorem 2.5 in ^[18] □.

3.4.   The fractionary approach

Finally, in this section explores the qualitative properties of the fractionary model (2.4). Let us note that if $\beta \dfrac{R}{1+\epsilon R}(N - R) - (\gamma + \mu)R = 0$, then the equilibrium points of (2.4) are as follows:

We start showing that $(0, \mathcal{K}_{0}^{f})$ is an invariant set. In this case, the invariance is more restrictive, since the fractional model has more limitations.

Theorem 3.10 (Invariance). For any $R_0 \in (0, \mathcal{K}_{0}^{f})$, there exists a unique global solution of (2.4) invariant in $(0, \mathcal{K}_{0}^{f}).$

Proof:

Let $R_0 \in \big(0, \mathcal{K}_{0}^{f}\big)$ and define $b(R): = \beta\dfrac{R}{1+ \epsilon R}(N-R)-(\gamma+\mu)R$, for $R \in \mathbb{R}^{+}.$

First, let us show that there exists a unique global positive solution of (2.4). Since $b$ is locally Lipschitz, there exists a unique local maximal solution $R(\cdot; R_{0}, \alpha):[0, T_{f}) \to \mathbb{R}$, where $T_{f}$ is the fractional explosion time. Since $R_{0} \in \mathbb{R}^{+}$, there exists $n_0 \in \mathbb{N}$ such that $\dfrac{1}{n_{0}} < R_{0}.$

Now, we define the following sequence of stopping times:

It is clear that $\{T_{n}\}_{n \ge n_0}$ is an increasing sequence. Moreover,

Let us show that $T_{\infty} = \infty.$ Suppose that $T_{\infty} < \infty$; then, there exists $T \in \mathbb{R}^{+}$ such that $T_{\infty} \le T.$

Moreover, there exists $n_1 \in \mathbb{N},$ with $n_1 > n_0$ such that $T_{n} \le T.$ On the other hand, consider the function $V: \mathbb{R}^{+} \to \mathbb{R}$ defined by the following:

Given $R \in \mathbb{R}^{+},$ by property Lemma 3.1 of ^[22], we have the following:

Let $c : = \beta N + \beta + \gamma + \mu.$ Thus, for all $R \in \mathbb{R}^{+}$, $L_{f}V(R) \le c.$ Applying the fractional integral to both sides from 0 to $T \wedge T_n$, for $n \in \mathbb{N}$, with $n > n_0$, we have the following:

Now, for all $n \in \mathbb{N}, n > n_1,$ we have $T \wedge T_n = T_n$ and $V(R(T \wedge T_n)) = V(R(T_n)).$

Note that:

hence,

Then, as $n \to \infty$, we have the following:

This is a contradiction, since the expression $V(R_0) + \dfrac{cT^{\alpha}}{\Gamma(\alpha+1)}$ is finite.

Now, let us show that $R(t; R_{0}, \alpha) \in (0, N),$ for all $t \in \mathbb{R}_{0}^{+}.$ Let us denote this solution by $R$.

Let $t \in \mathbb{R}^{+}.$ Note that:

Solving the inequality, we have the following:

where $E_{\alpha}$ is the Mittag-Leffler function of parameter $\alpha$ defined in ^[23]. Since $0 \le E_{\alpha}(-(\gamma+\mu)t^{\alpha}) \le 1$ and $R_{0} - \mathcal{K}_{0}^{f} < 0,$ then

□

Theorem 3.11 (Asymptotic stability). If $\mathcal{K}_{0}^{f} < 1$, then $R = 0$ is an asymptotically stable equilibrium point of (2.4). Moreover, if $\mathcal{K}_{0}^{f} > 1$, then $R = \xi_{f}$ is an asymptotically stable equilibrium point of (2.4).

Proof:

First, let us show that $R = 0$ is an asymptotically stable equilibrium point of (2.4). Consider $\delta > 0$ sufficiently small and let $D = (-\delta, \delta).$ Note that

for $R \in D.$ It is clear that $b$ is differentiable and

for $R \in D.$ Thus, $b'(0) = \beta N - (\gamma + \mu) < 0.$ By the fractional Lyapunov linearization theorem (Theorem 2 in ^[24]), we have that $R = 0$ is an asymptotically stable equilibrium point of (2.4).

Now, let us show that $R = \xi_{f}$ is an asymptotically stable equilibrium point of (2.4). Consider $\delta > 0$ sufficiently small and let $D = (\xi_{f} - \delta, \xi_{f} + \delta).$ Thus,

By the fractional Lyapunov linearization theorem (Theorem 2 in ^[24]), we have that $R = \xi_{f}$ is an asymptotically stable equilibrium point of (2.4). □

Theorem 3.12 (Extinction). If $\mathcal{K}_{0}^{f} < 1$, then for any $R_0 \in (0, \mathcal{K}_{0}^{f})$, we have the following

Proof:

Let $D : = [0, \mathcal{K}_{0}^{f}]$ and consider $V: D \to \mathbb{R}$ defined by the following:

Clearly, $V \in C^{1}(D; \mathbb{R})$ and $V$ are positive definite at $R = 0$. Additionally, given $R \in D$, using the Lemma 2.1 of ^[22], we have the following:

Thus, $L_{f}V(R) \le 0$ for all $R \in D;$ and moreover,

Then, by the fractional LaSalle theorem (Lemma 4.6 in ^[25]), we have that $R = 0$ is a globally asymptotically stable equilibrium point of (2.4). Then, $\lim\limits_{t \to \infty} R(t; R_0, \alpha) = 0$ for any $R_{0} \in (0, \mathcal{K}_{0}^{f})$.

□

Theorem 3.13 (Persistence). If $\mathcal{K}_{0}^{f} > 1$ and $\dfrac{\beta N-\gamma-\mu}{\beta N}-\dfrac{\beta}{\gamma+\mu} < \epsilon < 1$, then there exists $\xi_{f} \in (0, \mathcal{K}_{0}^{f})$ such that for any $R_0 \in (0, \mathcal{K}_{0}^{f})$,

where $\xi_{f}$ is the equilibrium point defined in (3.1).

Proof:

Let $D : = [0, \mathcal{K}_{0}^{f}]$ and consider $V: D \to \mathbb{R}$ defined by the following:

Clearly, $V \in C^{1}(D - \{0\}; \mathbb{R})$ and $V$ is positive definite at $R = \xi_{f}$. Additionally, given $R \in D - \{0\}$, using the Lemma 3.1 of ^[22], we have the following:

Thus, $L_{f}V(R) \le 0$ for all $R \in D - \{0\};$ and moreover,

Then, by the fractional LaSalle theorem (Lemma 4.6 in ^[25]), we have that $R = \xi_{f}$ is a globally asymptotically stable equilibrium point of (2.4). Then, $\lim\limits_{t \to \infty} R(t; R_0, \alpha) = \xi_{f}$ for any $R_{0} \in (0, \mathcal{K}_{0}^{f})$.

□

4.   Numerical simulations

Our three proposed approaches, (2.2)–(2.4), can be numerically validated, using data from various real-world epidemiological phenomena that align with our hypotheses. However, we have represented the phenomenon of AMR, specifically focusing on the resistance of Escherichia coli to colistin. Although this particular case serves as a representative study to explore our mathematical framework, it does not fully capture the complexities inherent in the AMR dynamics.

Colistin has been the final option for treating multidrug-resistant gram-negative bacteria (MDR-GNB), and it is also extensively used in veterinary practice ^[26]. The recent discovery of plasmid-mediated colistin resistance, exemplified by the mcr-1 gene in E. coli, has raised concerns regarding its use in food-producing animals and its potential to accelerate the spread of resistance. Studies have called for a reassessment of colistin usage and dosing in animal husbandry to safeguard its effectiveness in human healthcare ^[27,28]. The reversibility of AMR, especially for colistin, has been demonstrated in laboratory experiments using E. coli that carry the mcr-1 gene ^[29]. Under colistin-free conditions, there is a significant reduction in antibiotic resistance genes due to the elimination of MDR plasmids. This suggests that the high fitness costs associated with mobile genetic elements, such as plasmids, make resistance unstable, and consequently, strict antibiotic control can help reverse resistance driven by these genes ^[29]. Other studies have also demonstrated that E. coli adapts to prolonged antimicrobial exposure through the genetic regulation of porins and efflux pumps, which are crucial for bacterial resistance. Therefore, both transcriptional and post-translational regulation of membrane proteins are essential for the physiological adaptation of gram-negative bacteria to antibiotic stress ^[30].

We begin with a population of $N = 10^6$ bacteria and consider two initial conditions: $R(0) = N-1$ (indicating that almost all bacteria are resistant) to demonstrate extinction (or clearance of bacteria) and $R(0) = 1$ (where almost all bacteria are susceptible) to illustrate persistence. Under stressful conditions, the turnover rate of E. coli ($\mu$) is typically significantly lower than that under optimal conditions because of the adverse factors that affect bacterial growth. In our study, we assumed a natural turnover rate of $\mu = 0.1$, which imply that the bacterial population doubles or clears approximately every 10 days. The population-dependent rate of R-plasmid transfer through conjugation $\beta$ for E. coli under colistin exposure was derived from observational data in ^[29], where it was suggested that a susceptible bacterium acquires an R-plasmid approximately every two days. The rate of R-plasmid loss was estimated to be within the range of $[0, 2]$, thus indicating that resistant bacteria may lose plasmids every 12 hours ($\gamma = 2$), or never lose them ($\gamma = 0$). The saturation parameter $1/\epsilon$, which represents the number of resistant bacteria at which plasmid transfer becomes less efficient due to competition between sensitive and resistant bacteria, was assumed to be equal to the total population $N$ ($\epsilon = 1/N = 10^{-6}$).

Finally, considering that AMR largely stems from an evolutionary process in which bacteria gradually adapt to their environment under antibiotic stress ^[31], it is possible to capture AMR dynamics by highlighting the influence of historical exposure in the presence of the fitness costs associated with plasmids. This suggests that the resistance characteristics of E. coli are not only dependent on current conditions, but are also heavily influenced by their previous exposure to antibiotics (selective pressure) and the presence of plasmids. For $\alpha \in (0, 1)$, when $\alpha$ is closer to zero, the influence of the past states becomes more pronounced, indicating that the previous antibiotic pressures significantly contribute to the current state of resistance. Table 2 shows the values used for the validation experiments.

Figure 1 illustrates the extinction scenarios for antibiotic-resistant bacteria for the three different modeling approaches (deterministic, fractional, and stochastic). These scenarios were evaluated for varying values of the plasmid loss rate $\gamma$. In each of these plots, it can be observed that the resistant bacterial populations tended to decline over time, provided that the threshold $\mathcal{K}_0$ parameters remained below 1. In the deterministic and fractional models, the extinction of resistant bacteria occurs at a faster rate as $\gamma$ increases from 1 to 2. This behavior suggests a strong correlation between the plasmid loss rate and the speed of bacterial extinction, with higher values of $\gamma$ leading to a more rapid approach toward zero population levels for resistant strains. While exhibiting a more erratic pattern due to inherent randomness, the stochastic model, also shows a trend toward the extinction of resistant bacteria when the critical threshold $\mathcal{K}_{0}^{s}$ remains below 1. This threshold condition indicates that the resistant bacteria cannot sustain their population over time, the reinforcing the consistency of the extinction behavior across all three modeling approaches under these conditions.

Figure 2 illustrates the persistence scenario for resistant bacteria using the three modeling approaches. This condition, $\mathcal{K}_0 > 1$, signifies that the resistant population can sustain itself over time, thereby avoiding a clearance. In each model, as the plasmid loss rate $\gamma$ increased within the range of 0–1, resistant bacteria exhibited a stable persistence level. This indicates that lower values of $\gamma$ promoted plasmid retention, thereby supporting the maintenance of antibiotic resistance within the bacterial population. Notably, for values of $\gamma$ near one, the population of resistant bacteria approached a steady-state level early on and remained nearly constant over extended time periods. This steady-state persistence across all three models suggests that, when $\mathcal{K}_0 > 1$ and $\gamma \leq 1$, the conditions facilitate a resilient resistant population capable of persisting over time.

Figure 3 presents a comparison between the deterministic and stochastic models for different values of the white noise intensity $\sigma$. In these simulations, the plasmid loss rate was set to $\gamma = 2$ with an initial condition of $R(0) = N - 1$ for the extinction case, and $\gamma = 0$ with $R(0) = 1$ for the persistence case. The plots reveal how the stochastic noise intensity $\sigma$ influences the extinction and persistence behaviors in the stochastic model compared with the deterministic model. As the value of $\sigma$ increases, the stochastic model exhibits greater fluctuations, reflecting the impact of random perturbations on bacterial population dynamics. This variability contrasts with the deterministic model, which follows a smooth trajectory toward either extinction or persistence, depending on the initial conditions and plasmid loss rate. Specifically, under the extinction condition ($\gamma = 2$, $R(0) = N - 1$), the stochastic model shows a more erratic decline in resistant bacteria as $\sigma$ increases, thus suggesting that a higher noise intensity accelerates fluctuations toward extinction. Conversely, in the persistence case ($\gamma = 0$, $R(0) = 1$), the stochastic model demonstrated sustained fluctuations around a stable population level. Here, increasing the noise intensity $\sigma$ introduces variability, but does not prevent the resistant population from persisting. These results indicate that stochastic fluctuations modulated by $\sigma$ can significantly alter the dynamics of the resistant bacteria. Higher noise levels amplify random effects, thus leading to more pronounced deviations from the deterministic model outcomes, particularly affecting the stability of bacterial populations under both extinction and persistence conditions. In the stochastic case, the expression:

shows that $\mathcal{K}_{0}^{s}$ depends on $\sigma^2$ when the rest of the parameters are fixed. Note that when $\sigma^2$ is equal to zero, we have $\mathcal{K}_{0}^{s} = \mathcal{K}_{0}^{d}$, and $\mathcal{K}_{0}^{s}$ decreases as a function of $\sigma^2$. For this reason, it is clear that a deterministic extinction implies a stochastic extinction, and a stochastic persistence implies a deterministic persistence.

In Figure 4, the fractional model is analyzed for different values of the fractional derivative order $\alpha$. The simulations were conducted with $\gamma = 1.5$ and an initial condition of $R(0) = N - 1$ for the extinction scenario, and with $\gamma = 0.2$ and $R(0) = 1$ for the persistence case. The results show that as the order $\alpha$ of the fractional derivative increases, both the extinction and persistence behaviors in the fractional model become more pronounced. Specifically, in the extinction case, larger values of $\alpha$ led to a more rapid approach to zero for the resistant bacterial population. This suggests that a higher order of the fractional derivative accelerates the population decline by amplifying the memory effects inherent in the fractional systems. This reflects how the system's extended memory dynamically integrates past interactions, such as historical antibiotic usage or environmental pressures. These enhanced memory effects make the population more responsive to cumulative changes, thus intensifying the rate of decline as adverse factors compound over time. Conversely, in the persistence case, increasing $\alpha$ results in a faster stabilization at the persistence level. This stabilization implies that higher values of $\alpha$ facilitate more efficient regulation of the population dynamics, thus allowing the resistant bacteria to reach a steady state more quickly. This behavior is characteristic of fractional models, where memory effects become more pronounced as $\alpha$ increases, thus leading to a faster stabilization in response to lower plasmid loss rates.

The histograms in Figure 5 provide an approximation of the stationary distribution associated with the stochastic model (2.3). This stationary distribution corresponds to the probability measure induced by $R_{\infty}: = \lim_{t \to \infty}R(t; \cdot)$. Notably, this probability measure is independent of the initial condition. Moreover, if the system begins ($t = 0$) with this probability measure, it then remains the same for all $t > 0$, thus reflecting its invariance over time.

5.   Discussion

This study addressed antimicrobial resistance (AMR) reversal using three distinct modeling approaches: deterministic, stochastic, and fractional. By proposing a simplified epidemiological model hypothetically adapted to represent this phenomenon, this study aimed to analyze the qualitative properties of these approaches and compare their outcomes. Although the model provided an interesting framework for exploring AMR reversal, it did not incorporate all the complex dynamics that underlie the resistance reversal process. This mathematical exercise offered useful outputs to the strengths and limitations of each modeling perspective in capturing the phenomenon.

The results obtained from the analysis in this study showed that all three models (deterministic, stochastic, and fractional) are well posed, satisfying the conditions of existence, uniqueness, and invariance within its the domain, and $(0, N)$ refers to the number of resistant bacteria contained in the population, thus providing the model with coherence and biological relevance. For each of these models, three thresholds $\mathcal{K}_{0}^{d}, \mathcal{K}_{0}^{s},$ and $\mathcal{K}_{0}^{f}$, were identified as structurally similar to the basic reproduction number characteristic of epidemiological models. It was demonstrated that, if these thresholds are less than one, then the asymptotic stability of the trivial equilibrium point $R = 0$ is guaranteed. In the stochastic model, the condition $\mathcal{K}_{0}^{s} < 1$ also implies that the trivial equilibrium is stable in probability. In contrast, when these thresholds are greater than one, both the deterministic and fractional models guarantee the asymptotic stability of the nontrivial equilibrium point.

In relation to the phenomenon of extinction (bacterial clearance), it should be noted that the modeling in the deterministic and fractional cases is similar. The stochastic model was more detailed and presented more scenarios. The existence of the stationary measure in the stochastic case provides a probabilistic framework that allows for a better understanding of the evolution of AMR.

To validate the theoretical models, our numerical experiments focused on the AMR of Escherichia coli to colistin, a last-resort antibiotic for infections caused by multidrug-resistant gram-negative bacteria (MDR-GNB). By analyzing the extinction and persistence of resistant bacteria in relation to different plasmid loss rates ($\gamma$) and the intensity of stochastic noise ($\sigma$), we observed that the threshold parameter $\mathcal{K}_0^{\cdot}$ is critical to determine whether a resistant population goes extinct or persists. When $\mathcal{K}_0^{\cdot} < 1$ and $\gamma$ increase, plasmid loss becomes more frequent and accelerates the extinction of resistant bacteria, as the cost of maintaining the plasmid outweighs the adaptive benefits. In contrast, under persistent conditions ($\mathcal{K}_0^{\cdot} > 1$ and $\gamma \leq 1$), bacteria manage to maintain the plasmid, thus showing early stabilization in both deterministic and fractional models, as well as a persistent level of resistance in the stochastic model, even with small fluctuations. This suggests that plasmid retention is key to maintaining resistance in an environment with continuous antibiotic exposure, where a lower plasmid loss rate allows resistance to persist over time.

Furthermore, by incorporating the memory of prior exposure to antibiotics (evolutionary pressure) through the fractional approach, represented by the parameter $\alpha$, it was observed that higher values of $\alpha$ intensified the response of the bacterial population. In the extinction case, high values of $\alpha$ facilitate a more rapid decline in the resistant population due to memory effects in the system, which reinforces the response to plasmid loss pressure. In the case of persistence, high values of $\alpha$ allow for a more efficient and rapid stabilization of the resistant population, thus demonstrating that the history of antibiotic exposure and the regulation of resistance genes (such as the reduction of porins or increased expression of efflux pumps) have an accumulating and relevant effect on bacterial adaptation. These findings highlight the importance of biological memory in AMR, as the adaptation of E. coli to antibiotic stress is not solely dependent on the current conditions, but is strongly influenced by the history of prior exposure, thus implying that antibiotic control strategies could leverage this effect to reverse resistance in bacterial populations.

The results of this study provide an important basis for the qualitative analysis of SIS-type mathematical models in general and set the stage for future research aimed at applying these findings to broader and more complex epidemiological studies. Despite the relevance of these results, this study has some limitations. First, one of the most important limitations is the assumption of a constant population, whereas in a realistic scenario, these bacteria can grow or decline depending on the environmental conditions and other factors. Second, the function $g(R) = \dfrac{R}{1 + \epsilon R}$ is a functional response adapted from the literature and may not specifically represent the transfer (acquisition and loss) of plasmids. Nevertheless, it was proposed in ^[18] a general formulation for the stochastic case.

However, without significant challenges, similar results can be obtained with function $h$, which possess the properties determined in ^[18]. Third, the simplification of the model regarding the environmental conditions and population structure. In our approach, we assumed a homogeneous bacterial population and constant external conditions such as nutrient availability and antibiotic pressure. However, in real-world scenarios, factors such as environmental heterogeneity, competition between different bacterial strains, and fluctuations in antibiotic exposure can significantly influence the resistance dynamics. Including these factors in the model could offer a more detailed perspective and allow for a better representation of the complex interactions in bacterial systems exposed to antibiotics.

This study raises open questions for future research. A key area is the explicit calculation of the mean and variance of the stationary measure in the stochastic model (2.3). The other is to identify the specific conditions under which $R = 0$ becomes an unstable equilibrium point in the probability.

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 gratitude to the anonymous reviewers for their thoughtful and constructive feedback, which significantly contributed to improving the quality of this manuscript. This research was partially supported by the Project 2023-58830 from the University of Antioquia.

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix

Stability

Definition .1. Equilibrium point of an ordinary differential equation.

Consider the ordinary differential equation:

where $b$ is a continuous function on a non-empty interval $I \subset \mathbb{R}$, and let $\xi \in I$. We say that $\xi$ is an equilibrium point of (A.1), or that $\xi$ is a trivial solution of (A.1), if and only if $b(\xi) = 0$.

Definition .2. Stability of an equilibrium point of an ordinary differential equation.

Let $\xi$ be an equilibrium point of (A.1). We say that:

1) $\xi$ is a stable equilibrium point of (A.1) if and only if for every $\epsilon > 0$, there exists $\delta > 0$ such that for all $x_0 \in I$, with $|x_0 - \xi| < \delta$, we have:

Otherwise, we say that $\xi$ is an unstable equilibrium point of (A.1).

2) $\xi$ is an asymptotically stable equilibrium point of (A.1) if and only if $\xi$ is a stable equilibrium point of (A.1) and there exists $\delta > 0$ such that for all $x_0 \in I$ with $|x_0 - \xi| < \delta$, we have:

3) $\xi$ is a globally asymptotically stable equilibrium point of (A.1) if and only if $\xi$ is a stable equilibrium point of (A.1) and for all $x_0 \in I$, we have:

Here, $x(t; x_0)$ represents the global solution $x$ of (A.1) that starts at $x_0$ at time $t$.

Definition .3. Equilibrium point of a stochastic differential equation.

Consider the stochastic differential equation:

where $b$ and $\sigma$ are continuous functions on an interval $I \subset \mathbb{R}$, $B$ is a Brownian motion on a probability space $(\Omega, \mathcal{F}, \mathbb{P})$, and let $\xi \in I$. We say that $\xi$ is an equilibrium point of (A.2), or that $\xi$ is a trivial solution of (A.2), if and only if $b(\xi) = \sigma(\xi) = 0$.

Definition .4. Stability of an equilibrium point of a stochastic differential equation.

Let $\xi$ be an equilibrium point of (A.2). We say that:

1) $\xi$ is a stable equilibrium point in probability of (A.2) if and only if for every $\epsilon > 0$ and $r \in (0, 1)$ there exists a $\delta > 0$ such that for all $x_0 \in I$, $|x_0 - \xi| < \delta$, it holds that

In the contrary case, we say that $\xi$ is an unstable equilibrium point of (A.2).

2) $\xi$ is an asymptotically stable equilibrium point in probability of (A.2) if and only if $\xi$ is a stable equilibrium point in probability of (A.2) and for every $r \in (0, 1)$ there exists a $\delta > 0$ such that for all $x_0 \in I$, $|x_0 - \xi| < \delta$, it holds that

3) $\xi$ is a globally asymptotically stable equilibrium point in probability of (A.2) if and only if $\xi$ is a stable equilibrium point in probability of (A.2) and for all $x_0 \in I$, it holds that

Where $X(\cdot, t; x_0)$ represents the global solution $X$ of (A.2) that starts at $x_0$ at time $t.$

Definition .5. Equilibrium point of a fractional differential equation.

Consider the fractional differential equation:

where $b$ is a continuous function on a non-empty interval $I$ of $\mathbb{R},$ $\alpha \in (0, 1),$ and let $\xi \in I.$ We say that $\xi$ is an equilibrium point of (A.3) or that $\xi$ is a trivial solution of (A.3) if and only if $b(\xi) = 0.$

Definition .6. Stability of an equilibrium point of a fractional differential equation.

Let $\xi$ be an equilibrium point of (A.3). We say that:

1) $\xi$ is a stable equilibrium point of (A.3) if and only if for every $\epsilon > 0$, there exists a $\delta > 0$ such that for every $x_0 \in I$, $|x_0 - \xi| < \delta$, it holds that

Otherwise, we say that $\xi$ is an unstable equilibrium point of (A.3).

2) $\xi$ is an asymptotically stable equilibrium point of (A.3) if and only if $\xi$ is a stable equilibrium point of (A.3) and there exists $\delta > 0$ such that for every $x_0 \in I$, $|x_0 - \xi| < \delta$, it holds that

3) $\xi$ is a globally asymptotically stable equilibrium point of (A.3) if and only if $\xi$ is a stable equilibrium point of (A.3) and for every $x_0 \in I$, it holds that

where $x(t; x_0, \alpha)$ represents the global solution $x$ of (A.3) that starts at $x_0$ at time $t.$