Catalyst.jl API

@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.

@network_component

Equivalent to @reaction_network except the generated ReactionSystem is not marked as complete.

make_empty_network(; iv=DEFAULT_IV, name=gensym(:ReactionSystem))

Construct an empty ReactionSystem. iv is the independent variable, usually time, and name is the name to give the ReactionSystem.

@reaction

Macro for generating a single Reaction object using a similar syntax as the @reaction_network macro (but permitting only a single reaction). A more detailed introduction to the syntax can be found in the description of @reaction_network.

The @reaction macro is followed by a single line consisting of three parts:

• A rate (at which the reaction occurs).
• Any number of substrates (which are consumed by the reaction).
• Any number of products (which are produced by the reaction).

The output is a reaction (just like created using the Reaction constructor).

Examples: Here we create a simple binding reaction and store it in the variable rx:

rx = @reaction k, X + Y --> XY

The macro will automatically deduce X, Y, and XY to be species (as these occur as reactants) and k as a parameter (as it does not occur as a reactant).

The @reaction macro provides a more concise notation to the Reaction constructor. I.e. here we create the same reaction using both approaches, and also confirm that they are identical.

# Creates a reaction using the `@reaction` macro.
rx = @reaction k*v, A + B --> C + D

# Creates a reaction using the `Reaction` constructor.
t = default_t()
@parameters k v
@species A(t) B(t) C(t) D(t)
rx2 = Reaction(k*v, [A, B], [C, D])

# Confirms that the two approaches yield identical results:
rx1 == rx2

Interpolation of already declared symbolic variables into @reaction is possible:

t = default_t()
@parameters k b
@species A(t)
ex = k*A^2 + t
rx = @reaction b*$ex*$A, $A --> C

Notes:

• @reaction does not support bi-directional type reactions (using <-->) or reaction bundling

(e.g. d, (X,Y) --> 0).

• Interpolation of Julia variables into the macro works similarly to the @reaction_network

macro. See The Reaction DSL tutorial for more details.

struct Reaction{S, T}

One chemical reaction.

Fields

• rate: The rate function (excluding mass action terms).

• substrates: Reaction substrates.

• products: Reaction products.

• substoich: The stoichiometric coefficients of the reactants.

• prodstoich: The stoichiometric coefficients of the products.

• netstoich: The net stoichiometric coefficients of all species changed by the reaction.

• only_use_rate: false (default) if rate should be multiplied by mass action terms to give the rate law. true if rate represents the full reaction rate law.

• metadata: Contain additional data, such whenever the reaction have a specific noise-scaling expression for the chemical Langevin equation.

Examples

using Catalyst
t = default_t()
@parameters k[1:20]
@species A(t) B(t) C(t) D(t)
rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
       Reaction(k[2], [B], nothing),            # B -> 0
       Reaction(k[3],[A],[C]),                  # A -> C
       Reaction(k[4], [C], [A,B]),              # C -> A + B
       Reaction(k[5], [C], [A], [1], [2]),      # C -> A + A
       Reaction(k[6], [A,B], [C]),              # A + B -> C
       Reaction(k[7], [B], [A], [2], [1]),      # 2B -> A
       Reaction(k[8], [A,B], [A,C]),            # A + B -> A + C
       Reaction(k[9], [A,B], [C,D]),            # A + B -> C + D
       Reaction(k[10], [A], [C,D], [2], [1,1]), # 2A -> C + D
       Reaction(k[11], [A], [A,B], [2], [1,1]), # 2A -> A + B
       Reaction(k[12], [A,B,C], [C,D], [1,3,4], [2, 3]),          # A+3B+4C -> 2C + 3D
       Reaction(k[13], [A,B], nothing, [3,1], nothing),           # 3A+B -> 0
       Reaction(k[14], nothing, [A], nothing, [2]),               # 0 -> 2A
       Reaction(k[15]*A/(2+A), [A], nothing; only_use_rate=true), # A -> 0 with custom rate
       Reaction(k[16], [A], [B]; only_use_rate=true),             # A -> B with custom rate.
       Reaction(k[17]*A*exp(B), [C], [D], [2], [1]),              # 2C -> D with non constant rate.
       Reaction(k[18]*B, nothing, [B], nothing, [2]),             # 0 -> 2B with non constant rate.
       Reaction(k[19]*t, [A], [B]),                                # A -> B with non constant rate.
       Reaction(k[20]*t*A, [B,C], [D],[2,1],[2])                  # 2A +B -> 2C with non constant rate.
  ]

Notes:

• nothing can be used to indicate a reaction that has no reactants or no products. In this case the corresponding stoichiometry vector should also be set to nothing.
• The three-argument form assumes all reactant and product stoichiometric coefficients are one.
• Pass unit_checks = true to validate unit consistency at construction time. When enabled, checks that all substrates and products share the same units, and that the rate expression has internally consistent additive terms. Default is false.

struct ReactionSystem{V<:Catalyst.NetworkProperties} <: ModelingToolkitBase.AbstractSystem

A system of chemical reactions.

Fields

• eqs: The equations (reactions and algebraic/differential) defining the system.

• rxs: The Reactions defining the system.

• iv: Independent variable (usually time).

• sivs: Spatial independent variables

• unknowns: All dependent (unknown) variables, species and non-species. Must not contain the independent variable.

• species: Dependent unknown variables representing species

• ps: Parameter variables. Must not contain the independent variable.

• var_to_name: Maps Symbol to corresponding variable.

• observed: Equations for observed variables.

• name: The name of the system

• systems: Internal sub-systems

• initial_conditions: The initial values to use when initial conditions and/or parameters are not supplied to problem constructors.

• networkproperties: NetworkProperties object that can be filled in by API functions. INTERNAL – not considered part of the public API.

• combinatoric_ratelaws: Sets whether to use combinatoric scalings in rate laws. true by default.

• continuous_events: continuous_events: A Vector{SymbolicContinuousCallback} that model events. The integrator will use root finding to guarantee that it steps at each zero crossing.

• discrete_events: discrete_events: A Vector{SymbolicDiscreteCallback} that models events. Symbolic analog to SciMLBase.DiscreteCallback that executes an affect when a given condition is true at the end of an integration step.

• brownians: Brownian variables for non-reaction noise, created via @brownians.

• poissonians: Poissonian variables for Poisson jump noise, created via @poissonians.

• jumps: Non-reaction jumps (VariableRateJump, ConstantRateJump, MassActionJump).

• metadata: Metadata for the system, to be used by downstream packages.

• complete: complete: if a model sys is complete, then sys.x no longer performs namespacing.

