Model changes in DES

In the section of ODE and dosing, we saw that a finite infusion dose changes the ODE system by adding the infusion rate on the RHS. NONMEM provides a general way to extend this concept of system changes at given time, through random variable MTIME.

Consider the Enterohepatic Recycling (EHC) model below ("examples/mtimehill").

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
$PROB  EHC using MTIME for the flag

$INPUT ID DOSE=AMT TIME CP=DV WT SS=DROP II ADDL EVID
$DATA mtimess.dat

$SUBROUTINES  ADVAN6 TOL=6
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
COMP=(PERIPH) COMP=(GALL )

$PK
;# Save the value of II from the dose record.
   if (ii>0) inter=ii

   KA=THETA(1)*EXP(ETA(1))
   KE=THETA(2)*EXP(ETA(2))
   CL=THETA(3)*WT*EXP(ETA(3))
   S2=CL/KE/WT
   K41=THETA(4)*EXP(ETA(4))
   K23=THETA(5)*EXP(ETA(5))
   K32=THETA(6)*EXP(ETA(6))
   K1G=THETA(7)*EXP(ETA(7))

   if (newind<=1) then
     MTIME(1)=THETA(8)
     MTIME(2)=MTIME(1)+THETA(9)
   else
     MTIME(1)=MTIME(1)
     MTIME(2)=MTIME(2)
   endif

  ;# Update mtime variables after mtime(2) is reached.
  ;# New values are for the next dosing interval (TIME+II)
  IF (MNOW == 2) THEN
    MTIME(1) = MTIME(1)+inter
    MTIME(2) = MTIME(2)+inter
    MTDIFF=1
  ENDIF

$DES
   FLAG=MPAST(1)-MPAST(2)
   DADT(1)=-KA*A(1)+K41*A(4)*FLAG
   DADT(2)= KA*A(1)-KE*A(2)-K23*A(2)+K32*A(3)-K1G*A(2)
   DADT(3)= K23*A(2)-K32*A(3)
   DADT(4)= K1G*A(2)-K41*A(4)*FLAG

$ERROR
A1=A(1) ;# for display in table
A2=A(2)
A3=A(3)
A4=A(4)
   Y=F+EPS(1)

$THETA 1 1 1
$THETA 10  ;# rate of gall bladder emptying is large vs. other k's
$THETA 1 1 1
$THETA 3 5 ;# start emptying at T+theta(8) ; lasts till T+theta(9)
$OMEGA 1 1 1 1 1 1 1
$SIGMA  .4

$TABLE TIME A1 A2 A3 A4 NOAPPEND FILE=mtimess.tab NOPRINT

The model has drug transfer from compartment 4 to 1 that only appears during a certain period of time (to reflect EHC). To estimate the beginning and end of this period, we model them as MTIME(1) and MTIME(2), respectively. Then values of the internal indicator variable MNEXT, MPAST, and MNOW is shown in the follow scheme.

1
2
3
4
5
6
7
           MTIME(1)       MTIME(2)
  ............|..............|...........
  01111111111110000000000000000000000000 MNEXT(1)
  00000000000001111111111111110000000000 MNEXT(2)
  00000000000001111111111111111111111111 MPAST(1)
  00000000000000000000000000001111111111 MPAST(2)
  01111111111112222222222222220000000000 MNOW

In other words, MNEXT(i) and MPAST(i) flag the "before" and "after" of MTIME(i), and MNOW stores the counter of MNEXT. Therefore, we define FLAG on line 40 such that FLAG=1 during time interval (MTIME(1), MTIME(2)], indicating the transfer between the two compartments.

Note that MTIME variables are ignored during Steady-State calculation. This is because steady states are intended for a single repeated dose with no changes in the status of the system. Hence the $INPUT... SS=DROP option.

Smooth changes of DES

