Catalyst.jl
Catalyst.jl copied to clipboard
SciML
New functionality planning
Now that we have things working well again let's brainstorm what our next functionality efforts should be:
- ODE-jump correctness tests (then it is done)
- functional affects in the DSL
- database of models with solutions for correctness testing
- hybrid with leaping
- hybrid with SDEs
- dynamic hybrid
- expand / improve spatial models
- support for delays
- GPU Gillespie
- improved stability analysis
- explicit compartment support and API
Docs:
- Network docs @vyudu
- Gillespie docs
- callback docs
- LikelihoodProfiler docs
- FiniteStateProjection docs
- NeuralNetwork + CRN docs.
- Update the current parameter fitting example to be practically identifiable.
- Doc example with full parameter fitting workflow.
I believe it would be great to add optional compartment (e.g. cytosol and nucleus volumes). Manually handling compartment scaling is easy to get wrong, and compartments are among the few SBML features still missing in Catalyst.jl, so adding support would in the future make it easier to write a ReactionNetwork → SBML exporter in the future.
database of models with solutions for correctness testing
For this it could be worthwhile to glance at the stochastic SBML test-suite that covers many cases relevant for Catalyst.
Do you have any material / thoughts on how compartments should be included? What changes should ReactionSystems have, are compartments just metadata associated with a species and.or reaction, or do compartments store a ReactionSystem?
I think we had some discussion about compartments and volumes at some point a while ago. I don't remember exactly how it went, but I thin the conclusion was "This is something that we would want to have, but actually, how to implement it is not as obvious as one would immediately think". At some point we should sit down again and think through if there is a good solution (or if not, what is the best non-ideal solution).
I agree with Torkel that compartments is something that needs to be thought about carefully (from implementing SBMLImporter I encountered many edge cases). But, I think it would be best to implement them in a similar way to SBML:
- Compartments are in a sense a parameter class. We could enforce them to be constant, as dynamic compartments make things complicated.
- Each specie is optionally assigned a compartment. If one specie has a compartment, it probably makes sense to enforce it for each specie.
- If an event changes a compartment, species need to be re-scaled.
I think it is also worthwhile to consider if compartments make sense for discrete JumpProblems (that is, does JumpProblem make sense for systems where we are modeling concentrations?)
I thinkl one thing we discussed was how this ties into hierarchical reaction systems. I.e. it would make sense to enforce a single reaction system to be maximally one compartment, and then compose these to make multi-comaprtment models. But then you have cross compartment reactions. I have forgotten exactly how far we got, by I remember that it really wasn't obvious.
Here are some natural questions/complexities that compartments then introduce:
- Should reactions also be tied to compartments too, or just be possible in any compartment that contains the needed substrate species?
- Should parameters also be tied to compartments (i.e. a reaction rate could be different in two different compartments)?
- Are species allowed to be defined in more than one compartment, or do users need to introduce new species for each compartment (for example, to model cystolic vs. nuclear mRNA)?
- Rescalings are only needed in models that are based on concentrations. So do we need a way to let users specify they want species variables to be interpreted as concentrations somewhere in the workflow (even if they do not provide units)? In this case what should Catalyst be injecting vs. assuming a user has explicitly defined (i.e. right now can't one construct the appropriate volume rescalings for a reaction that represents transport between compartments when one wants to use concentration units)? Likewise, if one has an event changing volumes one can manually update concentrations in the event, but it seems like it would be hard for Catalyst to detect this and extend the event to change concentrations (i.e. I guess we would need to explicitly support some type of "change volume" operation that knows how to update everything that users are told to call instead of them manually changing a volume)?
Should parameters also be tied to compartments (i.e. a reaction rate could be different in two different compartments)?
Probably not, as this is not the SBML case
Are species allowed to be defined in more than one compartment, or do users need to introduce new species for each compartment (for example, to model cystolic vs. nuclear mRNA)?
No, each specie would have a compartment, and splitting RNA into species like into nuclear and cytosol is more clear. We would have to work out the scaling math also, as for besides import or export reactions, the scaling is a bit more tricky (e.g. when a nuclear protein aids importing RNA).
Compartments tied to species also means reactions per say are not tied to compartments (which I think would be the most easy solution)
Rescalings are only needed in models that are based on concentrations. So do we need a way to let users specify they want species variables to be interpreted as concentrations somewhere in the workflow (even if they do not provide units)? In this case what should Catalyst be injecting vs. assuming a user has explicitly defined (i.e. right now can't one construct the appropriate volume rescalings for a reaction that represents transport between compartments when one wants to use concentration units)? Likewise, if one has an event changing volumes one can manually update concentrations in the event, but it seems like it would be hard for Catalyst to detect this and extend the event to change concentrations (i.e. I gues s we would need to explicitly support some type of "change volume" operation that knows how to update everything that users are told to call instead of them manually changing a volume)?
How about, model with compartments is assumed to be concentration? As for events, a special class would likelly be best as here we also have the question priority, is volume changed before we also potentially change a specie?
I'll point out that if you have changing compartment volumes and the state variables are concentrations, you need to apply the quotient rule as d(amount/V)/dt = (V d(amount)/dt - amount dV/dt)/ V^2. I think this is why Roadrunner uses amounts instead of concentrations as state variables. I can imagine that generic solvers that read SBML files might get this wrong.
I'll point out that if you have changing compartment volumes and the state variables are concentrations, you need to apply the quotient rule as d(amount/V)/dt = (V d(amount)/dt - amount dV/dt)/ V^2. I think this is why Roadrunner uses amounts instead of concentrations as state variables.
This is a very good point, and something for example importers like SBMLImporter deals with. But I also think it should, given that each specie has a unit, should be something Catalyst can deal with.
This is a very good point, and something for example importers like SBMLImporter deals with. But I also think it should, given that each specie has a unit, should be something Catalyst can deal with.
I think there might also be a place for a special reaction symbol that denotes flow between compartments (i.e. the units of the rate constant are volume/time). Something like ~~> for unidirectional flow and <~~> for bidirectional flow.
I have a decent chunk of thoughts on compartments, but am travelling and mostly on my phone. Will write something when I get back home.
So, generally, I think it makes most sense for compartmental models to be separate ReactionSystem models, which are then composed together to from the final, full model. Here, to specific arrows are really needed for reactions between compartments, these are just normal reactions where substrate/products come from different compartments.
In most cases you would likely have a top-level model which contains things like inter-compartmental reactions, and each compartment is a sub-system.
One questions whether there should be a special "compartment" tag for a ReactionSystem, in which case one might be able to define some additional properties, but also disallow some stuff.
One questions whether there should be a special "compartment" tag for a
ReactionSystem, in which case one might be able to define some additional properties, but also disallow some stuff.
It is noteworthy that the way that Roadrunner and SBMLtoODEjax/SBMLtopyODE handle compartments is that by default the state variables are amounts, not concentrations. Each state variable is divided by the compartment volume at the current time step to yield a concentration, so that the rate constants can still have units with concentrations when appropriate. The trick is that each equation must then be multiplied by the compartment volume somewhere: dm1/dt = -k_2Om1/Vcompm2/VcompVcomp. I know it looks silly but it does make dynamic volumes easier. Flows and clearances naturally work out as: dm_origin/dt = -Qm_origin/V_origin_comp and dm_destination/dt = Q*m_origin/V_origin_comp. The trick is who multiplies or divides by Vcomp, is it the user or the interpreter. SBMLtoODEjax makes the user write the rate law as amount/time, but in the background divides any amounts in the rate law by volume to get a concentration. Antimony/Roadrunner does both steps in software. It might be good to have a zoom to discuss.
In most cases you would likely have a top-level model which contains things like inter-compartmental reactions, and each compartment is a sub-system.
Still think it would be worthwhile to follow SBML as closely as possible and not allow the same specie id in different compartments. Also, would the suggested approach work if import into, for example nucles, depends on concentration of a nuclear protein?
Yes, that approach should work fine, e.g.
@reaction_network begin
@compartment nucleus begin
@species nuc_prot(t)
(p,d), 0 <--> X
end
@compartment cytosol begin
d, X --> 0
end
k*nucleus.nuc_prot, nucleus.X --> cytosol.X
end
(something like this, but programmatic and no in the DSL already work).
I think species will always be tied to compartment though, so not sure what you mean. The question is where inra-compartment reactions live, what do SBML do here?
I think species will always be tied to compartment though, so not sure what you mean. The question is where inra-compartment reactions live, what do SBML do here?
SBML forces each specie, in each compartment to have a unique specie-id. So in your example it would be to have:
nucleus_X , cytosol_X instead of nucleus.X , cytosol.X
I think you get the same ids/names in Catalyst, but this is a nice way to reference them. I think this is nice since it someone how include the concept of the species bring the same one, but in different compartments
Some related docs https://docs.sciml.ai/Catalyst/stable/model_creation/compositional_modeling/ https://docs.sciml.ai/Catalyst/stable/model_creation/reactionsystem_content_accessing/#model_accessing_hierarchical
I think this is a fair point. It could also be dealt with with a potential SBMLExporter.jl package, if it is okay to not have a 1-1 mapping between exported and then imported model?
So I think in this case there would be? I.e. after complete the system is automatically flattened, and when this happens each species have their own name.
So I think in this case there would be? I.e. after complete the system is automatically flattened, and when this happens each species have their own name.
That would be great!
And yes, @elbert5770 It might be useful to discuss on Zoom some time. Do you have a link to how Antimony does it?
So @TorkelE basically you are just saying that compartments are really just sub-ReactionSystems?
I still think we likely need to a bit more than that, so that when we flatten we can keep track of the compartments the model was constructed from and their properties. (Like their dimensionality and size.)
Well, I think my point is that we already support sub system, and I think that we probably should base our support on compartments around this. But yes, we need some additional framework on it as well for volumes etc.
Do you have a link to how Antimony does it?
Antimony is just a text file or string that Tellurium converts to SBML before solving the SBML model with Roadrunner. Catalyst is different because its goal is to write a symbolic model. I'm going to meet with Lucian Smith soon to discuss compartments in Roadrunner/SBML. I want to first develop some tests that illustrate some of the more challenging cases. I suspect that there needs to be a tag in SBML that applies 'hasonlySubstanceunits' to reactions instead of species for compartments to be handled correctly in all cases. I will keep you posted as I think I will learn a great deal over the next week. There is an interesting problem in that Catalyst->SBML->Catalyst needs to give the same answer as Catalyst alone.
Antimony seems more inline with what @TorkelE is suggesting as it allows one to attach both species and reactions to compartments.
Each species is required to have a compartment in Tellurium, even if it is a 'default' compartment with volume = 1. Antimony allows attaching species to compartments to avoid the 'default' compartment. Reactions are not tied to compartments in SBML I think? Thus, reactions are not tied to compartments in Antimony, other than the fact that the species need to be in the same compartment. I'm not sure if this is checked though. Labeling species names with compartments is just good practice but not 100% required.
Antimony’s documentation seems to indicate reactions can be assigned to compartments too: https://tellurium.readthedocs.io/en/latest/antimony.html#id36
I met with Lucian Smith and Herb Sauro yesterday and an interesting insight was that the hasonlySubstanceunits tag is implemented in Antimony/roadrunner in a logical but non-intuitive way. The default value of hasonlySubstanceunits is false. Also, the rates (LHS) are ALWAYS amount/time. But similar to SBMLtoODEjax, the right hand side species are ALWAYS concentration if hasonlySubstanceunits is false. So in the background the amounts on the LHS are divided by compartment volume to yield a concentration for use on the RHS (to be honest, the actual algorithms are hard to decipher from the code because they write LLVM instead of equations so could be more complex). However, the RHS (rate eqtn) MUST deliver amount/time. It is up to the user to multiply by the compartment volume to convert the conc/time on the RHS to amount/time. You will otherwise get the wrong answer unless vol=1. If hasonlySubstanceunits is true, the amount/time on the LHS is used to update amounts on the RHS, so the equations can be written 'normally', i.e. without multiplying by compartment volume. However, the rate constants and IC need to be in amounts not concentrations. That's tricky if volume is changing. In terms of flows and clearances (vol/time) it is better to have hasonlySubstanceunits=false but of course the user then needs to be careful with amount vs. concentration and this means the rate equations in SBML need to be like k x conc x vol for 1st order reactions IN THE SBML. Catalyst is perhaps better positioned because only the rate constants are specified and rate equations are constructed given mass action or another function. Catalyst could write the equations correctly for the user, perhaps even in the background. However, Catalyst defaults units to whatever is in the rate constant/IC, I believe, which if true requires a graph of possible outcomes. Definitely worth discussing more. How should we schedule the proposed zoom? Perhaps on slack?
By default, Catalyst currently let's everything be unitless. I think this is a good default behaviour, and should work fine with compartments (i.e. it already does) if volumes are not specified.
I am unsure how unitless units should interact with volumes if this is defined (i.e. some units should be automatically assumed, error, or somehow carry through with it in some fashion anyway). If users want to define units they should use the system already in place, and from there we should be able to figure stuff out for the user (and attach units to any outputs).
I am pretty much happy to have Zoom meetings whenever, mostly free (on London time). I should mention though that we will unlikely be able to make any implementation anytime soon. E.g. the hybrid simulation stuff is ready, but we would want to add more tests and write docs for it. Then there is MTKv10 which we would want to update Catalyst to work with. Finally, both me and Sam have a decent amount of other stuff to work on in between which to squeeze in Catalyst development.