Design Evaluation when No. of subjects change towards the end of trial

[Anks2030]

Hi @andrewhooker
I am trying to optimize a PIII trial design using popED.
This trial is for one year and the number of subjects stays same (N=200) under week 52, but post-wk-52 the extensive PK is collected only for N=20 subjects
Before wk-52, trough samples are only collected
Do you know if this situation can be implemented in popED for design evaluation and optimization
Will appreciate your response
Thanks,
-Ankur

[andrewhooker]

Hi Ankur,

This can certainly be implemented! You just need to define different groups with different designs in the trial. See for example https://andrewhooker.github.io/PopED/articles/intro-poped.html. You can post example code here if you like and I can try and help.

Best regards,
Andy

[Anks2030]

Hi Andrew,
Thanks for looking into this.
In our trial, we have N=200 in treatment arm until wk-52. Post-52 week we have PK at D1, D2, D3, D4 and D7, D14. And post-52WK the collection is for only N=20.
I dont know how to implement this switch of N=200 to N=20 and have an overall influence of design?
So far I am just assuming 200 for entire duration and have even tried to evaluate the design only till WK-52 (using just N=200)
Please have code for entire study for N=200

#' Define the ODE system
PK.2.comp.sc.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {

CL=CL*(WT/76.8)^(WT_CL)
V1=V1*(WT/76.8)^(WT_V1)

dA1 <- -KA*A1 
dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
dA3 <- A2* Q/V1-A3* Q/V2
return(list(c(dA1, dA2, dA3)))

})
}

#' define the initial conditions and the dosing
ff.PK.2.comp.sc.md.ode <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0, A3=0)
times_xt <- drop(xt)
dose_times = c(seq(from=0,to=336, by=84),seq(from=336,to=8568,by=TAU))
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE), method = c("add"))
times <- sort(c(times_xt,dose_times))
times <- unique(times) # remove duplicates
out <- ode(A_ini, times, PK.2.comp.sc.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = (out[, "A2"]/(V1/Favail))+BL
y=y[match(times_xt,out[,"time"])]
y=cbind(y)
return(list(y=y,poped.db=poped.db))
})
}

#' parameter definition function
#' names match parameters in function ff
fg <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]),
V1=bpop[2]*exp(b[2]),
KA=bpop[3],
Q= bpop[4],
V2=bpop[5],
Favail=bpop[6],
WT_CL=bpop[7],
WT_V1=bpop[8],
BL=bpop[9],
DOSE=a[1],
TAU=a[2],
WT= a[3])
return( parameters )
}

-- Residual unexplained variablity (RUV) function

-- Additive + Proportional

feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]

y = y*(1+epsi[,1])+epsi[,2]

return(list( y= y,poped.db =poped.db ))
}

Original Design

#CL: L/h; V1,V2: L; Q: L/h; V2: L; BL: g/L, KA= h-1
#create poped database original design
poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.sc.md.ode",
fg_fun="fg",
fError_file="feps",

                              bpop=c(CL=0.00633,V1=8.75,KA=0.02108,Q= 0.0344, V2= 1.8, Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
                              notfixed_bpop=c(1,1,1,0,1,0,0,0,0),# decides which tested and which are not

                              d=c(CL=0.282,V1=0.569), # decides the variances of BSV
                              notfixed_d = c(1,1),

                              sigma=c(0.00662,0), # decides the variances of RV prop, additive
                              notfixed_sigma=c(1,0),

                              m=1,      #number of groups
                              groupsize=200,


                              #0,7,28,35,56,59,63,70,77,84
                              xt=   c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904),
                              minxt=c(0,312,648,1992,2664,3336,6024,7368,8544,8568,8592,8616,8640,8712,8880),
                              maxxt=c(0,360,696,2040,2712,3384,6072,7416,8592,8616,8640,8664,8688,8760,8928),

                              bUseGrouped_xt=1,
                              a=c(DOSE=0.15*76.8,TAU=168, WT=76.8))

