Modeling Drug Concentration Level in Blood Using Fractional Differential Equation Based on Psi-Caputo Derivative

This article studies a pharmacokinetics problem, which is the mathematical modeling of a drug concentration variation in human blood, starting from the injection time. Theories and applications of fractional calculus are the main tools through which we establish main results. The psi-Caputo fractional derivative plays a substantial role in the study. We prove the existence and uniqueness of the solution to the problem using the psi-Caputo fractional derivative. The application of the theoretical results on two data sets shows the following results. For the ﬁ rst data set, a psi-Caputo with the kernel ψ = x + 1 is the best approach as it yields a mean square error (MSE) of 0 : 04065 . The second best is the simple fractional method whose MSE is 0 : 05814 ; ﬁ nally, the classical approach is in the third position with an MSE of 0 : 07299 . For the second data set, a psi-Caputo with the kernel ψ = x + 1 is the best approach as it yields an MSE of 0 : 03482 . The second best is the simple fractional method whose MSE is 0 : 04116 and, ﬁ nally, the classical approach with an MSE of 0 : 048640 .


Introduction
To treat an infection from a human being or even from an animal, a suitable dose of medicine is substantial. Owing to the amount of the drug in the blood plasma decreasing with time, medicine must be given in multiple doses.
In phase I of clinical development, the time to achievement of steady state of a new regimen is routinely evaluated. The time to the achievability of the steady state is the time needed until the drug concentration is stable in the blood, i.e., does not display an increasing tendency by drug accumulation. If a drug is given at orderly dosing intervals, drug sediment from preceding doses is accumulated. Stabilization of the concentration occurs when the quantity of drug discarded during the dosing interval equals the amount that was given. In order to evaluate the time to achieve steady state, blood is sampled at a certain time point within each dosing interval see [1].
Following the drug concentrations at these time points [2] Jordan et al. proposed a model for the achievement of the steady state of the drug concentration in plasma. Authors of [3] proposed a model for the prediction of the "unbound brain-toplasma" drug concentration ratio. Zhang et al. [4] studied the ratio of the drug concentration between tissues and plasma.
A simple and novel sensor was developed for the analysis of clinical doxorubicin (DOX) concentration based on the screen-printed electrode by evaluating the DOX concentration, see [5]. In [6] the authors evaluated the drug concentrations in postmortem blood samples where the value of concentrations slightly differs depending on the sample site. For more works on drug concentration in the blood, we refer the reader to the references [7][8][9][10][11].
Fractional calculus (F.C) helps to describe models and natural phenomena problems. Many researches have dedicated their work in this branch (see, e.g., previous studies [12][13][14][15][16][17][18][19]). The results obtained were significantly positive in different fields of Medicine and Biology. The foundation of fractional calculus is laid on fractional integrals and derivatives. The efficiency of the fractional order model over the integer order is investigated by Bagley and Torvik [20]. Through fractional calculus, Djordjevic et al. [21] developed a rheological model of airway smooth muscle cells, which came as an alternative to the least square data fitting technique often used for this purpose. Recently, an application of fractional calculus to nanotechnology was proposed in [22]. These are just a few out of plenty examples of research works in which fractional calculus has proven its efficiency compared to existing classical approach.
This research article contributes to showing the power of mathematical modeling using fractional calculus. In particular, a class of fractional derivatives called ψ-Caputo, introduced by Almeida [23], has proven its efficiency in various applications including a recent study by Awadalla et al. [24]. A preliminary investigation on the topic of this study was carried out by us. The results of the said investigation are provided in this reference [25]. The main contribution of this work to the literature is the reduction or further minimization of the MSE in modeling the drug's concentration kinetics. The article starts with an introduction; then, some preliminaries of fractional calculus are covered. Elements of pharmacokinetics and the mathematical model of drug concentration in the blood are discussed in the third section. Main theoretical results are established in the fourth section followed by application examples in the fifth section. Finally, the last section provides concluding remarks on the overall study.
More generally, realization of this work was motivated by the aim to reduce modeling error in pharmacokinetics.