In the above example, the derivative DADT is discontinuous at MTIME(1) and MTIME(2), i.e., the changes to the ODEs is instataneous (though the states A stays continuous). To make DADT continuous, we can use a smooth function that changes between 0 and 1 as FLAG. Consider the following alternative implementation of the EHC model ("examples/mtimehill").

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
$PROB  EHC using Hill (sigmoid emax model) for the flags
$INPUT      ID DOSE=AMT TIME CP=DV WT SS II ADDL EVID
$DATA       hillss.dat

$SUBROUTINES  ADVAN6 TOL=6
$MODEL COMP=(DEPOT,INITIALOFF,DEFDOSE) COMP=(CENTRAL,DEFOBS,NOOFF)
COMP=(PERIPH) COMP=(GALL )

$PK
   KA=THETA(1)*EXP(ETA(1))
   KE=THETA(2)*EXP(ETA(2))
   CL=THETA(3)*WT*EXP(ETA(3))
   S2=CL/KE/WT
   K41=THETA(4)*EXP(ETA(4))
   K23=THETA(5)*EXP(ETA(5))
   K32=THETA(6)*EXP(ETA(6))
   K1G=THETA(7)*EXP(ETA(7))

$DES
;# Save the value of II from the dose record.
   if (ii>0) inter=ii
   mt1=inter*INT(T/inter)+THETA(8)
   mt2=mt1+THETA(9)
   hill1=exp(-THETA(10)*(T-mt1))
   hill2=exp(-THETA(10)*(T-mt2))
   flag1=1./(1+hill1)  ;# changes from 0 to 1 near t=mt1
   flag2=1./(1+hill2)  ;# changes from 0 to 1 near t=mt2
   FLAG=flag1-flag2

   DADT(1)=-KA*A(1)+K41*A(4)*FLAG
   DADT(2)= KA*A(1)-KE*A(2)-K23*A(2)+K32*A(3)-K1G*A(2)
   DADT(3)= K23*A(2)-K32*A(3)
   DADT(4)= K1G*A(2)-K41*A(4)*FLAG

$ERROR
A1=A(1) ;# for display in table
A2=A(2)
A3=A(3)
A4=A(4)
   Y=F+EPS(1)

$THETA 1 1 1
$THETA 10  ;# rate of gall bladder emptying is large vs. other k's
$THETA 1 1 1
$THETA 3 5 ;# start emptying at T+THETA(8) ; lasts till T+THETA(9)
$THETA 50 ;# a very large hill exponent
$OMEGA 1 1 1 1 1 1 1
$SIGMA  .4

$TABLE TIME A1 A2 A3 A4 NOAPPEND FILE=hillss.tab NOPRINT

Here we construct two Hill functions of the independent time variable T that change their values in the neighborhood of time mt1 and mt2, and use their differences as the smooth FLAG that rises from 0 to 1 around mt1 and back to 0 around mt2. THETA(10), the exponent in the Hill functions, represents how fast the derivatives change. The greater it is the closer the system approache the earlier model using MTIME for instataneous changes.

Note unlike the model using MTIME, stead-state dose is allowed here. But keep in mind the changes in differential equations do not act retroactively: the action of FLAG in the future does not affect the steady-state calculations. Specifically, with data (the first a few records of the data file for "examples/mtimehill")

1
2
3
4
5
6
7
8
#       ID DOSE=AMT      TIME        CP=DV     WT SS II ADDL EVID
         1    100.00      0.         .       79.6 1  12 100   1
         1       .        4.0         .      79.6 .  . .     2
         1       .        8.0         .      79.6 .  . .     2
         1       .       12.0         .      79.6 .  . .     2
         1       .       16.0         .      79.6 .  . .     2
         1       .       20.0         .      79.6 .  . .     2
         1       .       24.0         .      79.6 .  . .     2

that specifies II=12, the $DES changes at THETA(8) and THETA(8)+THETA(9) after the the start of the begin of intervals (0, 12), (12, 24), etc, when the system is steady-state. This is done by using INT (integer function)

1
mt1=II*INT(T/II)+THETA(8)

Periodic discontinuities in $DES