• parent: The hierarchical parent system before simplification that MTK now seems to require for hierarchical namespacing to work in indexing.

Example

Continuing from the example in the Reaction definition:

# simple constructor that infers species and parameters
@named rs = ReactionSystem(rxs, t)

# allows specification of species and parameters
@named rs = ReactionSystem(rxs, t, [A,B,C,D], k)

Keyword Arguments:

• observed::Vector{Equation}, equations specifying observed variables.
• systems::Vector{ReactionSystem}, vector of sub-ReactionSystems.
• name::Symbol, the name of the system (must be provided, or @named must be used).
• initial_conditions::SymmapT, a dictionary mapping parameters and species to their initial values.
• checks = true, boolean for whether to run structural checks at construction time.
• unit_checks = false, boolean for whether to perform unit validation at construction time. Uses Catalyst.validate_units / Catalyst.assert_valid_units. Units should be specified with symbolic units (us"...") from DynamicQuantities.
• networkproperties = NetworkProperties(), cache for network properties calculated via API functions.
• combinatoric_ratelaws = true, sets the default value of combinatoric_ratelaws used in conversion functions (ode_model, sde_model, jump_model, ss_ode_model, hybrid_model) or when calling problem constructors with the ReactionSystem.
• balanced_bc_check = true, sets whether to check that BC species appearing in reactions are balanced (i.e appear as both a substrate and a product with the same stoichiometry).
• brownians, vector of Brownian variables for non-reaction SDE noise (created via @brownians). Auto-discovered from equations in the two-argument constructor.
• poissonians, vector of Poissonian variables for Poisson jump noise (created via @poissonians). Auto-discovered from equations in the two-argument constructor.
• jumps, vector of non-reaction jump processes.

species(network)

Given a ReactionSystem, return a vector of all species defined in the system and any subsystems that are of type ReactionSystem. To get the species and non-species variables in the system and all subsystems, including non-ReactionSystem subsystems, uses unknowns(network).

Notes:

• If ModelingToolkitBase.get_systems(network) is non-empty will allocate.

get_species(sys::ReactionSystem)

Return the current dependent variables that represent species in sys (toplevel system only).

nonspecies(network)

Return the non-species variables within the network, i.e. those unknowns for which isspecies == false.

Notes:

• Allocates a new array to store the non-species variables.

reactions(network)

Given a ReactionSystem, return a vector of all Reactions in the system.

Notes:

• If ModelingToolkitBase.get_systems(network) is not empty, will allocate.

get_rxs(sys::ReactionSystem)

Return the system's Reaction vector (toplevel system only).

nonreactions(network)

Return the non-reaction equations within the network (i.e. algebraic and differential equations).

Notes:

• Allocates a new array to store the non-species variables.

numspecies(network)

Return the total number of species within the given ReactionSystem and subsystems that are ReactionSystems.

numparams(network)

Return the total number of parameters within the given system and all subsystems.

numreactions(network)

Return the total number of reactions within the given ReactionSystem and subsystems that are ReactionSystems.

speciesmap(network)

Given a ReactionSystem, return a Dictionary mapping species that participate in Reactions to their index within species(network).

paramsmap(network)

Given a ReactionSystem, return a Dictionary mapping from all parameters that appear within the system to their index within parameters(network).

isautonomous(rs::ReactionSystem)

Checks if a system is autonomous (i.e. no rate or equation depend on the independent variable(s)). Example:

rs1 = @reaction_system
    (p,d), 0 <--> X
end
isautonomous(rs1) # Returns `true`.

rs2 = @reaction_system
    (p/t,d), 0 <--> X
end
isautonomous(rs2) # Returns `false`.

has_alg_equations(sys::AbstractSystem)

For a system, returns true if it contain at least one algebraic equation (i.e. that does not contain any differentials).

Example:

using ModelingToolkitBase
using ModelingToolkitBase: t_nounits as t, D_nounits as D
@parameters p d
@variables X(t)
eq1 = D(X) ~ p - d*X
eq2 = 0 ~ p - d*X
@named osys1 = System([eq1], t)
@named osys2 = System([eq2], t)

has_alg_equations(osys1) # returns `false`.
has_alg_equations(osys2) # returns `true`.

alg_equations(sys::AbstractSystem)

For a system, returns a vector of all its algebraic equations (i.e. that does not contain any differentials).

Example:

using ModelingToolkitBase
using ModelingToolkitBase: t_nounits as t, D_nounits as D
@parameters p d
@variables X(t)
eq1 = D(X) ~ p - d*X
eq2 = 0 ~ p - d*X
@named osys = System([eq1, eq2], t)

alg_equations(osys) # returns `[0 ~ p - d*X(t)]`.

has_diff_equations(sys::AbstractSystem)

For a system, returns true if it contain at least one differential equation (i.e. that contain a differential).

Example:

using ModelingToolkitBase
using ModelingToolkitBase: t_nounits as t, D_nounits as D
@parameters p d
@variables X(t)
eq1 = D(X) ~ p - d*X
eq2 = 0 ~ p - d*X
@named osys1 = System([eq1], t)
@named osys2 = System([eq2], t)

has_diff_equations(osys1) # returns `true`.
has_diff_equations(osys2) # returns `false`.

diff_equations(sys::AbstractSystem)

For a system, returns a vector of all its differential equations (i.e. that does contain a differential).

Example:

using ModelingToolkitBase
using ModelingToolkitBase: t_nounits as t, D_nounits as D
@parameters p d
@variables X(t)
eq1 = D(X) ~ p - d*X
eq2 = 0 ~ p - d*X
@named osys = System([eq1, eq2], t)

diff_equations(osys) # returns `[Differential(t)(X(t)) ~ p - d*X(t)]`.

isspecies(s)

Tests if the given symbolic variable corresponds to a chemical species.

Catalyst.isconstant(s)

Tests if the given symbolic variable corresponds to a constant species.

Catalyst.isbc(s)

Tests if the given symbolic variable corresponds to a boundary condition species.

isvalidreactant(s)

Test if a species is valid as a reactant (i.e. a species variable or a constant parameter).

getname(x)::Symbol

Get the name of a symbolic variable as a Symbol. Acts as the identity function for x::Symbol.

ismassaction(rx, rs; rxvars = get_variables(rx.rate),
                              haveivdep = nothing,
                              unknownset = Set(unknowns(rs)),
                              ivset = nothing)

True if a given reaction is of mass action form, i.e. rx.rate does not depend on any chemical species that correspond to unknowns of the system, and does not depend explicitly on the independent variable (usually time).

Arguments