#' plot intial design with BSV and RUV in model
#plot_model_prediction(poped.db,IPRED=T,DV=T)
dat_Original <- model_prediction(poped.db)
plot_model_prediction(poped.db,DV=T,
# separate.groups=T,, groupsize_sim=100,PI=T
sample.times=F,alpha.sample.times.DV.points = 0.05,
model_num_points = 1000)+# , PI=T
labs(x = "Time from first dose (Weeks)", y="IgG Conc.(g/L)")+
theme_bw()+ geom_vline(xintercept=336,lty="dashed")+
geom_point(data = dat_Original, aes(Time, PRED ), size=2, alpha=1/2, colour="blue")+#colour="Before Optimization"
ggtitle("Original Design")+coord_cartesian(ylim=c(4,12))+
scale_y_continuous(breaks=sort(c(seq(0,12,1))))+
scale_x_continuous(breaks = seq(0,8928,168), labels=c(seq(1,54,1)))+
theme_bw()

Evaluate Original Design

(dso <- evaluate_design(poped.db))

[andrewhooker]

Hi Ankur,

Here is am implementation of what you are trying to do.

First define the model:

library(tidyverse)
library(deSolve)
library(PopED)
packageVersion("PopED") #0.5.0

#' Define the ODE system
PK.2.comp.sc.ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {
    
    
    
    CL=CL*(WT/76.8)^(WT_CL)
    V1=V1*(WT/76.8)^(WT_V1)
    
    dA1 <- -KA*A1
    dA2 <- KA*A1 + A3* Q/V2 -A2*(CL/V1+Q/V1)
    dA3 <- A2* Q/V1-A3* Q/V2
    return(list(c(dA1, dA2, dA3)))
  })
}

#' define the initial conditions and the dosing
ff.PK.2.comp.sc.md.ode <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0)
    times_xt <- drop(xt)
    dose_times = c(seq(from=0,to=336, by=84),seq(from=336,to=8568,by=TAU))
    eventdat <- data.frame(var = c("A1"),
                           time = dose_times,
                           value = c(DOSE), method = c("add"))
    times <- sort(c(times_xt,dose_times))
    times <- unique(times) # remove duplicates
    out <- ode(A_ini, times, PK.2.comp.sc.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = (out[, "A2"]/(V1/Favail))+BL
    y=y[match(times_xt,out[,"time"])]
    y=cbind(y)
    return(list(y=y,poped.db=poped.db))
  })
}

#' parameter definition function
#' names match parameters in function ff
fg <- function(x,a,bpop,b,bocc){
  parameters=c( CL=bpop[1]*exp(b[1]),
                V1=bpop[2]*exp(b[2]),
                KA=bpop[3],
                Q= bpop[4],
                V2=bpop[5],
                Favail=bpop[6],
                WT_CL=bpop[7],
                WT_V1=bpop[8],
                BL=bpop[9],
                DOSE=a[1],
                TAU=a[2],
                WT=  a[3])
  return( parameters )
}

## -- Residual unexplained variablity (RUV) function
## -- Additive + Proportional 
feps <- function(model_switch,xt,parameters,epsi,poped.db){
  returnArgs <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))
  y <- returnArgs[[1]]
  poped.db <- returnArgs[[2]]
  
  y = y*(1+epsi[,1])+epsi[,2]
  
  return(list( y= y,poped.db =poped.db ))
}

Add a design with design with extra observations after 52 wks in 20 patients.

poped_db_1 <- create.poped.database(
  ff_fun="ff.PK.2.comp.sc.md.ode",
  fg_fun="fg",
  fError_file="feps",
  
  bpop=c(CL=0.00633,V1=8.75,KA=0.02108,
         Q= 0.0344, 
         V2= 1.8, 
         Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
  notfixed_bpop=c(1,1,1,
                  0,
                  1,
                  0,0,0,0),# decides which tested and which are not
  
  d=c(CL=0.282,V1=0.569), # decides the variances of BSV
  notfixed_d = c(1,1),
  
  sigma=c(0.00662,0), # decides the variances of RV prop, additive
  notfixed_sigma=c(1,0),
  
  m=2,      #number of groups
  groupsize=c(20,180),
  
  
  #0,7,28,35,56,59,63,70,77,84
  # Here we have 20 individuals that have sampling like the rest until wk 52 and
  # then more samples post wk52.  
  # 180 have only sampling until wk 52.
  xt=   list(c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904),
             c(0,336,671,2016,2688,3360,6048,7392,8568)), 
  minxt=list(c(0,312,648,1992,2664,3336,6024,7368,8544,8568,8592,8616,8640,8712,8880),
             c(0,312,648,1992,2664,3336,6024,7368,8544)),
  maxxt=list(c(0,360,696,2040,2712,3384,6072,7416,8592,8616,8640,8664,8688,8760,8928),
             c(0,360,696,2040,2712,3384,6072,7416,8592)),
  
  #bUseGrouped_xt=1,
  a=list(c(DOSE=0.15*76.8,TAU=168, WT=76.8),
         c(DOSE=0.15*76.8,TAU=168, WT=76.8))
)

