Modeling Klebsiella pneumonia infections and antibiotic resistance dynamics with fractional differential equations: insights from real data in North Cyprus

This study presents an enhanced fractional-order mathematical model for analyzing the dynamics of Klebsiella pneumonia infections and antibiotic resistance over time. The model incorporates fractional Caputo derivative operators and kernel, to provide a more comprehensive understanding of the complex temporal dynamics. The model consists of three groups: Susceptible ( S ), Infected ( I ), and Resistant ( R ) individuals, each controlled by a fractional differential equation. The model represents the interaction between infection, recovery from infection, and the possible development of antibiotic resistance in susceptible individuals. The existence, uniqueness, stability, and alignment of the model’s prediction to the observed data were analyzed and buttressed with numerical simulations. The results show that imipenem has the highest efficacy compared with ertapenem and meropenem category drugs. The estimated reproduction number and reproduction coefficient illustrate the potential impact of this model in improving treatment strategies, while the memory effects highlight the advantages of fractional differentiation. The model predicts an increased possibility of antibiotic resistance despite effective treatment, suggesting a new treatment approach.


Introduction
Klebsiella pneumoniae is a Gram-negative bacterium that has posed a serious public health challenge in the world.It can cause different types of infections, from urinary tract infections to life-threatening pneumonia, and with its growing antibiotic resistance.This necessitates a deep understanding of how it causes infections as well as develops resistance [1][2][3][4][5].According to [6], K. pneumonia is by far the most significant contributor to global antibiotic resistance particularly within institutions where there are always high-risk (HiR) epidemic strains that are generally resistant to multiple antibiotics.Consequently, these strains cause limited treatment options leading to more deaths, complications, and increased costs of healthcare.Klebsiella pneumoniae, being a hospital-acquired pathogen, continuously receives exposure to different antibiotics leading to evolutionary pressure aiming for further positive mutations.The drug resistance pattern also comprises a wide range of antibiotic resistance genes (ARGs) that encompass both chromosomal and plasmid-encoded gene cassettes.Plasmid-mediated resistome and transposons facilitate horizontal transfer to other bacterial strains, contributing to the spread of resistance determinants.Authors of [7] emphasize that managing antimicrobial resistance in multi-drug-resistant K. pneumonia (MDR-KP) poses a significant challenge for clinicians.The optimal treatment approach for MDR-KP infections remains uncertain, prompting the exploration of various combination therapies involving antibiotics like meropenem, colistin, fosfomycin, tigecycline, and aminoglycosides.New antimicrobials targeting MDR-KP are in different stages of clinical research and development.They emphasize the need for coordinated strategies, infection control, and prudent antimicrobial use to limit MDR-KP spread and improve treatment outcomes.Combination therapies and newer drugs like ceftazidime-avibactam are explored for their potential in treating MDR-KP infections.Mathematical models have shown to be useful in describing complex phenomenons, which have led to the development and optimization of treatment strategies and public health interventions [8][9][10][11][12][13].Authors of [14] investigated the future clinical scenario of K. pneumonia, a significant pathogen causing healthcare-associated infections with antibiotic resistance.They used an integer mathematical modeling approach (SIS-type model) with retrospective medical data to predict the spread of extended-spectrum beta-lactamase (ESBL)-producing K. pneumonia.The study revealed that ESBL-producing strains will possibly exceed non-ESBL strains in around 70 months, indicating the urgent need for preventive measures.Sensitivity analysis highlighted the critical role of antibiotic use in the development of ESBL-producing strains.The study emphasizes proper antibiotic use, infection control, and surveillance to mitigate the rise of antibiotic-resistant K. pneumonia.Integer ordinary differential equation models (ODE), while widely used in epidemiology, have some limitations when applied to the study of complex infectious diseases such as K. pneumonia infections and antibiotic resistance.A significant limitation is their inability to fully represent the complex dynamics and differentials of non-integer order observed in real systems.Integer-order models typically assume constant rates of transmission, recovery, and resistance development and ignore possible differences that may be critical in the context of rapid pathogen evolution.Furthermore, these models can have difficulty describing population diversity because they often rely on averages and assumptions of homogeneous mixtures.Integer ODE methods can generalize dynamics, leading to inaccurate predictions of the increase and spread of antibiotic-resistant strains.The model proposed in this study incorporates fractional differential equations to better explain the complexity and dynamics of infectious disease dynamics by proposing a more flexible framework to address some of these limitations.Many scientists [15][16][17][18][19] have employed the Susceptible-Infected-Resistant (SIR) framework in studying infections and drug effectiveness.However, no author in the literature has narrowed their research to K. pneumonia using (SIR) fractional-order dynamics to ascertain and compare drug effectiveness as well as predict future scenarios of the susceptible, infected, and resistant population.
In this study, we introduce an advanced mathematical model that integrates fractional Caputo derivative operators into the traditional Susceptible-Infected-Resistant (SIR) framework.The use of fractional-order differentiation improves the model's competence to capture intricate chronological behaviors in the context of K. pneumonia infections and antibiotic resistance.Our model consists of three groups: Susceptible (S), Infected (I), and Resistant (R), each governed by fractional-order differential equations.Within the model, the Susceptible group accounts for individuals susceptible to infection and their recovery rates.In the Infected group, we study the dynamic interaction among susceptible individuals becoming infected, recuperating from infection, and possibly developing antibiotic resistance.The Resistant group demonstrates people who have become immune to antibiotics.
Our model uses fractional-order differentiation as well and this facilitates more sophisticated investigations of the sequential dynamics which may not be apparent using the traditional integer-order differential equations.It is a useful research tool for epidemiologists and medical mathematicians who aim at understanding K. pneumonia infections and antibiotic resistance better to improve disease management and treatment approaches.What distinguishes this study is the incorporation of fractional-order dynamics into the Susceptible-Infected-Resistant (SIR) framework for K. pneumonia.As a result, it explores more superior drug effectiveness future scenarios and new perspectives on how to combat antibiotic-resistant K. pneumonia.The method used in this research involves employing fractional differential equations within the Susceptible-Infected-Resistant (SIR) framework to simulate K. pneumonia infections and antibiotic resistance dynamics because it has several benefits.One significant advantage is its capacity to capture complex and intricate behaviors in the dynamics of the susceptible, infected, and resistant populations, such as memory effect [20][21][22].The fractional-order dynamics provide a refined representation of the temporal evolution of the system, allowing for a closer approximation of real-world complexities.Furthermore, the model's increased sensitivity to changes in parameters and initial conditions, especially as the fractional order approaches 1, has the potential to more accurately represent the impact of subtle changes on infection spread and resistance development.However, as with any modeling approach, there are inherent shortcomings.The complexity caused by fractional differential equations can lead to computational costs and interpretability challenges.The need for accurate parameter estimates and validation against empirical data is becoming increasingly important, and model performance depends on the availability of precise input parameters.Furthermore, adopting and applying fractional-order dynamics may require a more comprehensive understanding by practitioners and researchers, which may limit their accessibility for practical applications.Balancing the benefits of increased realism with the computational and interpretive requirements are key considerations for the successful application of this approach in understanding and predicting K. pneumonia infections and antibiotic resistance dynamics.To formulate a mathematical model to describe the dynamics of K. pneumonia infections and antibiotic resistance over time, we use a simple compartmental model as shown in Eq. ( 1), with the description of the variables and parameters in Table 1.