Preliminaries
This section is dedicated to preliminary tools that will enable a smooth study in upcoming sections. Indeed, F.C theories are built thanks to theorems, definitions, and lemmas. Similar to classical calculus, integrals and derivatives taken in the fractional sense are the foundation of F.C. Concerning the abovementioned sentences, selected definitions, theorems, and notations are discussed in this section. They play important roles in the rest of the paper.
Definition 1 (see [9]). The fractional integral of order ξ > 0 of a function H : ½a, b ⟶ R e taken in the Riemann-Liouville sense is defined as Equation (1) holds only if the right-hand side of the given integral is point-wise defined on 0, +∞½. Note that the function Γ is the commonly-known gamma function, which is defined as follows: Definition 2 (see [9]). The Caputo derivative of order ξ > 0 of a function H : ½a, b ⟶ R e is defined as where n = bξc + 1, bξc is the integer part of ξ.
Definition 3 (see [2]). Let ξ > 0, H ∈ L 1 ½a, b and ψ ∈ C 1 ½a , b is selected to be an increasing function with ψ ′ ðxÞ ≠ 0, ∀x ∈ ½a, b; then, the notation I ξ,ψ 0 + H ðϱÞ represents the fractional integral of H with respect to another function ψ, and it is defined by The idea of integrating a function definition is known as ψ − Caputo integral of the fractional integral is taken in the Caputo sense. This is used in the sequel for formulating the solution to the ψ − Caputo fractional model. Definition 4 (see [2]). Let ξ > 0 and H , ψ ∈ C n ½a, b where ψ is an increasing function and ψ ′ ðxÞ ≠ 0, ∀x ∈ ½a, b. Then, ð C D ξ,ψ 0 + H ÞðϱÞ denotes the fractional derivative of H with respect to ψ. ψ − Caputo if the fractional derivative is taken in the Caputo sense, and it is given by Lemma 5 (see [2]). Let ξ > 0 and n be a positive integer such that ξ ∈ n − 1, n½. For every H , ψ ∈ C n ½a, b, we have

Pharmacokinetics and Drug Concentration Model
A brief definition and overview of pharmacokinetics are proposed in this section. To treat an infection from a human body, a suitable dose of medicine is substantial. Once a drug is administrated to an individual through intravenous injection, it has an initial concentration that decreases over time.
The decrement appears as a result of metabolism and excretion. Pharmacokinetics is a branch of medicine that studies the dynamic (kinetics of drugs in a living body). Owing to the amount of the drug in the human body decreasing with time medicine must be given in multiple doses. Two main streams of study exist in pharmacokinetics, the theoretical approach and the experimental approach.

Journal of Mathematics
The former approach focuses on the development of a pharmacokinetics (mathematical) model that predicts drug concentration levels in the blood over time. The latter method involves empirical development based on a biological sample, during which analytical methods for drugs and their metabolites are measured. In this case, it is required to have an adequate experimental setting for data collection and handling. This article focuses in the sequel on the development of a mathematical pharmacokinetics model. The entire process of absorption, distribution, metabolism, and elimination (ADME) of a drug is illustrated in the next figure as described in [26,27]. Figure 1 represents the ADME process of a drug after it has been administrated to a human. The process is governed by a change in concentration over time. The change rate can be denoted by ±dðconcentrationÞ/dϱ. More generally, let us denote by Y the drug concentration in the body; then, the mathematical model describing the rate change is given by where k is a constant to be experimentally determined for each drug. If a patient is given an initial drug dose, Y 0 , then, the drug level in his body at any time is the solution of the differential equation defined by Eq. (6), that is, The objective of this work is to prove the advantages of modeling drug concentration kinetics using fractional differential equations. Indeed, we will show empirically that modeling results using F.C are better than those obtained from classical calculus. The starting point is to build a fractional counterpart of Eq. (6). The said equation is built as follows or simply The Lemma below is defined to introduce a general form of the solution to Eq. (9) representing drug concentration kinetics with initial condition. Lemma 6. Let us consider Q, with the assumption that it is an integrable function, which is defined over ½0, T . It follows that the solution of the fractional differential equation given by Eq. (9) has a general form which can be expressed by the following integral equation It is worth to highlight in Eq. (10) from Lemma 6 the presence of ψ-Caputo fractional integral.
In the application, Kernel functions are selected under the data distribution. There is not a steady rule for that. However, a linear kernel works well in many cases. Other commonly used kernel includes but is not limited to ffiffiffi u p , log u.

