On the contribution of qualitative analysis in mathematical modeling of plasmid-mediated ceftiofur resistance

1.   Introduction

In the area of mathematical modeling, it is a very common practice to use continuous dynamical systems defined through systems of ordinary differential equations to model population dynamics in different branches of knowledge. The purpose of this type of model is to predict the temporal evolution of interacting populations and characterize laws, properties or intrinsic patterns that govern their behavior ^[1].

There are multiple examples of mathematical models that have contributed to the solution of problems in real contexts and have generated theories that are still valid and applied today, such as the case of Hammer and his law of mass action, which expresses that the rate at which a disease spreads is proportional to the number of individuals susceptible to contracting it, multiplied by the number of infectious ^[2]. Ross developed a mathematical theory that correctly explained that to eradicate malaria it was not necessary to eliminate the total population of mosquitoes, but it was enough to maintain their population level below a threshold ^[3]. Kermack and McKendrick established the threshold theorem, which postulates that the introduction of an infectious individual in a community will not give result in an outbreak unless the density of the susceptible population exceeds a certain critical value ^[4]. Lotka and Volterra with their prey-predator model, described the interaction in which one organism, the predator, eats all or part of the body of another organism, the prey ^[5,6]. More recently, Stewart and Levin created the Stewart-Levin criterion, which determines the equilibrium frequencies of plasmid-carrying cells in terms of key modeling parameters: Population growth, conjugational transfer and segregation rate ^[7]. These are just some of the most cited cases in the literature, and their success is based on the complement of different factors such as:

$1)$ The robustness of the mathematical models; that is, well-defined models whose variables and parameters are consistent with the dynamics to be studied. In addition, the results of the qualitative analysis of the model are coherent with the phenomenon and determine the behavior in terms of the model parameters.

$2)$ The validation of the model through data obtained from experimentation in vivo, in vitro, or in field work.

$3)$ Parameter estimation

$4)$ Sensitivity analysis of parameters

$5)$ Numerical simulations

To carry out the first item, there is a broad theoretical analysis of autonomous and non-autonomous dynamical systems ^[8,9]. Items two through five require statistical, inferential, probabilistic, numerical and computational methods to be carried out ^[10]. As we can see, to carry out an investigation that contemplates the five previous points, interdisciplinary work is required. This implies that in areas such as microbiology this becomes a very difficult task because experimental data is not always available or the level of complexity of the models is very high. At present, in silico experimentation or artificial intelligence are areas that are supporting mathematical modeling. However, they are emerging sciences that are just being coupled and require an adaptation time to contribute in the best way to modeling ^[11,12].

At present, antimicrobial resistance is one of the main public health problems at the global level, and its sudden increase in the last decades has affected the health systems of both low-income and high-income countries ^[13]. In particular, it is of great relevance to advance in the understanding of the role played by plasmids in the acquisition of antibiotic resistance mediated by plasmids. In this regard, there are works focused on characterizing properties that determine the outcome of infection in the interaction dynamics between sensitive and resistant strains ^{[1,14,15,16,17,18,19]}, on pharmacokinetic (PK)/pharmacodynamic (PD) models that study the effects of the antibiotic ^{[20,21,22,23,24]}, among other topics.

Volkova et al. ^[24] established that antimicrobial use in food for animals may contribute to antimicrobial resistance in bacteria of animals and humans. Commensal bacteria of an animal intestine may serve as a reservoir of resistance-genes. In this sense, they developed a mathematical model on plasmid-mediated ceftiofur resistance in commensal enteric Escherichia coli (E. coli) of cattle, in which they implemented the last four literals mentioned above, except the qualitative analysis of the model in item one. In this work, we focus on analyzing the robustness of the model, and determining laws, patterns and properties that arise from the theoretical results, as well as relating the qualitative analysis with the results of Volkova et al. ^[24].

2.   Materials and methods

In order to study plasmid-mediated ceftiofur resistance in commensal enteric E. coli of cattle, Volkova et al. formulated systems of ordinary differential equations that measure both the therapeutic effect of parenteral treatment with ceftiofur on competition dynamics between ceftiofur-sensitive and resistant commensal enteric E. coli, as well as the impact of the absence of treatment on competition dynamics. In a way, the authors are analyzing the pharmaceutical response of ceftiofur in the outcome of the infection, which is closely related to the role played for both the PK and the PD of the antibiotic in the propagation of bacterial population within the host. Since the models formulated by Volkova et al. ^[24] incorporate both the PK and PD of ceftiofur, in this section we will present the concepts, as well as the models and data used by them.

2.1.   Pharmacokinetic modeling