Changes in $DES could be periodic. One classical example is circadian rhythm. Consider the following model that has a change in $DES every 24 hours ("examples/step_circexa").

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
;; Use mtime to model periodic discontinuouities in $DES,
;; suggested by Sebastien Bihorel, Luann Phillips and Jill Fiedler-Kelly.
$PROBLEM step_circadian
;; -----------------------------------------------------------------------------
;; PURPOSE:
;; 1- Illustrate how to turn a step function on at a certain time of the day,
;;    SHIFT, every 24h and back off after a certain duration, DUR, using model
;;    event (MTIME) variables.
;;
;;            ^
;;         1 -|  ---------      ---------      ---------
;;            |  |       |      |       |      |       |
;;            |  |       |      |       |      |       |
;;         0 -|---       --------       --------       -------- FLAG
;;            |
;;            ---------------------------------------------------> TIME
;;               |   |              |              |
;;        SHIFT-24   0              24             48
;;                   <-  SHIFT -><- DUR ->
;;
;;    FLAG is used as the rate of the differential equation dA(1)/dt integrated
;;    over TIME. A(1) increases by the amount DUR every day.
;;
;; 2- Illustrate how NONMEM updates MTIME, MPAST, and MNOW variables over TIME.
;;    $TABLE will only report the value of MTIME, MPAST, MNOW, and any variable
;;    dependent on them at the TIME of event records. Therefore, WRITE
;;    statements are implemented to document the update of these variables at
;;    and in between event records.
;;
;; -----------------------------------------------------------------------------

$INPUT ID TIME AMT CMT EVID DV MDV

$DATA circadian.csv IGNORE=@

$THETA
 (0.0001,0.5,0.9999)   ;-- typical shift as a fraction of 24h
 (0.0001,0.5,0.9999)   ;-- typical duration as a fraction of 24h

$OMEGA
  3                    ;-- IIV in shift (%CV=100*sqrt(1-th1)*eta1)
  3                    ;-- IIV in duration (%CV=100*sqrt(1-th2)*eta2)

$SUBROUTINE ADVAN13 TOL=9

$MODEL NCOMP=1
  COMP=(FLAG)

