A Stochastic Analysis of the One Compartment Pharmacokinetic Model Considering Optimal Controls

First-order 1-compartment pharmacokinetic model for extravascular administered drugs can be used to derive many useful quantities by comparing the predicted values with actual data. However, less research has been done in actually formulating them as optimal control problems. Moreover, real pharmacological processes are always exposed to influences that are not completely understood or not feasible to model explicitly. Ignoring these phenomena in the modeling may affect the estimation of PK/PD models’ (pharmacokinetic/pharmacodynamic models’) parameters and the derived conclusions.Therefore there is an increasing need to extend the deterministic models to models including a stochastic component. In our study, we modify the 1-compartment pharmacokinetic model to a stochastic differential equation model based on an optimal control problem. A schedule of optimal dosing and timing has been given from our proposed model.


I. INTRODUCTION
Pharmacokinetics (PK) is the study of what the body does to a drug. It studies the absorption, distribution, metabolism, and excretion of the medicine (ADME), as well as bioavailability. PK analysis forms a major part of the understanding and development of the Investigation Medicinal Product (IMP), and can also contribute heavily to the prescription once a drug has been approved.
There are three approaches that have been suggested for pharmacokinetic modeling: compartmental, physiological and model-independent. The first one is an empirical approach, which is based on simple compartmental models. These compartments have no strict physiological or anatomical basis. The compartment simply represents a body volume, or just as easily it could represent a chemical state, for example a metabolite of the drug. Usually this approach uses one or two compartments. Despite its simplistic nature, many useful quantities can be derived using this approach, and by comparing predicted values with actual data.
In [1], Sophie Donnet and Adeline Samson (2013) reviewed the examination of the pertinence of stochastic differential equations (SDEs) for pharmacokinetic/ pharmacodynamic models. A natural extension of deterministic The associate editor coordinating the review of this manuscript and approving it for publication was Feiqi Deng . differential equations model is a system of SDEs, where relevant parameters have been modeled as suitable stochastic processes, or stochastic processes have been added to the driving system equations [2]. The first papers encouraging the introduction of random fluctuations in PK/PD were authored by D'Argenio [3] and Ramanathan [4], [5]. The authors underline that PK/PD have contributions from both deterministic and stochastic components: drug concentrations follow determinable trends but the exact concentration at any given time is not completely determined. For example, Ramanathan [5] proposes a stochastic one-compartment PK model with a variable elimination rate. More sophisticated PK/PD models then have been proposed with multiple compartments, non-linear or time-inhomogeneous absorption or elimination (seen for example Ferrante et al. [6]; Tornøe et al. [7]; Ditlevsen and De Gaetano [8]; Ditlevsen et al. [9]; Picchini et al. [10]).
Although lots of mathematical models for pharmacokinetics has been build, less research has been done in actually formulating them as optimal control problems. Actually, suitable control such as food intake, special excercises, some vitamin, can affect the drug absorption. As in [11], the authors studied significant effect of infection and food intake on sirolimus pharmacokinetics and exposure in pediatric patients with acute lymphoblastic leukemia. A mathematical model for the depletion of bone marrow under cancer chemotherapy VOLUME 8, 2020 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ is analyzed as an optimal control problem in [12]. Reference [13] investigated multiscale tumor modeling that integrates drug pharmacokinetic and pharmacodynamic (PK/PD) information using stochastic hybrid system modeling framework. Reference [14] presented a new model to describe the single-dose pharmacokinetics of bevacizumab and predict its multiple-dose pharmacokinetics in beagle dogs. Reference [15] combined health record informatics and pharmacokinetic modeling and got a powerful translational approach to detect high dimensional drug-drug interactions.
Through clinical trial, [16] proposed a model based optimization of G-CSF treatment during cytotoxic chemotherapy and showed validity of model predictions regarding alternative G-CSF schedules.
In this study, we investigate to combine optimal control theory and stochastic methods to propose a novel pharmacokinetic model. Our aim is to find the optimal drug dosing schedule and predict the absorption rate and concentration rate. We added a control vector to the pharmacokinetic model and optimal control theory was used to analyze the modified model. We found the optimal dosing timing schedule and an equilibrium point of the dynamic system. Near the stationary point, white noise was added to the modified model, and the ensuing stochastic differential equations (SDEs) were presented. We proved existence, uniqueness and stability of the SDE system and found an explicit solution. Finally, the model was simulated and parameters of the model was estimated by using the R language and the stability of the numerical method was proven.