The PK describes the behavior of an administered drug in the body over time and is currently defined as the study of the absorption, distribution, metabolism and excretion of a drug ^[8]. We could say that the quantitative description of kinetic patterns of concentration due to the four drug properties mentioned above frame PKs. In general, these patterns are measured experimentally by means of kinetic parameters that allow determining the evolution of the drug concentration with respect to time. PK models are phenomenological or empirical approaches to describe drug concentration that allows us to compare their results with experimental evidence ^[25]. In ^[24], based on a literature review and experimental data, Volkova et al. developed an experimental protocol on the PKs of ceftiofur metabolities, which led them to the formulation of PK models for two different patterns. They began by defining concentration of ceftiofur equivalents (CE) as the total of ceftiofur and its active metabolites. Due to there being no published data on the pattern or inter-individual variability of the biliary excretion, they assumed for the protocol there was no enterohepatic ceftiofur circulation, and CE-concentration in the intestine was independent of that in systematic distribution. Let $C(t)$ be the CE-concentration per $g$ of feces in the large intestine at time $t$. From review, Volkova et al. concluded that $C(t)$ seemed to decay exponentially, and since the pattern of biliary excretion of ceftiofur is unknown, they explored two possibilities in a therapy with repeated injections of a non-sustained-release formulation of ceftiofur ^[24]. The parameters used in the pharmacokinetics of ceftiofur metabolites were

● Ceftiofur dose in one injection $D$.

● Fraction of ceftiofur dose in one injection excreted in bili $p$.

● Volume of the animal's large intestine $V$.

● Biodegradation at rate $\lambda$.

● Excretion fraction number $k$.

In the protocol, Volkova et al. assumed that, $D$, $p$ and $V$ were taken to be equal for every injection $n_j$, which occurred at time $T_j$ since start of treatment, $t = 0$. Based on the above information, they proposed the following concentration patterns.

Pattern 1: Amount $Dp$ was excreted at one hour post injection (p.i). After a passage time $T_{\delta}$, biliary metabolites entered the large intestine. At entry, for a given $n_j$, $C(t)$ per $g$ of feces (assuming weight-to-volume ratio of feces of one) was $Dp/V$, then it decayed exponentially due to the biodegradation at rate $\lambda$. In consequence, $C(t)$ is given by

where

being $n = 5$, $T_j = 24\left(n_j-1\right)$, and the injection number $j = 1, \dots, n$.

Note that the function $C(t)$ defined in (2.1) is a piecewise continuous function with a discontinuity at $t^* = T_j+1+T_{\delta}$. In ^[24], Volkova et al. modeled seven ceftiofur treatment scenarios for cattle and compared their results with experimental data. In this work we will use the data from the scenario $R_1$ which corresponds to a six-month-old dairy cattle weighing 180 kg and treated for the bovine respiratory disease named interdigital necrobacillosis. Table 1 shows the biliary ceftiofur metabolites parameters.

Using data from Table 1 we obtain the following expression for CE-concentration of Pattern 1.

Pattern 2: Six equal fractions of $Dp$ were excreted hourly at hour one to six p.i. The choice of hours one to six p.i. was based on the working hypothesis that $T_{\delta} = 6$ hours, thus the entire amount $Dp$ would reach the large intestine by 12 hours p.i and the initial concentration of $\mu_g/g$ feces was $Dp/6V$. The concentration C (CE $\mu_g/g$) at time $t$ was as in (2.1), $c_j$ for a given $n_j$ was

being $m = 6$, $k = 1, \dots, m$ and

being $n = 5$, $T_j = 24\left(n_j-1\right)$, and the injection number $j = 1, \dots, n$. In this case, the expression for CE-concentration of Pattern 2 is given by

where $c_{jk}(t)$ is defined in (2.3) and estimated using data from Table 1.

2.2.   Pharmacodynamics modeling

PD is the study of the molecular, biochemical, and physiologic effects or actions of drugs. The effect of the drug is associated with receptors, which are the most important target of therapeutic drugs. The main functions of the drug effect are the binding of the drug to the receptor, drug-induced receptor activation, and the propagation of this initial receptor activation into the observed drug effect. The drug-receptor interaction was developed ^[25] which deduced the most common function used to relate drug concentration with the pharmacological effect, the $E_{\max}$ model defined by

where $E_{\max}$ is the maximum effect, $Ec_{50}$ is the concentration at half the maximal observable in vivo effect, $\rho$ are drug molecules and the drug concentration $C$ is an equilibrium point of the differential equation that model drug–receptor complex interaction. In this case, (2.4) defines a static nonlinear model in which $C$ is constant. If we consider $C$ as a time course $C(t)$, we must implicitly assume that equilibrium is achieved rapidly throughout $C(t)\equiv C$ (see ^[25], chapter 12). When a baseline $E_0$ is introduced to (2.4), we obtain the $E_{\max}$ model describing either stimulation or inhibition of the effect by the concentration of the drug given by

In ^[24], Volkova et al. presented the protocol for the PK effect in which the changes in the net growth of ceftiofur-sensitive and resistant enteric E. coli depending on CE-concentration were modeled by means of the inhibitory version of $E_{\max}$ model defined in (2.5). Specifically, they defined $E_0 = 1$, $E_{\max} = 2$, $\rho = H$ as the Hill coefficient and $Ec_{50} = MIC$ as the minimum inhibitory concentration. Therefore, the fractional changes in net growth of ceftiofur-sensitive and resistant E. coli, respectively, at CE-concentration C are given by

