Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

[F] Better handling analytical PK models #34

Closed
kylebaron opened this issue Apr 25, 2016 · 16 comments
Closed

[F] Better handling analytical PK models #34

kylebaron opened this issue Apr 25, 2016 · 16 comments

Comments

@kylebaron
Copy link
Collaborator

library(knitr)
opts_chunk$set(fig.path="img/evid4-")

library(mrgsolve)
## mrgsolve: Community Edition

## www.github.com/metrumresearchgroup
options(mrgsolve_mread_quiet=TRUE)

Proposed change for analytical PK models

The idea is this: specify ADVAN and TRANS in $SUBROUTINES. Depending on what you pick for TRANS, mrgsolve will automatically capture the "fundamental" PK parameters for that ADVAN/TRANS.

For trans 2, the fundamental parameters are CL, V, and KA. We can utilize items in $PARAM here.

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
$SUBROUTINES advan=2, trans=2
'
mod <- mcode("advan2", code)
mod %>% ev(amt=1000, rate=25,cmt=2) %>% mrgsim(end=120)
## Model:  advan2.cpp 
## Dim:    122 x 4 
## Time:   0 to 120 
## ID:     1 
##      ID time GUT   CENT
## [1,]  1    0   0   0.00
## [2,]  1    0   0   0.00
## [3,]  1    1   0  24.65
## [4,]  1    2   0  48.60
## [5,]  1    3   0  71.88
## [6,]  1    4   0  94.50
## [7,]  1    5   0 116.48
## [8,]  1    6   0 137.85

Another trans (11) puts a lower case i at the end. We could draw from $PARAM again, but in this example, we use a derived variable.

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
$SUBROUTINES advan=2, trans=11

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
$CAPTURE ECL EV
'
mod <- mcode("advan2b", code)
set.seed(110102)
mod %>% ev(amt=1000, rate=25,cmt=2) %>% mrgsim(end=120)
## Model:  advan2b.cpp 
## Dim:    122 x 6 
## Time:   0 to 120 
## ID:     1 
##      ID time GUT   CENT     ECL     EV
## [1,]  1    0   0   0.00 -0.3357 -1.253
## [2,]  1    0   0   0.00 -0.3357 -1.253
## [3,]  1    1   0  24.13 -0.3357 -1.253
## [4,]  1    2   0  46.59 -0.3357 -1.253
## [5,]  1    3   0  67.50 -0.3357 -1.253
## [6,]  1    4   0  86.97 -0.3357 -1.253
## [7,]  1    5   0 105.10 -0.3357 -1.253
## [8,]  1    6   0 121.97 -0.3357 -1.253

A similar deal for ADVAN 4

code <- '
$PARAM CL = 1, V2 = 35, V3 = 200, Q = 12, KA = 1
$CMT GUT CENT PERIPH
$SUBROUTINES advan=4, trans=4

'
mod <- mcode("advan2c", code)
set.seed(110102)
mod %>% ev(amt=1000, rate=25,cmt=2) %>% mrgsim(end=120)
## Model:  advan2c.cpp 
## Dim:    122 x 5 
## Time:   0 to 120 
## ID:     1 
##      ID time GUT  CENT PERIPH
## [1,]  1    0   0  0.00   0.00
## [2,]  1    0   0  0.00   0.00
## [3,]  1    1   0 20.95   3.73
## [4,]  1    2   0 35.75  13.11
## [5,]  1    3   0 46.53  26.14
## [6,]  1    4   0 54.69  41.54
## [7,]  1    5   0 61.13  58.44
## [8,]  1    6   0 66.45  76.30

Just leave me alone

code <- '
$PARAM a = 1, b = 35, c = 12, d = 200, e = 1
$CMT GUT CENT PERIPH
$SUBROUTINES advan=4, trans=1
$MAIN
pred_CL = a;
pred_V2 = b;
pred_Q = c;
pred_V3 = d;
pred_KA = e;

'
mod <- mcode("advan2d", code)
set.seed(110102)
mod %>% ev(amt=1000, rate=25,cmt=2) %>% mrgsim(end=120)
## Model:  advan2d.cpp 
## Dim:    122 x 5 
## Time:   0 to 120 
## ID:     1 
##      ID time GUT  CENT PERIPH
## [1,]  1    0   0  0.00   0.00
## [2,]  1    0   0  0.00   0.00
## [3,]  1    1   0 20.95   3.73
## [4,]  1    2   0 35.75  13.11
## [5,]  1    3   0 46.53  26.14
## [6,]  1    4   0 54.69  41.54
## [7,]  1    5   0 61.13  58.44
## [8,]  1    6   0 66.45  76.30

