Re: [stan-dev/math] Feature: Torsten (pharmacometrics function library) (#314)

[Bob Carpenter]

[off issue]

We want to use issues only when the functionality is
well defined enough that unit tests can be written. We
have wikis and mailing list for general design discussions
before that.

We would also rather have issues broken down into bite-sized
pieces so that they're easy to code review for test coverage.

Software development should not be done by coding and then
adding unit tests and design as an afterthought. The API,
doc, and unit tests should all be designed together. A
good API design means something that's easy to test and doc.

The code should be designed top down and implemented
bottom up with unit tests on all the pieces written in tandem
with the code itself.

So it sounds like issue #314 is premature if you don't know
what the unit tests would look like.

- Bob

> On Aug 3, 2016, at 5:09 PM, Charles Margossian <notifi...@github.com> wrote:
>
> @syclik, no action required. I created this issue to be consistent with Stan practices when implementing a new feature. My understanding is that this precedes submitting a pull request. Btw, I could not create a new branch directly on stan-dev, so I first created a fork -- I don't believe this should pose any problem.
>
> I anticipate we will be revising the design for Torsten, and pharma applications in general, during the upcoming meeting with the Stan group and Sebastian Weber. So I'm holding off designing unit tests until that's settled.
>
> —
> You are receiving this because you are subscribed to this thread.
> Reply to this email directly, view it on GitHub, or mute the thread.
>

[Charles Margossian]

Hi Bob,

I apologize for creating an issue prematurely. My goal was to share the code on GitHub and I attempted to use the procedure specified on https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan. To do so, I created a branch, which I understood had to be named after an issue. 

My intent was to share a working prototype code to illustrate an approach to Stan functions for pharmacometrics modeling, and to clarify by example the desired functionality. How would you prefer such code be shared with the Stan community? 

Best,

Charles  

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Daniel Lee]

Hi Charles,

On Thu, Aug 4, 2016 at 9:48 AM, Charles Margossian <charles.ma...@gmail.com> wrote:

Hi Bob,

I apologize for creating an issue prematurely. My goal was to share the code on GitHub and I attempted to use the procedure specified on https://github.com/stan-dev/stan/wiki/Contributing-New-Functions-to-Stan. To do so, I created a branch, which I understood had to be named after an issue. 

That documentation was for developers that knew what function they were contributing to Stan. I'm not sure your contributions fall under that. More below.

 

My intent was to share a working prototype code to illustrate an approach to Stan functions for pharmacometrics modeling, and to clarify by example the desired functionality. How would you prefer such code be shared with the Stan community? 

I'm still missing a big piece of what you're doing. In particular, I don't understand how much of this is additions to:

- stan-dev/math

- stan-dev/stan

- stan-dev/cmdstan

It's hard to recommend how to share code without understanding what parts are changing. One option is to fork CmdStan and everything from there down as a prototype.

If I recall, I thought you were going to create a wiki page with a description and a design. Could you put that up somewhere?

[Charles Margossian]

Hi Daniel,

Yes, I indeed created  a wiki on Strategies to Facilitate Stan in Pharmacometrics, with a section on Torsten. I updated it so that it links to the example models I posted on a fork of example-model. 

I also wrote a doc on how to use the examples. Today I added a third section, "Under the Hood  Design", which I hope will answer your questions. 

I'm still missing a big piece of what you're doing. In particular, I don't understand how much of this is additions to:
- stan-dev/math
- stan-dev/stan
- stan-dev/cmdstan

At the end of the day, Torsten is four built-in Stan functions. They are grouped together because they have a similar architecture and functionality. Additions are made to stan-dev/math (see here) and stan-dev/stan (see here). stan-dev/cmdstan remains unchanged.  

I hope this sheds some light!

Best,

Charles 

On Thu, Aug 4, 2016 at 12:37 PM, Daniel Lee <bea...@alum.mit.edu> wrote:

Hi Charles,

 

Best,

Charles  

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Krzysztof Sakrejda]

Hey,

