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.

struct ReactionSystem{V<:Catalyst.NetworkProperties} <: AbstractTimeDependentSystem

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

• defaults: The default values to use when initial conditions and/or parameters are not supplied in ODEProblem.

• connection_type: Type of the system

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

• 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{AbstractSystems}, vector of sub-systems. Can be ReactionSystems, ODESystems, or NonlinearSystems.
• name::Symbol, the name of the system (must be provided, or @named must be used).
• defaults::Dict, a dictionary mapping parameters to their default values and species to their default initial values.
• checks = true, boolean for whether to check units.
• networkproperties = NetworkProperties(), cache for network properties calculated via API functions.
• combinatoric_ratelaws = true, sets the default value of combinatoric_ratelaws used in calls to convert or calling various problem types 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).

Notes:

• ReactionSystems currently do rudimentary unit checking, requiring that all species have the same units, and all reactions have rate laws with units of (species units) / (time units). Unit checking can be disabled by passing the keyword argument checks=false.

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 ModelingToolkit.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 ModelingToolkit.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 ModelingToolkit
using ModelingToolkit: 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 = ODESystem([eq1], t)
@named osys2 = ODESystem([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: ```julia using ModelingToolkit using ModelingToolkit: tnounits as t, Dnounits as D @parameters p d @variables X(t) eq1 = D(X) ~ p - dX eq2 = 0 ~ p - dX @named osys = ODESystem([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 ModelingToolkit
using ModelingToolkit: 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 = ODESystem([eq1], t)
@named osys2 = ODESystem([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: ```julia using ModelingToolkit using ModelingToolkit: tnounits as t, Dnounits as D @parameters p d @variables X(t) eq1 = D(X) ~ p - dX eq2 = 0 ~ p - dX @named osys = ODESystem([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).

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

setdefaults!(rn, newdefs)

Sets the default (initial) values of parameters and species in the ReactionSystem, rn.

For example,

sir = @reaction_network SIR begin
    β, S + I --> 2I
    ν, I --> R
end
setdefaults!(sir, [:S => 999.0, :I => 1.0, :R => 1.0, :β => 1e-4, :ν => .01])

# or
t = default_t()
@parameter β ν
@species S(t) I(t) R(t)
setdefaults!(sir, [S => 999.0, I => 1.0, R => 0.0, β => 1e-4, ν => .01])

gives initial/default values to each of S, I and β

Notes:

• Can not be used to set default values for species, variables or parameters of subsystems or constraint systems. Either set defaults for those systems directly, or flatten to collate them into one system before setting defaults.
• Defaults can be specified in any iterable container of symbols to value pairs or symbolics to value pairs.

extend(
    sys::ModelingToolkit.AbstractSystem,
    basesys::ModelingToolkit.AbstractSystem;
    name,
    description,
    gui_metadata
) -> ReactionSystem{Catalyst.NetworkProperties{Int64, V}} where V<:SymbolicUtils.BasicSymbolic{Real}

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

See also compose.

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

Extends the indicated ReactionSystem with another AbstractSystem.

Notes:

• The AbstractSystem being added in must be an ODESystem, NonlinearSystem, or ReactionSystem currently.
• 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. The resulting system would inherit the first system's name.

See also extend.

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

Compose the indicated ReactionSystem with one or more AbstractSystems.

Notes:

• The AbstractSystem being added in must be an ODESystem, NonlinearSystem, or ReactionSystem currently.
• Returns a new ReactionSystem and does not modify rs.
• By default, the new ReactionSystem will have the same name as sys.

ModelingToolkit.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.
• Currently only ReactionSystems, NonlinearSystems and ODESystems are supported as sub-systems when flattening.
• rs.networkproperties is reset upon flattening.
• The default value of combinatoric_ratelaws will be the logical or of all ReactionSystems.

==(rx1::Reaction, rx2::Reaction)

Tests whether two Reactions are identical.

Notes:

• Ignores the order in which stoichiometry components are listed.
• Does not currently simplify rates, so a rate of A^2+2*A+1 would be considered different than (A+1)^2.

isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true,
    debug = false)

Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystems. Ignores the names of the systems in testing equality.

Notes:

• Does not currently simplify rates, so a rate of A^2+2*A+1 would be considered different than (A+1)^2.
• ignorenames = false is used when checking equality of sub and parent systems.
• Does not check that parent systems are the same.
• Pass debug = true to print out the field that caused the two systems to be considered different.

==(rn1::ReactionSystem, rn2::ReactionSystem)

Tests whether the underlying species, parameters and reactions are the same in the two ReactionSystems. Requires the systems to have the same names too.

Notes:

• Does not currently simplify rates, so a rate of A^2+2*A+1 would be considered different than (A+1)^2.
• Does not include defaults in determining equality.

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.

Base.convert(::Type{<:ODESystem},rs::ReactionSystem)

Convert a ReactionSystem to an ModelingToolkit.ODESystem.

Keyword args and default values:

• combinatoric_ratelaws=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!. Set combinatoric_ratelaws=false for a ratelaw of k*S^2, i.e. the scaling factor is ignored. Defaults to the value given when the ReactionSystem was constructed (which itself defaults to true).
• remove_conserved=false, if set to true will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
• expand_catalyst_funs = true, replaces Catalyst defined functions like hill(A,B,C,D) with their rational function representation when converting to another system type. Set to false` to disable.

Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)

Convert a ReactionSystem to an ModelingToolkit.NonlinearSystem.

Keyword args and default values:

• combinatoric_ratelaws = 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!. Set combinatoric_ratelaws=false for a ratelaw of k*S^2, i.e. the scaling factor is ignored. Defaults to the value given when the ReactionSystem was constructed (which itself defaults to true).
• remove_conserved = false, if set to true will calculate conservation laws of the underlying set of reactions (ignoring coupled ODE or algebraic equations). For each conservation law one steady-state equation is eliminated, and replaced with the conservation law. This ensures a non-singular Jacobian.
• conseqs_remake_warn = true, set to false to disable warning about remake and conservation laws. See the FAQ entry for more details.
• expand_catalyst_funs = true, replaces Catalyst defined functions like hill(A,B,C,D) with their rational function representation when converting to another system type. Set to false` to disable.

Base.convert(::Type{<:SDESystem},rs::ReactionSystem)

Convert a ReactionSystem to an ModelingToolkit.SDESystem.

Notes:

• combinatoric_ratelaws=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!. Set combinatoric_ratelaws=false for a ratelaw of k*S^2, i.e. the scaling factor is ignored. Defaults to the value given when the ReactionSystem was constructed (which itself defaults to true).
• remove_conserved=false, if set to true will calculate conservation laws of the underlying set of reactions (ignoring constraint equations), and then apply them to reduce the number of equations.
• expand_catalyst_funs = true, replaces Catalyst defined functions like hill(A,B,C,D) with their rational function representation when converting to another system type. Set to false` to disable.

Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)

Convert a ReactionSystem to an ModelingToolkit.JumpSystem.

Notes:

• combinatoric_ratelaws=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_ratelaws=false then the ratelaw is k*S*(S-1), i.e. the rate law is not normalized by the scaling factor. Defaults to the value given when the ReactionSystem was constructed (which itself defaults to true).
• Does not currently support ReactionSystems that include coupled algebraic or differential equations.
• Does not currently support continuous events as these are not supported by ModelingToolkit.JumpSystems.
• expand_catalyst_funs = true, replaces Catalyst defined functions like hill(A,B,C,D) with their rational function representation when converting to another system type. Set to false` to disable.
• save_positions = (true, true), indicates whether for any reaction classified as a VariableRateJump to save the solution before and/or after the jump occurs. Defaults to true for both.

struct JumpInputs{S<:JumpSystem, T<:SciMLBase.AbstractODEProblem}

Inputs for a JumpProblem from a given ReactionSystem.

Fields

• sys: The JumpSystem to define the problem over

• prob: The problem the JumpProblem should be defined over, for example DiscreteProblem

structural_simplify(sys; ...)
structural_simplify(
    sys,
    io;
    additional_passes,
    simplify,
    split,
    allow_symbolic,
    allow_parameter,
    conservative,
    fully_determined,
    kwargs...
)

Structurally simplify algebraic equations in a system and compute the topological sort of the observed equations in sys.

Optional Arguments:

• optional argument io may take a tuple (inputs, outputs). This will convert all inputs to parameters and allow them to be unconnected, i.e., simplification will allow models where n_unknowns = n_equations - n_inputs.

Optional Keyword Arguments:

• When simplify=true, the simplify function will be applied during the tearing process.
• allow_symbolic=false, allow_parameter=true, and conservative=false limit the coefficient types during tearing. In particular, conservative=true limits tearing to only solve for trivial linear systems where the coefficient has the absolute value of $1$.
• 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.
• sort_eqs=true controls whether equations are sorted lexicographically before simplification or not.

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 one each 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{Pari{Symbol,Int64}}, listing a compounds species (created using e.g. the @compound macro) all the coefficients and their stoichiometric coefficients.

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

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(rs::ReactionSystem, info::String="")

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

Notes:

• Does not check subsystems, constraint equations, or non-species variables.

symmap_to_varmap(sys, symmap)

Given a system and map of Symbols to values, generates a map from corresponding symbolic variables/parameters to the values that can be used to pass initial conditions and parameter mappings.

For example,

sir = @reaction_network sir begin
    β, S + I --> 2I
    ν, I --> R
end
subsys = @reaction_network subsys begin
    k, A --> B
end
@named sys = compose(sir, [subsys])

gives

Model sys with 3 equations
Unknowns (5):
  S(t)
  I(t)
  R(t)
  subsys₊A(t)
  subsys₊B(t)
Parameters (3):
  β
  ν
  subsys₊k

to specify initial condition and parameter mappings from symbols we can use

symmap = [:S => 1.0, :I => 1.0, :R => 1.0, :subsys₊A => 1.0, :subsys₊B => 1.0]
u0map  = symmap_to_varmap(sys, symmap)
pmap   = symmap_to_varmap(sys, [:β => 1.0, :ν => 1.0, :subsys₊k => 1.0])

u0map and pmap can then be used as input to various problem types.

Notes:

• Any Symbol, sym, within symmap must be a valid field of sys. i.e. sys.sym must be defined.

struct LatticeReactionSystem{Q, R, S, T} <: AbstractTimeDependentSystem

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 number of species that can move spatially that a LatticeReactionSystem contains.

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