• rx, the Reaction.
• rs, a ReactionSystem containing the reaction.
• Optional: rxvars, Variables which are not in rxvars are ignored as possible dependencies.
• Optional: haveivdep, true if the Reactionrate field explicitly depends on any independent variable (i.e. t or for spatial systems x,y,etc). If not set, will be automatically calculated.
• Optional: unknownset, set of unknowns which if the rxvars are within mean rx is non-mass action.
• Optional: ivset, a Set of the independent variables of the system. If not provided and the system is spatial, i.e. isspatial(rs) == true, it will be created with all the spatial variables and the time variable. If the rate expression contains any element of ivset, then ismassaction(rx,rs) == false. Pass a custom set to control this behavior.

Notes:

• Non-integer stoichiometry is treated as non-mass action. This includes symbolic variables/terms or floating point numbers for stoichiometric coefficients.

dependents(rx, network)

Given a Reaction and a ReactionSystem, return a vector of the non-constant species and variables the reaction rate law depends on. e.g., for

k*W, 2X + 3Y --> 5Z + W

the returned vector would be [W(t),X(t),Y(t)].

Notes:

• Allocates
• Does not check for dependents within any subsystems.
• Constant species are not considered dependents since they are internally treated as parameters.
• If the rate expression depends on a non-species unknown variable that will be included in the dependents, i.e. in

  t = default_t()
  @parameters k
  @variables V(t)
  @species A(t) B(t) C(t)
  rx = Reaction(k*V, [A, B], [C])
  @named rs = ReactionSystem([rx], t)
  issetequal(dependents(rx, rs), [A,B,V]) == true

dependents(rx, network)

See documentation for dependents.

substoichmat(rn; sparse=false)

Returns the substrate stoichiometry matrix, $S$, with $S_{i j}$ the stoichiometric coefficient of the ith substrate within the jth reaction.

Note:

• Set sparse=true for a sparse matrix representation
• Note that constant species are not considered substrates, but just components that modify the associated rate law.

prodstoichmat(rn; sparse=false)

Returns the product stoichiometry matrix, $P$, with $P_{i j}$ the stoichiometric coefficient of the ith product within the jth reaction.

Note:

• Set sparse=true for a sparse matrix representation
• Note that constant species are not treated as products, but just components that modify the associated rate law.

netstoichmat(rn, sparse=false)

Returns the net stoichiometry matrix, $N$, with $N_{i j}$ the net stoichiometric coefficient of the ith species within the jth reaction.

Notes:

• Set sparse=true for a sparse matrix representation
• Caches the matrix internally within rn so subsequent calls are fast.
• Note that constant species are not treated as reactants, but just components that modify the associated rate law. As such they do not contribute to the net stoichiometry matrix.

reactionrates(network)

Given a ReactionSystem, returns a vector of the symbolic reaction rates for each reaction.

hasnoisescaling(reaction::Reaction)

Returns true if the input reaction has the noise_scaing metadata field assigned, else false.

Arguments:

• reaction: The reaction we wish to check for the noise_scaing metadata field.

Example:

reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
hasnoisescaling(reaction)

getnoisescaling(reaction::Reaction)

Returns noise_scaing metadata field for the input reaction.

Arguments:

• reaction: The reaction we wish to retrieve the noise_scaing metadata field.

Example:

reaction = @reaction k, 0 --> X, [noise_scaling=0.0]
getnoisescaling(reaction)

hasdescription(reaction::Reaction)

Returns true if the input reaction has the description metadata field assigned, else false.

Arguments:

• reaction: The reaction we wish to check for the description metadata field.

Example:

reaction = @reaction k, 0 --> X, [description="A reaction"]
hasdescription(reaction)

getdescription(reaction::Reaction)

Returns description metadata field for the input reaction.

Arguments:

• reaction: The reaction we wish to retrieve the description metadata field.

Example:

reaction = @reaction k, 0 --> X, [description="A reaction"]
getdescription(reaction)

hasmisc(reaction::Reaction)

Returns true if the input reaction has the misc metadata field assigned, else false.

Arguments:

• reaction: The reaction we wish to check for the misc metadata field.

Example:

reaction = @reaction k, 0 --> X, [misc="A reaction"]
hasmisc(reaction)

getmisc(reaction::Reaction)

Returns misc metadata field for the input reaction.

Arguments:

• reaction: The reaction we wish to retrieve the misc metadata field.

Example:

reaction = @reaction k, 0 --> X, [misc="A reaction"]
getmisc(reaction)

Notes:

• The misc field can contain any valid Julia structure. This mean that Catalyst cannot check it

for symbolic variables that are added here. This means that symbolic variables (e.g. parameters of species) that are stored here are not accessible to Catalyst. This can cause troubles when e.g. creating a ReactionSystem programmatically (in which case any symbolic variables stored in the misc metadata field should also be explicitly provided to the ReactionSystem constructor).

U0Map

Metadata key for storing initial condition / species value mappings on a ReactionSystem. Intended for use by file parsers that need to preserve mappings in a format different from initial_conditions.

See also: has_u0_map, get_u0_map, set_u0_map

ParameterMap

Metadata key for storing parameter value mappings on a ReactionSystem. Intended for use by file parsers that need to preserve parameter mappings in a format different from initial_conditions.

See also: has_parameter_map, get_parameter_map, set_parameter_map

has_u0_map(rs::ReactionSystem)

Returns true if the ReactionSystem has a U0Map metadata entry.

get_u0_map(rs::ReactionSystem)

Returns the U0Map metadata from the ReactionSystem, or nothing if no U0Map has been set.

set_u0_map(rs::ReactionSystem, u0map)

Returns a newReactionSystem with the U0Map metadata set to u0map. The original system is not modified.

has_parameter_map(rs::ReactionSystem)

Returns true if the ReactionSystem has a ParameterMap metadata entry.

get_parameter_map(rs::ReactionSystem)

Returns the ParameterMap metadata from the ReactionSystem, or nothing if no ParameterMap has been set.

set_parameter_map(rs::ReactionSystem, pmap)

Returns a newReactionSystem with the ParameterMap metadata set to pmap. The original system is not modified.

extend(
    sys::ModelingToolkitBase.AbstractSystem,
    basesys::ModelingToolkitBase.AbstractSystem;
    name,
    description,
    gui_metadata
) -> ReactionSystem{Catalyst.NetworkProperties{Int64, SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{SymReal}}}

Extend basesys with sys. This can be thought of as the merge operation on systems. Values in sys take priority over duplicates in basesys (for example, initial conditions).

By default, the resulting system inherits sys's name and description.

The & operator can also be used for this purpose. sys & basesys is equivalent to extend(sys, basesys).

