Model(reaction, diffusion, boundary_conditions, initial_conditions)

An object containing a mathematical description of a reaction diffusion system to be simulated, independent of parameter values.

Fields

• reaction::ReactionSystem
• diffusion::DiffusionSystem
• boundary_flux::(ReactionSystem, ReactionSystem)
• initial_conditions::SpeciesValues

@reaction_network

Macro for generating chemical reaction network models (Catalyst ReactionSystems). See the (DSL introduction and advantage usage) sections of the Catalyst documentation for more details on the domain-specific language (DSL) that the macro implements. The macro's output (a ReactionSystem structure) is central to Catalyst and its functionality. How to e.g. simulate these is described in the Catalyst documentation.

Returns:

• A Catalyst ReactionSystem, i.e. a symbolic model for the reaction network. The returned

system is marked complete. To obtain a ReactionSystem that is not marked complete, for example to then use in compositional modelling, see the otherwise equivalent @network_component macro.

Examples: Here we create a basic SIR model. It contains two reactions (infection and recovery):

sir_model = @reaction_network begin
    c1, S + I --> 2I
    c2, I --> R
end

Next, we create a self-activation loop. Here, a single component (X) activates its own production with a Michaelis-Menten function:

sa_loop = @reaction_network begin
    mm(X,v,K), 0 --> X
    d, X --> 0
end

This model also contains production and degradation reactions, where 0 denotes that there are either no substrates or no products in a reaction.

Options: In addition to reactions, the macro also supports "option" inputs (permitting e.g. the addition of observables). Each option is designated by a tag starting with a @ followed by its input. A list of options can be found here.

@diffusion_system L begin D, species;... end

Define a spatial domain of length L and a set of diffusion rates and Neumann boundary conditions for the given species. The boundary conditions uₓ(0)=a and uₓ(L)=b default to 0 if ommitted.

Example

@diffusion_system L begin
    0.5,             U
    Dᵥ, (0.0, 0.5),  V
    Dᵣ/k, (a, a*s),  R
end

@initial_conditions begin IC, species;... end

Define a set of initial conditions for the given species. IC may depend on arbitrary parameters and additionally the spatial variable x.

Example

@inital_conditions begin
    U0,             U
    V0 + exp(x),    V
end

function parameter_set(model, params)

Create a set of parameter values and initial conditions for model. Defaults are used for values missing from params.

simulate(model, params; output_func=nothing, full_solution=false, alg=ETDRK4(), num_verts=64, dt=0.1, max_attempts = 4, tol=1e-4, kwargs...)

Simulate model for the parameters and initial conditions given in params, stopping when a steady state is reached. Returns (u,t) with the solution values and time.

Arguments

• model: Model object containing the system to be simulated.
• params: Either a single parameter set or a vector of parameter sets to be solved as an ensemble. Parameter sets can be created manually with parameter_set or supplied as a dict or collection of pairs in which case defaults will be used for any missed values and low-level noise added to initial conditions. Parameters values may be either single numbers which are replicated homogenously over the domain, or functions mapping the interval [0.0,1.0] to values for the corresponding point in space.
• output_func(u, t): Function to transform output values.
• full_solution: Return a vector of values at each time point if true, instead of just the steady-state solution.
• max_attempts: Number of times to retry with reduced dt before giving up if the solution fails to converge.
• num_verts: Number of points in spatial discretisation.

For other keyword arguments see https://docs.sciml.ai/DiffEqDocs/stable/basics/commonsolveropts/.

function simulate(model; output_func=nothing, full_solution=false, alg=ETDRK4(), num_verts=64, dt=0.1, max_attempts = 4, tol=1e-5, noise=1e-4, kwargs...)

Partially applied version of simulate to avoid repeating expensive setup when simulating the same model reapeatedly.

turing_wavelength(model, params; k=logrange(0.01,1000,1000))

Compute dominant wavelengths of Turing instabilities for each of params. Returns 0.0 for parameter sets for which Turing instability does not occur.

is_turing(model,params)

Test whether params exhibit Turing instability.

filter_turing(model,params)

Return only params which demonstrate Turing instability. Multithreaded.

function steady_state_plot(model, params; normalise=true, hide_y=normalise, kwargs...)

Run a simulation and plot the final values for each species. If normalise is true, values for different species will be normalised to a common scale.

function steady_state_plot(model, sol; normalise=true, hide_y=normalise, kwargs...)

Plot the final value in the solution object for each species. If normalise is true, values for different species will be normalised to a common scale.

