INTERNATIONAL JOURNAL OF COMPUTING SCIENCE AND APPLIED MATHEMATICS VOL.
NO.
MARCH 2024
Integration-Based Method as an Alternative Way to Estimate Parameters in the IV Bolus Compartment Model Diny Zulkarnaen*.
Fadilah Ilahi.
Muhammad Syifa Irfani, and Dani Suandi AbstractAiAn alternative method of integration-based parameter estimation applied in pharmacokinetics problems is proposed The method, introduced by Holder and Rodrigo .
, is used to estimate the rate of drug elimination and distribution when it enters the body via intravenous bolus.
The estimation results are then compared with the classical method, the least squares method for the one-compartment model, and the residual method for the two-compartment model.
Graphical simulations of drug concentration versus time are also performed in this article to view not only the dynamics of drug delivery in the body, but also the comparisons between the approximate solutions and the arbitrarily generated data points.
Comparisons are also presented when the data points take into account noise in the form of random values.
Based on the estimation and simulation results, the integration-based method gives good results and even better than the classical method although when noise is applied to the data points.
Index TermsAiPharmacokinetics.
Compartment Model.
Least Squares.
Residual.
Integration-based Method.
I NTRODUCTION
HE Pharmacokinetic compartment model describes the dynamics of the drug delivery system when drugs enter the body by certain routes such as oral.
IV bolus, intramuscular, and rectal.
When a drug enters the body, it undergoes four fundamental processes, namely absorption, distribution, metabolism and excretion.
The last two processes can also be referred to elimination process.
The process of drug delivery can be modeled mathematically by an ordinary differential equation (ODE) which represents the rate of change of drug concentration in the body .
, .
The number of ODEs used in the model tells us about the number of compartment that the body has.
When all organs or body tissues are considered to have homogenous rapid drug perfusion, then a one-compartment model with an ODE can be constructed.
However, when the body is supposed to have two different classeses of drug perfusion: fast and slow, then a twocompartment model is more suitable, and so on.
From the model, we can find the a solution to see the dynamics of drug concentration either analytically, for example using the superposition method .
or the Laplace transform .
, .
, .
, or numerically using finite difference method .
, .
Another Zulkarnaen.
Ilahi, and M.
Irfani are with the Department of Mathematics.
Universitas Islam Negeri Sunan Gunung Djati.
Bandung, 40614 Indonesia e-mail*: dzulkarnaen@uinsgd.
Suandi is with the School of Computer Science.
Bina Nusantara University.
Jakarta, 11480 Indonesia e-mail: dani.
suandi@binus.
Manuscript received August 24, 2023.
accepted November 28, 2023.
way to represent the drug delivery process is to construct a mathematical model that takes the form of a fractional-order system .
, .
, .
, .
In addition to finding the exact solution, the parameters that appear in the model can be estimated based on the data collection of drug concentration with various times.
The least squares method .
and residual method .
are very well-kown that are commonly used to estimate parameters, especially for intravenous bolus (IV bolu.
Meanwhile the Wagner-Nelson .
and Loo-Riegelman .
, .
methods are also the most common methods used to estimate parameters for oral administration for one and twocompartment models, respectively.
In this paper, we propose an alternative method for parameter estimation called integration-based method, introduced by Holder and Rodrigo .
, for one and two-compartment models, where the drug is given via IV bolus route.
The method is implemented using data in the form of drug concentration versus time so that the parameters in the model can be There are two types of data points provided in the implementation: the data points with smooth pattern, and the other is the data points with noise.
In the end this method is compared with the classical method of least squares and II.
T HE L EAST S QUARES AND R ESIDUAL M ETHODS
In this section the original one and two-compartment models are presented .
, .
, .
, in which the drug is given by IV bolus route.
From the given model, we use classical methods to estimate the parameters that appear in the model.
The least squares method for one-compartment model, and the residual method for two-compartment model.
One-Compartment IV Bolus Model Here, we assume that the body has homogeneous rapid blood flow, so a one-compartment model containing the rate of drug concentration C.
is utilized to investigate the drug delivery process, specified by = Oeke C.
= C0 .
As can be seen in the model, there is only one parameter appears, that is ke which indicates the elimination rate.
Meanwhile C0 , the constant initial drug concentration in the blood, is assumed to be known.
To estimate this parameter.
INTERNATIONAL JOURNAL OF COMPUTING SCIENCE AND APPLIED MATHEMATICS VOL.
NO.
MARCH 2024
the classical method of least squares fitting is employed.
The first step that needs to be done is to find the exact solution of the model .
By integrating both sides we can easily obtain the solution as C.
= C0 eOeke t .
From here, we linearize the solution to become ln C.
= Oeke t ln C0 with the aim that the least squares method can be applied to find the slope of .
By having the slope, we can estimate the parameter value of the elimination rate by the formula m i=1 ti ln Ci Oe i=1 ti i=1 ln Ci ke = Oe m i=1 t2i Oe ( i=1 ti )2 where m is the number of data points.
Two-Compartment IV Bolus Model When the body is considered to consist of two classes of blood flow: rapid and slow, then a two-compartment model is an appropriate choice to describe the drug delivery process.
The compartments with rapid and slow blood perfusion are called the central and peripheral compartments, respectively.
Thus, the model can be set up as a system of two differential equations of the rate of drug concentration in the central C1 .
and peripheral C2 .
compartments, governed as = k21 C2 Oe .
12 ke )C1 .
C1 .
= C0 , .
= k12 C1 Oe k21 C2 .
C2 .
= 0.
Notice that there are three parameters that appear in this model, namely k12 , k21 , and ke , where the first two parameters denote the respective rate of drug distribution from central to peripheral and from peripheral to central.
Therefore, these three parameters will be estimated, the residual method is employed here.
Like in the one-compartment model, the solution of .
also needs to be found.
One way to produce such a solution is using the Laplace transform method, and we obtain C0 C1 .
= (Oe1 k21 )eOe1 t Oe (Oe2 k21 )eOe2 t , 2 Oe 1 C0 k12 Oe1 t C2 .
= Oe eOe2 t , 2 Oe 1 where 1 < 2 and 1 2 = k12 k21 ke , 1 2 = k21 ke .
Now we have the supporting parameters that need to be obtained, i.
1 and 2 , before the main parameters can be Starting from the solution C1 .
we can estimate 1 as well as k21 using large time data points, say the last m1 data In other words, when time is large we can consider that the first term in the solution .
of C1 .
is close to zero since 2 > 1 .
This means the solution of C1 .
becomes C0 (Oe1 k21 ) Oe1 t C1 .
= 2 Oe 1 which can be transformed into the linear equation as C0 (Oe1 k21 ) ln C1 .
= Oe1 t ln 2 Oe 1 From this equation we can calculate the slope and the yintercept by the least squares method, so that we can estimate 1 as Pm1 Pm1 Pm1 ln C1,i ti ln C1,i Oe i=1 m1 i=1 Pm1i i=1 Pm1 2 1 = Oe m1 i=1 ti Oe ( i=1 ti ) Since the slope and the y-intercept of the equation .
has been estimated, we can approximate the solution in .
Consequently, the deviation between the original exact solution in .
and the approximate solution in .
can be calculated C0 .
Oe k21 ) Oe2 t R.
= 2 Oe 1 We can also call R.
as the residual equation.
Linearizing this equation we can calculate the slope to obtain 2 , and the y-intercept to obtain 2 and k21 by using the first m2 data points with m2 < m as Pm2 Pm2 Pm2 m2 i=1 ti ln Ri Oe i=1 ln Ri Pm2 2 Pm2i i=1 2 = Oe , .
m2 i=1 ti Oe ( i=1 ti )2 .
Oe 1 ) y k21 = 2 Oe Pm2 i=1 ti Pm2 Pm2 Pm2 ln R Oe i=1 t ln R Pm2 i 2 Pm2i 2 i i=1 i .
m2 i=1 ti Oe ( i=1 ti ) Observe that k21 can be found by the formula .
since C0 is given and 1 as well as 2 can be calculated by the respective formulas given in .
Then, the parameters ke and k12 can be estimated sequentially using the formula .
as 1 2 k12 = 1 2 Oe k21 Oe ke .
T HE INTEGRATION - BASED METHOD
The integration-based method, proposed by Holder and Rodrigo .
, is utilized here as an alternative way to estimate the parameters that appear in the ODE.
It can be used as long as the parameters to be estimated is linear in ODE.
The general step that needs to be done is to define the weight function followed by performing the integration to both sides of the ODE.
One-compartment model To apply the integration-based method to the onecompartment model .
, the weight function w.
with s > 0 is multiplied to both sides of .
, then integration is performed between time 0 to T .
As a result the equation .
becomes Z T Z T w.
= Oeke Then, applying the integration by parts on the left-hand side of the latter equation, we now have INTERNATIONAL JOURNAL OF COMPUTING SCIENCE AND APPLIED MATHEMATICS VOL.
NO.
MARCH 2024
Z T
w(T .
C(T ) Oe w.
Oe wA .
Z T
= Oeke w.
Note that we choose the weight function which takes the form w.
= eOest since this form has been used by .
, .
From the latter equation, we can estimate the elimination rate by the formula w(T .
s0 )C(T ) Oe w.
s0 )C.
Oe 0 wA .
s0 )C.
dt ke = Oe w.
s0 )C.
for a given value s0 > 0.
Here the integral term in .
is calculated numerically .
n scilab, we use intspline functio.
Remark 1.
The idea of using and choosing the weight function w.
= eOest is that, based on the original article .
, when this weight function is multiplied and integrated to an equation, itR will be analogous to applying the finite Laplace transform 0 eOest f .
But one can use another type of weight function.
Since this paper focuses on comparing the alternative method with the classical method, so to keep on that focus, we restrict this alternative method to only use one type of weight function, that is w.
= eOest .
Thus, k12 and k21 can be obtained provided the inverse of the leftmost matrix exists.
Next step, the first ODE in .
is employed to estimate ke .
Like the previous step, we obtain the equation Z T wA .
C1 .
dt w(T .
C1 (T ) Oe w.
C1 .
Z T
Z T
= k21 w.
C2 dt Oe k12 w.
C1 dtOe
Z T
C1 dt.
Since k12 and k21 are now known, the last equation can be represented by A.
= B.
ke D.
Z T wA .
C1 .
= w(T .
C1 (T ) Oe w.
C1 .
Z T
Z T
= k21 w.
C2 dt.
= Oek12 w.
C1 dt.
Z T D.
= Oe w.
C1 dt.
Thus, the parameter ke can be estimated by Two-compartment model As previously stated that there are three parameters to be estimated from the two-compartment model given in .
The way to estimate these parameters based on the integrationbased method is estimating the parameters in each ODE sequentially, not simultaneously.
It should be noted that we use the same weight function for each ODE for consistency, which is the same function with the one-compartment model, namely w.
= eOest .
We first select the second ODE to estimate k21 and k12 .
Like the one-compartment model, we can transform the second ODE in .
to become Z T w(T .
C2 (T ) Oe w.
C2 .
Oe wA .
C2 .
Z T
Z T
= k12 w.
C1 dt Oe k21 w.
C2 dt, .
or simply rewrite the equation as P .
= k12 Q.
k21 R.
Z T P .
= w(T .
C2 (T ) Oe w.
C2 .
Oe wA .
C2 .
Z T
Z T
= w.
C1 dt.
= Oe w.
C2 dt.
Since there are two parameters to be estimated, we have to choose two different values of s, say s1 and s2 .
As a result, we can construct a matrix from the latter equation, written as Q.
1 ) R.
1 ) k12 P .
1 ) .
2 ) R.
2 ) k21 P .
2 ) A.
0 ) Oe B.
0 ) Oe C.
0 ) D.
0 ) .
for given s0 > 0.
IV.
S IMULATION
In this simulation, the parameters that appear in the one and two-compartment models are estimated not only using the alternative integral-based method but also the classical method of least squares and residual for the one and twocompartment models, respectively.
The estimation requires time versus drug concentration data points, where these data points are generated by the exact solutions given in .
with arbitrarily selected parameter values.
Then the generated data points are used to estimate the parameters using the integrate-based method as well as the classical methods to compare them with the parameter values that were set in the beginning.
Furthermore, the generated data points are given noise and investigate whether there is a significant difference in estimation between the two methods.
Finally, several graphical simulations are presented to see not only the dynamics of drug concentration in the body, but also the comparison between the generated data points and with the approximate solutions.
One-compartment model In this part, we do not use randomly selected timeconcentration data points, but rather follow the trend of the data points based on equation .
that decreases exponentially.
In other words, we generate data points from the model solution in .
by choosing the parameter values .
r we can the exact value.
that appear in this equation, namely ke and INTERNATIONAL JOURNAL OF COMPUTING SCIENCE AND APPLIED MATHEMATICS VOL.
NO.
MARCH 2024
C0 .
The reason we use this way of generating data points is that from the data points we have generated, we will obtain parameter estimates which are calculated in Section II and Section i.
Subsequently, we compare these estimated values with the pre-selected .
The closer the estimated values to the exact parameter values using a certain method, the more robust the method is.
Here we choose the elimination rate ke = 0.
2184/hour and the initial concentration C0 = 24 mg/ml to generate the data points by inserting these two values to the exact solution presented in .
We use time interval 1 < t < 20 with a one-hour time step, which means we have hourly drug concentration data points for up to 20 hours.
The generated data points can be seen in Table I.
TABLE I: Generated data points for the one-compartment t .
C .
g/m.
As we can see in the table, the drug concentration is large in the beginning and then decreases to zero.
This happens since there is no absorption process so that the drug enters the blood directly, then undergoes elimination process.
By the data points in Table (I), we can now estimate the elimination rate using the formula given in .
and in .
by choosing s0 = 0.
01 for the respective least squares and integration-based methods.
The result for these calculations is ke = 0.
2184 for both methods, which is the same as the original values.
data points integrationOebased method least squares fitting generated as shown in Table II.
It is important to note that the rule of generating data points here is the same as the case of the one-compartment model.
TABLE II: Generated data points for the two-compartment t .
C1 .
g/m.
C2 .
g/m.
Notice that, unlike the central compartment, where the concentration C1 decreases over time, the drug concentration in the peripheral compartment C2 increases first, then at some point starts to decrease towards zero.
This occurs because the amount of drug in the peripheral compartment yield drug supply from the central compartment.
Then when the drug concentration in the central compartment starts to decrease, the supply also decreases which results in the reduction of drug concentration in the peripheral compartment.
TABLE i: Exact and estimated parameter values for the twocompartment model.
Methods Exact Integration-based Residual The generated data points given in Table II are then implemented to estimate the parameters using the residual and the integration-based methods.
The estimation values can be seen in Table i.
These estimations involve the formulas presented in .
for the residual method, and the formulas in .
for the integration-based method.
When we Fig.
1: Comparison between data points and the approximate solutions for the one-compartment model.
As a result the graphs for the generated data points and the approximate solutions give the overlapping lines as depicted in Figure 1.
Two-compartment model Recall that there are three parameters involved in the twocompartment model stated in .
Here, we input the parameters k12 = 0.
93/hour, k21 = 0.
32/hour and ke = 0.
22/hour together with the initial concentration C0 = 20 mg/ml to the exact solution .
so that the generated data points can be Fig.
2: Comparison between data points and the approximate solutions for the two-compartment model.
compare the estimates of k12 , the residual method gives the better approximation.
While for the other two parameters, the integration-based method gives better results than the residual Now, when a graphical simulation is made and it can be observed in Figure 2, the least squares method gives a better approximation for the drug concentration in the peripheral INTERNATIONAL JOURNAL OF COMPUTING SCIENCE AND APPLIED MATHEMATICS VOL.
NO.
MARCH 2024
In contrast, the integration-based method provides better approximation for the drug concentration in the central compartment.
TABLE IV: Exact and estimated parameter values for twocompartment model with noise data set.
Methods Exact Integration-based Residual In actual experiments, the data collection of drug concentrations over time may have an unsmooth graph due to measurement errors or other factors that make data collection values Therefore, we make a noise from the generated data .
Fig.
3: Comparisons between data points with noise and the approximate solutions for the two-compartment model.
points by adding random values to the drug concentrations for C1 and C2 .
Then, both methods are reapplied using the new data points to perform parameter estimate.
The results can be seen in Table IV and Figure 3.
From Table IV, it can be seen clearly that the integrationbased method gives a much better estimate than the residual The cause of the inconsistent results that we can analyze from the residual method is the semilog graph which forces the data points to be linear as given in equation .
, whereas the truth is the data points are now no longer smooth so that for the noise case they no longer form a linear trend as shown in Figure 4b.
We can compare the semilog graph when the noise is not present as depicted in Figure 4a which is linear since the data points are smooth.
This forced linearity produces unexpected results from the residual method, where the approximation of drug concentration in both the central and peripheral compartments lie too far from the data points, pictured as broken lines in Figure 3.
Fig.
4: Semilog graph of large-time data points, .
without noise and .
with noise, for the two-compartment model.
least squares fitting in estimating parameters for the onecompartment model.
As for two-compartment model, the method integration-based method provides better estimates than the classical residual.
This is based in the closeness of the estimated parameter values to their original values as given in Table II, that is the integration-based method produces two better parameter estimates, while the residual method only produces one better parameter estimate.
When the data points are perturbed by a noise, it turns out that the integration-based method gives consistent estimates, unlike the residual method which gives estimates that are much different than the data points without noise .
ee the comparison in Table II and Table IV).
As a result, the drug concentration dynamics are far from the generated data points as depicted in Figure 3.
With these statements we can conclude that the alternative method, namely the integrationbased method gives good approximation and even better than the classical method despite the perturbation in the data points.
C ONCLUSION
The integration-based method can be used as an alternative way to estimate parameters in pharmacokinetic problems specifically in one and two-compartment model.
This method gives as same precision as the classic method of ACKNOWLEDGMENT