Main Result of Psi-Caputo Drug Concentration Model
This section aims to investigate theoretical study around Eq. (9). The final goal is to build prove of the existence and the uniqueness of a solution to Eq. (9). Let C½0, T be the space of real valued functions that are continuous on ½0, T endowed with the norm of the uniform convergence: kYk ∞ = sup ϱ∈½0,T jYðϱÞj for every Y ∈ C½0, T . Then, Φ ≔ ðC½0, T , k•k ∞ Þ is a Banach space.
An operator T : Φ ⟶ Φ defined and attached to the problem introduced by Eq. (9) can be built as In what follows, the existence of a solution to Eq. (10) is proved followed by a proof of the uniqueness of the said solution. Before establishing the proof of the main results, let us first establish the following statements. Indeed, the statements are mathematical hypotheses that are used in sections dedicated to proofs.

Journal of Mathematics
Existence of at least one solution to the problem Eq. (9) is proved in the theorem below Theorem 7. Let us assume that the three conditions ðA 1 Þ, ð A 3 Þ, and ðA 4 Þ hold. Then, there exists at least one solution to Eq. (9). The said solution is in the interval ½0, T .
Proof. The proof of heorem 7 will be divided into several steps. The first step consists of showing that the operator Applying the supremum on t on both sides of Eq. (17) leads to The second step of this proof is to show that the operator Let ϱ 1 , ϱ 2 ∈ ½0, T with ϱ 1 < ϱ 2 and Y ∈ B λ . The relation below holds as results of the assumptions ðA1Þ-ðA4Þ The right-hand side of Eq. (19) tends to zero as ϱ 1 ⟶ ϱ 2 . That is jTYðϱ 2 Þ − TYðϱ 1 Þj ⟶ 0 as ϱ 1 ⟶ ϱ 2 . It is worth observing that the right hand part of Eq. (19) does not depend on Y ∈ B λ , this implies by Arzela-Ascoli theorem that TðB λ Þ is completely continuous.
The third step of the proof requires a last intermediate step to complete the assumptions of the Leray-Schauder nonlinear alternative theorem. This consists of showing the boundedness of the set of all solutions to equation Y = δT ðYÞ. Assume that Y is a solution Eq. (9), then, it follows from Eq. (10) that Inverting both sides of Eq. (20) and dividing them by R.H.S of Eq. (20) leads to the following relation   Journal of Mathematics Recalling ðA 4 Þ, there exists a constant W > 0, which is indeed such that W ≠ Y. Moreover, let us construct the set Ω = fY ∈ Φ ; kYk ∞ < Wg. It is obvious that the operator T : Ω ⟶ Φ is continuous and completely continuous. Based on the constructed Ω, there exists Y ∈ ∂Ω such that YδTðYÞ for some δ ∈ 0, 1½. Consequently, by the nonlinear alternative of Leray-Schauder type, we deduce that T has a fixed point Y ∈ Ω which is a solution to the problem defined by Eq. (10). Theorem 8. Let us assume that conditions ðA 1 Þ and ðA 2 Þ hold. Moreover, if the condition holds; then, the problem defined by Eq. (9) has a solution which is unique. The said solution belongs to the interval ½0, T .
Proof. Let us consider the operator T defined in Eq. (12). Let us also define a ball where M Q ≔ sup 0≤ϱ≤T jQðϱ, 0Þj.
First, let us show that TB ε ⊂ B ε . For any Y ∈ B ε , ρ ∈ ½0, T , and using Eq. (12), we have the following relation On the other hand, Substituting the appropriate fragment of equation in Eq. (24) by Eq. (25) implies a new relation which is given by Equation (26) is sufficient to conclude that TB ε ⊂ B ε . The second step of the proof is to show that the considered operator is a contraction mapping. For every Y 1 , Y 2 ∈ Φ, the following relation holds.
From Eq. (27), we deduce that T is a contraction. By the Banach contraction mapping theorem, the problem defined by Eq. (9) has a unique solution. The said solution belongs to the interval ½0, T .

Application Example
In this section, application examples are provided to support the theoretical work developed above. Data set obtained from real-life experiment were used. The classical method, simple fractional method, and kernel fractional method were used to fit the data set. Finally, a comparison of results is done to support theoretical findings.

