A Pathwise Fractional one Compartment Intra-Veinous Bolus Model

Extending deterministic compartments pharmacokinetic models as diffusions seems not realistic on biological side because paths of these stochastic processes are not smooth enough. In order to extend one compartment intra-veinous bolus models, this paper suggests to modelize the concentration process $C$ by a class of stochastic differential equations driven by a fractional Brownian motion of Hurst parameter belonging to $]1/2,1[$. The first part of the paper provides probabilistic and statistical results on the concentration process $C$ : the distribution of $C$, a control of the uniform distance between $C$ and the solution of the associated ordinary differential equation, an ergodic theorem for the concentration process and its application to the estimation of the elimination constant, and consistent estimators of the driving signal's Hurst parameter and of the volatility constant. The second part of the paper provides applications of these theoretical results on simulated concentration datas : a qualitative procedure for choosing parameters on small sets of observations, and simulations of the estimators of the elimination constant and of the driving signal's Hurst parameter. The relationship between the estimations quality and the size/length of the sample is discussed.


Introduction
Compartments pharmacokinetic models describe the way an administered drug is transmitted among the body's compartments. The concentration of the drug in each compartment can be modeled by ordinary differential equations (cf. Y. Jacomet [13]). In particular, in one compartment models, the concentration is classically modeled by a linear (deterministic) differential equation with negative constant coefficient, taking in account the absorption and elimination steps. Only one compartment models are studied in this paper. By D. D'Argenio and K. Park [5], the elimination process has both deterministic and random components. A natural way to take in account these components is to add a stochastic noise in the linear differential equation that classically models the concentration. That has been well studied in the Itô stochastic calculus framework by many authors (cf. S. Donnet and A. Samson [10]). However, as mentioned in M. Delattre and M. Lavielle [8], since the standard Brownian motion has α-Hölder continuous paths with α ∈]0, 1/2[, the extension of the deterministic model as a diffusion is not realistic on biological side. M. Delattre and M. Lavielle force the paths regularity of the concentration process C by putting where D is the diffusion that extends the deterministic model. As mentioned in N. Marie [16], another way to increase the regularity of the concentration process paths is to replace the standard Brownian motion by a fractional Brownian motion B H of Hurst parameter H ∈]1/2, 1[ as driving signal. Since the signal is not a semi-martingale anymore, the stochastic integral is taken pathwise, in the sense of Young (cf. A. Lejay [15]). The Young integral keeps the regularity of the driving signal, therefore the concentration process has α-Hölder continuous paths with α ∈]0, H[. In both Itô and pathwise stochastic calculus frameworks, an interesting volatility function is x ∈ R + → σx β with σ ∈ R and β ∈ [0, 1]. It covers classical models : • β = 0, σ = 0 : Langevin equation. Its solution is the so-called Ornstein-Uhlenbeck process. • β = 1/2, σ = 0 : Cox-Ingersoll-Ross model. • β = 1, σ = 0 : Linear stochastic differential equation.
• σ = 0 : Linear ordinary differential equation. In the Itô stochastic calculus framework, that concentration model has been studied on statistical side in K. Kalogeropoulos et al. [14]. In the pathwise stochastic calculus framework, it has been studied on probabilistic side in N. Marie [16]. This paper is devoted to the probabilistic and statistical study of the special case of the one compartment intra-veinous (i.v.) bolus model with fractional Brownian signal : where τ 0 := inf {t ∈ R + : C t = 0} , the exponent β belongs to [0, 1[, υ > 0 is the rate of elimination describing the removal of the drug by all elimination processes including excretion and metabolism, and C 0 := A 0 /V with A 0 > 0 the administered dose and V > 0 the volume of the elimination compartment.
Since its vector field is C ∞ on bounded sets of R * + , equation (1) admits a unique continuous pathwise solution defined on [0, τ 0 ] and satisfying C . = X γ+1 . , where γ := β/(1 − β) and X is the solution of the following fractional Langevin equation : That equation is obtained by applying the rough change of variable formula to the process C and to the map x ∈ R + → x 1−β on [0, τ 0 ]. For details, the reader can refer to N. Marie [16]. The fractional Langevin equation is deeply studied in P. Cheridito et al. [3].
Since the concentration process has to be positive and stop when it hits zero, it can be defined as the solution of equation (1) For the sake of simplicity, even if the following equality only holds on [0, τ 0 ], throughout this paper, C is defined on R + by Note that for H = 1 and β = 0, the fractional Brownian motion is matching with t ∈ R + → ξt such that ξ N (0, 1), and That limit case illustrates that the Hurst parameter H is continuously controlling the regularity of the concentration process paths, but also that σ and H provides two complementary ways to control the impact of the random component on the elimination process with respect to its deterministic component.
In mathematical finance, the semi-martingale property of the prices process is crucial in order to ensure the market's completeness. The Itô stochastic calculus is then tailor-made to model prices in finance. In pharmacokinetic, the semi-martingale property of the concentration process seems not crucial on biological side. To replace the standard Brownian motion by a fractional Brownian motion in the pathwise stochastic calculus framework implies that the concentration process doesn't satisfy the Markov property anymore. In general, it makes the estimation of parameters υ, σ and H difficult, but the relationship between C and X mentioned above allows to bypass these difficulties by using results coming from Y. Hu and D. Nualart [11], J. Istas and G. Lang [12] and, A. Brouste and S. Iacus [2].
The second section is devoted to probabilistic and statistical properties of processes X and C. The first part concerns the distribution of the concentration process C and a control, in probability, of the uniform distance between the fractional Ornstein-Uhlenbeck process X and the solution of the associated ordinary differential equation. The second part provides a strongly consistent estimator of the elimination constant υ, and an extension of existing ergodic theorems for the fractional Ornstein-Uhlenbeck process X is established. The third part provides a strongly consistent estimator of (H, σ). A weakly consistent estimator of υ is deduced for unknown values of H and σ.
The third section is devoted to the application of the second subsection's theoretical results on simulated concentrations. For small sets of observations, the first part provides a qualitative procedure for choosing parameters H, σ and β. The cornerstone of the procedure is the control of the uniform distance between X and the solution of the associated ordinary differential equation mentioned above. The second part illustrates the convergence of estimators provided at Section 2. The relationship between the estimations quality and the size/length of the sample is discussed. Appendices A and B provide respectively useful definitions and results on fractional Brownian motions, and proofs of results stated at Section 2.

Probabilistic and statistical properties of the concentration process
The first subsection concerns the distribution of the concentration process C by using that C t = |X t | γ+1 ; t ∈ R + . Lemma 2.1 provides the covariance function of the fractional Ornstein-Uhlenbeck process X. Proposition 2.3 allows to control, in probability, the uniform distance between the process X and the solution of the associated ordinary differential equation. Refer to Appendix B for proofs of these results.
Essentially inspired by Y. Hu and D. Nualart [11], the second subsection provides an ergodic theorem for the concentration process and its application to the estimation of the parameter υ. An extension of existing ergodic theorems for the fractional Ornstein-Uhlenbeck process X is established. Essentially inspired by J. Istas and G. Lang [12] and, A. Brouste and S. Iacus [2], the third subsection provides a strongly consistent estimator of the parameter (H, σ). A weakly consistent estimator of υ is deduced for unknown values of H and σ.
2.1. Distribution of the concentration process and related topics. In order to provide the distribution of the concentration process C by using that C . = |X . | γ+1 , the following lemma provides first the covariance function of the fractional Ornstein-Uhlenbeck process X.
is a centered Gaussian process of covariance function R H,ϑ such that : for every s, t ∈ R + . Then, the covariance function R X of the fractional Ornstein-Uhlenbeck process X satisfies : for every s, t ∈ R + .
The following proposition provides the finite-dimensional distributions of the concentration process C.
Proposition 2.2. For every n ∈ N * and t 1 , . . . , t n ∈ R + , the distribution of the random vector (C t1 , . . . , C tn ) admits a density χ n with respect to the Lebesgue measure on (R n , B(R n )) such that : It is a straightforward application of Lemma 2.1 together with N. Marie [16], Proposition 5.1.
The following proposition allows to control, in probability, the uniform distance between the process X and the solution X det of the associated ordinary differential equation For every x > 0 and T > 0, For a level λ ∈]0, 1[, in order to ensure with probability greater than 1 − λ that An illustration of how one can use Proposition 2.3 and Corollary 2.4 in practice is provided at Subsection 3.1.

2.2.
Ergodic theorem and estimator of the elimination constant. The following proposition is an extension of existing ergodic theorems for the fractional Ornstein-Uhlenbeck process X. Let Y be the stochastic process defined by : Proposition 2.5. Let f : R → R be a continuous function such that : where a i := n − i and Then, by Proposition 2.5 : by Y. Hu and D. Nualart [11], Lemma 5.1. For n = 2, (3) coincides with Y. Hu and D. Nualart [11], Lemma 3.3.
Assume that values of parameters H and σ are known. For T > 0 arbitrarily chosen, consider . Proposition 2.6. υ T is a strongly consistent estimator of υ.

Estimators of the Hurst parameter and of the volatility constant.
Assume that the concentration process C is discretely observed at times t 0 , . . . , t n , where n ∈ N * , t k := kδ n for every k ∈ {0, . . . , n}, and (δ n , n ∈ N) is a R * + -valued sequence such that lim n→∞ δ n = 0 and lim n→∞ nδ n = ∞.
Proposition 2.7 provides a strongly consistent estimator, easy to implement, of the Hurst parameter H coming from J. Istas and G. Lang [12]. Proposition 2.8 provides an associated consistent estimator of σ. Proposition 2.9 provides a weakly consistent estimator of υ for unknown values of parameters H and σ.

Proposition 2.7. Consider
H n is a strongly consistent estimator of H.
. σ n is a strongly consistent estimator of σ.

Numerical simulations and pharmacokinetics
For small sets of observations, the first subsection provides a qualitative procedure for choosing parameters H, σ and β. Proposition 2.3 is the cornerstone of the procedure. The second subsection illustrates the convergence of estimators provided at Section 2. The relationship between the estimations quality and the size/length of the sample is discussed.

Consider the following values of the other parameters involved in equation
Parameters Values T 3h υ 3.5h −1 C 0 1g n 500 In order to choose H, σ and β, a qualitative procedure is provided by using these values as an example. That method is simple and doesn't require a lot of observations of the concentration process.
On one hand, as mentioned at Proposition 2.3, for a level λ ∈]0, 1[, in order to ensure with probability greater than Moreover, by Corollary 2.4 : On the other hand, as mentioned in introduction, the Hölder regularity of the concentration process paths is continuously controlled by the Hurst parameter H.    the concentration process paths seem not significantly perturbed with respect to the associated deterministic model. Then, to take β = 0.9 ensures that the value of the parameter σ can be chosen such that the following realistic condition is satisfied : On the observed concentrations c 1 , . . . , c n , the following procedure allows to choose H, σ and β qualitatively : • Step 1. Take H ∈]0.5, 1[ sufficiently close to 1 as 0.9.
If the value of υ is unknown, since paths of the concentration process have to be moderately perturbed with respect to the associated deterministic model, it can be approximated by linear regression as in Y. Jacomet [13] (see Subsection 3.2). • Step 4. Take σ 2 ∈]0; M (λ, x, H)] such that the concentration process paths seem regular enough locally and globally to model the elimination of the administered drug.
If the concentration process paths are not significantly perturbed with respect to the associated deterministic model for usual levels λ ∈]0, 1[, then go to the second step and choose a greater value of the parameter β. If the concentration process paths are not globally regular enough for standard levels λ ∈]0, 1[, then go to the second step and choose a smaller value of the parameter β.
Consider the following values of parameters involved in equation (1) : Parameters Values T nδ n ; n = 10, . . . , 10 The two following figures illustrate the convergence of estimators υ n and H n provided at propositions 2.6 and 2.7 respectively. For every n belonging to {10, . . . , 10 3 }, the concentration process C is simulated at times t 0 , . . . , t n and estimators υ n and H n are computed with these simulated observations denoted by c 1 , . . . , c n . Estimations are plotted in black and parameters values are plotted in red : The estimator υ n converges slowly to the elimination constant υ. Then, if the number n of observations is insufficient, since paths of the concentration process have to be moderately perturbed with respect to the associated deterministic model, one can take as in Y. Jacomet [13] : . . , t n ; log(c 1 ), . . . , log(c n )] var(t 1 , . . . , t n ) .
The qualitative procedure provided at Subsection 3.1 is also an alternative for choosing H and σ with few observed concentrations.

