and
could be
written as
The user may define their own models for fitting or simulation by writing the appropriate subroutines CP, MODEL, and MODELOUT. To do this the user will need to have a copy of Absoft MACFORTRAN/020 (Macintosh), MS FORTRAN 5.1(MS DOS), or VAX FORTRAN (VAX VMS). The subroutines are written, compiled and then linked with the application MULTI-PART. A dummy PART subroutine should also be included or an external reference note will be displayed by the linker. The program should still run however without the subroutine PART (if integration method 4 is never used or error 75 will result).
The subroutine CP is required to define the integrated equations, differential equations, and/or initial conditions as necessary for each of the models included in the program. A number of models can be included and selected from within the subroutine CP. See the MODEL.DEMO file on this disk for the complete listing of an example subroutine CP.
SUBROUTINE CP(Y,CC,X,INIT)
C
REAL*4 CC(30),X,Y(30)
* REAL*4 T1,T2,T3,T4,t5,t6,t7,t8,TMAX
INTEGER INIT,ndose,i,j
include modcom.for
save
c
c call MFinder for MultiFinder access
c
call MFinder
The initial section of subroutine CP is shown above. Data is interfaced with the subroutine by an argument list and two COMMON areas which are specified in the file modcom.for. The statement labeled '*' is optional depending on the users model statements. The '*' should not be included. The variables are:-
Y(30) Real. On exit this array should return i) the dependent value, ii) the differential value, or iii) the initial condition for up to 30 differential equations to the main program.
CC(30) Real. On entry into CP this array contains the current values for the dependent variable, for up to 30 differential equations. This is only used with differential equation models. It is undefined for integrated models.
X Real. On entry this variable contains the current value of the independent variable (e.g. time).
INIT Integer. This variable is used as a flag with differential equation models. INIT = 0 means return initial conditions in the Y array, INIT = 1 means return the differential equation values in the Y array.
In the modcom.for file.
C(100) Real. In common, contains the values for the constants of the model(s)
P(30) Real. In common, contains the current values for the parameters
MODNUM Integer. In common, contains the number of the model selected by the user during analysis. This allows the subroutines, CP, MODEL, MODELOUT, and PART to contain information about a collection of models.
The equations:-
i) and
could be
written as
y(1) = p(1)*exp(-p(2)*x)+p(3)*exp(-p(4)*x) y(2) = c(1)*exp(-p(5)*t)
when A, B, C, D, and F are parameters and E is a constant.
ii) dCp/dt = -kel.Cp and dU/dt = V.ke.Cp, with Cp(0) = dose/V could be written as
if (init.eq.0) then
y(1) = c(1)/p(1)
y(2) = 0.0
else
y(1) = -p(2)*cc(1)
y(2) = p(1)*p(3)*cc(1)
endif
with one constant and three parameters.
1. Non-Continuous Functions
One class of pharmacokinetic model described by differential equations in subroutine CP needs special attention if the more efficient differential equation solvers (RKF45 and GEAR) are to work properly. This class of problem is the non-continuous function. This non-continuity may arise from any number of reasons, but most commonly it is the result of giving a second dose, or starting or stopping an infusion. Any time dependent change of parameter may result in a non-continuous function. Hopefully, I have taken care of this properly in Boomer. (It seems to work most of the time). If you write your own CP subroutine you will need to take a little care. Examples are shown in the source code for MODEL.DEMO.
The calculation routine will now call CP with INIT = -1 once for each time point specified in the data file if differential equations are being used. This will allow you to check the time (variable X) and compare it with the value of a particular constant to determine if a dose is to be given or an infusion is to be changed. If the time is one at which a change is made, make the change and set INIT = -2 to let the main program know of the change. This method means that a data time point equal to the dosing time change must be specified. If you don't have data at this time you can add a dummy point with concentration equal to zero. The program will automatically assign a weight of zero to this point.
Example: Multiple IV bolus. If c(2) is the time of the second dose, c(3) is the value of the second dose put into component 1 and p(2) is the apparent volume of distribution the following code could be added to subroutine CP
if (init.eq.-1) then
if (x.eq.c(2)) then
cc(1) = cc(1) + c(3)/p(2)
init = -2
return
endif
endif
Example: IV Infusion. If c(2) is the time to stop the infusion, which has started at time zero and c(3) is value of the infusion rate and p(3) is the apparent volume of distribution.
if (init.eq.-1) then
if (x.eq.c(2)) then
c(nc+1) = 0.0
init = -2
return
endif
endif
..
..
if (init.eq.0) then
cc(1) = 0.0
c(nc+1) = c(3)
return
else
y(1) = c(nc+1)/p(3) - ...
Note the use of c(nc+1) to store the value of the infusion rate. This value can be turned on or off at each appropriate time with the suitable code.
2. Mixed Models
Sometimes it is useful to be able to specify a model in terms of both differential equations and in integrated form as well. If the number of lines is set to a value greater than the number of differential equations, the calculation routines will call CP with an INIT value of -2. It is during this call that values for lines greater than the number of differential equations are calculated.
This subroutine is called during the data entry process to define a number of variables.
SUBROUTINE MODEL C * INTEGER I,maxmod,ndose * character*2 nds include modcom.for save
The initial section of subroutine MODEL is shown above. Again the lines marked '*' are optional, dependent on the model definition. Variables of interest include:-
MODNUM Integer. The main reason for the call to MODEL is to define the model number, MODNUM, and set the model specification parameters LN, M, NC, and NDEQ. The user should be presented with a selection menu and asked to choose which model. If there is only one model don't change MODNUM (i.e. don't set it equal to anything) just define the other model specification parameters. See the demo subroutines in the MODEL.DEMO text file.
LN, NL(20), M, N, NC, and NDEQ. Integer. In common, contain the number of lines, data per line, parameters, data, constants, and differential equation for the model to be analyzed.
Note that values for the variables LN, M, NC and NDEQ must be specified in the user's subroutine MODEL.
PS(30) Character (length 20), LINEST(20) Character (length 15), and CST(100) Character (length 20). On exit, these contain the name for each parameter, line, or constant of each model to aid in the data entry process.
MODST Character (length 40). On exit, this contains the name of the selected model for later printout.
IWRITE, IREAD, IDISK, and IBATCH Integer. In common, are the current unit numbers for writing, reading, saving to disk, and batch file control. Writing will be to the screen, reading maybe from the keyboard or a batch file, and idisk is usually not used. The value of ibatch depends on where the data comes from (keyboard or batch file) or whether a batch file is being produced.
This subroutine is optional, although if not used a dummy subroutine is required. That is, the following start with a return and an end statement are required. This subroutine is used to do any post-fitting processing. In some of the demo models the area under the curve is calculated in this subroutine.
SUBROUTINE MODELOUT (cc)
C
* REAL*4 AUC,czero,beta
* real*4 k10,k12,k21,v1,vbeta,vss,tbc
* INTEGER I
include modcom.for
save
The only additional variables, not previously described, are:-
CY(200) and TX(200) Real. These are the observed data fitted by the program. This information is necessary for AUC calculations.
CC(200) Real. This is calculated data, used at times to calculate AUC etc.
If Gear's method is selected for differential equations with the optional subroutine PEDERV, then a subroutine PART must be written. If this option is not used then it is still a good idea to include a dummy PART subroutine in the form of the following statements followed by a return and an end statement.
SUBROUTINE part(t,cc,pw)
real*4 t,cc(30),pw(30*30)
* real*4 t1
* integer i
include modcom.for
save
The variables used or defined in this subroutine are:-
T Real. One entry, this is the current value of the independent variable, time. This value should not be changed in subroutine PART.
CC(30) Real. On entry this is the current value of the dependent variable for up to 30 differential equations.
PW(30*30) Real. On exit, this array of ndeq*ndeq elements should contain the partial derivative of each equation with respect to each dependent variable. The order is equation 1 with respect to variable 1, eq 2 with var 1, eq 3 with var 1,...,eq ndeq with var 1, eq 1 with var 2, eq 2 with var 2,...,eq 1 with var ndeq,...eq ndeq with var ndeq. For demo model 13 the model equations are:-
diff eq #1 y(1) = -c(1)*cc(1)+c(2)*cc(2)*cc(3) diff eq #2 y(2) = c(1)*cc(1)-c(2)*cc(2)*cc(3)-c(3)*cc(2)*cc(2) diff eq #3 y(3) = c(3)*cc(2)*cc(2)
and the equations for PW are:-
c partials with respect to variable cc(1)
c
diff eq #1 pw(1) = -c(1)
diff eq #2 pw(2) = c(1)
diff eq #3 pw(3) = 0.0
c
c partials with respect to variable cc(2)
c
diff eq #1 pw(4) = c(2)*cc(3)
diff eq #2 pw(5) = -c(2)*cc(3)-2.0*c(3)*cc(2)
diff eq #3 pw(6) = 2.0*c(3)*cc(2)
c
c partials with respect to variable cc(3)
c
diff eq #1 pw(7) = c(2)*cc(2)
diff eq #2 pw(8) = -c(2)*cc(2)
diff eq #3 pw(9) = 0.0
Material on this website should be used for Educational or Self-Study Purposes Only
Copyright © 2001-2008 David W. A. Bourne (david@boomer.org)