Supplemental material
Example4: Lifespan indirect response model with a precursor pool
Theoretical
SemiDDE/Category 2.5
The LIDR model (43)-(45) is a Semi-DDE/Category 2.5 model and reads
/ / (S01)/ / (S02)
/ / (S03)
Category 2.5 ODE formulation
The equivalent ODE Category 2.5 form reads for
/ / (S04)/ / (S05)
/ / (S06)
/ / (S07)
/ / (S08)
and for the right hand side is
/ (S09)/ (S10)
/ (S11)
/ (S12)
/ (S13)
Category 2.5 ODE formulation with lag-time
The system to implement in PKPD software without a DDE solver reads:
/ / (S14)/ / (S15)
/ / (S16)
/ / (S17)
/ / (S18)
wherethe Dirac delta function is modeled by the lag-time feature from the PKPD software.
Implementation
All parameter estimates are similar from the PKPD programs, except for small differences based on tolerance settings and error models.
acslX
The acslX code for the DDE formulation (S01)-(S03) reads:
PROGRAM
INITIAL
constant dose = 1000
schedulestartdose .at. 0.0
scheduleenddose .at. 0.1
END
DYNAMIC
ALGORITHM IALG = 15
NSTEPS NSTP = 100
MAXTERVAL MAXT = 1
MINTERVAL MINT = 1.0e-9
CINTERVAL CINT = 0.01
DERIVATIVE
constantkel = 0.25, V = 1
constantSmax = 50, SC50 = 1
constant k0 = 0.5, k1 = 0.05, Tau = 20
x10 = 0
x20 = k0/k1
x30 = k0*Tau
c = x1/V
x2_del = delay(x2,x20,Tau,100000)
x1 = integ(ivbolus-kel*x1,x10)
x2 = integ(k0*(1+(Smax*c)/(SC50+c))-k1*x2,x20)
x3 = integ(k1*x2-k1*x2_del,x30)
response = log(x3+0.0001)
END
DISCRETE startdose
ivbolus = dose/0.1
END
DISCRETE enddose
ivbolus = 0.0
END
CONSTANT TSTOP = 100.0
TERMT (T .GE. TSTOP, 'REACHED TSTOP')
END
END
ADAPT-5
The subroutines of ADAPT-5 code for the Category 2.5 ODE formulation (S14)-(S18) reads:
Subroutine SYMBOL
Implicit None
Include 'globals.inc'
Include 'model.inc'
NDEqs = 5
NSParam = 7
NVparam = 1
NSecPar = 0
NSecOut = 0
Ieqsol = 1
Descr = ''
Psym(1) = 'Smax'
Psym(2) = 'SC50'
Psym(3) = 'k0'
Psym(4) = 'k1'
Psym(5) = 'Tau'
Psym(6) = 'kel'
Psym(7) = 'V'
PVsym(1) = 'sig1'
Return
End
Subroutine DIFFEQ(T,X,XP)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Real*8 T,X(MaxNDE),XP(MaxNDE)
Real*8 c,cdel
Real*8 Smax,SC50,k0,k1,Tau,kel,V
Smax = P(1)
SC50 = P(2)
k0 = P(3)
k1 = P(4)
Tau = P(5)
kel = P(6)
V = P(7)
c = x(1)
cdel = x(3)
xp(1) = -kel*x(1)
xp(2) = k0*(1+(Smax*c)/(SC50+c)) - k1*(x(2)+k0/k1)
if (t .lt. Tau) then
xp(3) = 0
xp(4) = 0
xp(5) = k1*(x(2)+k0/k1) - k0
else
xp(3) = -kel*x(3)
xp(4) = k0*(1+(Smax*cdel)/(SC50+cdel)) - k1*(x(4)+k0/k1)
xp(5) = k1*(x(2)+k0/k1) - k1*(x(4)+k0/k1)
end if
Return
End
Subroutine OUTPUT(Y,T,X)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Real*8 Y(MaxNOE),T,X(MaxNDE)
Real*8 k0,Tau
k0 = P(3)
Tau = P(5)
y(1) = log(x(5) + k0*Tau + 0.0001)
Return
End
Subroutine VARMOD(V,T,X,Y)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Real*8 V(MaxNOE),T,X(MaxNDE),Y(MaxNOE)
v(1) = (pv(1)*y(1))**2
Return
End
Berkeley Madonna
The Berkley Madonna code for the DDE formulation (S01)-(S03) reads:
METHOD STIFF
STARTTIME = 0
STOPTIME = 100
DTMIN = 1e-6
DTMAX = 1e-2
TOLERANCE = 1e-6
DT = 1e-6
Dose = 1000
kel = 0.25
V = 1
Smax = 50
SC50 = 1
k0 = 0.5
k1 = 0.05
Tau = 20
INIT x1 = dose/V
INIT x2 = k0/k1
INIT x3 = k0*Tau
d/dt(x1) = - kel*x1
d/dt(x2) = k0*(1+(Smax*x1)/(SC50+x1)) - k1*x2
d/dt(x3) = k1*x2 - k1*delay(x2,tau)
response = logn(x3+0.0001)
DISPLAY response
DISPLAY Smax,SC50,k0,k1,Tau
MATLAB
The Matlab code for the DDE formulation (S01)-(S03) reads:
clear all;
clear global all;
global g_params;
opt_dde = ddeset('RelTol',1e-8,'AbsTol',1e-8);
kel = 0.25;
V = 1;
Smax = 86.74;
SC50 = 9.888;
k0 = 0.3135;
k1 = 0.052198;
Tau = 33.71;
g_params = [kel,V,Smax,SC50,k0,k1,Tau];
sol = dde23(@fun_PreLSM_Exa4,[Tau],@fun_PreLSM_Exa4_hist,[0, 100],opt_dde);
t_plot = sol.x;
y_plot = sol.y;
figure(1); plot(t_plot,log(y_plot(3,:)+0.0001),'black'); hold on;
function dx = fun_PreLSM_Exa4(t,x,Z)
globalg_params;
kel = g_params(1);
V = g_params(2);
Smax = g_params(3);
SC50 = g_params(4);
k0 = g_params(5);
k1 = g_params(6);
dx(1,1) = -kel*x(1);
dx(2,1) = k0*(1+(Smax*x(1))/(SC50+x(1))) - k1*x(2);
dx(3,1) = k1*x(2) - k1*Z(2,1);
function erg = fun_PreLSM_Exa4_hist(t)
globalg_params;
dose = 1000;
V = g_params(2);
k0 = g_params(5);
k1 = g_params(6);
T = g_params(7);
erg = [dose/V, k0/k1, k0*T];
MONOLIX
The MONOLIX code for the DDE formulation (S01)-(S03) reads:
DESCRIPTION:
PreLSM_II
INPUT:
parameter= {k,V,Smax,SC50,k0,k1,Tau}
PK:
compartment(cmt=1, volume=V, concentration=c)
elimination(cmt=1, k)
iv(cmt=1)
EQUATION:
t0 = 0
y1_0 = k0/k1
y2_0 = k0*Tau
ddt_y1 = k0*(1+(Smax*c)/(SC50+c)) - k1*y1
ddt_y2 = k1*y1 - k1*delay(y1,Tau)
out = log(y2+0.00001)
OUTPUT:
output = {out}
NONMEM
The NONMEM code for the Category 2.5 ODE formulation (S14)-(S18) reads:
$PROB
$INPUT ID TIME AMT CMT DV
$DATA data_log_v2.txt IGNORE #
$SUB
ADVAN13 TOL=8
$MODEL
COMP = (Comp1)
COMP = (Comp2)
COMP = (Comp3)
COMP = (Comp4)
COMP = (Comp5)
$PK
kel = Theta(1)
V = Theta(2)
Smax = Theta(3)
SC50 = Theta(4)
k0 = Theta(5)
k1 = Theta(6)
Tau = Theta(7)
A_0(1) = 0
A_0(2) = k0/k1
A_0(3) = 0
A_0(4) = k0/k1
A_0(5) = k0*Tau
ALAG3 = Tau
$DES
c = A(1)/V
cdel = A(3)/V
DADT(1) = -kel*A(1)
DADT(2) = k0*(1+(Smax*c)/(SC50+c)) - k1*A(2)
DADT(3) = -kel*A(3)
DADT(4) = k0*(1+(Smax*cdel)/(SC50+cdel)) - k1*A(4)
DADT(5) = k1*A(2) - k1*A(4)
$ERROR
PRED = A(5)
Y = log(PRED*(1+err(1))+0.0001)
$THETA
0.25 FIX
1 FIX
50 ; Smax
1 ; SC50
0.5 ; k0
0.05 ; k1
20.0 ; Tau
$OMEGA
1e-2
$EST METHOD = 0 NSIG=8 MAXEVAL=9999 NOABORT PRINT = 1
$TABLE ID TIME FILE=xptab1.txt
NONMEM data set (fragment):
#IDTIMEAMTCMTDV
1010001.
1010003.
10.52.367657649
14.52.971483993
18.54.030494105
Phoenix PLM
The PLM code for the ODE formulation (S14)-(S18) reads:
PHX_C25ODE(){
deriv(x = -kel*x)
deriv(P = k0*(1+(Smax*x/V)/(SC50+x/V))-k1*P)
deriv(xdel = -kel*xdel)
deriv(Pdel = k0*(1+(Smax*xdel/V)/(SC50+xdel/V))-k1*Pdel)
deriv(R = k1*P - k1*Pdel)
dosepoint(x)
dosepoint(xdel,tlag=Tau)
sequence{P = k0/k1}
sequence{Pdel = k0/k1}
sequence{R = k0*Tau}
error(REps = 1)
observe(RObs = log(R * (1 + REps) + 0.0001))
stparm(kel = 0.25, V = 1)
stparm(Smax = tvSmax, SC50 = tvSC50, k0 = tvk0, k1 = tvk1, Tau = tvTau)
fixef(tvSmax = c(0, 50, 100))
fixef(tvSC50 = c(0, 1, 20))
fixef(tvk0 = c(0, 0.5, 1))
fixef(tvk1 = c(0, 0.05, 0.1))
fixef(tvTau = c(0, 20, 40))
}
S-ADAPT
The S-ADAPT code for the DDE formulation (S01)-(S03) reads:
Subroutine DEL_DIFFEQ(T,X,XP,RD)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Include 'tdelay.inc'
Real*8 T,X(*),XP(*),RD(*)
Real*8 kel,V,k0,k1,Smax,SC50
Real*8 c
TD(1) = P(1)
kel = P(2)
V = P(3)
k0 = P(4)
k1 = P(5)
Smax = P(6)
SC50 = P(7)
c = X(1)/V
XP(1) = -kel*X(1)
XP(2) = k0*(1 + Smax*c/(SC50+c)) - k1*X(2)
XP(3) = k1*X(2) - k1*XD(2,1)
Return
End
Subroutine DEL_OUTPUT(Y,T,X)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Include 'tdelay.inc'
Real*8 Y(*),T,X(*)
Real*8 kel,V,k0,k1,Smax,SC50
Real*8 c
NOEQS = 1
TD(1) = P(1)
kel = P(2)
V = P(3)
k0 = P(4)
k1 = P(5)
Smax = P(6)
SC50 = P(7)
X0(1) = 0
X0(2) = k0/k1
X0(3) = k0*TD(1)
Y(1)=log(X(3)+0.0001)
Return
End
Subroutine SYMBOL
Implicit None
Include 'globals.inc'
Include 'model.inc'
Include 'tdelay.inc'
character*60 descr
common /descr/ Descr
NDE = 3
NDEL = 1
NDEqs = NDE*(NTD+1)
NSParam = 7
NVparam = 1
Ieqsol = 1
NTPARAM = 5
Descr = 'Example 4 - PreLSM II'
PSYM(1) = 'T'
PSYM(2) = 'kel'
PSYM(3) = 'V'
PSYM(4) = 'k0'
PSYM(5) = 'k1'
PSYM(6) = 'Smax'
PSYM(7) = 'SC50'
PVSYM(1) = 'RESV'
Return
End
Subroutine DEL_VARMOD(V,T,X,Y,J)
Implicit None
Include 'globals.inc'
Include 'model.inc'
Include 'tdelay.inc'
Real*8 V(*),T,X(*),Y(*)
V(1) = PV(1)**2
Return
End
WINONLIN CLASSIC
The WINONLIN CLASSIC code for the ODE formulation (S14)-(S18) reads:
MODEL LIDR_Pre
COMMANDS
NFUNCTIONS 1
NPARAMETERS 5
PNAMES 'Smax','SC50','k0','k1','Tau'
NDERIVATIVES 3
END
TEMPORARY
T = X
kel = 0.25
V = 1
dose = 1000
END
START
Z(1) = k0/k1
Z(2) = k0/k1
Z(3) = k0*Tau
END
DIFFERENTIAL
IF T >= 0 THEN
c = (dose/V)*exp(-kel*T)
ELSE
c = 0
ENDIF
IF T >= Tau THEN
cdel = (dose/V)*exp(-kel*(T-Tau))
ELSE
cdel = 0
ENDIF
kin = k0*(1+(Smax*c)/(SC50+c))
kindel = k0*(1+(Smax*cdel)/(SC50+cdel))
DZ(1) = kin - k1*Z(1)
DZ(2) = kindel - k1*Z(2)
DZ(3) = k1*(Z(1)-Z(2))
END
FUNCTION 1
F = DLOG(Z(3)+0.0001)
END
EOM