What's included?

  1. ADVAN 1 through 4 (as defined by NONMEM)
  2. TRANS
  • 1 = nothing
  • 2 = CL/V/KA or CL/V2/Q/V3/KA
  • 11 = CLi/Vi/KAi or CLi/V2i/Qi/V3i/KAi
  1. You can still interact directly with pred_CL, pred_V3 etc ... the TRANS stuff is just convenience to save you typing.
  2. I'm not planning a ton of error checking to make sure you have valid symbols defined given your selection for TRANS. But if you don't, the compiler will give a pretty clear error that Vi (for example) wasn't found.
  3. I'd like to eventually deprecate $ADVAN2 and $ADVAN4. For now, they should work as they always have.
@dpastoor
Copy link
Contributor

Looks like a nice set of ideas. I'm a little confused regarding the CL vs CLi distinction. For example, given this trans11 example:

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
$SUBROUTINES advan=2, trans=11

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
$CAPTURE ECL EV
'

Wouldn't this be the same as just using trans2 and explicitly capturing? :

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
-$SUBROUTINES advan=2, trans=11
+$SUBROUTINES advan=2, trans=2

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
-$CAPTURE ECL EV
+$CAPTURE ECL EV CLi Vi KAi
'

aka the 'benefit' is your <PARAM>i values are automagically captured rather than explicitly defining them? Am I missing some other distinction between the 2?

Secondly, it seems trans1 could just be a default?

code <- '
$PARAM a = 1, b = 35, c = 12, d = 200, e = 1
$CMT GUT CENT PERIPH
-$SUBROUTINES advan=4, trans=1
+$SUBROUTINES advan=4

$MAIN
pred_CL = a;
pred_V2 = b;
pred_Q = c;
pred_V3 = d;
pred_KA = e;
'

@kylebaron
Copy link
Collaborator Author

Yes trans = 1 is the default. I explicitly set it, but you are right ... it wasn't necessary to do.

The main benefit / motivation for doing this is getting the user out of having to write stuff like this:

pred_CL = a;
pred_V2 = b;
pred_Q = c;
pred_V3 = d;
pred_KA = e;

Once you commit to, say, trans 2, you are saying that you will have CL, V and KA defined and they contain the values you want to use to calculate the system. mrgsolve basically connects CL with pred_CL, V with pred_V, etc .... Not only is it less coding for the user but hopefully more familiar for those coming from NONMEM.

It's kind of a fix for a self-inflicted wound: I fooled around with 3 or 4 different ways to get CL, V, and KA from the model dll back to the mrgsolve dll. Also, toyed around with putting the machinery to calculate the solution in the model dll, but thought it would be more trouble than it was worth. I came up with (I think) something safe and easier to keep track of, but potentially required the user to do something like this:

$PARAM CL=1, V=2, KA=3
$ADVAN2
$MAIN
pred_CL = CL;
pred_V = V;
pred_KA = KA;

which seems really wasteful when you should be able to do this:

$PARAM CL=1, V=2, KA=3
$SUB advan=2, trans=2

The trans 1, 2, and 4 are as NONMEM does. I invented trans 13 because I tend to code like that and wanted to have it around.

Posted the proposal to see if any other unforeseen consequences. I want to deprecate $ADVAN2 and $ADVAN4 if switching to $SUB makes sense.

@kylebaron
Copy link
Collaborator Author

Another try: drop $SUBROUTINES (too NONMEM). Introduce $PKMODEL which for now does the same thing, but you specify the number of compartments (ncmt) and a logical about whether you want a depot compartment or not.

This still has trans = 1 (nothing) trans=2 (CL/V/KA) and trans=4 (CL/V2/Q/V3/KA).

