New functionality planning

[isaacsas]

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.

[sebapersson]

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.

[isaacsas]

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?

[TorkelE]

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).

[sebapersson]

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:

1. Compartments are in a sense a parameter class. We could enforce them to be constant, as dynamic compartments make things complicated.
2. Each specie is optionally assigned a compartment. If one specie has a compartment, it probably makes sense to enforce it for each specie.
3. 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?)

[TorkelE]

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.

[isaacsas]

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)?

[sebapersson]

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?

[elbert5770]

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.

[sebapersson]

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.

[elbert5770]

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.

[TorkelE]

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.

[TorkelE]

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.

[elbert5770]

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.

[sebapersson]

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?

[TorkelE]

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?