$PK
  REPID=(IREP-1)*10 + ID

  ; Define shift and duration parameters

  TVSHIFT = 24*THETA(1)                   ; typical shift on 24h scale; (0,24)
                                          ; where THETA(1) represents this typical shift as a fraction of 24h; (0,1)
  LGSHI = LOG(THETA(1)/(1-THETA(1)))      ; logit transformation of the typical fraction of 24h; (-INF,+INF)
  ILGSHI = LGSHI+ETA(1)                   ; addition of IIV; (-INF,+INF)
  SHIFT = 24*EXP(ILGSHI)/(1+EXP(ILGSHI))  ; back-transformation and scale-up of individual shift to 24h scale; (0,24)

  TVDUR = 24*THETA(2)                     ; typical duration on 24h scale; (0,24)
                                          ; where THETA(1) represents this typical duration as a fraction of 24h; (0,1)
  LGDUR = LOG(THETA(2)/(1-THETA(2)))      ; logit transformation of the typical fraction of 24h; (-INF,+INF)
  ILGDUR = LGDUR+ETA(2)                   ; addition of IIV; (-INF,+INF)
  DUR = 24*EXP(ILGDUR)/(1+EXP(ILGDUR))    ; back-transformation and scale-up of individual duration to 24h scale; (0,24)

  ; Define offset for the start of the circadian rhythm parameter (in case the
  ; step function "starts" before TIME=0
  IF ((SHIFT+DUR) > 24) THEN
    OFFSET = 24
  ELSE
    OFFSET = 0
  END IF

  TX = SHIFT-OFFSET

  IF (NEWIND <= 1) THEN
    MTIME(1) = TX
    MTIME(2) = TX+DUR
  ELSE
    MTIME(1) = MTIME(1)
    MTIME(2) = MTIME(2)
  ENDIF

  ; Defines of mtime(1) to use for integration during this interval
  INTMTIME = MTIME(1)

  ; Update mtime variables after mtime(2) is reached
  IF (MNOW == 2) THEN
    MTIME(1) = MTIME(1)+24
    MTIME(2) = MTIME(2)+24
    MTDIFF=1
  ENDIF

  ; Save mtime(1), mtime(2), and tstate for the TABLE statement
  MT1=MTIME(1)
  MT2=MTIME(2)
  MP1=MPAST(1)
  MP2=MPAST(2)
  STS=TSTATE

  ; Save mnow as a real variable for the WRITE statement
  RMNOW=MNOW

  ; Initialize the file used to output variables at and in-between observation
  ; records
  IF (ICALL == 1) THEN
    "OPEN (FILE='step_circexa.txt',UNIT=96)
    "WRITE (96,'(10A14)') 'REPID','RMNOW','TIME','MT1','MT2','MP1','MP2','TSTATE','INTMTIME'
  ENDIF

  ; Output variables in-between observation records using WRITE statement
  WRITE (96,991) REPID,RMNOW,TIME,MT1,MT2,MP1,MP2,TSTATE,INTMTIME

$DES
  ; Define FLAG variable:
  ; FLAG=1 when T is in [SHIFT+n*24, SHIFT+DUR+n*24], where n=0,1,2,3,...
  ; FLAG=0 otherwise
  FLAG = MPAST(1)-MPAST(2)

  ; Monitor FLAG:
  DADT(1) = FLAG

$ERROR
  ; Output variables at observation records using WRITE statement and flagging
  ; the output row with MNOW=9.
  ; Uncomment the following 3 lines of verbatim code when using nm7.1 or nm7.2,
  ; as TSTATE is not a reserved variable in $ERROR prior to nm7.3
    "FIRST
    "USE PROCM_REAL, ONLY: TSTATE
    "WRITE (96,991) REPID,9.0,TIME,MT1,MT2,MP1,MP2,TSTATE,INTMTIME
  ;
  ; Uncomment the following line of code when using nm7.3
  ;  WRITE (96,991) REPID,9.0,TIME,MT1,MT2,MP1,MP2,TSTATE,INTMTIME

  ; Define ouput for TABLE statement
  FLAGCMT = A(1)
  IF (ICALL == 4) THEN
    DV = FLAGCMT
  ENDIF

$SIMULATION(12345) NSUB=10 ONLYSIM NOPREDICTION

$TABLE REPID TIME FLAG FLAGCMT NOAPPEND NOPRINT ONEHEADER FILE=step_circexa.tab

$TABLE REPID RMNOW TIME MT1 MT2 MP1 MP2 STS INTMTIME SHIFT DUR
  NOAPPEND NOPRINT ONEHEADER FILE=step_circexa_debug.tab
  ; Uncomment the following line of code when using nm7.1 or nm7.2
  FORMAT=s34F13.4
  ; Uncomment the following line of code when using nm7.3
  ;RFORMAT="(10F10.4)"  LFORMAT="(3x,10A10)"

$SCAT FLAG VS TIME BY REPID OBSONLY

Simlar to the first model, FLAG is the difference of MPAST variables. The MTIME(1) and MTIME(2) variables respectively define the times at which FLAG is set to 1 and 0, such as MTIME(2)-MTIME(1)=DUR. After MTIME(2) is reached, both variables are incremented by 24, and the FLAG update process perpetuates itself every 24 hours. Compartment 1 is a dummy compartment intended to monitor the FLAG variable: the amount A1 starts at 0 and is incremented by DUR every day. A related information is the time at which the current state vector was computed, i.e. the time to which the system was most recently advanced. It is reported in the reserved variable TSTATE.

Both $TABLE and WRITE are used in to report the values of the variables of interest. WRITE statements are used because $TABLE functions only for event records, but not for nonevent times such as Model times.

For an extension of the above example, consider the following biological response

\begin{equation} R(t) = R_{\text{min}} + a\sin(z(t)), \quad t\in [t_{\text{shift}}+24n, t_{\text{shift}}+t_{\text{dur}} + 24n]. ) \end{equation}