Collection of microbiological data
The study comprised 937 strains of K. pneumonia isolated in the university hospital microbiology laboratory between January 2018 and July 2023.Antibiotic susceptibilities of these strains, obtained from hospitalized patients, were determined using the Biomérieux Vitek ® 2 automated system.The effectiveness of carbapenem group antibiotics, namely ertapenem, imipenem, and meropenem, on these strains was evaluated using a mathematical model.

The K. pneumonia fractional-order model
A fractional-order model of the Caputo type is employed to describe the disease dynamics, schematically shown in Fig. 1, and mathematically in system (1).The proof of the existence and uniqueness of the system is established using the Banach contraction principle, and the stability of the solution is shown using the Jacobian matrix.Numerical analysis was simulated using MATLAB implementing the Predictor Corrector Scheme [23,24].All variables and parameters are explained in Table 1.
where P(t) = S(t) + I(t) + R(t) and c 0 D α t is the Caputo fractional differential operator which is defined in Definition 3, with the following initial conditions: (1)  Fractional-order Caputo differential operator c 0 D α t is used in the model (1).It presents a more comprehensive view of how K. pneumonia infections and antibiotic resistance develop over time.This model is divided into three groups: Susceptible (S), Infected (I), and Resistant (R) which are governed by fractional differential equations.The model explores how susceptible individuals may get infected at ( β α ) rate, recover from the infection at ( δ α ) rate, and possibly develop antibiotic resistance over time at ( µ α ) rate.In the case of the Susceptible group, the rate of change c 0 D α t S(t) shows the number of people who can be infected and then recover concerning the number of infected ones (I).The Resistant group consists of individuals who have developed resistance towards antibiotics and recover at a rate ( θ α ) from infection.Complex dynamics are better stud- ied using this Caputo derivative-based model for K. pneumonia infections and antibiotic resistance, which takes into account non-integer order differentiation to capture more non-linear and complex behaviors.

