What is the best way to extract system information from roxde2 object

[john-harrold]

I need to programmatically find information about the names of different components of a system. I wanted to confirm what I'm currently doing is correct. I also wanted to ask about some I'm not sure of. Here are two examples one with odes and another using an analytical solution:

library(rxode2)

fun_odes  <- function() {
  description <- "Dupilumab PK model (Kovalenko 2020)"
  reference <- "Kovalenko P, Davis JD, Li M, et al. Base and Covariate Population Pharmacokinetic Analyses of Dupilumab Using Phase 3 Data. Clinical Pharmacology in Drug Development. 2020;9(6):756-767. doi:10.1002/cpdd.780"
  # Model 1 from table 1 and supplementary Table 2 in the publication and its
  # supplement.
  covariateData <-
    list(
      WT = "Body weight in kg"
    )
  ini({
    lvc <- log(2.48); label("central volume (L)")
    lke <- log(0.0534); label("elimination rate (1/d)")
    lkcp <- log(0.213); label("central-to-peripheral rate (1/d)")
    Mpc <- 0.686; label("ratio of kcp and kpc (kpc is peripheral to central rate with units of 1/d)")
    lka <- log(0.256); label("absorption rate (1/d)")
    lMTT <- log(0.105); label("mean transit time (d)")
    lVm <- log(1.07); label("maximum target-mediated rate of elimination (mg/L/d)")
    Km <- fixed(0.01); label("Michaelis-Menten constant (mg/L)")
    lfdepot <- log(0.643); label("Bioavailability (fraction)")
    e_wt_vc <- 0.711; label("Exponent of weight on central volume (unitless)")

    etalvc ~ 0.192
    etalke ~ 0.285
    etalka ~ 0.474
    etalvm ~ 0.236
    etamtt ~ 0.525 # etamtt is assumed to be on log-scale MTT to prevent negative values; this is a difference relative to Supplementary Table 2

    cppropSd <- 0.15; label("Proportional residual error (fraction)")
    cpaddSd <- fixed(0.03); label("Additive residual error (mg/L)")
  })
  model({
    # Weight normalization to 75 kg is assumed based on prior publication.  It
    # is not specified in the current publication:
    # Kovalenko P, DiCioccio AT, Davis JD, et al. Exploratory Population PK
    # Analysis of Dupilumab, a Fully Human Monoclonal Antibody Against
    # IL-4Ralpha, in Atopic Dermatitis Patients and Normal Volunteers. CPT
    # Pharmacometrics Syst Pharmacol. 2016;5(11):617-624. doi:10.1002/psp4.12136
    vc <- exp(lvc + etalvc)*(WT/75)^e_wt_vc
    ke <- exp(lke + etalke)
    kcp <- exp(lkcp)
    ka <- exp(lka + etalka)
    MTT <- exp(lMTT + etamtt)
    Vm <- exp(lVm + etalvm)

    # Derived parameters
    kpc <- kcp/Mpc
    ktr <- (3 + 1)/MTT

    d/dt(depot) <- -ktr*depot
    d/dt(transit1) <- ktr*(depot - transit1)
    d/dt(transit2) <- ktr*(transit1 - transit2)
    d/dt(transit3) <- ktr*transit2 - ka*transit3
    # Linear and Michaelis-Menten clearance
    d/dt(central) <-                 ka*transit3 - ke*central - kcp*central + kpc*periph - central*(Vm/(Km + central/vc))
    d/dt(periph) <-                                             kcp*central - kpc*periph

    f(depot) <- exp(lfdepot)
    # No unit conversion is required to change mg/L (dosing amount/central
    # volume unit) to mg/L (measurement unit)
    Cc <- central/vc
    Cc ~ add(cpaddSd) + prop(cppropSd)
  })
}

