“All models are wrong, some models are useful” George Box 1979 (1)

A model is a representation of observed data or an observed system. Models may be useful to describe a collection of observations and their relationship with one another. More importantly, in attempting to explain and predict observations and their relationships, modelling can be the mechanism to connect description and understanding. Pharmacometrics involves the analysis and interpretation of data and seeks to bridge stamp collecting and physics (to paraphrase Ernest Rutherford quote (2)). 

Metrics is a term used to denote a system of measurement. It involves a certain way of assessing and interpreting parameters, and is used in a variety of fields including economics, education, politics as well as pharmacology. However in pharmacometrics, as with any metrics system, quantitative measurement and estimation does not automatically imply validity and truth.

Models vary in complexity; but inherently modelling involves simplification. Not only does the model have to deal with the known knowns, the known unknowns, but perhaps also the unknown knowns, and the unknown unknowns (paraphrasing a confused/confusing quote from Donald Rumsfield (3)). So while no model can ever speak the truth, a model may be more or less valid depending on the perspective of the modeller and the problem to be solved. There is error inherent in the process of observation and collection of information. Error and variability exists in every step of model building and testing. More useful models attempt to quantify and classify error and variability.

“In pharmacology models are ultimately expressed by mathematical equations” (Bourne 1995 (4)). Equations can encapsulate and describe large amounts of information. Differential equations are used to describe rate of change of a variables over time. There are different numerical methods of solving differential equations such as the Runge-Kutta method (RK4) developed by  mathematicians several decades ago. The robustness of this and other forms of iterative solution has been demonstrated by their successful use for satellite exploration of the solar system. Improvements in accuracy with which a model describes observations and estimatesparameters, may come from more powerful and robust mathematical functions and repeated observations.

Pharmacometrics uses mathematical modelling enabled by software applications to estimate parameters and define and predict relationships between variables. Several different computer software programs can be used for this purpose, all with unique features and limitations. These software programs can render complex mathematics useful to the nonmathematician in solving pharmacometric problems.

In order to use pharmacometric modelling, it is necessary to understand the components of model equations. Equations summarise large amounts of data with a small number of parameter values. A parameter is a constant that is estimated, e.g. volume of distribution, clearance, rate constants for elimination/absorption/equilibration, EC50. Parameters are adjusted during the modelling process to obtain the best fit to the data. Variables that are fixed and do not change in the modelling process are true constants e.g. dose, infusion time. The dependant variable is the variable whose value the equation describes, e.g. concentration in pharmacokinetic equations, effect in pharmacodynamic equations. Parameters and the dependant variable are related by mathematical operators. Time is the only truly independent variable, as it is the only variable we cannot control.

In pharmacology models can not only describe data such as dose-concentration, concentration-effect, but more usefully predict these relationships. This is most useful if the model can be used to define parameters not just for the data observed for particular individuals at particular times, but for different individuals at different times. Population modelling has been used since the 1970s to describe and predict pharmacokinetic, pharmacodynamic and linked PKPD relationships for a population rather than just individuals. NONMEM is the most widely used software tool for population modelling, and uses Fortran source code.

Simulation involves the generation of data from a proposed model equation. Simulation has a role in both model generation and model evaluation. Software programs such as Excel and Berkeley Madonna involve simulation of data. Software programs such as NONMEM, Monolix and WinNonlin use nonlinear regression analysis to simulate data and estimate parameters. Observed values of collected data can be compared with predicted values for the dependant variables.

Differences in predicted and observed values may be due to error and variation. The goal of model development is to find parameter values that decrease the differences between predicted and observed values, in other words to develop a model that fits the data and describes it well. Models are evaluated by tabular and graphical representations and statistical methods. Error modelling is used to describe the variability in how well the data is described by a parameter. In fact measurement of error involves estimation. Statistical error describes the difference between observed and expected values. Residual error is an observable estimate of unobservable statistical error. Residuals is the term given to describe the difference between group mean and individual values in a sample. Standard error (SE) is a measure of precision. It is used to estimate a population parameter from a sample. Standard error is used to define a range in which the true population mean value should lie. Confidence intervals are calculated from standard error. Coefficient of variation is another measure of variability, and is equal to the standard deviation divided by the mean (multiplied by 100 to express as a percentage). All modelling solutions involve assumption. A common assumption is that the independent variable can be measured without significant error, however this is unlikely to be always true.

