Torsten implementation of IV

General
1

I don’t know is this is the right place but while trying to implement multi dose IV I happened to look at torsten codes and was confused by IV rate implementation. It applies both to 1cmp and 2cmp models.

In the Pred1_one function there is a statement which updates concentration in central compartment:
if ((init[1] != 0) || (rate[1] != 0)) {
a[0] = 1;
pred(0, 1) += PolyExp(x=dt, dose=init[1], rate=0, xinf=0, tau=0, ss=false, a, alpha, n=1) +
PolyExp(dt, 0, rate[1], dt, 0, false, a, alpha, 1);
}

which translates into:
bolusResult = init[1] * exp (-k10 * dt);
result = (1 - exp (-k10 * dt)) / k10;
pred(0, 1) = bolusResult + rate[1] * result;

However units of rate in mg/h and concentration in mg/L become inconsistent. However, if we replace k10 with CL (L/h) then units become consistent.

Please pardon for my ignorance.

2

@charlesm93 is the one to ask. I don’t know if Metrum set up a mailing list or forum for Torsten.

3

Hi Linas,

When amt (and init) has units of mass and rate has units of mass/time, the pred functions of Torsten all calculate mass, not concentration. The user is responsible for converting mass to concentration in the Stan model code, Thus the units are consistent.

Bill

4

Regarding a mail list for Torsten, there is none yet. My guess is there’ll be one once we expect higher traffic, and v1.0 gets released. I’d say it’s appropriate to post issues on github, github.com/metrumresearchgroup/stan. If you post on the Stan forum, we’ll respond + other Stan users and developers can step in, which is if course a plus.

5

I assume I can divide the output by Vd then?

6

Yes. You may want to look at examples at https://github.com/metrumresearchgroup/example-models to see some use cases for the Torsten functions.

7

Thanks. The next item I am struggling with is the parameters that change
over the time. The signature of one cpt model requires that theta (the
parameter next to biovar) has to be real[] or real[,]. While the attached
code runs for real[] case it doesn’t for real[,] case (commented code). Any
idea where is the problem?

onecmpaddpopt.stan (4.2 KB)

8

Hi Linas,

Sorry for the delayed reply.

What error messages do you get? Does it indicate a function signature problem or something else?

Thanks,
Bill

9

Hi Linas,

One possible source of the problem is the following loop:

  for (i in 1:N_rest) {
  	theta[i, 1] = thetaM[j, 1] * (WT_data_rest[i] / NORM_WT) ^ 0.75 * 
  		thetaM[j, 3] * SCR_data_rest[i] / NORM_SCR; // CL
  	theta[i, 2] = thetaM[j, 2] * (WT_data_rest[i] / NORM_WT); // V2
  	theta[i, 3] = 0; // ka
  }

I think your intention is to only populate theta with parameters for individual j, but the loop index corresponds to the entire data set since N_rest appears to be the total number of event times for all individuals. Maybe you want something like “for(i in 1:nt[j])” for the loop index where nt[j] is the number of event times for the jth individual, i.e.,

for(i in 1:nSubjects)
ntj[i] = end_rest[j] - start_rest[j] + 1;

Since you’re probably dealing with a ragged array, you will also need to specify the same index range when you pass theta as an argument, i.e.,

  x[start_rest[j]:end_rest[j]] = PKModelOneCpt(
  	time_rest[start_rest[j]:end_rest[j]], 
  	AMT_data_rest[start_rest[j]:end_rest[j]],
  	RATE_data_rest[start_rest[j]:end_rest[j]],
  	ii_rest[start_rest[j]:end_rest[j]],
  	evid_rest[start_rest[j]:end_rest[j]],
  	cmt_rest[start_rest[j]:end_rest[j]],
  	addl_rest[start_rest[j]:end_rest[j]],
  	ss_rest[start_rest[j]:end_rest[j]],
  	theta[1:nt[j], ], biovar, tlag);

Finally declare theta with a size corresponding to the largest number of event times for a single individual, i.e.,
real<lower = 0> theta[max(ntj), 3];

Cheers,
Bill

10

Bill,

Thanks. I have corrected this. There is still relatively small mismatch
between torsten and my implementation of one compartmental model. Both are
attached in case you are interested. Torsten model gives much more
divergencies even if I use Matt’s trick and sd of estimates is much larger.

Linas
Linas

onecmpaddpop.stan (4.5 KB)

onecmpaddpopt.stan (5.58 KB)

11

Hi Linas,

Thanks for these. Are you able to share the data? Otherwise it’s hard for me to diagnose the discrepancy. The fact the sd estimates are much larger makes me suspect the PK model calculations are different event given the same parameter values. Otherwise I would expect them to converge to the same posterior.

Thanks,
Bill

12

Bill,

One possibility for mismatch is the way IV is handled. In my
implementation whenever IV occurs I add two dosing events: one with real
dose and real rate and another 1h later with 0 dose and 0 rate. Since it
is possible that some labs will occur within this 1h window I sort based
on time. I don’t know how this is handled in torsten.

Linas

13

One clarification. For my model I use two events - start and stop of IV.
For torsten model I drop stop IV event. Maybe I shouldn’t do it?

14

If you can’t share the data set, would you share an excerpt of the data showing how the dosing is specified? I would really like to understand if there is a bug I need to fix or just a misunderstanding about how to specify dosing.

Torsten uses NONMEM conventions for dose specifications. If I understand correctly, you are trying to specify a constant rate input with a 1 hour duration. Is that correct? If so, that is specified with a single dosing event record in which rate = the input rate, amt = rate * (input duration – 1 h in this case), evid = 1, cmt = compartment receiving input (2 for PKModelOneCpt if it is IV input).