See also compose.

extend(
    sys,
    basesys::Array{T<:ModelingToolkitBase.AbstractSystem, 1}
) -> Any

Extend sys with all systems in basesys in order.

ModelingToolkitBase.extend(sys::ReactionSystem, rs::ReactionSystem; name::Symbol=nameof(sys))

Extends the indicated ReactionSystem with another ReactionSystem.

Notes:

• Only ReactionSystems can be used to extend a ReactionSystem.
• Returns a new ReactionSystem and does not modify rs.
• By default, the new ReactionSystem will have the same name as sys.

compose(sys, systems; name)

Compose multiple systems together. This adds all of systems as subsystems of sys. The resulting system inherits the name of sys by default.

The ∘ operator can also be used for this purpose. sys ∘ basesys is equivalent to compose(sys, basesys).

See also extend.

compose(syss...; name) -> Any

Syntactic sugar for adding all systems in syss as the subsystems of first(syss).

ModelingToolkitBase.compose(sys::ReactionSystem, systems::AbstractArray; name = nameof(sys))

Compose the indicated ReactionSystem with one or more ReactionSystems.

Notes:

• Only ReactionSystems can be composed with a ReactionSystem.
• Returns a new ReactionSystem and does not modify sys.
• By default, the new ReactionSystem will have the same name as sys.
• Brownians and jumps from subsystems are collected at flatten time via recursive accessors.

flatten(sys::System) -> System
flatten(sys::System, noeqs) -> System

Flatten the hierarchical structure of a system, collecting all equations, unknowns, etc. into one top-level system after namespacing appropriately.

ModelingToolkitBase.flatten(rs::ReactionSystem)

Merges all subsystems of the given ReactionSystem up into rs.

Notes:

• Returns a new ReactionSystem that represents the flattened system.
• All Reactions within subsystems are namespaced and merged into the list of Reactions of rs. The merged list is then available as reactions(rs).
• All algebraic and differential equations are merged in the equations of rs.
• Only ReactionSystems are supported as subsystems when flattening.
• rs.networkproperties is reset upon flattening.
• The default value of combinatoric_ratelaws will be the logical or of all ReactionSystems.

plot_network(rn::ReactionSystem; kwargs...)

Converts a ReactionSystem into a GraphMakie plot of the species reaction graph (or Petri net representation). Reactions correspond to small green circles, and species to blue circles.

Notes:

• Black arrows from species to reactions indicate reactants, and are labelled with their input stoichiometry.
• Black arrows from reactions to species indicate products, and are labelled with their output stoichiometry.
• Red arrows from species to reactions indicate that species is used within the rate expression. For example, in the reaction k*A, B --> C, there would be a red arrow from A to the reaction node. In k*A, A+B --> C, there would be red and black arrows from A to the reaction node.

For a list of accepted keyword arguments to the graph plot, please see the GraphMakie documentation.

plot_complexes(rn::ReactionSystem; show_rate_labels = false, kwargs...)

Creates a GraphMakie plot of the Catalyst.ReactionComplexs in rn. Reactions correspond to arrows and reaction complexes to blue circles.

Notes:

• Black arrows from complexes to complexes indicate reactions whose rate is a parameter or a Number. i.e. k, A --> B.
• Red arrows from complexes to complexes indicate reactions whose rate constants

depends on species. i.e. k*C, A --> B for C a species.

• The show_rate_labels keyword, if set to true, will annotate each edge

with the rate constant for the reaction.

For a list of accepted keyword arguments to the graph plot, please see the GraphMakie documentation.

oderatelaw(rx; combinatoric_ratelaw=true)

Given a Reaction, return the symbolic reaction rate law used in generated ODEs for the reaction. Note, for a reaction defined by

k*X*Y, X+Z --> 2X + Y

the expression that is returned will be k*X(t)^2*Y(t)*Z(t). For a reaction of the form

k, 2X+3Y --> Z

the expression that is returned will be k * (X(t)^2/2) * (Y(t)^3/6).

Notes:

• Allocates
• combinatoric_ratelaw=true uses factorial scaling factors in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S^2/2!. If combinatoric_ratelaw=false then the ratelaw is k*S^2, i.e. the scaling factor is ignored.

jumpratelaw(rx; combinatoric_ratelaw=true)

Given a Reaction, return the symbolic reaction rate law used in generated stochastic chemical kinetics model SSAs for the reaction. Note, for a reaction defined by

k*X*Y, X+Z --> 2X + Y

the expression that is returned will be k*X^2*Y*Z. For a reaction of the form

k, 2X+3Y --> Z

the expression that is returned will be k * binomial(X,2) * binomial(Y,3).

Notes:

• Allocates
• combinatoric_ratelaw=true uses binomials in calculating the rate law, i.e. for 2S -> 0 at rate k the ratelaw would be k*S*(S-1)/2. If combinatoric_ratelaw=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor.

mm(X,v,K) = v*X / (X + K)

A Michaelis-Menten rate function.

mmr(X,v,K) = v*K / (X + K)

A repressive Michaelis-Menten rate function.

hill(X,v,K,n) = v*(X^n) / (X^n + K^n)

A Hill rate function.

hillr(X,v,K,n) = v*(K^n) / (X^n + K^n)

A repressive Hill rate function.

hillar(X,Y,v,K,n) = v*(X^n) / (X^n + Y^n + K^n)

An activation/repressing Hill rate function.

mtkcompile(
    sys;
    additional_passes,
    inputs,
    outputs,
    disturbance_inputs,
    split,
    kwargs...
)

Compile the given system into a form that ModelingToolkitBase can generate code for. Also performs order reduction for ODEs and handles simple discrete/implicit-discrete systems.

Keyword Arguments

• fully_determined=true controls whether or not an error will be thrown if the number of equations don't match the number of inputs, outputs, and equations.
• inputs, outputs and disturbance_inputs are passed as keyword arguments.All inputs get converted to parameters and are allowed to be unconnected, allowing models where n_unknowns = n_equations - n_inputs.

setdefaultnoisescaling(rs::ReactionSystem, noisescaling)

Creates an updated ReactionSystem. This is the old ReactionSystem, but each Reaction that does not have a noise_scaling metadata have its noisescaling metadata updated. The input ReactionSystem is not mutated. Any subsystems of rs have their `noisescaling` metadata updated as well.

Arguments:

• rs::ReactionSystem: The ReactionSystem which you wish to remake.
• noise_scaling: The updated noise scaling terms

@compound

Macro that creates a compound species, which is composed of smaller component species.

Example:

t = default_t()
@species C(t) O(t)
@compound CO2(t) ~ C + 2O