where $MIC_s$ is the $MIC$ for ceftiofur-sensitive and $MIC_r$ is the $MIC$ for ceftiofur-resistant. In the protocol, they consider the following possibilities

PD-effect A Multiplicative PD effect on E. coli net growth with a constant minimum inhibitory concentration (MIC = 1 $\mu$g/mL) and changing ceftiofur concentration expressed as multiples of MIC; Hill coefficient = 1.5.

PD-effect B Multiplicative PD effect on E. coli net growth with changing MIC and a constant ceftiofur concentration (4 $\mu$g/mL); Hill coefficient = 1.5.

Under the assumption that the CE-concentration defined in (2.2) reaches an equilibrium $C$ that satisfies the conditions of the dynamics exposed in PD-effect A, we have that $E_s$ and $E_r$ defined in (2.6) are rewritten as

where $\phi\in(0, \infty)$. Observe that $E_s(0) = 1 = E_r(0)$, which indicates that for $\phi = 0$ the dynamics in absence and during treatment coincide. On the other hand,

Similarily, for PD-effect B we obtain

From (2.9) we obtain the following properties

From (2.8) and (2.10) we observe that the PD- effect A defined in (2.7) is completely opposite to PD-effect B defined in (2.9). Those $PD-effects$ will be substituted in (2.11) to perform numerical simulations of the bacterial dynamics during parenteral ceftiofur treatment.

2.3.   Volkova-Lanzas-Lu-Gröhn model

In ^[24], V. Volkova et al. formulated a system of ordinary differential equations to explore the dynamics in vivo of ceftiofur-resistant and sensitive commensal enteric E. coli in the absence and during parenteral ceftiofur therapy within the gut of cattle. The model is given by

where $N_s(t)$ and $N_r(t)$ denote the number of sensitive and resistant E. coli cells at time $t$, respectively, $N(t) = N_s(t)+N_r(t)$ is the total number of E. coli cells a time $t$, $N_{\max}$ is the upper limit for total E. coli per g of feces, $r$ is the maximum net growth rate (in exponential growth phase) in numbers of enteric E. coli and $\beta$ is frequency-dependent transmission for blaCMY-2-carrying plasmids from resistant donor to sensitive cells; that is, it is the plasmmid transfer rate, $v$ represents the fraction of the ingested bacteria carrying plasmids with blaCMY-2, $\gamma$ is both the rates of hourly fractional in-flow and outflow of enteric E. coli, $\alpha$ is the reduction rate of the net growth rate due to the fitness cost, and $E_s$ and $E_r$ are defined in (2.6), and measure the PD effect on the net growth rate, $r$, for sensitive and resistant E. coli, respectively. Note that $E_s = E_r\equiv1$ corresponds to the bacterial dynamics without treatment, while $E_s\neq1$ or $E_r\neq1$ corresponds to the dynamics during treatment.

Remark 2.1. It is important to highlight that in the context of the biological phenomenon and the definition of the parameters these are within the following ranges: $r > 0$, $0 < \beta\leq1$, $N_{\max} > 1$, $0\leq v\leq1$, $0 < \alpha < 1$ and $0\leq\gamma\leq1$.

3.   Qualitative analysis of the model without treatment

By substituting $E_s = E_r = 1$ in the system of differential equation (2.11) we obtain the dynamical system of ceftiofur-sensitive and resistant commensal enteric E. coli in the absence of immediate ceftiofur pressure. In this section we will do the qualitative analysis of the above system.

3.1.   Invariant set

The matrix form of (2.11) is given by

where

and

Let $\mathbb{R}_+ = \{p\in\mathbb{R}:p\geq0\}$. Since $f\in\mathcal{C}^1(\mathbb{R}_+^2\setminus(0, 0))$, there exists an $\varepsilon > 0$ such that the initial value problem (ivp) defined by (3.1) and $x(0) = x^0\in\mathbb{R}^2$ has a unique solution $x(t)$ on the interval $[0, \varepsilon]$ (fundamental theorem of existence and uniqueness ^[9]). In addition, for a compact set $\Omega\subset \mathbb{R}_+^2$ such that

it follows that $\varepsilon = \infty$. In consequence, the ivp has a solution $x(t)$ for all $t\geq0$ (Corollary 2, p91 ^[9]). Now, by adding the two equations of (2.11) we obtain

where $N = N_s+N_r$. In consequence, from (3.4) we obtain

The solution $N(t) = N_s(t)+N_r(t)$ of inequality (3.5) is globally bounded. In effect, let us suppose that there exists a time $t$ with $N(t) = N_s(t)+N_r(t) = \xi > N_{\max}$. Let $t^*$ be the infimum of those times where this inequality holds for a fixed choice of $\xi$, then $dN(t^*)/dt = dN_s(t^*)/dt+dN_r(t^*)/dt\geq0$. On the other hand, from (3.5) it is followed that