Here \(R_{\text{min}}\) is the flat baseline, \(n\) the counter, \(z(t)\) the circadian function of time \(t\) that scales the sinusoidal response interval to \([0, \pi]\). This model can be used to describe the response baseline of a drug that influences circadian rhythms related biological processes, e.g. metabolic effects of corticoids. If the effect of the drug is not direct, one can express the biological response and the drug effect in terms of an indirect response model. For example, assume an agonist that stimulates the formation of the said biological processes, we can model the response as

$$ \frac{dR}{dt} = k_{\text{in}}(t)(1+s(C)) - k_{\text{out}}R. $$ where \(k_{\text{out}}\) is the constant elimination rate of the response, \(k_{\text{in}}\) the time-varying formation rate, and \(s(C)\) the stimulation that depends on drug concentration \(C\). The Eq.(1) then describes the response baseline in the absence of drug, $$ \frac{dR}{dt} = k_{\text{in}}(t), $$ where \(K_{\text{in}}\) is sinusoidal inside certain intervals and constant otherwise. Again we use a FLAG variable to model it ("examples/idr_circexa").

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
; Use mtime to model periodic discontinuouities in $DES,
; suggested by Sebastien Bihorel, Luann Phillips and Jill Fiedler-Kelly.
$PROBLEM idr_circadian
;; -----------------------------------------------------------------------------
;; PURPOSE:
;; 1- Illustrate an indirect response model describing a baseline response with
;;    circadian rhythm. The model implements a time-varying kin(t)=dR/dt+Kout*R(t)
;;    where R(t) is a modified sine function. Every 24h, this function increases
;;    from a minimum value, Rmin, at a certain time, SHIFT, reaches a maximum,
;;    Rmin+Amp, then decreases to return to its minimum at time SHIFT+DUR
;;
;;            ^
;;  Rmin+Amp -|     ---            ---            ---
;;            |    |   |          |   |          |   |
;;            |   |     |        |     |        |     |
;;            |  |       |      |       |      |       |
;;  Rmin     -|---       --------       --------       ------ Response R
;;            |
;;            ---------------------------------------------------> TIME
;;               |   |              |              |
;;        SHIFT-24   0              24             48
;;                   <-  SHIFT -><- DUR ->
;;
;;    Response must oscillate between 0.5 and 4.5
;;
;; 2- Illustrate how to create a 24h clock time variable T24 in $DES by
;;    transforming the T variable ranging from 0 to infinity into repeating
;;    intervals of continuous time between 0 and 24.
;;
;; -----------------------------------------------------------------------------

$INPUT ID TIME AMT CMT EVID DV MDV

$DATA circadian.csv IGNORE=@

$THETA
 (0.0001,0.5,0.9999)   ;-- typical scaled-down shift
 (0.0001,0.5,0.9999)   ;-- typical scaled-down duration
 (0,4)                 ;-- typical amplitude of response
 (0,0.5)               ;-- typical minimum response
 (0,0.1)               ;-- typical elimination rate

$OMEGA
  3                    ;-- IIV in shift (%CV=100*sqrt(1-th1)*eta1)
  3                    ;-- IIV in duration (%CV=100*sqrt(1-th2)*eta2)

$SUBROUTINE ADVAN13 TOL=9

$MODEL NCOMP=1
  COMP=(RESP)