Notes:

• The component species must be defined before using the @compound macro.

@compounds

Macro that creates several compound species, which each is composed of smaller component species. Uses the same syntax as @compound, but with one compound species per line.

Example:

t = default_t()
@species C(t) H(t) O(t)
@compounds
    CH4(t) = C + 4H
    O2(t) = 2O
    CO2(t) = C + 2O
    H2O(t) = 2H + O
end

Notes:

• The component species must be defined before using the @compound macro.

iscompound(s)

Returns true if the input is a compound species (else false).

components(s)

Returns a vector with a list of all the components of a compound species (created using e.g. the @compound macro).

coefficients(s)

Returns a vector with a list of all the stoichiometric coefficients of the components of a compound species (created using e.g. the @compound macro).

component_coefficients(s)

Returns a Vector{Pair{Symbol,Int64}} listing, for a compound species (created using e.g. the @compound macro), all the components and their stoichiometric coefficients.

validate_units(rx::Reaction; info::String = "", warn::Bool = true)

Check that all substrates and products within the given Reaction have the same units, and that the units of the reaction's rate expression are internally consistent (i.e. if the rate involves sums, each term in the sum has the same units).

validate_units(rs::ReactionSystem; info::String="", warn::Bool = true)

Check that all species in the ReactionSystem have the same units, and that the rate laws of all reactions reduce to units of (species units) / (time units). Also validates unit consistency of non-reaction equations.

Uses catalyst_get_unit for SymbolicDimensions-preserving unit inference, avoiding the floating-point precision loss that occurs with MTKBase's get_unit when using non-SI units like M or μM.

Notes:

• Correctly handles only_use_rate=true reactions (does not multiply substrate units into the rate).
• Assumes reaction-local rate-expression checks (e.g. additive-term consistency) were already performed on each Reaction (for example via unit_checks = true at Reaction construction time, or by calling validate_units(rx) separately).
• If all species/time/parameters are unitless, reaction-rate dimensional checks are skipped. This mode assumes rate/equation expressions do not include literal dimensional quantities (for example us"..." constants), which are currently unsupported model inputs.
• Does not check subsystems, use flatten(rs) and then call validate_units on the flattened system if you want to check the full composed system.
• Does not require that non-species variables have consistent units (outside of the equations in which they appear).
• Does not handle events or user-provided jumps.

assert_valid_units(rx::Reaction; info::String = "")

Run strict unit validation on a Reaction. Throws UnitValidationError if any unit inconsistency is detected.

assert_valid_units(rs::ReactionSystem; info::String = "")

Run strict unit validation on a ReactionSystem. Throws UnitValidationError if any unit inconsistency is detected.

unit_validation_report(rx::Reaction; info::String = "")

Run unit validation on a Reaction and return a UnitValidationReport containing both overall validity and structured issue diagnostics.

unit_validation_report(rs::ReactionSystem; info::String = "")

Run unit validation on a ReactionSystem and return a UnitValidationReport containing both overall validity and structured issue diagnostics.

UnitValidationIssue

A single unit-validation diagnostic entry.

UnitValidationReport

Structured output from unit validation.

UnitValidationError

Exception thrown by strict unit-validation entrypoints (e.g. assert_valid_units). Wraps the full validation report.

struct LatticeReactionSystem{Q, R, S, T} <: ModelingToolkitBase.AbstractSystem

A representation of a spatial system of chemical reactions on a discrete (lattice) space.

Fields

• reactionsystem: The (non-spatial) reaction system within each vertex.

• spatial_reactions: The spatial reactions defined between individual vertices.

• lattice: The lattice on which the (discrete) spatial system is defined.

• num_verts: The number of vertices (compartments).

• num_edges: The number of edges.

• num_species: The number of species.

• spatial_species: List of species that may move spatially.

• parameters: All parameters related to the lattice reaction system (both those whose values are tied to vertices and edges).

• vertex_parameters: Parameters which values are tied to vertices, e.g. that possibly could have unique values at each vertex of the system.

• edge_parameters: Parameters whose values are tied to edges (adjacencies), e.g. that possibly could have unique values at each edge of the system.

• edge_iterator: An iterator over all the lattice's edges. Currently, the format is always a Vector{Pair{Int64,Int64}}. However, in the future, different types could potentially be used for different types of lattice (E.g. for a Cartesian grid, we do not technically need to enumerate each edge)

Arguments:

• rs: The non-spatial ReactionSystem model that is expanded to a spatial model.
• srs: A vector of spatial reactions. These provide the rules for how species may move spatially.
• lattice: Either a Cartesian grid, a masked grid, or a graph. This describes the discrete space

to which the non-spatial model is expanded.

Keyword Arguments:

• diagonal_connections = false: Only relevant for Cartesian and masked lattices. If true,

diagonally adjacent compartments are considered adjacent, and spatial reactions in between these are possible.

Example:

# Fetch packages.
using Catalyst, OrdinaryDiffEqDefault
import CairoMakie

# Creates the `LatticeReactionSystem` model.
rs = @reaction_network begin
    (p,d), 0 <--> X
end
diffusion_rx = @transport_reaction D X
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(rs, [diffusion_rx], lattice)

# Simulates the model (using ODE and jumps).
u0 = [:X => rand(5,5)]
tspan = (0.0, 1.0)
ps = [:p => 1.0, :d => 0.5, :D => 0.1]
oprob = ODEProblem(lrs, u0, tspan, ps)
osol = solve(oprob)

# Saves an animation of the solution to the file "lattice_animation.mp4".
lattice_animation(osol, :X, lrs, "lattice_animation.mp4")

Notes:

• Spatial modelling in Catalyst is still a work in progress, any feedback (or contributions) to this

is highly welcome.

• LatticeReactionSystems are primarily intended to model systems in discrete space. Modelling

continuous space systems with them is possible, but requires the user to determine the discretisation (the lattice). Better support for continuous space models is a work in progress.

• Catalyst contains extensive documentation on spatial modelling, which can be found here.

reactionsystem(lrs::LatticeReactionSystem)

Returns the non-spatial ReactionSystem stored in a LatticeReactionSystem.

spatial_reactions(lrs::LatticeReactionSystem)

Returns a vector with all the spatial reactions stored in a LatticeReactionSystem.

lattice(lrs::LatticeReactionSystem)

Returns the lattice stored in a LatticeReactionSystem.

num_verts(lrs::LatticeReactionSystem)

Returns the number of vertices (i.e. compartments) in the lattice stored in a LatticeReactionSystem.

num_edges(lrs::LatticeReactionSystem)