Since $N(t^*)/N_{\max} > 1$, then from (3.6) we obtain $dN(t^*)/dt < 0$, which is a contradiction. On the other hand, since $N'(0) = 0$ and $N''(0) = r > 0$ then $N'(t)$ is nonnegative increasing function at $t = 0$, and since $N(0)\geq 0$, then $N(t)$ is nonnegative increasing function at $t = 0$. In consequence, the solution of inequality (3.5) satisfies $0\leq N(t) = N_s(t)+N_r(t)\leq N_{\max}$ for all $t\geq0$. Therefore, the invariant set of (2.11) is given by the compact set

3.2.   Equilibrium points

The equilibria of (2.11) are given by the solutions of the algebraic equation

Now, we will find the equilibrium points for which $N_s = 0$. In this case, (3.8) is reduced to

Observe that $N_r = 0$ is a solution of (3.9), which implies the existence of the equilibrium solution $x_0 = (0, 0)^T$ of (2.11). However, the function $z$ defined in (3.2) is not defined in the point $x_0$. Therefore, $x_0$ is not an equilibrium point of (2.11). When $N_r > 0$, we obtain $(1-v)\gamma = 0$ from the first equation of (3.9), and $N_r = N_{\max}$ from the second ones. In consequence, if $(1-v)\gamma = 0$ there exists an equilibrium point $x_1 = (0, N_{\max})^T$. Now, when $N_r = 0$, (3.8) is reduced to

The solutions of (3.10) are $N_s = 0$, and if $v\gamma = 0$ then $N_s = N_{\max}$. Therefore, if $v\gamma = 0$ there exists an equilibrium point $x_2 = (N_{\max}, 0)^T$. Now, we will find the equilibrium points for which $N_s > 0$ and $N_r > 0$. By adding the two equations of (3.8) we obtain

Since $\alpha < 1$, then the solution of (3.11) is given by $N = N_s+N_r = N_{\max}$. From the above equation we obtain

By substituting $N = N_{\max}$ in (3.8) we obtain the following system

By replacing (3.12) in the second equation of (3.13) we obtain the following quadratic equation

The only positive solution of (3.14) is

By substituting (3.15) in (3.12) we obtain

The existence results are summarized in the following proposition.

Proposition 3.1. The systen (2.11) always has the equilibrium $x_+ = (N_s^+, N_r^+)^T$. If $(1-v)\gamma = 0$ there exists the resistant bacteria-equilibrium $x_1 = (0, N_{\max})^T$, and if $v\gamma = 0$ there exists the sensitive bacteria-equilibrium $x_2 = (N_{\max}, 0)^T$.

3.2.1.   Existence analysis of equilibrium $x_0$

As we could observe, the point $x_0$ is not an equilibrium solution of (2.11). However, if we define

Then $x_0 = (0, 0)^T$ is an equilibrium solution of the system

Since

then the function $\tilde{f}$ is continuous in $x_0$. Therefore, $\tilde{f}\in\mathcal{C}(\mathbb{R}_+^2)$. Now, the natural derivative of $\bar f$ could be

However, $D\tilde{f}$ is not a continuous function in $x_0$. Therefore $\tilde{f}$ is differentiable but not continuously differentiable, which prevents us from using the Picard's Theorem to guarantee the uniqueness of the equilibrium solution $x_0$ of (3.18).

Remark 3.2. Note that in the absence of the ceftiofur pressure ($E_s = E_r = 1$), all the equilibrium solutions of (2.11), including the infection-free point $x_0$ are on the boundary of the set $\Omega$.

3.3.   Stability analysis

In this section we determine the stability of the equilibrium points of (2.11). The linearization of (2.11) around an equilibrium point $\bar x$ is given by $x' = J(\bar x)x$, where the Jacobian matrix $J$ evaluated at $x$ is given by

where

3.3.1.   Stability of the resistant-bacteria equilibrium $x_1$

In the case $(1-v)\gamma = 0$, the Jacobian matrix evaluated at $x_1 = (0, N_{\max})^T$ is given by

The eigenvalues of $J(x_1)$ are $\lambda_1 = -(v\gamma +\beta) < 0$ and $\lambda_2 = -r(1-\alpha) < 0$, which implies that $x_1$ is locally asymptotically stable (l.a.s) in $\Omega$. This result is summarized in the following proposition.

Proposition 3.3. The resistant-bacteria equilibrium $x_1$ is l.a.s. in $\Omega$.

3.3.2.   Stability of the sensitive-bacteria equilibrium $x_2$

In the case $v\gamma = 0$, the Jacobian matrix evaluated at $x_2 = (N_{\max}, 0)^T$ is given by