fun_ana <- function() {
  description <- "Two compartment PK model with linear clearance for average monoclonal antibodies (Davda 2014)"
  reference <- "Davda JP, Dodds MG, Gibbs MA, Wisdom W, Gibbs JP. A model-based meta-analysis of monoclonal antibody pharmacokinetics to guide optimal first-in-human study design. MAbs. 2014;6(4):1094-1102. doi:10.4161/mabs.29095"
  ini({
    lfdepot <- log(0.744) ; label("Subcutaneous bioavailability (fraction)")
    lka <- log(0.282) ; label("Absorption rate (Ka, 1/day)")
    lcl <- log(0.200) ; label("Clearance (CL, L/day)")
    lv  <- log(3.61) ; label("Central volume of distribution (V, L)")
    lvp  <- log(2.75) ; label("Peripheral volume of distribution (Vp, L)")
    lq  <- log(0.747) ; label("Intercompartmental clearance (Q, L/day)")

    allocl <- 0.865 ; label("Allometric exponent on clearance and intercompartmental clearance (unitless)")
    allov <- 0.957 ; label("Allometric exponent on volumes of distribution (unitless)")

    etafdepot ~ 0
    etaka ~ 0.416
    etacl + etav + etavp ~ c(0.0987,
                             0.0786, 0.116,
                             0.0377, 0.0619, 0.0789)
    etaq ~ 0.699

    prop.err <- 0.144 ; label("Proportional residual error (fraction)")
  })
  model({
    # WT is body weight in kg
    fdepot <- exp(lfdepot + etafdepot)
    ka <- exp(lka + etaka)
    wtnorm <- log(WT/70)
    cl <- exp(lcl + allocl*wtnorm + etacl)
    q  <- exp(lq + allocl*wtnorm + etaq)
    v <- exp(lv + allov*wtnorm + etav)
    vp <- exp(lvp + allov*wtnorm + etavp)
    Cc <- linCmt()
    
    f(depot) <- fdepot  # Units are dosing units/L (typically mg/L = ug/mL)
    Cc ~ prop(prop.err)
    
  })
}

obj_ana = rxode2(fun_ana)
obj_odes = rxode2(fun_odes)

For covariates, parameters, and iiv terms I'm using this which seems to work. Is there anything wrong with it?

# Covariates
covariates_ana  = obj_ana$allCovs
covariates_odes = obj_odes$allCovs

# System, input, assignment values in ini
parameters_ana  = names(obj_ana$theta)
parameters_odes = names(obj_odes$theta)

# IIV terms
iiv_ana         = unique(names(obj_ana$etaLhs))
iiv_odes        = unique(names(obj_odes$etaLhs))

I also want a list of states or compartments in the model. For ODE based systems this seems to work:

states_odes     = obj_odes$state

But I don't know how to get states/compartments for analytical solutions. What I"m looking for here is compartments that can be dosed into. So if there is a test where I can do something like if length(obj_ana$state) == 0 then here are default dosing compartments. But I need to know how to do that test and what dosing compartments there are.

Next what I'm looking for is secondary parameters. These are values assigned in the model() chunk of the code. I'm using this but it seems a little tenuous because it looks like a part of the object tracking only the relationship between typical values and individual parameters.

sec_parameters_ana = unique(names(obj_ana$lhsTheta))

Lastly I need a list of outputs. Basically things assigned with a ~ in the model() chunk of the code.

[billdenney]

I think that a general function that gets the information about a model would be very helpful. That way, it could be kept as a consistent API while the internals that you're digging into could change.

[mattfidler]

rxModelVars() is the most basic function that gets model variables, but I think mostly John is working with the UI, which has some additional hurdles.

I hope to comment on this today or tomorrow.

[john-harrold]

There isn't any rush on this.

I'm using the rxode2 object that is built from the functions. I left it out above but I'm building it with:

obj_ana = rxode2(fun_ana)
obj_odes = rxode2(fun_odes)

I looked at rxModelVars but params seems to contain both the system parameters as well as etas. So I'm trying to break them out up there.

I think lhs will give me all of the secondary parameters.

It also doesn't really tell me what to do about the states of solved systems and defined outputs as well.

I agree with Bill. It would be great to have a simple function that would do this :).

[mattfidler]

I think that everything should be in the UI. I am struggling with maintaining an additional function to reapply what is already in the UI.

If things are not easily accessible, simply make a way to access them easily and keep them the same for the API.

[mattfidler]

Perhaps simply a way to access them from the UI

That would be something like $params which would list the relevant parameters.

[mattfidler]

Can you clairfy what you mean by:

Next what I'm looking for is secondary parameters. These are values assigned in the model() chunk of the code. I'm using this but it seems a little tenuous because it looks like a part of the object tracking only the relationship between typical values and individual parameters.

