How to capture time of dosing events

[sbihorel]

Hi,

To implement an approximated transit compartment absorption model for repeating dosing study design, I would need to capture the time of dosing events (which might vary from simulation/id to simulation/id) in an array then make the array available in [ ODE ].

Can this be implemented in mrgsolve?

[kylebaron]

Here is a single dose; I just captured the value for demonstration purposes; you can use as you wish in ODE

dose_time.cpp:

[ MAIN ] 

if(NEWIND <=1) double told = 0;

if(EVID==1) {
  told = TIME;
}

[ ODE ]
// USE told and SOLVERTIME in this 
dxdt_COMPARTMENT = -0.1*COMPARTMENT;

[ CMT ] COMPARTMENT

[ CAPTURE ] told

library(mrgsolve)

mod <- mread("dose_time")
#> Building dose_time ... done.

mod %>% 
  ev(amt = 100, ii = 6, addl = 4) %>% 
  mrgsim(end = 72) %>% 
  plot()

^{Created on 2019-02-03 by the reprex package (v0.2.1)}

[kylebaron]

We can do something more complex too (track several doses etc); Just let me know the specifics and I can give you the example code.

[sbihorel]

Thanks Kyle

I did not reply right away to your post because I want to take the information you gave and try to find the solution by myself... I think I got it.

The following is a pseudo direct translation of a NONMEM code described by Alison Boeckmann (NONMEM examples/simdosetn.ctl). It gives virtually identical prediction in mrgsolve, NONMEM, and Berkeley Madonna for a deterministic simulation in 1 subject.

Do you recommend any change to way it is coded or anticipate any issue when used in a stochastic multi-subject scenario?

[ GLOBAL ]
  #define CP (A2/VC)
  double dosetime [99];
  double dose [99];

[ CMT ] @annotated
  A1 : DEPOT
  A2 : CENTRAL

[ PARAM ] @annotated
  CL  : 3    : Elimination Clearance (-)
  VC  : 10   : Central Volume (-)
  MTT : 3    : Mean Transit Time (-)
  NN  : 2.5  : Number of Transit Compartment (-)

[ MAIN ]

  // Capture dose information
  if (NEWIND <= 1) {
    int ndose = -1;
  }
  if (EVID == 1){
    ndose = ndose + 1;
    dosetime[ndose] = TIME;
    dose[ndose] = self.amt;
  }

  // Derived parameters
  double K20 = CL/VC;
  double KTR = NN/MTT;
  double KINPT =  pow(KTR,NN)/exp(lgamma(NN));

  // Bioavailability
  F_A1 = 0;

[ ODE ]
  double INPT = 0;
  double IPT = 0;
  int i = 0;
  while (i <= ndose) {
    if (SOLVERTIME >= dosetime[i]){
      IPT = dose[i]*pow(SOLVERTIME - dosetime[i], NN - 1)*exp(-KTR*(SOLVERTIME - dosetime[i]));
    }
    INPT = INPT + IPT;
    i = i + 1;
  }

  dxdt_A1 = KINPT*IPT - KTR*A1;
  dxdt_A2 = KTR*A1 - K20*A2;

[ CAPTURE ]
  CP

[sbihorel]

I spoke too fast. The dose time capture does not work in a multiple dose settings.

require(mrgsolve)

mod <- mread('mrg-1cmt-transit.cpp')

out <- mod %>%
  ev(amt = 1000, addl=3, ii = 24, cmt = 1) %>%
##  ev(amt = 1000, time = 24, addl=3, ii = 24, cmt = 1) %>%
  mrgsim(
    end = 48*2,
    delta = 0.01)

print(plot(out))

[kylebaron]

I can take a look with you.

[kylebaron]

Your code is right; the problem is related to this: https://github.com/metrumresearchgroup/mrgsolve/issues/368

It should have gotten resolved in the last release. I will post a work-around.

[kylebaron]

This works:

library(mrgsolve)

mod <- mread('mrg-1cmt-transit.cpp')

e <- ev(amt = 1000, addl=3, ii = 24, cmt = 1) %>% realize_addl()

out <- mod %>% mrgsim_e(e, end = 48*2, delta = 0.01)

plot(out)

[sbihorel]

Awesome ! Thanks !

Out of curiosity, what is your release cycle on CRAN?

[kylebaron]

We have an internal change control process that needs to happen and we can't release to CRAN until then. So it's pretty easy to push this stuff to CRAN but more work (people, paperwork) to do the internal stuff.

I'll see about doing something at the end of Feb (another undocumented feature that I wanted to tweak a bit too).

I'll get #368 addressed and at least available in the dev branch.

Thanks for your patience.

BTW: this is another workaround that is functionally the same as what you originally proposed

• mrgsolve doesn't mind of you have dose amount zero on evid 1 record
• so just do zero dose and include Amt as the dose amount you want
• be sure to list Amt as a parameter.
• and then don't set bioavailability to zero

library(mrgsolve)
library(dplyr)

mod <- mread('mrg-1cmt-transit-2.cpp')

e <- ev(amt = 0, Amt = 1000, addl=3, ii = 24, cmt = 1) 

out <- mod %>% mrgsim_e(e, end = 48*2, delta = 0.01)

plot(out)

mrg-1cmt-transit-2.cpp

[ GLOBAL ]
#define CP (A2/VC)
double dosetime [99];
double dose [99];

[ CMT ] @annotated
A1 : DEPOT
A2 : CENTRAL

[ PARAM ] @annotated
CL  : 3    : Elimination Clearance (-)
VC  : 10   : Central Volume (-)
MTT : 3    : Mean Transit Time (-)
NN  : 2.5  : Number of Transit Compartment (-)
Amt : 0    : Dose amount

[ MAIN ]

// Capture dose information
if (NEWIND <= 1) {
  int ndose = -1;
}

if(EVID == 1){
  ndose = ndose + 1;
  dosetime[ndose] = TIME;
  dose[ndose] = Amt;
}

// Derived parameters
double K20 = CL/VC;
double KTR = NN/MTT;
double KINPT =  pow(KTR,NN)/exp(lgamma(NN));

[ ODE ]
double INPT = 0;
double IPT = 0;
int i = 0;
while (i <= ndose) {
  if (SOLVERTIME >= dosetime[i]){
    IPT = dose[i]*pow(SOLVERTIME - dosetime[i], NN - 1)*exp(-KTR*(SOLVERTIME - dosetime[i]));
  }
  INPT = INPT + IPT;
  i++;
}

dxdt_A1 = KINPT*IPT - KTR*A1;
dxdt_A2 = KTR*A1 - K20*A2;

[ CAPTURE ] CP

[sbihorel]

No rush. I am still at the stage of playing with mrgsolve, so I can wait.

[kylebaron]

Ok; good. I'd like to capture this example here too. I could do it with attribution or I'd be happy to accept a pull request too.

[sbihorel]

No problem. I am not yet up to speed on the whole pull request thing... github is still a bit of a mistery to me (I come from SVN). Attribution is fine (Alison Boeckmann is the main source, IMO).

[kylebaron]

Thanks. Closing this as the fundamental issue was addressed in https://github.com/metrumresearchgroup/mrgsolve/commit/e697edd4aadf129132e853d2b573c4217d7dd1e1 . See #368