One-compartment stochastic pharmacokinetic model

In this work, we consider a pharmacokinetic (PK) model with first-order drug absorption and first-order elimination that represent the concentration of drugs in the body, including both the absorption and elimination parts, and we also add a random factor to describe the variability between patients and the environment. Using Itô’s lemma and the Laplace transform, we obtain the solutions in integral form for a single and constant dosage regimen in time. Moreover, formulas for the expected value and the variance for each case of study are presented, which allows the statistical assessment of the proposed models, as well as predicting the ideal path of drug concentration and its uncertainty. These results are important in the long-term analysis of drug concentration and the persistence of therapeutic level. Further, a numerical method for the solution of the stochastic differential equation (SDE) is introduced and developed.


Introduction
Drug concentration levels in the body vary among different patients according to their weight, age, stress, or genetic factors [1]. In addition, it has also been observed that the food intake, special excercises and some vitamin can affect the drug absorption [2]. Due to the variability and uncertainty of such factors, several researches have introduced stochastic corrections on the well known compartmental pharmacokinetics deterministic models [3].
The mathematical modeling of stochastic compartmental systems has received great attention in the literature and has produced many useful mathematical models. For example, an analysis of a one-compartment model to describe the spontaneous erratic variations of drug concentration decay with a single dose administration [4,5], where the elimination rate fluctates according a Wiener motion. More recently the same case is studied in [6] where stochastic differential equations (SDE) driven by the Liu process (see [7]) are proposed and this model presents the uncertainty that the drug concentration is larger than the minimum effective concentration. In contrast to [6], in this paper we include the drug absorption phenomena under two kinds of dosage forms and the elimination using a Wiener motion.
Another approach is presented in [4] and [8] using a Vasicek model to study the stochastic variability for continuous dosing. Some recent studies have proposed more complex compartments models [9,10,11,12,13,14]. These works consider the noise term as a constant and focus on parameter estimation techniques. However, in such works closed expressions for the solution, the expected value, and the time-dependent variance are not derivated. The aim of this article is to determine such formulas.
In addition, the problem to determine an optimal dosing timing schedule to control the blood drug concentration is studied in [2]. The authors combine optimal control theory and stochastic methods, prove existence, uniqueness of the solution and present the corresponding stability analisys. In such work the goal is to determine the source term of the SDEs system. On the other hand, in this paper we are interested in determining the terapeutic range of the drug when dosification function is known.
In this article, we consider an one-compartment PK model where drug is absorbed according to a first-order process and first order elimination. Such situation is described by a differential equation coupled with a SDE, that takes into account the variability among individuals under dose regimen. In order to solve such system, we employ Itô's calculus and Laplace transformation to obtain a solution in Riemann integral form. The advantage is that one can approximate numerically the realistic trajectory of the drug concentration by standard quadrature rules. To determine the therapeutic range of the drug, it is necessary to know the expected value and the time-dependent variance of the drug concentration. In order to find them, it is possible to solve the differential equations for the first moment (i.e. the expected value) and for the second moment (and so find the variance), this is usefull when an exact solution of the corresponding SDE is not known. However, in this work we obtain an explicit formula for the solution of SDE in an integral form, this allow us to directly calculate the expected value and variance for this solution.

Mathematical model
In order to describe mathematically the absorption and elimination of a drug in the body, we consider the following hypotheses in the deterministic part of the model [3]: a) The absorbed drug concentration decreases proportionally to the amount of drug at site of administration at the time t .
b) The drug is rapidly and uniformly distributed throughout the body, i.e., when the body behaves as a single central compartment.
c) The drug elimination is considered as a first-order process, ie, is proportional to the amount of drug in the body The classical model, which describes drug absorption in a one-compartment and relates the changes in drug concentration in the blood with time to the absorption and the elimination rates is based on the Wagner-Nelson method [15]: where dX= dt is the rate (mg h 1 ) of change of drug concentration in the blood, X.t / is the concentration of drug in the blood or body at time t , X a is the concentration of absorbable drug at the absorption site at time t , r a and r e are the first-order absorption and elimination constants rates, respectively (e.g., h 1 ), r a X a is the first-order rate of absorption (mg h 1 , g h 1 , etc.), r e X.t / is the first-order rate of elimination (e.g., mg h 1 ), X a0 is the initial dose of drug and for the dosage regimen of the drug per unit of time f .t /. We consider the following two cases: • Single dose: f .t / D 0 with X a0 ¤ 0.
• Constant dosage: f .t / D r a X p , with X p positive constant.
However, model (1) does not take into account that drug concentration levels vary among different patients according to their weight, age, stress or genetic factors ( [16] and [1]). Therefore, we consider that the elimination rate r e is not constant in time but randomly fluctuates around a mean value as r e t ; where t is a Gaussian white noise process. Then t dt can be written as c dW .t /; where W .t / is a standard Wiener motion and the positive constant c is a diffusion coefficient.
Incorporating this assumption into the deterministic model (1), we obtain the following SDE system 8 < : Theorem 2.1. The solution to the initial value problems (2) with f .t / D r a X p and X p a positive constant is given by where L 1 is the inverse Laplace transform operator, s is the Laplace transform parameter and X a .t / is the solution of first equation from (2) which is given by the indefinite integral form: Proof.
To determine the solution to problem (2), we use the integration factor method (see [17]) as follows: 1. We solve the non-deterministic part of equation (2), i.e dX.t / D cX.t / dW .t /: Using Itô formula [18] with F .X; t / D ln.X / and after some calculations, we get 2. Applying the method of variation of parameters, we make the constant K of the nondeterministic solution (6) vary as a function of time, so X.t / is of the form where G.t / satisfies equation (5). Then, By substituting (7) into the above equation, we have By cancelling the similar terms and dividing by G.t / in both sides of above equation, we get the following initial value problem where X a .t / is given by (4). To solve the initial value problem (8) we use the Laplace transform, then Now, we apply the inverse Laplace transformation on both sides: We plug equation (9) into (7) to obtain the general solution of (2) where u is an auxiliary variable.