Modelling is used to quantify and summarise data and furthermore to explain and explore mechanism and predict values. The usefulness of a model depends on the problem to be solved, the population to which it is applied, and the perspective of modeller.

There are different types of structural models used in pharmacology: compartment models, mechanistic or physiological models, empirical models, and hybrid approaches. Traditional compartmental models describe a number of well stirred compartments. The number of compartments and structure of the model is determined by the data and route of administration. Noncompartmental modelling can be described as including nonparametric data, and a more general approach rather than a structural system. Mechanistic models are based on physical and physiological principles. For example physiological models may include factors describing blood flow, elimination rate, partition coefficients, diffusion, and the kinetics of receptor binding. Empiric modelling however requires few assumptions about the data generating mechanisms and is useful when little is known about underlying physiological processes (black boxes).

Model equations and parameter estimates are useful to summarise large amounts of data. Excel & Berkeley Madonna simulate data using initial estimates and a model equation. This is a potentially useful step for developing a model. MONOLIX and NONMEM are more powerful software packages that perform nonlinear regression analysis and can be used to estimate parameters, in addition to simulation.

Introduction by Anita Sumpter (2008).

References

1. Box GEP. Robustness in the strategy of scientific model building. In: Launer RL, Wilkinson GN, editors. Robustness in Statistics; 1979. p. 202.
2. Reference to Ernest Rutherford quote “Science is either stamp collecting or physics”. http://en.wikiquote.org/wiki/Ernest_Rutherford. All science is either physics or stamp collecting. As quoted in "Rutherford at Manchester" (1962) by J. B. Birks.
3. Paraphrased Donald Rumsfield quote “ http://en.wikiquote.org/wiki/Donald_Rumsfield. “Reports that say that something hasn't happened are always interesting to me, because as we know, there are known knowns; there are things we know we know. We also know there are known unknowns; that is to say we know there are some things we do not know. But there are also unknown unknowns -- the ones we don't know we don't know. And if one looks throughout the history of our country and other free countries, it is the latter category that tend to be the difficult ones. “ Department of Defense news briefing, February 12, 2002 [1].
4. Bourne D. Mathematical Modeling of Pharmacokinetic Data. CRC Press, Boca Raton, 1995. Preface x.

Note: All files should be loaded from and saved to your Pharmacometrics Data\User Defined Models folder for this assignment. You may find it useful to create a separate folder for each program inside the User Defined Models folder as you will create several files during this assignment.

Some hints (shown in grey boxes) refer to changes needed for subsequent problems.

You may find it easier to read  this page if your click on the "Printer Friendly version" on the right hand side and then click "Cancel" when prompted for a printer

1. Open Berkeley Madonna shortcut in folder "Pharmacometrics Programs".
2. Add code to the Equations window defining the STARTTIME and STOPTIME, the output interval (DTOUT), the model parameters and the model equation (code shown in Figure 1).
3. Click on Run to run the model.
4. Save your model with the name bo1.mmd (in your Pharmacometrics Data\User Defined Models\Berkely Madonna folder).
5. View a table of times and concentrations by clicking on the Table icon in the Run 1 graph window.

Repeat the above steps for the k01 and emaxc models.

 Figure 1. Code for bo1.mmd (Berkeley Madonna)

1. Open the Monolix shortcut in folder "Pharmacometrics Programs".
2. Click on the New Project icon (top left of Monolix GUI).
3. To export the simulated bo1 data from 3models.xls to a format that Monolix can read, create a new Excel workbook and add the following headings in the top row. It is important to put the # before ID so that the same file can be used for NONMEM.
   
   #ID TIME AMT DV
   
   
   For k01, you will also need a RATE column.
   For emaxc you will only need #ID CONC and DV columns.
