RXR-2378 bug with Fi defined in ode producing NA concentrations

[roninsightrx]

A model where F (bioav) was defined in ODE block as:

    F1_avg = F1 * pow(prv_dose/1000.0, -0.15) ;\
    F1i = F1_avg * exp(kappa_F1) ;\

resulted in NAs in the simulated concentrations. The code however runs fine when used in PK code.

This was because in the internal C++ code, the first row is not actually a dosing event, but an initalization row that does not have a dose. However, prv_dose is initialized as: prv_dose = doses[0];, which then results in prv_dose being zero. Then in the above definition a value of prv_dose=0 will results in Inf for F1_avg.

We can avoid this by instead just reading the first non-zero value in the dose column of the event table and initializing prv_dose with that.

This change shouldn't affect anything else.

[roninsightrx]

Unrelated. This line was actually not used, the grep function was looking for // insert bioav definition, so never picked up on this one. A little down in the code it was written correctly. So it seems this line was not functional, therefore can be removed.

[jasmineirx]

thanks for looking into this. i didn't realize this was the cause of the issue. a couple thoughts on if it makes sense to code around this or go with the fix.

[jasmineirx]

Suggested change

	prv_dose = doses[idx];
	t_prv_dose = times[0];
	prv_dose = doses[idx];
	t_prv_dose = times[idx];

[jasmineirx]

with this change, we are building in a rather quiet assumption that this variable is only correct if doses start at t = 0.

maybe for now this is okay and we can rewrite the model code to something like

F1_avg = F1 \
if (t > 0) { F1_avg = F1 * pow(prv_dose/1000.0, -0.15) } \

[jasmineirx]

On the other hand, other models also could use prv_dose in the ode_code block, and it is maybe also unexpected for it to be 0 momentarily.

[roninsightrx]

yes, my thinking was that if we would keep PKPDsim as-is, and have to make some kind of modification in each model, that would be more annoying, and we might forget for future models. I was thinking this modification, but it's similar to yours.

FACT = 1
if(prv_dose > 0) FACT = pow(prv_dose/1000.0, -0.15)
F1_avg = F1 * FACT

Neither solution is perfect. Let me look into one other solution I was thinking about though, i.e. we could also make sure that doses[0] is not zero to begin with (at least as long as at t=0 there is a dose).

[roninsightrx]

Update: there is actually another problem that this PR should address. This is that t_prv_dose (the time of previous dose, accessible in ODE block) is currently calculated at the time of the dose + lag time. This is unintuitive, and different from how one would implement it in NONMEM.
TODO: the changes I just put up fixes this behavior somewhat, but is still not entirely correct.

[roninsightrx]

Another issue that investigations into this PRs revealed is that lag-times, when specified as an estimated parameter, stay fixed because the design going into the core simulation function does not dynamically update. I.e. it is impossible to estimate individual estimates for lag time at the moment. It's not really a problem for our models (we only have a handful of models for which a lag time is estimated, and even for those we rarely have information in the absorption phase), but we should fix.