Comparison to popED
Introduction
This page shows a comparison of the results of mlxDesignEval and popED. To do so, we re-use the examples presented in the popED documentation.
popED calculates the RSE on the variance of the random effects and error model, while mlxDesignEval works with standard deviations by default. To request RSEs on the variances with mlxDesignEval, we will use the option rse_on_variance=T. For each example, the RSEs are displayed side by side and the percentage difference is calculated. If the difference is less than 1%, we consider that the results are the same (i.e is_identical=T). Note that the OFV and FIM will differ between mlxDesignEval and popED and cannot be compared directly.
The RSE are identical between popED and mlxDesignEval on all examples.
In order to allow an easy comparison of the models syntax, we will use the inlineModel() function to create model files. The results would be the same using a Monolix or Simulx project (containing the same model).
Ex1a: PK,1-comp,oral,MD
This example uses the analytical solution of a 1-compartment model with first-order absorption. The bioavailability is set as fixed as it cannot be estimated, and the error model parameters are also fixed. There are two groups with different dose amounts.
popED
CODE
##-- Model: One comp first order absorption
## -- Analytic solution for both mutiple and single dosing
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
N = floor(xt/TAU)+1
y=(DOSE*Favail/V)*(KA/(KA - CL/V)) *
(exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
return(list( y=y,poped.db=poped.db))
})
}
## -- parameter definition function
## -- names match parameters in function ff
sfg <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU=a[2])
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 ))
}
## -- Define design and design space
poped.db <- create.poped.database(ff_fun="ff",
fg_fun="sfg",
fError_fun="feps",
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
notfixed_bpop=c(1,1,1,0),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(0.04,5e-6),
notfixed_sigma=c(0,0),
m=2,
groupsize=20,
xt=c( 1,2,8,240,245),
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
bUseGrouped_xt=1,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24))
res1a_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
With mlxDesignEval:
DOSE and TAU do not need to be specified in the model
the model is defined using the
pkmodel()macro, which will be replaced by the corresponding analytical solution in the backgroundthe bioavailability
Fhas no random effects, it can be set as eithernormalorlogNormal(same results)y = y*(1+epsi[,1])+epsi[,2]in popED corresponds to acombined2error model in mlxtran language
CODE
model1a <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=normal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, Cl}
PK:
Cc = pkmodel(ka, V, Cl, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")When defining the population_parameters, remember that omega parameters represent a standard deviation while the popED parameter are variances. We thus use the square root of the value given in popED.
CODE
treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20))
treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
res1a_mlx <- mlxDesignEval::evaluate_design(
model_file=model1a,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","a","b"),
group=list(g1,g2),
output=list(output="DV", data=data.frame(time=c(1,2,8,240,245))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 8.215338 V_pop 8.215339 1.661207e-05 TRUE
#> 2 KA 10.090955 ka_pop 10.090937 1.835064e-04 TRUE
#> 3 CL 4.400304 Cl_pop 4.400307 6.276719e-05 TRUE
#> 4 d_V 39.833230 var_omega_V 39.833208 5.578078e-05 TRUE
#> 5 d_KA 60.089601 var_omega_ka 60.089600 1.920909e-06 TRUE
#> 6 d_CL 27.391518 var_omega_Cl 27.391561 1.557890e-04 TRUEEx1b: PK,1-comp,oral,MD
Same example as before but parametrized with elimination rate instead of clearance.
popED
CODE
## -- names match parameters in function defined in ff_file
fg.PK.1.comp.oral.md.param.2 <- function(x,a,bpop,b,bocc){
## -- parameter definition function
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
KE=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU=a[2])
return( parameters )
}
## -- Define design and design space
poped.db <- create.poped.database(ff_file="ff.PK.1.comp.oral.md.KE",
fg_file="fg.PK.1.comp.oral.md.param.2",
fError_file="feps.add.prop",
groupsize=20,
m=2,
sigma=c(0.04,5e-6),
bpop=c(V=72.8,KA=0.25,KE=3.75/72.8,Favail=0.9),
d=c(V=0.09,KA=0.09,KE=0.25^2),
notfixed_bpop=c(1,1,1,0),
notfixed_sigma=c(0,0),
xt=c( 1,2,8,240,245),
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
bUseGrouped_xt=1,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=40),
mina=c(DOSE=0,TAU=2))
## evaluate initial design
res1b_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
CODE
model1b <- inlineModel("
[INDIVIDUAL]
input = {k_pop, omega_k, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
k = {distribution=logNormal, typical=k_pop, sd=omega_k}
F = {distribution=normal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, k}
PK:
Cc = pkmodel(ka, V, k, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
CODE
treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20))
treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
res1b_mlx <- mlxDesignEval::evaluate_design(
model_file = model1b,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, k_pop=3.75/72.8,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_k=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","a","b"),
group = list(g1,g2),
output = list(output="DV", data=data.frame(time=c( 1,2,8,240,245))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 8.215338 V_pop 8.215339 1.980025e-05 TRUE
#> 2 KA 10.090955 ka_pop 10.090937 1.838245e-04 TRUE
#> 3 KE 7.566975 k_pop 7.566974 4.300219e-06 TRUE
#> 4 d_V 31.220520 var_omega_V 31.220506 4.337963e-05 TRUE
#> 5 d_KA 44.677836 var_omega_ka 44.677859 5.190823e-05 TRUE
#> 6 d_KE 38.005067 var_omega_k 38.005052 3.844144e-05 TRUEEx1c: PK,1-comp,oral,MD
Same example as 1a but written with an ODE system instead of an analytical solution.
popED
The ODE system can be solved in R with deSolve (slow) or in C++ using Rcpp (faster).
CODE
## define the ODE
PK.1.comp.oral.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
dA1 <- -KA*A1
dA2 <- KA*A1 - (CL/V)*A2
return(list(c(dA1, dA2)))
})
}
## define the initial conditions and the dosing
ff.ode <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0)
times_xt <- drop(xt) #xt[,,drop=T]
dose_times = seq(from=0,to=max(times_xt),by=TAU)
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE*Favail), method = c("add"))
times <- sort(c(times_xt,dose_times))
out <- ode(A_ini, times, PK.1.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = out[, "A2"]/(V)
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
sfg <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU=a[2])
return( parameters )
}
poped.db <- create.poped.database(ff_fun=ff.ode,
fError_fun=feps.add.prop,
fg_fun=sfg,
groupsize=20,
m=2, #number of groups
sigma=c(0.04,5e-6),
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9),
d=c(V=0.09,KA=0.09,CL=0.25^2),
notfixed_bpop=c(1,1,1,0),
notfixed_sigma=c(0,0),
xt=c( 1,2,8,240,245),
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
discrete_xt = list(0:248),
bUseGrouped_xt=1,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24))
# calculations are noticeably slower than with the analytic solution
res1c_popED_Rode <- PopED::evaluate_design(poped.db)
# using Rcpp
cppFunction('List one_comp_oral_ode(double Time, NumericVector A, NumericVector Pars) {
int n = A.size();
NumericVector dA(n);
double CL = Pars[0];
double V = Pars[1];
double KA = Pars[2];
dA[0] = -KA*A[0];
dA[1] = KA*A[0] - (CL/V)*A[1];
return List::create(dA);
}')
ff.ode.rcpp <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0)
times_xt <- drop(xt) #xt[,,drop=T]
dose_times = seq(from=0,to=max(times_xt),by=TAU)
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE*Favail), method = c("add"))
times <- sort(c(times_xt,dose_times))
out <- ode(A_ini, times, one_comp_oral_ode, c(CL,V,KA), events = list(data = eventdat))
y = out[, "A2"]/(V)
y=y[match(times_xt,out[,"time"])]
y=cbind(y)
return(list(y=y,poped.db=poped.db))
})
}
## -- Update poped.db with compiled function
poped.db.compiled.rcpp <- create.poped.database(poped.db,ff_fun=ff.ode.rcpp)
# calculations are much faster than with the pure R solution but slower than .dll compiled solution in desolve
res1c_popED_Rcpp <- PopED::evaluate_design(poped.db.compiled.rcpp)mlxDesignEval
The model in mlxtran language is always and automatically converted to C++ code to run fast. There is nothing to do on the user side. The solver for stiff ODEs has a better precision and is usually preferred. It can be set with odeType=stiff in the structural model.
CODE
model1c <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=logNormal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, Cl}
PK:
depot(target=Ad)
EQUATION:
odeType=stiff
ddt_Ad = -ka*Ad
ddt_Ac = ka*Ad - (Cl/V)*Ac
Cc = Ac/V
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
CODE
treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20))
treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
res1c_mlx <- mlxDesignEval::evaluate_design(
model_file = model1c,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","a","b"),
group = list(g1,g2),
output = list(output="DV", data=data.frame(time=c( 1,2,8,240,245))),
rse_on_variance = T)Comparison
Comparison to solution with deSolve:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 8.215334 V_pop 8.212722 0.031791907 TRUE
#> 2 KA 10.090963 ka_pop 10.084245 0.066577753 TRUE
#> 3 CL 4.400304 Cl_pop 4.400000 0.006913958 TRUE
#> 4 d_V 39.833193 var_omega_V 39.821267 0.029938234 TRUE
#> 5 d_KA 60.089706 var_omega_ka 60.030178 0.099065990 TRUE
#> 6 d_CL 27.391516 var_omega_Cl 27.387705 0.013911217 TRUEComparison to solution with Rcpp:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 8.215334 V_pop 8.212722 0.031791907 TRUE
#> 2 KA 10.090963 ka_pop 10.084245 0.066577753 TRUE
#> 3 CL 4.400304 Cl_pop 4.400000 0.006913958 TRUE
#> 4 d_V 39.833193 var_omega_V 39.821267 0.029938234 TRUE
#> 5 d_KA 60.089706 var_omega_ka 60.030178 0.099065990 TRUE
#> 6 d_CL 27.391516 var_omega_Cl 27.387705 0.013911217 TRUEEx2: warfarin
This example uses the Warfarin model used in the design evaluation software comparison Nyberg et al., “Methods and software tools for design evaluation for population pharmacokinetics-pharmacodynamics studies”, Br. J. Clin. Pharm., 2014.. It uses the analytical solution of a 1-compartment model with first-order absorption with a single dose.
popED
CODE
ff <- function(model_switch,xt,parameters,poped.db){
##-- Model: One comp first order absorption
with(as.list(parameters),{
y=xt
y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
return(list(y=y,poped.db=poped.db))
})
}
sfg <- function(x,a,bpop,b,bocc){
## -- parameter definition function
parameters=c(CL=bpop[1]*exp(b[1]),
V=bpop[2]*exp(b[2]),
KA=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1])
return(parameters)
}
feps <- function(model_switch,xt,parameters,epsi,poped.db){
## -- Residual Error function
## -- Proportional
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])
return(list(y=y,poped.db=poped.db))
}
## -- Define initial design and design space
poped.db <- create.poped.database(ff_file="ff",
fg_file="sfg",
fError_file="feps",
bpop=c(CL=0.15, V=8, KA=1.0, Favail=1),
notfixed_bpop=c(1,1,1,0),
d=c(CL=0.07, V=0.02, KA=0.6),
sigma=0.01,
groupsize=32,
xt=c( 0.5,1,2,6,24,36,72,120),
minxt=0,
maxxt=120,
a=70)
## evaluate initial design
res2a_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
The model is the same as for the example 1, except that a proportional error model is used.
CODE
model2 <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=logNormal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {b}
input = {F, ka, V, Cl}
PK:
Cc = pkmodel(ka, V, Cl, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=proportional(b)}")
CODE
res2a_mlx <- mlxDesignEval::evaluate_design(
model_file = model2,
population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15,
omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07),
b=sqrt(0.01)),
fixed_parameters = c("F_pop"),
group = list(size=32),
treatment = list(data=data.frame(time=0, amount=70)),
output =list(output="DV", data=data.frame(time=c(0.5,1,2,6,24,36,72,120))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 4.738266 Cl_pop 4.738270 6.837835e-05 TRUE
#> 2 V 2.756206 V_pop 2.756208 4.766023e-05 TRUE
#> 3 KA 13.925829 ka_pop 13.925839 7.575014e-05 TRUE
#> 4 d_CL 25.627205 var_omega_Cl 25.627221 6.124399e-05 TRUE
#> 5 d_V 30.344316 var_omega_V 30.344315 3.414725e-06 TRUE
#> 6 d_KA 25.777327 var_omega_ka 25.777328 3.518511e-06 TRUE
#> 7 SIGMA[1,1] 11.170784 var_b 11.170783 9.912263e-06 TRUEEx3: PKPD,1-comp,oral,MD
This example uses a model with two outputs (PK and PD), both using an analytical solution. The PK model has 1 compartment and first order absorption. The PD is an Imax model.
popED
CODE
##-- Model: One comp first order absorption + inhibitory imax
## -- works for both mutiple and single dosing
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
MS <- model_switch
# PK model
N = floor(xt/TAU)+1
CONC=(DOSE*Favail/V)*(KA/(KA - CL/V)) *
(exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
# PD model
EFF = E0*(1 - CONC*IMAX/(IC50 + CONC))
y[MS==1] = CONC[MS==1]
y[MS==2] = EFF[MS==2]
return(list( y= y,poped.db=poped.db))
})
}
## -- parameter definition function
sfg <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1],
TAU = a[2],
E0=bpop[5]*exp(b[4]),
IMAX=bpop[6],
IC50=bpop[7])
return( parameters )
}
## -- Residual Error function
feps <- function(model_switch,xt,parameters,epsi,poped.db){
returnArgs <- ff(model_switch,xt,parameters,poped.db)
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
MS <- model_switch
pk.dv <- y*(1+epsi[,1])+epsi[,2]
pd.dv <- y*(1+epsi[,3])+epsi[,4]
y[MS==1] = pk.dv[MS==1]
y[MS==2] = pd.dv[MS==2]
return(list( y= y,poped.db =poped.db ))
}
poped.db <- create.poped.database(ff_fun="ff",
fError_fun="feps",
fg_fun="sfg",
groupsize=20,
m=3,
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,E0=1120,IMAX=0.807,IC50=0.0993),
notfixed_bpop=c(1,1,1,0,1,1,1),
d=c(V=0.09,KA=0.09,CL=0.25^2,E0=0.09),
sigma=c(0.04,5e-6,0.09,100),
notfixed_sigma=c(0,0,0,0),
xt=c( 1,2,8,240,240,1,2,8,240,240),
minxt=c(0,0,0,240,240,0,0,0,240,240),
maxxt=c(10,10,10,248,248,10,10,10,248,248),
discrete_xt = list(0:248),
G_xt=c(1,2,3,4,5,1,2,3,4,5),
bUseGrouped_xt=1,
model_switch=c(1,1,1,1,1,2,2,2,2,2),
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24),c(DOSE=0, TAU=24)),
maxa=c(DOSE=200,TAU=40),
mina=c(DOSE=0,TAU=2),
ourzero=0)
## evaluate initial design
res3a_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
CODE
model3 <- inlineModel("
[INDIVIDUAL]
input = {V_pop, omega_V, KA_pop, omega_KA, CL_pop, omega_CL, Favail_pop, E0_pop, omega_E0, IMAX_pop, IC50_pop}
DEFINITION:
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
KA = {distribution=logNormal, typical=KA_pop, sd=omega_KA}
CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL}
Favail = {distribution=logNormal, typical=Favail_pop, no-variability}
E0 = {distribution=logNormal, typical=E0_pop, sd=omega_E0}
IMAX = {distribution=logNormal, typical=IMAX_pop, no-variability}
IC50 = {distribution=logNormal, typical=IC50_pop, no-variability}
[LONGITUDINAL]
input = {aCONC, bCONC, aEFF, bEFF}
input = {V,KA,CL,Favail,E0,IMAX,IC50}
EQUATION:
Cc = pkmodel(ka=KA,V,Cl=CL,p=Favail)
Eff = E0*(1 - Cc*IMAX/(IC50 + Cc))
OUTPUT:
output = {Cc,Eff}
DEFINITION:
yCONC = {distribution=normal, prediction=Cc, errorModel=combined2(aCONC, bCONC)}
yEFF = {distribution=normal, prediction=Eff, errorModel=combined2(aEFF, bEFF)}")The popED example defines two PK measurements at the same time t=240. If the same is done with mlxDesignEval, the two times are merged into a single one, and the information that two measurements are done at the same time is lost. To capture the double measurements, we give two slightly different times: 240 and 240.001.
CODE
treatment1 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20))
treatment2 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40))
treatment3 <- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=0))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
g3 <- list(size=20, treatment=treatment3)
# it is not possible to give twice the same time (will be merged into one by simulx) so giving slightly different value
outPK <- list(output="yCONC", data=data.frame(time=c(1,2,8,240,240.001)))
outPD <- list(output="yEFF", data=data.frame(time=c(1,2,8,240,240.001)))
res3a_mlx_option1 <- mlxDesignEval::evaluate_design(
model_file = model3,
population_parameters = data.frame(Favail_pop=0.9, KA_pop=0.25, V_pop=72.8, CL_pop=3.75,
E0_pop=1120,IMAX_pop=0.807,IC50_pop=0.0993,
omega_KA=sqrt(0.09), omega_V=sqrt(0.09),
omega_CL=sqrt(0.25^2),omega_E0=sqrt(0.09),
bCONC=sqrt(0.04),aCONC=sqrt(5e-6),bEFF=sqrt(0.09),aEFF=sqrt(100)),
fixed_parameters = c("Favail_pop","bCONC","aCONC","bEFF","aEFF"),
group = list(g1,g2,g3),
output = list(outPK,outPD),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 8.119842 V_pop 8.119667 2.153595e-03 TRUE
#> 2 KA 9.968612 KA_pop 9.968355 2.580165e-03 TRUE
#> 3 CL 4.304635 CL_pop 4.304765 3.011526e-03 TRUE
#> 4 E0 7.076883 E0_pop 7.076857 3.718358e-04 TRUE
#> 5 IMAX 9.895340 IMAX_pop 9.894347 1.003892e-02 TRUE
#> 6 IC50 39.478269 IC50_pop 39.477335 2.365973e-03 TRUE
#> 7 d_V 38.960998 var_omega_V 38.960436 1.443299e-03 TRUE
#> 8 d_KA 58.523188 var_omega_KA 58.521487 2.905635e-03 TRUE
#> 9 d_CL 25.832775 var_omega_CL 25.833807 3.996265e-03 TRUE
#> 10 d_E0 22.036110 var_omega_E0 22.036132 9.915972e-05 TRUEEx4: PKPD,1-comp,Emax
This example uses a model with two outputs (PK and PD), both using an analytical solution. The PK model has 1 compartment and bolus administration. The PD is an Emax model.
popED
CODE
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
MS <- model_switch
# PK model
CONC = DOSE/V*exp(-CL/V*xt)
# PD model
EFF = E0 + CONC*EMAX/(EC50 + CONC)
y[MS==1] = CONC[MS==1]
y[MS==2] = EFF[MS==2]
return(list( y= y,poped.db=poped.db))
})
}
## -- parameter definition function
sfg <- function(x,a,bpop,b,bocc){
parameters=c(
CL=bpop[1]*exp(b[1]) ,
V=bpop[2]*exp(b[2]) ,
E0=bpop[3]*exp(b[3]) ,
EMAX=bpop[4],
EC50=bpop[5]*exp(b[4]) ,
DOSE=a[1]
)
return( parameters )
}
## -- Residual Error function
## -- Proportional PK + additive PD
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]]
MS <- model_switch
prop.err <- y*(1+epsi[,1])
add.err <- y+epsi[,2]
y[MS==1] = prop.err[MS==1]
y[MS==2] = add.err[MS==2]
return(list( y= y,poped.db =poped.db ))
}
poped.db <- create.poped.database(ff_fun=ff,
fError_fun=feps,
fg_fun=sfg,
groupsize=20,
m=3,
sigma=diag(c(0.15,0.015)),
bpop=c(CL=0.5,V=0.2,E0=1,EMAX=1,EC50=1),
d=c(CL=0.09,V=0.09,E0=0.04,EC50=0.09),
xt=c( 0.33,0.66,0.9,5,0.1,1,2,5),
bUseGrouped_xt=1,
model_switch=c( 1,1,1,1,2,2,2,2),
minxt=0,
maxxt=5,
ourzero = 0,
a=list(c(DOSE=0),c(DOSE=1),c(DOSE=2)),
maxa=c(DOSE=10),
mina=c(DOSE=0))
## evaluate initial design
res4_popED <- PopED::evaluate_design(poped.db)
#> Problems inverting the matrix. Results could be misleading.mlxDesignEval
CODE
model4 <- inlineModel("
[INDIVIDUAL]
input = {V_pop, omega_V, CL_pop, omega_CL, E0_pop, omega_E0, EC50_pop, omega_EC50, EMAX_pop}
DEFINITION:
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL}
E0 = {distribution=logNormal, typical=E0_pop, sd=omega_E0}
EC50 = {distribution=logNormal, typical=EC50_pop, sd=omega_EC50}
EMAX = {distribution=logNormal, typical=EMAX_pop, no-variability}
[LONGITUDINAL]
input = {bCONC, aEFF}
input = {V,CL,E0,EMAX,EC50}
EQUATION:
Cc = pkmodel(V=V,Cl=CL)
Eff = E0 + Cc*EMAX/(EC50 + Cc)
OUTPUT:
output = {Cc,Eff}
DEFINITION:
yCONC = {distribution=normal, prediction=Cc, errorModel=proportional(bCONC)}
yEFF = {distribution=normal, prediction=Eff, errorModel=constant(aEFF)}")
CODE
treatment1<- list(data=data.frame(time=0, amount=0))
treatment2<- list(data=data.frame(time=0, amount=1))
treatment3<- list(data=data.frame(time=0, amount=2))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
g3 <- list(size=20, treatment=treatment3)
outPK <- list(output="yCONC", data=data.frame(time=c(0.33,0.66,0.9,5)))
outPD <- list(output="yEFF", data=data.frame(time=c(0.1,1,2,5)))
res4_mlx <- mlxDesignEval::evaluate_design(
model_file = model4,
population_parameters = data.frame(CL_pop=0.5,V_pop=0.2,E0_pop=1,EMAX_pop=1,EC50_pop=1,
omega_V=sqrt(0.09), omega_CL=sqrt(0.09),omega_E0=sqrt(0.04),omega_EC50=sqrt(0.09),
bCONC=sqrt(0.15),aEFF=sqrt(0.015)),
group = list(g1,g2,g3),
output = list(outPK,outPD),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 6.032896 CL_pop 6.032893 4.680235e-05 TRUE
#> 2 V 6.288412 V_pop 6.288410 2.768071e-05 TRUE
#> 3 E0 2.802303 E0_pop 2.802303 1.030050e-05 TRUE
#> 4 EMAX 3.555580 EMAX_pop 3.555587 1.938459e-04 TRUE
#> 5 EC50 14.508492 EC50_pop 14.508480 7.788633e-05 TRUE
#> 6 d_CL 30.322296 var_omega_CL 30.322269 8.930289e-05 TRUE
#> 7 d_V 32.964398 var_omega_V 32.964358 1.224145e-04 TRUE
#> 8 d_E0 20.298320 var_omega_E0 20.298325 2.401653e-05 TRUE
#> 9 d_EC50 142.920289 var_omega_EC50 142.920335 3.266862e-05 TRUE
#> 10 SIGMA[1,1] 14.728239 var_bCONC 14.728240 1.091985e-05 TRUE
#> 11 SIGMA[2,2] 11.723581 var_aEFF 11.723583 2.183795e-05 TRUEEx5: PD,Emax-Hill
This example uses a model with PD only, with an analytical solution of an Emax model with sigmoidicity depending on the dose.
popED
Note that the time xt is used as a surrogate for the dose (DOSE=xt).
CODE
ff.emax.hill <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
DOSE = xt
y=BASE + EMAX*DOSE^(GAMMA)/(ED50^(GAMMA) + DOSE^(GAMMA))
return(list( y= y,poped.db=poped.db))
})
}
## -- parameter definition function
sfg.emax.hill <- function(x,a,bpop,b,bocc){
parameters=c( EMAX=bpop[1]*exp(b[1]),
ED50=bpop[2]*exp(b[2]),
GAMMA=bpop[3],
BASE=bpop[4]+b[3])
return( parameters )
}
poped.db <- create.poped.database(ff_fun=ff.emax.hill,
fError_fun=feps.add.prop,
fg_fun=sfg.emax.hill,
groupsize=100,
m=1,
bpop=c(EMAX=100,ED50=20,GAMMA=4.5,BASE=1),
d=c(EMAX=0.0625,ED50=0.0625,BASE=0.0625),
sigma=diag(c(0.01,.1)),
xt=seq(0,50,length.out=8),
minxt=0,
maxxt=50,
ourzero=0)
## evaluate initial design
res5_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
Emax and EC50 have a lognormal distribution while Base has a normal distribution. Gamma has no variability and can be set as either logNormal or normal. In order to consider the dose on the x axis instead of time, we rename the internal time variable t as DOSE (in the same way as done with popED).
CODE
model5 <- inlineModel("
[INDIVIDUAL]
input = {EMAX_pop, BASE_pop, omega_BASE, ED50_pop, omega_ED50, GAMMA_pop, omega_EMAX}
DEFINITION:
EMAX = {distribution=logNormal, typical=EMAX_pop, sd=omega_EMAX}
BASE = {distribution=normal, typical=BASE_pop, sd=omega_BASE}
ED50 = {distribution=logNormal, typical=ED50_pop, sd=omega_ED50}
GAMMA = {distribution=normal, typical=GAMMA_pop, no-variability}
[LONGITUDINAL]
input = {aCONC, bCONC}
input = {BASE,EMAX,ED50,GAMMA}
EQUATION:
DOSE=t
Pred=BASE + EMAX*DOSE^(GAMMA)/(ED50^(GAMMA) + DOSE^(GAMMA))
OUTPUT:
output = {Pred}
DEFINITION:
yDV = {distribution=normal, prediction=Pred, errorModel=combined2(aCONC, bCONC)}")The different doses tested are defined in the output element, which data frame column should stay time (even if it now represents the dose).
CODE
out <- list(output="yDV", data=data.frame(time=seq(0,50,length.out=8)))
res5_mlx <- mlxDesignEval::evaluate_design(
model_file = model5,
population_parameters = data.frame(EMAX_pop=100,ED50_pop=20,GAMMA_pop=4.5,BASE_pop=1,
omega_EMAX=sqrt(0.0625),omega_ED50=sqrt(0.0625),omega_BASE=sqrt(0.0625),
bCONC=sqrt(0.01),aCONC=sqrt(0.1)),
group = list(size=100),
output = out,
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 EMAX 2.588804 EMAX_pop 2.588805 1.742681e-05 TRUE
#> 2 ED50 2.554668 ED50_pop 2.554668 9.715213e-06 TRUE
#> 3 GAMMA 1.112444 GAMMA_pop 1.112444 1.946455e-05 TRUE
#> 4 BASE 3.922208 BASE_pop 3.922206 5.208619e-05 TRUE
#> 5 d_EMAX 14.825482 var_omega_EMAX 14.825492 6.249867e-05 TRUE
#> 6 d_ED50 14.375072 var_omega_ED50 14.375067 3.504967e-05 TRUE
#> 7 d_BASE 33.257104 var_omega_BASE 33.257108 1.172602e-05 TRUE
#> 8 SIGMA[1,1] 7.072694 var_bCONC 7.072680 1.946146e-04 TRUE
#> 9 SIGMA[2,2] 18.487445 var_aCONC 18.487444 1.021925e-06 TRUEEx6: PK,1-comp,oral,SD
This example uses a simple PK model with 1 compartment and first-order absorption, receiving a single dose. The ODE system solution and analytical solution are compared.
popED
CODE
##-- Model: One comp first order absorption, analytic solution
PK.1.comp.oral.sd.ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
y=(DOSE*Favail*KA/(V*(KA-KE)))*(exp(-KE*xt)-exp(-KA*xt))
return(list( y= y,poped.db=poped.db))
})
}
# ODE solution
PK.1.comp.oral.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
dA1 <- -KA*A1
dA2 <- KA*A1 - KE*A2
return(list(c(dA1, dA2)))
})
}
PK.1.comp.oral.sd.ff.ode <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
A_ini <- c(A1 = DOSE, A2 = 0)
times <- drop(xt)
times <- sort(times)
times <- c(0,times) # add extra time for start of integration
out <- ode(A_ini, times, PK.1.comp.oral.ode, parameters) #,atol=1e-13,rtol=1e-13)
y = out[,"A2"]/(V/Favail)
y=y[-1] # remove initial time for start of integration
y = cbind(y) ## must be a row vector
return(list( y= y,poped.db=poped.db))
})
}
## -- parameter definition function
PK.1.comp.oral.sd.fg.param.1 <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1])
parameters["KE"]=parameters["CL"]/parameters["V"]
return( parameters )
}
PK.1.comp.oral.sd.fg.param.2 <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
KE=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1])
return( parameters )
}
poped.db.1 <- create.poped.database(ff_fun="PK.1.comp.oral.sd.ff",
fg_fun="PK.1.comp.oral.sd.fg.param.1",
fError_fun="feps.add.prop",
groupsize=32,
m=1,
sigma=diag(c(0.01,0.25)),
bpop=c(V=8,KA=1,CL=0.15,Favail=1),
d=c(V=0.02,KA=0.6,CL=0.07),
notfixed_bpop=c(1,1,1,0),
xt=c(1,2,3,6,24,36,72,120),
minxt=0,
maxxt=c(25,25,25,120,120,120,120,120),
discrete_xt = list(1:120),
a=cbind(c(70)),
bUseGrouped_xt=1,
maxa=c(200),
mina=c(0))
poped.db.2 <- create.poped.database(poped.db.1,
fg_fun="PK.1.comp.oral.sd.fg.param.2",
bpop=c(V=8,KA=1,KE=0.15/8,Favail=1),
d=c(V=0.02,KA=0.6,KE=0.07))
poped.db.3 <- create.poped.database(poped.db.1,
ff_fun="PK.1.comp.oral.sd.ff.ode")
## evaluate initial designs
res6a_popED <- PopED::evaluate_design(poped.db.1) # analytical solution with Cl
res6b_popED <- PopED::evaluate_design(poped.db.2) # analytical solution with KE
res6c_popED <- PopED::evaluate_design(poped.db.3) # ODEsmlxDesignEval
CODE
# analytical solution with pkmodel() parametrized with Cl
model6a <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=normal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, Cl}
PK:
Cc = pkmodel(ka, V, Cl, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
# analytical solution with pkmodel() parametrized with k
model6b <- inlineModel("
[INDIVIDUAL]
input = {k_pop, omega_k, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
k = {distribution=logNormal, typical=k_pop, sd=omega_k}
F = {distribution=normal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, k}
PK:
Cc = pkmodel(ka, V, k, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
# ODE model parametrized with Cl
model6c <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=normal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, Cl}
PK:
depot(target=Ad)
EQUATION:
odeType=stiff
ddt_Ad = -ka*Ad
ddt_Ac = ka*Ad - (Cl/V)*Ac
Cc = Ac/V
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")The same design is evaluated for the three models.
CODE
# analytical solution with Cl
res6a_mlx <- mlxDesignEval::evaluate_design(
model_file = model6a,
population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15,
omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07),
a=sqrt(0.25),b=sqrt(0.01)),
fixed_parameters = c("F_pop"),
group = list(size=32),
treatment = list(data=data.frame(time=0,amount=70)),
output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))),
rse_on_variance = T)
# analytical solution with k
res6b_mlx <- mlxDesignEval::evaluate_design(
model_file = model6b,
population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, k_pop=0.15/8,
omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_k=sqrt(0.07),
a=sqrt(0.25),b=sqrt(0.01)),
fixed_parameters = c("F_pop"),
group = list(size=32),
treatment = list(data=data.frame(time=0,amount=70)),
output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))),
rse_on_variance = T)
# ODE system
res6c_mlx <- mlxDesignEval::evaluate_design(
model_file = model6c,
population_parameters = data.frame(F_pop=1, ka_pop=1, V_pop=8, Cl_pop=0.15,
omega_ka=sqrt(0.6), omega_V=sqrt(0.02), omega_Cl=sqrt(0.07),
a=sqrt(0.25),b=sqrt(0.01)),
fixed_parameters = c("F_pop"),
group = list(size=32),
treatment = list(data=data.frame(time=0,amount=70)),
output = list(output="DV", data=data.frame(time=c(1,2,3,6,24,36,72,120))),
rse_on_variance = T)Comparison
Analytical solution with Cl:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 2.946875 V_pop 2.946876 4.858601e-05 TRUE
#> 2 KA 14.572054 ka_pop 14.572062 5.500982e-05 TRUE
#> 3 CL 5.096568 Cl_pop 5.096567 1.770528e-05 TRUE
#> 4 d_V 34.369463 var_omega_V 34.369453 3.158243e-05 TRUE
#> 5 d_KA 27.846026 var_omega_ka 27.846004 7.837298e-05 TRUE
#> 6 d_CL 29.782073 var_omega_Cl 29.782069 1.165045e-05 TRUE
#> 7 SIGMA[1,1] 26.863004 var_b 26.862986 6.770186e-05 TRUE
#> 8 SIGMA[2,2] 25.642019 var_a 25.642028 3.589427e-05 TRUEAnalytical solution with k:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 2.946875 V_pop 2.946876 5.191786e-05 TRUE
#> 2 KA 14.572054 ka_pop 14.572062 5.524263e-05 TRUE
#> 3 KE 5.437768 k_pop 5.437766 3.526612e-05 TRUE
#> 4 d_V 33.258184 var_omega_V 33.258192 2.252980e-05 TRUE
#> 5 d_KA 27.727936 var_omega_ka 27.727963 9.708115e-05 TRUE
#> 6 d_KE 32.722277 var_omega_k 32.722299 6.570502e-05 TRUE
#> 7 SIGMA[1,1] 26.509409 var_b 26.509395 5.598068e-05 TRUE
#> 8 SIGMA[2,2] 25.306924 var_a 25.306916 3.139526e-05 TRUEODE system with Cl:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 2.947023 V_pop 2.946876 4.989732e-03 TRUE
#> 2 KA 14.572062 ka_pop 14.572062 5.305491e-06 TRUE
#> 3 CL 5.097250 Cl_pop 5.096567 1.339638e-02 TRUE
#> 4 d_V 34.368092 var_omega_V 34.369453 3.957485e-03 TRUE
#> 5 d_KA 27.845897 var_omega_ka 27.846004 3.848681e-04 TRUE
#> 6 d_CL 29.780671 var_omega_Cl 29.782069 4.694171e-03 TRUE
#> 7 SIGMA[1,1] 26.864337 var_b 26.862986 5.030488e-03 TRUE
#> 8 SIGMA[2,2] 25.645051 var_a 25.641997 1.190959e-02 TRUEEx7: PK,1-comp,maturation
This example uses a 1-comp PK model with weight and post-menstrual age (PMA) as a covariate on Cl and V. The weight is not an input, instead it is calculated from the PMA value with different formulas for males and females. Four groups with different PMA values are defined in the design.
popED
The function to calculate the weight from the post-menstrual age is defined as a separate function which is then called in the model.
CODE
# typical WT prediction from Sumpter & Holford 2011
# PMA in years
get_typical_weight <- function(PMA,SEX=2){
WTmax1 <- 2.76
TM50wt1 <- 38.5
HILLwt1 <- 12.9
HILL2wt1 <- 2.74
PMA1 <- PMA*52
WTmax2 <- 16.4
TM50wt2 <- 2.1
HILLwt2 <- 2.04
HILL2wt2 <- 1
WTmax3 <- 40.2
TM50wt3 <- 12.4
HILLwt3 <- 2.87
HILL2wt3 <- 0
TLAGwt4 <- 12.4
WTmax4 <- 33.6
THALFwt4 <- 3.61
FFEM <- 1
if(SEX==1) FFEM <- 0.884
wt1 <- FFEM*WTmax1/(1+(TM50wt1/PMA1)^((PMA1<TM50wt1)*HILLwt1+(PMA1>=TM50wt1)*HILL2wt1))
wt2 <- WTmax2/(1+(TM50wt2/PMA)^((PMA<TM50wt2)*HILLwt2+(PMA>=TM50wt2)*HILL2wt2))
wt3 <- WTmax3/(1+(TM50wt3/PMA)^((PMA<TM50wt3)*HILLwt3+(PMA>=TM50wt3)*HILL2wt3))
wt4 <- (PMA > TLAGwt4)*(FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4))))
return(wt1 + wt2 + wt3 + wt4)
}
PK.1.comp.maturation.ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
WT <- get_typical_weight(PMA,SEX=SEX)
CL=CL*(WT/70)^(3/4)*(PMA^HILL)/(TM50^HILL+PMA^HILL)
V=V*(WT/70)
DOSE=1000*(WT/70)
y = DOSE/V*exp(-CL/V*xt)
return(list( y= y,poped.db=poped.db))
})
}
PK.1.comp.maturation.fg <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]),
V=bpop[2]*exp(b[2]),
TM50=bpop[3]*exp(b[3]),
HILL=bpop[4],
PMA=a[1],
SEX=a[2])
return( parameters )
}
poped.db <- create.poped.database(ff_file="PK.1.comp.maturation.ff",
fError_file="feps.add.prop",
fg_file="PK.1.comp.maturation.fg",
groupsize=rbind(50,20,20,20),
m=4,
sigma=c(0.015,0.0015),
notfixed_sigma = c(1,0),
bpop=c(CL=3.8,V=20,TM50=60/52,HILL=3),
d=c(CL=0.05,V=0.05,TM50=0.05),
xt=c( 1,2,4,6,8,24),
minxt=0,
maxxt=24,
bUseGrouped_xt=1,
a=list(c(PMA=25,SEX=2),
c(PMA=15,SEX=2),
c(PMA=10,SEX=2),
c(PMA=5,SEX=2)),
maxa=list(c(PMA=30,SEX=2)),
mina=list(c(PMA=1,SEX=2)))
## evaluate initial design
res7_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
PMA and SEX are passed as regressors such that they can be handled in a flexible way in the structural model. The calculation of the typical weight for each PMA is done directly in the model.
In this example, the dose is defined directly in the model and the explicit analytical formula is used. Note that it would also be possible to use the pkmodel() macro and define the treatment as part of the evaluate_design input. In that case, the scaling by weight could be done via the p= argument of pkmodel.
CODE
model7 <- inlineModel("
[INDIVIDUAL]
input = {CL_pop, omega_CL, V_pop, omega_V, TM50_pop, omega_TM50, HILL_pop}
DEFINITION:
CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
TM50 = {distribution=logNormal, typical=TM50_pop, sd=omega_TM50}
HILL = {distribution=logNormal, typical=HILL_pop, no-variability}
[LONGITUDINAL]
input = {a, b}
input = {CL,V,TM50,HILL,PMA,SEX}
PMA = {use=regressor}
SEX = {use=regressor}
EQUATION:
WTmax1 = 2.76
TM50wt1 = 38.5
HILLwt1 = 12.9
HILL2wt1 = 2.74
PMA1 = PMA*52
WTmax2 = 16.4
TM50wt2 = 2.1
HILLwt2 = 2.04
HILL2wt2 = 1
WTmax3 = 40.2
TM50wt3 = 12.4
HILLwt3 = 2.87
HILL2wt3 = 0
TLAGwt4 = 12.4
WTmax4 = 33.6
THALFwt4 = 3.61
if(SEX==1)
FFEM = 0.884
else
FFEM = 1
end
if (PMA1<TM50wt1)
if (PMA1>=TM50wt1)
wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILLwt1+HILL2wt1))
else
wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILLwt1))
end
else
if (PMA1>=TM50wt1)
wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(HILL2wt1))
else
wt1 = FFEM*WTmax1/(1+(TM50wt1/PMA1)^(0))
end
end
if (PMA<TM50wt2)
if (PMA>=TM50wt2)
wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILLwt2+HILL2wt2))
else
wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILLwt2))
end
else
if (PMA>=TM50wt2)
wt2 = WTmax2/(1+(TM50wt2/PMA)^(HILL2wt2))
else
wt2 = WTmax2/(1+(TM50wt2/PMA)^(0))
end
end
if (PMA<TM50wt3)
if (PMA>=TM50wt3)
wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILLwt3+HILL2wt3))
else
wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILLwt3))
end
else
if (PMA>=TM50wt2)
wt3 = WTmax3/(1+(TM50wt3/PMA)^(HILL2wt3))
else
wt3 = WTmax3/(1+(TM50wt3/PMA)^(0))
end
end
if (PMA > TLAGwt4)
wt4 = (FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4))))
else
wt4 = 0*(FFEM*WTmax4*(1-exp(-log(2)/THALFwt4*(PMA-TLAGwt4))))
end
WT = wt1 + wt2 + wt3 + wt4
CL1 = CL*(WT/70)^(3/4)*(PMA^HILL)/(TM50^HILL+PMA^HILL)
V1 = V *(WT/70)
DOSE=1000*(WT/70)
Cc = DOSE/V1*exp(-CL1/V1*t)
OUTPUT:
output = {Cc}
DEFINITION:
y1 = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")Each group has its own regressor element defining the PMA and SEX. There is no treatment definition in evaluate_design as the dose is defined as part of the structural model.
CODE
reg1 <- data.frame(time=0,PMA=25, SEX=2)
reg2 <- data.frame(time=0,PMA=15, SEX=2)
reg3 <- data.frame(time=0,PMA=10, SEX=2)
reg4 <- data.frame(time=0,PMA=5, SEX=2)
g1 <-list(size=50, regressor=reg1)
g2 <-list(size=20, regressor=reg2)
g3 <-list(size=20, regressor=reg3)
g4 <-list(size=20, regressor=reg4)
res7_mlx <- mlxDesignEval::evaluate_design(
model_file = model7,
population_parameters = data.frame(CL_pop=3.8,V_pop=20,TM50_pop=60/52,HILL_pop=3,
omega_CL=sqrt(0.05), omega_V=sqrt(0.05), omega_TM50=sqrt(0.05),
a=sqrt(0.0015),b=sqrt(0.015)),
fixed_parameters = c("a"),
group = list(g1,g2,g3,g4),
output = list(output="y1", data=data.frame(time=c(1,2,4,6,8,24))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 3.724808 CL_pop 3.724809 1.139376e-05 TRUE
#> 2 V 2.251799 V_pop 2.251799 9.947249e-06 TRUE
#> 3 TM50 3048.507973 TM50_pop 3048.507832 4.620080e-06 TRUE
#> 4 HILL 2117.356641 HILL_pop 2117.356843 9.528602e-06 TRUE
#> 5 d_CL 15.751133 var_omega_CL 15.751139 3.949789e-05 TRUE
#> 6 d_V 14.997327 var_omega_V 14.997333 3.759456e-05 TRUE
#> 7 d_TM50 27968.366291 var_omega_TM50 27968.367565 4.556077e-06 TRUE
#> 8 SIGMA[1,1] 6.950869 var_b 6.950875 7.528790e-05 TRUEEx8: TMDD
This example uses a TMDD model with quasi steady-state (QSS) approximation, encoded as an ODE system.
popED
The system is solved using deSolve in R (results with C++ using Rcpp are identical). The results are quite sensitive to the ODE solver tolerance. We fave set the tolerance in the popED script to be the same as the default Monolix tolerances. Starting from version 2024, the tolerances can also be changed in Monolix.
CODE
tmdd_qss_one_target_model_ode <- function(Time,State,Pars){
with(as.list(c(State, Pars)), {
RTOT = A4
CTOT= A2/V1
CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT))
dA1 = -KA*A1
dA2 = FAVAIL*KA*A1+(Q/V2)*A3-(CL/V1+Q/V1)*CFREE*V1-RTOT*KINT*CFREE*V1/(KSSS+CFREE)
dA3 = (Q/V1)*CFREE*V1 - (Q/V2)*A3
dA4 = R0*KDEG - KDEG*RTOT - (KINT-KDEG)*(RTOT*CFREE/(KSSS+CFREE))
return(list(c(dA1,dA2,dA3,dA4)))
})
}
sfg <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]) ,
V1=bpop[2]*exp(b[2]) ,
Q=bpop[3]*exp(b[3]) ,
V2=bpop[4]*exp(b[4]) ,
FAVAIL=bpop[5]*exp(b[5]) ,
KA=bpop[6]*exp(b[6]) ,
VMAX=bpop[7]*exp(b[7]) ,
KMSS=bpop[8]*exp(b[8]) ,
R0=bpop[9]*exp(b[9]) ,
KSSS=bpop[10]*exp(b[10]) ,
KDEG=bpop[11]*exp(b[11]) ,
KINT=bpop[12]*exp(b[12]) ,
DOSE=a[1] ,
SC_FLAG=a[2])
return(parameters)
}
tmdd_qss_one_target_model <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
#The initialization vector for the compartment
A_ini <- c(A1=DOSE*SC_FLAG,
A2=DOSE*(1-SC_FLAG),
A3=0,
A4=R0)
#Set up time points for the ODE
times_xt <- drop(xt)
times <- sort(times_xt)
times <- c(0,times) ## add extra time for start of integration
# solve the ODE
out <- ode(A_ini, times, tmdd_qss_one_target_model_ode, parameters, atol=1e-9,rtol=1e-6) # as monolix default
# extract the time points of the observations
out = out[match(times_xt,out[,"time"]),]
# Match ODE output to measurements
RTOT = out[,"A4"]
CTOT = out[,"A2"]/V1
CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT))
COMPLEX=((RTOT*CFREE)/(KSSS+CFREE))
RFREE= RTOT-COMPLEX
y[model_switch==1]= RTOT[model_switch==1]
y[model_switch==2]= CFREE[model_switch==2]
#y[model_switch==3]=RFREE[model_switch==3]
return(list( y=y,poped.db=poped.db))
})
}
tmdd_qss_one_target_model_ruv <- 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[model_switch==1] = log(y[model_switch==1])+epsi[,1]
y[model_switch==2] = log(y[model_switch==2])+epsi[,2]
#y[model_switch==3] = log(y[model_switch==3])+epsi[,3]
return(list(y=y,poped.db=poped.db))
}
#################################################
# for study 1 in gibiansky,JPKPD,2012 table 2
#################################################
# for study 1 in gibiansky,JPKPD,2012 table 2
poped.db.1 <- create.poped.database(ff_fun=tmdd_qss_one_target_model,
fError_fun=tmdd_qss_one_target_model_ruv,
fg_fun=sfg,
groupsize=6,
m=4, #number of groups
sigma=c(0.04,0.0225),
bpop=c(CL=0.3,V1=3,Q=0.2,V2=3,FAVAIL=0.7,KA=0.5,VMAX=0,
KMSS=0,R0=0.1,KSSS=0.015,KDEG=10,KINT=0.05),
d=c(CL=0.09,V1=0.09,Q=0.04,V2=0.04,FAVAIL=0.04,KA=0.16,VMAX=0,
KMSS=0,R0=0.09,KSSS=0.09,KDEG=0.04,KINT=0.04),
notfixed_bpop=c( 1,1,1,1,1,1,0,0,1,1,1,1),
notfixed_d=c( 1,1,1,1,1,1,0,0,1,1,1,1),
xt=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56,
0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56),
model_switch=c(1,1,1,1,1,1,1,1,1,1,1,1,1,
2,2,2,2,2,2,2,2,2,2,2,2,2),
bUseGrouped_xt=1,
G_xt=c(1,2,3,4,5,6,7,8,9,10,11,12,13,
1,2,3,4,5,6,7,8,9,10,11,12,13),
a=list(c(DOSE=100, SC_FLAG=0),
c(DOSE=300, SC_FLAG=0),
c(DOSE=600, SC_FLAG=0),
c(DOSE=1000, SC_FLAG=1)),
discrete_a = list(DOSE=seq(100,1000,by=100),
SC_FLAG=c(0,1)))
res8a_popED <- PopED::evaluate_design(poped.db.1)mlxDesignEval
The model in mlxtran language is always and automatically converted to C++ code to run fast. There is nothing to do on the user side. The solver for stiff ODEs has a better precision and is usually preferred. It can be set with odeType=stiff in the structural model. IV ans sub-cutaneous (SC) administrations are distinguished by the adm keyword.
CODE
model8 <- inlineModel("
[INDIVIDUAL]
input = {Q_pop, omega_Q, R0_pop, omega_R0, V2_pop, omega_V2, CL_pop, omega_CL, FAVAIL_pop, omega_FAVAIL, KA_pop, omega_KA, KDEG_pop, omega_KDEG, KINT_pop, omega_KINT, KSSS_pop, omega_KSSS, V1_pop, omega_V1}
DEFINITION:
Q = {distribution=logNormal, typical=Q_pop, sd=omega_Q}
R0 = {distribution=logNormal, typical=R0_pop, sd=omega_R0}
V2 = {distribution=logNormal, typical=V2_pop, sd=omega_V2}
CL = {distribution=logNormal, typical=CL_pop, sd=omega_CL}
FAVAIL = {distribution=logNormal, typical=FAVAIL_pop, sd=omega_FAVAIL}
KA = {distribution=logNormal, typical=KA_pop, sd=omega_KA}
KDEG = {distribution=logNormal, typical=KDEG_pop, sd=omega_KDEG}
KINT = {distribution=logNormal, typical=KINT_pop, sd=omega_KINT}
KSSS = {distribution=logNormal, typical=KSSS_pop, sd=omega_KSSS}
V1 = {distribution=logNormal, typical=V1_pop, sd=omega_V1}
[LONGITUDINAL]
input = {aCfree, aRtot}
input = {KA, FAVAIL, V1, KINT, KSSS, KDEG, R0, CL, Q, V2}
PK:
; IV
depot(target = A2, adm = 1)
; SC
depot(target = A1, adm = 2)
EQUATION:
odeType = stiff
A1_0 = 0
A2_0 = 0
A3_0 = 0
A4_0 = R0
; Parameter transformations
RTOT = A4
CTOT = A2/V1
CFREE = 0.5*((CTOT-RTOT-KSSS)+sqrt((CTOT-RTOT-KSSS)^2+4*KSSS*CTOT))
ddt_A1 = -KA*A1
ddt_A2 = FAVAIL*KA*A1+(Q/V2)*A3-(CL/V1+Q/V1)*CFREE*V1-RTOT*KINT*CFREE*V1/(KSSS+CFREE)
ddt_A3 = (Q/V1)*CFREE*V1 - (Q/V2)*A3
ddt_A4 = R0*KDEG - KDEG*RTOT - (KINT-KDEG)*(RTOT*CFREE/(KSSS+CFREE))
OUTPUT:
output = {CFREE, RTOT}
DEFINITION:
yCfree = {distribution=logNormal, prediction=CFREE, errorModel=constant(aCfree)}
yRtot = {distribution=logNormal, prediction=RTOT, errorModel=constant(aRtot)}")We define 3 IV groups with ascending doses and one SC group.
CODE
trt1 <- list(data=data.frame(time=0, amount=100), admID=1)
trt2 <- list(data=data.frame(time=0, amount=300), admID=1)
trt3 <- list(data=data.frame(time=0, amount=600), admID=1)
trt4 <- list(data=data.frame(time=0, amount=1000), admID=2)
g1 <- list(size=6, treatment=trt1)
g2 <- list(size=6, treatment=trt2)
g3 <- list(size=6, treatment=trt3)
g4 <- list(size=6, treatment=trt4)
outCfree <- list(output="yCfree", data=data.frame(time=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56)))
outRtot <- list(output="yRtot", data=data.frame(time=c(0.0417,0.25,0.5,1,3,7,14,21,28,35,42,49,56)))
res8_mlx <- mlxDesignEval::evaluate_design(
model_file = model8,
population_parameters = data.frame(KA_pop=0.5, V1_pop=3, CL_pop=0.3, FAVAIL_pop=0.7,
KINT_pop=0.05, KSSS_pop=0.015, Q_pop=0.2, V2_pop=3,
R0_pop=0.1, KDEG_pop=10,
omega_KA=sqrt(0.16), omega_V1=sqrt(0.09), omega_CL=sqrt(0.09),
omega_FAVAIL=sqrt(0.04),omega_KINT=sqrt(0.04), omega_KSSS=sqrt(0.09),
omega_Q=sqrt(0.04), omega_V2=sqrt(0.04),
omega_R0=sqrt(0.09), omega_KDEG=sqrt(0.04),
aCfree=sqrt(0.0225), aRtot=sqrt(0.04)),
group = list(g1,g2,g3,g4),
output = list(outCfree,outRtot),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 7.125817 CL_pop 7.125686 1.844242e-03 TRUE
#> 2 V1 6.866975 V1_pop 6.867008 4.873102e-04 TRUE
#> 3 Q 7.409078 Q_pop 7.408239 1.132901e-02 TRUE
#> 4 V2 10.435157 V2_pop 10.436235 1.032650e-02 TRUE
#> 5 FAVAIL 11.471389 FAVAIL_pop 11.471443 4.658107e-04 TRUE
#> 6 KA 19.592448 KA_pop 19.592448 3.111509e-06 TRUE
#> 7 R0 8.408518 R0_pop 8.408258 3.089987e-03 TRUE
#> 8 KSSS 11.031286 KSSS_pop 11.028962 2.107175e-02 TRUE
#> 9 KDEG 8.176590 KDEG_pop 8.175696 1.094110e-02 TRUE
#> 10 KINT 7.344912 KINT_pop 7.343644 1.726866e-02 TRUE
#> 11 d_CL 32.695336 var_omega_CL 32.695158 5.452809e-04 TRUE
#> 12 d_V1 33.548081 var_omega_V1 33.548108 7.942589e-05 TRUE
#> 13 d_Q 74.871951 var_omega_Q 74.864811 9.535884e-03 TRUE
#> 14 d_V2 84.900937 var_omega_V2 84.901826 1.046666e-03 TRUE
#> 15 d_FAVAIL 98.434249 var_omega_FAVAIL 98.434293 4.446352e-05 TRUE
#> 16 d_KA 78.750482 var_omega_KA 78.750556 9.337863e-05 TRUE
#> 17 d_R0 36.358020 var_omega_R0 36.359746 4.746907e-03 TRUE
#> 18 d_KSSS 49.621583 var_omega_KSSS 49.607974 2.742569e-02 TRUE
#> 19 d_KDEG 66.473682 var_omega_KDEG 66.478869 7.803029e-03 TRUE
#> 20 d_KINT 49.673985 var_omega_KINT 49.675044 2.131655e-03 TRUE
#> 21 SIGMA[1,1] 8.935594 var_aRtot 8.935620 2.905946e-04 TRUE
#> 22 SIGMA[2,2] 9.676886 var_aCfree 9.676859 2.780330e-04 TRUEEx9: PK,2-comp,oral,ODE
This is a two compartment model with oral absorption using ODEs.
popED
The system is solved using deSolve in R.
CODE
PK.2.comp.oral.ode <- function(Time, State, Pars){
with(as.list(c(State, Pars)), {
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.oral.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 = seq(from=0,to=max(times_xt),by=TAU)
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE), method = c("add"))
times <- sort(c(times_xt,dose_times))
out <- ode(A_ini, times, PK.2.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = out[, "A2"]/(V1/Favail)
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],
KA=bpop[3]*exp(b[2]),
Q=bpop[4],
V2=bpop[5],
Favail=bpop[6],
DOSE=a[1],
TAU=a[2])
return( parameters )
}
#' create poped database
poped.db <- create.poped.database(ff_fun="ff.PK.2.comp.oral.md.ode",
fError_fun="feps.add.prop",
fg_fun="fg",
groupsize=20,
m=1, #number of groups
sigma=c(prop=0.1^2,add=0.05^2),
bpop=c(CL=10,V1=100,KA=1,Q= 3.0, V2= 40.0, Favail=1),
d=c(CL=0.15^2,KA=0.25^2),
notfixed_bpop=c(1,1,1,1,1,0),
xt=c( 48,50,55,65,70,85,90,120),
minxt=0,
maxxt=144,
discrete_xt = list(0:144),
a=c(DOSE=100,TAU=24),
maxa=c(DOSE=1000,TAU=24),
mina=c(DOSE=0,TAU=8),
discrete_a = list(DOSE=seq(0,1000,by=100),TAU=8:24))
res9_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
CODE
model9 <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, Q_pop, V1_pop, V2_pop, ka_pop, omega_ka}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
Q = {distribution=normal, typical=Q_pop, no-variability}
V1 = {distribution=normal, typical=V1_pop, no-variability}
V2 = {distribution=normal, typical=V2_pop, no-variability}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {ka, Cl, V1, Q, V2}
PK:
depot(target=Ad)
EQUATION:
odeType=stiff
; Parameter transformations
V = V1
k12 = Q/V1
k21 = Q/V2
k = Cl/V1
ddt_Ad = -ka * Ad
ddt_Ac = ka*Ad - k*Ac - k12*Ac + k21*Ap
ddt_Ap = k12*Ac - k21*Ap
Cc=Ac/V
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
CODE
res9_mlx <- mlxDesignEval::evaluate_design(
model_file = model9,
population_parameters = data.frame(ka_pop=1, V1_pop=100, Cl_pop=10, Q_pop=3, V2_pop =40,
omega_ka=0.25, omega_Cl=0.15,
a=0.05, b=0.1),
group = list(size=20),
treatment = list(data=data.frame(start=0, interval=24, nbDoses=6, amount=100)),
output = list(output="DV", data=data.frame(time=c(48,50,55,65,70,85,90,120))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 4.409665 Cl_pop 4.40966 1.204776e-04 TRUE
#> 2 V1 12.071097 V1_pop 12.07108 1.341462e-04 TRUE
#> 3 KA 32.825506 ka_pop 32.82545 1.773338e-04 TRUE
#> 4 Q 65.765415 Q_pop 65.76515 3.965393e-04 TRUE
#> 5 V2 77.649537 V2_pop 77.64977 3.058298e-04 TRUE
#> 6 d_CL 36.561733 var_omega_Cl 36.56174 2.154027e-05 TRUE
#> 7 d_KA 191.947572 var_omega_ka 191.94783 1.324427e-04 TRUE
#> 8 sig_prop 64.576067 var_b 64.57616 1.385793e-04 TRUE
#> 9 sig_add 20.574939 var_a 20.57494 1.310480e-05 TRUEEx10: PKPD,HCV
This example uses the hepatitis C virus (HCV) model from the publication Nyberg et al., “Methods and software tools for design evaluation for population pharmacokinetics-pharmacodynamics studies”, Br. J. Clin. Pharm., 2014..
popED
The popED example uses the compiled code. It is therefore not possible to see the equations.
CODE
sfg <- function(x,a,bpop,b,bocc){
## -- parameter definition function
parameters=c(p=bpop[1],
d=bpop[2],
e=bpop[3],
s=bpop[4],
KA=bpop[5] + b[1],
KE=bpop[6] + b[2],
VD=bpop[7] + b[3],
EC50=bpop[8] + b[4],
n=bpop[9] + b[5],
delta=bpop[10] + b[6],
c=bpop[11] + b[7],
DOSE=a[1],
TINF=a[2],
TAU=a[3])
return(parameters)
}
ff_ODE_compiled <- function(model_switch,xt,parameters,poped.db){
parameters[5:11] <- exp(parameters[5:11])
with(as.list(parameters),{
A_ini <- c(A1 = 0, A2 = 0, A3=c*delta/(p*e),
A4=(s*e*p-d*c*delta)/(p*delta*e),
A5=(s*e*p-d*c*delta)/(c*delta*e))
#Set up time points for the ODE
times_xt <- drop(xt)
times <- c(0,times_xt) ## add extra time for start of integration
times <- sort(times)
times <- unique(times) # remove duplicates
# compute values from ODEs
out <- ode(A_ini, times, func = "derivs", parms = parameters,
#method="daspk",
#jacfunc = "jac",
dllname = "HCV_ode",
initfunc = "initmod", #nout = 1, outnames = "Sum",
atol=1e-12,rtol=1e-12)
# grab timepoint values
out = out[match(times_xt,out[,"time"]),]
y <- xt*0
pk <- out[,"A2"]/VD
pd <- out[,"A5"]
y[model_switch==1] <- pk[model_switch==1]
y[model_switch==2] <- log10(pd[model_switch==2])
y = cbind(y) # must be a column matrix
return(list( y= y,poped.db=poped.db))
})
}
feps_ODE_compiled <- function(model_switch,xt,parameters,epsi,poped.db){
## -- Residual Error function
MS<-model_switch
y <- ff_ODE_compiled(model_switch,xt,parameters,poped.db)[[1]]
pk.dv <- y + epsi[,1]
pd.dv <- y + epsi[,2]
y[MS==1] = pk.dv[MS==1]
y[MS==2] = pd.dv[MS==2]
return(list(y=y,poped.db=poped.db))
}
## -- Define initial design and design space
poped_db_compiled <- create.poped.database(ff_file="ff_ODE_compiled",
fg_file="sfg",
fError_file="feps_ODE_compiled",
bpop=c(p=100,
d=0.001,
e=1e-7,
s=20000,
KA=log(0.8),
KE=log(0.15),
VD=log(100),#VD=log(100000),
EC50=log(0.12), #EC50=log(0.00012),
n=log(2),
delta=log(0.2),
c=log(7)),
notfixed_bpop=c(0,0,0,0,1,1,1,1,1,1,1),
d=c(KA=0.25,
KE=0.25,
VD=0.25,
EC50=0.25,
n=0.25,
delta=0.25,
c=0.25),
sigma=c(0.04,0.04),
groupsize=30,
xt=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28,
0,0.25,0.5,1,2,3,4,7,10,14,21,28),
model_switch=c(rep(1,12),rep(2,12)),
a=c(180,1,7))
res10_popED <- PopED::evaluate_design(poped_db_compiled)mlxDesignEval
CODE
model10 <- inlineModel("
[INDIVIDUAL]
input = {lEC50_pop, omega_lEC50, lVd_pop, omega_lVd, lc_pop, omega_lc, ldelta_pop, omega_ldelta, lka_pop, omega_lka, lke_pop, omega_lke, ln_pop, omega_ln}
DEFINITION:
lEC50 = {distribution=normal, typical=lEC50_pop, sd=omega_lEC50}
lVd = {distribution=normal, typical=lVd_pop, sd=omega_lVd}
lc = {distribution=normal, typical=lc_pop, sd=omega_lc}
ldelta = {distribution=normal, typical=ldelta_pop,sd=omega_ldelta}
lka = {distribution=normal, typical=lka_pop, sd=omega_lka}
lke = {distribution=normal, typical=lke_pop, sd=omega_lke}
ln = {distribution=normal, typical=ln_pop, sd=omega_ln}
[LONGITUDINAL]
input = {aC, aW}
input = {lka,lke,lVd,lEC50,ln,ldelta,lc}
PK:
depot(target = X, p = 1, adm = 1)
EQUATION:
odeType = stiff
ka=exp(lka)
ke=exp(lke)
Vd=exp(lVd)
EC50=exp(lEC50)
n=exp(ln)
delta=exp(ldelta)
c=exp(lc)
q=100
d=0.001
e=0.0000001
s=20000
; Initial conditions
t_0 = 0
X_0 = 0
T_0 = (c*delta)/(q*e)
A_0 = 0
I_0 = (s*e*q-d*c*delta)/(q*delta*e)
W_0 = (s*e*q-d*c*delta)/(c*delta*e)
; Ordinary Differential Equations
C=A/Vd
ddt_X = -ka*X
ddt_A = ka*X-ke*A
ddt_T = s-T*(e*W+d)
ddt_I = e*W*T-delta*I
ddt_W = q*(1-(max(C,0)^n)/(max(C,0)^n+EC50^n))*I-c*W
log10W = log10(W)
OUTPUT:
output = {C,log10W}
DEFINITION:
yC = {distribution=normal, prediction=C, errorModel=constant(aC)}
yW = {distribution=normal, prediction=log10W, errorModel=constant(aW)}")
CODE
outPK <- list(output="yC", data=data.frame(time=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28)))
outPD <- list(output="yW", data=data.frame(time=c(0,0.25,0.5,1,2,3,4,7,10,14,21,28)))
res10_mlx <- mlxDesignEval::evaluate_design(
model_file = model10,
population_parameters = data.frame(lka_pop=log(0.8), lke_pop=log(0.15), lVd_pop=log(100), lEC50_pop=log(0.12),
ln_pop=log(2), ldelta_pop=log(0.2), lc_pop=log(7),
omega_lka=sqrt(0.25), omega_lke=sqrt(0.25), omega_lVd=sqrt(0.25),
omega_lEC50=sqrt(0.25), omega_ln=sqrt(0.25),
omega_ldelta=sqrt(0.25), omega_lc=sqrt(0.25),
aC=sqrt(0.04), aW=sqrt(0.04)),
group = list(size=30),
output = list(outPK,outPD),
treatment = list(data=data.frame(start=0, interval=7, nbDoses=6, amount=180, tinf=1)),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 KA 54.103880 lka_pop 54.103938 1.079865e-04 TRUE
#> 2 KE 5.520538 lke_pop 5.520530 1.466914e-04 TRUE
#> 3 VD 2.163313 lVd_pop 2.163313 9.106420e-06 TRUE
#> 4 EC50 7.432080 lEC50_pop 7.432087 8.158977e-05 TRUE
#> 5 n 15.079243 ln_pop 15.079235 5.244492e-05 TRUE
#> 6 delta 5.833940 ldelta_pop 5.833941 1.394196e-05 TRUE
#> 7 c 5.656054 lc_pop 5.656058 5.658128e-05 TRUE
#> 8 d_KA 39.856357 var_omega_lka 39.856362 1.173282e-05 TRUE
#> 9 d_KE 30.738289 var_omega_lke 30.738302 4.195786e-05 TRUE
#> 10 d_VD 28.610608 var_omega_lVd 28.610599 3.133398e-05 TRUE
#> 11 d_EC50 60.518673 var_omega_lEC50 60.518691 3.001940e-05 TRUE
#> 12 d_n 28.979133 var_omega_ln 28.979137 1.595303e-05 TRUE
#> 13 d_delta 27.179856 var_omega_ldelta 27.179845 3.854860e-05 TRUE
#> 14 d_c 32.751496 var_omega_lc 32.751501 1.418497e-05 TRUE
#> 15 SIGMA[1,1] 8.373173 var_aC 8.373171 1.446238e-05 TRUE
#> 16 SIGMA[2,2] 9.164741 var_aW 9.164742 8.385020e-06 TRUEEx11: PK, prior FIM
This example shows the usage of a prior FIM.
popED
CODE
# This example shows how to include a prior FIM into the design evaluation.
# We look at PK assessment in pediatrics, where we are mainly interested
# in assessing if there is a substantial difference of more than 20% in
# clearance (CL) between children and adults.
# First we setup the general model, then the designs for adults and pediatrics,
# and finally we can evaluate the separate designs and the pooled data.
##-- Model: One comp first order absorption
## -- Analytic solution for both multiple and single dosing
ff <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
N = floor(xt/TAU)+1
y=(DOSE*Favail/V)*(KA/(KA - CL/V)) *
(exp(-CL/V * (xt - (N - 1) * TAU)) * (1 - exp(-N * CL/V * TAU))/(1 - exp(-CL/V * TAU)) -
exp(-KA * (xt - (N - 1) * TAU)) * (1 - exp(-N * KA * TAU))/(1 - exp(-KA * TAU)))
return(list( y=y,poped.db=poped.db))
})
}
## -- parameter definition function
## -- names match parameters in function ff
## -- note, covariate on clearance for pediatrics (using isPediatric x[1])
sfg <- function(x,a,bpop,b,bocc){
parameters=c( V=bpop[1]*exp(b[1]),
KA=bpop[2]*exp(b[2]),
CL=bpop[3]*exp(b[3])*bpop[5]^a[3], # add covariate for pediatrics
Favail=bpop[4],
isPediatric = a[3],
DOSE=a[1],
TAU=a[2])
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 ))
}
## -- Define design and design space for adults (isPediatric = 0)
## Two arms, 5 time points
poped.db <- create.poped.database(ff_fun="ff",
fg_fun="sfg",
fError_fun="feps",
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8),
notfixed_bpop=c(1,1,1,0,1),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(0.04,5e-6),
notfixed_sigma=c(0,0),
m=2,
groupsize=20,
xt=c( 1,8,10,240,245),
bUseGrouped_xt=1,
a=list(c(DOSE=20,TAU=24,isPediatric = 0),
c(DOSE=40, TAU=24,isPediatric = 0)))
# Note, to be able to use the adults FIM to combine with the pediatrics,
# both have to have the parameter "pedCL" defined and set notfixed_bpop to 1.
## Define pediatric model/design (isPediatric = 1)
## One arm, 4 time points only
poped.db.ped <- create.poped.database(ff_fun="ff",
fg_fun="sfg",
fError_fun="feps",
bpop=c(V=72.8,KA=0.25,CL=3.75,Favail=0.9,pedCL=0.8),
notfixed_bpop=c(1,1,1,0,1),
d=c(V=0.09,KA=0.09,CL=0.25^2),
sigma=c(0.04,5e-6),
notfixed_sigma=c(0,0),
m=1,
groupsize=6,
xt=c( 1,2,6,240),
bUseGrouped_xt=1,
a=list(c(DOSE=40,TAU=24,isPediatric = 1)))
# We can evaluate the adult design without warning, by setting the pedCL
# parameter to be fixed (i.e., not estimated)
res11_popED_adult <- PopED::evaluate_design(create.poped.database(poped.db, notfixed_bpop=c(1,1,1,0,0)))
# One obtains good estimates for all parameters for adults (<60% RSE for all).
## evaluate design of pediatrics only - insufficient
# Similarly as before with only pediatrics we cannot estimate the covariate effect,
# so we fix it.
res11_popED_ped <- PopED::evaluate_design(create.poped.database(poped.db.ped, notfixed_bpop=c(1,1,1,0,0)))
# Due to having less subjects, less samples per subject, and only one dose level
# the variability in pediatrics cannot be estimated well.
## Add adult prior
# Now we combined the two sutdies, where we assume that we only assess a difference
# between adults and pediatrics in CL.
# We can set the prior FIM to the adult one:
outAdult <- res11_popED_adult
poped.db.all <- create.poped.database(
poped.db.ped,
prior_fim = outAdult$fim,
notfixed_bpop=c(1,1,1,0,0)
)
## evaluate design using prior FIM from adults
res11_popED_both <- PopED::evaluate_design(poped.db.all)
# Obviously, the pooled data leads to much higher precision in parameter estimates
# compared to the pediatrics only.mlxDesignEval
The isPediatric indicator is passed as a regressor to adjust the clearance value between adults and children.
CODE
model11 <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka, pedCL_pop}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=logNormal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
pedCL = {distribution=logNormal, typical=pedCL_pop, no-variability}
[LONGITUDINAL]
input = {a, b}
input = {F, ka, V, Cl, pedCL, isPediatric}
isPediatric = {use=regressor}
PK:
Clnew = Cl * pedCL^isPediatric
Cc = pkmodel(ka, V, Cl=Clnew, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")
CODE
treatment1<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=20))
treatment2<- list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
res11_mlx_adult <- mlxDesignEval::evaluate_design(
model_file = model11,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","pedCL_pop","a","b"),
group = list(g1,g2),
regressor = data.frame(time=0,isPediatric=0),
output = list(output="DV", data=data.frame(time=c( 1,8,10,240,245))),
rse_on_variance = T)
res11_mlx_ped <- mlxDesignEval::evaluate_design(
model_file = model11,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","pedCL_pop","a","b"),
group = list(size=6),
treatment = list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)),
regressor = data.frame(time=0,isPediatric=1),
output = list(output="DV", data=data.frame(time=c( 1,2,6,240))),
rse_on_variance = T)
res11_mlx_both <- mlxDesignEval::evaluate_design(
model_file = model11,
population_parameters = data.frame(F_pop=0.9, ka_pop=0.25, V_pop=72.8, Cl_pop=3.75, pedCL_pop=0.8,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09), omega_Cl=0.25,
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("F_pop","pedCL_pop","a","b"),
group = list(size=6),
treatment = list(data=data.frame(start=0, interval=24, nbDoses=11, amount=40)),
regressor = data.frame(time=0,isPediatric=1),
output = list(output="DV", data=data.frame(time=c( 1,2,6,240))),
rse_on_variance = T,
prior_FIM = res11_mlx_adult$FIM)Comparison
Adult trial:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 6.634931 V_pop 6.634930 3.545338e-06 TRUE
#> 2 KA 8.587203 ka_pop 8.587206 3.117019e-05 TRUE
#> 3 CL 4.354792 Cl_pop 4.354790 2.426530e-05 TRUE
#> 4 d_V 33.243601 var_omega_V 33.243613 3.433354e-05 TRUE
#> 5 d_KA 55.689432 var_omega_ka 55.689417 2.726942e-05 TRUE
#> 6 d_CL 27.133255 var_omega_Cl 27.133212 1.586636e-04 TRUEPediatric trial:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 24.72088 V_pop 24.72088 5.977031e-06 TRUE
#> 2 KA 30.84953 ka_pop 30.84954 1.908690e-05 TRUE
#> 3 CL 11.94767 Cl_pop 11.94768 8.210502e-05 TRUE
#> 4 d_V 116.23095 var_omega_V 116.23100 4.371910e-05 TRUE
#> 5 d_KA 181.19778 var_omega_ka 181.19774 2.271892e-05 TRUE
#> 6 d_CL 77.29188 var_omega_Cl 77.29188 2.266350e-06 TRUEPediatric trial taking into account knowledge from adult trial:
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 V 6.379075 V_pop 6.379074 3.848623e-06 TRUE
#> 2 KA 8.219929 ka_pop 8.219932 2.814061e-05 TRUE
#> 3 CL 4.086855 Cl_pop 4.086855 1.122750e-05 TRUE
#> 4 d_V 31.808871 var_omega_V 31.808884 3.972277e-05 TRUE
#> 5 d_KA 52.858399 var_omega_ka 52.858388 2.094961e-05 TRUE
#> 6 d_CL 25.601551 var_omega_Cl 25.601515 1.415043e-04 TRUEEx12: Covariates
This example shows the usage of covariate distributions. The model includes the effect of weight on Cl and V.
popED
CODE
## Introduction
# Lets assume that we have a model with a covariate included in the model description.
# Here we define a one-compartment PK model that has weight on both clearance and volume of distribution.
mod_1 <- function(model_switch,xt,parameters,poped.db){
with(as.list(parameters),{
y=xt
CL=CL*(WT/70)^(WT_CL)
V=V*(WT/70)^(WT_V)
DOSE=1000*(WT/70)
y = DOSE/V*exp(-CL/V*xt)
return(list( y= y,poped.db=poped.db))
})
}
par_1 <- function(x,a,bpop,b,bocc){
parameters=c( CL=bpop[1]*exp(b[1]),
V=bpop[2]*exp(b[2]),
WT_CL=bpop[3],
WT_V=bpop[4],
WT=a[1])
return( parameters )
}
# distribution of covariates: We set
#`groupsize=1`, the number of groups to be 50 (`m=50`) and assume that WT is
#sampled from a normal distribution with mean=70 and sd=10
set.seed(123456)
weigths <- rnorm(50,mean = 70,sd=10)
poped_db_2 <- create.poped.database(ff_fun=mod_1,
fg_fun=par_1,
fError_fun=feps.add.prop,
groupsize=1,
m=50,
sigma=c(0.015,0.0015),
notfixed_sigma = c(1,0),
bpop=c(CL=3.8,V=20,WT_CL=0.75,WT_V=1),
d=c(CL=0.05,V=0.05),
xt=c( 1,2,4,6,8,24),
minxt=0,
maxxt=24,
bUseGrouped_xt=1,
a=as.list(weigths)
)
res12b_popED <- PopED::evaluate_design(poped_db_2)mlxDesignEval
The covariate effect is defined in the [INDIVIDUAL] block with a covariate transformation in the [COVARIATE] block to obtain a power-law relationship between the parameter and the covariate. Less standard parameter-covariate relationships can also be coded in the structural model using regressors, see our dedicated documentation.
CODE
model12 <- inlineModel("
[COVARIATE]
input = WT
EQUATION:
tWT = log(WT/70)
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, V_pop, omega_V, tWT, beta_V_tWT, beta_Cl_tWT}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, covariate=tWT, coefficient=beta_Cl_tWT, sd=omega_Cl}
V = {distribution=logNormal, typical=V_pop, covariate=tWT, coefficient=beta_V_tWT, sd=omega_V}
[LONGITUDINAL]
input = {a, b}
input = {V, Cl}
PK:
Cc = pkmodel(V, Cl)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")In R seed is fixed to use the same weights table in popED and mlxDesignEval. Instead of using one group per covariate value (as done in popED), we can simply define a data frame of weight values as covariate argument.
CODE
# table of weight values and weight-based dosing
set.seed(123456)
covweigths <- rnorm(50,mean = 70,sd=10)
res12b_mlx <- mlxDesignEval::evaluate_design(
model_file = model12,
population_parameters = data.frame(V_pop=20, Cl_pop=3.8,
beta_V_tWT=1,beta_Cl_tWT=0.75,
omega_V=sqrt(0.05), omega_Cl=sqrt(0.05),
a=sqrt(0.0015),b=sqrt(0.015)),
fixed_parameters = c("a"),
group = list(size=50),
covariate = data.frame(id=1:50, WT=covweigths),
treatment = list(data=data.frame(time=0, amount=1000/70),
scale=list(covariate="WT", intercept=0)),
output = list(output="DV", data=data.frame(time=c(1,2,4,6,8,24))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 3.249976 Cl_pop 3.249969 1.930655e-04 TRUE
#> 2 V 3.319989 V_pop 3.319985 1.253870e-04 TRUE
#> 3 WT_CL 29.001417 beta_Cl_tWT 29.001425 2.897006e-05 TRUE
#> 4 WT_V 22.228796 beta_V_tWT 22.228788 3.603646e-05 TRUE
#> 5 d_CL 21.026128 var_omega_Cl 21.026136 3.565315e-05 TRUE
#> 6 d_V 21.953885 var_omega_V 21.953879 2.690191e-05 TRUE
#> 7 SIGMA[1,1] 10.064637 var_b 10.064643 5.809050e-05 TRUEEx13: Shrinkage
It is not possible to calculate shrinkage with mlxDesignEval.
Ex14: PK, IOV
This example shows the usage of inter-occasion variability (IOV) on the clearance parameter.
popED
With popED, we need to define a different clearance parameter for each occasion and choose to use the first or second one using an if/else statement.
CODE
cppFunction(
'List one_comp_oral_ode(double Time, NumericVector A, NumericVector Pars) {
int n = A.size();
NumericVector dA(n);
double CL_OCC_1 = Pars[0];
double CL_OCC_2 = Pars[1];
double V = Pars[2];
double KA = Pars[3];
double TAU = Pars[4];
double N,CL;
N = floor(Time/TAU)+1;
CL = CL_OCC_1;
if(N>6) CL = CL_OCC_2;
dA[0] = -KA*A[0];
dA[1] = KA*A[0] - (CL/V)*A[1];
return List::create(dA);
}'
)
ff.ode.rcpp <- function(model_switch, xt, parameters, poped.db){
with(as.list(parameters),{
A_ini <- c(A1=0, A2=0)
times_xt <- drop(xt) #xt[,,drop=T]
dose_times = seq(from=0,to=max(times_xt),by=TAU)
eventdat <- data.frame(var = c("A1"),
time = dose_times,
value = c(DOSE), method = c("add"))
times <- sort(c(times_xt,dose_times))
out <- ode(A_ini, times, one_comp_oral_ode, c(CL_OCC_1,CL_OCC_2,V,KA,TAU),
events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
y = out[, "A2"]/(V)
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
sfg <- function(x,a,bpop,b,bocc){
parameters=c( CL_OCC_1=bpop[1]*exp(b[1]+bocc[1,1]),
CL_OCC_2=bpop[1]*exp(b[1]+bocc[1,2]),
V=bpop[2]*exp(b[2]),
KA=bpop[3]*exp(b[3]),
DOSE=a[1],
TAU=a[2])
return( parameters )
}
## -- Define design and design space
poped.db <- create.poped.database(ff_fun=ff.ode.rcpp,
fError_fun=feps.add.prop,
fg_fun=sfg,
bpop=c(CL=3.75,V=72.8,KA=0.25),
d=c(CL=0.25^2,V=0.09,KA=0.09),
sigma=c(0.04,5e-6),
notfixed_sigma=c(0,0),
docc = matrix(c(0,0.09,0),nrow = 1),
m=2,
groupsize=20,
xt=c( 1,2,8,240,245),
minxt=c(0,0,0,240,240),
maxxt=c(10,10,10,248,248),
bUseGrouped_xt=1,
a=list(c(DOSE=20,TAU=24),c(DOSE=40, TAU=24)),
maxa=c(DOSE=200,TAU=24),
mina=c(DOSE=0,TAU=24))
## evaluate initial design
res14_popED <- PopED::evaluate_design(poped.db)mlxDesignEval
The IOV is defined in the [INDIVIDUAL] block. Omega_Cl represents the IIV and gamma_Cl the IOV. Given the lognormal distribution, the formula for Cl is the same as in popED.
CODE
model14 <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, V_pop, omega_V, ka_pop, omega_ka, gamma_Cl}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, varlevel={id, id*occ}, sd={omega_Cl, gamma_Cl}}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
[LONGITUDINAL]
input = {a, b}
input = {ka, V, Cl}
PK:
Cc = pkmodel(ka, V, Cl)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=combined2(a, b)}")The occasions are defined in the occasion argument and repeated in the other elements as well.
CODE
treatment1<- list(data=data.frame(occ=c(1,1,1,1,1,1,2,2,2,2,2) ,
time=c(0,24,48,72,96,120,144,168,192,216,240),
amount=20))
treatment2<- list(data=data.frame(occ=c(1,1,1,1,1,1,2,2,2,2,2) ,
time=c(0,24,48,72,96,120,144,168,192,216,240),
amount=40))
g1 <- list(size=20, treatment=treatment1)
g2 <- list(size=20, treatment=treatment2)
res14_mlx <- mlxDesignEval::evaluate_design(
model_file = model14,
population_parameters = data.frame(ka_pop=0.25, V_pop=72.8, Cl_pop=3.75,
omega_ka=sqrt(0.09), omega_V=sqrt(0.09),
omega_Cl=0.25, gamma_Cl=sqrt(0.09),
a=sqrt(5e-6), b=sqrt(0.04)),
fixed_parameters = c("a","b"),
group = list(g1,g2),
occasion = data.frame(time=c(0,144),occ=c(1,2)),
output = list(output="DV", data=data.frame(occ=c(1,1,1,2,2),time=c(1,2,8,240,245))),
rse_on_variance = T)Comparison
The RSEs are identical, except a small difference for the IOV variance.
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 6.225029 Cl_pop 6.244430 0.31166092 TRUE
#> 2 V 8.817098 V_pop 8.826273 0.10406706 TRUE
#> 3 KA 10.672172 ka_pop 10.680824 0.08106658 TRUE
#> 4 d_CL 106.438935 var_omega_Cl 106.309736 0.12138269 TRUE
#> 5 d_V 42.915831 var_omega_V 42.942908 0.06309316 TRUE
#> 6 d_KA 62.884047 var_omega_ka 62.911048 0.04293874 TRUE
#> 7 D.occ[1,1] 79.419819 var_gamma_Cl 78.406349 1.27609228 FALSEEx15: full covariance matrix
This example shows the inclusion of a full covariance matrix, i.e correlations between the random effects.
popED
With popED, the covariance terms are defined in the covd argument of the database.
CODE
ff <- function(model_switch,xt,parameters,poped.db){
##-- Model: One comp first order absorption
with(as.list(parameters),{
y=xt
y=(DOSE*Favail*KA/(V*(KA-CL/V)))*(exp(-CL/V*xt)-exp(-KA*xt))
return(list(y=y,poped.db=poped.db))
})
}
sfg <- function(x,a,bpop,b,bocc){
## -- parameter definition function
parameters=c(CL=bpop[1]*exp(b[1]),
V=bpop[2]*exp(b[2]),
KA=bpop[3]*exp(b[3]),
Favail=bpop[4],
DOSE=a[1])
return(parameters)
}
feps <- function(model_switch,xt,parameters,epsi,poped.db){
## -- Residual Error function
## -- Proportional
returnArgs <- ff(model_switch,xt,parameters,poped.db)
y <- returnArgs[[1]]
poped.db <- returnArgs[[2]]
y = y*(1+epsi[,1])
return(list(y=y,poped.db=poped.db))
}
## -- With covariances
poped.db_with <- create.poped.database(ff_file="ff",
fg_file="sfg",
fError_file="feps",
bpop=c(CL=0.15, V=8, KA=1.0, Favail=1),
notfixed_bpop=c(1,1,1,0),
d=c(CL=0.07, V=0.02, KA=0.6),
covd = c(.03,.1,.09),
sigma=0.01,
groupsize=32,
xt=c( 0.5,1,2,6,24,36,72,120),
minxt=0,
maxxt=120,
a=70)
res15_popED <- PopED::evaluate_design(poped.db_with)mlxDesignEval
With mlxDesignEval, we define correlations instead of covariances, and this is done in the [INDIVIDUAL] block.
CODE
model15 <- inlineModel("
[INDIVIDUAL]
input = {Cl_pop, omega_Cl, F_pop, V_pop, omega_V, ka_pop, omega_ka, corr_ka_V, corr_V_Cl, corr_ka_Cl}
DEFINITION:
Cl = {distribution=logNormal, typical=Cl_pop, sd=omega_Cl}
F = {distribution=logNormal, typical=F_pop, no-variability}
V = {distribution=logNormal, typical=V_pop, sd=omega_V}
ka = {distribution=logNormal, typical=ka_pop, sd=omega_ka}
correlation = {level=id, r(V, Cl)=corr_V_Cl, r(ka, Cl)=corr_ka_Cl, r(ka, V)=corr_ka_V}
[LONGITUDINAL]
input = {b}
input = {F, ka, V, Cl}
PK:
Cc = pkmodel(ka, V, Cl, p=F)
OUTPUT:
output = {Cc}
DEFINITION:
DV = {distribution=normal, prediction=Cc, errorModel=proportional(b)}")To calculate the correlation value corresponding to the covariance value in popED, use:
corr_i_j = covar_i_j/(sqrt(var_i) * sqrt(var_j))
CODE
res15_mlx <- mlxDesignEval::evaluate_design(
model_file = model15,
population_parameters = data.frame(F_pop=1, Cl_pop=0.15,V_pop=8, ka_pop=1,
omega_Cl=sqrt(0.07),omega_V=sqrt(0.02),omega_ka=sqrt(0.6),
corr_V_Cl=0.8017837,corr_ka_Cl=0.4879500, corr_ka_V=0.8215838,
b=sqrt(0.01)),
fixed_parameters = c("F_pop"),
group = list(size=32),
treatment = list(data=data.frame(time=0, amount=70)),
output = list(output="DV", data=data.frame(time=c(0.5,1,2,6,24,36,72,120))),
rse_on_variance = T)Comparison
CODE
#> popED_name popED_RSE mlx_names mlx_RSE %_difference is_identical
#> 1 CL 4.738266 Cl_pop 4.738270 6.837835e-05 TRUE
#> 2 V 2.756206 V_pop 2.756208 4.766023e-05 TRUE
#> 3 KA 13.925829 ka_pop 13.925839 7.575014e-05 TRUE
#> 4 d_CL 25.660169 var_omega_Cl 25.660200 1.211581e-04 TRUE
#> 5 d_V 30.482032 var_omega_V 30.482018 4.601698e-05 TRUE
#> 6 d_KA 25.860002 var_omega_ka 25.860007 1.764386e-05 TRUE
#> 7 D[2,1] 30.823749 covar_V_Cl 30.823745 1.409248e-05 TRUE
#> 8 D[3,1] 41.469557 covar_ka_Cl 41.469561 1.042334e-05 TRUE
#> 9 D[3,2] 30.732069 covar_ka_V 30.732062 2.328111e-05 TRUE
#> 10 SIGMA[1,1] 11.180340 var_b 11.180340 5.004789e-12 TRUE