4. The next row will be a dose record, which is different to an observation record. Enter 1 into the ID column, 0 into the TIME column, the dose into the AMT column and put a period (.) into the DV column. A period is used in programs such as Monolix and  NONMEM when a value is not required.
   
   
   For k01, the AMT for the dose record is the total dose (Rate * infusion duration) and the RATE is the Rate. For emaxc, the CONC and DV for the dose record are both 0.
5. Rows 3 to 13 are observation records. Fill the ID column with the value 1. Copy the times from 3models.xls into the TIME column. Fill the AMT column with periods and copy the concentrations from 3models.xls into the DV column. Your worksheet should look similar to Table 1, but your DV values will be different.
   
   

   	A	B	C	D
   1	#ID	TIME	AMT	DV
   2	1	0	100	.
   3	1	0	.	10.13898
   4	1	1	.	7.346313
   5	1	2	.	5.270305
   6	1	3	.	4.082883
   7	1	4	.	3.342263
   8	1	5	.	2.244587
   9	1	6	.	1.689803
   10	1	7	.	1.534163
   11	1	8	.	0.951986
   12	1	9	.	0.612068
   13	1	10	.	0.601922

   Table 1. Monolix input file for bo1.
6. Now save the data file using 'Save As…' and choose CSV (Comma delimited, *.csv) format. Name the file bo1.csv.
7. In Monolix , click 'The data', locate your bo1.csv file. This can be done by copying the path to the folder from Windows Explorer into the Monolix "Select datafile File name:" box.   e.g. P:\My Pharmacometrics\Pharmacometrics Data\User Defined Models
8. Select bo1.csv and click open. The data information window will appear.
9. Monolix should identify that your input file has a header row and your data should appear under the headings: ID DOSE TIME Y. If the 'Use header' box is not checked then click on it so that Monolix uses the headers in the input file. Click accept if Monolix has assigned the correct headings.
   
   
   You can select the type of data associated with the headings manually by clicking the drop down menus. If the independent variable is not time then select X for the type.
10. Click 'The structural model' and the model library will appear. Here we can select from prebuilt models or write our own using MLXTRAN. First you will use a prebuilt bolus dose model. Select the bolus_1cpt_VCl model. A description of the model appears under Model Infos. The parameters to be estimated are listed under Variable names. Click accept.
    
    
    The Monolix library models for the k01 and emaxc problems have these names:
    k01: infusion_1cpt_VCl
    emaxc: immed_emax_null (you will need to change to the PD library by clicking 'Other list').
11. Change the initial parameter estimates under Fixed effects. Holding the mouse over each box will tell you, which parameter the value is for.
12. Set the 'Stand. dev. of the random effects' to 0 for each parameter and fix it by right clicking on each box and selecting fix (it will turn pink and will not be estimated).
    
    
    Because these are data from a single individual, there is no between subject variability (random effects). The SD of random effects is therefore 0.
13. Set the residual error model to 'exp'.
14. Set the 'Residual error parameters' to 0.1 (CV=10%).
15. Click 'Check initial fixed effects'. A plot of predictions based on the model and initial parameter estimates will display along with the observed values.
    You can visualise the effect of changing your parameter estimates by adjusting the values in the bottom left of the window. When you have chosen initial estimates that form a prediction that is similar to the observations, click 'Set as initial values' to apply these values and close the window.
16. Save the project as bo1.mat in your User Defined Models\Monolix folder.
17. Run the model by clicking 'Estimate the Population Parameters' at the top of the Monolix window. This will take a while depending on the complexity of the model. During the estimation process you can see how the parameter estimates are being changed and settle down towards the final value.
18. When the estimation finishes click on 'Settings Results'. Choose 'Individual Fits' then at 'VPC'. Look at the graphs to see how well the model and the final parameters fit the data.
19. View the parameter estimates by clicking 'Last Results'. The parameter estimates are also displayed in the command prompt that is behind the main Monolix window.