The eigenvalues of $J(x_2)$ are $\bar\lambda_1 = -r$ and $\bar\lambda_2 = \gamma\left(R_0-1\right)$, where

Since $\lambda_1 < 0$ and $\lambda_2 < 0$ if and only if $R_0 < 1$, we have the following result.

Proposition 3.4. If $R_0 < 1$, then sensitive-bacteria equilibrium $x_2$ is l.a.s. in $\Omega$.

Remark 3.5. In the case $\gamma > 0$, it implies that Proposition 3.4 makes sense when $v = 0$.

Remark 3.6. Since $\beta$ is the plasmid transfer rate and $\gamma$ is the outflow rate, then $R_0 = \beta/\gamma$ can be interpreted as the offspring ratio of bacteria measuring the new bacteria generated by resistant bacteria population.

3.3.3.   Stability of the coexistence-bacteria equilibrium $x_+$

The Jacobian matrix evaluated at $x_+ = (N^+_{s}, N_r^+)^T$ is given by

The characteristic polynomial of $J(x_+)$ is

where

Observe that if $N^+_r\geq N^+_s$, then $a_1 > 0$ and $a_2 > 0$. Now, we will determine the conditions for which the above inequality is satisfied. Observe that

If

then $N^+_r- N^+_s\geq0$, which implies that $a_1 > 0$ and $a_2 > 0$. Therefore, from the Routh-Hurwitz criterion we conclude that the eigenvalues of $J(x_+)$ have negative real part. In consequence, we have the following result.

Proposition 3.7. If the condition (3.28) or equivalently $R_0 > 2(1-2v)$ is satisfied, then the equilibrium point $x_+$ is l.a.s. in $\Omega$.

Remark 3.8. Note that (3.28) establishes a relation between the fraction of the ingested bacteria carrying plasmids with blaCMY-2, $v$, and the number of new bacteria generated by them, $R_0$, suggests that the larger the fraction of ingested bacteria, the lower the number of bacteria generated by a bacterium carrying plasmids with resistance genes.

3.3.4.   Stability analysis of the free-bacteria equilibrium $x_0$

The linearization of (3.18) around of an equilibrium point $\bar x$ is given by $x' = \tilde{J}(\bar x)x$, where the Jacobian matrix $\tilde{J}$ evaluated at $x$ is given by

From (3.29) we obtain $\tilde{J}(x_0) = A$, and its characteristic polynomial is given by

where

If $\tilde{a}_1 > 0$ and $\tilde{a}_2 > 0$, then $x_0$ is locally asymptotically stable. Above implies the following result.

Proposition 3.9. If

then $x_0$ is locally asymptotically stable in the closure of $\Omega$ given by

Remark 3.10. The product $\gamma v\alpha$ can be interpreted as the fraction of ingested bacteria carrying blaCMY-2 plasmids that acquired resistance to ceftiofur and were subsequently secreted.

4.   Qualitative analysis of the model during treatment

The system of differential equations (2.11) describes the dynamical system of ceftiofur-sensitive and resistant commensal enteric E. coli during parenteral ceftiofur treatment. In this section we will do the qualitative analysis of the above system. We begin by highlighting that in the PD modeling carried out by Volkova et al. ^[24], the PD effect defined in (2.6) does not vary with respect to time $t$. This is due to the fact that both $MIC$ and the stationary ceftiofur concentration $C$ were taken constant with respect to time $t$. However, in scenario A, $MIC$ is constant and $C$ is a multiple of $MIC$. In this sense, the PD effect varies with respect to a parameter but not with respect to time. Scenario B is interpreted in a similar way. Consequently, in this case the nonlinear system of ordinary differential equations (2.11) is autonomous. For this reason, we will use the theory of autonomous systems to determine their equilibrium solutions.

Carrying out the same procedure as in the subsection 3.1 (invariant set), it is verified that the set $\Omega$ defined in (3.7) is the invariant set of system (3.4).

4.1.   Equilibrium points

In a similar way to the developed $N_s = 0$ o $N_r = 0$ in subsection 3.2, the following result is verified.

Proposition 4.1. If $(1-v)\gamma = 0$ there exists the resistant bacteria-equilibrium $\bar x_1 = (0, N_{\max})^T$, and if $v\gamma = 0$ there exists the sensitive bacteria-equilibrium $\bar x_2 = (N_{\max}, 0)^T$.

Now, for $N_s > 0$ and $N_r > 0$, by adding the two equations of (3.8) we obtain

Since $-1\leq E_s\leq1$ and $-1\leq E_r\leq1$, then the solutions of (4.1) are

From the first one solution $\tilde{N}_r$ of (4.2) we obtain the coexistence equilibrium $\bar x_+ = x_+ = (N_s^+, N_r^+)^T$, which is the same equilibrium defined in Proposition 3.1. Before determining the second coexistence equilibrium, let us note that the second solution $\bar{N}_r$ of (4.2) has biological sense if and only if $E_s$ and $E_r$ have opposite signs. Now, by substituting $\bar{N}_r$ into the total bacterial population $N$ we obtain $\bar{N} = a\bar{N}_s$ where $a = 1+E_s/(\alpha-1)E_r$. Taking the righthand side of the first equation of (2.11) equal to zero and substituting $\bar{N}_r$, $\bar{N}_s$ and $\bar{N}$ in the resulting equation we obtain the following equation