Returns the number of edges (i.e. connections between vertices) in the lattice stored in a LatticeReactionSystem.

num_species(lrs::LatticeReactionSystem)

Returns the number of species that a LatticeReactionSystem contains.

spatial_species(lrs::LatticeReactionSystem)

Returns the species that can move spatially in a LatticeReactionSystem.

vertex_parameters(lrs::LatticeReactionSystem)

Returns all the parameters of a LatticeReactionSystem whose values are tied to vertices.

edge_parameters(lrs::LatticeReactionSystem)

Returns all the parameters of a LatticeReactionSystem whose values are tied to edges.

edge_iterator(lrs::LatticeReactionSystem)

Returns an iterator over all of the edges in the lattice stored in a LatticeReactionSystem. Each edge is a Pair{Int64, Int64}, taking the source vertex to the destination vertex.

is_transport_system(lrs::LatticeReactionSystem)

Returns true if all spatial reactions in lrs are TransportReactions.

has_cartesian_lattice(lrs::LatticeReactionSystem)

Returns true if lrs was created using a cartesian grid lattice (e.g. created via CartesianGrid(5,5)). Otherwise, returns false.

has_masked_lattice(lrs::LatticeReactionSystem)

Returns true if lrs was created using a masked grid lattice (e.g. created via [true true; true false]). Otherwise, returns false.

has_grid_lattice(lrs::LatticeReactionSystem)

Returns true if lrs was created using a cartesian or masked grid lattice. Otherwise, returns false.

has_graph_lattice(lrs::LatticeReactionSystem)

Returns true if lrs was created using a graph grid lattice (e.g. created via path_graph(5)). Otherwise, returns false.

grid_size(lrs::LatticeReactionSystem)

Returns the size of lrs's lattice (only if it is a cartesian or masked grid lattice). E.g. for a lattice CartesianGrid(4,6), (4,6) is returned.

grid_dims(lrs::LatticeReactionSystem)

Returns the number of dimensions of lrs's lattice (only if it is a cartesian or masked grid lattice). The output is either 1, 2, or 3.

make_edge_p_values(lrs::LatticeReactionSystem, make_edge_p_value::Function)

Generates edge parameter values for a lattice reaction system. Only works for (Cartesian or masked) grid lattices (without diagonal adjacencies).

Input:

• lrs: The lattice reaction system for which values should be generated.
• make_edge_p_value: a function describing a rule for generating the edge parameter values.

Output: - ep_vals: A sparse matrix of size (numverts,numverts) (where numverts is the number of vertices in lrs). Here, eps[i,j] is filled only if there is an edge going from vertex i to vertex j. The value of eps[i,j] is determined by `makeedgepvalue`.

Here, make_edge_p_value should take two arguments, src_vert and dst_vert, which correspond to the grid indices of an edge's source and destination vertices, respectively. It outputs a single value, which is the value assigned to that edge.

Example: In the following example, we assign the value 0.1 to all edges, except for the one leading from vertex (1,1) to vertex (1,2), to which we assign the value 1.0.

using Catalyst
rn = @reaction_network begin
    (p,d), 0 <--> X
end
tr = @transport_reaction D X
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(rn, [tr], lattice)

function make_edge_p_value(src_vert, dst_vert)
    if src_vert == (1,1) && dst_vert == (1,2)
        return 1.0
    else
        return 0.1
    end
end

D_vals = make_edge_p_values(lrs, make_edge_p_value)

make_directed_edge_values(lrs::LatticeReactionSystem, x_vals::Tuple{T,T}, y_vals::Tuple{T,T} = (undef,undef),
                 z_vals::Tuple{T,T} = (undef,undef)) where {T}

Generates edge parameter values for a lattice reaction system. Only works for (Cartesian or masked) grid lattices (without diagonal adjacencies). Each dimension (x, and possibly y and z), and direction has assigned its own constant edge parameter value.

Input: - lrs: The lattice reaction system for which values should be generated. - x_vals::Tuple{T,T}: The values in the increasing (from a lower x index to a higher x index) and decreasing (from a higher x index to a lower x index) direction along the x dimension. - y_vals::Tuple{T,T}: The values in the increasing and decreasing direction along the y dimension. Should only be used for 2 and 3-dimensional grids. - z_vals::Tuple{T,T}: The values in the increasing and decreasing direction along the z dimension. Should only be used for 3-dimensional grids.

Output: - ep_vals: A sparse matrix of size (numverts,numverts) (where numverts is the number of vertices in lrs). Here, eps[i,j] is filled only if there is an edge going from vertex i to vertex j. The value of eps[i,j] is determined by the `xvals,yvals, andzvals` Tuples, and vertices i and j's relative position in the grid.

It should be noted that two adjacent vertices will always be different in exactly a single dimension (x, y, or z). The corresponding tuple determines which value is assigned.

Example: In the following example, we wish to have diffusion in the x dimension, but a constant flow from low y values to high y values (so not transportation from high to low y). We achieve it in the following manner:

using Catalyst
rn = @reaction_network begin
    (p,d), 0 <--> X
end
tr = @transport_reaction D X
lattice = CartesianGrid((5,5))
lrs = LatticeReactionSystem(rn, [tr], lattice)

D_vals = make_directed_edge_values(lrs, (0.1, 0.1), (0.1, 0.0))

Here, since we have a 2d grid, we only provide the first two Tuples to make_directed_edge_values.

lat_getu(sim_struct, sp, lrs::LatticeReactionSystem)

For a problem or integrators, retrieves its u values. For non-lattice models, this is can be done through direct interfacing (e.g. prob[X]). However, for LatticeReactionSystem-based problems and integrators, this function must be used instead. The output format depends on the lattice (a dense array for cartesian grid lattices, a sparse array for masked grid lattices, and a vector for graph lattices). This format is similar to which is used to designate species initial conditions.

Arguments:

• sim_struct: The simulation structure which u value we wish to retrieve. Can be either a ODEProblem, JumpProblem, or an integrator derived from either of these.
• sp: The species which value we wish to update. Can be provided either in its symbolic form (e.g. X) or as a symbol (e.g. :X).
• lrs: The LatticeReactionSystem which was used to generate the structure we wish to modify.

Notes:

• Even if the species is spatially uniform, a full array with its values across all vertices will be retrieved.

Example:

# Prepare `LatticeReactionSystem`s.
using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
lrs = LatticeReactionSystem(rs, [tr], CartesianGrid((2,3)))

# Prepares a corresponding ODEProblem.
u0 = [:X1 => [1.0 2.0 3.0; 4.0 5.0 6.0], :X2 => 2.0]
tspan = (0.0, 50.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.01]
oprob = ODEProblem(lrs, u0, tspan, ps)