Repeat the above steps for the k01 and emaxc models.

1. Open the NONMEM shortcut in folder "Pharmacometrics Programs".
2. Change directory to the User Defined Models folder by typing this command in the NONMEM window then press <ENTER>.
   
   
   cd User*
   
   
   
   All commands shown in the green boxes should end by pressing <Enter>.
   
   
3. Make a NONMEM directory then change to the NONMEM folder by typing this command in the NONMEM window.
   
   
   md NONMEM
   
   
   
   You only have to make the directory once. You dont have to give the md NONMEM command when you are working on the K01 and emaxc problems.
   
   
4. Change to the NONMEM directory if you are in the User Defined Models directory..
   
   
   cd NONMEM
   
   
   
5. Use the Windows Notepad to create the following code in a file named bo1.ctl e.g. by using the command:
   
   style="border-right: dimgray 1px dashed; padding-right: 5px; border-top: dimgray 1px dashed; padding-left: 5px; padding-bottom: 5px; border-left: dimgray 1px dashed; width: 213px; padding-top: 5px; border-bottom: dimgray 1px dashed; background-color: #00cc00"> notepad bo1.ctl
   
6. Enter the code for bo1.ctl shown in Figure 2.
   
   
   $PROB one compartment bolus
   $INPUT ID TIME AMT DV
   $DATA bo1.csv
   IGNORE #
   IGNORE (AMT.GT.0)
   $ESTIM METHOD=CONDITIONAL
   $COV
   $THETA (0,3,) ; POP_CL
   $THETA (0,10,) ; POP_V
   $OMEGA 0 FIX ; ETA_CL
   $OMEGA 0 FIX ; ETA_V
   $SIGMA 0.01 ; ERR_CV
   $PRED
   DOSE=100
   CL=THETA(1)*EXP(ETA(1))
   V=THETA(2)*EXP(ETA(2))
   F= DOSE/V*EXP(-CL/V*TIME)
   Y=F + F*ERR(1)
   $TABLE ID TIME Y
   NOPRINT ONEHEADER FILE=bo1.fit
   

   Figure 2. NM-TRAN code for bo1.ctl (NONMEM)

   The $OMEGA parameters are FIXed to 0. This is because there is only one individual being modelled. These random effects parameters are used to describe between subject variability when there is more than one individual.
   
   
7. Save the file. Check using Windows Explorer that you can find the file “bo1.ctl” in the User Defined Models\NONMEM folder.
8. Copy your bo1.csv file from your User Defined Models\Monolix folder to your NONMEM folder. The data file name must match in the $DATA record of the bo1.ctl file.
9. Open each .csv file and type # before the the ID column values of rows 1,2 and 3 (header and dose record and first observation at TIME=0).
   
   
   NONMEM ignores rows that start with # (defined by IGNORE # in the control stream [.ctl file]). You could just delete these rows, but if you open the file again at a later date, you may not remember what your data represent.
   
   NONMEM can also ignore rows based on the values in the column. In this example the AMT row is dropped. This is because we are defining the dose directly in this $PRED model. More commonly the PREDPP library is used and this requires the dose to be defined in the AMT column.
   
   
10. Execute NONMEM with this command in the NONMEM window:
    
    
    nmgo bo1
    
    
    
    When you get errors from NONMEM with the nmgo command then please read the error message carefully and try to understand what it is telling you. The usual errors that occur with these example problems will give you some clues to what you might need to change in your ctl file.
    
    
    Commands can be recalled by using the Up arrow on your keyboard. You can easily repeat the command without more typing or edit it to save the amount of typing you do.
    
    
11. Use the Windows notepad or EditPlus to open this file and look at the saved results e.g. with the following command:
    
    
    notepad bo1.std\bo1.smr
    
    
12. Repeat the above steps for the K01 and emaxc models. For Monolix you will need to find a suitable model in the Monolix library of models. For NONMEM you should change the bo1.ctl file and save it with a suitable name such as K01.ctl or emaxc.ctl..