Pharmacokinetic/Pharmacodynamic Anesthesia Model Incorporating psi-Caputo Fractional Derivatives

We present a novel Pharmacokinetic/Pharmacodynamic (PK/PD) model for the induction phase of anesthesia, incorporating the $\psi$-Caputo fractional derivative. By employing the Picard iterative process, we derive a solution for a nonhomogeneous $\psi$-Caputo fractional system to characterize the dynamical behavior of the drugs distribution within a patient's body during the anesthesia process. To explore the dynamics of the fractional anesthesia model, we perform numerical analysis on solutions involving various functions of $\psi$ and fractional orders. All numerical simulations are conducted using the MATLAB computing environment. Our results suggest that the $\psi$ functions and the fractional order of differentiation have an important role in the modeling of individual-specific characteristics, taking into account the complex interplay between drug concentration and its effect on the human body. This innovative model serves to advance the understanding of personalized drug responses during anesthesia, paving the way for more precise and tailored approaches to anesthetic drug administration.


Introduction
Pharmacokinetic/Pharmacodynamic (PK/PD) modeling is a mathematical approach used in pharmacology to study the relationship between drug concentrations (Pharmacokinetics) and their effects on the body (Pharmacodynamics) [31,37].The PK/PD models help researchers and clinicians to understand how drugs are absorbed, distributed, metabolized, and eliminated from the body [7].
The PK/PD models integrate Pharmacokinetic and Pharmacodynamic data to characterize the time course of drug action [29] .These models can be simple or complex, depending on the drug's characteristics and the purpose of the modeling.The parameters of these models were fitted by Schnider et al. in [36].Some common types of PK/PD models include: 1. Linear PK/PD models.The basic structure of a linear PK/PD model consists of two main components: the Pharmacokinetic component, which describes the drug concentration-time profile in the body, and the Pharmacodynamic component, which relates to the drug concentration to the observed effect [16].
2. Non-linear PK/PD models.These models can incorporate various components, such as saturable drug elimination, receptor binding kinetics, and indirect response models.These models provide a more accurate representation of the concentration-effect relationship and can be used to optimize dosing regimens and predict drug responses in different populations [11].
3. Mechanistic PK/PD models.These models incorporate known biological mechanisms, such as drug-receptor interactions, enzyme kinetics, or signal transduction pathways.They provide a more detailed representation of drug action but require more data and knowledge about the underlying biology [33].
4. Population Pharmacokinetic models.These models are used to describe the time course of drug exposure in patients and to investigate sources of variability in patient exposure.They can be used to simulate alternative dose regimens, allowing for informed assessment of dose regimens before study conduct [10,32].
In recent years, the field of fractional derivatives has emerged as a promising approach to model and understand complex biological processes characterized by non-integer order dynamics [21,27,35].This unique mathematical framework has found diverse applications in various areas of biology, where traditional integer-order calculus falls short in capturing the intricacies of these systems [6,25].
One prominent field where fractional derivatives have made significant contributions is Neurobiology.By employing fractional calculus, researchers have been able to delve into the dynamics of neural systems with a greater level of realism.This includes modeling the behavior of neurons, synaptic transmission, and the propagation of nerve impulses.The incorporation of fractional derivatives enables the consideration of memory effects and non-local behavior, providing a more accurate representation of neural processes [28].
Another area where fractional derivatives have shown their efficacy is in Biomedical Signal Processing.Biomedical signals, such as electroencephalograms (EEG), electrocardiograms (ECG), and blood pressure signals, are often complex and exhibit non-integer order dynamics.By utilizing fractional order filters and operators, meaningful information can be extracted from these signals, leading to improved analysis and interpretation.Fractional calculus has proven to be a valuable tool in enhancing our understanding of these vital physiological signals [15].
Furthermore, the field of Cancer Modeling has witnessed the application of fractional derivatives [34].Tumor growth and the intricate interactions between cancer cells and the immune system present complex dynamics that can be effectively captured using fractional order models.By incorporating non-local effects and memory into the modeling process, fractional derivatives provide a comprehensive framework for studying cancer progression.This approach has the potential to shed light on the underlying mechanisms and aid the development of novel therapeutic strategies [2,39].
The field of Pharmacokinetics and Pharmacodynamics plays a crucial role in understanding the behavior of drugs in biological systems [37].Traditional PK/PD models have predominantly relied on integer-order derivatives to describe various processes involved in drug absorption, distribution, metabolism, and elimination [13,31,43].However, these models often fall short in accurately capturing the complexity of Pharmacokinetic behaviors [8,19].
Fractional calculus offers a promising alternative by providing a more precise representation of these intricate dynamics [38].A number of studies have shown that certain drugs follow an anomalous kinetics that can hardly be represented by classical models.Indeed, fractional-order pharmacokinetics models have proved to be better suited to represent the time course of these drugs in the body and also they can describe memory effects and a power-law terminal phase.Therefore, they give rise to more complex kinetics that better reflects the complexity of the human body.In [19], a fractional one-compartment model with a continuous intravenous infusion is considered, where it allows to determine how the infusion rate influences the total amount of drug in the compartment.Moreover, in the case of multiple dosing administration, recurrence relations for the doses and the dosing times that also prevent drug accumulation are presented.Hence, in [9], a PK model was introduced employing a fractional-order approach akin to mammillary dynamics.This model was specifically designed to incorporate considerations of tissue entrapment, thereby altering the anticipated drug concentration profiles.The proposed model shows a limitation in data fit profiles, without transgressing the fundamental principles of mass balance and physiological states.The mathematical study of the amount of drug administered as a continuous intravenous infusion or oral dose for fractional-order mammillary type models is investigated in [30].
For pharmacokinetic and pharmacodynamic models, the first fractional one was introduced in [40].By utilizing the Caputo fractional derivative, the authors presented a fresh perspective on the intricate relationships between drug dosages, absorption rates, and therapeutic outcomes.Since then, some few applications of fractional PK/PD models have appeared in the literature [3,20].In [3], the authors study the discrete and discrete fractional representation of a PK/PD model describing tumor growth and anti-cancer effects in continuous time considering a specific times scale, while in [20] a fractional PK/PD model in anesthesia is developed to describe the nonlinear characteristics of the PK/PD patient models.
In this article, we propose a novel fractional four-compartmental PK/PD model for the induction phase of anesthesia, employing the ψ-Caputo fractional derivative.This innovative approach involves replacing each ordinary derivative in the classical PK/PD model investigated in [41,43] with the ψ-Caputo fractional derivative of order α ∈ (0, 1] [1].
By incorporating fractional derivatives into the model, we aim to capture the intricate dynamics of anesthesia induction more accurately.Furthermore, based on the Picard iterative process, we establish a solution for the nonhomogeneous ψ-Caputo fractional differential system of equations in terms of matrix Mittag-Leffler functions.The solution provides a deeper understanding of the dynamics and behavior of the anesthesia process, taking into account the complex interplay between drug concentration and its effect on the body.
To validate the effectiveness of our proposed model, we conduct numerical simulations.Specifically, we determine the optimal anesthesia dosage for a 53-year-old male weighing 77 Kg and measuring 177 cm, while considering the minimum treatment time [41].These simulations enable us to evaluate the performance and applicability of the fractional PK/PD model in a practical context.
Through the article, we aim to demonstrate the potential of fractional derivatives in advancing our understanding of biology, particularly in the field of PK/PD modeling.By incorporating fractional calculus into these models, we can unlock new insights and improve our ability to predict and optimize drug behavior within biological systems.
The paper is organized as follows.In Section 2, we provide a review of several definitions and properties of fractional calculus that are essential for our subsequent discussions (Section 2.1); and we recall the Pharmacokinetic and Pharmacodynamic models proposed by Bailey et al. [5] and Schnider [36] (Section 2.2), discussing their bispectral index (Section 2.3) and the equilibrium points (Section 2.4).Our original contributions are then given in Section 3: we obtain a solution to a general linear nonhomogeneous ψ-Caputo fractional system (Section 3.1); we introduce a novel PK/PD model for the induction phase of anesthesia based on the ψ-Caputo fractional derivative (Section 3.2); and finally we compute the model parameters using the Schnider model [36], presenting the numerical results of the fractional PK/PD model corresponding to different ψ functions and fractional orders (Section 3.3).We conclude with Section 4, discussing the implications and limitations of our results, and with Section 5, summarizing our findings and outlining potential directions for future research.