Preliminaries
Definition 1 Caputo derivative [25] The Caputo derivative of order α ∈ (0, 1) of a sufficiently differentiable function f(t) is defined as follows: where Ŵ is the gamma function.
Definition 2 Gamma function [26] The gamma function Ŵ(z) is defined for Re(z) > 0 by the integral and f(t) be absolutely continuous on [0, ∞) .The Caputo fractional derivative of order α is defined as [27,28]: where f(t) is n times differentiable and Ŵ(x) is the Gamma function:

Definition 4
The Riemann-Liouville fractional integral of order α > 0 of a function f(t) is defined as [27,28]: This integral exists almost everywhere for any integrable function f(t).
The Riemann-Liouville integral and the Caputo fractional derivative operators satisfy the following property: Definition 5 The Mittag-Leffler function E σ (w) is defined as: Its generalized form E σ ,ς (w) is defined as: These functions are called Mittag-Leffler functions [29].

Definition 6 Laplace transform of the Caputo derivative [30]
The Laplace transform of the Caputo derivative D α t of order α ∈ (0, 1) of a function f(t) is defined as: where Lf (t)(s) is the Laplace transform of f(t) and f (0 + ) denotes the right-sided limit of f(t) at t = 0. Definition 7 Banach contraction principle [31] Let (X, d) be a metric space, and let T : X → X be a function.Then T is a Banach con- traction if there exists a constant 0 ≤ k < 1 such that for all x, y ∈ X,

Model analysis
Theorem 1 The solution to system (1) exists and is unique.
Proof We first show that system (1) is bounded [32,33]: The system solution (1) exists and is unique.

Proof
Using the Laplace transform method in Defintion 6 to solve Gronwall's inequality [34] with initial condition P(t 0 ) ≥ 0 , we get: where the series of the Mittag-Leffler functions E α (−θ α t α ) and E (α,i+1) (−θ α t α ) converge (see Definition 5).Therefore, the system (1) has a bounded solution and the biologically feasible region is given as: , Next, we rewrite system (1) as: where y(t) is the vector of dependent variables S(t), I(t), and R(t), and H(t, y(t)) is the vector of corresponding right-hand sides of the differential equations where L = (�A� + 1) , and Hence, H is uniformly Lipschitz continuous and bounded.We proceed to complete the proof for the uniqueness of the system (1).Let 0 < α < 1 , φ = [0, z * ] ⊆ R and ψ = �y(t) − y(0)� ≤ K , and h : φ × ψ → R be a continuous bounded function, such that |h(t, y)| ≤ M , since H is Lipschitz continuous.
Let LK < M , then there exists a unique y ∈ C α [0, z * ] for the initial value problem (2), where and hence a complete metric space.
Transforming the system (2) to the equivalent Volterra integral equation: Defining an operator H in E , as H : E → E: (2) It follows that: M �y − y * � , and from hypothesis LK M < 1.Therefore, as a consequence of the Banach contraction principle [35], E is a contrac- tion and has a unique fixed point [36].Hence, from the Picard-Lindelöf theorem [37], the system (1) has a unique solution.