Discussion and perspectives
The stochastic model studied in this paper is a natural extension of usual deterministic models used in pharmacokinetics, it has smooth enough paths to take realistically in account the random component of the elimination process, and its explicit expression together with Decreusefond-Lavaud method allow to simulate it easily. As mentioned at Section 3, estimators of parameters υ, H and σ provide good estimations for large sets of observed concentrations. For small sets of observations, the qualitative procedure described at Subsection 3.1 is simple and seems quite efficient. For these reasons, the model could be used in clinical applications.
Assume that the therapeutic response R t to the administered drug at time t ∈ [0, τ 0 ] and O is a stochastic process with R-valued paths that doesn't depend on the initial concentration C 0 = A/V . The random variable C t is derivable with respect to C 0 > 0 and Differential calculus arguments could then allow to compute the dose that maximises the therapeutic response R t for some well chosen functions F and well chosen stochastic processes O.
Since the stochastic process C seems to model the elimination process more realistically than the deterministic function C det , the perspective of clinical applications described above could be interesting for potentially toxic drugs. For instance, the elimination of the ketamine, that can be neurotoxic but more effective than classic antidepressant in the treatment of major depressive disorders (cf. G.E. Correll and G.E. Futter [4]), could be modeled by the stochastic process studied in that paper. To choose F and O such that R t models the Hamilton rating scale or the Beck depression inventory at time t could allow to compute the dose of ketamine maximizing its antidepressant effect and minimizing its neurotoxic effect.