timeseries_plot(model, params; normalise=true, hide_y=normalise, autolimits=true, kwargs...)

Simulate and display the results with an interactive slider to move through time. The remaining kwargs are passed to simulate. If normalise is true, values for different species will be normalised to a common scale.

timeseries_plot(model, sol; normalise=true, hide_y=normalise, autolimits=true, kwargs...)

Display the solution with an interactive slider to move through time. If normalise is true, values for different species will be normalised to a common scale.

timeseries_plot(model, u, t; normalise=true, hide_y=normalise, autolimits=true, kwargs...)

Display the solution with an interactive slider to move through time. If normalise is true, values for different species will be normalised to a common scale.

interactive_plot(model, param_ranges; normalise=false, hide_y=normalise, autolimits=true, num_verts=32, kwargs...)

Generate an interactive plot of the steady state solution with sliders to adjust each of the parameters within param_ranges. param_ranges should be a dictionary mapping parameter names to either Range objects or collections of possible values.

Construct a SplitODEProblem to solve a reaction diffusion system with reflective boundaries.

Returns the SplitODEProblem with solutions in the frequency (DCT-1) domain and a FFTW plan to transform solutions back to the spatial domain.

PseudoSpectralProblem(species, reaction_rates, diffusion_rates, boundary_conditions, initial_conditions, num_verts; p=nothing, noise=1e-4, rng=default_rng(), kwargs...)

Construct a PsuedoSpectralProblem object representing a reaction diffusion system of the form uₓₓ(x,t) = Duₜ(x,t) + f(u(x,t)).

Arguments

PseudoSpectral expects Symbolics.jl expressions as inputs. The special variable 'x' ∈ [0,1] represents the spatial coordinate. Any variables other than 'x' and those supplied in species will be interpreted as parameters.

• species: Vector of variables corresponding to the components of u.
• reaction_rates: Vector of expressions representing f(u).
• diffusion_rates: Vector of Expressions representing diag(D).
• boundary_conditions: 2xn matrix of expressions representing Neumann boundary conditions. The two rows correspond to uₓ at the left and right boundaries.
• initial_conditions: Vector of expressions representing u(x,0).
• num_verts: Number of points in spatial discretisation.
• p=nothing: Dictionary associating parameters with numerical values.
• noise=1e-4: Guassian noise with σ²=noise is added to the initial conditions.
• rng=default_rng(): Random number generator for noise.
• kwargs...: Keyword arguments passed on to SciML's solve. For details see https://docs.sciml.ai/DiffEqDocs/stable/basics/commonsolveropts/.

PseudoSpectralSolution

Solution object for PsuedoSpectralProblem.

Indexing

• By time-step sol[3].
• By species sol[U].

EnsembleProblem(prob::PseudoSpectralProblem, params; output_func=nothing)

Construct an ensemble problem to solve the system in parallel for each of the supplied parameter sets.

EnsembleProblem(prob::PseudoSpectralProblem; prob_func, output_func=nothing)

Construct an ensemble problem to run the solver in parallel. For details see https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/.

steady_state_callback(tol=1e-4)

Callback function to be passed to solve to detect steady state. Terminates solver when |uₜ| ≤ tol.

solve(prob::PseudoSpectralProblem, alg=ETDRK4(); kwargs...)

See https://docs.sciml.ai/DiffEqDocs/stable/basics/commonsolveropts/. Algorithm defaults to ETDRK4.

PseudoSpectralIntegrator(prob::PseudoSpectralProblem; alg=ETDRK4(), kwargs...)

Initialize an integrator for the problem. See https://docs.sciml.ai/DiffEqDocs/stable/basics/integrator/. Algorithm defaults to ETDRK4.

remake(prob::PseudoSpectralProblem; p=nothing, rng=nothing, kwargs...)

Return a new problem with updated parameters, random number generator, and/or solver options.

remake(integrator::PseudoSpectralIntegrator; kwargs...)

Return a new integrator updated with remake(integrator.prob; kwargs...).

get_u(integrator::PseudoSpectralIntegrator)

Return solution values at time t, stepping the integrator as necessary.

get_sol(integrator::PseudoSpectralIntegrator)

Return a solution object for the current integrator state.

step!(integrator::PseudoSpectralIntegrator, dt=nothing, stop_at_tdt=false)

Advance the iterator by dt.

step_to!(integrator::PseudoSpectralIntegrator, t, stop_at_tdt=false)

Advance the iterator to time t.

Spatial variable x∈[0,1].