Take an example of the pk.3.comp.olal.ode model

[liangqiongyue]

I have encountered a problem and need to use the 3-compartment oral model for modeling, but I did not find any example of this model on the Internet, so I would like to ask you if there is any problem with my code. If there is a problem, I hope you can help me modify it, thank you very much.
Here is my code:

#' An implementation of a two compartment model with oral absorption using ODEs.
library(PopED)
library(deSolve)

#' Define the ODE system
PK.3.comp.oral.ode <- function(Time, State, Pars){
  with(as.list(c(State, Pars)), {    
    dA1 <- -KA*A1 
    dA2 <- KA*A1 + A3* Q3/V3 -A2*(CL/V2+Q3/V2+Q4/V2)+A4*Q4/V4
    dA3 <- A2* Q3/V2-A3* Q3/V3
    dA4 <- A2* Q4/V2-A4* Q4/V4
    return(list(c(dA1, dA2, dA3, dA4)))
  })
}

#' define the initial conditions and the dosing
ff.PK.3.comp.oral.md.ode <- function(model_switch, xt, parameters, poped.db){
  with(as.list(parameters),{
    A_ini <- c(A1=0, A2=0, A3=0, A4=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.3.comp.oral.ode, parameters, events = list(data = eventdat))#atol=1e-13,rtol=1e-13)
    y = out[, "A2"]/(V2/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]),
                V2=bpop[2]*exp(b[2]),
                V3=bpop[3]*exp(b[3]),
                V4=bpop[4]*exp(b[4]),
                KA=bpop[5]*exp(b[5]),
                Q3=bpop[6]*exp(b[6]),
                Q4=bpop[7]*exp(b[7]),
                Favail=bpop[8]*exp(b[8]),
                ALAG1=bpop[9]*exp(b[9]),
                DOSE=a[1],
                TAU=a[2])
  return( parameters ) 
}

#' create poped database
poped.db <- create.poped.database(ff_fun="ff.PK.3.comp.oral.md.ode",
                                  fError_fun="feps.add.prop",
                                  fg_fun="fg",
                                  groupsize=20,
                                  m=1,      #number of groups
                                  sigma=c(prop=0.0714^2,add=6.13^2),                                
                                  bpop=c(CL=42.9, V2=102, V3=2340, V4=421, KA=3.54, Q3=105, Q4=397, Favail=1, ALAG1=0.0626),    (bpop)
                                  d=c(CL=0.487^2, V2=0.385^2, V3=0.196^2, V4=0.427^2, KA=0, Q3=0.319^2, Q4=0.450^2, Favail=0, ALAG1=0.101^2),  
                                  notfixed_bpop=c(1,1,1,1,1,1,1,0,1),
                                  xt=c(0, 7*24, 7*24+5/60, 7*24+20/60, 7*24+1, 7*24+12, 7*24+24),
                                  minxt=c(-1, 7*24-1, 7*24+0, 7*24+10/60, 7*24+0.5, 7*24+10, 7*24+20),
                                  maxxt=c(0, 7*24, 7*24+10/60, 7*24+0.5, 7*24+1.5, 7*24+14, 7*24+28),
                                  a=c(DOSE=75,TAU=24),
                                  mina=c(DOSE=37.5,TAU=24),
                                  maxa=c(DOSE=150,TAU=24),
                                  discrete_a = list(DOSE=seq(0,150,by=37.5),TAU=24))

#' plot intial design just PRED
plot_model_prediction(poped.db,model_num_points = 500)