Discrepancy in cmax estimation under default settings

[JohannesStarp]

Hi MRGsolve team,

I’ve been working on replicating the workflow outlined in your recent blog post "Simulate true cmax with ode" (https://mrgsolve.org/blog/posts/true-cmax.html), and I noticed some unexpected behavior related to cmax estimation.

Summary

In line with the approach described in the blog post, I implemented code to track PK parameters such as cmax and tmax directly in the $DES block. However, when running simulations with the default solver settings, I observed that the predicted cmax values sometimes precede the actual peak in the concentration-time curve, especially during steep absorption phases.

When I limited the solver step size using the hmax (hmax = 0.001) option, the predicted cmax aligned correctly with the true maximum concentration. This suggests the issue may be tied to how the ODE solver handles internal step sizes versus output time points.

Hypothesis

It seems like the ODE solver uses lager step sizes in order to then calculate the actual step sizes. During this proces higher concentrations are used in the $DES block, which are then saved as Cmax. These values safed for Cmax are not readjusted when a finer step size is used to calculate the concentrations in the central compartment.

Reproducible Example

To better understand and document this issue, I created a GitHub repository with a Quarto report that reproduces the findings using both one- and two-compartment models and different dosing routes:

📁 Repository: https://github.com/JohannesStarp/mrgsolve_pk_parameters
📄 Quarto Document (in repo): Quarto pk parameters mrgsolve.qmd

The report includes side-by-side plots showing the impact of solver settings on cmax predictions.

Feedback

I’d appreciate your thoughts on whether this behavior is expected under the current implementation. Is there a more robust way to capture cmax and tmax values internally without relying on modifying solver step sizes (which, as noted in the docs, is discouraged for regular use)?

Thanks for your time and for continuing to develop such a powerful simulation tool!

Best regards,
Johannes

[kylebaron]

Hi @JohannesStarp -

Thanks for reaching out and for interacting with the Blog post. I think you're exactly right: using this ODE method to calculate Cmax probably needs the stricter settings to prevent this behavior that is nicely demonstrated by the quarto content in your GitHub. The blog post focused on the infusion, but the ODE solver used in mrgsolve will always show this behavior. As we noted in the Blog post, this is feature of the ODE solver ... to help it go faster. But we'll see this outcome when trying to find Cmax unless we take another approach.

Ruben is right in this post here: the right way to do this is with a root finder. I think that will properly handle the numerical considerations here.

In the majority of my day to day work, I don't use either of these approaches. I tend to do a very simple approximation of Cmax by just simulating out a "fine" grid of points and picking off the maximum concentration. It's definitely an approximation but IMO it tends to work pretty well and it can handle a lot of different "odd" scenarios. It's not satisfying for people who love exact solutions to questions, but it's a practical solution, IMO. Especially when you're simulating Cmax and limiting the number of significant digits when you report the result. Obviously, you can easily come up with a model where this might not work well ... or where you need a lot of points to get close to Cmax. In my experience, it doesn't happen that often. And when absorption happens over an extended period of time, you will be hard-pressed to tell the difference in the Cmax methods. When the drug is infused, you don't have to do any finding ... you know when Tmax is and you know when to look for Cmax. Obviously, when this approximation method is used, I'd report exactly what you did so the reader understands how the numbers were generated.

Anyhow - thanks for your message; the Quarto stuff is really great and I very much appreciate the fact that you are digging into this and helping us understand the software and the problem better.

Best regards,
Kyle

PS: Ruben's stuff shows I think that the "fine grid" approach isn't the fastest computationally. The root finder is faster. For me, part of the practical speed of "find grid" is not computational efficiency, but faster time to final result by minimizing set up time etc.

PSS: this was one of my most controversial posts ever on LinkedIn. I'm sure there will continue to be disagreements and I'm happy to keep the discussion going. I'm als happy to acknowledge the shortcomings of my approaches and the strengths of other more sophisticated approaches.

[JohannesStarp]

Hi @kylebaron,

Thank you so much for your detailed and thoughtful response - really appreciated!

Your explanation and the link to your post helped me a lot in understanding not just the behavior in MRGsolve, but also how this issue is rooted in the ODE solver's design, and how it appears in other tools like NONMEM as well (https://www.mail-archive.com/nmusers@globomaxnm.com/msg07043.html).

For my upcoming projects, I’ll most likely follow the "fine grid" strategy to approximate PK/PD parameters like cmax and tmax, and be transparent in reporting how they were derived. This is transperent and makes it easier to spot potential abnormalities in the output file directly.

Thanks again for your time, the blog and your respond!

Best regards,
Johannes

[kylebaron]

Thanks for following up, @JohannesStarp. Just connected with you on LinkedIn. Not sure how much longer you have in your PhD program, but this might be an interesting question to explore. The fine grid approach is an approximate Cmax and it will be lower than the "true" Cmax most of the time (it will never be higher than the true). So there should be some bias in there. But how much? You could explore this through data simulated from models with different absorption profiles. Maybe this is just a generic model (like 2 compartment with first-order absorption or zero-order or sequential zero-first or transit etc) with different parameter values that might affect how flat (or not flat) the profile is around Cmax. Assess with the root finder or with the "embedded" Cmax code or fine grid. If embedded Cmax code, when do you need start worrying about ODE solver settings? If fine grid, what should be the grid? can you predict a sufficient grid? what is the bias? I think if you could answer some of these questions, it would be very interesting to our community.

[kylebaron]

If any of this is interesting to you and something you might want to pursue, I'd be glad to collaborate. If not ... I totally understand!

[JohannesStarp]

Thanks again @kylebaron for the thoughtful follow-up and for connecting on LinkedIn - I really appreciate it!

You're absolutely right, the idea of systematically exploring the bias introduced by the fine grid approach across different absorption profiles is a very interesting research question!

At the moment, though, I’m still fully occupied with other projects that are more practically oriented and not focused on simulation-based methodology work. So unfortunately, I won’t be able to dive into this topic in the near term. That said, I’ll definitely keep this idea in mind for future work - it’s a great example of how understanding what the program is doing is essential for interpreting what it outputs.

I really appreciate your offer to collaborate, and I’d be very happy to come back to you if and when I start working on this topic in more depth. Thanks again for the engaging discussion and the generous support!

Best regards,
Johannes