A. A BRIEF REVIEW OF FIRST-ORDER 1-COMPARTMENT UNCONTROLLED MODEL (Extravascular Administration)
This approach models the entire body as a single compartment into which a drug is added by a rapid single dose, or bolus. It is assumed that the drug concentration is uniform in the body compartment at all times, and is eliminated by a first order process that is described by a first order rate constant K 10 : where A a = Amount of drug absorption deposit, A c = Amount of drug in central compartment, K a = Absorption rate constant, K 10 = first-order elimination rate, indicating elimination of drug out of the central compartment into urine, feces, etc.
It's easy to get the solution of the system (1): where A a (0) is the initial amount of drug in the gastrointestinal tract.

B. STABILITY OF THE SDEs
Definition 1: [17] Let ( , F, P) be a complete probability space with a filtration {F t } t≥0 satisfying the usual conditions.
with initial value x(t 0 ) = x 0 . Assume that for any initial value x(t 0 ) = x 0 ∈ d , equation (2) has a unique global solution which is denoted by x(t; t 0 , x 0 ). We know that the solution has continuous sample paths and its every moment is finite. Assume furthermore that f (0, t) = 0 and g(0, t) = 0 for all t ≥ t 0 .
So equation (2) has the solution x(t) ≡ 0 corresponding to the initial value x(t 0 ) = 0. This solution is called the trivial solution or equilibrium position. Definition 2: [17] Let K denote the family of all continuous nondecreasing functions µ : + → + such that is said to be decrescent (i.e. to have an arbitrarily small upper bound) if for some µ ∈ K, Definition 3: [17] The trivial solution of equation (2) is said to be (i) stochastically stable or stable in probability if for every pair of ∈ (0, 1) and r > 0, there exists a δ = δ( , r, t 0 ) > 0 such that Otherwise, it is said to be stochastically unstable.
Let 0 < h ≤ ∞. Denote by C 2,1 S h × + ; + the family of all nonnegative functions V (x, t) defined on S h × + such that they are continuously twice differentiable in x and once in t. Define the differential operator L associated with equation (2) by If L acts on a function V ∈ C 2,1 S h × + ; + , then Theorem 1: [17] If there exists a positive-definite (2) is stochastically asymptotically stable in the large. Definition 4: [17] The Euler-Maruyama approximate solutions are defined as follows: For every integer n ≥ 1, define x n (t 0 ) = x 0 , and then for t 0 + (k − 1)/n < t ≤

. Here A and B i s are all the d ×d matrices while a and b i s are d-dimensional vectors.) is stable in distribution. In particular, for a scalar linear SDE
Let t be the set of all possible outcomes (or realisations) at the point t, and define the random variable Y t as the function Y t : t → . Define the set of possible outcomes over all time as = ⊗ ∞ t t , and the random variables X t : → , where for every ω ∈ , with ω = (ω 0 , ω 1 , · · · ), we have X t (ω) = Y t (ω t ). Hence we have a sequence of random variables {X t } t (which we call a random process). When we observe {x t } t , this means there exists an ω ∈ , such that X t (ω) = x t . To complete things we have a sigma-algebra F whose elements are subsets of and a probability measure Definition 5: We say that the sequence {X t } converges almost sure to µ, if there exists a set M ⊂ , such that P(M ) = 1 and for every ω ∈ M , we have In other words for every > 0, there exists an N (ω) such that