I read through (much of) your code and examples and I'm trying to understand the difference between the standard Stan ODE solver and your compartmental model functions. I had a few thoughts:

- Seems like PKModelTwoCpt (and the OneCpt version) are convenience functions in the sense that they wrap up those specific ODE systems for easy use. They are not really pharma-specific in the sense that other modeling could use them (e.g.-geochemical systems). Someday it might be nice to have a library of pre-implemented ODE's for a variety of model types like this and starting with pharma-related examples might be a nice thing but for Stan itself it might make more sense to cast them more broadly

- The GeneralCptModel_* functions seem like they are just adapters to use ODE systems with the NONMEM data format---is that right or do they serve another function?

- Then there are a lot of helper functions that make it really slick to use the NONMEM data format (processing multiple events and restarts of the solver associated with the events). It seems like having
helper functions to turn NONMEM data into a more ODE-compatible format could live separately from the code required to deal with multiple solver restarts associate with events. That might make it easier for both pieces to get adopted into Stan (and/or the interfaces).

So if I'm mostly correct here, my advice would be to break these things out more if they are going to stand a chance of getting integrated into Stan (even as libraries). Maybe something like:

1) Some constructs to make it more convenient to use the Stan ODE solvers with multiple restarts (associated with events in your terminology).

2) Some functions to handle converting NONMEM data into whatever more generic format is required by 1)

3) A library of ODE's that implement pharma-specific models and accept data in the NONMEM format.

As is all three seem to be packaged together and that might be hard for the development process to swallow.

Krzysztof

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

>
> For more options, visit https://groups.google.com/d/optout.
>
>
>
>
>
>
>
> --
>
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
>

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

>
> For more options, visit https://groups.google.com/d/optout.
>
>
>
>
>
>
>
>
>
> --
>
> You received this message because you are subscribed to the Google Groups "stan development mailing list" group.
>

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+u...@googlegroups.com.

[Sebastian Weber]

Hi!

Great to have this thread! As I have already done projects with Stan in this field of applications, I wrote my own version of Torsten. However, as we "only" have vanilla Stan in our production I used only Stan functionality. I shared a condensed set of this here

https://github.com/stan-dev/example-models/tree/feature/issue-72-stan-pkpdlib/misc/pkpd

I think it needs to be balanced as to what approach we want in going forward. My idea behind the above code is to hide all the boiler plate code for the index fiddling (and ODE restarting) from the user by using the include functionality of rstan. I could have made more use of it. Retrospectively, the constructs which would have made my life easier to program this is

- ragged array structures
- tuples (and ragged arrays of tuples)
- looping constructs over ragged arrays
- a more convenient way to have data variables in my functions (it is a pain to pass from the transformed data block the "x_r" ODE argument through all functions)

Other than that my little library supports analytic and ODE based systems. Features missing in comparison to Torsten are (which I am aware of)

- steady-state conditions (although I have regular dosing); but I am not sure as to how this is implemented in Torsten as it needs an implicit equation solver which we do not have yet in Stan - sundials has such a solver with sensitivity support

- infusions; this can be added, just more index index fiddling and often there are ways to avoid this

I did want to add to the wiki's Charles has started. My apologies for not yet doing so, instead I have written an extensive README which is placed in the branch I point to above.

Best,
Sebastian

[Charles Margossian]

Hi Krzysztoph, 

Thank you for taking the time to read through the code and for the very constructive comments. 

I agree with you that Torsten should be broken up, especially since many of its features appear to be extendable to fields other than pharmacometrics. The initial idea was to take a set of functions present in NONMEM and BUGS, and see if we could implement them in Stan. Fitting everything into one function was meant to meet a desired and expected simplicity for the user. I believe one of Torsten's strengths is how accessible and convenient it is for pharmacometricians. 

 Evidently, there seems to be a more ambitious and better approach. 

They are not really pharma-specific in the sense that other modeling could use them (e.g.-geochemical systems). 

I'm trying to learn about ODE based models in other fields. Right now, my focus in on system biology, but if you have any papers or documents on geochemistry, I am very interested. My hope is to get a better sense of what is pharma specific and what is generalizable in Torsten. 

