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

Declaring agent names

The first step is declaring the names of agents in the system.

This system will contain three diffusing species.

agent_names = MEDYAN.AgentNames(;
    diffusingspeciesnames=[
        :R, # Rock
        :P, # Paper
        :S, # Scissors
    ],
)
MEDYAN.AgentNames([:R, :P, :S], Symbol[], Symbol[], Symbol[], Tuple{Symbol, Vector{Symbol}}[], Symbol[], Symbol[], Symbol[])

Defining System

After the agent names are declared, the system can be defined.

This is done by constructing a MEDYAN.SysDef object from the agent_names and then mutating it to add parameters, reactions, and callbacks.

s = MEDYAN.SysDef(agent_names)
MEDYAN.SysDef
 Diffusing species:
  R: 0.0 nm²/s
  P: 0.0 nm²/s
  S: 0.0 nm²/s

Add a diffusion coefficient of 1E6 nm²/s to all species. All units are based on nm, s, pN

add_diffusion_coeff!

D = 1E6
add_diffusion_coeff!(s, :R, D) # nm²/s
add_diffusion_coeff!(s, :P, D) # nm²/s
add_diffusion_coeff!(s, :S, D) # nm²/s
MEDYAN.SysDef
 Diffusing species:
  R: 1.0e6 nm²/s
  P: 1.0e6 nm²/s
  S: 1.0e6 nm²/s

Add 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
 Compartment 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 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["diffusingcounts"]);dims=2)
    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()
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
└ @ Makie ~/.julia/packages/Makie/aG6xp/src/scenes.jl:227

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