Theorem 3 System (1) is locally asymptotically stable at the given positive equilibrium.
Proof To prove the local asymptotic stability of the system described by Eq. ( 1), we will analyze the system's equilibrium points and the stability of those points using linearization.The equilibrium points are those where the derivatives are zero, i.e., where Equilibrium points: To find the equilibrium points, we set c 0 D α t S(t) = 0 , c 0 D α t I(t) = 0 , and c 0 D α t R(t) = 0 .This leads to the following equations: From Eq. (3), we can solve for S eq : ( Substituting this into Eq.( 4), we can solve for I eq : Finally, substituting the values of I eq and S eq into Eq.( 5), we can solve for R eq : So, the equilibrium points are given by: Stability analysis [38]: To analyze the stability of the equilibrium points, we will linearize the system around these points.We will calculate the Jacobian matrix J of the system evaluated at the equi- librium point (S eq , I eq , R eq ): where Ṡ, İ, and Ṙ are the right-hand sides of your equations.Then, we will evaluate the eigenvalues of the Jacobian matrix J at the equilibrium point.
If all eigenvalues of J have negative real parts, then the equilibrium point is locally asymptotically stable.
Calculation of Jacobian matrix: Let us calculate the Jacobian matrix J for the given system.We will differentiate each equation with respect to S, I, and R: where 1, 2, 3 are the right-hand sides of your equations.Differentiating Eq. ( 1) with respect to S, I, and R , we get: Differentiating Eq. ( 4) with respect to S, I, and R , we get: S eq = β δ α I eq .
(S eq , I eq , R eq ) = β δ α I eq , Differentiating Eq. ( 5) with respect to S, I, and R , we get: Now, we can assemble the Jacobian matrix J: Next, we will evaluate the eigenvalues of this matrix at the equilibrium point (S eq , I eq , R eq ) to determine the stability.To compute the eigenvalues of the Jacobian matrix J at the equilibrium point (S eq , I eq , R eq ) , we will use the Jacobian matrix J that we derived previously.Now, we can calculate the eigenvalues by solving the characteristic equation: where I is the identity matrix, and is the eigenvalue we are solving for.
Substituting the values from the Jacobian matrix: Now, we can simplify and solve for by expanding the equation: So we get: Since all three eigenvalues have negative real parts, it indicates that the equilibrium point of the system (1) is locally asymptotically stable.Therefore, the system (1) is stable.

Reproduction number and coefficient
Reproduction number (R 0 ) [39]: The basic reproduction number, denoted as (R 0 ) , por- trays the expected number of secondary cases produced by one infected individual when introduced into a completely susceptible population.For our model, the linearized equations can be written as follows: Now, let us find the Laplace transforms of I(t) and R(t) by solving these equations: Since we are interested in the behavior near the disease-free equilibrium, we consider the case where I(s) and R(s) are both zero.Now, we can compute (R 0 ) using the next- generation matrix approach: Reproduction coefficient (R) [40]: The reproduction coefficient, denoted as R, portrays the expected number of secondary cases produced by one infected individual when introduced into a partially immune population.It is related to(R 0 ) as follows:

Sensitivity analysis
We want to compute the sensitivity of the final susceptible population S(T) to a change in the parameters β α and δ α .To compute the sensitivity coefficients, we first need to find the solution to the fractional differential equation for the baseline parameters and then for the perturbed parameters.
Baseline solution: Let us denote the baseline solution as S b (T ) for the baseline param- eters β b and δ α b .
R 0 = Sum of products of contact rates and fractions of the population that is susceptible Recovery rate , Perturbed solution (for β α ): We perturb the parameter β α by a small amount δ α β α to obtain a new equation: Solve this equation to obtain the solution S ′ (β)(T ) with the perturbed parameter Perturbed solution (for δ α ): Similarly, we perturb the parameter δ α by δ α β α to obtain a new equation: Solve this equation to obtain the solution S ′′ (δ α )(T ) with the perturbed parameter δ α b + δ α β α .Compute sensitivity coefficients: For β α : For δ α : These sensitivity coefficients S β α and S δ α represent how changes in the parameters β α and δ α influence the final susceptible population S(T).

Sensitivity of group I(t):
For group I(t), the fractional differential equation is: Sensitivity to β α : We perturb the parameter β α by δ α β α : The sensitivity coefficient for β α is calculated as: Sensitivity to δ α : We perturb the parameter δ α by δ α β α : The sensitivity coefficient for δ α is calculated as: .