II. USING OPTIMAL CONTROL THEORY TO STUDY THE FIRST-ORDER 1-COMPARTMENT MODEL
Pharmacokinetic (PK) equations model the time evolution of the drug's concentration in the body/plasma. Let u denote the drug dosage with u = 1 corresponding to a maximal dose and u = 0 denoting no treatment. Simple models considered in the literature (for example, see [19], [20]) use a first-order linear systemċ where f and h are positive constants to represent the dynamics for the drug concentration c in the plasma. The model itself is one of exponential growth/decay as it is commonly used as model for continuous infusions. As in [12], a more generally bilinear system of the form was proposeḋ with an additional parameter g added.
Here we add a control vector to equation (1) of the first-order 1-compartment model, and derive the following first-order 1-compartment dynamic system, where α and β are positive constants representing the dynamics for the drug and concentration in absorption deposit A a and central compartment A c in the plasma, respectively. VOLUME 8, 2020 Let A := a 11 a 12 a 21 a 22 Then, we can write (3) in matrix notation asẋ = Ax + Bu. By Von Karman controllability, the necessary condition for controllability is that the matrix B AB has full rank, i.e., B AB = 0, i.e., α β The fact is that maximizing the Hamiltonian function H with respect to u, is only possible if u is bounded, i.e. |u| ≤ u 0 for some u 0 .
As we mentioned above, in our system, the drug dosage u is between 0 and 1 which surely satisfies the above bounded condition.
In order to avoid drug side effects, we hope to achieve the best possible therapeutic effect with as few medications as possible. So in this paper we consider to a performance index in the form (4) where T represents the time betweeen two treatments, r 1 , r 2 > 0, q 1 , q 2 ≥ 0, and b > 0 are constants.
In the objective (4), we have incorporated a term q 1 A a (t) + q 2 A c (t) in the Lagrangian in an effort to achieve the best possible therapeutic effect. At the same time, in order to avoid drug side effects, we also want to schedule as few drug dosage as possible by including a term b(1 − u(t)) in the objective. In addition we have added a terminal term r 1 A a (T ) + r 2 A c (T ) which represents a weighted average of the total drug concentration at the end of an assumed fixed therapy interval [0,T] in order to prevent that the drug concentration would be too low to keep certain effectiveness towards the end of the therapy interval.
The mathematical problem therefore can be formulated as to maximize (4) over all Lebesgue measurable functions u which take values in [0, 1] subject to the dynamics (1) respectively (3) and given initial conditions.
By (4) and (3) respectively, we denote f 0 (x 1 , where ψ 1 and ψ 2 are the adjoint functions associated with constraints (3) and the given initial conditions. We will obtain them from the following discussion.

Define the Hamiltonian
By the Pontryagin's maximum principle [21] [22], if u * is the optimal control (optimal controls u * maximize the Hamiltonian H, i.e., (ψ 1 α +ψ and x * is the state trajectory, then there exists an absolutely continuous function ψ(t) := (ψ 1 (t), ψ 2 (t)) defined on [0, T ] satisfying the adjoint equations with transversality condition, such that the following condition is satisfied: the optimal control u * maximizes the Hamiltonian H (see (5)) over the control set [0, 1] along ψ 1 (t), ψ 2 (t), A * a , A * c , u * (t) . We call a pair ((A a , A c ), u) consisting of an admissible control u with corresponding trajectory (A a , A c ) for which there exist multipliers (ψ 1 , ψ 2 ) such that the conditions of the Maximum Principle are satisfied an extremal (pair) and the triple ((A a , A c ), u, (ψ 1 , ψ 2 )) is an extremal lift (to the cotangent bundle).
From Pontryagin's Maximum Principle [21] [22], the optimal control u * maximizes H as a function of u. Since H is linear in u and u ∈ [0, 1], the maximum value of H is Since u = 1 represents maximal dosage and u = 0 corresponds to no drug dosed, we have that, in principle, optimal controls alternate between sessions of ''full control'' or rest periods and partial controls are not optimal.
By solving the above ODE, we get where C 1 is a constant. By the transversality condition ψ 2 (T ) = r 2 , we get C 2 = r 2 e −K 10 T . So ψ 2 = r 2 e K 10 (T −t) .
By the transversality condition ψ 1 (T ) = r 1 , and C 2 = r 2 e −K 10 T , we obtain It implies C 1 = − r 1 + r 2 K a K 10 − K a r 2 K a e (K a −K 10 )T . The necessary condition for t being a switching time is that t satisfies i.e., the swithching function (t) = αψ 1 (t) + βψ 2 (t) − b is strictly monotonous on [0, T ], there's at most one switch. Therefore, by the above two conclusions, there is one and only one switch in the control.
The switching time can be obtained by solving − C 2 K a e K 10 t K 10 − K a − C 1 C 2 K a e K a t α + C 2 e K 10 t β − b = 0, using numerical method.
Optimal trajectories satisfẏ The equilibrium point P is Since both eigenvectors of matrix A, λ 1 = −K a and λ 2 = −K 10 , are negative, the positive-plane has stable node at P = αu * K a (α+β)u * K 10 .

III. THE FIRST-ORDER 1-COMPARTMENT STOCHASTIC MODEL
From the above discussion, we see that the optimal control occurs at u * . Thus, we revise the first-order compartment optimal control model (3) as We assume that the stochastic perturbations of the variables around their values given in (8), are white noise type, which are proportional to the distances of A a , A c from the values A * a , A * c . We then arrive to the system where σ 1 and σ 2 are real constants, and can be defined as the intensities of stochasticity, and ξ t = (ξ 1 t , ξ 2 t ) is a 2-dimensional white noise process. We wonder whether the dynamical behavior of model (9) is robust with respect to such a kind of stochastic perturbations by investigating the asymptotic stochastic stability behavior of equilibrium P for (10), and comparing the results with those obtained from the system (9). 10 , then A a = X 1 + αu * K a and A c = X 2 + (α+β)u * K 10 . Substituting them into the above system (10), and after simplification, we arrive to a first order 1-compartment SDE model: Using the Itô s formula [23], it's easy to solve dX 1 = −K a X 1 dt + σ 1 X 1 dξ 1 t and get which is a so-called geometric Brownian motion (GBM). If ξ 1 t is independent of X 1 (0), we have that and Var(X 1 (t)) = X 2 1 (0)e −2K a t (e σ 2 1 − 1).
The probability density function of X 1 (t) is f X 1 (x, K a , σ 1 , t) If K a + σ 2 1 2 = 0, then X 1 (t) will fluctuate between arbitrary large and arbitrary small values as t → ∞, P − a.s.

Theorem 3: The solution of the stochastic model (11) is unique.
Since the stochastic system model (11) is a linear system, the uniqueness of its solution can been obtained directly from Theorem 2.1 in Reference [17].
Proof: Let the Lyapunov function be V (x, t) = By the application, we know K a > 0 and K 10 > 0, hence Moreover, Thus, Q is a symmetric positive-definite matrix as required.

V. QUANTITATIVE ANALYSIS OF THE FIRST-ORDER 1-COMPARTMENT MODEL (EXTRAVASCULAR ADMINISTRATION) A. PARAMETER ESTIMATION OF THE FIRST-ORDER 1-COMPARTMENT MODEL (Extravascular Administration)
In the original ODE model, in order to estimate the parameters K a and K 10 in the model, firstly, we generated sampled data with K a = 2, K 10 = 1 and with initial concentrations being A a = 0.957, A c = 0.031. Then we set the initial values for the optimizer as K a = K 10 = 0.5, and we specify the coefficients drift and diffusion as expressions. We can now use the Levenberg-marquardt routine in package minpack.lm to estimate the parameters K a and K 10 of the model. The estimated coefficients are extracted from the output object fitmod as follows:  For u * = 0, we have done above. For u * = 1, in order to estimate the parameters K a , K 10 , α and β in the model, firstly, we generated sampled data with K a = 2, K 10 = 1, α = 0.16, β = 0.2, and with initial concentrations being A a = 0.957, A c = 0.031. Then we set the initial values for the optimizer as K a = K 10 = 0.5, α = 0.12, β = 0.03, and we specify the coefficients drift and diffusion as expressions. We can now use the Levenberg-marquardt routine in package minpack.lm to estimate the parameters K a , K 10 , α and β of the model.
The estimated coefficients are extracted from the output object fitmod as follows:

C. SIMULATIONS OF THE THREE FIRST-ORDER 1-COMPARTMENT MODELS
Using the above estimated parameters in V .A and V .B, we solved the three models numerically and plotted the simulation curves (Figure 1, 2, 3), then compared the three models visually. For 0 < u(t) < 1, all of the curves of (A a , A c ) are between curves on Figure 1 and Figure 2. It's clear that our revised ODE model with optimal control increased the amount of drug absorption deposit (A a ) and drug in central compartment (A c ), and SDE with optimal control model did the same improvement but with consideration of the influences that are not completely understood or not feasible to model explicitly. So the SDE with optimal control model improves the model and is more reasonable than the original ODE model and our revised ODE model with optimal control. If we keep the same (A a , A c ), then we need less drug dosage.

1) USING NUMERICAL METHOD TO VERIFY THE EXPLICIT SOLUTION OF THE FIRST ORDER 1-COMPARTMENT SDE MODEL
In this section, we'll give a theorem for the stability of the EM method [24], [25] for our SDE with optimal control model (11). Then we'll verify our explicit solutions of our SDE with optimal control model through comparing the explicit solution with the numerical solution from Euler-Maruyama(E-M) method.