Wondering about changing trans to be either "CLV" ... "clearances and volumes" ... mrgsolve will pick the right parameter set based on the number of compartments chosen (most users doing this) or "CLVi" ... clearances and volumes with i appended at the end (still think this would be useful.

code <- '
$CMT GUT CENT 
$PKMODEL ncmt=1, trans=2, depot=TRUE
$PARAM CL=1, V = 20, KA=1
'

mod <-mcode("model", code)
## Compiling model.cpp.cpp ...

## done.

## Loading: model5fdc1214e14.so
mod %>% ev(amt=100) %>% mrgsim
## Model:  model.cpp 
## Dim:    26 x 4 
## Time:   0 to 24 
## ID:     1 
##      ID time      GUT  CENT
## [1,]  1    0   0.0000  0.00
## [2,]  1    0 100.0000  0.00
## [3,]  1    1  36.7879 61.41
## [4,]  1    2  13.5335 81.00
## [5,]  1    3   4.9787 85.36
## [6,]  1    4   1.8316 84.25
## [7,]  1    5   0.6738 81.27
## [8,]  1    6   0.2479 77.72

@dpastoor
Copy link
Contributor

haha I was already thinking about pkmodel - do you think we could even completely drop the code block, and make some assumptions about what they can do based on the data

pkmodel(list(CL = 1, V = 20, KA = 1)) %>% ev(amt = 100) %>% mrgsim

pkmodel could quickly check for invalid param combos so

pkmodel(list(pred_CL, pred_V, pred_ka)) %>% ...

could fail with like

Error: invalid param combination, pred_ka not a valid parameter name

This is of course, assuming that the flexibility of these analytical models is still low enough to handle this.

@vjd
Copy link

vjd commented Apr 26, 2016

The problem with dropping code blocks I think has to deal with
customization of the models once the base model is built. Any further
complexity in terms of covariates or variability models has to be coded
anyway.

On Tue, Apr 26, 2016 at 5:37 PM Devin Pastoor notifications@github.com
wrote:

haha I was already thinking about pkmodel - do you think we could even
completely drop the code block, and make some assumptions about what they
can do based on the data

pkmodel(list(CL = 1, V = 20, KA = 1)) %>% ev(amt = 100) %>% mrgsim

pkmodel could quickly check for invalid param combos so

pkmodel(list(pred_CL, pred_V, pred_ka)) %>% ...

could fail with like

Error: invalid param combination, pred_ka not a valid parameter name

This is of course, assuming that the flexibility of these analytical
models is still low enough to handle this.


You are receiving this because you are subscribed to this thread.
Reply to this email directly or view it on GitHub
#34 (comment)

~Vijay

@dpastoor
Copy link
Contributor

dpastoor commented Apr 26, 2016

and to your previous post - absolutely get (and really like) the idea of providing those shorthands - I still am just trying to wrap my head around CL vs CLi distinction that trans11 provides. Does it provide additional error checking or something? Wouldn't CLi Vi etc essentially always need to be defined 'customly' given the expectation of a random effects structure being added (else they'd just be using trans2)

As in, what is the benefit of this:

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
+$SUBROUTINES advan=2, trans=11

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
$CAPTURE ECL EV
'

over

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
+$SUBROUTINES advan=2, trans=2

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
+$CAPTURE ECL EV CLi Vi KAi
'

Other than not having to explicitly capture CLi Vi? Is there addtional error checking optimization etc?

@kylebaron
Copy link
Collaborator Author

@dpastoor wrote:

pkmodel(list(CL = 1, V = 20, KA = 1)) %>% ev(amt = 100) %>% mrgsim

We could totally do this; you wouldn't need to compile anything, but let odeproblem::advan2 do all of the calculating back in the mrgsolve dll. As @vjd noted you couldn't do anything complicated like covariate model (or population stuff) ... but it could be useful / cool do to for teaching "two-compartment model" or whatever.

@kylebaron
Copy link
Collaborator Author

Sorry for confusion on: trans 2 / trans 11 ... or ... CL/CLi

All it does is give the user more flexibility around the names of the variables. Whether you use trans 2 CL or trans 11 CLi, they both contain the value of if the parameter that you want to use to advance the system. Either one may or may not have some random variability with it.

@dpastoor wrote:

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
+$SUBROUTINES advan=2, trans=2

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
+$CAPTURE ECL EV CLi Vi KAi
'

Since trans is 2, mrgsolve will advance the system with CL, not CLi. Probably not what you want.

I was showing Vijay this today ... When you use trans = 2, mrgsolve quitely adds this to the very end of $MAIN:

pred_CL = CL;  pred_V  = V;   pred_KA = KA;

as if you yourself wrote it. It uses CL/V/KA because trans is 2. This will be wrong if using the model at the top of this post.

If, instead, you used trans=11, you'd get this at the bottom of $MAIN

pred_CL = CLi;  pred_V  = Vi;   pred_KA = KAi;