Sensitivity of group R(t):
For group R(t), the fractional differential equation is: Sensitivity to β α : We perturb the parameter β α by δ α β: The sensitivity coefficient for β α is calculated as: Sensitivity to δ α : We perturb the parameter δ α by δ α β α : The sensitivity coefficient for δ α is calculated as: These sensitivity coefficients represent how changes in the parameters β α and δ α influ- ence the final populations of groups I(t) and R(t).

Drug effectiveness analysis
To integrate a drug effectiveness analysis with model (1), we will introduce variables to represent the effectiveness indices of each drug category and modify the infection and recovery rates based on these indices.Let us denote the effectiveness indices as follows: We can update the model equations in (1) as follows to incorporate the drug effectiveness: We can compute the efficiency indices for each drug class as follows: In Eq. ( 22), the recovery rate ( δ α ) is adapted based on the efficiency catalogs of each drug category.To regulate the most effective drug category over time, we can calculate the effectiveness indices ( E 1 , E 2 , E 3 ) based on the data.Then, we can use these efficiency indices in the adapted model equations to simulate the dynamics of K. pneumonia (17) .
infections and antibiotic resistance over time.The drug category with the highest effectiveness index at a given time will be the most effective in controlling the infection.This method allows us to implement the analysis of drug efficiency into the original model, providing an all-inclusive understanding of the disease dynamics while considering the influence of diverse treatments.

Numerical analysis
To numerically solve systems ( 1) and ( 22), we consider the initial value problem in (1): Employing the Riemann-Liouville integral operator in Definition 4, we get: Substituting t = t m and t = t m+1 into Eq.( 23) and subtracting the obtained equa- tions, we get: where t j = jh , j = 0, 1, . . ., N , and h = T f /N portrays the step size.We approximate f (s, y(s)) on [t m , t m+1 ] using two-step Lagrange polynomial interpolation: where y k = y(t k ) .Using (24) and ( 25), we have: Employing integration by parts, (26) becomes: Since y m+1 appears on the right side of (27), this equation is implicit, requiring the prediction of y m+1 as y p m+1 .Therefore, (27) serves as a corrector formula.In (24), we use the rectangle rule for the integral part, resulting in the following predictor formula: where Therefore, the numerical formula for system (1) is as follows: The predictor formula: The corrector formula:

Results and discussion
The analytic computations show that the model (1) has a unique and stable solution, and the numerical analysis further buttresses these results.According to Fig. 2, our computations demonstrate a strong fit between model predictions and the empirical data.Our model correctly predicts the average of the observed annual data.Based on our model predictions, Meropenem's effectiveness index is negative in Figs. 3 and 6.For instance, this could indicate that Meropenem would not be best suited to treating K. pneumonia according to such assumptions about what works within such models.The efficiency ratios obtained from the modeling process give insight into which drug category are most effective against K. pneumonia infections.Imipenem emerged as the most efficacious antibiotic, showcasing a positive effectiveness index as depicted in Figs. 4 and 6, (28) while Ertapenem is deemed the least effective.It suggests that Imipenem should be prioritized for treatment.The combination of its lower transmission rate, faster recovery, reduced mortality, and the development of resistance positions Imipenem as a potent tool in managing and controlling the disease.Figures 5 and 6 indicate that Ertapenem  has a significantly negative effectiveness index, thus showing that it is not effective in reducing infection.Therefore, Ertapenem may be inappropriate for treating this infection under these circumstances because its efficacy is very low (Figs. 7, 8, 9, and 10).Conversely, the model shows a positive effectiveness index for Imipenem (as compared with other antibiotics), which implies that it will reduce infections in 2023 as per the model hypothesis.Thus, there exists a better choice of Imipenem than Meropenem and Ertapenem having favorable contributions toward reducing infections.The reproduction number and coefficient ( R 0 and R) depicted in Fig. 11 displays the reproduction number and coefficient as R 0 and R, which provide valuable insights into possible K. pneumonia infections and antibiotic resistance within a population.The baseline estimate of transmission potential is R 0 while R considers immunity plus interventions.These measures are essential in evaluating whether public health controls, vaccination programs, and antibiotic strategies are working effectively against this organism.A value of between 1.02 and 1.25 for R 0 indicates that the illness has low to moderate potential for sus- tained transmission within the population (Fig. 11).It can cause local outbreaks but is not highly contagious.Therefore, it cannot spread easily and quickly across populations.Other calculated values of R ranged from between 0.02 to 0.1 demonstrating that there could be additional infections despite partial immune response or drug-resistant populations; however, the infection rates are extremely slow with limited spread pathways implying reduced fast trackability among people.
Sensitivity coefficients calculated for each group, (S) Sensitive, (I) Infected, and (R) Resistant, provide detailed information about the response of the system It was observed for the sensitive group that infection rate ( β α ) and intensity ( δ α ) were  (10).b Sensitivity coefficient of δ α to S(t) in (11).c Sensitivity coefficient plot of β α and δ α to S(t) in ( 10) and (11) affected which stands out as shown in Fig. 12.A sensitivity coefficient close to 0.3 for β α indicated that a small increase in infection rate resulted in a corresponding increase in the number of susceptible individuals, emphasizing significance if the spread of the disease is suppressed.On the other hand, a sensitivity coefficient of 1.8 for δ α indicated that small perturbations have a more pronounced effect of increasing sensitivity in the group a for the infected case; changes in β α and in δ α had a notice- able effect, with sensitivity coefficients of − 0.2 and − 2.8, respectively, as shown in Fig. 13.This indicated that infection rates resulted in a corresponding decrease in infected individuals, emphasizing the importance of monitoring disease transmission, whereas higher recovery rates resulted in human infection rates decreasing dramatically, emphasizing the importance of effective treatment in the management of the disease.Differences in factors such as the development of antimicrobial resistance ( µ α ) and recovery rates had a significant effect on individual infection rates of resist- ance ( θ α ).The rate of development of antibiotic resistance against the resistant group µ α was found to reduce the number of resistant individuals (R) with a sensitivity coef- ficient of − 0.2 as depicted in Fig. 14.This highlights the importance of controlling the development of antibiotic resistance to maintain the effectiveness of treatments.These parameters' sensitivity can give insight into how K. pneumonia infections and antibiotic resistance change with model parameters.This helps to improve or even re-engineer disease control strategies of antibiotic resistance, leading to a better understanding of the system's behavior as well as possible impacts caused by parameter fluctuations in time.

