2. Rock Paper Scissors Stochastic Reaction Diffusion

Author

Nathan Zimmerberg

Published

February 13, 2023

2. Rock Paper Scissors Stochastic Reaction Diffusion

Inspired by:

https://www.youtube.com/watch?v=TORwMc2AaRE&t=0s

Setup

Open a julia REPL or jupyter notebook using the medyan-tutorial environment you created in tutorial 1.

Load MEDYAN with MEDYANVis for visualization

using MEDYAN
using MEDYANVis
using CairoMakie
using Random
using SmallZarrGroups
Random.seed!(1234);

Defining System

The first step is to create a system definition (MEDYAN.SysDef).

This system will contain three diffusing species: Rock, Paper, and Scissors.

We create the species using add_diffusing_species!, which takes the species name and diffusion coefficient.

All units are based on nm, s, pN.

s = MEDYAN.SysDef()
D = 1E6 # nm²/s
add_diffusing_species!(s, :R; coeff=D) # Rock
add_diffusing_species!(s, :P; coeff=D) # Paper
add_diffusing_species!(s, :S; coeff=D) # Scissors

MEDYAN.SysDef
 Diffusing species:
  :R: 1.0e6 nm²/s
  :P: 1.0e6 nm²/s
  :S: 1.0e6 nm²/s

Adding Reactions

Paper “beats” rock, scissors “beats” paper, rock “beats” scissors. In this simulation “beats” means convert to self type.

In this case because the reaction is between two diffusing species, it has units of nm³/s. You can think of this as the volume where if exactly one of each reactant existed the rate would be 1 per second.

addreaction!

rate = 50.0*(500.0^3) # nm³/s
addreaction!(s,
    "diffusing.P + diffusing.R --> 2diffusing.P",
    rate,
    1,
)
addreaction!(s,
    "diffusing.S + diffusing.P --> 2diffusing.S",
    rate,
    1,
)
addreaction!(s,
    "diffusing.R + diffusing.S --> 2diffusing.R",
    rate,
    1,
)

MEDYAN.SysDef
 Diffusing species:
  :R: 1.0e6 nm²/s
  :P: 1.0e6 nm²/s
  :S: 1.0e6 nm²/s
 Chem voxel reactions without callbacks:
  "diffusing.P + diffusing.R --> 2diffusing.P" 6.25e9 nm³/s
  "diffusing.S + diffusing.P --> 2diffusing.S" 6.25e9 nm³/s
  "diffusing.R + diffusing.S --> 2diffusing.R" 6.25e9 nm³/s

Creating a grid

All simulations run in a grid.

This is a box that contains the simulation.

The grid is composed of chemistry voxels, small cubes that can have different counts of the various diffusing species.

Create a 50 by 50 by 1 grid of 500 nm side length voxels.

L = 50
grid = CubicGrid((L,L,1),500.0)

CubicGrid([50, 50, 1], 500.0)

Creating a Context

The Context is the object the handles the state of a running simulation.

Create a Context using s and grid

c = MEDYAN.Context(s, grid)

MEDYAN.Context at time 0.0s in MEDYAN.CubicGrid([50, 50, 1], 500.0)

Add agents to the context

The context starts empty.

Distribute 2000 of each species randomly to the chem voxels by mutating the context.

adddiffusingcount_rand!

adddiffusingcount_rand!(c, s.diffusing.R, 2000)
adddiffusingcount_rand!(c, s.diffusing.P, 2000)
adddiffusingcount_rand!(c, s.diffusing.S, 2000)

Run chemistry

Run chemistry for 1.0s and visualize the results.

vis = Visualizer()
setvisible!(vis["/Grid"], false)
setvisible!(vis["/Axes"], false)
setvisible!(vis["/Background"], false)

If you open the visualizer in a browser with open(vis) and run the following you should get an animated version of below.

for i in 1:100
    run_chemistry!(c, 0.01)
    MEDYANVis.drawdiffusing!(vis["diffusing"], c.grid, c.chemistryengine, s; size = 0.8)
    sleep(0.02)
end

Snapshots and plotting

Next we will see how to plot quantities from the simulation.

If you want to analyze a trajectory in detail, or load a system state, one method is to use snapshots. For details see snapshot group reference

MEDYAN.load_snapshot! can be used to reload snapshots, this is very useful for restarting simulations, or changing the simulation state in a different programming language, but not needed in this tutorial.

SmallZarrGroups.jl can be used to save/load snapshot groups to/from disk, in zip files, or in directories.

Run chemistry for another 1.0s and store a snapshot every 0.01s in a vector.

using SmallZarrGroups
snapshots = ZGroup[]
for i in 1:100
    run_chemistry!(c, 0.01)
    push!(snapshots, MEDYAN.snapshot(c))
end

Now that the simulation is done the results can be analyzed. For longer simulations, these snapshots should be saved to disk with for example:

SmallZarrGroups.save_dir("snapshot$i.zarr.zip",MEDYAN.snapshot(c))

to avoid needed to rerun the simulation and to avoid using too much RAM.

total_rocks = Int[]
total_papers = Int[]
total_scissors = Int[]
times = Float64[]
for group in snapshots
    totals = sum(collect(group["chem/dc"]);dims=(2,3,4))
    push!(times, attrs(group)["time (s)"])
    push!(total_rocks, totals[s.diffusing.R])
    push!(total_papers, totals[s.diffusing.P])
    push!(total_scissors, totals[s.diffusing.S])
end


figure, axis, lineplot = lines(times, total_rocks; label="rock")
lines!(times, total_papers; label="paper")
lines!(times, total_scissors; label="scissors")
lines!(
    times,
    total_scissors.+total_papers.+total_rocks;
    label="combined")
axislegend()
current_figure()

If you are using VSCode or Jupyter the figure will be displayed.

If you are using the REPL, save the figure as a PDF to view it.

save("tutorial2-plot.pdf", current_figure())