Theorem 5: If
then for any sufficiently small step size , the E-M approximate solution of the SDE model (11) with the optimal control is stable in distribution.  i.e., then U is negative-definite, so by Theorem 2 for any sufficiently small step size , the E-M approximate solution of the SDE with optimal control model (11) is stable in distribution. In our SDE with optimal control model (11), K a = 2.07126, K 10 = 1.03649, σ 1 = 0.5, and σ 2 = 0.1, so Thus for any sufficiently small step size , the E-M approximate solution of our SDE with optimal control model (11) is stable in distribution.
In order to compare the explicit A a and the numerical A a , we draw them on the same set of axes ( Figure 5), and to show the explicit A c and the numerical A c overlap each other, we also draw them on the same set of axes (Figure 7). In order to compare the explicit A a and the numerical A a , we draw them on the same set of axes ( Figure 5), and to show the explicit A c and the numerical A c overlap each other, we also draw them on the same set of axes (Figure 7).
From Figure 4 and Figure 5, we can see the two curves of explicit A a and the numerical A a completely overlap each other, thus the explicit solution of A a is correct. From Figure 6   and Figure 7, we can see the two curves of explicit totally A c and the numerical A c also overlap each other, thus the explicit solution of A c is also correct.

VI. CONCLUSION
In this paper, we combined optimal control theory and stochastic differential equation and proposed an optimal controlled pharmacokinetic SDE model which provides an optimal dosing and timing schedule. The original ODE PK model was augmented through optimal control analysis and moreover, in order to make the model more reasonable, the optimal control model further was augmented into SDE model by considering the stochastic perturbations of the variables around their values at optimal control equilibrium point being white noise type, which are proportional to the distances of A a and A c from values at the optimal control equilibrium point.
We got the explicit solution of our SDE model with optimal control and proved the uniqueness and stability of explicit solution. In part V, we discussed our model quantitatively. We estimated the parameters of the model and compared the simulation curves of the three model and found that our our SDE model with optimal control increased the amount of drug absorption deposit (A a ) and drug in central compartment (A c ) and SDE more reasonable. Then we provided the condition of stability of E-M method and used it to get the numerical solutions. Then we verified our explicit solution is correct through comparing the curves of simulated the explicit and numerical solution.
Reference [12] analyzed a model for cancer chemotherapy that aims at minimizing the damage done to bone marrow cells during the chemotherapy, and concluded that partial doses are not optimal and in principle optimal contrals alternate between chemotheropy sessions of ''full dose'' and rest-periods. Reference [16] established a biomathematical model of human granulopoiesis under chemotherapy which allows predictions of yet untested G-CSF schedules, comparisons between them, and with it, optimization of flgrastim and pegfilgrastim treatment, and showed validity of model predictions regarding alternative G-CSF schedules by clinical trial results. From our study, we conclude theoretically that optimal dosing alternate between sessions of ''full dose'' or rest periods and partial dosage are not optimal through analyzing the PK model as an optimal control problem. Moreover, we found the switching time of ''full dose'' and rest-periods which can be used for drug dosage schedule easily. We can see that appropriate control can achieve the best possible therapeutic effect with as few medications as possible so that we can avoid drug side effects. And the stochastic differential equation makes the model more reasonable and the predictions about drug concentration more accurate.