Conclusion
We proposed a fractional-order model to simulate K. pneumonia infections and antibiotic resistance using real data from Northern Cyprus.Comparatively, this is a model that has several advantages over traditional integer-order models.Firstly, the fractionalorder models allow for more long-term memory of bacterial infection and antibiotic resistance dynamics accurately represented than traditional integer order models do.This is much closer to reality given the fact that such system behaviors change with time.Also, anomalous dynamical phenomena observed in K. pneumonia infections are effectively described by the model through the incorporation of fractional derivatives.Furthermore, the complex spatial dynamics and processes that the model can cover are improved.It also helps in simulating non-local interactions between bacterial strains, antibiotics, and host immune responses which contributes to a more comprehensive understanding of system dynamics.Moreover, as regards real data obtained from North Cyprus, the fractional-order model demonstrates an improved fit with the specific parameters or trends for K. pneumonia infections including antibiotic resistance patterns which indicate high efficiency for monitoring the spread out process of this infection in a region like North Cyprus.However, unlike other studies that modeled bacterial infection and antibiotic resistance dynamics using ordinary differential equations, our approach uses fractional differential equations with real-world data that completely revolutionize how modeling disease transmission using ODEs looks Fig. 14 a Sensitivity coefficient of β α to R(t) in (19).b Sensitivity coefficient of δ α to R(t) in (21).c Sensitivity coefficient plot of β α and δ α to R(t) in ( 19) and (21) like today.Future recommendations would include a further study designed to expand its applicability to other geographic regions and infectious diseases.Additionally, the model can be expanded to include factors such as environmental influences, treatment strategies, and evolutionary dynamics, which will provide valuable insights for disease surveillance and public health intervention.

Fig. 2 aFig. 3 a
Fig. 2 a Observed data versus model's prediction.b Observed data plot over model's prediction mean projection

Fig. 4 a
Fig. 4 a Imipenem sensitive and resistance trend.b Imipenem effectiveness trend

Fig. 11 aFig. 12 a
Fig. 11 a Reproduction number (6) trend for various values of α. b Reproduction coefficient (7) trend for various values of α

Table 1
Variables and parameters in the K. pneumonia infection and antibiotic resistance model