Back to the Top
Hello,
I am trying to model the PK of a drug and its metabolite
simultaneously. The parent drug is modeled with 2-compartments, with
elimination from the central compartment. The metabolite is modeled
as a 3rd compartment with its own elimination. We have concentration
data for the parent drug and the metabolite in the central compartment.
I am attempting to write a user compiled model (using differential
equations) using Fortran for WinNonMix, but for some reason the
software is unable to solve the differential equations, and I keep
getting errors. I don't think there is a problem with the equations
but I am not sure. I have placed the fortran code below. I suppose
I can try to get numerical solutions for the differential equations
using Laplace transforms, but that is very cumbersome.
So here are my questions:
1. Do you see any problems with the fortran code or way the
differential equations are written that would cause WinNonMix to be
unable to solve them?
[What errors are you getting? Is it a stiff system - db]
2. Does anyone have a numerical solution to these differential
equations?
* This way I could write the fortran code using algebraic functions,
rather than differential equations.
[See http://www.boomer.org/c/p4/c02/c07/c0707.html - db]
Thanks in advance for your help. The code is below:
--
subroutine usrmod(mode, funcno, dta, f, p, dz, s, con, x, z, wt, y,
objfn)
implicit real*8(a-h,o-z)
implicit integer(i-n)
!ms$attributes dllexport :: usrmod
integer mode, funcno
real*8 p(1), con(1), x, f, s(1), z(1), dz(1), dta(1), wt(1), y, objfn
!
************************************************************************
******
!
*
*
!* Block 1: declare any model specific variables (if any),
e.g. *
!
*
*
!* real*8
var *
!
*
*
!
************************************************************************
******
real*8 k10, k12, k21, v, Dose, TSTART, TSTOP, Q, cl1, v2, cl2
real*8 k13, k30, v3, cl3, k31
integer i, k
!---- end of block (1)
--
!
************************************************************************
******
!
*
*
!* Block 2: initialize variables (if any),
e.g. *
!
*
*
!* var =
1.0d0 *
!
*
*
!
************************************************************************
******
Dose = con(1)
TSTART = con(2)
TSTOP = con(3)
Q = Dose/(TSTOP - TSTART)
v = p(1)
k10 = p(2)
k12 = p(3)
k13 = p(4)
k21 = p(5)
k31 = p(6)
k30 = p(7)
!---- end of block (2)
--
if (mode == 2) then
!
************************************************************************
******
!
*
*
!* Block 3: set function value (required). It must include the
following *
!*
statement: *
!
*
*
!* f = ;value of the model
function *
!
*
*
!
************************************************************************
******
if (funcno == 1) f = z(1)
if (funcno == 2) f = z(3)
!---- end of block (3)
--
else if (mode == 5) then
!
************************************************************************
******
!
*
*
!* Block 4: set secondary parameter values (if any). Depending
on *
!* the number of secondary parameters, one or more of
the *
!* following statements may appear in this
block *
!
*
*
!* s(1) = ;value of the first secondary
parameter *
!
*
*
!
************************************************************************
******
!---- end of block (4)
--
else if (mode == 4) then
!
************************************************************************
******
!
*
*
!* Block 5: Set initial values of differential equations (if any).
Depending *
!* on the number of differential equations involved in the
model, *
!* one or more of the following statements may appear in
this *
!*
block: *
!
*
*
!* z(1) = ;initial value of the 1st differential
equation *
!
*
*
!
************************************************************************
******
z(1) = 0.0d0
z(2) = 0.0d0
z(3) = 0.0d0
!---- end of block (5)
--
else if (mode == 3) then
!
************************************************************************
******
!
*
*
!* Block 6: Describe differential equations (if any). Depending
on *
!* the number of differential equations involved in the
model, *
!* one or more of the following statements may appear in
this *
!*
block: *
!
*
*
!* dz(1) = ;value of the 1st differential
equation *
!
*
*
!
************************************************************************
******
if (x <= TSTOP) then
dz(1) = Q/v - (k13 + k10 + k12) * z(1) + k21 * z(2) + k31 * z
(3)
dz(2) = k12 * z(1) - k21 * z(2)
dz(3) = k13 * z(1) - k31 * z(3) - k30 * z(3)
else
dz(1) = - (k13 + k10 + k12) * z(1) + k21 * z(2) + k31 * z(3)
dz(2) = k12 * z(1) - k21 * z(2)
dz(3) = k13 * z(1) - k31 * z(3) - k30 * z(3)
end if
!---- end of block (6)
--
end if
!
************************************************************************
******
!
*
*
!* Block 7: description of local subroutines/functions
(optional) *
!
*
*
!* Note: if any local subroutine/function is to be included in this
block, *
!* the following statement must appear before any of the
subroutine/ *
!*
function: *
!
*
*
!*
contains
*
!
*
*
!
************************************************************************
******
!---- end of block (7)
--
end subroutine usrmod
--
Tarek A. Leil, Ph.D.
Department of Oncology, Guggenheim 1311
Mayo Clinic
200 First Street SW
Rochester, MN 55905
ph: 507-538-4227
fax: 507-266-5146
Back to the Top
The following message was posted to: PharmPK
It is hard to suggest to you why you are having errors when you don't
tell us what they are!
It is unlikely that the metabolite returns reversibly to the parent.
If this reaction is not reversible then you do not need k31 * z(3) in
the differential equations for parent (compartment 1) and metabolite
(compartment 3)
You have a start and stop time for the infusion which you use to
calculate the input rate Q but you do not use the start time in the
differential equations. If TSTART is not 0 then the input duration
will be TSTOP not (TSTOP-TSTART).
You also need to be careful about which parameters you are trying to
estimate. Because you do not give metabolite directly but only
observe the concentration derived from the parent you cannot identify
k13 separately from k10. Because you cannot estimate k10 directly you
can simply make the assumption that all the parent is converted to
the metabolite ie. k10=0.
A simpler way of coding your system which includes these suggestions
would be:
Dose = con(1)
TSTART = con(2)
TSTOP = con(3)
v = p(1)
if (x>=TSTART and x<=TSTOP) then
Q = Dose/(TSTOP - TSTART)/v
else
Q = 0
end if
dz(1) = Q - (k13 + k12) * z(1) + k21 * z(2)
dz(2) = k12 * z(1) - k21 * z(2)
dz(3) = k13 * z(1) - k30 * z(3)
--
Nick Holford, Dept Pharmacology & Clinical Pharmacology
University of Auckland, 85 Park Rd, Private Bag 92019, Auckland, New
Zealand
email:n.holford.-a-.auckland.ac.nz tel:+64(9)373-7599x86730 fax:373-7556
http://www.health.auckland.ac.nz/pharmacology/staff/nholford/
Back to the Top
The following message was posted to: PharmPK
Dear Tarek A. Leil, Ph.D.
As an alternative to a compartment method I am suggesting a
system-approach based method for modeling formation of a metabolite from
a parent drug. The use of the latter method is exemplified in: Dedik
L., \0x00uri\0x00ova M.: System approach to modeling metabolite
formation from
parent drug: A working example with methotrexate. Methods Find Exp Clin
Pharmacol 24, 2002, 481-486.
With best regards,
Maria Durisova, DrSc (Math/Phys)
http://www.uef.sav.sk/
http://www.uef.sav.sk/advanced.htm
http://www.klinickafarmakologie.cz/artkey/far-20030
PharmPK Discussion List Archive Index page
Copyright 1995-2010 David W. A. Bourne (david@boomer.org)