class: center, middle, inverse, title-slide # A gentle introduction of PopED functions ### <a href="https://www.linkedin.com/in/han-liu-a8824b14a/">Han</a> ### 2020/07/16 --- # Optimal design background PK/PD studies should be designed in such a way that model parameters can be estimated with adequate precision and bias. This can be assessed by simulation, but depending on the study and model(s) involved, it can be impractical to evaluate many combinations of design variables. Optimal design tools allow us to quickly (however approximately) evaluate designs, and even search over a design space for the best possible design. --- # Content - Not going to talk about the theory - [Stephen Duffull - Optimal Design of PKPD Studies](https://link.springer.com/chapter/10.1007%2F978-1-4419-7937-7_8) - Andy's seminar PPT - [Joakim Nyberg - PhD thesis](http://uu.diva-portal.org/smash/get/diva2:451173/FULLTEXT01.pdf) -- - Multiple examples can be find at - [Metrum group webinar](https://www.metrumrg.com/course/a-gentle-introduction-to-optimaldesign-for-pharmacometric-models/) - [Code for 15 examples with different feature by Andy](https://andrewhooker.github.io/PopED/articles/examples.html#analytic-solution-of-pkpd-model-multiple-study-arms-1) -- ### Get started ! ```r library(PopED) ``` --- # First thing to do is … - Define the **model** and **design**  -- - Functions - *ff ()*, the structure and covariate model; - *fg ()*, the IIV and IOV model; - *feps ()*, the residual error model; - *create.poped.database ()*, collect all models, define parameter values, define study design and save them as a database object **poped.db** --- ## Structure and covariate model ff() - Analytical solution - one compartment PK model - allometric scaling with a weight effect on clearance ```r ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ CL = CL*(WT/70)^(WT_CL) y = DOSE/V*exp(-CL/V*xt) return(list(y=y,poped.db=poped.db)) }) } ``` --- ## Structure and covariate model ff() - Analytical solution - one compartment PK model + direct Emax PD model - allometric scaling with a weight effect on clearance ```r ff <- function(model_switch,xt,parameters,poped.db){ with(as.list(parameters),{ # PK model CL=CL*(WT/70)^(WT_CL) CONC = DOSE/V*exp(-CL/V*xt) # PD model EFF = E0 + CONC*EMAX/(EC50 + CONC) y[MS==1] = CONC[model_switch==1] y[MS==2] = EFF[model_switch==2] return(list(y=y,poped.db=poped.db)) }) } ``` --- ## Structure and covariate model ff() - ODEs-based structure model - PopED + deSolve/ PKPDsim/ RxODE/ **mrgsolve** - [Tutorial by Andy](https://andrewhooker.github.io/PopED/articles/model_def_other_pkgs.html) ```r library(mrgsolve) code <- ' $PARAM CL=3.75, V=72.8, KA=0.25, BIO=0.9 $CMT DEPOT CENT $ODE dxdt_DEPOT = -KA*DEPOT; dxdt_CENT = KA*DEPOT - (CL/V)*CENT; $TABLE double CP = CENT/V*BIO; $CAPTURE CP ' moda <- mcode("optim", code, atol=1e-8, rtol=1e-8,maxsteps=5000) ``` --- ```r ff <- function(model_switch, xt, p, poped.db){ times_xt <- drop(xt) dose_times <- seq(from=0,to=max(times_xt),by=p[["TAU"]]) time <- sort(unique(c(0,times_xt,dose_times))) is.dose <- time %in% dose_times data <-tibble::tibble(ID = 1, time = time, amt = ifelse(is.dose,p[["DOSE"]], 0), cmt = ifelse(is.dose, 1, 0), evid = cmt, CL = p[["CL"]], V = p[["V"]], KA = p[["KA"]], BIO = p[["BIO"]]) out <- mrgsim_q(moda, data=data) y <- out$CP y <- y[match(times_xt,out$time)] return(list(y=matrix(y,ncol=1),poped.db=poped.db)) } ``` --- ## Stochastic model fg() - Define individual parameters with - **bpop**: typical value - **b**: IIV etas - **bocc**: IOV eta - **a**: covariates ```r fg <- function(x,a,bpop,b,bocc){ parameters=c( KA=bpop[1]*exp(b[1]), CL=bpop[2]*exp(b[2]), #CL=bpop[2]*exp(b[2]+bocc[1,1]), V=bpop[3]*exp(b[3]), BIO=bpop[4], DOSE=a[1], TAU=a[2]) return( parameters ) } ``` --- ## Residual model feps() ```r feps <- function(model_switch,xt,parameters,epsi,poped.db){ y <- do.call(poped.db$model$ff_pointer,list(model_switch,xt,parameters,poped.db))[[1]] y = y*(1+epsi[,1])+epsi[,2] return(list(y=y,poped.db=poped.db)) } ``` --- ## Collect all models and combined with design ```r poped.db <- create.poped.database( # Model ff_fun =ff, fg_fun =fg, fError_fun=feps, bpop=c(KA=0.25,CL=3.75,V=72.8,BIO=0.9), notfixed_bpop=c(1,1,1,0), d=c(V=0.09,KA=0.09,CL=0.15^2), # covd=c(0.48,0,0,0,0.45,0.53), covariance of IIV sigma=c(prop=0.04,add=0.0025), ds_index=c(0,0,0,0,0,0,0,0), ofv_calc_type=6, # Design m=2, groupsize=10, xt=c(1,8,240,245), a=cbind(DOSE=c(20,40),TAU=c(24,24)), # Design space discrete_xt = list(seq(0,250,1)), maxa=c(DOSE=200,TAU=24), mina=c(DOSE=0,TAU=24)) ``` --- # After Poped.db is ready … - Model prediction - Data simulation - Design evaluation - expected parameter uncertainty - **SSE** is recommended to validate the parameter uncertainty and calculate the bias - Design optimization - decrease expected parameter uncertainty - other stuff by define your own **criteria function** --- # Check design and parameter ```r poped.db$design$xt ``` ``` ## obs_1 obs_2 obs_3 obs_4 ## grp_1 1 8 240 245 ## grp_2 1 8 240 245 ``` ```r poped.db$parameter$param.pt.val$bpop ``` ``` ## [,1] ## KA 0.25 ## CL 3.75 ## V 72.80 ## BIO 0.90 ``` ```r model_prediction(poped.db, IPRED = T, DV = T) ``` ``` ## ID Time DV IPRED PRED Group Model DOSE TAU ## 1 1 1 0.091317306 0.05782968 0.05325024 1 1 20 24 ## 2 1 8 0.099595422 0.15927616 0.16409609 1 1 20 24 ## 3 1 240 0.116416047 0.10817642 0.12671376 1 1 20 24 ## 4 1 245 0.192054145 0.23579831 0.24980320 1 1 20 24 ## 5 2 1 0.039848348 0.04092116 0.05325024 1 1 20 24 ## 6 2 8 0.126890828 0.14536814 0.16409609 1 1 20 24 ## 7 2 240 0.062260946 0.11569698 0.12671376 1 1 20 24 ## 8 2 245 0.149253829 0.21580174 0.24980320 1 1 20 24 ## 9 3 1 0.120789758 0.06460762 0.05325024 1 1 20 24 ## 10 3 8 0.170583209 0.17885527 0.16409609 1 1 20 24 ## 11 3 240 0.093557498 0.07826940 0.12671376 1 1 20 24 ## 12 3 245 0.197761601 0.22682453 0.24980320 1 1 20 24 ## 13 4 1 0.003209404 0.06945290 0.05325024 1 1 20 24 ## 14 4 8 0.144675805 0.16242492 0.16409609 1 1 20 24 ## 15 4 240 0.097424146 0.15854318 0.12671376 1 1 20 24 ## 16 4 245 0.376104914 0.29341577 0.24980320 1 1 20 24 ## 17 5 1 -0.012132702 0.03252013 0.05325024 1 1 20 24 ## 18 5 8 0.121912986 0.12126602 0.16409609 1 1 20 24 ## 19 5 240 0.116388410 0.10170285 0.12671376 1 1 20 24 ## 20 5 245 0.164464113 0.18225046 0.24980320 1 1 20 24 ## 21 6 1 0.180486040 0.12065562 0.05325024 1 1 20 24 ## 22 6 8 0.252908496 0.20423375 0.16409609 1 1 20 24 ## 23 6 240 0.031083031 0.07911781 0.12671376 1 1 20 24 ## 24 6 245 0.227314515 0.28591386 0.24980320 1 1 20 24 ## 25 7 1 -0.025799750 0.04331601 0.05325024 1 1 20 24 ## 26 7 8 0.170419335 0.15128205 0.16409609 1 1 20 24 ## 27 7 240 0.079994340 0.16284803 0.12671376 1 1 20 24 ## 28 7 245 0.243312971 0.26607710 0.24980320 1 1 20 24 ## 29 8 1 0.013715264 0.05256555 0.05325024 1 1 20 24 ## 30 8 8 0.122920249 0.18612825 0.16409609 1 1 20 24 ## 31 8 240 0.173786224 0.17434819 0.12671376 1 1 20 24 ## 32 8 245 0.295036667 0.30138862 0.24980320 1 1 20 24 ## 33 9 1 0.033364356 0.02626024 0.05325024 1 1 20 24 ## 34 9 8 0.148485112 0.10108519 0.16409609 1 1 20 24 ## 35 9 240 0.151613797 0.16309355 0.12671376 1 1 20 24 ## 36 9 245 0.205484108 0.22642091 0.24980320 1 1 20 24 ## 37 10 1 0.054589540 0.03320253 0.05325024 1 1 20 24 ## 38 10 8 0.075765897 0.11300603 0.16409609 1 1 20 24 ## 39 10 240 0.221520557 0.09665134 0.12671376 1 1 20 24 ## 40 10 245 0.184237693 0.17618611 0.24980320 1 1 20 24 ## 41 11 1 0.092531052 0.10002743 0.10650048 2 1 40 24 ## 42 11 8 0.261185093 0.30920887 0.32819217 2 1 40 24 ## 43 11 240 0.185474240 0.17596131 0.25342752 2 1 40 24 ## 44 11 245 0.487548376 0.41320486 0.49960640 2 1 40 24 ## 45 12 1 0.088802940 0.08506654 0.10650048 2 1 40 24 ## 46 12 8 0.449303857 0.32349856 0.32819217 2 1 40 24 ## 47 12 240 0.156803376 0.26610620 0.25342752 2 1 40 24 ## 48 12 245 0.467291574 0.47828686 0.49960640 2 1 40 24 ## 49 13 1 0.040397668 0.13493558 0.10650048 2 1 40 24 ## 50 13 8 0.390010430 0.37447291 0.32819217 2 1 40 24 ## 51 13 240 0.252918201 0.26029164 0.25342752 2 1 40 24 ## 52 13 245 0.388962192 0.55870151 0.49960640 2 1 40 24 ## 53 14 1 0.087033812 0.03710848 0.10650048 2 1 40 24 ## 54 14 8 0.121029601 0.17266226 0.32819217 2 1 40 24 ## 55 14 240 0.213253244 0.24038342 0.25342752 2 1 40 24 ## 56 14 245 0.452376981 0.33353210 0.49960640 2 1 40 24 ## 57 15 1 0.097694920 0.14869083 0.10650048 2 1 40 24 ## 58 15 8 0.422879078 0.38429975 0.32819217 2 1 40 24 ## 59 15 240 0.148421270 0.21304006 0.25342752 2 1 40 24 ## 60 15 245 0.646284833 0.53563269 0.49960640 2 1 40 24 ## 61 16 1 0.182627573 0.16165331 0.10650048 2 1 40 24 ## 62 16 8 0.212255309 0.38155445 0.32819217 2 1 40 24 ## 63 16 240 0.167329489 0.23204445 0.25342752 2 1 40 24 ## 64 16 245 0.586922148 0.56167997 0.49960640 2 1 40 24 ## 65 17 1 0.104111487 0.08051749 0.10650048 2 1 40 24 ## 66 17 8 0.233330269 0.31644615 0.32819217 2 1 40 24 ## 67 17 240 0.068145605 0.29496294 0.25342752 2 1 40 24 ## 68 17 245 0.664197899 0.49599496 0.49960640 2 1 40 24 ## 69 18 1 0.212913669 0.13220228 0.10650048 2 1 40 24 ## 70 18 8 0.310307004 0.33284181 0.32819217 2 1 40 24 ## 71 18 240 0.317866331 0.26001573 0.25342752 2 1 40 24 ## 72 18 245 0.468182652 0.53396676 0.49960640 2 1 40 24 ## 73 19 1 0.141179438 0.18516444 0.10650048 2 1 40 24 ## 74 19 8 0.525965718 0.48640964 0.32819217 2 1 40 24 ## 75 19 240 0.149454425 0.19930271 0.25342752 2 1 40 24 ## 76 19 245 0.729465288 0.61594117 0.49960640 2 1 40 24 ## 77 20 1 0.143585388 0.16023821 0.10650048 2 1 40 24 ## 78 20 8 0.233084871 0.40034961 0.32819217 2 1 40 24 ## 79 20 240 0.156554823 0.13676966 0.25342752 2 1 40 24 ## 80 20 245 0.493405407 0.49341875 0.49960640 2 1 40 24 ``` --- # Plot the model ```r plot_model_prediction(poped.db, model_num_points = 300) ``` <!-- --> --- # Plot the model ```r plot_model_prediction(poped.db, PI=TRUE, separate.groups=T, model_num_points = 300, sample.times = FALSE) ``` <!-- --> ```r # By default 95% prediction interval of varibility # FO approximation to and normality assumption of the variance ``` --- # Evaluate the design ```r evaluation <- evaluate_design(poped.db) round(evaluation$rse) # in percent ``` ``` ## KA CL V d_KA d_CL d_V sig_prop sig_add ## 30 5 17 273 86 95 59 39 ``` --- # Optimize the design - optimize for parameter precision ```r #output <- poped_optim(poped.db, opt_xt=TRUE) ``` - optimize for parameter precision of non-nested models ```r #output <- poped_optim_MA(list(poped.db1, poped.db2, opt_xt=TRUE) ``` --- # Optimize the dose - optimize the dose for target concentration - typical population trough concentrations of 0.2 and 0.35 for the two groups of patients at 240 hours. ```r crit_fcn <- function(poped.db,...){ pred_df <- model_prediction(poped.db) sum((pred_df[pred_df["Time"]==240,"PRED"] - c(0.2,0.35))^2) } output_cost <- poped_optim(poped.db, opt_a = TRUE, ofv_fun=crit_fcn, maximize = FALSE) ``` --- ```r summary(output_cost) ``` ``` ## =============================================================================== ## FINAL RESULTS ## ## Optimized Covariates: ## Group 1: 31.5672 : 24 ## Group 2: 55.2426 : 24 ## ## OFV = 6.56342e-15 ## ## Efficiency: ## (ofv_final / ofv_init) = 4.4658e-13 ## ## Expected relative standard error ## (%RSE, rounded to nearest integer): ## Parameter Values RSE_0 RSE ## KA 0.25 30 24 ## CL 3.75 5 5 ## V 72.8 17 15 ## d_KA 0.09 273 209 ## d_CL 0.0225 86 70 ## d_V 0.09 95 78 ## sig_prop 0.04 59 48 ## sig_add 0.0025 39 54 ## ## Total running time: 14.5 seconds ``` --- # Optimize the dose - optimize the dose to maximize the percentage of IPRED falls into a therapeutic window ```r crit_fcn <- function(poped.db,...){ poped.db$design$groupsize[2,1] <- 1000 pred_df <- model_prediction(poped.db, IPRED = T, groups_to_use=2) Cmin <- pred_df[pred_df["Time"]==240,"IPRED"] length(Cmin[Cmin<0.4 & 0.25<Cmin])/length(Cmin) } # output_cost <- poped_optim(poped.db, opt_a = T, ofv_fun=crit_fcn, maximize = T) # summary(output_cost) ```