Someday it might be nice to have a library of pre-implemented ODE's for a variety of model types like this

Agreed. It might even be possible to code the gradients by hand for certain models, and I want to try it out with the one and two compartment models. 

The GeneralCptModel_* functions seem like they are just adapters to use ODE systems with the NONMEM data format---is that right or do they serve another function? 

I think you have the right idea.

It seems like having helper functions to turn NONMEM data into a more ODE-compatible format could live separately from the code required to deal with multiple solver restarts associate with events. 

An "internal format", to use Sebastian's term, would indeed be the way to go to make the code more general, and is what has been implemented in Sebastian's code. We should keep in mind that some components of the event schedule could also be parameters. So we would need a wrapper at both ends of the process, one that converts NONEM into the internal format, and one that converts the result back into NONMEM format. 

So if I'm mostly correct here, my advice would be to break these things out more if they are going to stand a chance of getting integrated into Stan (even as libraries).  Maybe something like: 
1) Some constructs to make it more convenient to use the Stan ODE solvers with multiple restarts (associated with events in your terminology). 
2) Some functions to handle converting NONMEM data into whatever more generic format is required by 1) 
3) A library of ODE's that implement pharma-specific models and accept data in the NONMEM format. 

This strikes me as a very sensible approach.

(1) would be the "Event Handler" -- currently the pred function in Torsten -- which augments the event schedule and sequentially applies to each event the "solution operator" (to use another very good term Sebastian coins in his ReadMe). It might be a good idea to split pred into two functions that individually handle (i) augmenting the data and (ii) applying the solution operator; mostly because (i) often only involves data and could thus happen in the transformed data block for efficiency purposes. 

I would like to find an easy way for the user to specify the solution operator, by letting him: 

(a) write the operator in the function block of Stan

(b) write out an ODE system and pick an ODE integrator

(c) call a built-in analytical operator -- currently pred1_one, pred1_two, predSS_one, and predSS_two in Torsten. I am not sure whether these should work with NONMEM data or the internal format in order for them to work with the Event Handler. 

Currently, Torsten handles (b) and (c), but it constrains the user to pass in NONMEM data. 

As is all three seem to be packaged together and that might be hard for the development process to swallow.

I would agree with you, though in the end, I think it comes down to only testing 4 functions (currently finishing up unit tests for OneCptModel). Though admittedly, I am not sure on that one. 

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Charles Margossian]

Hi Sebastian,

Great to have this thread!

Agreed. This discussion should make next week's meeting significantly more productive -- Krzysztof, I cordially invite you to participate (Friday the 12th @ 10 am). Thank you for posting your code again. I very much enjoyed the README you wrote, and I've started looking at the examples you shared.   

Other than that my little library supports analytic and ODE based systems. Features missing in comparison to Torsten are (which I am aware of) 
- steady-state conditions (although I have regular dosing); but I am not sure as to how this is implemented in Torsten as it needs an implicit equation solver which we do not have yet in Stan - sundials has such a solver with sensitivity support
- infusions; this can be added, just more index index fiddling and often there are ways to avoid this  

 

I need to learn more about your Torsten, and probably Bill's expertise to comment on that. Steady-state approximations are only available for analytical one and two compartment model, but currently not for the generalCptModel_* functions.

I did want to add to the wiki's Charles has started. My apologies for not yet doing so, instead I have written an extensive README which is placed in the branch I point to above.  

I myself gave up a little bit on this wiki, because I expected we were going to create a more general one on ODE based models -- in fact it would be a good idea to do that, and include a pharma section in the new wiki. I'll set it up once I finish reading about ODE models outside pharma, but feel free to go ahead (and if you do, delete/merge the existing wiki). 

This wiki might also be a place to discuss ODE solving techniques, such as root solvers and matrix exponentials. 

Andrew, at a recent Columbia - Metrum meeting, suggested using the user mail list to get a better sense of what the community needs in terms of ODE based models, so I'll share the wiki there once it is set up. 

[Krzysztof Sakrejda]

