The basic premise for the clinical utility of pharmacokinetics is that there is a clearly defined relationship between drug concentration in readily available samples and drug response. Application of pharmacokinetic methods allows us to account for the variability in an individual's ability to absorb, distribute, metabolize and excrete drugs. The objective of these methods is to control the drug concentration in blood or plasma and potentially other fluids and tissues. For this approach to be effective there needs to be a relationship between these concentrations and the response to the drug. In some cases the response is direct and reversible and the Hill equation or some variation of it may be be applied.
In general the receptor might be 'mathematically' within the central compartment, a peripheral compartment or a separate effect compartment. In each case the relationship between concentration at the receptor and the response might be described with a Sigmoid Emax model (Hill equation)
In Equations 1 and 2, Emax is the maximum response, CR, 50% is the concentration which produces a 50% of maximum response (often called EC50%) and γ is a slope factor. For responses which increases from zero, Equation 1 is more appropriate. In other cases where the response is an increase or decrease from a baseline value Equation 2 may be a better choice. Another version of this equation, Equation 1, is the Emax model but this is just the same except that γ is set to 1.
A two compartment PK/PD model with an effect compartment
| Parameter | Value | Units |
|---|---|---|
| Dose | 200 | mg |
| k10 | 0.2 | hr-1 |
| k12 | 1.3 | hr-1 |
| k21 | 1.1 | hr-1 |
| V1 | 12.0 | L |
| k1e | 0.01 | hr-1 |
| -k1e | -0.01 | hr-1 |
| keo | 0.8 | hr-1 |
| Emax | 100 | percent |
| CR, 50% | 0.25 | mg/L |
| γ | 1.5 |
A two compartment PK/PD model with a response compartment
>cd /Directory/For/This/Analysis
>boomer
Boomer v3.4.5
Ref: Comput.Meth.Prog.Biomed.,29 (1989) 191-195
David W.A. Bourne
Copyright 1986-2018 D.W.A. Bourne
See http://www.boomer.org or email david@boomer.org for more info
Based on the original MULTI by K. Yamaoka, et.al.
J. Pharmacobio-Dyn., 4,879 (1981); ibid, 6,595 (1983); ibid, 8,246 (1985)
DATA ENTRY
0) From KEYBOARD
1) From .BAT file -1) From .BAT file (with restart)
2) From KEYBOARD creating .BAT file
3) From .BAT file (quiet mode) -3) From .BAT file (quiet mode-with restart)
4) to enter data only 5) to calculate AUC from a .DAT file
-9) to quit -8) Registration Information
Enter choice (0-5, -1, -3, -8 or -9) 1
Enter .BAT filename
fit
The input .BAT filename is fit.BAT
METHOD OF ANALYSIS
0) Normal fitting
1) Bayesian
2) Simulation only
3) Iterative Reweighted Least Squares
4) Simulation with random error or sensitivity analysis
5) Grid Search
-5) To perform Monte Carlo run (Only once at the start of BAT file)
-4) To perform multi-run (End of BAT file only)
-3) To run random number test subroutine
-2) To close (or open) .BAT file
-1) To finish
Enter choice (-3 to 5) 0
Where do you want the output?
0) Terminal screen
1) Disk file
Enter choice (0-2) 1
Enter .OUT filename
The output .OUT filename is fit.OUT
Enter the parameters for the run.
Parameters currently available in Boomer for model building
Enter type# for parameter 1 (-5 to 51) 1
Enter parameter name Dose
Enter Dose value 200.0
0) fixed, 1) adjustable, 2) single dependence
or 3) double dependence 0
Enter component to receive dose 1
Enter component for F-dependence ( 1 to - 1 or 0 for no dependence) 0
Input summary for Dose (type 1)
Fixed value is 200.0
Dose/initial amount added to 1
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter pharmacokinetic parameters (k10, k12, k21 and V1) for the two compartment model. Note the line name, Cp, for compartment one. These parameters are specified as adjustable and a lower and upper limit are provided. Output should be checked for a hit limit condition and changes made before the analysis is re-run.
Enter type# for parameter 2 (-5 to 51) 2
Enter parameter name k10
Enter k10 value 0.2000
0) fixed, 1) adjustable, 2) single dependence
or 3) double dependence 1
Enter lower limit 0.1000
Enter upper limit 0.4000
Enter component to receive flux 0
Enter component to lose flux 1
Input summary for k10 (type 2)
Initial value 0.2000 float between 0.1000 and 0.4000
Transfer from 1 to 0
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter type# for parameter 3 (-5 to 51) 2
Enter parameter name k12
Enter k12 value 1.300
0) fixed, 1) adjustable, 2) single dependence
or 3) double dependence 1
Enter lower limit 0.6000
Enter upper limit 2.600
Enter component to receive flux 2
Enter component to lose flux 1
Input summary for k12 (type 2)
Initial value 1.300 float between 0.6000 and 2.600
Transfer from 1 to 2
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter type# for parameter 4 (-5 to 51) 2
Enter parameter name k21
Enter k21 value 1.100
0) fixed, 1) adjustable, 2) single dependence
or 3) double dependence 1
Enter lower limit 0.5000
Enter upper limit 2.200
Enter component to receive flux 1
Enter component to lose flux 2
Input summary for k21 (type 2)
Initial value 1.100 float between 0.5000 and 2.200
Transfer from 2 to 1
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter type# for parameter 5 (-5 to 51) 18
Enter parameter name V1
Enter V1 value 12.00
0) fixed, 1) adjustable, 2) single dependence
or 3) double dependence 1
Enter lower limit 6.000
Enter upper limit 24.00
Enter data set (line) number 1
Enter line description Cp
Enter component number (0 for obs x) 1
Input summary for V1 (type 18)
Initial value 12.00 float between 6.000 and 24.00
Component 1 added to line 1
Enter 0 if happy with input, 1 if not, 2 to start over 0
Next we enter the parameters for the effect compartment. Note inclusion of the -k1e term means that the effect compartment amount/concentration doesn't influence the pharmacokinetic results. The model described here places the drug receptor for the effect/response measured in a hypothetical effect compartment.
Enter type# for parameter 6 (-5 to 51) 2
Enter parameter name k1e/V1
Enter k1e/V1 value 0.3000
Enter component to receive flux 3
Enter component to lose flux 1
Input summary for k1e/V1 (type 2)
Fixed value is 0.3000
Transfer from 1 to 3
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter type# for parameter 7 (-5 to 51) 2
Enter parameter name -k1e/V1
Enter -k1e/V1 value -0.3000
Enter component to receive flux 0
Enter component to lose flux 1
Input summary for -k1e/V1 (type 2)
Fixed value is -0.3000
Transfer from 1 to 0
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter type# for parameter 8 (-5 to 51) 2
Enter parameter name ke0
Enter ke0 value 0.8000
Enter component to receive flux 0
Enter component to lose flux 3
Input summary for ke0 (type 2)
Fixed value is 0.8000
Transfer from 3 to 0
Enter 0 if happy with input, 1 if not, 2 to start over 0
Finally we can enter the Hill equation parameters, eMax, EC50% and gamma (here S or i for slope term).
Enter type# for parameter 9 (-5 to 51) 11
Enter parameter name Effect
Enter Emax or a Effect value 100.0
Enter data set (line) number 2
Enter line description Effect
Enter component number (0 for obs x) 3
Input summary for Emax or a Effect (type 11)
Fixed value is 100.0
Component 3 added to line 2
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter Ec(50%) or b Effect value 8.000
Input summary for Ec(50%) or b Effect (type 12)
Fixed value is 8.000
Enter 0 if happy with input, 1 if not, 2 to start over 0
Enter S or i Effect value 1.500
Input summary for S or i Effect (type 13)
Fixed value is 1.500
Enter 0 if happy with input, 1 if not, 2 to start over 0
Next the numerical integration method is chosen and a description entered before entering the data. Compartmental pharmacokinetic models are typically not stiff systems so the Fehlberg RKF45 method is usually a good choice. I like to start the fitting with the Simplex method to get close and finish off fitting with the Damping Gauss-Newton method. The Simplex method start at random parameter spaces so it is useful to use this method with multiple run to avoid local minima.
Method of Numerical Integration
0) Classical 4th order Runge-Kutta
1) Runge-Kutta-Gill
2) Fehlberg RKF45
3) Adams Predictor-Corrector with DIFSUB
4) Gears method for stiff equations with PEDERV
5) Gears method without PEDERV
Enter choice (0-5) 2
Enter Relative error term for
Numerical integration (0.0001) 0.000
Enter Absolute error term for
Numerical integration (0.0001) 0.000
FITTING METHODS
0) Gauss-Newton
1) Damping Gauss-Newton
2) Marquardt
3) Simplex
4) Simplex->Damping GN
Enter Choice (0-4) 4
Enter PC for convergence (0.00001) 0.000
Enter description for this analysis: Fit to PK PD Model
Enter data from
0) Disk file 2) ...including weights
1) Keyboard 3) ...including weights
Enter Choice (0-3) 1
Enter data for Cp
Enter x-value (time) = -1 to finish data entry
X-value (time) 0.000
Y-value (concentration) 0.000
After entering values for each data set exit by entering -1 for the last X-value. Check for errors, save or not and move to the next data set. After all the data sets are entered you can choose a weighting scheme for each data set. Here I've selected 1/Obs value squared. This assumes a constant coefficient of variation for all the data. After the fitting has converged you can choose to calculate AUC values or other options such as graphs.
...
X-value (time) 36.00
Y-value (concentration) 8.200
X-value (time) -1.000
Data for Effect
DATA # Time Concentration
1 0.000 0.000
2 0.1000E-01 2.200
3 0.2000E-01 5.800
4 0.5000E-01 19.00
5 0.1000 38.00
6 0.2000 59.00
7 0.5000 79.00
8 1.000 85.00
9 2.000 87.00
10 6.000 82.00
11 12.00 68.00
12 18.00 49.00
13 24.00 28.00
14 30.00 16.00
15 36.00 8.200
Do you want to
0) Accept data
1) Correct data point
2) Delete data point
3) Insert new data point
4) Add offset to x-value
Enter choice (0-3) 0
Save Observed Data to Disk Module
0) Continue without saving
1) Save data for on disk
Enter choice 0
Weighting function entry for Cp
0) Equal weights
1) Weight by 1/Cp(i)
2) Weight by 1/Cp(i)^2
3) Weight by 1/a*Cp(i)^b
4) Weight by 1/(a + b*Cp(i)^c)
5) Weight by 1/((a+b*Cp(i)^c)*d^(tn-ti))
Data weight as a function of Cp(Obs)
Enter choice (0-5) 2
Weighting function entry for Response
0) Equal weights
1) Weight by 1/Cp(i)
2) Weight by 1/Cp(i)^2
3) Weight by 1/a*Cp(i)^b
4) Weight by 1/(a + b*Cp(i)^c)
5) Weight by 1/((a+b*Cp(i)^c)*d^(tn-ti))
Data weight as a function of Response(Obs)
...
After convergence with the Simplex method and the damping Gauss-Newton method.
...
Calculation of AUC and AUMC section
0) Exit this section
1) Cp
2) Effect
3) Receptor
4) Tissue
Enter line # for required AUC (0- 4) 0
Additional Output
Enter -> 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
Save Data x x x x x x x x
Graphs x x x x x x x x
Extra files x x x x x x x x
Sensitivity x x x x x x x x
Save calculated data files, Produce linear, semi-log and
weighted residual plots, save extra data files or perform
sensitivity or optimal sampling analysis
Enter choice (0-15) 2
The run parameters are displayed and final parameter values tabulated. The integration and fitting methods are listed. The parameter values with S.D. and C.V. % are provided. CV% less than 10% is usually a good result. Various other parameters are provided.
** FINAL OUTPUT FROM Boomer (v3.4.5) ** 5 March 2020 --- 1:06:06 pm
Title: Fit to PK PD Model
Input: From fit.BAT
Output: To fit.OUT
Data for Cp came from keyboard (or ?.BAT)
Data for Response came from keyboard (or ?.BAT)
Fitting algorithm: DAMPING-GAUSS/SIMPLEX
Weighting for Cp by 1/Cp(Obs )^2
Weighting for Response by 1/Cp(Obs )^2
Numerical integration method: 2) Fehlberg RKF45
with 3 de(s)
With relative error 0.1000E-03
With absolute error 0.1000E-03
DT = 0.1000E-02 PC = 0.1000E-04 Loops = 3
Damping = 4
** FINAL PARAMETER VALUES ***
# Name Value S.D. C.V. % Lower <-Limit-> Upper
1) k10 0.19707 0.305E-02 1.5 0.10 0.40
2) k12 1.1945 0.737E-01 6.2 0.60 2.6
3) k21 1.0382 0.479E-01 4.6 0.50 2.2
4) V1 12.225 0.178 1.5 6.0 24.
5) ke0 0.82833 0.151E-01 1.8 0.20 3.2
6) Ec(50%) or b Emax 0.25075 0.333E-02 1.3 0.10 0.50
7) S or i Emax 1.4951 0.103E-01 0.69 0.70 3.0
Final WSS = 0.871169E-02 R^2 = 0.9997 Corr. Coeff = 0.9999
AIC = -123.550 AICc = -118.550
Log likelihood = 76.5 Schwarz Criteria = -113.979
R^2 and R - jp1 0.9996 0.9998
R^2 and R - jp2 0.9996 0.9998
RMSE = 0.4835 or 1.733 % RMSE
MAE = 0.2775 ME = -0.7631E-01
Among other output information Boomer provides a summary of the model as entered. Here we can check for errors. The BAT file can be edited and re-run as needed.
Model and Parameter Definition # Name Value Type From To Dep Start Stop 1) Dose = 200.0 1 0 1 0 0 0 2) k10 = 0.1971 2 1 0 0 0 0 3) k12 = 1.194 2 1 2 0 0 0 4) k21 = 1.038 2 2 1 0 0 0 5) V1 = 12.23 18 1 1 0 0 0 6) k1e = 0.1000E-01 2 1 3 0 0 0 7) -k1e = -0.1000E-01 2 1 0 0 0 0 8) ke0 = 0.8283 2 3 0 0 0 0 9) Emax or a Emax = 100.0 11 3 2 0 0 0 10) Ec(50%) or b Emax = 0.2507 12 0 0 0 0 0 11) S or i Emax = 1.495 13 0 0 0 0 0The next table presents the observed and calculated values along with weighted residual. Check this table carefully for errors and non-random weighted residuals.
Data for Cp :-
DATA # Time Observed Calculated (Weight) Weighted residual
1 0.000 0.00000 16.3599 0.00000 -0.00000
2 0.5000E-01 15.0000 15.2841 0.666667E-01 -0.189425E-01
3 0.1000 14.0000 14.3241 0.714286E-01 -0.231484E-01
4 0.1500 14.0000 13.4669 0.714286E-01 0.380753E-01
5 0.2000 13.0000 12.7014 0.769231E-01 0.229708E-01
6 0.2500 12.0000 12.0173 0.833333E-01 -0.143957E-02
7 0.5000 9.40000 9.53657 0.106383 -0.145292E-01
8 1.000 7.20000 7.23043 0.138889 -0.422683E-02
9 2.000 5.90000 5.87982 0.169492 0.342022E-02
10 4.000 4.90000 4.86490 0.204082 0.716287E-02
11 9.000 3.20000 3.14295 0.312500 0.178287E-01
12 12.00 2.40000 2.41845 0.416667 -0.768781E-02
13 18.00 1.40000 1.43198 0.714286 -0.228405E-01
14 24.00 0.850000 0.847885 1.17647 0.248811E-02
15 30.00 0.500000 0.502038 2.00000 -0.407684E-02
16 36.00 0.300000 0.297256 3.33333 0.914802E-02
WSS data set 1 = 0.4171E-02 R^2 = 0.9986 Corr. Coeff. = 0.9993
R^2 and R - jp1 0.9996 0.9998
R^2 and R - jp2 0.9996 0.9998
RMSE = 0.1975 or 1.668 % RMSE
MAE = 0.1196 ME = -0.6927E-02
Data for Response :-
DATA # Time Observed Calculated (Weight) Weighted residual
1 0.000 0.00000 0.00000 0.00000 0.00000
2 0.1000E-01 2.20000 2.19383 0.454545 0.280424E-02
3 0.2000E-01 5.80000 5.85533 0.172414 -0.953921E-02
4 0.5000E-01 19.0000 18.8989 0.526316E-01 0.531980E-02
5 0.1000 38.0000 37.7362 0.263158E-01 0.694275E-02
6 0.2000 59.0000 59.3721 0.169492E-01 -0.630757E-02
7 0.5000 79.0000 78.8822 0.126582E-01 0.149073E-02
8 1.000 85.0000 85.1443 0.117647E-01 -0.169785E-02
9 2.000 87.0000 86.6743 0.114943E-01 0.374323E-02
10 6.000 82.0000 81.3412 0.121951E-01 0.803403E-02
11 12.00 68.0000 66.6963 0.147059E-01 0.191723E-01
12 18.00 49.0000 47.7762 0.204082E-01 0.249752E-01
13 24.00 28.0000 29.4720 0.357143E-01 -0.525708E-01
14 30.00 16.0000 16.0282 0.625000E-01 -0.176525E-02
15 36.00 8.20000 8.01970 0.121951 0.219884E-01
WSS data set 2 = 0.4540E-02 R^2 = 0.9996 Corr. Coeff. = 0.9998
R^2 and R - jp1 0.9995 0.9998
R^2 and R - jp2 0.9995 0.9998
RMSE = 0.2045 or 1.709 % RMSE
MAE = 0.1279 ME = -0.7225E-02
Linear and semi-log printer-style plots can be output for a quick review of the results. Note that you can select (above) to save the calculated data for plotting with more advanced graphing software. Below some of these data were plotted using a spreadsheet program. Here we see observed and calculated data on linear and semi-log plots. Check for outliers which may be caused by data entry error, model or weight misspecification. We also have plots of the weighted residuals. Patterns might also indicate data entry error, model or weight misspecification.
Plots of observed (*) and calculated values (+)
versus time for Cp. Superimposed points (X)
16.36 Linear 16.36 Semi-log
|+ |+
| |+
|+ |X
|* |X
|+ |
|* |X
|+ |
|X | X
|X | X
| |
| | X
| |
| |
|X | X
| |
| | X
| |
| X |
| |
| | X
| X |
| |
| X |
| | X
| |
| X |
| X |
| | X
| X |
| X |
|* X X | X
|_____________________________________ |*____________________________________
0.000 0.2973
0 <--> 36. 0 <--> 36.
Plot of Std Wtd Residuals (X) Plot of Std Wtd Residuals (X)
versus time for Cp versus log(calc Cp(i)) for Cp
1.998 1.998
|X | X
| |
| |
| |
|X | X
| |
| X | X
| |
| X |X
| X | X
0X=X=====================X============ 0=========X================X=========X
|X | X
| X X X | X X X
| |
|X | X
|X | X
|X X | X X
| |
-1.215 -1.215
0.0 <--> 36. 0.30 <--> 16.
Plots of observed (*) and calculated values (+)
versus time for Response. Superimposed points (X)
87.00 Linear 87.00 Semi-log
| * | *
| X+ |XX+ X
| X |
|X | X
| |X
| | X
| |
| * |X
| + |
| | +
|X | *
| |
| |
| |X
| X | X
| |
| |
|X |
| |
| |
| + | X
| * |
| |+
| |*
|X |
| X |
| |
| |
|X X |
| |
|X |X
|_____________________________________ |X____________________________________
0.000 2.194
0 <--> 36. 0 <--> 36.
Plot of Std Wtd Residuals (X) Plot of Std Wtd Residuals (X)
versus time for Response versus log(calc Cp(i)) for Response
1.311 1.311
| X | X
| X | X
| X | X
| |
|X X | X X
|X X |X X XX
0XX============================X====== 0===================X===============X=
|X | X
|X | X
| |
| |
| |
| |
| |
| |
| |
| X | X
| |
-2.759 -2.759
0.0 <--> 36. 2.2 <--> 87.
Plot of Std Wtd Residuals (X) Plot of Std Wtd Residuals (X)
versus time for all data versus log(calc Cp(i)) for all data
1.998 1.998
|X | X
| |
| |
|X X X | X X X
| X X | X X
| |
|X X X X |X X X X X
0X=X=====================X============ 0======X=====X=====X====X=X=========XX
|XX X | X X X X X
|X X | X X
|X | X
|X X | X X
| |
| |
| |
| |
| X | X
| |
-2.759 -2.759
0.0 <--> 36. 0.30 <--> 87.
It might interesting to view some of the plots that could be prepared from this output.
Concentration versus Time. Observed and Calculated Data.
The BAT and OUT files can be dowloaded.
The BAT file can be saved as sim.BAT and run with Boomer. With the current version of Boomer for Windows the lineEnter component for F-dependence ( 1 to - 1 or 0 for no dependence) 0will need to be removed.