Lili,
you create a discontinuity in the first derivative and that can be
difficult to handle numerically.
Other than that, nonmem's failure to execute the covariance step is not
necessarily a sign of anything.
I assume that your two halves, Tmin and Tmax, are not 12 hours apart and
that is why you do not use a standard trigonometric function (sine or
cosine) for modeling the circadian rhythm.
Here is an alternative: Tweak the standard cosine function to split into
two parts of different lengths.
Basically you stretch the x-axis differently in the two parts.
L2=24-L1 ; L1=length of first period, L2=length of second
period
T24 = (TIME-CSHIFT) - INT((TIME-CSHIFT)/24); modulo 24
T1=T24/L1*12 ; part 1: map from (0, L1) to (0, 12)
T2=(T24-L1)/L2*12+12 ; part 2: map from (L1, L2) to (12, 24)
IF(T24.LT.L1)THEN
TX=T1
ELSE
TX=T2
ENDIF ; step-function but at L1, T1=T2=12 and gradient=0
PI=3.141592653
CIRC = CAMP*COS(TX/24*2*PI)
Now TMAX equals CSHIFT and TMIN equals CSHIFT+L1, the length of the first
half-period.
Here is an R function for illustration:
f <- function(
camp=10, # circadian amplitude
cshift=3, # circadian shift, time of occurrence of maximum
time=seq(0, 100, length=100), # time axis
L1=5 # length of first interval or time of occurrence of the
minimum (time after cshift)
)
{
L2 <- 24-L1 # derived length of second period
# Split time axis from (0, L1) and (L1, L2) to (0, 12) and (12,
24),
# then stretch the two intervals individually (note that time gets
stretched)
t24 <- (time-cshift)%%24 # remainder of division by 24
t1 <- t24/L1*12 # part 1 from (0, L1) to (0, 12)
t2 <- (t24-L1)/L2*12 + 12 # part 2 from (L1, L2) to (12, 24)
tx <- ifelse(t24 < L1, t1, t2) # step-function but at L1,
t1=t2=12 and gradient 0
y <- camp*cos(tx/24*2*pi)
# output some of the data and a graph
data <- cbind(time, t24, t1, t2, tx, y)
print(data[1:min(25, nrow(data)), ])
str <- paste("L1=", L1, ", L2=", L2, ", cshift=", cshift, ",
camp=", camp, sep="")
plot (time, y, type="l", lwd=5, main=str)
abline(v=seq( cshift, max(time), by=24), col="darkblue") # daily
max
abline(v=seq(L1+cshift, max(time), by=24), col="blue") # daily
min
}
# Call the function:
f() # default arguments
f(cshift=10, L1=12) # standard with equal lengths
f(cshift=10, L1=6)
f(cshift=10, L1=6, time=seq(0, 25, length=1000)) # zoom in
-----
Andreas Krause, PhD
Director, Lead Scientist Modeling and Simulation
Dept. of Clinical Pharmacology
Actelion Pharmaceuticals Ltd, Hegenheimermattweg 91, CH-4123 Allschwil,
Switzerland
From:
Li Li <[email protected]>
To:
[email protected]
Date:
25.05.2011 19:57
Subject:
[NMusers] Coding INTEGER Function in NONMEM
Sent by:
[email protected]
Hi,
I would like to build a model to describe a circadian rhythm. Basically I
assume that the release rate is linearly increasing from Tmin to Tmax, and
linearly decreasing from Tmax to Tmin in 24 hr period. My code is as
follows (I tried to put integer operation in $DES, to convert every time
point after 24hr to 0-24hr period). My run failed at the covariance step,
and nonmem gave me a warning:
(WARNING 68) THE INTEGER FUNCTION IS BEING USED OUTSIDE OF A SIMULATION
BLOCK. IF THE INTEGER VALUE AFFECTS THE VALUE OF THE OBJECTIVE FUNCTION,
THEN AN ERROR WILL PROBABLY OCCUR.
How the warning will affect my nonmem output? Why the run failed? Any
suggestions on the coding will be highly appreciated. Thank you.
Best wishes,
lili
NONMEM CODE:
;PD Model
KOUT = 0.4 ;Elimination rate constant (1/hr)
RMAX = THETA(1)*EXP(ETA(1)) ;Maximum release rate (ng/ml/hr)
TMIN = THETA(2)*EXP(ETA(2)) ;Time at minimum release rate (hr)
TGAP = THETA(3)*EXP(ETA(3)) ;Time from minimum rate to maximum release
rate (hr)
BSL = THETA(4)*EXP(ETA(4)) ;Concentration at time zero (ng/ml)
SD = THETA(5) ;CV% of residual
TMAX = TMIN+TGAP
A_0(1) = BSL
$DES
T1= T-24*INT(T/24) ;Convert time to 0-24 period, the MOD function is not
working in NONMEM
IF (T1.LE.TMIN) THEN ;Make the time start from Tmin
T2=T1+24
ELSE
T2=T1
ENDIF
IF (T2.LE.TMAX) THEN
KIN=RMAX*(T2-TMIN)/TGAP
ELSE
KIN=RMAX*(TMIN+24-T2)/(24-TGAP)
ENDIF
DADT(1)= KIN-KOUT*A(1)
$ERROR
IPRED=F
Y=F*(1+SD*EPS(1))
IRES=DV-IPRED
IWRES=IRES/SD/F
The information of this email and in any file transmitted with it is strictly
confidential and may be legally privileged.
It is intended solely for the addressee. If you are not the intended recipient,
any copying, distribution or any other use of this email is prohibited and may
be unlawful. In such case, you should please notify the sender immediately and
destroy this email.
The content of this email is not legally binding unless confirmed by letter.
Any views expressed in this message are those of the individual sender, except
where the message states otherwise and the sender is authorised to state them
to be the views of the sender's company. For further information about Actelion
please see our website at http://www.actelion.com