The solution of (4.3) is

From (4.4) we conclude that $\bar{N}_s > 0$ if and only if

By substituting, $\bar{N}_s$ in the second equation of (4.2) we obtain

The following proposition summarized the existence result for coexistence equilibria.

Proposition 4.2. The system (2.11) always has the equilibrium $\bar x_+ = (N_s^+, N_r^+)^T$. If $E_s$ and $E_r$ have opposite signs and the inequality (4.5) is satisfied then there exists the coexistence equilibrium $\bar{x} = (\bar{N}_s, \bar{N}_r)^T$.

4.2.   Stability analysis

In this subsection we will determine the stability of the equilibrium points of (2.11) during the treatmet. The linearization of (2.11) around an equilibrium point $\hat{x}$ is given by $x' = \bar{J}(\hat x)x$, where the Jacobian matrix $\bar J$ evaluated at $x$ is given by

where

and

Following a procedure similar to the one carried out in subsection 3.3.1, we verify that the eigenvalues of the Jacobian $\bar J(\bar x_1)$ are given by $\bar \lambda_1 = -(v\gamma+\beta)$ and $\bar \lambda_2 = -r(1-\alpha)E_r$. In consequence, if ($\alpha > 1\, \wedge\, E_r < 0$) or ($\alpha < 1\, \wedge\, E_r > 0$), then $\bar x_1$ is locally asymptotically stable in $\Omega$. Similarly, we verify that the eigenvalues of the Jacobian $\bar J(\bar x_2)$ are given by $\bar \lambda_3 = -rE_s$ and $\bar \lambda_4 = \beta-\gamma$. Therefore, if $E_s > 0$ and $\beta < \gamma$, then $\bar x_2$ is locally asymptotically stable in $\Omega$. We can also verify that the characteristic polynomial of $\bar{J}(\bar{x}_+)$ is given by

where

Observe that if $E_s > 0$, $E_r > 0$ and the inequality (3.28) is satisfied, then $g_1 > 0$ and $g_2 > 0$, which implies the local stability of $\bar{x}_+$. Since $\alpha < 1$, then the above results are summarized in the following proposition

Proposition 4.3. If $E_r > 0$, then $\bar x_1$ is l.a.s. in $\Omega$. If $E_s > 0$ and $R_0 < 1$, then $\bar x_2$ is l.a.s. in $\Omega$, and if $E_s > 0$, $E_r > 0$ and the inequality (3.28) is satisfied, then $\bar{x}_+$ is $\textrm{l.a.s}$ in $\Omega$.

By substituting $b = E_s/(\alpha-1)E_r$ in (4.2) we obtain $\bar x = (\bar N_s, b\bar N_s)^T$ where $\bar N_s$ is defined in (4.4). From direct calculations we verify that the Jacobian $\bar J$ defined in (4.6) evaluated at $\bar x$ is given by

where

After some simplifications we verify that the trace and the determinant of $\bar J(\bar x)$ are given by

As we can see, $trace\left(\bar J(\bar x)\right)$ and $det\left(\bar J(\bar x)\right)$ do not have defined sign. In fact, to determine ranges of the parameters for which they have a defined sign is a difficult task. However, in the section on numerical simulations we will determine values of the parameters for which $trace\left(\bar J(\bar x)\right) < 0$ and $det\left(\bar J(\bar x)\right) > 0$, which implies the stability of the equilibrium $\bar x$. Finally, in the following proposition we verified that (2.11) does not have periodic solutions in $\Omega$.

Proposition 4.4. System (2.11) does not have periodic orbits in $\Omega$.

Proof. To prove this result, we will use the Theorem 2 of ^[26], which established that if there exists a function $c\in\mathcal{F}_{\Omega}$ where

such that $h$ is a solution of the system

being $F = (f_1, f_2)^T$ the righthand side of (2.11), with $h\in\mathcal{F}_{\Omega}$, then $h$ is a Dulac function for (2.11) on $\Omega$. For

is satisfied by (4.7) it implies that $h$ is a Dulac function of (2.11). Therefore, (2.11) does not have periodic orbits in $\Omega$. □

5.   Numerical simulation

In this section we use data of one of the scenarios modeled by Volkova et al. in ^[24] (six month beef) to perform numerical simulations of (2.11). The data in Table 2 was used for the simulations of the dynamics without treatment. In this case, the stability condition (3.28) is satisfied, which implies that solutions tend to the equilibrium $x_+ = (310535, 5692)$ on the boundary of the $\Omega$ (See Figure 1). From Figure 1 we see that if the ingested bacteria contains both sensitive and resistant strains, in the absence of the antibiotic the population of resistant bacteria will persist.