Appendix A. Fractional Brownian motion
Essentially inspired by D. Nualart [19] and, L. Decreusefond and A.S. Ustünel [7], this section gives basics on the fractional Brownian motion B H of Hurst parameter H ∈]1/2, 1[, its reproducing kernel Hilbert space and the fractional Young/Wiener integral with respect to B H . On Gaussian processes, the reader can refer to J. Neveu [18].
For a time T > 0 arbitrarily chosen, consider The process B H is a semi-martingale if and only if H = 1/2 (cf. [19], Proposition 5.1.1). Then, it is not possible to integrate with respect to B H in the sense of Itô. However, since E |B H t − B H s | 2 = |t − s| 2H for every s, t ∈ [0, T ], the Kolmogorov continuity criterion ensures that B H has α-Hölder continuous paths with α ∈]0, H[. Therefore, for any stochastic process X with β-Hölder continuous paths such that α + β > 1, it is possible to integrate X with respect to B H in the sense of Young. About the Young integral, that extends the classic Riemann-Stieljès integral, the reader can refer to A. Lejay [15].
That achieves the proof.
On the other hand, let ω be an element of Ω such that X(ω) − X det ∞,T x. In other words, for every t ∈ [0, T ], and so, by Jensen's inequality : Therefore, That achieves the proof by inequality (4).
Proof of Proposition 2.5. The ergodic theorem provided by A. Neuenkirch and S. Tindel at [17], Proposition 2.3 allows to conclude if f is in addition continuously differentiable. However, in the particular case of the fractional Ornstein-Uhlenbeck process, let show that the condition of Proposition 2.5 is sufficient.
Since Y is a centered, stationary and ergodic Gaussian process (cf. P. Cheridito et al. [3]), by the Birkhoff-Chintchin ergodic theorem together with the Fernique theorem : 1 T In order to conclude, it is sufficient to show that For T > 0 arbitrarily chosen : where p i (x) := c i (1 + |x|) bi ; x ∈ R, i ∈ {1, . . . , n}.