That's better ... we want to advance with the individual, not population, parameters. Either way, the i at the end doesn't have anything to do with it; it's just that you have the individual parameter in CLi not in CL in this case. You can either set trans to 11 or re-write your model so that CL has the individual parameter value.

I wouldn't say there is and advantage of 11 over 2. Only that, if you had CL in $PARAM AND trans = 2 AND some covariate model or random effect on clearance, there would be no way to advance with the individual parameter.

@kylebaron
Copy link
Collaborator Author

Here's the latest stab:

code <- '
$CMT GUT  CENT
$PKMODEL ncmt=1,depot=TRUE

$PARAM CL=1, V = 20, KA=1
'

Still calling this advan/trans, but only under the hood. Hope to spare the user with having to think about it.

##' Parse data from \code{$PKMODEL}.
##'
##' @param ncmt number of compartments
##' @param depot logical indicating whether to add depot/dosing compartment
##' @param trans the parameterization for the PK model
PKMODEL <- function(ncmt=1, depot=FALSE, trans = pick_trans(ncmt,depot), ...) {
    stopifnot(ncmt %in% c(1,2))
    advan <- pick_advan(ncmt,depot)
    return(list(advan=advan, trans=trans, n=ncmt))
}

We want 1-compartment model with depot dosing; trans is assumed to be 2 for one-compartment model we need to supply CL, V and KA. The default is trans=4 when ncmt is 2. Hopefully, people don't need to even know what trans is and they just accept the default. But it's there

@kylebaron kylebaron changed the title PROPOSAL for better handling analytical PK models [F] Better handling analytical PK models Apr 27, 2016
kylebaron added a commit that referenced this issue Apr 27, 2016
@dpastoor
Copy link
Contributor

I "kind of" see - pardon my naivety - if it's not too much trouble can you expand on what the difference in advancing the system between those 2.

** I had typed out a longer response, and have now deleted it since I think something clicked in my head, hopefully my understanding below is correct **

are you saying, that if this code was used

code <- '
$PARAM CL = 1, V = 35, KA = 1
$CMT GUT CENT
+$SUBROUTINES advan=2, trans=2

$MAIN 
double CLi = CL*exp(ETA(1));
double Vi =  V*exp(ETA(2));
double KAi = KA;

$OMEGA 1 1
labels=s(ECL,EV)
+$CAPTURE ECL EV CLi Vi KAi
'

the conc values, which are ultimately derived from pred_CL, etc. would be referenced to CL still, so our CLi derivation would draw from an eta distribution and be output, but that would just be seen as another user defined variable, nothing 'special' so wouldn't actually impact underlying concentration calculations, which would still be driven by the pop CL/V since they are mapped to pred_CL = CL etc.

If this is the case, then it of course makes absolute sense why you need that explicit distinction.

@kylebaron
Copy link
Collaborator Author

kylebaron commented Apr 27, 2016

@dpastoor 's explanation matches up with what I was thinking.

@kylebaron
Copy link
Collaborator Author

Hopefully the specification has been simplified. But I can still see confusion ahead.

I'm about to check in some code that searches $MAIN for tokens matching the trans (even if the user didn't specify one. Informative error if something not right.

library(dplyr)
library(magrittr)
library(mrgsolve)

PROBLEM

code <- '
$CMT CENT
$PKMODEL ncmt=1
$PARAM V1 = 2

$MAIN
double CL = 1;
double VC = 2;

'


mod <- try(mcode("foo1", code))

mod
## [1] "Error : Required PK parameters not found: V\n"
## attr(,"class")
## [1] "try-error"
## attr(,"condition")
## <simpleError: Required PK parameters not found: V>

FINE

code <- '
$CMT CENT
$PKMODEL ncmt=1
$PARAM V = 2

$MAIN
double CL = 1;
double VC = 2;

'

mod <- try(mcode("foo2", code))
## Compiling foo2.cpp.cpp ...

## done.

## Loading: foo21cda2e16405b.so
mod
## 
## 
## -------- mrgsolve model object (unix) --------
##   Project: /private/var/fol...dbsj0x4637ptj4hnd7kg900000gp/T/RtmpdSgpKr
##   source:        foo2.cpp
##   shared object: 1cda2e16405b (loaded)
## 
##   compile date:  04/28 16:28
##   Time:          start: 0 end: 24 delta: 1
##   >              add: <none>
##   >              tscale: 1
## 
##   Compartments:  CENT [1]
##   Parameters:    V [1]
##   Omega:         0x0 
##   Sigma:         0x0 
## 
##   Solver:        atol: 1e-08 rtol: 1e-08
##   >              maxsteps: 2000 hmin: 0 hmax: 0

