ReactionDiffusion.jl for modelling pattern formation in biological systems
Reaction-diffusion dynamics are present across many areas of the physical and natural world, and allow complex spatiotemporal patterns to self-organize de novo. ReactionDiffusion.jl aims to be an easy-to-use and performant pipeline to simulate reaction-diffusion PDEs of arbitrary complexity, with a focus on pattern formation in biological systems. Using this package, complex, biologically-inspired reaction-diffusion models can be:
- specified using an intuitive syntax
- screened across millions of parameter combinations to identify pattern-forming networks (i.e., those that undergo a Turing instability)
- rapidly simulated to predict spatiotemporal patterns
Quick start guide
Here we show how ReactionDiffusion.jl can be used to simulate a biologically-inspired reaction-diffusion system, responsible for generating evenly spaced joints along the length of your fingers and toes (from Grall et el, 2024).
We begin by specifying the reaction-diffusion dynamics via the intuitive syntax developed in Catalyst.jl, which naturally mirrors biochemical feedbacks and interactions.
using ReactionDiffusion
model = @reaction_network begin
# complex formation
(k₊, k₋), GDF5 + NOG <--> COMPLEX
# degradation
δ₁, GDF5 --> ∅
δ₂, NOG --> ∅
δ₃, pSMAD --> ∅
# transcriptional feedbacks (here: repressive hill functions)
hillr(pSMAD,μ₁,K₁,n₁), ∅ --> GDF5
hillr(pSMAD,μ₂,K₂,n₂), ∅ --> NOG
# signalling
μ₃*GDF5, ∅ --> pSMAD
end\[ \begin{align*} \mathrm{\mathtt{GDF5}} + \mathrm{\mathtt{NOG}} &\xrightleftharpoons[\mathtt{k_-}]{\mathtt{k.}} \mathrm{\mathtt{COMPLEX}} \\ \mathrm{\mathtt{GDF5}} &\xrightarrow{\mathtt{\delta_1}} \varnothing \\ \mathrm{\mathtt{NOG}} &\xrightarrow{\mathtt{\delta_2}} \varnothing \\ \mathrm{\mathtt{pSMAD}} &\xrightarrow{\mathtt{\delta_3}} \varnothing \\ \varnothing &\xrightarrow{\frac{\mathtt{K_1}^{\mathtt{n_1}} \mathtt{\mu_1}}{\mathtt{pSMAD}^{\mathtt{n_1}} + \mathtt{K_1}^{\mathtt{n_1}}}} \mathrm{\mathtt{GDF5}} \\ \varnothing &\xrightarrow{\frac{\mathtt{K_2}^{\mathtt{n_2}} \mathtt{\mu_2}}{\mathtt{K_2}^{\mathtt{n_2}} + \mathtt{pSMAD}^{\mathtt{n_2}}}} \mathrm{\mathtt{NOG}} \\ \varnothing &\xrightarrow{\mathtt{GDF5} \mathtt{\mu_3}} \mathrm{\mathtt{pSMAD}} \end{align*} \]
We can then specify values for each parameter:
params = model_parameters()
# constant parameters
params.reaction["δ₃"] = [1.0]
params.reaction["μ₁"] = [1.0]
params.reaction["μ₃"] = [1.0]
params.reaction["n₁"] = [8]
params.reaction["n₂"] = [2]
# varying parameters
num_params = 5
params.reaction["δ₁"] = screen_values(min = 0.1,max = 10, number=num_params)
params.reaction["δ₂"] = screen_values(min = 0.1,max = 10,number=num_params)
params.reaction["μ₂"] = screen_values(min = 0.1,max = 10, number=num_params)
params.reaction["k₊"] = screen_values(min = 10.0,max = 100.0, number=num_params)
params.reaction["k₋"] = screen_values(min = 10.0,max = 100.0, number=num_params)
params.reaction["K₁"] = screen_values(min = 0.01,max = 1,number=num_params)
params.reaction["K₂"] = screen_values(min = 0.01,max = 1, number=num_params)
params.reactionDict{Any, Any} with 12 entries:
"μ₁" => [1.0]
"μ₂" => [0.1, 2.575, 5.05, 7.525, 10.0]
"K₂" => [0.01, 0.2575, 0.505, 0.7525, 1.0]
"k₋" => [10.0, 32.5, 55.0, 77.5, 100.0]
"δ₂" => [0.1, 2.575, 5.05, 7.525, 10.0]
"n₁" => [8]
"n₂" => [2]
"K₁" => [0.01, 0.2575, 0.505, 0.7525, 1.0]
"μ₃" => [1.0]
"δ₁" => [0.1, 2.575, 5.05, 7.525, 10.0]
"k₊" => [10.0, 32.5, 55.0, 77.5, 100.0]
"δ₃" => [1.0]We must then also specify the diffusion coefficients for each variable:
params.diffusion = Dict(
"NOG" => [1.0],
"GDF5" => screen_values(min = 0.1,max = 30, number=10),
"COMPLEX" => screen_values(min = 0.1,max = 30, number=10)
)
params.diffusionDict{String, Vector{Float64}} with 3 entries:
"COMPLEX" => [0.1, 3.42222, 6.74444, 10.0667, 13.3889, 16.7111, 20.0333, 23.3…
"NOG" => [1.0]
"GDF5" => [0.1, 3.42222, 6.74444, 10.0667, 13.3889, 16.7111, 20.0333, 23.3…Then, with a single line of code, we can perform a Turing instability analysis across all combinations of parameters:
julia> turing_params = returnTuringParams(model, params);Screening parameter sets: 25%[============> ] ETA: 0:03:28 Screening parameter sets: 38%[==================> ] ETA: 0:03:00 Screening parameter sets: 50%[=========================> ] ETA: 0:02:10 Screening parameter sets: 62%[===============================> ] ETA: 0:01:38 Screening parameter sets: 75%[=====================================> ] ETA: 0:01:00 Screening parameter sets: 88%[===========================================> ] ETA: 0:00:29 Screening parameter sets: 100%[==================================================] Time: 0:03:38 3405/7812500 parameters are pattern forming
This returns all parameter combinations that can break symmetry from a homogeneous initial condition. We take advantage of the highly performant numerical solvers in DifferentialEquations.jl to be able to simulate millions of parameter sets per minute on a standard laptop.
We may then take a single parameter set and simulate its spatiotemporal dynamics directly, using Plots.jl to visualize the resulting pattern:
param1 = get_params(model, turing_params[1000])
sol = simulate(model,param1)
using Plots
plot(endpoint(),model,sol)
We may also view the full spatiotemporal dynamics:
@gif for t in 0.0:0.01:1
plot(timepoint(),model,sol,t)
end fps=20
Support, citation and future developments
If you find ReactionDiffusion.jl helpful in your research, teaching, or other activities, please star the repository and consider citing this paper:
@article{TBD,
doi = {TBD},
author = {Muzatko, Daniel AND Daga, Bijoy AND Hiscock, Tom W.},
journal = {biorXiv},
publisher = {TBD},
title = {TBD},
year = {TBD},
month = {TBD},
volume = {TBD},
url = {TBD},
pages = {TBD},
number = {TBD},
}We are a small team of academic researchers from the Hiscock Lab, who build mathematical models of developing embryos and tissues. We have found these scripts helpful in our own research, and make them available in case you find them helpful in your research too. We hope to extend the functionality of ReactionDiffusion.jl as our future projects, funding and time allows.
This work is supported by ERC grant SELFORG-101161207, and UK Research and Innovation (Biotechnology and Biological Sciences Research Council, grant number BB/W003619/1)
Funded by the European Union. Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union or the European Research Council Executive Agency. Neither the European Union nor the granting authority can be held responsible for them