The model simulation currently calculates the relationship between typical and population parameters, Did you want a list of output parameters without this calculation?

Lastly I need a list of outputs. Basically things assigned with a ~ in the model() chunk of the code.

Do you mean the endpoints in the model (ie the derived variables that go into the endpoints)?

[john-harrold]

To me secondary parameters are any parameters defined in terms of your typical value/etas (anything in the ini section) which would include population parameters but also any other values that are defined in terms of secondary parameters. So for example kpc and qtr below:

kpc <- kcp/Mpc
ktr <- (3 + 1)/MTT

To me outputs are anything you're comparing to data. So I think that would be anything in the model section that is associated with error using a ~. Like Cc below:

    Cc ~ prop(prop.err)

[mattfidler]

I hope I understood correctly.

A current pull request seems to do what you wish with $params:

library(rxode2)
#> rxode2 2.1.1.9000 using 8 threads (see ?getRxThreads)
#>   no cache: create with `rxCreateCache()`

fun_odes  <- function() {
  description <- "Dupilumab PK model (Kovalenko 2020)"
  reference <- "Kovalenko P, Davis JD, Li M, et al. Base and Covariate Population Pharmacokinetic Analyses of Dupilumab Using Phase 3 Data. Clinical Pharmacology in Drug Development. 2020;9(6):756-767. doi:10.1002/cpdd.780"
  # Model 1 from table 1 and supplementary Table 2 in the publication and its
  # supplement.
  covariateData <-
    list(
      WT = "Body weight in kg"
    )
  ini({
    lvc <- log(2.48); label("central volume (L)")
    lke <- log(0.0534); label("elimination rate (1/d)")
    lkcp <- log(0.213); label("central-to-peripheral rate (1/d)")
    Mpc <- 0.686; label("ratio of kcp and kpc (kpc is peripheral to central rate with units of 1/d)")
    lka <- log(0.256); label("absorption rate (1/d)")
    lMTT <- log(0.105); label("mean transit time (d)")
    lVm <- log(1.07); label("maximum target-mediated rate of elimination (mg/L/d)")
    Km <- fixed(0.01); label("Michaelis-Menten constant (mg/L)")
    lfdepot <- log(0.643); label("Bioavailability (fraction)")
    e_wt_vc <- 0.711; label("Exponent of weight on central volume (unitless)")

    etalvc ~ 0.192
    etalke ~ 0.285
    etalka ~ 0.474
    etalvm ~ 0.236
    etamtt ~ 0.525 # etamtt is assumed to be on log-scale MTT to prevent negative values; this is a difference relative to Supplementary Table 2

    cppropSd <- 0.15; label("Proportional residual error (fraction)")
    cpaddSd <- fixed(0.03); label("Additive residual error (mg/L)")
  })
  model({
    # Weight normalization to 75 kg is assumed based on prior publication.  It
    # is not specified in the current publication:
    # Kovalenko P, DiCioccio AT, Davis JD, et al. Exploratory Population PK
    # Analysis of Dupilumab, a Fully Human Monoclonal Antibody Against
    # IL-4Ralpha, in Atopic Dermatitis Patients and Normal Volunteers. CPT
    # Pharmacometrics Syst Pharmacol. 2016;5(11):617-624. doi:10.1002/psp4.12136
    vc <- exp(lvc + etalvc)*(WT/75)^e_wt_vc
    ke <- exp(lke + etalke)
    kcp <- exp(lkcp)
    ka <- exp(lka + etalka)
    MTT <- exp(lMTT + etamtt)
    Vm <- exp(lVm + etalvm)

    # Derived parameters
    kpc <- kcp/Mpc
    ktr <- (3 + 1)/MTT

    d/dt(depot) <- -ktr*depot
    d/dt(transit1) <- ktr*(depot - transit1)
    d/dt(transit2) <- ktr*(transit1 - transit2)
    d/dt(transit3) <- ktr*transit2 - ka*transit3
    # Linear and Michaelis-Menten clearance
    d/dt(central) <-                 ka*transit3 - ke*central - kcp*central + kpc*periph - central*(Vm/(Km + central/vc))
    d/dt(periph) <-                                             kcp*central - kpc*periph

    f(depot) <- exp(lfdepot)
    # No unit conversion is required to change mg/L (dosing amount/central
    # volume unit) to mg/L (measurement unit)
    Cc <- central/vc
    Cc ~ add(cpaddSd) + prop(cppropSd)
  })
}