First Experimental Data
Set. This data set was retrieved from [28]. The author carried out an experiment in which he measured a drug concentration in (mg/L) over 6 hours of an antibiotic. A single dose of the said antibiotic was administered intravenously to a 50-kilogram woman. The dose level was 20 mg/kg. A scatter plot of the concentration data over time shows a decay. Three deterministic approaches were used to fit the data set. The first approach is what is referred to as the classical approach, in which the general solution is defined by Eq. (7). The second and third approaches are fractional differential method and kernel fractional differential method, respectively, which general solutions are given by Eq. (11), respectively, with ψðxÞ = x, trivial kernel, and ψðxÞ ≠ x, pure kernel. The selected pure kernel here is linear ψðxÞ = x + 1. It was observed empirically that using any linear kernel ψðxÞ = x + a, with a ≠ 0, would produce a similar result to the case where ψðxÞ = x + 1 is used. Table 1 summarises one hand best estimates of the parameters for both classical and fractional approaches. On the other hand, it displays the MSE of each method. It is observed that the fractional kernel method with ψðxÞ = x + 1 performed the best, followed by the fractional method with ψðxÞ = x and lastly the classical method. It is worth highlighting consistency in the results. Indeed, the solution to the classical approach is a first order differential equation; therefore, one would expect the solution to its fractional counterparts to be such that ξ ∈ ð0, 1Þ ∪ ð1, 2Þ. Figure 2 is the graph of the original data set, the fitted line is obtained from the classical model, and the fitted line is obtained from the fractional method with ψðxÞ = x. Both fitted lines overlap over each other at the beginning of the plot but show a difference toward the end. The fractional 5 Journal of Mathematics method with ψðxÞ = x seems to be closer to the true data. MSE in Table 1 is an evidence. Figure 3 is the graph of the original data set, the fitted line is obtained from the classical model, and the fitted line is obtained from the fractional method with ψðxÞ = x + 1. Similar to Figure 2, a close look at the figure reveals that the fractional method with ψðxÞ = x + 1 does the job better than the classical method. Moreover, recalling Figures 2 and 3 as well as Table 1, the ordinal classification (first : fractional method with ψðxÞ = x + 1; second : fractional method with ψðxÞ = x; and third : classical method) of the three methods used in this work becomes explicit.

Second Experimental Data Set.
In this experiment, a newly developed drug was administrated to a patient. The administration was done through an IV injection. A sample of blood was taken regularly, and the drug plasma concentration was determined. The data set was retrieved from [29].
Similar experimental procedures to those from the first example are used. The best values of parameters as well as MSE of each method are consigned in the table below. Table 2 summarises and displays the results of the second experiment. It is observed that the fractional kernel method with ψðxÞ = x + 1 performed the best, followed by the fractional method with ψðxÞ = x and lastly the classical method. Moreover, the fractional order of derivatives is always such that ξ ∈ ð0, 1Þ ∪ ð1, 2Þ, proving that the results are in line with the first-order differential equation. Hence, the results are consistent. Figure 4 displays the original data, the classical solution, and the fractional with the kernel ψðxÞ = x. The results in Table 2 are reflected by the fact that the fractional approach fits the original data points better than the classical approach does.     Journal of Mathematics Figure 5 is the graph of the original data set, the fitted line is obtained from the classical model, and the fitted line is obtained from the fractional method with ψðxÞ = x + 1. A close check of the figure reveals that the fractional method with ψðxÞ = x + 1 does the job better than the classical method.

Conclusion and Future Work
In this work, we studied the ψ-Caputo type of fractional differential equation. This derivative is the fractional analog of the so-called ðfogÞ′ derivative in classical calculus. The existence and uniqueness of the proposed method were discussed before the application examples. Experiment results show that the ψ-Caputo method which uses a pure kernel function performed the best, followed by a simple fractional approach and finally the classical method. The fractional order of derivative that allows to best fit the data is always such that ξ ∈ ð0, 1Þ ∪ ð1, 2Þ, which is in line with the setup of the studied problem since the classical approach solution is a first-order differential equation. The experimental section has revealed the following results: For the first data set, a psi-Caputo with the kernel ψ = x + 1 is the best approach as it yields a mean square error (MSE) of 0:04065. The second best is the simple fractional method whose MSE is 0:05814; finally, the classical approach is in the third position with an MSE of 0:07299.
For the second data set, a psi-Caputo with the kernel ψ = x + 1 is the best approach as it yields an MSE of 0:03482 . The second best is the simple fractional method whose MSE is 0:04116 and, finally, the classical approach with an MSE of 0:048640.
In future works, we aim to investigate if the obtained results hold for all or most of existing fractional derivatives. We may also study properties that can help in the selection of a suitable kernel function. In the current study, the selection ψ function was done randomly, on a try and error basis, until we found out that a family of linear functions could well do the job.

Data Availability
The data set used in application section is available through the url provided in [27].