$PK

  ;------------------------
  ; Indirect response model
  ;------------------------
  PI = 3.14159263

  ; Define shift and duration parameters

  TVSHIFT = 24*THETA(1)                   ; typical shift on 24h scale; (0,24)
                                          ; where THETA(1) represents this typical shift as a fraction of 24h; (0,1)
  LGSHI = LOG(THETA(1)/(1-THETA(1)))      ; logit transformation of the typical fraction of 24h; (-INF,+INF)
  ILGSHI = LGSHI+ETA(1)                   ; addition of IIV; (-INF,+INF)
  SHIFT = 24*EXP(ILGSHI)/(1+EXP(ILGSHI))  ; back-transformation and scale-up of individual shift to 24h scale; (0,24)

  TVDUR = 24*THETA(2)                     ; typical duration on 24h scale; (0,24)
                                          ; where THETA(1) represents this typical duration as a fraction of 24h; (0,1)
  LGDUR = LOG(THETA(2)/(1-THETA(2)))      ; logit transformation of the typical fraction of 24h; (-INF,+INF)
  ILGDUR = LGDUR+ETA(2)                   ; addition of IIV; (-INF,+INF)
  DUR = 24*EXP(ILGDUR)/(1+EXP(ILGDUR))    ; back-transformation and scale-up of individual duration to 24h scale; (0,24)

  AMP = THETA(3)
  RMIN = THETA(4)
  KOUT = THETA(5)
  RMAX = RMIN + AMP

  ; Circadian rhythm
  ; Determine initial value of the response and time offset based upon whether
  ; SHIFT and DUR are such that the response is flat or within the sinusoid
  ; phase at time 0
  IF ((SHIFT+DUR).GT.24) THEN
    OFFSET = 24
    INIT1 = RMIN+AMP*SIN(PI*(OFFSET-SHIFT)/DUR)
  ELSE
    OFFSET = 0
    INIT1 = RMIN
  END IF
  TX = SHIFT-OFFSET

  IF (NEWIND <= 1) THEN
    MTIME(1) = TX
    MTIME(2) = TX+DUR
  ELSE
    MTIME(1) = MTIME(1)
    MTIME(2) = MTIME(2)
  ENDIF

  ; Defines of mtime(1) to use for integration during this interval
  INTMTIME = MTIME(1)

  ; Update mtime variables after mtime(2) is reached
  IF (MNOW == 2) THEN
    MTIME(1) = MTIME(1)+24
    MTIME(2) = MTIME(2)+24
    MTDIFF=1
  ENDIF

  ; Initialize the response compartment
  IF (A_0FLG == 1) THEN
    A_0(1) = INIT1
  ENDIF

  ;---------------
  ; 24h clock time
  ;---------------

  ; Use mtime(3) to scale the T variable in $DES from 0 to infinity to repeating
  ; 0-24h intervals
  IF (NEWIND<=1) THEN
    MTIME(3)=TIME+24
  ELSE
    MTIME(3)=MTIME(3)
  ENDIF

  INTMTIME3=MTIME(3)

  IF (MNOW==3) THEN
    MTIME(3)=MTIME(3)+24
    MTDIFF=1
  ENDIF

$DES
  ; Define FLAG variable:
  ; FLAG=1 when T is in [SHIFT+n*24,SHIFT+DUR+n*24], where n=0,1,2,3,...
  ; FLAG=0 otherwise
  FLAG = MPAST(1)-MPAST(2)

  ; Define time-varying kin:
  ; * kin is constant (KOUT*RMIN) when FLAG=0
  ; * kin is time-varying when FLAG=1, and function of TS, a transformation
  ;   function of T scaling time from 0 to infinity to repeating intervals [0,1]
  TS = (T-INTMTIME)/DUR
  KIN = KOUT*RMIN + FLAG*((AMP*PI/DUR)*COS(PI*TS) + KOUT*AMP*SIN(PI*TS))

  ; Response
  DADT(1) = KIN-KOUT*A(1)

  ; 24h clock time
  T24 = T-(INTMTIME3-24)

$ERROR
  REPID = (IREP-1)*10 + ID
  RESPONSE = A(1)
  IF (ICALL == 4) THEN
    DV = RESPONSE
  ENDIF

$SIMULATION(12345) NSUB=10 ONLYSIM NOPREDICTION

$TABLE REPID TIME RESPONSE SHIFT DUR ONEHEADER NOAPPEND NOPRINT
  FILE=idr_circexa.tab

$SCAT FLAG VS TIME BY REPID
$SCAT RESPONSE VS TIME BY REPID

The $DES block describes the response as the state variable. Line 120 defines a continuous 24h clock time, through the circadian reset of the MTIME(3) variable in $PK. Note that here \(z(t)\) is implemented using the INTMTIME variable that records MTIME(1).