On Friday, August 5, 2016 at 8:46:54 AM UTC-4, Charles Margossian wrote:
> I agree with you that Torsten should be broken up, especially since many of its features appear to be extendable to fields other than pharmacometrics. The initial idea was to take a set of functions present in NONMEM and BUGS, and see if we could implement them in Stan. Fitting everything into one function was meant to meet a desired and expected simplicity for the user. I believe one of Torsten's strengths is how accessible and convenient it is for pharmacometricians. 

That's ultimately a great approach to an interface for scientists. I'm only suggesting making some under-the-hood changes to make that layer thinner and less specialized.

>  Evidently, there seems to be a more ambitious and better approach. 

My suggestion is to go for something less ambitious first. The Stan development process for bringing in small well-specified chunks of functionality is pretty clear and well greased at this point. At the same time I'm not sure I've ever seen a big ambitious set of changes that went in smoothly. Maybe Bob/Daniel will just work their magic and see how to integrate all this but

> They are not really pharma-specific in the sense that other modeling could use them (e.g.-geochemical systems). 
>
> I'm trying to learn about ODE based models in other fields. Right now, my focus in on system biology, but if you have any papers or documents on geochemistry, I am very interested. My hope is to get a better sense of what is pharma specific and what is generalizable in Torsten. 

Bob has some examples in Stan----I can't find them at the moment but maybe he'll pipe up. I only did some thesis stuff by simulation as an undergrad (it's on floppies somewhere...) so I'm out of touch there. The other uses I'm familiar with are compartmental models in hydrology and SIR-type compartmental models in infectious disease but again it's not my speciality.

> Someday it might be nice to have a library of pre-implemented ODE's for a variety of model types like this
>
>
> Agreed. It might even be possible to code the gradients by hand for certain models, and I want to try it out with the one and two compartment models. 
>
> The GeneralCptModel_* functions seem like they are just adapters to use ODE systems with the NONMEM data format---is that right or do they serve another function? 
>
>
> I think you have the right idea.
>
> It seems like having helper functions to turn NONMEM data into a more ODE-compatible format could live separately from the code required to deal with multiple solver restarts associate with events. 
>
>
> An "internal format", to use Sebastian's term, would indeed be the way to go to make the code more general, and is what has been implemented in Sebastian's code. We should keep in mind that some components of the event schedule could also be parameters. So we would need a wrapper at both ends of the process, one that converts NONEM into the internal format, and one that converts the result back into NONMEM format. 
>
> So if I'm mostly correct here, my advice would be to break these things out more if they are going to stand a chance of getting integrated into Stan (even as libraries).  Maybe something like: 
> 1) Some constructs to make it more convenient to use the Stan ODE solvers with multiple restarts (associated with events in your terminology). 
> 2) Some functions to handle converting NONMEM data into whatever more generic format is required by 1) 
> 3) A library of ODE's that implement pharma-specific models and accept data in the NONMEM format. 
>
>
> This strikes me as a very sensible approach.
>
>
> (1) would be the "Event Handler" -- currently the pred function in Torsten -- which augments the event schedule and sequentially applies to each event the "solution operator" (to use another very good term Sebastian coins in his ReadMe). It might be a good idea to split pred into two functions that individually handle (i) augmenting the data and (ii) applying the solution operator; mostly because (i) often only involves data and could thus happen in the transformed data block for efficiency purposes. 
>
>
> I would like to find an easy way for the user to specify the solution operator, by letting him: 
> (a) write the operator in the function block of Stan
> (b) write out an ODE system and pick an ODE integrator
> (c) call a built-in analytical operator -- currently pred1_one, pred1_two, predSS_one, and predSS_two in Torsten. I am not sure whether these should work with NONMEM data or the internal format in order for them to work with the Event Handler. 
>
>
> Currently, Torsten handles (b) and (c), but it constrains the user to pass in NONMEM data. 

All of this makes sense and I'm looking forward to the meeting next week. I'll give Sebastien's code a better read as well.

