- Email: [email protected]
Simulation of complex pharmacokinetic models in Microsoft EXCEL
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Simulation of complex pharmacokinetic models in Microsoft EXCEL ¨ ¨ Ingolf Meineke ∗ , Jurgen Brockmoller ¨ ¨ Department of Clinical Pharmacology, University of Gottingen, Robert-Koch-Street 40, 37075 Gottingen, Germany
a r t i c l e
i n f o
a b s t r a c t
Article history:
With the arrival of powerful personal computers in the office numerical methods are acces-
Received 14 March 2007
sible to everybody. Simulation of complex processes therefore has become an indispensible
Received in revised form
tool in research and education. In this paper Microsoft EXCEL is used as a platform for a
12 September 2007
universal differential equation solver. The software is designed as an add-in aiming at a
Accepted 12 September 2007
minimum of required user input to perform a given task. Four examples are included to demonstrate both, the simplicity of use and the versatility of possible applications. While
Keywords:
the layout of the program is admittedly geared to the needs of pharmacokineticists, it can
Differential equations
be used in any field where sets of differential equations are involved. The software package
Pharmacokinetics
is available upon request. © 2007 Elsevier Ireland Ltd. All rights reserved.
VBA Simulation
1.
Introduction
Simulation of physical, chemical and biological processes can provide substantial insight into underlying mechanisms in research and academic teaching. In particular, the simulation of concentration–time or effect vs. time profiles is a frequently recurring task in drug research, toxicology or biochemistry. An increasing number of programs is presently available, both, commercially distributed and as freeware, which allow to perform such tasks at any level of sophistication [1–10]. The references are provided to illustrate the point under discussion, but do not give a complete collection of such programs. For occasional users, or for those who want to obtain an almost instantaneous answer, many software packages cannot be intuitively applied to a given problem or require so much overhead (e.g. writing of scripts, preparation of input data files) that frequent use is not encouraged. Moreover, only the more basic situations like monoexponential radioactive decay or simple pharmacokinetic models can be easily described in algebraic form or be solved by linearization. To
∗
apply the method of Laplace transforms [11] for the solution of complex systems used in pharmacokinetic analysis is tedious and not efficient. If saturable processes, metabolic or transport, have to be taken into account then this way is barred, too. There is no explicit solution to the simple Michaelis–Menten equation [12]. Coupled differential equations on the other hand represent a universal approach to describe complex models in a direct and efficacious manner. A universal solver is presented in this work that was developed for non-standard problems as a tool to provide fast answers to what–if questions. This approach gives priority to the speed of preparation (writing of the model) rather than the speed of execution (calculating the concentrations). Nonstandard problems are for instance catenary models, models with non-central elimination, multicompartment systems with parent compound, several metabolites and interconversion. It was our intention to develop a set of routines which allow the direct translation of a typical model blueprint with its compartments and transfer routes into an equivalent mathe-
Corresponding author. Tel.: +49 551 395796; fax: +49 551 396970. E-mail address: [email protected] (I. Meineke). 0169-2607/$ – see front matter © 2007 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2007.09.007
240
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
matical form with the absolute minimum of user input. With such a tool we furthermore hope to give an incentive for more intensive use of simulation especially outside those fields where mathematical modelling is traditionally found. Microsoft® EXCEL under Windows XP was chosen as the programming environment since this software is a de facto standard for spreadsheet applications with widespread availability and offers many built-in mathematical and graphical routines which can be called in user defined functions. It is gratefully acknowledged that important concepts used in this work have been borrowed from the NONMEM project which is an outstanding exemplary.
2.
Computational methods
The differential equation solver is implemented in two files named pkengine 2.xla and transmission.dll. Pkengine 2.xla defines a simple interface which is called by the user in order to make the model known to the solver and then to perform the calculations. The subroutine calls (Fig. 1) are made from a EXCEL workbook designed for this purpose. In this model workbook the user provides the necessary information to specify a model as a set of differential equations. It is further required that a reference to the pk-engine is set in the model workbook. Since circle references are not allowed in EXCEL for obvious reasons, the call from the pk-engine to the subroutine model() is made with the aid of the support DLL (Fig. 1). The Runge-Kutta method (fourth order with variable step size) is implemented to numerically integrate the differential equations [13]. The main loop of the pk-engine advances from time zero to the final time in steps defined by the time interval and keeps track of dosing events. Dosing events are handled with priority and must be specified in a separate list. Especially with zero order infusions additional calls to the solver can be required within time intervals. Results are collected into a variable array which is written to the respective worksheet with a simple assignment at the end of the calculations. EXCEL chart objects are used to automatically provide graphic output. User defined additional calculations can be implemented in a
Fig. 1 – Program structure. The model independent routines are encapsulated in the pk-engine. The user defined model definitions remain in the current workbook.
Fig. 2 – Minimal model. The example shows a one-compartment model with bolus drug administration. For each compartment the differential equation must be defined and the variables used must be declared.
subroutine that is called at each time point after the call to the model subroutine. The additional subroutine is present as a stub in the default situation.
3.
Program description
A model workbook (see examples below) consists of four mandatory worksheets. Additional worksheets can de added at the users own discretion. Two of the mandatory worksheets (Model, Dose History) require input from the user, the remaining two (Concentrations, Diagrams) receive the output from the program. In the model-sheet the set of differential equations is defined within the body of a subroutine called model(). The declaration of this subroutine should not be altered. Within the body any user defined variables and expressions are admissible. In addition to the differential equations the corresponding transfer constants must be specified by the user as illustrated in Fig. 2 in the model specific mandatory information section. No assumptions are made whatsoever in terms of dimensions of the variables. In a pharmacokinetic context, however, ti would typically denote time, ai() amount and dadt() amount per time. An easy and effective way to implement a new model is a draft sketch of the intended structure with consecutively numbered compartments with a list of transfer constants beneath the sketch. The constants can be indexed to indicate transfer direction, e.g. k12 denotes transfer from entity 1 to 2. While the size of user input will vary with the complexity of the model to be simulated, only four global constants must be adapted by the user. These are arranged in the mandatory general information section and denote the total number of compartments (NN), the time interval (tint) between two consecutive observations and the total time (ttotal). More information with respect to advanced features is provided in Section 4 of his manuscript.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 3 – Dose history. In the dose history sheet arbitrary complex dose schemes can be defined. For each dose event amount, route of administration and the input compartment can be chosen. For oral doses an additional dose compartment must be provided in the model.
Certain constants (see Appendix A for details) can be defined in the model specific facultative information section. They may also be assigned empty strings (“”) by the user. In the dose history sheet one single line of input is required which must specify the start of the time profile at time t = 0. The dose entered at this time can be even zero, which would of course result in a set of flat lines for all compartments. In addition to the dose the rate can be entered for zero order infusions together with the compartment for which the dose is intended. However, more complex dosing histories will be investigated as a rule and the number of entries is only limited by the number of rows in an EXCEL sheet (Fig. 3). The few lines in Fig. 3 state that in a model with at least two compartments a bolus dose (dose rate 0) is introduced into compartment 2 at time zero. At the same time a zero order infusion is started into the same compartment which will last 50 time units. After 30 time units a second infusion is started which ends 10 time units later, i.e. before the first infusion is complete. Finally at time 60 a bolus dose is placed into compartment 1. Depending on the model this could be an oral dose. From the dose history sheet the model is compiled and calculated. For these two steps command buttons have been provided in the worksheet. They trigger calls to pkengine.xla, where the routines needed for the compilation of the model given in the modelsheet and the calculation of the resulting concentrations are located (Fig. 1). The concentrations sheet is simply a container which receives the concentration–time information for all compartments. The compartment names defined in the model specific facultative section are used as column headers to ease association of the results. It is important to note that the concentrations sheet will be rebuilt each time the model is recalculated. Any input from the user will then be overwritten. The diagrams sheet finally presents fundamental plots of all concentration–time profiles for the convenience of the user provided the plot option for a compartment is set. While the model building steps only require the absolute minimum input from the user, the module offers additional optional functionality. As explained in Section 4 in full detail, in the facultative information section of the model-sheet names for kinetic and auxiliary result columns can be specified and it is possible to select the profiles to be plotted. The default for plot is 1, i.e. the corresponding profile is shown. Moreover, the entries in the string variable vi are used for scaling purposes. It is important to note that the
241
values entered here are used internally to obtain concentrations, which are then tabulated in the concentrations sheet and handed to the auxiliary subroutine as an [in] parameter. Therefore, if the vi are set to 1 (the default) concentrations and amounts as reported by the model subroutine are numerically identical. If the volumes of distribution are used for scaling in vi, true concentrations (amount/volume) are listed. Otherwise the concentrations calculated will be scaled by volume/vi. Some intelligence has also been placed into the dose history sheet. This code is used to call the pk-engine routines as shown in Fig. 1 The user is advised to create new models by saving an existing model workbook under a different name and then modifying this new instance. This way the code in the dose history sheet and the pivotal reference to pkengine.xla are preserved.
4.
Example program runs
The first example (minimal model.xls) (Fig. 2) describes the basic requirements for a one-compartment body model. Hence, in the subroutine in the model-sheet one differential equation describes the transfer from the single compartments with one model specific constant. Three general constants are also needed. Of these NN (the number of differential equations in model()) is consequently set to 1, the total time (ttotal) is assumed to be 10 time units with a time interval (tint) of 1 units, resulting in 11 observations. The history sheet contains a single dose event of 100 units at time zero. No further input is needed for this bare-bones example. The unspectacular result diagram shows the well-known exponential decay. The main purpose of this example is to demonstrate the concept of modelling with the pk-engine and to point out how little overhead the program requires. The second example (diazepam model.xls) is based on ¨ a diazepam clinical study [Brockmoller, unpublished data] which produced concentration–time data for diazepam and three of its metabolites, namely desmethyldiazepam, temazepam and oxazepam. The assumed relationship between the compounds is given in Fig. 4 which was used as the blueprint for the model definition. The model was implemented as follows (Fig. 5). In the first step the model() subroutine was coded. For this purpose the transfer into and from each compartment was described, beginning with the dose compartment (#1) to the oxazepam compartment (#6). Two additional compartments were then defined, one effect compartment (#7) for diazepam and a term to calculate the oxazepam AUC (#8). The total number of differential equations in NN was therefore 8. In the second step the auxiliary() subroutine was coded. In the example a simple EMAX effect model was implemented with two parameters EMAX and C50, the effect was linked to the concentration in the diazepam compartment (#2) and for comparison also to the effect compartment (#7). The number of auxiliary equations is saved in NA. It is noteworthy that EMAX and C50 are declared without the public keyword. They are local to auxiliary(). The subroutine makes use of the concentration–time information provided in ti and ci() and returns its calculations in results(). In the example a
242
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 4 – Diazepam metabolic pathway. All rate transfer constants are h−1 . All volumes are L. The flowchart can be directly translated into differential equations as shown in the implementation (see Fig. 5).
pharmacokinetic–pharmacodynamic link via an effect compartment is thus compared with a direct link. The model specific additional constants section was then completed, defining compartment names, scaling factors and, via the plot option, whether a compartment is shown in the model plot. Finally, the transfer constants used in the model() subroutine were defined in the model specific mandatory constants section. Model parameterization can be made in term of clearances and volumes or transfer constants and volumes or any other combination. Here clearance and volume are the primary parameters and the transfer constants calculated in the model-worksheet. Without the advanced features the model would appear more lean, but probably less informative. The results for the diazepam model are presented in Fig. 6. It is important to note that the dose compartment is not plotted and that no scaling was performed. Therefore the vertical axis in the model graph is in amount rather than concentration units and the diagram shows how the dose is distributed among compartments. If this is not desired, individual volumes for each compartment can be provided in the variable vi. The second diagram (plot AUXILIARY) in the diazepam example compares the effect-time course of a direct effect (located in the central compartment) with that of an effect mediated through an effect compartment. This approach can be pushed further for instance by asking what effect-time course would result if the contribution of the metabolites to the effect is taken into consideration. A wealth of scenarios is imaginable and can be easily implemented.
The third example (pb-pk-model.xls) was taken from literature [14]. The example shows the distribution and elimination of styrene in the rat during and after the exposure to 80 ppm styrene in the air. In contrast to mammillary models physiological based models attempt to mimic the mass transport to important organs via the circulation (Fig. 7) and include additional a priori information like distribution factors into the model. Moreover, saturable elimination cannot be algebraically described in a closed form [12]. For these reason it was deemed instructive to include such a model in the set of examples. The model implementation is not depicted, since all the relevant information is available in the structural model. The results plot has some additional features of interest (Fig. 8). METAB and PULMO are amounts of styrene metabolized and eliminated via the pulmonal route, respectively. The contents of the other compartments are given as concentrations. FAT/10 indicates that FAT has been downscaled by a factor of 10. Otherwise this compartment would dominate the plot and the other graphs would be squeezed against the x-axis. The last example (lotka-volterra.xls) contains a non-linear multiplicative model (Fig. 9). The oscillating prey–predator populations described by the model are particularly sensitive to changes in the start conditions (Fig. 10). The differential equations used in the model represent only a coarse approximation and more complex and better models can be easily devised.
4.1.
Specifications
The pk-engine package was developed and tested in EXCEL 2003 under Microsoft Windows XP service pack 2. The models were also successfully run in EXCEL 2000. Other EXCEL versions or Windows versions were not tested. The results produced by the software were numerically identical with those of corresponding model simulations in NONMEM.
5.
Limitations
The Runge-Kutta algorithm has the reputation to perform badly with stiff differential equations. We have tested models with transfer constants ranging from 0.001 to 10 and have so far observed no restrictions in terms of accuracy. Graphic representations of the compartments in the model by default use only one ordinate. Therefore if large concentration or amount differences are encountered the graph will be scaled to the highest value and smaller value will occasionally be difficult to discern. There is however nothing that prevents the user to produce better adapted and better looking graphs from the concentration–time data on a case-to-case basis. One simple remedy would be the use of the secondary ordinate EXCEL offers via axgrp (see Appendix A). Additional worksheets are an effective solution if the user wants to keep all columns or selected columns from a run in order to compare results obtained under different conditions or between models. Data can either be copied to the new sheet or the concentrations sheet is renamed. In that case a new concentrations worksheet must be generated or a range error
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
243
Fig. 5 – Diazepam model implementation Following the principles explained with the minimal model the differential equations are an equivalent representation to the flowchart (see Fig. 4).
Fig. 6 – Diazepam model results The plot titled Model shows the results calculated in the model() subroutine. Calculations from the auxiliary() subroutine are presented in a separate plot called auxiliary.
244
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 8 – PB-PK model results. the appropriate alterations in the model-sheet and rerun the model.
6.
Fig. 7 – Physiological based pharmacokinetic model: structural model QF , QL , QR , QS partial flow rates and Qp air flow; Ca and Cv arterial and venous concentrations; respectively; Cin styrene concentration in the inhaled air.
will occur. The pk-engine in its present form does not handle such errors and code execution will be stopped by the Visual Basic (VB) run time environment. Similarly, syntax errors and the like in the module code or elsewhere will be detected by the VB compiler and the VB editor will appear with the offending code highlighted. This does not mean that the users need to learn VB syntax or are required to understand details about the workings of the VB editor. It is recommended to take a note of the line containing the error and then to simply close the VB editor. This will abort the current run and the user can make
Availability
The pk-engine package is available as a compressed archive via electronic mail from the authors upon request. The request should therefore include an e-mail address. The package contains the pk-engine, the support DLL and the four examples presented above. More examples have been developed by the authors including experimental alternative absorption models. We will gladly support scientists who wish to implement a particular model. Likewise we are interested to receive any bug reports describing unexpected behaviour of the software. The software is mostly harmless insofar as it contains no malign code and performs only the tasks of solving ordinary differential equations and plotting the results. Since EXCEL is also an automation server, a variety of scenarios is imaginable that run the pk-engine package from within other applications like Microsoft Powerpoint, Microsoft Word or literally any other automation client. Web based projects as front-ends could further add to user-friendliness.
Fig. 9 – Lotka-Volterra model: example of a non-linear multiplicative model. Small changes in the start values (PREY = 5, PREDATOR = 10) have pronounced effect.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 10 – Lotka-Volterra type oscillations. Stable attractors can only be obtained with suitable combinations of start values and model parameter values. We are looking forward to receiving the suggestions of our readers.
Conflicts of interest
245
comp names
Array of user defined compartment names, facultative
dadt
Array of 1..NN differential equations
model
Subroutine defining the set of differential equations
NA
Constant defining the number of additional results in auxiliary
NN
Constant defining the number of differential equations
plot
String of flags (0,1) determining whether a profile is plotted (1) or not (0), facultative.
results
Array of 1..NA user defined equations, output from auxiliary
ti
Constant stating the time interval between observations
ttotal
Constant stating the total time for profiling
vi
String of compartmental distribution volumes (comma separated), facultative
The authors report no conflicts of interest.
references
Appendix A Quick start: 1. Extract the contents of the pk-engine archive to a directory of your choice. 2. With the helper application setref.exe check the presence of a reference to PKENGINE in the file minimum model.xls or add this reference. 3. Start EXCEL and open the file minimum model.xls. 4. From the Dose history sheet compile and calculate the model. 5. Change the dose and recalculate. 6. Change k10 on the model-sheet, then recompile and recalculate. 7. Close minimal model.xls and repeat from step2 with diazepam model.xls or pb-pk-model.xls. Reserved variables: names and functions ai Array 1..NN describing the compartment amounts auxiliary
Optional subroutine with parameters ti, ci, results
aux names
Array of user defined names for optional data columns, facultative
axgrp
String of flags (1,2) stating the axis (primary or secondary) for the plot, facultative
ci
Array 1..NN describing compartment concentrations, input to auxiliary
[1] PK Solutions [online], Available from the World Wide Web: http://www.summitpk.com/pksolutions/pksolutions.htm. [2] GastroPlus [online], Available from the World Wide Web: http://www.simulations-plus.com. [3] SAAM II [online], Available from the World Wide Web: http://depts.washington.edu/saam2/. [4] PCCAL [online], Available from the World Wide Web: http://www.coacs.com/PCCAL/. [5] ADAPT II [online], Available from the World Wide Web: http://bmsr.usc.edu/. [6] MODELMAKER [online], Available from the World Wide Web: http://www.modelkinetix.com/. [7] G.M. Hood, PopTools version 2.6.9, Available from the World Wide Web: http://www.cse.sciro.au/poptools. [8] WinNonLin Professional Edition 2.1, Pharsight Corporation, Palo Alto, California. [9] A. Boeckmann, L. Sheiner, S. Beal, NONMEM User’s Guide, San Francisco, University of California, San Francisco, 1994. [10] R.F. Muzic, S. Cornelius, COMKAT: compartment model kinetic analysis tool, J. Nucl. Med. 42 (2001) 636–644. [11] M. Gibaldi, D. Perrier, Pharmacokinetics, Marcel Dekker, New York, 1982, p. 419. [12] S.L. Beal, Computation of the explicit solution to the Michaelis–Menten equation, J. Pharmacok. Biopharm. 11 (1983) 641–657. [13] R.G. Davies, Computer Programming In Quantitative Biology, Academic Press, London/New York, 1971, p. 388. [14] S. Haddad, M. Pelekis, K. Krishnan, A methodology for solving physiologically based pharmacokinetic models without the use of simulation software, Toxicol. Lett. 85 (1996) 113–126.
journal homepage: www.intl.elsevierhealth.com/journals/cmpb
Simulation of complex pharmacokinetic models in Microsoft EXCEL ¨ ¨ Ingolf Meineke ∗ , Jurgen Brockmoller ¨ ¨ Department of Clinical Pharmacology, University of Gottingen, Robert-Koch-Street 40, 37075 Gottingen, Germany
a r t i c l e
i n f o
a b s t r a c t
Article history:
With the arrival of powerful personal computers in the office numerical methods are acces-
Received 14 March 2007
sible to everybody. Simulation of complex processes therefore has become an indispensible
Received in revised form
tool in research and education. In this paper Microsoft EXCEL is used as a platform for a
12 September 2007
universal differential equation solver. The software is designed as an add-in aiming at a
Accepted 12 September 2007
minimum of required user input to perform a given task. Four examples are included to demonstrate both, the simplicity of use and the versatility of possible applications. While
Keywords:
the layout of the program is admittedly geared to the needs of pharmacokineticists, it can
Differential equations
be used in any field where sets of differential equations are involved. The software package
Pharmacokinetics
is available upon request. © 2007 Elsevier Ireland Ltd. All rights reserved.
VBA Simulation
1.
Introduction
Simulation of physical, chemical and biological processes can provide substantial insight into underlying mechanisms in research and academic teaching. In particular, the simulation of concentration–time or effect vs. time profiles is a frequently recurring task in drug research, toxicology or biochemistry. An increasing number of programs is presently available, both, commercially distributed and as freeware, which allow to perform such tasks at any level of sophistication [1–10]. The references are provided to illustrate the point under discussion, but do not give a complete collection of such programs. For occasional users, or for those who want to obtain an almost instantaneous answer, many software packages cannot be intuitively applied to a given problem or require so much overhead (e.g. writing of scripts, preparation of input data files) that frequent use is not encouraged. Moreover, only the more basic situations like monoexponential radioactive decay or simple pharmacokinetic models can be easily described in algebraic form or be solved by linearization. To
∗
apply the method of Laplace transforms [11] for the solution of complex systems used in pharmacokinetic analysis is tedious and not efficient. If saturable processes, metabolic or transport, have to be taken into account then this way is barred, too. There is no explicit solution to the simple Michaelis–Menten equation [12]. Coupled differential equations on the other hand represent a universal approach to describe complex models in a direct and efficacious manner. A universal solver is presented in this work that was developed for non-standard problems as a tool to provide fast answers to what–if questions. This approach gives priority to the speed of preparation (writing of the model) rather than the speed of execution (calculating the concentrations). Nonstandard problems are for instance catenary models, models with non-central elimination, multicompartment systems with parent compound, several metabolites and interconversion. It was our intention to develop a set of routines which allow the direct translation of a typical model blueprint with its compartments and transfer routes into an equivalent mathe-
Corresponding author. Tel.: +49 551 395796; fax: +49 551 396970. E-mail address: [email protected] (I. Meineke). 0169-2607/$ – see front matter © 2007 Elsevier Ireland Ltd. All rights reserved. doi:10.1016/j.cmpb.2007.09.007
240
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
matical form with the absolute minimum of user input. With such a tool we furthermore hope to give an incentive for more intensive use of simulation especially outside those fields where mathematical modelling is traditionally found. Microsoft® EXCEL under Windows XP was chosen as the programming environment since this software is a de facto standard for spreadsheet applications with widespread availability and offers many built-in mathematical and graphical routines which can be called in user defined functions. It is gratefully acknowledged that important concepts used in this work have been borrowed from the NONMEM project which is an outstanding exemplary.
2.
Computational methods
The differential equation solver is implemented in two files named pkengine 2.xla and transmission.dll. Pkengine 2.xla defines a simple interface which is called by the user in order to make the model known to the solver and then to perform the calculations. The subroutine calls (Fig. 1) are made from a EXCEL workbook designed for this purpose. In this model workbook the user provides the necessary information to specify a model as a set of differential equations. It is further required that a reference to the pk-engine is set in the model workbook. Since circle references are not allowed in EXCEL for obvious reasons, the call from the pk-engine to the subroutine model() is made with the aid of the support DLL (Fig. 1). The Runge-Kutta method (fourth order with variable step size) is implemented to numerically integrate the differential equations [13]. The main loop of the pk-engine advances from time zero to the final time in steps defined by the time interval and keeps track of dosing events. Dosing events are handled with priority and must be specified in a separate list. Especially with zero order infusions additional calls to the solver can be required within time intervals. Results are collected into a variable array which is written to the respective worksheet with a simple assignment at the end of the calculations. EXCEL chart objects are used to automatically provide graphic output. User defined additional calculations can be implemented in a
Fig. 1 – Program structure. The model independent routines are encapsulated in the pk-engine. The user defined model definitions remain in the current workbook.
Fig. 2 – Minimal model. The example shows a one-compartment model with bolus drug administration. For each compartment the differential equation must be defined and the variables used must be declared.
subroutine that is called at each time point after the call to the model subroutine. The additional subroutine is present as a stub in the default situation.
3.
Program description
A model workbook (see examples below) consists of four mandatory worksheets. Additional worksheets can de added at the users own discretion. Two of the mandatory worksheets (Model, Dose History) require input from the user, the remaining two (Concentrations, Diagrams) receive the output from the program. In the model-sheet the set of differential equations is defined within the body of a subroutine called model(). The declaration of this subroutine should not be altered. Within the body any user defined variables and expressions are admissible. In addition to the differential equations the corresponding transfer constants must be specified by the user as illustrated in Fig. 2 in the model specific mandatory information section. No assumptions are made whatsoever in terms of dimensions of the variables. In a pharmacokinetic context, however, ti would typically denote time, ai() amount and dadt() amount per time. An easy and effective way to implement a new model is a draft sketch of the intended structure with consecutively numbered compartments with a list of transfer constants beneath the sketch. The constants can be indexed to indicate transfer direction, e.g. k12 denotes transfer from entity 1 to 2. While the size of user input will vary with the complexity of the model to be simulated, only four global constants must be adapted by the user. These are arranged in the mandatory general information section and denote the total number of compartments (NN), the time interval (tint) between two consecutive observations and the total time (ttotal). More information with respect to advanced features is provided in Section 4 of his manuscript.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 3 – Dose history. In the dose history sheet arbitrary complex dose schemes can be defined. For each dose event amount, route of administration and the input compartment can be chosen. For oral doses an additional dose compartment must be provided in the model.
Certain constants (see Appendix A for details) can be defined in the model specific facultative information section. They may also be assigned empty strings (“”) by the user. In the dose history sheet one single line of input is required which must specify the start of the time profile at time t = 0. The dose entered at this time can be even zero, which would of course result in a set of flat lines for all compartments. In addition to the dose the rate can be entered for zero order infusions together with the compartment for which the dose is intended. However, more complex dosing histories will be investigated as a rule and the number of entries is only limited by the number of rows in an EXCEL sheet (Fig. 3). The few lines in Fig. 3 state that in a model with at least two compartments a bolus dose (dose rate 0) is introduced into compartment 2 at time zero. At the same time a zero order infusion is started into the same compartment which will last 50 time units. After 30 time units a second infusion is started which ends 10 time units later, i.e. before the first infusion is complete. Finally at time 60 a bolus dose is placed into compartment 1. Depending on the model this could be an oral dose. From the dose history sheet the model is compiled and calculated. For these two steps command buttons have been provided in the worksheet. They trigger calls to pkengine.xla, where the routines needed for the compilation of the model given in the modelsheet and the calculation of the resulting concentrations are located (Fig. 1). The concentrations sheet is simply a container which receives the concentration–time information for all compartments. The compartment names defined in the model specific facultative section are used as column headers to ease association of the results. It is important to note that the concentrations sheet will be rebuilt each time the model is recalculated. Any input from the user will then be overwritten. The diagrams sheet finally presents fundamental plots of all concentration–time profiles for the convenience of the user provided the plot option for a compartment is set. While the model building steps only require the absolute minimum input from the user, the module offers additional optional functionality. As explained in Section 4 in full detail, in the facultative information section of the model-sheet names for kinetic and auxiliary result columns can be specified and it is possible to select the profiles to be plotted. The default for plot is 1, i.e. the corresponding profile is shown. Moreover, the entries in the string variable vi are used for scaling purposes. It is important to note that the
241
values entered here are used internally to obtain concentrations, which are then tabulated in the concentrations sheet and handed to the auxiliary subroutine as an [in] parameter. Therefore, if the vi are set to 1 (the default) concentrations and amounts as reported by the model subroutine are numerically identical. If the volumes of distribution are used for scaling in vi, true concentrations (amount/volume) are listed. Otherwise the concentrations calculated will be scaled by volume/vi. Some intelligence has also been placed into the dose history sheet. This code is used to call the pk-engine routines as shown in Fig. 1 The user is advised to create new models by saving an existing model workbook under a different name and then modifying this new instance. This way the code in the dose history sheet and the pivotal reference to pkengine.xla are preserved.
4.
Example program runs
The first example (minimal model.xls) (Fig. 2) describes the basic requirements for a one-compartment body model. Hence, in the subroutine in the model-sheet one differential equation describes the transfer from the single compartments with one model specific constant. Three general constants are also needed. Of these NN (the number of differential equations in model()) is consequently set to 1, the total time (ttotal) is assumed to be 10 time units with a time interval (tint) of 1 units, resulting in 11 observations. The history sheet contains a single dose event of 100 units at time zero. No further input is needed for this bare-bones example. The unspectacular result diagram shows the well-known exponential decay. The main purpose of this example is to demonstrate the concept of modelling with the pk-engine and to point out how little overhead the program requires. The second example (diazepam model.xls) is based on ¨ a diazepam clinical study [Brockmoller, unpublished data] which produced concentration–time data for diazepam and three of its metabolites, namely desmethyldiazepam, temazepam and oxazepam. The assumed relationship between the compounds is given in Fig. 4 which was used as the blueprint for the model definition. The model was implemented as follows (Fig. 5). In the first step the model() subroutine was coded. For this purpose the transfer into and from each compartment was described, beginning with the dose compartment (#1) to the oxazepam compartment (#6). Two additional compartments were then defined, one effect compartment (#7) for diazepam and a term to calculate the oxazepam AUC (#8). The total number of differential equations in NN was therefore 8. In the second step the auxiliary() subroutine was coded. In the example a simple EMAX effect model was implemented with two parameters EMAX and C50, the effect was linked to the concentration in the diazepam compartment (#2) and for comparison also to the effect compartment (#7). The number of auxiliary equations is saved in NA. It is noteworthy that EMAX and C50 are declared without the public keyword. They are local to auxiliary(). The subroutine makes use of the concentration–time information provided in ti and ci() and returns its calculations in results(). In the example a
242
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 4 – Diazepam metabolic pathway. All rate transfer constants are h−1 . All volumes are L. The flowchart can be directly translated into differential equations as shown in the implementation (see Fig. 5).
pharmacokinetic–pharmacodynamic link via an effect compartment is thus compared with a direct link. The model specific additional constants section was then completed, defining compartment names, scaling factors and, via the plot option, whether a compartment is shown in the model plot. Finally, the transfer constants used in the model() subroutine were defined in the model specific mandatory constants section. Model parameterization can be made in term of clearances and volumes or transfer constants and volumes or any other combination. Here clearance and volume are the primary parameters and the transfer constants calculated in the model-worksheet. Without the advanced features the model would appear more lean, but probably less informative. The results for the diazepam model are presented in Fig. 6. It is important to note that the dose compartment is not plotted and that no scaling was performed. Therefore the vertical axis in the model graph is in amount rather than concentration units and the diagram shows how the dose is distributed among compartments. If this is not desired, individual volumes for each compartment can be provided in the variable vi. The second diagram (plot AUXILIARY) in the diazepam example compares the effect-time course of a direct effect (located in the central compartment) with that of an effect mediated through an effect compartment. This approach can be pushed further for instance by asking what effect-time course would result if the contribution of the metabolites to the effect is taken into consideration. A wealth of scenarios is imaginable and can be easily implemented.
The third example (pb-pk-model.xls) was taken from literature [14]. The example shows the distribution and elimination of styrene in the rat during and after the exposure to 80 ppm styrene in the air. In contrast to mammillary models physiological based models attempt to mimic the mass transport to important organs via the circulation (Fig. 7) and include additional a priori information like distribution factors into the model. Moreover, saturable elimination cannot be algebraically described in a closed form [12]. For these reason it was deemed instructive to include such a model in the set of examples. The model implementation is not depicted, since all the relevant information is available in the structural model. The results plot has some additional features of interest (Fig. 8). METAB and PULMO are amounts of styrene metabolized and eliminated via the pulmonal route, respectively. The contents of the other compartments are given as concentrations. FAT/10 indicates that FAT has been downscaled by a factor of 10. Otherwise this compartment would dominate the plot and the other graphs would be squeezed against the x-axis. The last example (lotka-volterra.xls) contains a non-linear multiplicative model (Fig. 9). The oscillating prey–predator populations described by the model are particularly sensitive to changes in the start conditions (Fig. 10). The differential equations used in the model represent only a coarse approximation and more complex and better models can be easily devised.
4.1.
Specifications
The pk-engine package was developed and tested in EXCEL 2003 under Microsoft Windows XP service pack 2. The models were also successfully run in EXCEL 2000. Other EXCEL versions or Windows versions were not tested. The results produced by the software were numerically identical with those of corresponding model simulations in NONMEM.
5.
Limitations
The Runge-Kutta algorithm has the reputation to perform badly with stiff differential equations. We have tested models with transfer constants ranging from 0.001 to 10 and have so far observed no restrictions in terms of accuracy. Graphic representations of the compartments in the model by default use only one ordinate. Therefore if large concentration or amount differences are encountered the graph will be scaled to the highest value and smaller value will occasionally be difficult to discern. There is however nothing that prevents the user to produce better adapted and better looking graphs from the concentration–time data on a case-to-case basis. One simple remedy would be the use of the secondary ordinate EXCEL offers via axgrp (see Appendix A). Additional worksheets are an effective solution if the user wants to keep all columns or selected columns from a run in order to compare results obtained under different conditions or between models. Data can either be copied to the new sheet or the concentrations sheet is renamed. In that case a new concentrations worksheet must be generated or a range error
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
243
Fig. 5 – Diazepam model implementation Following the principles explained with the minimal model the differential equations are an equivalent representation to the flowchart (see Fig. 4).
Fig. 6 – Diazepam model results The plot titled Model shows the results calculated in the model() subroutine. Calculations from the auxiliary() subroutine are presented in a separate plot called auxiliary.
244
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 8 – PB-PK model results. the appropriate alterations in the model-sheet and rerun the model.
6.
Fig. 7 – Physiological based pharmacokinetic model: structural model QF , QL , QR , QS partial flow rates and Qp air flow; Ca and Cv arterial and venous concentrations; respectively; Cin styrene concentration in the inhaled air.
will occur. The pk-engine in its present form does not handle such errors and code execution will be stopped by the Visual Basic (VB) run time environment. Similarly, syntax errors and the like in the module code or elsewhere will be detected by the VB compiler and the VB editor will appear with the offending code highlighted. This does not mean that the users need to learn VB syntax or are required to understand details about the workings of the VB editor. It is recommended to take a note of the line containing the error and then to simply close the VB editor. This will abort the current run and the user can make
Availability
The pk-engine package is available as a compressed archive via electronic mail from the authors upon request. The request should therefore include an e-mail address. The package contains the pk-engine, the support DLL and the four examples presented above. More examples have been developed by the authors including experimental alternative absorption models. We will gladly support scientists who wish to implement a particular model. Likewise we are interested to receive any bug reports describing unexpected behaviour of the software. The software is mostly harmless insofar as it contains no malign code and performs only the tasks of solving ordinary differential equations and plotting the results. Since EXCEL is also an automation server, a variety of scenarios is imaginable that run the pk-engine package from within other applications like Microsoft Powerpoint, Microsoft Word or literally any other automation client. Web based projects as front-ends could further add to user-friendliness.
Fig. 9 – Lotka-Volterra model: example of a non-linear multiplicative model. Small changes in the start values (PREY = 5, PREDATOR = 10) have pronounced effect.
c o m p u t e r m e t h o d s a n d p r o g r a m s i n b i o m e d i c i n e 8 8 ( 2 0 0 7 ) 239–245
Fig. 10 – Lotka-Volterra type oscillations. Stable attractors can only be obtained with suitable combinations of start values and model parameter values. We are looking forward to receiving the suggestions of our readers.
Conflicts of interest
245
comp names
Array of user defined compartment names, facultative
dadt
Array of 1..NN differential equations
model
Subroutine defining the set of differential equations
NA
Constant defining the number of additional results in auxiliary
NN
Constant defining the number of differential equations
plot
String of flags (0,1) determining whether a profile is plotted (1) or not (0), facultative.
results
Array of 1..NA user defined equations, output from auxiliary
ti
Constant stating the time interval between observations
ttotal
Constant stating the total time for profiling
vi
String of compartmental distribution volumes (comma separated), facultative
The authors report no conflicts of interest.
references
Appendix A Quick start: 1. Extract the contents of the pk-engine archive to a directory of your choice. 2. With the helper application setref.exe check the presence of a reference to PKENGINE in the file minimum model.xls or add this reference. 3. Start EXCEL and open the file minimum model.xls. 4. From the Dose history sheet compile and calculate the model. 5. Change the dose and recalculate. 6. Change k10 on the model-sheet, then recompile and recalculate. 7. Close minimal model.xls and repeat from step2 with diazepam model.xls or pb-pk-model.xls. Reserved variables: names and functions ai Array 1..NN describing the compartment amounts auxiliary
Optional subroutine with parameters ti, ci, results
aux names
Array of user defined names for optional data columns, facultative
axgrp
String of flags (1,2) stating the axis (primary or secondary) for the plot, facultative
ci
Array 1..NN describing compartment concentrations, input to auxiliary
[1] PK Solutions [online], Available from the World Wide Web: http://www.summitpk.com/pksolutions/pksolutions.htm. [2] GastroPlus [online], Available from the World Wide Web: http://www.simulations-plus.com. [3] SAAM II [online], Available from the World Wide Web: http://depts.washington.edu/saam2/. [4] PCCAL [online], Available from the World Wide Web: http://www.coacs.com/PCCAL/. [5] ADAPT II [online], Available from the World Wide Web: http://bmsr.usc.edu/. [6] MODELMAKER [online], Available from the World Wide Web: http://www.modelkinetix.com/. [7] G.M. Hood, PopTools version 2.6.9, Available from the World Wide Web: http://www.cse.sciro.au/poptools. [8] WinNonLin Professional Edition 2.1, Pharsight Corporation, Palo Alto, California. [9] A. Boeckmann, L. Sheiner, S. Beal, NONMEM User’s Guide, San Francisco, University of California, San Francisco, 1994. [10] R.F. Muzic, S. Cornelius, COMKAT: compartment model kinetic analysis tool, J. Nucl. Med. 42 (2001) 636–644. [11] M. Gibaldi, D. Perrier, Pharmacokinetics, Marcel Dekker, New York, 1982, p. 419. [12] S.L. Beal, Computation of the explicit solution to the Michaelis–Menten equation, J. Pharmacok. Biopharm. 11 (1983) 641–657. [13] R.G. Davies, Computer Programming In Quantitative Biology, Academic Press, London/New York, 1971, p. 388. [14] S. Haddad, M. Pelekis, K. Krishnan, A methodology for solving physiologically based pharmacokinetic models without the use of simulation software, Toxicol. Lett. 85 (1996) 113–126.