Drug dose regimens
In this section, we analyze the behavior of the drug concentration in the body under constant dosage (infusion) regimen and the single dose (push) administration.

Continuous dosage (infusion) of the drug
A constant amount of the drug is supplied continuously for long periods of time, in the treatment of some chronic diseases. Such types of administration are: intravenous infusion, certain oral formulations based on the phenomenon of osmosis and certain transdermal patches.

Remark 3.1.
When the source function f .t / D r a X p with X p a positive constant, the solution of the first equation of (2) given by (4) is Theorem 3.1. Let f .t / D r a X p be the source function with X p a positive constant that represents the concentration of the drug that is administered at all time t > 0. Then, the concentration of the drug is given by the integral form: where and Proof. By replacing (10) in (3) we get using the convolution theorem, we obtain Finally, we get the integral form where Ä 1 .t u/ and Ä 2 .t u/ are defined as in (12) and (13) respectively.

Remark 3.2.
In conclusion, the solution of the system (2) (12) and (13). (10) and the dose is supplied at the initial time .t D 0/. Then the concentration is given by the integral form Proof. As X p D 0; it follows from (12) and (13) that Thus, immediately from (11) the result is obtained.

Expected value and variance for constant dose
In this section we obtain explicit formulas for the expected value and the variance of the drug concentration X.t /. Having these formulas, we can identify the model parameters from observed data and to establish therapeutic ranges of the drug for each dosing regimen. At this point, it is important to differentiate the case when the absorption constant rate r a is equal to the elimination constant rate r e to determine the relevant pharmacokinetic parameters (see 3.3). In the deterministic model similar results have been obtained in [23] and [24].
Proposition 3.1. The expected value of the stochastic process X.t / given by (11) is if r e D r a : • If r e ¤ r a the expected value of X.t / in (11) is given by D r a r e X p 1 e r e t C r a r a r e X a0 X p e r e t e r a t : • For the case that r e D r a , the result is obtained easily from equation (16).

Proposition 3.2.
If r e ¤ fr a ; r a Cc 2 ; r a Cc 2 =2; c 2 ; c 2 =2g; the variance of the stochastic process X.t / given by (11) is Proof. Details of the calculation of Var X.t / can be seen in Appendix B.

Proposition 3.3.
If r e D r a ; the variance of the stochastic process X.t / given by (11) is Proof. The result can be obtained from the variance definition and equation (23) in Appendix B.
Remark 3.3. The expected value given in (15) satisfies the ordinary differential equation given in (1) when f .t / D r a X p , since if we differentiate (15) with respect to t and substitute into (1), it is satisfied. Furthermore, if X p D 0 in (15), this expression coincides with the exact solution of (1), that is, in the absence of the stochastic term.

Remark 3.4.
To compute the variance of process X.t / in presence of flip-flop kinetics (r e > r a ), we consider the cases r e D r a C c 2 and r e D r a C c 2 =2 (see Appendix B).

Remark 3.5.
In pharmacology is very important to determine the therapeutic range of a drug, which is the range in which the drug can be used without causing toxic or lethal effects on the individual. From equations (15) and (17) we obtain that in the stationary state (t ! 1), the minimum effective concentration X min and the concentration maximum admissible X max must be such that where X.t/ D c r a X p r e p 2r e c 2 when t ! 1.