fun_ana <- function() {
  description <- "Two compartment PK model with linear clearance for average monoclonal antibodies (Davda 2014)"
  reference <- "Davda JP, Dodds MG, Gibbs MA, Wisdom W, Gibbs JP. A model-based meta-analysis of monoclonal antibody pharmacokinetics to guide optimal first-in-human study design. MAbs. 2014;6(4):1094-1102. doi:10.4161/mabs.29095"
  ini({
    lfdepot <- log(0.744) ; label("Subcutaneous bioavailability (fraction)")
    lka <- log(0.282) ; label("Absorption rate (Ka, 1/day)")
    lcl <- log(0.200) ; label("Clearance (CL, L/day)")
    lv  <- log(3.61) ; label("Central volume of distribution (V, L)")
    lvp  <- log(2.75) ; label("Peripheral volume of distribution (Vp, L)")
    lq  <- log(0.747) ; label("Intercompartmental clearance (Q, L/day)")

    allocl <- 0.865 ; label("Allometric exponent on clearance and intercompartmental clearance (unitless)")
    allov <- 0.957 ; label("Allometric exponent on volumes of distribution (unitless)")

    etafdepot ~ 0
    etaka ~ 0.416
    etacl + etav + etavp ~ c(0.0987,
                             0.0786, 0.116,
                             0.0377, 0.0619, 0.0789)
    etaq ~ 0.699

    prop.err <- 0.144 ; label("Proportional residual error (fraction)")
  })
  model({
    # WT is body weight in kg
    fdepot <- exp(lfdepot + etafdepot)
    ka <- exp(lka + etaka)
    wtnorm <- log(WT/70)
    cl <- exp(lcl + allocl*wtnorm + etacl)
    q  <- exp(lq + allocl*wtnorm + etaq)
    v <- exp(lv + allov*wtnorm + etav)
    vp <- exp(lvp + allov*wtnorm + etavp)
    Cc <- linCmt()

    f(depot) <- fdepot  # Units are dosing units/L (typically mg/L = ug/mL)
    Cc ~ prop(prop.err)

  })
}

fun_ana2 <- function() {
  description <- "Two compartment PK model with linear clearance for average monoclonal antibodies (Davda 2014)"
  reference <- "Davda JP, Dodds MG, Gibbs MA, Wisdom W, Gibbs JP. A model-based meta-analysis of monoclonal antibody pharmacokinetics to guide optimal first-in-human study design. MAbs. 2014;6(4):1094-1102. doi:10.4161/mabs.29095"
  ini({
    lfdepot <- log(0.744) ; label("Subcutaneous bioavailability (fraction)")
    lka <- log(0.282) ; label("Absorption rate (Ka, 1/day)")
    lcl <- log(0.200) ; label("Clearance (CL, L/day)")
    lv  <- log(3.61) ; label("Central volume of distribution (V, L)")
    lvp  <- log(2.75) ; label("Peripheral volume of distribution (Vp, L)")
    lq  <- log(0.747) ; label("Intercompartmental clearance (Q, L/day)")

    allocl <- 0.865 ; label("Allometric exponent on clearance and intercompartmental clearance (unitless)")
    allov <- 0.957 ; label("Allometric exponent on volumes of distribution (unitless)")

    etafdepot ~ 0
    etaka ~ 0.416
    etacl + etav + etavp ~ c(0.0987,
                             0.0786, 0.116,
                             0.0377, 0.0619, 0.0789)
    etaq ~ 0.699

    prop.err <- 0.144 ; label("Proportional residual error (fraction)")
  })
  model({
    # WT is body weight in kg
    fdepot <- exp(lfdepot + etafdepot)
    ka <- exp(lka + etaka)
    wtnorm <- log(WT/70)
    cl <- exp(lcl + allocl*wtnorm + etacl)
    q  <- exp(lq + allocl*wtnorm + etaq)
    v <- exp(lv + allov*wtnorm + etav)
    vp <- exp(lvp + allov*wtnorm + etavp)
    linCmt() ~ prop(prop.err)
    f(depot) <- fdepot  # Units are dosing units/L (typically mg/L = ug/mL)

  })
}

