Including global variables that participate in MassActionJumps in spatial systems

[alexeychizh]

I have a grid (membrane) where species (molecules) can hop around, react, and interact/exchange with a global, well-mixed pool (cytoplasm). Species A, B and C can transition between the global pool and grid sites with rates k_on, k_off.

Suppose the mass-action reaction we consider is A + B <--> C. Alongside reactions between A and B molecules at specific grid sites, how can I include a MassActionJump that consumes, say, a B in a specific grid site and an A in the global pool (that is equally accessible to all the grid sites), to create +1 C in the grid site? Likewise, the reverse case, where a grid site C can dissociate to form a global pool A and grid site B.

Very similarly, I would like to include recruitment, where a global pool A can be recruited to a grid site at a rate proportional to the number of A at that site (or, from the perspective of the grid site, A --> 2A at rate k_fb * A_pool, where A_pool is the population of A in the pool, which is depleted by 1 with each recruitment).

I would ideally like to do this with the existing spatial solvers.

[Vilin97]

Say, you want a reaction that consumes species s in site i (assuming the species are in the [A,B,C] order, s=2 would be species B). You would construct the stoichiometry vectors the usual way:

netstoch = [[s => -1]] # indicates that species s goes down by 1 when this reaction happens
reactstoch = [[s => 1]] # this reaction has a single reactant -- species s

And would would pass in a matrix for reaction rates, whose r, i entry is the reaction rate of reaction r at site i.

rates = zeros()
rates[r,i] = 1.0 # or whatever you want your rate to be.
majump = SpatialMassActionJump(rates, reactstoch, netstoch)

To change the reaction from B -> 0 to A+B -> C just change the stoichiometry vectors.

You might wonder how you can know what site has what index. If you are using a CartesianGrid, then it's just the index of the array, presented as a vector. For example, CartesianGrid( (2,3) ) has site 2 in the lower-left corner (Julia is column-major). But in your case you want to use a Graph because you want not just a grid but a grid with an extra node. You can construct it using g = Graphs.grid( (a,b) ) and then adding a vertex with Graphs.add_vertex!(g) and adding all edges with Graphs.add_edge!(g, s, d), where you loop over all nodes s in g and d is the newly added special node. In this case i is the index of the node in the graph.

Hopefully, this helps. You might also want to look at the tests (e.g. this one), this tutorial and this comment

[alexeychizh]

This doesn't address my question. Perhaps I was unclear: I want to include a MassActionJump or SpatialMassActionJump that consumes, or produces, species in two different 'sites' in one event.

So, let me be clearer and describe species as s[i], where s is species and i is its location in the Graph. Let's use the topology you suggest - my grid site positions i run from 1 to a*b, and I attach an extra vertex (position a*b + 1), corresponding to my global pool, connected to all the nodes in the grid.

I want to include reactions where A[a*b + 1] + B[i] --> C[i] and its reverse C[i] --> A[a*b + 1] + B[i], where i \in 1:a*b. In other words, an A is recruited from the global pool by B at a specific node to form C at the same node. The rates of these reactions would similarly scale with the populations of A[a*b + 1] and B[i].

Alternatively, I could solve this problem without adding this extra vertex if I could have a global variable that refers to the population of A in the global pool - let's call it A_pool. I could then have a reaction B[i] --> C[i] whose rate scales with A_pool, but, crucially, A_pool -= 1 with each event. Similarly, I have a reaction C[i] --> B[i], whose rate is independent of A_pool, but each time it fires, A_pool += 1.

I don't know how I could solve this problem while also being able to use the specialised spatial solvers NSM and DirectCRDirect. Let me know if anything is unclear in my explanation.

[isaacsas]

I don't think this can be done currently unfortunately, but @Vilin97 can hopefully better comment on this. Allowing such reactions is definitely on our TODO list, but I'm not sure we can say it will be available anytime soon.

Basically we ultimately would like to allow

• Coupled well-mixed reactions and species
• Multiple sub-domains (with coupling transitions between them allowed)
• Non-local reactions i.e. of the form A_i + B_j -> C_k.

But I don't think any of these are supported currently.

[Vilin97]

Non local reactions are not implemented and are not currently in the works AFAIK, although it is a long term TODO. But can this special case perhaps be done with a callback that triggers when B[i] -> C[i] fires and decrements A_pool by 1, @isaacsas ?

[isaacsas]

How is the information in the change in the rate due to the change in the number of A going to be propagated across all spatial locations?

[isaacsas]

@Vilin97 I was thinking about non-local reactions. What do you think about formalizing our "update" API, and then making it available to users via affect functions. The idea being that they have to take responsibility for registering the "updates" they've made, and then we handle the efficient updating using that information inside the internal data structures?

It seems like something like that might give the most flexibility in the types of reactions we can accept, since then users can tell us the updates that have been made (i.e. species and locations).

[Vilin97]

@isaacsas, I missed your comment about the change in number of A propagating across locations. I guess a hacky way to do this would be to use constant rate jumps from https://github.com/SciML/JumpProcesses.jl/pull/343. It seems that you are suggesting something similar when you talk about the "update" API, although I did not fully understand what you meant.

[isaacsas]

But the constant rate PR doesn't allow users to register where they are making updates right? It assumes they are made at a specific node only I thought? I was suggesting we consider/design a way to register that one has made updates of specific species at specific locations within an affect function (and use this ourselves in the solvers), and then we process this list automatically.

This would allow a user to say that they've updated species 10 at location 1 and species 2 locations 10,15 for example, and we would know where our data structures need updating. It seems like this is a more general way to work towards allowing non-local reactions.

With that said, I don't have a good feel if this can be done in a way that doesn't have a significant impact on speed.

For @alexeychizh's problem, what would probably be the most efficient way to do it would be if he could have one constant rate jump for selecting the reaction across the whole domain, and an affect function he writes that samples the location of the reaction when it occurs somewhere (which would require access to our internal information probably, and a way to tell us where he is making updates). This would then only need to update the one "global" reaction, and the one site where the reaction occurs, within our data structures.

I don't see how in our current constant rate setup we wouldn't require recalculating rates everywhere, since this reaction would have to be represented as A + B[i] --> C[i], and changing the amount of A would require recalculating this rate at all sites. (In contrast to what I say above.)

[Vilin97]

Now I see what you are saying. After merging the three outstanding PRs let's think what the interface should look like and then I'll try implementing it.

[isaacsas]

Sounds good. I'm working on Torkel's spatial PR again finally, hopefully we will merge it soon even it is not really in a good final form, and then I will circle back to your PRs.