Before presenting simulations for (2.11) during the treatment, let us note that in this case the model is a nonlinear system of non-autonomous differential equations. In fact, since the concentration of ceftiofur at time $t$ defined in (3.4) and $C(t)$ is a piecewise continuous function, then the functions of the PD effects with respect to time $t$, $E_s (t)$ and $E_r(t)$ are also discontinuous functions, which implies that (2.11) is a nonsmooth dynamical system. However, due to properties of pharmacotherapy, it is assumed that equilibrium $C$ is rapidly reached by $C(t)$, $C(t)\equiv C$ (see ^[25], chapter 12). In this sense, both $E_s(t)$ and $E_r(t)$ will reach the equilibria $E_s$ and $E_r$ defined in (2.6). Figure 2 shows numerical simulations of both ceftiofur-sensitive and resistant E. coli under the PD-effect A defined in (2.7). Since both $E_s$ and $E_r$ depend on the ratio $\phi$ between CE-concentration and MIC, we perform the simulations with different values of $\phi$. the lefthand side of Figure 2 shows graphs of the time course of E. coli population sensitive to ceftiofur and the righthand side shows graphs of the time course of E. coli population resistant to ceftiofur. For $\phi = 0$ we observe that the population of resistant bacteria presents a small increase, remaining almost constant, while the population of sensitive bacteria grows in such a way that the sum of both bacterial quantities reaches the carrying capacity, $N_{ \max}$. This dynamic coincides with the dynamics in the absence of the antibiotic. For $\phi\in(0, 2]$ the bacterial population reaches an equilibrium of coexistence in the interior set of $\Omega$ with the characteristic that as $\phi$ grows, the population of sensitive bacteria decreases and increases the population of resistant bacteria. For $\phi\in(2, 4]$ both populations tend to decrease as $\phi$ approaches the value $\phi = 4$ where both populations cancel out. Figures 3 and 4 show numerical simulations of both ceftiofur-sensitive and resistant E. coli under the PD-effect B defined in (2.9). These simulations were performed for nine equidistant values of $MIC_s$ in the interval $[0, 4]$, and the fixed value $MIC_r = 10$ for Figure 3 and $MIC_r = 0.001$ for Figure 4. the lefthand side of Figure 3 shows graphs of the time course of E. coli population sensitive to ceftiofur and the righthand side shows graphs of the time course of E. coli population resistant to ceftiofur. For $MIC_s = 4$ we observe that the population of resistant bacteria presents a slight increase, and then decreases, reaching an equilibrium close to its initial population, while the sensitive bacteria maintains a sustained increase until reaching equilibrium in such a way that the sum of both bacterial populations reaches the capacity bacterial load, $N_{ \max}$. In this dynamic it is observed that as $MIC_s$ increases, the population of sensitive bacteria decreases, while the population of resistant bacteria increases. In Figure 4 both bacterial populations tend to be cleared. We have also verified numerically that when $MIC_r$ decreases the behavior shown in Figure 4 becomes a trend.

6.   Discussion

Qualitative analysis of the model without treatment corresponds to the scenario of ceftiofur-sensitive and resistant commensal enteric E. coli in the absence of immediate ceftiofur pressure. In this case, we verify that there is always an equilibrium of coexistence of both bacterial populations on the boundary of $\Omega$, $x_+ = (N_s^+, N_r^+)^T$, whereas if there is no inflow and no outflow of E. coli ($\gamma = 0$), no ingestion of bacteria carrying plasmids with blaCMY-2 ($v = 0$), or if $100%$ of ingested bacteria carries blaCMY-2 plasmids, then there are two equilibrium on the boundary of $\Omega$, the ceftiofur-resistant equilibrium $x_1 = (0, N_{\max})^T$ and the ceftiofur-sensitive equilibrium $x_2 = (N_{\max}, 0)^T$. We have also proved that if ($v = 1$ and $\alpha < 1$) or ($\gamma = 0$ and $\alpha < 1$) then $x_1$ is l.a.s. in $\Omega$. Under the assumption that the reduction rate of the net growth rate due to the fitness cost is less than one, the previous results suggest that if $100%$ of ingested bacteria are ceftiofur-resistant or if there is no inflow and no outflow of commensal E. coli, then the population of ceftiofur-sensitive commensal enteric E. coli will be cleared while ceftiofur-resistant commensal enteric E. coli will reach carrying capacity $N_{\max}$. In other words, the infection will persist, but only with resistant bacteria. On the other hand, if $v = 0$ and $\beta < \gamma$ then $x_2$ is l.a.s. in $\Omega$, which suggests that if $100%$ of ingested bacteria are ceftiofur-sensitive and the plasmid transfer rate is less than in-flow and outflow rates then the population of ceftiofur-resistant commensal enteric E. coli will be cleared while ceftiofur-sensitive commensal enteric E. coli will reach carrying capacity $N_{\max}$. Now when, the condition (3.28) is satisfied, then $x_+$ is l.a.s. in $\Omega$ which means that ceftiofur-sensitive and resistant commensal enteric E. coli will proliferate in a controlled way until they reach a steady state. From (3.27) it is observed that the size of the bacterial populations and the speed at which they reach equilibrium will depend on the parameters $\beta$, $v$ and $\gamma$. However, in (3.28) a sufficient condition for stability of $x_+$ is established that only depends on $\beta$ and $v$.