Preliminaries
In this section, we recall several definitions, properties of fractional calculus, and the classical PK/PD model, that will be used in the sequel.

Fundamental definitions and results
Throughout the paper, ψ designates a function of class Definition 1 (See [26]).The left ψ-Riemann-Liouville fractional integral of a function f of order α ∈ (0, 1) is defined by where Γ(•) is the Euler Gamma function.
Definition 2 (See [26]).The ψ-Caputo fractional derivative of a function f of order α ∈ (0, 1) can be defined as follows: We have the following properties of the fractional operators with respect to function ψ.
Definition 3 (See [17]).The Mittag-Leffler function of one parameter, of a matrix A, is defined as , Re(α) > 0. ( Definition 4 (See [17]).The Mittag-Leffler function of two parameters, of a matrix A, is defined as Remark 2. The matrix exponential function is a special case of the matrix Mittag-Leffler function [24].For We now recall the notion of generalized convolution integral.
Definition 5 (See [23]).Let f and g be two functions which are piecewise continuous at any interval [a, b] and of exponential order.The generalized convolution of f and g is defined by

The classical PK/PD model: state of the art
The Pharmacokinetic/Pharmacodynamic (PK/PD) model comprises four compartments: intramuscular blood (y 1 (t)), muscle (y 2 (t)), fat (y 3 (t)), and the effect site (y 4 (t)).The inclusion of the effect site compartment (representing the brain) is necessary to account for the finite equilibration time between the central compartment and the central nervous system concentrations [5].This model is employed to characterize the distribution of drugs within a patient's body and can be mathematically described by a four-dimensional dynamical system: subject to the initial conditions where y 1 (t), y 2 (t), y 3 (t) and y 4 (t) represent, respectively, the masses of the propofol in the compartments of blood, muscle, fat, and effect site at time t.The control u(t) represents the continuous infusion rate of the anesthetic.The parameters a 1 0 and a e 0 represent, respectively, the rate of clearance from the central compartment and the effect site.The parameters a 1 2 , a 1 3 , a 2 1 , a 3 1 and a e 0 /v 1 are the transfer rates of the drug between compartments.All these parameters depend on age, weight, height and gender, and the relations can be found in Table 1.
A schematic diagram of the dynamical control system (3) is given in Figure 1.Following Schnider et al. [36], the lean body mass (LBM) is calculated using the James formula, which performs satisfactorily in normal and moderately obese patients, but not so well for severely obese cases [22].The James formula calculates LBM as follows: for Female, LBM =

The bispectral index (BIS)
The bispectral index (BIS) serves as an indicator of anesthesia depth, obtained by analyzing the electroencephalogram (EEG) signal and reflecting the effect site concentration of y 4 (t).It provides a quantitative measure of a patient's level of consciousness, ranging from 0 (indicating no cerebral activity) to 100 (representing a fully awake patient).According to clinical guidelines, maintaining BIS values within the range of 40 to 60 is considered essential to ensure adequate general anesthesia during surgical procedures [4].BIS values below 40 indicate a deep hypnotic state, while values above 60 may increase the risk of awareness under anesthesia.Thus, it is crucial to closely monitor and regulate the BIS value within the optimal range of 40 to 60 to prevent unintended consciousness during anesthesia and ensure patient safety [14].
Empirically, the BIS can be described by a decreasing sigmoid function, as outlined by Bailey et al. [5]: where the parameters in the BIS model have specific interpretations.Function BIS 0 gives the BIS value of an awake patient that is typically set to 100; EC 50 corresponds to the drug concentration at which 50% of the maximum effect is achieved, while γ is a parameter that captures the degree of nonlinearity in the model.According to Haddad et al. [18], typical values for these parameters are EC 50 = 3.4 mg/l, and γ = 3.

The equilibrium point
The equilibrium points are computed by putting the right-hand side of the four equations given in (3) equal to zero with the condition It results that the equilibrium point y e = (y e 1 , y e 2 , y e 3 , y e 4 ) is given by and the value of the continuous infusion rate for this equilibrium is The fast state is defined by For more information on the classical PK/PD model we refer the interested reader to [42].

Main Results
We begin by using the Picard iterative process to prove a series solution to a linear nonhomogeneous ψ-Caputo fractional system: see Theorem 2, in Section 3.1.Then, we generalize the state-of-theart PK/PD model (3) by introducing in Section 3.2 a more general ψ-Caputo fractional PK/PD model that is covered by our Theorem 2. We finish our new results in Section 3.3, by investigating numerically the new fractional model and comparing the efficacy of function ψ.

Solution of linear nonhomogeneous ψ-Caputo fractional systems
Consider the following linear nonhomogeneous fractional equation: subject to the initial condition y(a) = y 0 , where Proof.Follows by using the change of variable z = ψ −1 (ψ(t) + ψ(a) − ψ(s)), Definition 5, and performing direct calculations.
Lemma 3. Let α ∈ (0, 1] and C be a constant.Then, one has Proof.From Definition 1, we have and the proof is complete.Now, we shall utilize the Picard iterative process [12] to formulate a series solution to ( 12)-( 13).
Theorem 2. The solution of the initial value problem (12)-( 13) can be given in series form as Proof.Applying the fractional integration operator I α,ψ a to both sides of equation ( 12), and using Theorem 1, we obtain the following expression: y(t) = y(a) + AI α,ψ a y(t) + I α,ψ a u(t).Let ϕ k be the kth approximate solution with the initial one given by ϕ 0 (a) = y(a) and, for k ≥ 1, the recurrent formula being satisfied.From formula (15) and Lemma 3, one has By virtue of Lemma 2 and by taking the limit k −→ ∞ for ϕ k (•), we obtain the series formula (14) for the solution of ( 12)-( 13).
Note that, in terms of the matrix Mittag-Leffler functions (1) and ( 2), the solution ( 14) may be written as

A fractional PK/PD model
Motivated by system (3), we introduce here our ψ-Caputo fractional Pharmacokinetic/Pharmacodynamic model, which is obtained by replacing each ordinary derivative in the system (3) by the ψ-Caputo fractional derivative of order α ∈ (0, 1].Then, our proposed PK/PD model can be expressed by the following four-dimensional fractional dynamical system: subject to the initial conditions According to the dynamical system (12), one may write system ( 17)-( 18) in a matrix form as follows: with One mentions that the continuous infusion rate u 1 (t) is to be chosen in such a way to transfer the system (17) from the initial state (wake state) to the fast final state (anesthetized state).

Numerical simulations
To administer anesthesia to a 53-year-old man weighing 77 Kg and measuring 177 cm, we utilize our proposed fractional PK/PD system described by: where, according with Table 1 and [41], the matrix A is taken as with From Theorem 2 of Section 3.1, written in form ( 16), the solution of system ( 20) is given by with u(t) = Bu 1 (t) = [u 1 (t), 0, 0, 0] T .Figure 2 showcases the solutions derived from the fractional PK/PD model (20), considering the function ψ(t) = t and exploring different fractional order values: α = 1, α = 0.95, α = 0.9, and α = 0.85.In Figure 3, the curves represent the controlled BIS (Bispectral Index) associated with the optimal continuous infusion rate of the administered anesthetic u(t).It is noteworthy that when the function ψ(t) = t and the fractional order is set to α = 1, then the obtained results resemblance those derived from the classical PK/PD model (3).However, altering the fractional orders introduces variations in the degree of anesthesia.The recorded values for all fractional orders fell within the range of 40 to 50 (corresponding to the classical model), thus ensuring the condition of anesthesia.Nevertheless, it is crucial to acknowledge that lower fractional order values entail a higher risk of awareness during anesthesia.
, and ψ(t) = t + 0.2, when considering a fractional order of α = 1.The graphs shown in Figure 5 depict the controlled BIS corresponding to a specific value of the fractional order α = 1, under functions ψ(t) = t, ψ(t) = √ t, ψ(t) = t 2 and ψ(t) = t + 0.2.It is observed that selecting functions ψ(t) = √ t and ψ(t) = t 2 does not yield satisfactory anesthesia results.On the other hand, employing the functions ψ(t) = t and ψ(t) = t + 0.2 leads to favorable anesthesia outcomes.In subsequent simulations, we will maintain the functions ψ(t) = t and ψ(t) = t + 0.2 while altering the fractional orders.

Discussion
The ψ-Caputo fractional PK/PD model represents a notable advancement in modeling due to its ability to capture intricate drug response dynamics, considering the complex relationship between drug concentration and physiological effects.In comparison to conventional modeling approaches,  the incorporation of the ψ-Caputo fractional derivative and the fractional order α enable a more comprehensive representation of non-local and memory-dependent effects, offering a more nuanced understanding of the induction phase of anesthesia.However, it should be noted that the ψ-Caputo fractional PK/PD model is more complex compared to conventional modeling approaches [41,43] because it relies on fractional calculus and an unknown ψ function, introducing additional parameters, and requires specialized mathematical knowledge.Therefore, its use should be reserved for situations where traditional models are insufficient for capturing the intricacies of the system being modeled with conventional models, based on integer-order calculus, remaining the standard for most PK/PD modeling tasks due to their ease of use and the availability of established methodologies and software [42].
Complex systems often require sophisticated modeling techniques that consider intricate interactions mechanisms.Additionally, the availability of high-quality data, incorporating a wide range of variables, significantly enhances the reliability of model predictions.Given this, the accuracy and predictability of ψ-Caputo fractional PK/PD models are contingent upon the complexity of the research domain and the quality of available data, presenting some limitations: (i) additional parameters; and (ii) data requirements.Indeed, (i) the ψ-Caputo model introduces extra parameters associated with the fractional orders, which need to be estimated from data.These additional parameters can make the model more complex and require more data for accurate estimation.(ii) Fractional PK/PD models may require more extensive data to accurately estimate parameters and capture the complexities of the system.Conventional models, in some cases, might work with fewer data points.Moreover, the ψ-Caputo fractional PK/PD model inherently carries uncertainties related to the selection of the function ψ and the fractional order α, which can impact the model's reliability in predicting anesthesia dosage.Uncertainties in this model can arise from various sources, e.g.(i) measurement uncertainty (errors in drug concentration measurements can propagate into the model, affecting the accuracy of estimated parameters and predictions), (ii) inter-individual variability (patients can have varying physiological characteristics, which can lead to inter-individual variability in drug response, with the model not capturing all aspects of this variability, potentially resulting in suboptimal dosing for some patients), (iii) intra-individual variability (even within the same individual, the response to anesthesia drugs can vary due to factors like changes in organ function, health status, or concurrent medications, and in that case the ψ-Caputo fractional model may not account for all these factors, leading to uncertainty in dosing recommendations for the same patient at different times), (iv) parameter estimation uncertainty (the model relies on various parameters such as clearance, volume of distribution, and rate constants, and estimating these experimental data can introduce uncertainty), (v) data uncertainty (the quality and quantity of data used to estimate model parameters can vary, and this can lead to uncertainties in the model; sparse or noisy data can result in less reliable model predictions, affecting the accuracy of anesthesia dosage recommendations), (vi) model structure uncertainty (the ψ-Caputo fractional PK/PD model itself is a simplification of the complex processes occurring in the body, not fully capturing all relevant mechanisms, leading to uncertainties in predictions).Consequences of these uncertainties in the case of anesthesia (dosage) can be significant: underdosing, if the model underestimates the required dosage, in which case patients may not achieve the desired level of anesthesia, leading to inadequate pain control or awareness during surgery; overdosing, if the model overestimates the required dosage, with patients receiving an excessive anesthesia, leading to complications such as extended recovery times, respiratory depression, or other side effects.
While the ψ-Caputo fractional PK/PD model shows advancements in capturing complex drug responses, there might be certain types of drugs or treatments for which this modeling approach may not be as suitable.Some reasons why fractional PK/PD modeling may not be ideal in certain cases include: (i) non-linear pharmacokinetics; (ii) complex drug interactions; (iii) short-acting anesthetics.
Finally, it should be remarked that patient safety must be always of primary importance, especially concerning regulatory aspects like the Bispectral Index (BIS).Any modeling approach should promote patient safety by accurately predicting drug responses and aiding in the development of safer and more effective treatment protocols.Some key aspects on patient safety and the regulation of PKPD models, with a focus on BIS monitoring, are: • Balancing Technology and Clinical Expertise.While PKPD models and BIS monitoring are valuable tools in anesthesia, they should not replace the judgment and expertise of skilled anesthesiologists.The safe administration of anesthesia requires a combination of technology and clinical acumen.
• Evidence-Based Practice.The development and use of PKPD models and monitoring devices like BIS should be based on rigorous scientific evidence and clinical studies.Regulatory bodies should require strong evidence of safety and efficacy before approving or endorsing these technologies.
• Regulatory Oversight.Regulatory agencies play a crucial role in ensuring that medical devices and models are safe and effective.They should establish and enforce guidelines for the development and use of PKPD models and monitoring devices, including BIS.Regular assessments and updates should be conducted to account for evolving scientific knowledge.
• Continuous Monitoring.Real-time monitoring of patients during surgery is essential for patient safety.Monitoring parameters, including BIS values, should be continuously observed and interpreted by trained personnel to respond promptly to any adverse events or deviations.
• Individualized Care.Each patient is unique, and anesthesia management should be tailored to their specific needs and responses.PKPD models can provide guidance, but individualized care remains central to patient safety.

Conclusion
The incorporation of the ψ-Caputo fractional derivative in Pharmacokinetics and Pharmacodynamics modeling represents a significant advancement in the field.Indeed, by utilizing fractional-order derivatives, researchers can more accurately capture the complex and non-local behavior observed in drugs within biological systems.The choice of the function ψ and the fractional order α holds critical importance in modeling the relationship between drug concentrations and pharmacological effects.This approach provides a more realistic representation of drug efficacy and dose-response relationships, allowing for a deeper understanding of the intricate dynamics involved in drug-target interactions.
While our study successfully demonstrates the potential of the ψ-Caputo fractional PK/PD model in capturing complex drug responses during the induction phase of anesthesia, it is crucial to acknowledge certain limitations.The current model's applicability may be constrained by the specific parameters chosen and the simplifications employed to facilitate numerical analysis.Furthermore, the availability of comprehensive clinical data for model validation remains a challenge, which could affect the generalizability of the findings.Future research efforts should focus on addressing these limitations to further enhance the applicability and robustness of the proposed model.Moreover, further research is necessary to explore the impact of the chosen function ψ and the fractional order α on time-delayed responses.This area remains open for investigation, and future studies can delve into understanding how different choices of ψ and α influence the temporal aspects of drug responses.
In summary, the incorporation of ψ-Caputo fractional derivatives in Pharmacokinetics and Pharmacodynamics modeling offers valuable insights and advancements.By refining the choice of function ψ and fractional order α, researchers can enhance the accuracy and realism of drug modeling, paving the way for a more comprehensive understanding of drug behavior in biological systems.

Figure 1 :
Figure 1: Schematic diagram of the PK/PD model with the effect site compartment of Bailey and Haddad [5].

Figure 4
Figure 4 illustrates the solutions of the fractional PK/PD model (20) associated with the functions ψ(t) = t, ψ(t) = √ t, ψ(t) = t 2, and ψ(t) = t + 0.2, when considering a fractional order of α = 1.The graphs shown in Figure5depict the controlled BIS corresponding to a specific value of the fractional order α = 1, under functions ψ(t) = t, ψ(t) = √ t, ψ(t) = t 2 and ψ(t) = t + 0.2.It is observed that selecting functions ψ(t) = √ t and ψ(t) = t 2 does not yield satisfactory anesthesia results.On the other hand, employing the functions ψ(t) = t and ψ(t) = t + 0.2 leads to favorable anesthesia outcomes.In subsequent simulations, we will maintain the functions ψ(t) = t and ψ(t) = t + 0.2 while altering the fractional orders.In Figure6, we present the solutions of the fractional PK/PD model(20) corresponding to the functions ψ(t) = t and ψ(t) = t + 0.2, under the fractional orders α = 1, α = 0.9, and α = 0.8.The curves representing the controlled BIS are displayed in Figure7.It is worth noting that the
. It is worth noting that the