obj_ana = rxode2(fun_ana)
obj_ana2 = fun_ana2()
obj_odes = rxode2(fun_odes)
#> ℹ parameter labels from comments will be replaced by 'label()'


obj_ana$params
#> $pop
#> [1] "lfdepot" "lka"     "lcl"     "lv"      "lvp"     "lq"      "allocl" 
#> [8] "allov"  
#> 
#> $resid
#> [1] "prop.err"
#> 
#> $group
#> $group$id
#> [1] "etafdepot" "etaka"     "etacl"     "etav"      "etavp"     "etaq"     
#> 
#> 
#> $linCmt
#> [1] TRUE
#> 
#> $cmt
#> [1] "depot"   "central"
#> 
#> $output
#> $output$primary
#> [1] "fdepot" "ka"     "wtnorm" "cl"     "q"      "v"      "vp"    
#> 
#> $output$secondary
#> character(0)
#> 
#> $output$endpoint
#> [1] "Cc"
#> 
#> $output$state
#> character(0)
obj_ana2$params
#> $pop
#> [1] "lfdepot" "lka"     "lcl"     "lv"      "lvp"     "lq"      "allocl" 
#> [8] "allov"  
#> 
#> $resid
#> [1] "prop.err"
#> 
#> $group
#> $group$id
#> [1] "etafdepot" "etaka"     "etacl"     "etav"      "etavp"     "etaq"     
#> 
#> 
#> $linCmt
#> [1] TRUE
#> 
#> $cmt
#> [1] "depot"   "central"
#> 
#> $output
#> $output$primary
#> [1] "fdepot" "ka"     "wtnorm" "cl"     "q"      "v"      "vp"    
#> 
#> $output$secondary
#> character(0)
#> 
#> $output$endpoint
#> character(0)
#> 
#> $output$state
#> character(0)
obj_odes$params
#> $pop
#>  [1] "lvc"     "lke"     "lkcp"    "Mpc"     "lka"     "lMTT"    "lVm"    
#>  [8] "Km"      "lfdepot" "e_wt_vc"
#> 
#> $resid
#> [1] "cppropSd" "cpaddSd" 
#> 
#> $group
#> $group$id
#> [1] "etalvc" "etalke" "etalka" "etalvm" "etamtt"
#> 
#> 
#> $linCmt
#> [1] FALSE
#> 
#> $cmt
#> [1] "depot"    "transit1" "transit2" "transit3" "central"  "periph"  
#> 
#> $output
#> $output$primary
#> [1] "vc"  "ke"  "kcp" "ka"  "MTT" "Vm"  "kpc"
#> 
#> $output$secondary
#> [1] "ktr"
#> 
#> $output$endpoint
#> [1] "Cc"
#> 
#> $output$state
#> [1] "depot"    "transit1" "transit2" "transit3" "central"  "periph"

^{Created on 2024-01-12 with reprex v2.0.2}

There are a few things to keep in mind:

• ktr is the only secondary parameter for the ode system. I define secondary parameters by if they depend directly on input population or between subject variabilities directly. Since kpc is calculated in terms of the population parameter Mpc, it does not qualify.
• The endpoints may not be defined (as in the case of linCmt() ~
• States may be empty though dosing compartment $cmt can contain items
• In models without linCmt() or ODEs, there will be no dosing compartments

[john-harrold]

I think this is what I"m looking for. I'm scanning through the examples above and I only see one issue. Shouldn't obj_odes$params$secondary also contain kpc as well?

[mattfidler]

No. It is defined in terms of an estimated population parameter

[john-harrold]

Are you thinking of kcp?

   kcp <- exp(lkcp)
   
    # Derived parameters
    kpc <- kcp/Mpc
    ktr <- (3 + 1)/MTT

[mattfidler]

No. As documented in the original fix, kpc depends on the population parameter Mpc

You can see it in the ini block:

ini({
...
Mpc <- 0.686; label("ratio of kcp and kpc (kpc is peripheral to central rate with units of 1/d)")
...
})

[john-harrold]

I see so if it's defined in terms of something in ini then it's found in obj_odes$params$output$primary. Would everything else be in obj_odes$params$output$secondary. For example bob and sam below:

bob = kpc*2
sam = Cc*kpc

[mattfidler]

Unless bob or sam are a endpoint variable, Yes.

[john-harrold]

Cool this should work then. Thanks.