Skip to content

Commit d15a28b

Browse files
Vilin97Vilin97
committed
added a tutorial for constructing, modifying and visualizing spatial reversible binding model
1 parent 2cc7ef3 commit d15a28b

File tree

2 files changed

+229
-0
lines changed

2 files changed

+229
-0
lines changed

‎tutorials/jumps/assets/ABC_anim_275frames_5fps.gif‎

1.29 MB
Loading

‎tutorials/jumps/spatial.jmd‎

Lines changed: 229 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,229 @@
1+
---
2+
title: "Tutorial on using spatial SSAs in DiffEqJump"
3+
author: Vasily Ilin"
4+
---
5+
6+
This tutorial shows how to use spatial solvers added to [`DiffEqJump`](https://github.com/SciML/DiffEqJump.jl) in summer 2021. See [the documentation](https://diffeq.sciml.ai/latest/types/jump_types/) for a tutorial on getting started with `DiffEqJump`.
7+
8+
## Installing `DiffEqJump`
9+
10+
Once in REPL, do `] add DiffEqJump`. After the installation finishes, you will be able to use all the functionality described below.
11+
12+
## Reversible binding model on a grid
13+
14+
A 5 by 5 Cartesian grid:
15+
16+
| <!-- --> | <!-- --> | <!-- --> | <!-- --> | <!-- --> |
17+
|---|---|---|---|---|
18+
| . | . | . | . | B |
19+
| . | . | . | . | . |
20+
| . | . | . | . | . |
21+
| . | . | . | . | . |
22+
| A | . | . | . | . |
23+
24+
Suppose we have a reversible binding system described by $A+B \xleftrightarrow[k_2]{k_1} C$, where $k_1$ is the forward rate and $k_2$ is the backward rate. Further suppose that all $A$ molecules start in the lower left corner, while all $B$ molecules start in the upper right corner of a 5 by 5 grid. There are no $C$ molecules at the start.
25+
26+
We first create the grid:
27+
28+
```julia
29+
using DiffEqJump
30+
dims = (5,5)
31+
num_nodes = prod(dims) # number of sites
32+
grid = CartesianGrid(dims) # or use LightGraphs.grid(dims)
33+
```
34+
35+
Now we set the initial state of the simulation. It has to be a matrix with entry $(s,i)$ being the number of species $s$ at site $i$ (with the standard column-major ordering of the grid).
36+
37+
```julia
38+
num_species = 3
39+
starting_state = zeros(Int, num_species, num_nodes)
40+
starting_state[1,1] = 25
41+
starting_state[1,2] = 25
42+
```
43+
44+
We now set the time-span of the simulation and the reaction rates. These can be chosen arbitrarily.
45+
46+
```julia
47+
tspan = (0.0, 2.0)
48+
rates = [3.0, 0.05] # k_1 = rates[1], k_2 = rates[2]
49+
```
50+
51+
Now we can create the `DiscreteProblem`:
52+
53+
```julia
54+
prob = DiscreteProblem(starting_state, tspan, rates)
55+
```
56+
57+
Since both reactions are [massaction reactions](https://en.wikipedia.org/wiki/Law_of_mass_action), we put them together in a `MassActionJump`. In order to do that we create two stoichiometry vectors. The net stoichiometry vector describes which molecules change in number and how much after each reaction; for example, `[1 => -1]` is the first molecule disappearing. The reaction stoichiometry vector describes what the reactants of each reaction are; for example, `[1 => 1, 2 => 1]` would mean that the reactants are one molecule of type 1 and one molecule of type 2.
58+
59+
```julia
60+
netstoch = [[1 => -1, 2 => -1, 3 => 1],[1 => 1, 2 => 1, 3 => -1]]
61+
reactstoch = [[1 => 1, 2 => 1],[3 => 1]]
62+
majumps = MassActionJump(rates, reactstoch, netstoch)
63+
```
64+
65+
The last thing to set up is the hopping constants -- the probability per time of an individual molecule of each species hopping from one site to another site. In practice this parameter, as well as reaction rates, are obtained empirically. Suppose that molecule $C$ cannot diffuse, while molecules $A$ and $B$ diffuse at probability per time 1 (i.e. the time of the diffusive hop is exponentially distributed with mean 1). Entry $(s,i)$ of `hopping_constants` is the hopping rate of species $s$ at site $i$ to any of its neighboring sites (diagonal hops are not allowed).
66+
67+
```julia
68+
hopping_constants = ones(num_species, num_nodes)
69+
hopping_constants[3, :] .= 0.0
70+
```
71+
72+
We are now ready to set up the `JumpProblem` with the Next Subvolume Method.
73+
74+
```julia
75+
alg = NSM()
76+
jump_prob = JumpProblem(prob, alg, majumps, hopping_constants=hopping_constants, spatial_system = grid, save_positions=(true, false))
77+
```
78+
79+
The `save_positions` keyword tells the solver to save the positions just before the jumps. To solve the jump problem do
80+
81+
```julia
82+
solution = solve(jump_prob, SSAStepper())
83+
```
84+
85+
## Animation
86+
87+
Visualizing solutions of spatial jump problems is best done with animations.
88+
89+
~~~
90+
<figure>
91+
<img src="/assets/ABC_anim_275frames_5fps.gif" width="100%" />
92+
</figure>
93+
~~~
94+
95+
This animation was produced by the following script.
96+
97+
```julia
98+
using Plots
99+
is_static(spec) = (spec == 3) # true if spec does not hop
100+
"get frame k"
101+
function get_frame(k, sol, linear_size, labels, title)
102+
num_species = length(labels)
103+
h = 1/linear_size
104+
t = sol.t[k]
105+
state = sol.u[k]
106+
xlim=(0,1+3h/2); ylim=(0,1+3h/2);
107+
plt = plot(xlim=xlim, ylim=ylim, title = "$title, $(round(t, sigdigits=3)) seconds")
108+
109+
species_seriess_x = [[] for i in 1:num_species]
110+
species_seriess_y = [[] for i in 1:num_species]
111+
CI = CartesianIndices((linear_size, linear_size))
112+
for ci in CartesianIndices(state)
113+
species, site = Tuple(ci)
114+
x,y = Tuple(CI[site])
115+
num_molecules = state[ci]
116+
sizehint!(species_seriess_x[species], num_molecules)
117+
sizehint!(species_seriess_y[species], num_molecules)
118+
if !is_static(species)
119+
randsx = rand(num_molecules)
120+
randsy = rand(num_molecules)
121+
else
122+
randsx = zeros(num_molecules)
123+
randsy = zeros(num_molecules)
124+
end
125+
for k in 1:num_molecules
126+
push!(species_seriess_x[species], x*h - h/2 + h*randsx[k])
127+
push!(species_seriess_y[species], y*h - h/2 + h*randsy[k])
128+
end
129+
end
130+
for species in 1:num_species
131+
scatter!(plt, species_seriess_x[species], species_seriess_y[species], label = labels[species], marker = 6)
132+
end
133+
xticks!(plt, range(xlim...,length = linear_size+1))
134+
yticks!(plt, range(ylim...,length = linear_size+1))
135+
xgrid!(plt, 1, 0.7)
136+
ygrid!(plt, 1, 0.7)
137+
return plt
138+
end
139+
140+
"make an animation of solution sol in 2 dimensions"
141+
function animate_2d(sol, linear_size; species_labels, title, verbose = true)
142+
num_frames = length(sol.t)
143+
anim = @animate for k=1:num_frames
144+
verbose && println("Making frame $k")
145+
get_frame(k, sol, linear_size, species_labels, title)
146+
end
147+
anim
148+
end
149+
# animate
150+
anim=animate_2d(solution, 5, species_labels = ["A", "B", "C"], title = "A + B <--> C", verbose = false)
151+
fps = 5
152+
name = "ABC_anim_$(length(solution.u))frames_$(fps)fps.gif"
153+
path = joinpath(name)
154+
gif(anim, path, fps = fps)
155+
```
156+
157+
## Making changes to the model
158+
159+
Now suppose we want to make some changes to the reversible binding model above. There are three "dimensions" that can be changed: the topology of the system, the structure of hopping rates and the solver. The supported topologies are `CartesianGrid` -- used above, and any `AbstractGraph` from `LightGraphs`. The supported forms of hopping rates are $D_{s,i}, D_{s,i,j}, D_s * L_{i,j}$, and $D_{s,i} * L_{i,j}$, where $s$ denotes the species, $i$ -- the source site, and $j$ -- the destination. The supported solvers are `NSM`, `DirectCRDirect` and any of the standard non-spatial solvers.
160+
161+
### Topology
162+
163+
If our mesh is a grid (1D, 2D and 3D are supported), we can create the mesh as follows.
164+
165+
```julia
166+
dims = (2,3,4) # can pass in a 1-Tuple, a 2-Tuple or a 3-Tuple
167+
grid = CartesianGrid(dims)
168+
```
169+
170+
The interface is the same as for [`LightGraphs.grid`](https://juliagraphs.org/LightGraphs.jl/latest/generators/#LightGraphs.SimpleGraphs.grid-Union{Tuple{AbstractVector{T}},%20Tuple{T}}%20where%20T%3C:Integer). If we want to use an unstructured mesh, we can simply use any `AbstractGraph` from `LightGraphs` as follows:
171+
172+
```julia
173+
using LightGraphs
174+
graph = cycle_digraph(5) # directed cyclic graph on 5 nodes
175+
```
176+
177+
Now either `graph` or `grid` can be used as `spatial_system` in creation of the `JumpProblem`.
178+
179+
### Hopping rates
180+
181+
The most general form of hopping rates that is supported is $D_{s,i,j}$ -- each (species, source, destination) triple gets its own independent hopping rate. To use this, `hopping_constants` must be of type `Matrix{Vector{F}} where F <: Number` (usually `F` is `Float64`) with `hopping_constants[s,i][j]` being the hopping rate of species $s$ at site $i$ to neighbor at index $j$. Note that neighbors are in ascending order, like in `LightGraphs`. Here is an example where only hopping up and left is allowed.
182+
183+
```julia
184+
hopping_constants = Matrix{Vector{Float64}}(undef, num_species, num_nodes)
185+
for ci in CartesianIndices(hopping_constants)
186+
(species, site) = Tuple(ci)
187+
hopping_constants[species, site] = zeros(outdegree(grid, site))
188+
for (n, nb) in enumerate(neighbors(grid, site))
189+
if nb < site
190+
hopping_constants[species, site][n] = 1.0
191+
end
192+
end
193+
end
194+
```
195+
196+
To pass in `hopping_constants` of form $D_s * L_{i,j}$ we need two vectors -- one for $D_s$ and one for $L_{i,j}$. Here is an example.
197+
198+
```julia
199+
species_hop_constants = ones(num_species)
200+
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
201+
for site in 1:num_nodes
202+
site_hop_constants[site] = ones(outdegree(grids[1], site))
203+
end
204+
hopping_constants=Pair(species_hop_constants, site_hop_constants)
205+
```
206+
207+
We must combine both vectors into a pair as in the last line above.
208+
209+
Finally, to use in `hopping_constants` of form $D_{s,i} * L_{i,j}$ we construct a matrix instead of a vector for $D_{s,j}$.
210+
211+
```julia
212+
species_hop_constants = ones(num_species, num_nodes)
213+
site_hop_constants = Vector{Vector{Float64}}(undef, num_nodes)
214+
for site in 1:num_nodes
215+
site_hop_constants[site] = ones(outdegree(grids[1], site))
216+
end
217+
hopping_constants=Pair(species_hop_constants, site_hop_constants)
218+
```
219+
220+
We can use either of the four versions of `hopping_constants` to construct a `JumpProblem` with the same syntax as in the original example. The different forms of hopping rates are supported not only for convenience but also for better memory usage and performance. So it is recommended that the most specialized form of hopping rates is used.
221+
222+
### Solvers
223+
224+
There are currently two specialized "spatial" solvers: `NSM` and `DirectCRDirect`. The former stands for Next Subvolume Method \citep{elf2004spontaneous}. The latter employs Composition-Rejection to sample the next site to fire, similar to the ordinary DirectCR method. For larger networks `DirectCRDirect` is expected to be faster. Both methods can be used interchangeably.
225+
226+
Additionally, all standard solvers are supported as well, although they are expected to use more memory and be slower. They "flatten" the problem, i.e. turn all hops into reactions, resulting in a much larger system. For example, to use the Next Reaction Method (`NRM`), simply pass in `NRM()` instead of `NSM()` in the construction of the `JumpProblem`. Importantly, you *must* pass in `hopping_constants` in the `D_{s,i,j}` or `D_{s,i}` form to use any of the non-specialized solvers.
227+
228+
## References
229+
* \biblabel{elf2004spontaneous}{Elf and Ehrenberg (2004)} Elf, Johan and Ehrenberg, Mäns. “Spontaneous separation of bi-stable biochemical systems into spatial domains of opposite phases”. In: _Systems biology_ 1.2 (2004), pp. 230–236.

0 commit comments

Comments
 (0)