GitHub

API for ReactionDiffusion.jl

ReactionDiffusion aims to be an easy-to-use and computationally-efficient pipeline to simulate biologically-inspired reaction-diffusion models. It is our hope that models can be built with just a few lines of code and solved without the user having any knowledge of PDE solver methods.

This is achieved by drawing from a range of numerical routines from the SciML packages, including Catalyst.jl, Symbolics.jl, ModelingToolkit.jl, and DifferentialEquations.jl.

For users more familiar with PDE solvers, it is possible to specify optional arguments to e.g., control the specific solver algorithm used (see the API below). However, ReactionDiffusion does not aim to be a fully customizable PDE solving package that covers all bases; our focus is to make something that is relatively easy-to-use and performant but only for a restricted number of use cases.

If you require more customization or a more flexible PDE-solving framework, we highly recommend ModelingToolkit.jl and/or Catalyst.jl.

Data structures

We introduce two structures

  • save_turing is an object that records parameter sets that undergo a Turing instability. It has the following fields:

    • steady_state_values: The computed steady state values of each variable
    • reaction_params: The reaction parameters
    • diffusion_constants: The diffusion constants
    • initial_conditions: The initial conditions of each variable used to compute the steady state values
    • pattern_phase: The predicted phase of each of the variables in the final pattern.[1 1]would be in-phase,[1 -1]` would be out-of-phase
    • wavelength: The wavelength that is maximally unstable
    • max_real_eigval: The maximum real eigenvalue associated with the Turing instability
    • non_oscillatory: If true, this parameter set represents a stationary Turing pattern. If false, the unstable mode has complex eigenvalues and thus may be oscillatory.

  • model_parameters is an object that records parameter values for model simulations. It has the following fields:

    -reaction: reaction parameters -diffusion: diffusion constants -initial_condition: values of each of the variables used for homogeneous initial conditions -initial_noise: level of (normally distributed) noise added to the homogeneous initial conditions -domain_size: size of 1D domain -random_seed: seed associated with random number generation; only set this when you need to reproduce exact simulations each run

Functions

ReactionDiffusion.screen_values — Function
screen_values(;min=0,max=1,number=10, mode="linear")

Returns a series of num_params parameter values that are linearly spaced between the min and max limits. The argument mode="log" may be used to sample in log-space instead.

source
ReactionDiffusion.get_params — Function
get_params(model, turing_param)

For a given model and a single pattern-forming parameter set, turing_param, this function creates a corresponding model_parameters variable. This sets, by default, the model_parameters fields:

  • domain_size: chosen to be 3x the computed pattern wavelength, i.e., 3*turing_param.wavelength
  • initial_condition: chosen to be the computed steady state values, i.e., turing_param.steady_state_values
  • initial_noise = 0.01: the magnitude of noise (normally distributed random numbers) added to the steady state values to define the initial conditions.
source
ReactionDiffusion.get_param — Function
get_param(model, turing_params, name, type)

For a given model and a (potentially large) number of pattern-forming parameter sets, turing_params, this function extracts the parameter values prescribed by the input name. For reaction parameters, used type="reaction", for diffusion constants, use type="diffusion".

Example:

δ₁ = get_param(model, turing_params,"δ₁","reaction")
D_COMPLEX = get_param(model, turing_params,"COMPLEX","diffusion")
source
ReactionDiffusion.returnTuringParams — Function
returnTuringParams(model, params; maxiters = 1e3,alg=Rodas5(),abstol=1e-8, reltol=1e-6, tspan=1e4,ensemblealg=EnsembleThreads(),batch_size=1e4)

Return a save_turing object of parameters that are predicted to be pattern forming.

Required inputs:

  • model: specified via the @reaction_network macro
  • params: all reaction and diffusion parameters, in a model_parameters object

Optional inputs:

  • batch_size: the number of parameter sets to consider at once. Increasing/decreasing from the default value may improve speed.

Inputs carried over from DifferentialEquations.jl; see here for further details:

  • maxiters: maximum number of iterations to reach steady state (otherwise simulation terminates)
  • alg: ODE solver algorithm
  • abstol and reltol: tolerance levels of solvers
  • tspan: maximum time allowed to reach steady state (otherwise simulation terminates)
  • ensemblealg: ensemble simulation method
source
ReactionDiffusion.simulate — Function
simulate(model,param;alg=KenCarp4(),reltol=1e-6,abstol=1e-8, dt = 0.1, maxiters = 1e3, save_everystep = true)

Simulate model for a single parameter set param.

Required inputs:

  • model: specified via the @reaction_network macro
  • param: all reaction and diffusion parameters, in a model_parameters object. This must be a single parameter set only

Inputs carried over from DifferentialEquations.jl; see here for further details:

  • maxiters: maximum number of iterations to reach steady state (otherwise simulation terminates)
  • alg: solver algorithm
  • abstol and reltol: tolerance levels of solvers
  • dt: initial value for timestep
  • save_everystep: controls whether all timepoints are saved, defaults to true
source