# Updates the `ODEProblem`.
lat_getu(oprob, :X1, lrs) # Retrieves the value of `X1`.

lat_getu(sol, sp, lrs::LatticeReactionSystem; t = nothing)

A function for retrieving the solution of a LatticeReactionSystem-based simulation on various desired forms. Generally, for LatticeReactionSystems, the values in sol is ordered in a way which is not directly interpretable by the user. Furthermore, the normal Catalyst interface for solutions (e.g. sol[:X]) does not work for these solutions. Hence this function is used instead.

The output is a vector, which in each position contains sp's value (either at a time step of time, depending on the input t). Its shape depends on the lattice (using a similar form as heterogeneous initial conditions). I.e. for a NxM cartesian grid, the values are NxM matrices. For a masked grid, the values are sparse matrices. For a graph lattice, the values are vectors (where the value in the n'th position corresponds to sp's value in the n'th vertex).

Arguments:

• sol: The solution from which we wish to retrieve some values.
• sp: The species which value we wish to update. Can be provided either in its symbolic form (e.g. X) or as a symbol (e.g. :X).
• lrs: The LatticeReactionSystem which was simulated to generate the solution.
• t = nothing: If nothing, we simply return the solution across all saved time steps (default). If t instead is a vector (or range of values), returns the solution interpolated at these time points.

Notes:

• The lat_getu is not optimised for performance. However, it should still be quite performant, but there might be some limitations if called a very large number of times.
• Long-term it is likely that this function gets replaced with a sleeker interface.

Example:

using Catalyst, OrdinaryDiffEqDefault

# Prepare `LatticeReactionSystem`s.
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
lrs = LatticeReactionSystem(rs, [tr], CartesianGrid((2,2)))

# Create problems.
u0 = [:X1 => 1, :X2 => 2]
tspan = (0.0, 10.0)
ps = [:k1 => 1, :k2 => 2.0, :D => 0.1]

oprob = ODEProblem(lrs1, u0, tspan, ps)
osol = solve(oprob)
lat_getu(osol, :X1, lrs) # Returns the value of X1 at each time step.
lat_getu(osol, :X1, lrs; t = 0.0:10.0) # Returns the value of X1 at times 0.0, 1.0, ..., 10.0

lat_setu!(sim_struct, sp, lrs::LatticeReactionSystem, u)

For a problem or integrators, update its u vector with the input u. For non-lattice models, this is can be done through direct interfacing (e.g. prob[X] = 1.0). However, for LatticeReactionSystem-based problems and integrators, this function must be used instead.

Arguments:

• sim_struct: The simulation structure which u value we wish to update. Can be either a ODEProblem, JumpProblem, or an integrator derived from either of these.
• sp: The species which value we wish to update. Can be provided either in its symbolic form (e.g. X) or as a symbol (e.g. :X).
• lrs: The LatticeReactionSystem which was used to generate the structure we wish to modify.
• u: The species's new values. Must be given in a form which is also a valid initial input to the ODEProblem/JumpProblem.

Example:

# Prepare `LatticeReactionSystem`s.
using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
lrs = LatticeReactionSystem(rs, [tr], CartesianGrid((2,3)))

# Prepares a corresponding ODEProblem.
u0 = [:X1 => [1.0 2.0 3.0; 4.0 5.0 6.0], :X2 => 2.0]
tspan = (0.0, 50.0)
ps = [:k1 => 2.0, :k2 => 1.0, :D => 0.01]
oprob = ODEProblem(lrs, u0, tspan, ps)

# Updates the `ODEProblem`.
lat_setu!(oprob, :X1, lrs, 0.0) # Sets `X1` to uniformly 0 across the lattice.
lat_setu!(oprob, :X2, lrs, [1.0 0.0 0.0; 0.0 0.0 0.0]) # Sets `X2` to `1.0` in one vertex, and 0 elsewhere.

lat_getp(sim_struct, p, lrs::LatticeReactionSystem)

For a problem or integrators, retrieves its p values. For non-lattice models, this is can be done through direct interfacing (e.g. prob[p]). However, for LatticeReactionSystem-based problems and integrators, this function must be used instead. The output format depends on the lattice (a dense array for cartesian grid lattices, a sparse array for masked grid lattices, and a vector for graph lattices). This format is similar to what is used to designate parameter initial values.

Arguments:

• sim_struct: The simulation structure which p value we wish to retrieve. Can be either a ODEProblem,

JumpProblem, or an integrator derived from either of these.

• p: The species which value we wish to update. Can be provided either in its symbolic form (e.g. k) or as a symbol (e.g. :k).
• lrs: The LatticeReactionSystem which was used to generate the structure we wish to modify.

Notes:

• Even if the parameter is spatially uniform, a full array with its values across all vertices will be retrieved.

Example:

# Prepare `LatticeReactionSystem`s.
using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
lrs = LatticeReactionSystem(rs, [tr], CartesianGrid((2,3)))

# Prepares a corresponding ODEProblem.
u0 = [:X1 => 1.0, :X2 => 2.0]
tspan = (0.0, 50.0)
ps = [:k1 => [1.0 2.0 3.0; 4.0 5.0 6.0], :k2 => 1.0, :D => 0.01]
oprob = ODEProblem(lrs, u0, tspan, ps)

# Updates the `ODEProblem`.
lat_getp(oprob, :k1, lrs) # Retrieves the value of `k1`.

lat_setp!(sim_struct, p, lrs::LatticeReactionSystem, p_val)

For a problem or integrators, update its p vector with the input p_val. For non-lattice models, this is can be done through direct interfacing (e.g. prob[p] = 1.0). However, for LatticeReactionSystem-based problems and integrators, this function must be used instead.

Arguments:

• sim_struct: The simulation structure which u value we wish to update. Can be either a ODEProblem, JumpProblem, or an integrator derived from either of these.
• p: The species which value we wish to update. Can be provided either in its symbolic form (e.g. k) or as a symbol (e.g. :k).
• lrs: The LatticeReactionSystem which was used to generate the structure we wish to modify.
• p_val: The parameter's new values. Must be given in a form which is also a valid initial input to the ODEProblem/JumpProblem.

Example:

# Prepare `LatticeReactionSystem`s.
using Catalyst
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
lrs = LatticeReactionSystem(rs, [tr], CartesianGrid((2,3)))

# Prepares a corresponding ODEProblem.
u0 = [:X1 => 1.0, :X2 => 2.0]
tspan = (0.0, 50.0)
ps = [:k1 => [1.0 2.0 3.0; 4.0 5.0 6.0], :k2 => 1.0, :D => 0.01]
oprob = ODEProblem(lrs, u0, tspan, ps)