Evaluate design and plot design.

(p1 <- plot_model_prediction(poped_db_1,separate.groups = T,model_num_points = 1000))
(des1 <- evaluate_design(poped_db_1))

Create a design without the extra observations after 52 weeks:

poped_db_2 <- create.poped.database(
  ff_fun="ff.PK.2.comp.sc.md.ode",
  fg_fun="fg",
  fError_file="feps",
  
  bpop=c(CL=0.00633,V1=8.75,KA=0.02108,
         Q= 0.0344, 
         V2= 1.8, 
         Favail=0.668,WT_CL=0.796, WT_V1=1.1, BL=4),# TV 0.007475 L/h, L
  notfixed_bpop=c(1,1,1,
                  0,
                  1,
                  0,0,0,0),# decides which tested and which are not
  
  d=c(CL=0.282,V1=0.569), # decides the variances of BSV
  notfixed_d = c(1,1),
  
  sigma=c(0.00662,0), # decides the variances of RV prop, additive
  notfixed_sigma=c(1,0),
  
  m=1,      #number of groups
  groupsize=c(200),
  
  
  #0,7,28,35,56,59,63,70,77,84
  # Here we have 20 individuals that have sampling like the rest until wk 52 and
  # then more samples post wk52.  
  # 180 have only sampling until wk 52.
  xt= c(0,336,671,2016,2688,3360,6048,7392,8568), 
  minxt=c(0,312,648,1992,2664,3336,6024,7368,8544),
  maxxt=c(0,360,696,2040,2712,3384,6072,7416,8592),
  
  #bUseGrouped_xt=1,
  a=c(DOSE=0.15*76.8,TAU=168, WT=76.8)
)

Evaluate design and plot design.

(p2 <- plot_model_prediction(poped_db_2,separate.groups = T,model_num_points = 1000))
(des2 <- evaluate_design(poped_db_2))

There is a big difference in parameter RSE between the two designs

tibble::tibble(names(des1$rse),des1$rse,des2$rse)

There is also a big difference in efficiency (52% more individuals in the design without extra observations after 52 wks are needed to achieve same information as the extended design)

efficiency(ofv_init = des2$ofv, ofv_final = des1$ofv, poped_db = poped_db_1)

[Anks2030]

Hi Andy,
Thanks for taking time to go over this challenge
I have a question.
I am trying hard to understand if the code is capturing the design appropriately

I am attaching a single design figure (total duration is 1year) in which we have N=180 only till WK-52 and after WK-52 we have ONLY N=20

I trying to understand if this design is applied in below code you shared?

m=2, #number of groups
groupsize=c(20,180),

#0,7,28,35,56,59,63,70,77,84

Here we have 20 individuals that have sampling like the rest until wk 52 and

then more samples post wk52.

180 have only sampling until wk 52.

xt= list(c(0,336,671,2016,2688,3360,6048,7392,8568,8592,8616,8640,8664,8736,8904),
c(0,336,671,2016,2688,3360,6048,7392,8568)),

The output is Group 1 has N=20 only till WK54 and Group 2 is N=180 till WK-52
I wish to see one group which has overall impact of design (N=180 patients only till WK52 and after WK52 I have only 20 subjects for rich PK)

I am not sure if I am interpreting it correctly please advice

Really appreciate you spreading the word for PopED. Very useful tool.
Thanks,
-Ankur

[andrewhooker]

Hi

The design I implemented has 180 individuals with measurements between 0 and 52 weeks, and 20 individuals with those measurements PLUS additional measurements after 52 weeks. I assumed that the 20 individuals would be a part of the 200 total that are studies under 52 weeks and then 20 of those individuals are studied for an additional period.

[Anks2030]

Hi Andy,
Thank you very much it is clear now. Very useful and I really appreciate for your time and effort.
Do you have the analytical solution for the 2-compartment model.
I will try to do the optimization for this design and seems it is taking forever for ODE's
Thanks,
-Ankur