Feature Request: Providing more reaction information in types generated by @reaction_network

[isaacsas]

Many SSAs (Gibson-Bruck, partial propensity methods, RSSA...) require more reaction network info than is easily available from the type generated by the @reaction_network macro. It would be great if the generated type could also include fields that provide a mapping from reaction id to the coefficients of reactant species, and a mapping from the reaction id to the coefficients of product species. Such info is needed to determine which reactions need their propensities updated when another reaction occurs or a species changes its population size.

For a bit more detail, many of these other SSAs require knowing things like

Reactants( rx id ) = set of reactants for a given reaction Products( rx id) = set of products for a given reaction DependentRxs( rx id) = set of reactions with jump rates depending on species populations that "rx id" modifies.

All these types of sets can be built from the two mappings I mentioned.

Sam

[TorkelE]

Yes, to be honest I think the part of the DSL used for SSAs needs some serious attention. I have been meaning to sit down (likely over a weekend) and have a look at these kinds of issues (admittedly no great speciality of mine) and how the reaction networks creates jump problems.

To begin with ConstantRateJumps are currently only created, while if e.g. functions of t are used there should be VariableRateJumps (?). This is a plain issue that I need to fix.

Adding fields like you are suggesting is no problem, but in the interest of my own learning it would be interesting to know exactly how you plan to use this knowledge?

[isaacsas]

I'm starting to implement some different SSAs (as both a Julia learning exercise, and to get a better understanding about implementation details that influence performance). There are lots of SSAs that have been published in the past 10-15 years, with a number that claim substantial performance gains over the direct method for all but the simplest systems. I'd like to get a better feel for how they perform on different types of networks beyond just the theoretical asymptotic work per event. Ultimately, I'm planning to then use / adapt / extend some of these for my own research (which uses spatial SSAs where molecules diffuse and react).

Recent SSA methods are almost uniformly designed/optimized to avoid re-evaluating as many propensities (i.e. jump rates) as possible after an event occurs (as this updating is usually the rate limiting step per event). For this reason, they often rely on data structures that do things like map the ID of a reaction to a list of other "dependent" reactions whose propensities should be recalculated when the first reaction occurs (say because the first reaction changes the number of reactants that the other reactions would need, hence modifying their propensity/jump rate). The idea being that one wants to minimize reevaluating as many jump rates per event as possible and/or update associated data structures as quickly as possible.

I would also like to maintain input compatibility with DiffEqs.jl, since that would allow me to simultaneously look at ODE solutions or SDE solutions to the chemical system and/or to use DiffEqs.jl SSA methods too as they expand/improve.

-Sam

[ChrisRackauckas]

I would also like to maintain input compatibility with DiffEqs.jl, since that would allow me to simultaneously look at ODE solutions or SDE solutions to the chemical system and/or to use DiffEqs.jl SSA methods too as they expand/improve.

DiffEq is a distributed set of solvers from a bunch of different packages. If you add a package which adds an algorithm to the common interface, it can be called like any other. A package of SSA methods would be great to have! When the time comes we can talk about migrating it to DiffEq and giving you membership if you're interested, but being in DiffEq isn't required to be callable from DiffEq (see LSODA.jl).

@alanderos91 wrote a bunch of SSA methods in BioSimulator.jl so maybe he can weigh in.

https://github.com/alanderos91/BioSimulator.jl

[isaacsas]

If down the line I end up with a good package of SSAs that would certainly be nice! For now I'm still learning Julia, and don't really yet understand the general program flow used internally by DiffEqs.jl when setting up and running an SSA simulation. (I have trouble understanding where all the various dispatches route me in the overall DiffEqs.jl codebase.)

[alanderos91]

It looks like both ReactantStruct and ReactionStruct already have everything needed to construct a dependency graph for the propensities (see here). It looks like the macro passes around aReactionStruct vector so you can just use the indices for the IDs. Are these data structures passed down to the generated type?

Here's an example of how I construct dependency graphs in BioSimulator.jl. All the reaction and substrate information is stored in a dictionary.

function make_dependency_graph!(dg, reactions_dict)
  reactions = values(reactions_dict)

  for (i, reaction) in enumerate(reactions)
    for (j, other) in enumerate(reactions)
      affected_species = union(keys(reaction.reactants), keys(reaction.products))
      other_reactants  = keys(other.reactants)

      if !isempty(intersect(affected_species, other_reactants))
        push!(dg[i], j)
      end
    end
  end

  return dg
end

[isaacsas]

Related features have been in DiffEqBiological for a while now.