>
> As is all three seem to be packaged together and that might be hard for the development process to swallow.
>
>
> I would agree with you, though in the end, I think it comes down to only testing 4 functions (currently finishing up unit tests for OneCptModel). Though admittedly, I am not sure on that one. 

So far Stan is tested all the way down to the math library so I'm pretty sure we'll need tests for all the pieces here.

Krzysztof

[Bill Gillespie]

Hi Sebastian,

Thanks for the comments and your code. For now I just want to comment on the following regarding the steady-state case.

> - steady-state conditions (although I have regular dosing); but I am not sure as to how this is > implemented in Torsten as it needs an implicit equation solver which we do not have yet in
> Stan - sundials has such a solver with sensitivity support

The way I have formulated the periodic steady-state case resulting from multiple dosing is in terms of:
y'(t) = f(t, y(t)),
where
y(t) = y(t + tau).
It is sufficient to collapse the latter condition to
y(0) = y(tau).
In other words it is a boundary value problem. In the past I handled that using a numerical root solver to solve for the initial conditions required to satisfy the last equation. Thus the initial value problem is embedded in the root solving problem.

You mention IDAS from the Sundials collection. It solves DAE systems of the form:
F(t, y(t), y'(t)) = 0
where
y(t0) = y0, y'(t0) = y'0.
I haven't figured out how to apply this to the aforementioned boundary value problem. I would be curious to know if you have formulated the steady-state problem in a way that IDAS could handle it. Otherwise I expected we would need to use some combination of the Stan ODE solvers with a root finder like KINSOL (also in Sundials).

Cheers,
Bill

[Sebastian Weber]

Hi!

I have never worked with IDAS, but you are right, I had that one in mind. So no, I have not attempted to solve this problem yet, but I think that IDAS should solve this problem (or at least I would try it with that).

As far as I remember, Ben also wanted to have such a functionality in Stan, so this could be of general interest?

Sebastian

 

You received this message because you are subscribed to a topic in the Google Groups "stan development mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-dev/xGlyeLkdSiE/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-dev+u...@googlegroups.com.

[Bob Carpenter]

No problem---we're establishing process as the
project grows.

We'd like to share with Wiki for design and branch for
development. Sometimes it's necessary to create a pull request
just to do code commenting. We're open to suggestions on better
ways to do things.

We just don't want issues stacking up if they're not ready to
act on, because it adds overhead whenever we go through and triage.

I gave you permission to create branches on stan-dev in all the repos.
Please don't push without pull requests to any of
the core repos (stan, cmdstan, pystan, rstan, math).

- Bob

[Bob Carpenter]

> On Aug 4, 2016, at 9:20 PM, Charles Margossian <charles.ma...@gmail.com> wrote:
>
> Hi Daniel,
>
> Yes, I indeed created a wiki on Strategies to Facilitate Stan in Pharmacometrics, with a section on Torsten. I updated it so that it links to the example models I posted on a fork of example-model.
>
> I also wrote a doc on how to use the examples.

OK --- is the source for the doc there somewhere? We like
to keep all this stuff open and editable by everyone in
branches.

> Today I added a third section, "Under the Hood Design", which I hope will answer your questions.
>
>
> I'm still missing a big piece of what you're doing. In particular, I don't understand how much of this is additions to:
> - stan-dev/math
> - stan-dev/stan
> - stan-dev/cmdstan
>
> At the end of the day, Torsten is four built-in Stan functions. They are grouped together because they have a similar architecture and functionality. Additions are made to stan-dev/math (see here) and stan-dev/stan (see here). stan-dev/cmdstan remains unchanged.

Are they literally functions or are they custom expressions
like integrate_ode()? The differences is that functions go into
the general function-handling syntax whereas things like integrate_ode()
are custom and allow specifying that some arguments are data rather
than parameters, etc.

- Bob

[Bob Carpenter]

> On Aug 5, 2016, at 9:34 AM, Krzysztof Sakrejda <krzysztof...@gmail.com> wrote:
>
> On Friday, August 5, 2016 at 8:46:54 AM UTC-4, Charles Margossian wrote:
>> I agree with you that Torsten should be broken up, especially since many of its features appear to be extendable to fields other than pharmacometrics. The initial idea was to take a set of functions present in NONMEM and BUGS, and see if we could implement them in Stan. Fitting everything into one function was meant to meet a desired and expected simplicity for the user. I believe one of Torsten's strengths is how accessible and convenient it is for pharmacometricians.
>
> That's ultimately a great approach to an interface for scientists. I'm only suggesting making some under-the-hood changes to make that layer thinner and less specialized.
>
>> Evidently, there seems to be a more ambitious and better approach.
>
> My suggestion is to go for something less ambitious first. The Stan development process for bringing in small well-specified chunks of functionality is pretty clear and well greased at this point. At the same time I'm not sure I've ever seen a big ambitious set of changes that went in smoothly. Maybe Bob/Daniel will just work their magic and see how to integrate all this but
>
>> They are not really pharma-specific in the sense that other modeling could use them (e.g.-geochemical systems).
>>
>> I'm trying to learn about ODE based models in other fields. Right now, my focus in on system biology, but if you have any papers or documents on geochemistry, I am very interested. My hope is to get a better sense of what is pharma specific and what is generalizable in Torsten.
>
> Bob has some examples in Stan----I can't find them at the moment but maybe he'll pipe up.

I just have the very first case study on the doc page of our
web page under case studies.

How easy things are going to be get merged is more a matter
of how simple they are than how big they are, but all else
being equal, big pull requests are hard to review. So before
doing four functions, it might be good to do one. We got into
trouble with vectorization trying to do 40+ functions at once.

But I'd rather everyone agree on at least a high-level design
goal before spending too much time coding. The worst possible outcome
is someone putting lots of work into something we can't figure out
how to merge safely.

- Bob

[Bob Carpenter]

Didn't we propose to add a DAE solver? Ben, Michael and
maybe Sebastian are going to need to be involved heavily with
that.

- Bob

[Charles Margossian]

Hi Bob,

 

We'd like to share with Wiki for design and branch for
development.

I created the wiki on Complex ODE Based Models and got rid of the pharmacometrics page. Feel free to edit. Sebastian, could you include a link to the README you wrote for your PKPD library? I would do it, but I want to make sure I have your permission. I'll share the page on the user-list tomorrow, so we can get feedback as to how users would like to see Stan tackle ODE models.

I gave you permission to create branches on stan-dev in all the repos.

Nice! Should I move the code to a branch of stan-dev to make it more visible?

OK --- is the source for the doc there somewhere?  We like 
to keep all this stuff open and editable by everyone in 
branches. 

Just uploaded the source code in the source directory.  

Are they literally functions or are they custom expressions 
like integrate_ode()?

They are custom expressions. Some of the arguments are required to be data. 

So before doing four functions, it might be good to do one

Understood. If one works, the other ones should not be too hard to add. OneCptPKModel() might be a good first case. If we break down Torsten, there may be even simpler "intermediate" expressions we can start with. 

But I'd rather everyone agree on at least a high-level design 
goal before spending too much time coding.  

 

Understood. It seems we are heading towards breaking down and generalizing Torsten, which should not be too difficult to do. 

Thank you for taking the time to clarify development practices!

best,

Charles 

> To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

> For more options, visit https://groups.google.com/d/optout.

--
You received this message because you are subscribed to the Google Groups "stan development mailing list" group.

To unsubscribe from this group and stop receiving emails from it, send an email to stan-dev+unsubscribe@googlegroups.com.

[Sebastian Weber]

Hi!

Just add a link to the README - of course! Whatever is on github/stan-dev is in the public.

Sebastian

[Bill Gillespie]

On Saturday, August 6, 2016 at 8:39:57 PM UTC-4, Bob Carpenter wrote:
> Didn't we propose to add a DAE solver? Ben, Michael and
> maybe Sebastian are going to need to be involved heavily with
> that.

Yes, but as I commented in my post on 8/5 I am not sure the periodic steady-state problem commonly encountered in PK can be handled by the usual DAE solvers, e.g., IDA from SUNDIALS.

Bill