# Updates the `ODEProblem`.
lat_setp!(oprob, :k1, lrs, 0.0) # Sets `k1` to uniformly 0 across the lattice.
lat_setp!(oprob, :k2, lrs, [1.0 0.0 0.0; 0.0 0.0 0.0]) # Sets `k2` to `1.0` in one vertex, and 0 elsewhere.

rebuild_lat_internals!(sciml_struct)

Rebuilds the internal functions for simulating a LatticeReactionSystem. Whenever a problem or integrator has had its parameter values updated, this function should be called for the update to be taken into account. For ODE simulations, rebuild_lat_internals! needs only to be called when

• An edge parameter has been updated.
• When a parameter with spatially homogeneous values has been given spatially heterogeneous values (or vice versa).

Arguments:

• sciml_struct: The problem (e.g. an ODEProblem) or an integrator which we wish to rebuild.

Notes:

• Currently does not work for DiscreteProblems, JumpProblems, or their integrators.
• The function is not built with performance in mind, so avoid calling it multiple times in performance-critical applications.

Example:

# Creates an initial `ODEProblem`
rs = @reaction_network begin
    (k1,k2), X1 <--> X2
end
tr = @transport_reaction D X1
grid = CartesianGrid((2,2))
lrs = LatticeReactionSystem(rs, [tr], grid)

u0 = [:X1 => 2, :X2 => [5 6; 7 8]]
tspan = (0.0, 10.0)
ps = [:k1 => 1.5, :k2 => [1.0 1.5; 2.0 3.5], :D => 0.1]

oprob = ODEProblem(lrs, u0, tspan, ps)

# Updates parameter values.
oprob.ps[:ks] = [2.0 2.5; 3.0 4.5]
oprob.ps[:D] = 0.05

# Rebuilds `ODEProblem` to make changes have an effect.
rebuild_lat_internals!(oprob)

lattice_plot(sol, sp, lrs::LatticeReactionSystem, filename::String; t = sol.tspan[2], kwargs...)

Creates a plot of a LatticeReactionSystem simulation. The plot is created at the time point specified by t (defaults to the simulation's final time point).

Arguments (all lattices):

• sol: The simulation we wish to plot.
• sp: The species whose values we wish to plot. Can be provided either in its symbolic form or as a symbol.
• lrs: The LatticeReactionSystem which was simulated.
• t = sol.t[end]: The time point at which we wish to plot the solution

In addition, depending on the type of lattice used, the following optional arguments might be relevant.

Arguments (1d lattices):

• markersize = 20: The size of the markers marking each compartment's value.

Arguments (Graph & 2d lattices):

• colormap = :BuGn_7: The colour map with which we display the species amounts in the animation.
• plot_min = nothing: The minimum value for the colour scale (values less than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value (across the entire simulation, not just at the plotted time value).
• plot_max = nothing: The maximum value for the colour scale (values more than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value (across the entire simulation, not just at the plotted time value).

Arguments (Graph lattices):

• node_size = 50: The size of the compartments in the plot.
• layout = Spring(): The layout for the graph nodes in the plot. Can be provided as a vector, where the i'th element is a 2-valued tuple (determining the i'th compartment's y and x positions, respectively).

Notes:

• For masked lattices, there are no value displayed for grid points which do not correspond to a compartments.
• The current plotting interface is a work in progress, and modifications are expected. if you have any feedback, please contact the package authors.
• Additional arguments can be passed to lattice_plot, which then will be passed to Makie's lines plotting command.

lattice_animation(sol, sp, lrs::LatticeReactionSystem, filename::String; kwargs...)

Creates an animation of a LatticeReactionSystem simulation. The animation is saved to a file, whose name is provided in the filename argument.

Arguments (all lattices):

• sol: The simulation we wish to animate.
• sp: The species which values we wish to animate. Can be provided either in its symbolic form or as a symbol.
• lrs: The LatticeReactionSystem which was simulated.
• filename: The name of the file to which we wish to save the animation.
• nframes = 200: The number of frames in the animation (these are evenly samples across the simulation).
• framerate = 20: The frame rate of the animation.
• ttitle = true: Whether to add a title showing the simulation's time throughout the animation.

In addition, depending on the type of lattice used, the following optional arguments might be relevant.

Arguments (1d lattices):

• markersize = 20: The size of the markers marking each compartment's value.
• plot_min = nothing: The y-scale's minimum. If nothing, use the simulation's minimum value.
• plot_max = nothing: The y-scale's maximum. If nothing, use the simulation's maximum value.

Arguments (Graph & 2d lattices):

• colormap = :BuGn_7: The colour map with which we display the species amounts in the animation.
• plot_min = nothing: The minimum value for the colour scale (values less than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value.
• plot_max = nothing: The maximum value for the colour scale (values more than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value.

Arguments (Graph lattices):

• node_size = 50: The size of the compartments in the plot.
• layout = Spring(): The layout for the graph nodes in the plot. Can be provided as a vector, where the i'th element is a 2-valued tuple (determining the i'th compartment's y and x positions, respectively).

Notes:

• For masked lattices, there are no value displayed for grid points which do not correspond to a compartments.
• The current animation interface if a work in progress, and modifications are expected. if you have any feedback, please contact the package authors.
• Additional arguments can be passed to lattice_animation, which then will be passed to Makie's heatmap plotting command.

lattice_kymograph(sol, sp, lrs::LatticeReactionSystem, kwargs...)

Creates a kymograph of a LatticeReactionSystem simulation based on a Cartesian or masked lattice. The plot shows the compartments on the y-axis, and the time development of the system's state along the x-axis. Species amounts are shown as a heatmap.

Arguments (all lattices):

• sol: The simulation we wish to plot.
• sp: The species whose values we wish to plot. Can be provided either in its symbolic form or as a symbol.
• lrs: The LatticeReactionSystem which was simulated.
• colormap = :BuGn_7: The colour map with which we display the species amounts in the kymograph.
• nframes = 200: The number of time samples which the time series is sampled with.
• plot_min = nothing: The minimum value for the colour scale (values less than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value.
• plot_max = nothing: The maximum value for the colour scale (values more than this will be set at this value when the colour scale is computed). If nothing, use the simulation's minimum value.

Notes:

• For masked lattices, there are no value displayed for grid points which do not correspond to a compartments.
• The current plotting interface is a work in progress, and modifications are expected. if you have any feedback, please contact the package authors.
• Additional arguments can be passed to lattice_plot, which then will be passed to Makie's heatmap plotting command.