Numerical simulation
Several numerical schemes have been used to solve SDEs (for example [19] and [20]). However, in this work we propose an alternative method to simulate the evolution of the drug concentration, when a single dose or constant dose are administrated. To compute numerically an approximation of the solution (11), let us denote by˛WD r e C c 2 =2,ˇWD˛ r a , when X p ¤ 0 then h.t / WD r a X p e ˛tCcw.t/ and g.t / WD .X a0 =X p 1/h.t /. Hence the solution (11) can be written as follows Note that the integrals in (18) are just Riemann Integrals, then in order to evaluate X.t /. First, we divide the interval OE0; t max into N sub-intervals of equal length t WD t max =N . This defines a set of discrete times t i D i t; i D 1; : : : ; N . Next, we discretize the Wiener process with a time step t and interpolate linearly the term e cw.u/ on the interval OEt i 1 ; t i . Then, an approximation of X.t k / for k D 1; : : : ; N is where X.t 0 / D X.0/ D 0 and the integrals are computed exactly.

Single dosage (push) administration
We consider the experimental data of Theophylline concentrations (in mg/L) for 12 subjects following a single oral dose of 320 mg. The data is reported in [21] and its time series graphs are shown in Figure 2(a). From this data we identify the parameters r e ; r a ; X a0 and c (see Figure 2(a)) by the method of moments (see [22]). We found that elimination rate, absorption rate, coefficient of variance, and initial absorption concentration are r e D 0:0782, r a D 1:7603, c D 0:1004 and X a0 D 9:5206 respectively. Figure 2(b) shows a simulation of the drug concentration decay. As we can observe, computational simulations are consistent with the experimental measurements. Furthermore, in both figures, the sample paths are in between the strip formed by E X.t / ˙2 X.t/ .

Continuous dosage (infusion) of the drug
We will now study the Propofol concentration behavior during 60 minutes infusion dose administration with an infusion rate of 25 g kg 1 min 1 . Experimental data are taken from [25]. We estimate the parameters of the model by the method of moments, we find that r a D 0:8261; r e D 0:1247; c D 0:1435, X p D 0:0029 and X a0 D 0:0137.
From Figure 3(a) we see that the average mean of data (red dots) almost coincides with the expected value (black curve). Furthermore, the experimental data lie within a band around the expected value with a width of two standard deviations. Figure 3(b) illustrates a simulation of the drug concentration in five individuals when the dosage is continuous (infusion). Here, solution of the differential equation (2) with f .t / D r a X p , was evaluated by the numerical approximation (19).

Conclusions
We study an one compartment stochastic differential model that describe absorption and elimination of a drug under a single and constant dosage regimens. We derived an integral equation for the solution, that allow to determine an explicit formulas for the expected value and variance. The model is identified using the formulas of the expected value and the variance, which allows to predict the realistic path of the solution and its uncertainty, as well as, to determine the therapeutic range of the drug. Further, a numerical method for the solution of the SDE is introduced and developed, in this case, using a quadrature rule; instead of solving the SDE system by standard numerical methods. We leave the convergence analysis of the numerical method for future work.
a i a j :

Identity
Proof. In the left hand we make change of variablé

B. Calculation of variance for a constant dose
From equation (11) we get where the integrals are and Ä j .t u/ are given by the equations (12) and (13) respectively. Using the identity (20) with n D 2, we obtain To compute E I 2 2 .t / of process is similar to previous one, then If r e ¤ r a and r e ¤ r a C c 2 then By Fubini's Theorem we have that: Substituting (13) in the above expression, it results Thus, from (15), (22), (24) and (26) we obtain r e X a0 r a X p / r e .r e r a / Ã 2 e 2rt Â 2c 2 r 2 a X p r 2 e Ã Â r e X a0 r a X p .r e r a /.r e c 2 / Ã e r e t C Â c 2 r 2 a .X a0 X p / 2 .r a r e / 2 .2.r e r a / c 2 / Ã e 2r a t 2c 2 r 2 a X p .X a0 X p / r e .r a r e /.r e c 2 / e r a t 2r a .X a0 X p / r e .r e r a / 2 Â c 2 .r e X a0 r a X p / .r e r a / 2 X p r e r a c 2 Ã e .r e Cr a /t C 2r 2 a X p .X p X a0 / .r e r a /.r e c 2 / e .r e Cr a c 2 /t C c 2 r 2 a X 2 p r 2 e .2r e c 2 / C Â 2r 2 a .X p X a0 / 2 .r e r a c 2 /.2.r e r a / c 2 / C 2r 2 a X p .c 2 r e r a /X p C .2r e c 2 /X a0 .r e r a /.r e c 2 /.2r e c 2 / ! e 2.r e 1