Model descriptionA system of ordinary differential equations was developed to track the dynamics of the antimicrobial drug concentrations and the bacterial populations in the lungs and gut. The model flowchart is presented in (Fig. 1).Fig. 1Schematic representation of the mathematical model depicting enrofloxacin concentrations in a steer’s gastrointestinal and respiratory tracts, and its effects on both E. coli and P. multocida bacteria populations. The compartments S, P, L, and C represent antimicrobial drug concentrations in subcutaneous tissue, plasma, lung, and colon, respectively. \({\text{S}}_{\text{e}}\), \({\text{R}}_{\text{e}}\), \({\text{S}}_{\text{m}}\), and \({\text{R}}_{\text{m}}\) represent the susceptible and resistant subpopulation of E. coli and P. multocida. Parameters are described in the text.Enrofloxacin is administered to the subcutaneous compartment S with a dose denoted as \({\text{S}}_{0}\), given i times at regular intervals T. This input process is mathematically represented by the Dirac function \(\delta\)(t-iT). The drug is eliminated from the subcutaneous compartment S at a rate of k. The antimicrobial drug enters the plasma at a rate of \(\beta\) k. Within the plasma, the total concentration of both bound and unbound antimicrobial drug is denoted as P. The drug is then transferred from the plasma to the colon and lungs at a net transfer rate of \(\alpha \) and \(\gamma\), respectively. Finally, the drug levels in the colon C and the lungs L are eliminated at rates of \(\vartheta\) and \(\delta\), respectively. All the mentioned rates apply to both enrofloxacin and its metabolized form, ciprofloxacin, combined.The bacterial populations, E. coli and P. multocida, grow logistically with net growth rates of \(\sigma\) and \(\lambda\), respectively. The maximum carrying capacity for the E. coli and P. multocida are denoted as \({\text{N}}_{\text{emax}}\) and \({\text{N}}_{\text{mmax}}\). The net growth rate of a bacterial strain is considered to represent the bacterial fitness at the within host level15. We modeled the fitness costs associated with resistance by including a fractional reduction in the net growth rates, represented by c and p for the resistant E. coli (\({\text{R}}_{\text{e}}\)) and resistant P. multocida (\({\text{R}}_{\text{m}}\)) subpopulations, respectively.The pharmacodynamic effects of enrofloxacin are modeled by a sigmoidal function representing the saturation of the death rate as a function of the antimicrobial drug concentration5. We assume that the concentration required to produce half of the maximum death effect is greater for resistant bacteria than sensitive bacteria, hence \({\text{C}}_{\text{r50}}\) > \({\text{C}}_{\text{s50}}\) and \({\text{L}}_{\text{r50}}\) > \({\text{L}}_{\text{s50}}\).The model described above is represented by the system of differential equations given below:$$\frac{dS}{dt}= \sum_{i=0}^{n}{S}_{0 }{\delta }_{0} \left(t-iT\right)-KS$$$$\frac{dP}{dt}=\beta KS-\left(\alpha +\gamma \right)P$$$$\frac{dC}{dt}=\alpha P- \vartheta C$$$$\frac{dL}{dt}= \gamma P- \delta L$$$$\frac{d{S}_{e}}{dt}= \sigma \left(1-\frac{{S}_{e} + {R}_{e}}{{N}_{emax}}\right){S}_{e}-d \frac{C}{C+ {C}_{s50}} {S}_{e}$$
(1)
$$\frac{d{R}_{e}}{dt}=\left(1-c\right) \sigma \left(1- \frac{{S}_{e}+{R}_{e}}{{N}_{emax}}\right) {R}_{e}-d \frac{C}{C+{C}_{r50}} {R}_{e}$$$$\frac{d{S}_{m}}{dt}= \lambda \left(1- \frac{{S}_{m}+{R}_{m}}{{N}_{mmax}}\right){R}_{m}- \eta \frac{L}{L+{L}_{s50}} {S}_{m}- \phi {S}_{m}$$$$\frac{d{R}_{m}}{dt}=\left(1-p\right)\lambda \left(1-\frac{{S}_{m}+ {R}_{m}}{{N}_{mmax}}\right){R}_{m}- \eta \frac{L}{L+{L}_{r50}} {R}_{m}- \phi {R}_{m}$$with initial conditions S (0) = 0 mg/kg, P (0) = 0 mg/kg, C (0) = 0 μ/ml, L (0) = 0 μ/ml, \({\text{S}}_{\text{e}}\) = 400,000 colony forming units (CFUs), \({\text{R}}_{\text{e}}\) = 1000 CFUs, \({\text{S}}_{\text{m}}\) = 40,000 CFUs, \({\text{R}}_{\text{m}}\) = 1000 CFUs.Data and model parameterizationThe experimental data used to partially parameterize the model were previously described in Foster et al.16,17 and was approved by the North Carolina State University Institutional Animal Care and Use Committee. Moreover, Foster et al.16,17 affirmed adherence to animal welfare and ARRIVE guidelines in their study, demonstrating meticulous attention to ethical standards and methodological rigor. The data for fitting the PK components of the model were collected from two experimental studies where a cohort of twelve steers were monitored following two approved dosing scenarios: a single dose of enrofloxacin (12.5 mg/kg), and a single dose of enrofloxacin (7.5 mg/kg), both administered subcutaneously. After the subcutaneous administration, enrofloxacin concentration in plasma, colon, and interstitial fluid were measured over time. Additionally, E. coli concentration in the gut, measured as CFU/ml, and its minimum inhibitory concentration (MIC) was determined according to the established guidelines of the Clinical and Laboratory Standards Institute6. To capture the changes of overall bacteria population in the model, we denoted the sub-population of E. coli and P. multocida with MICs above the epidemiological cut-off (ECOFF) value as resistant (\({\text{R}}_{\text{e}}\) and \({\text{R}}_{\text{m}}\)), and those below the cut-off point were defined as susceptible (\({\text{S}}_{\text{e}}\) and \({\text{S}}_{\text{m}}\)). The ECOFF value for enrofloxacin on E. coli and P. multocida is the same, 0.125 \({\upmu}\)g/ml18,19.Parameters that can be derived uniquely from the data are considered identifiable20,21. Structural identifiability is a theoretical way to determine whether the within-host ODE model’s parameters are identifiable from the noise-free observations without the actual data22. A structurally unidentifiable parameter obtained from fitting is not valid for further analysis of the ODE model20,21,23. In our model, we investigated the structural identifiability of the parameters before estimating their values. We used a user friendly and universally accessible web application COMBOS24 for checking structural identifiability. COMBOS uses a Grobner-based computation to check the structural identifiability of the parameters24.We used the software Monolix 2020 R1 (Lixoft, Antony, France) to estimate the PK parameters from the data. We independently fitted the model to the concentration of antimicrobial drugs in the plasma, colon, and lung, and susceptible and resistant E. coli population in the feces for each steer in the treatment groups. The parameter estimation was done using stochastic approximation expectation maximization (SAEM). SAEM is a technique that combines maximum likelihood estimation with stochastic approximation to estimate conditional expectations25. We performed a maximum of 200 Monte Carlo runs, each with 10,000 iterations. The parameters k, \({\text{C}}_{\text{s50}}\), \({\text{L}}_{\text{s50}}\), \({\text{L}}_{\text{r50}}\) were fixed based on previous studies, as shown in Table 1. In our PK-PD model, we assumed all the parameters had log-normal distributions. Also, we let \(\sigma\), d, and \({\text{N}}_{emax}\) vary across study populations to capture the dynamics of E. coli at the individual level. The reason to estimate those parameters is that commensal E. coli populations typically vary across individuals and can be highly dynamic26. Using the data from Foster et al.17, we estimated \(\beta\), \(\alpha \), \(\vartheta\), \(\sigma\), and \({\text{N}}_{emax}\). Also, by using data from Foster et al.16, we estimated γ, d and δ. Then \({\text{N}}_{mmax}\), \(\lambda\), c, \({\eta}\), p, and \({\phi}\) values were calibrated to a scenario in which the disease-causing P. multocida bacteria was cleared in the model.Table 1 Parameter values from literature, calibration, and estimation. Here, aindicates that parameters value used to clear the infection scenario.The model was run in R software version 4.3.1. The integrator used in our model is lsoda from the deSolve package (version 1.40) in R software version 4.3.1. lsoda is a robust numerical solver that automatically switches between non-stiff and stiff integration methods.Simulated scenariosIn our study, we simulated the impact of three approved enrofloxacin dosing scenarios and some non-approved proposed scenarios on the treatment efficacy against one of the causative agents of BRD, P. multocida. Specifically, in the case of approved scenarios, we simulated 12.5 mg/kg and 7.5 mg/kg as a single dose, and doses of 5 mg/kg given 24 h apart for three days, all administered subcutaneously. In case of proposed non-approved dosing scenarios, we simulated doses of 12.5 mg/kg, 7.5 mg/kg, 6.25 mg/kg, and 3.75 mg/kg, administered 24 h apart for 2 days. Additionally, we examined doses of 4.15 mg/kg, given subcutaneously 24 h apart for three consecutive days. Our investigation aimed to simulate the effects of these drug dosing regimens and gain insights into their therapeutic efficacy, treatment cost, and bacterial resistance outcomes.Our study assessed the treatment expenses associated with different dosing scenarios. To quantify each treatment scenario’s cost, we considered the average cost of medications and labor and the productivity losses resulting from the treatment process28. If the infection was not fully resolved (defined in the model as presence of bacteria in the lung after treatment), we added an additional treatment in our analysis.Uncertainty and sensitivity analysisTo identify the PK and PD parameters that have the most influence on the bacterial resistance within the lung and gut sites, \({\text{R}}_{\text{m}}\) and \({\text{R}}_{\text{e}}\), we performed a variance-based sensitivity analysis. This Sobol sensitivity approach uses model output uncertainty to describe the variance29,30. We have chosen to vary the following PK and PD parameter set, q = (k, \(\beta\), \(\alpha \), \(\gamma\), \(\vartheta\), \(\delta\), \({\text{C}}_{\text{s50}}\), \({\text{C}}_{\text{r50}}\), \({\text{L}}_{\text{s50}}\), \({\text{L}}_{\text{r50}}\)). As a model output, we evaluated the uncertainty of both the susceptible and resistant of E. coli and P. multocida bacterial populations. We calculated the first-order and total-order Sobol indices of the parameter set q for the resistant bacteria31,32,33. The basic algorithm for the calculation of the indices is as follows:
Providing lower and upper bound values for each parameter of q as an input, we obtain N possible values uniformly distributed for each of the parameters of q.
Taking the distributions of the set of parameters q as an input variable, the within-host model (1) generates N = 2000 set of solutions, also called uncertainty distribution of solutions. In this study, we focus on the uncertainty distribution of states \({\text{S}}_{\text{e}}\)(t), \({\text{R}}_{\text{e}}\)(t), \({\text{S}}_{\text{m}}\)(t), and \({\text{R}}_{\text{m}}\)(t).
As mentioned, we are only concerned with calculating the Sobol indices of \({\text{R}}_{\text{m}}\) and \({\text{R}}_{\text{e}}\). Here, we calculated the overall variance of the resistant using deviation from the average value of the time-dependent uniform distribution of the resistant. Let the overall variance of \({\text{R}}_{\text{m}}\) be denoted by \({\text{V}}_{{\text{R}}_{\text{m}}}\).
The first-order variance \({\text{V}}_{{\text{S}}_{\text{i}}}\) of \({\text{R}}_{\text{m}}\) is calculated with respect to an individual parameter \({\text{q}}_{\text{i}}\) of the set q.
The total-order variance \({\text{V}}_{{\text{T}}_{\text{i}}}\) of \({\text{R}}_{\text{m}}\) is calculated with respect to all parameters except the individual parameter \({\text{q}}_{\text{i}}\) of the set q.
The time-dependent first-order Sobol index of the parameter \({\text{q}}_{\text{i}}\) for \({\text{R}}_{\text{m}}\) is calculated with the formula \({\text{S}}_{\text{i}}\) = \(\frac{{\text{V}}_{{\text{S}}_{\text{i}}}}{{\text{V}}_{{\text{R}}_{\text{m}}}}\) and the time-dependent total-order Sobol index of the parameter qi for \({\text{R}}_{\text{m}}\) is calculated with the formula \({\text{T}}_{\text{i}}\) = 1 − \(\frac{{\text{V}}_{{\text{T}}_{\text{i}}}}{{\text{V}}_{{\text{R}}_{\text{m}}}}\).
The sum of the first-order indices of all parameters q such as \({\sum }_{\text{i}}{{\text{S}}}_{\text{i}}\) is less than 1. Similarly, \({\sum }_{\text{i}}{{\text{T}}}_{\text{i}}\) also is less than 1. If \({\text{S}}_{\text{i}}\) or \({\text{T}}_{\text{i}}\) of a parameter is close to 1, then that parameter significantly influences the state \({\text{R}}_{\text{m}}\). On the other hand, if \({\text{S}}_{\text{i}}\) or \({\text{T}}_{\text{i}}\) of a parameter is close to 0, that parameter has low influence on the state \({\text{R}}_{\text{m}}\).
As the indices \({\text{S}}_{\text{i}}\)(t) and \({\text{T}}_{\text{i}}\)(t) of each parameter of set q are time dependent, we plotted the curves of the indices over time T = 150 h
The curve thus created by the indices of each parameter of q over time 150 h is integrated to get area using the Riemann sum. These areas help us to compare the level of sensitivity of parameters of set q for \({\text{R}}_{\text{m}}\) for the complete simulation period.
Similarly, we repeated the similar process of finding Sobol indices and integrated area for \({\text{R}}_{\text{e}}\text{.}\)
We used the “Sensobal package” (version 1.1.5) of R software version 4.3.1 to calculate the Sobol indices33,34.The input bounds for q, the generated uniform distribution of states \({\text{S}}_{\text{e}}\), \({\text{R}}_{\text{e}}\), \({\text{S}}_{\text{m}}\), and \({\text{R}}_{\text{m}}\) as well as the results of Sobol indices of the model parameters q for resistant bacteria \({\text{R}}_{\text{m}}\) and \({\text{R}}_{\text{e}}\) are presented in the results section. We utilized the proposed value for each parameter in Table 1 and select upper and lower bounds not more than 3 standard deviations from each parameter of set q. Each parameter’s upper and lower bounds of q are in Table 2Table 2 The lower and upper bound of the set of parameters q, which are applied to perform the Sobol sensitivity analysis in the ‘Sensobal package’ (version 1.1.5) of R software version 4.3.1. We iterated the simulation 2000 times and used a uniform distribution of each parameter of set q.