New learner here. Thanks to help of several generous contributors I have a working “toy” of a two-pool compartmental mechanistic biological model, produced in Rust, using the “ode_solvers” crate. When I get this completed I will expand it to an actual 18-20 pool model.
My current problem is an inability to pull out the values of some non-state variables (results) at each iteration of the ode solver. I want to export these variable values to both a text file and to a simple onscreen graph. I can export and graph the State variable values, just not the non-state ones.
On “the bigger picture”, I might be able to attack this problem if I had a way to interrogate the the inner workings (options) of the ode_solver crate but I don’t yet know how to do that. Any advice on how to proceed in this area and learn myself most welcome!
Here is a sample of code from a working mwe, happy to supply the github repo address if anyone is interested.
//declare system function, used by the stepper above
struct TwoPool;
impl ode_solvers::System<f64, State> for TwoPool {
fn system(&self, time: Time, y: &State, dy: &mut State) {
/*calculate the concentration of metabolite in each pool at
each iteration*/
let con_a = y[0] / I[0];
let con_b = y[1] / I[1];
/*calculate fluxes corresponding to arrows on diagram, using a
HMM equation*/
let fab = C[0] / (1.0 + (C[3] / con_a));
let fba = C[1] / (1.0 + (C[4] / con_b));
let fbo = C[2] / (1.0 + (C[5] / con_b));
/* specify the differential equations for each of the state
variables*/
dy[0] = I[4] + fba - fab;
dy[1] = fab - fba - fbo;
// y[2] = con_a;
If I uncomment the last line and attempt to store the iterative values of “con_a” in the vector y (state variables) it blows up with a non-mutable vector error (not surprisingly!, my bad)
Main issue: Any suggestions on how to extract non-state variable values, eg. con_a and con_b and add them to my graph, which currently graphs the state variable values very nicely. Thx. J
I encountered a similar limitation with ode_solvers crate, but on a different domain (mechanical system). The easiest way I found around it was to use a different ODE solver. I found in particular Russel’s approach more to my liking. In particular, Russel_ODE allows you to pass parameters to the ODE solver call (i.e., your non-states).
@pzometa Thanks very much for that, when learning a new language it is difficult to ascertain whether the evolving language (or crates) has shortcomings or is it my inability. This helps a lot. I don’t suppose you a public example or repo showing how you used the ode solver? Thx. J
I just simplified an old code that I had for a mechanical system (first commit). Using an AI-agent, it was adapted to your example (biological system, second commit). Note that I use a constant time step to advance the solver. Not sure if that is acceptable in your case. I hope it helps.
@pzometa Thanks for this …. unfortunately I just get a 404 error. For interest my model is located at git@github.com:jamaas/jam_two_pool_no_total.git
The diagram is in pdf … might save some time. Please note that this example is trivial but when working correctly can be used as a template for much more complex biological models. Thx. J.
@pzometa That appears to work just great ! I need to do a couple of checks and then figure out how to draw rudimentary graphs to screen in Linux. Will be in touch, should I publish the location of the public github repo on here? Best, J.
Just very simple runtime ones, in past I have found ... for me, the combination of numerical data and graphs helps identify and correct simple errors much more quickly. Probably more importantly, are great to demo the model for seminars to people who are not familiar with the model or the practice. Something like attached would be great! It might be difficult in my case to graph both Quantities and Concentrations on the same graph because the scales are so different. Perhaps two y axis? Thx. J.
if interested clone or copy git@github.com:jamaas/jam_two_pool_russell_zometa.git
Maybe to address the initial question more directly:
That may be a use case for RefCell to store a trace of your con_a variables into the system itself.
struct TwoPool {
// … w/e you have
pub trace: RefCell<Vec<f64>>,
}
// Later:
impl ode_solvers::System<f64, State> for TwoPool {
fn system(&self, time: Time, y: &State, dy: &mut State) {
// … previous code
// not entirely sure about the types, `State` isn't listed in your code
self.trace.borrow_mut().push(con_a);
}
}
This works as long as the system runs on a single thread. The introduction of RefCell will mark the struct as !Sync, meaning you can't share it across two threads at the same time. You can swap out RefCell for Mutex which will be slower but work across threads (if necessary). I don't think it is, the ode_solvers implementations that I see from a quick look across the crate run on one thread.