PROBLEM

code <- '
$CMT CENT PERIPH
$PKMODEL ncmt=2
$PARAM V = 2

$MAIN
double CL = 1;
double VC = 2;

'

mod <- try(mcode("foo3", code))
mod
## [1] "Error : Required PK parameters not found: V2,Q,V3\n"
## attr(,"class")
## [1] "try-error"
## attr(,"condition")
## <simpleError: Required PK parameters not found: V2,Q,V3>

FINE

code <- '
$CMT CENT PERIPH
$PKMODEL ncmt=2

$PARAM V = 2, KM = 2, V3 = 2

$GLOBAL
bool foo = false;

$MAIN 
double CL = 1;
double V2 = 2;
 double Q = 3;
double KA = 1;
'


mod <- try(mcode("foo4", code))
## Compiling foo4.cpp.cpp ...

## done.

## Loading: foo41cda4e2517ef.so
mod
## 
## 
## -------- mrgsolve model object (unix) --------
##   Project: /private/var/fol...dbsj0x4637ptj4hnd7kg900000gp/T/RtmpdSgpKr
##   source:        foo4.cpp
##   shared object: 1cda4e2517ef (loaded)
## 
##   compile date:  04/28 16:28
##   Time:          start: 0 end: 24 delta: 1
##   >              add: <none>
##   >              tscale: 1
## 
##   Compartments:  CENT PERIPH [2]
##   Parameters:    V KM V3 [3]
##   Omega:         0x0 
##   Sigma:         0x0 
## 
##   Solver:        atol: 1e-08 rtol: 1e-08
##   >              maxsteps: 2000 hmin: 0 hmax: 0

PROBLEM

code <- '
$CMT GUT CENT PERIPH
$PKMODEL ncmt=2, depot=TRUE, trans=11

$PARAM V = 2, KM = 2, V3 = 2, Qi = 2

$GLOBAL
bool foo = false;

$MAIN 
double CL = 1;
double V2 = 2;
double Q = 3;
double KA = 1;
'

mod <- try(mcode("foo5",code))

mod
## [1] "Error : Required PK parameters not found: CLi,V2i,V3i,KAi\n"
## attr(,"class")
## [1] "try-error"
## attr(,"condition")
## <simpleError: Required PK parameters not found: CLi,V2i,V3i,KAi>

FINE

code <- '
$CMT GUT CENT PERIPH
$PKMODEL ncmt=2, depot=TRUE, trans=11

$PARAM V = 2, KM = 2, V3 = 2, Qi = 2

$GLOBAL
bool foo = false;

$MAIN 
double CLi = 1;
double V2i = 2;
double Q = 3;
double V3i = 200;
double KAi = 1;
'

mod <- try(mcode("foo6",code))
## Compiling foo6.cpp.cpp ...

## done.

## Loading: foo61cda450777ef.so
mod
## 
## 
## -------- mrgsolve model object (unix) --------
##   Project: /private/var/fol...dbsj0x4637ptj4hnd7kg900000gp/T/RtmpdSgpKr
##   source:        foo6.cpp
##   shared object: 1cda450777ef (loaded)
## 
##   compile date:  04/28 16:28
##   Time:          start: 0 end: 24 delta: 1
##   >              add: <none>
##   >              tscale: 1
## 
##   Compartments:  GUT CENT PERIPH [3]
##   Parameters:    V KM V3 Qi [4]
##   Omega:         0x0 
##   Sigma:         0x0 
## 
##   Solver:        atol: 1e-08 rtol: 1e-08
##   >              maxsteps: 2000 hmin: 0 hmax: 0

@kylebaron
Copy link
Collaborator Author

@dpastoor If you're using this ... in addition to any bugs, can you report any awkwardness in the use or options? The latest commit 9c5130a fills in some of the details via ?PKMODEL.

@dpastoor
Copy link
Contributor

dpastoor commented May 1, 2016

Yep would have already brought it up if it felt too weird. So far been very pleasant (beyond that one tricky error)

@kylebaron
Copy link
Collaborator Author

Third or fourth iteration on this ... I like the way this turned out too. Thanks @vjd and @dpastoor .

@kylebaron
Copy link
Collaborator Author

This seems pretty stable so far. Closing.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

3 participants