$SIZES      NO=800 LIM6=800                  ; needed for big datasets
$PROB Models for survival and drop out

; Reference: Hansson et. al. (2013) PKPD Modeling of VEGF, sVEGFR-2, sVEGFR-3, and sKIT as Predictors of Tumor Dynamics and Overall Survival Following Sunitinib Treatment in GIST. CPT Pharmacometrics Syst Pharmacol, 20;2:e84.
; For information on imprecision (defined as standard errors) in model parameter estimates kindly refer to the bootstrap results in the published article.


$INPUT
ID                                           ; Patient identification
TIME                                         ; Time in weeks
DAYS                                         ; Time in days
FLAG                                         ;  3. survival ; 10. censoring from survival
DV
CENS                                         ; If FLAG=3:  0. no event ; 1. event (death)
                                             ; If FLAG=10: 0. censored ; 1. not censored
DOS                                          ; dose
PLA                                          ; Placebo: 1. untreated, 0. treated
CL                                           ; posthoc total plasma clearance
EVID
BAS3 MRT3 EC53                               ; posthoc parameters for VEGFR3 timecourse
TBAS                                         ; observed tumor size at baseline

$DATA   Simulated_Survival_GIST.csv
      IGN=#
      IGN=(FLAG.EQ.2)     ; drop out for tumor measurement
      IGN=(FLAG.EQ.4)     ; SLD
      IGN=(FLAG.EQ.5)     ; VEGF
      IGN=(FLAG.EQ.6)     ; VEGFR2
      IGN=(FLAG.EQ.7)     ; VEGFR3
      IGN=(FLAG.EQ.8)     ; SKIT


$SUBROUTINE ADVAN6 TOL=5
$MODEL      NCOMP=3

$PK
; Verbatim code: changes the iteration maximum (IMAX) (default value 100000)
"FIRST
" COMMON /PRCOMG/ IDUM1,IDUM2,IMAX,IDUM4,IDUM5
" INTEGER IDUM1,IDUM2,IMAX,IDUM4,IDUM5
" IMAX=100000000

;-----VEGFR3---------------------------
  BM03    =   BAS3                 ; model predicted parameters for VEGFR3 timecourse
  IMAX1   =   1
  IC503   =   EC53*24*7            ; scale EC50 on AUC from .hour to .week
  MRT33   =   MRT3/24/7            ; scale from hour to week
  KOUT3   =   1/MRT33

;-----Compartment initialization-------
A_0(1)    =   BM03                 ; VEGFR3

AUC=(DOS/CL)*24*7                  ; scale from .hour to .week

LAMBH     =   THETA(1)*EXP(ETA(1)) ; Scale parameter in the Weibull probability density function for the survival model (the ETA is a placeholder here)
ALPHH     =   THETA(2)             ; Shape parameter in the Weibull probability density function for the survival model
LAMBD     =   THETA(5)             ; Scale parameter in the Weibull probability density function for the drop out model
ALPHD     =   THETA(6)             ; Shape parameter in the Weibull probability density function for the drop out model


$DES

;-----VEGFR-3--------------------------
  EFF3    =   IMAX1*AUC/(IC503+AUC)  ; inhibitory Emax drug effect on VEGFR3
  KIN3    =   BM03*KOUT3

DADT(1)   =   KIN3*(1-EFF3)-KOUT3*A(1)     ; VEGFR3 timecourse
BM3       =   (A(1)-BM03)/BM03             ; relative change in VEGFR3 from baseline


;-----Weibull model for survival-------
  DEL=1E-6
  IF(T.GT.4)DADT(2) = LAMBH*ALPHH*(T+DEL)**(ALPHH-1)* EXP(-THETA(3)*BM3+(-THETA(4))*TBAS)
; THETA(3): parameter relating sVEGFR3 to the hazard
; THETA(4): parameter relating   TBAS  to the hazard
; patients who died before the first follow-up (week 4) are not included in the study


;-----Weibull model for dropout
  IF(T.GT.4)DADT(3) = LAMBD*ALPHD*(T+DEL)**(ALPHD-1)


$ERROR

SA1       =   A(1)            ;VEGFR3 timecourse
BBM3      =   (SA1-BM03)/BM03 ; relative change in VEGFR3 from baseline

;-----TTD (Time to death)--------------

  CHZ     =   A(2)           ; cumulative hazard for survival
  DCHZ    =   A(3)           ; cumulative hazard for drop out


  SUR     =   EXP(-CHZ)      ; survival probability
  SURD    =   EXP(-DCHZ)     ; drop out probability
  DELX    =   1E-6

  HAZNOW  =   0
  HAZDNOW =   0
  IF(TIME.GT.4)  HAZNOW  = LAMBH*ALPHH*(TIME+DELX)**(ALPHH-1)*EXP(-THETA(3)*BBM3+(-THETA(4))*TBAS)   ;BASE FOR DROP OUT
  IF(TIME.GT.4)  HAZDNOW = LAMBD*ALPHD*(TIME+DELX)**(ALPHD-1)


  IF(EVID.EQ.0.AND.CENS.EQ.0)     Y = SUR                  ; probability of survival (censored event)
  IF(EVID.EQ.0.AND.CENS.EQ.1)     Y = SUR*HAZNOW           ; probability of event (death) at time=TIME

  IF(FLAG.EQ.10.AND.CENS.EQ.1)    Y = SURD                 ; probability of not dropping out
  IF(FLAG.EQ.10.AND.CENS.EQ.0)    Y = SURD*HAZDNOW         ; probability of dropping out at time=TIME



$THETA  (0,0.00596) ; LAMBH
$THETA  (0,1.226)   ; ALPHH
$THETA  -3.77       ; parameter relating sVEGFR3 to the hazard
$THETA -0.002371    ; parameter relating  TBAS   to the hazard
$THETA (0,0.00196)  ; LAMBD
$THETA (0,1.273)    ; ALPHD

$OMEGA  0  FIX      ; placeholder

;$COV
$EST    PRINT=1 MAXEVAL=0 METHOD=1 LAPLACE LIKE MSFO=MSF3       ; Laplace estimation method: Y is a likelihood
$TABLE  ID TIME DV FLAG EVID CENS SA1 CHZ NOPRINT ONEHEADER FILE=sdtab3