We have verified that for an extension of the dynamical system (3.1) there exists an infection-free equilibrium $x_0$ on the boundary of $\Omega$, which is l.a.s when (3.32), depending on $\alpha, \, \gamma, \, r$ and $v$ is satisfied. The qualitative analysis of (2.11) in the absence of treatment correlates the results of Volkova et al. ^[24] in which the authors verified that reported fractions of enteric commensal E. coli carrying blaCMY-2 in cattle could persist in the absence of immediate ceftiofur pressure. Even more, these fractions were more sensitive to $\beta$ and $v$. In fact, if $v$ held constant, then they would be more sensitive to $\beta$ and $\gamma$. By taking $\beta\geq0.01$ (somewhat unrealistic for authors) allowed them to reproduce the reported fraction if there was no resistance to ceftiofur among ingested E. coli, $v = 0$. In addition, through qualitative analysis, specific conditions were determined for which bacterial populations are eliminated, or persist only with sensitve bacteria, resitanst bacter or with both populations. It is also noted that all equilibrium solutions are in the boundary set of $\Omega$, which implies that the bacterial populations reach equilibrium in the carrying capacity $N_{\max}$.

Qualitative analysis of the model during treatment corresponds to the scenario of ceftiofur-sensitive and resistant commensal enteric E. coli during parenteral ceftiofur treatment. In this case, the existence of the same equilibrium points on the boundary of $\Omega$, ($x_0, \, \bar{x}_1 = x_1, \, \bar{x}_2 = x_2, \, \bar{x}_+ = x_+$ in $\partial\Omega$) under the same existence condition of the previous case was proved. Furthermore, if the fractional changes in net growth of ceftiofur-sensitive and resistant E. coli, ($E_s$ and $E_r$, respectively), have opposite signs and the inequality (4.5) is satisfied, the system (2.11) has another coexistence equilibrium $\bar x$ in the interior of the set $\Omega$. In addition, we have shown that if $E_r > 0$ then $\bar x_1$ is l.a.s. in $\Omega$, which suggests that if the fractional changes in net growth of ceftiofur-resistant E. coli are positive then the population of ceftiofur-sensitive commensal enteric E. coli will be cleared while ceftiofur-resistant commensal enteric E. coli will reach carrying capacity $N_{\max}$. If the fractional changes in the net growth of ceftiofur-sensitive E. coli are positive and the plasmid transfer rate is less than in-flow and outflow rates, then the population of ceftiofur-resistant commensal enteric E. coli will be cleared while ceftiofur-sensitive commensal enteric E. coli will reach carrying capacity $N_{\max}$. If both fractional changes in the net growth of ceftiofur-sensitive and resistant E. coli are positive and the condition (3.28) is satisfied, then $\bar{x}_+$ is l.a.s. In other words, the bacterial populations will reach a steady state in $\partial\Omega$. We could not determinate the stability of $\bar{x}$. However, for a certain set of parameters we have numerically verified its stability (see Figure 3). Newly, not only were we able to corroborate the results of Volkova et al. ^[24], but we also supplemented them by determining conditions that govern the outcome of the infection. In addition, in certain cases the ranges of the parameters were determined. Finally, the fact that it has been proven that (2.11) does not have periodic solutions implies that the solutions will always reach an equilibrium solution, which is a required condition in PD.

7.   Conclusions

This work provides the mathematical theoretical framework for modeling developed by Volkova et al. ^[24] complementing the robustness of the model (2.11) and allowing:

● To verify the properties obtained by Volkova et al., but also evidencing other properties of bacterial dynamics that had not been previously shown.

● To develop specific relationships between the parameters that determine when a) ceftiofur-sensitive enteric E. coli tends to be cleared, b) ceftiofur-resistant enteric E. coli tends to be cleared, c) ceftiofur-sensitive and resistant commensal enteric E. coli tends to be cleared, and d) neither of the two populations are eliminated.

● To assess the impact of PD and PK on infection outcome.

Use of AI tools declaration

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

Acknowledgments

Eduardo Ibargüen-Mondragón and Miller Cerón-Gómez acknowledge support from project No 34-08/03/2022 (2466) financed by Vicerectoría de Investigaciones e Interacción Social (VIIS)-Universidad de Nariño. M. Victoria Otero-Espinar is partially supported by the Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación (Spain), grant PID2020-115155GB-I00.

Conflict of interest

The authors declare there are no conflicts of interest.