Skip to content

Latest commit

 

History

History
394 lines (351 loc) · 14.1 KB

api.md

File metadata and controls

394 lines (351 loc) · 14.1 KB

[Catalyst.jl API](@id api)

CurrentModule = Catalyst

Reaction network generation and representation

Catalyst provides the @reaction_network macro for generating a complete network, stored as a ReactionSystem, which in turn is composed of Reactions. ReactionSystems can be converted to other ModelingToolkit.AbstractSystems, including a ModelingToolkit.ODESystem, ModelingToolkit.SDESystem, or ModelingToolkit.JumpSystem.

When using the @reaction_network macro, Catalyst will automatically attempt to detect what is a species and what is a parameter. Everything that appear as a substrate or product in some reaction will be treated as a species, while all remaining symbols will be considered parameters (corresponding to those symbols that only appear within rate expressions and/or as stoichiometric coefficients). I.e. in

rn = @reaction_network begin
    k*X, Y --> W
end

Y and W will all be classified as chemical species, while k and X will be classified as parameters.

The ReactionSystem generated by the @reaction_network macro is a ModelingToolkit.AbstractSystem that symbolically represents a system of chemical reactions. In some cases it can be convenient to bypass the macro and directly generate a collection of symbolic Reactions and a corresponding ReactionSystem encapsulating them. Below we illustrate with a simple SIR example how a system can be directly constructed, and demonstrate how to then generate from the ReactionSystem and solve corresponding chemical reaction ODE models, chemical Langevin equation SDE models, and stochastic chemical kinetics jump process models.

using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
t = default_t()
@parameters β γ
@species S(t) I(t) R(t)

rxs = [Reaction(β, [S,I], [I], [1,1], [2])
       Reaction(γ, [I], [R])]
@named rs = ReactionSystem(rxs, t)
rs = complete(rs)

u₀map    = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => 1/10000, γ => 0.01]
tspan    = (0.0, 250.0)

# solve as ODEs
odesys = convert(ODESystem, rs)
odesys = complete(odesys)
oprob = ODEProblem(odesys, u₀map, tspan, parammap)
sol = solve(oprob, Tsit5())
p1 = plot(sol, title = "ODE")

# solve as SDEs
sdesys = convert(SDESystem, rs)
sdesys = complete(sdesys)
sprob = SDEProblem(sdesys, u₀map, tspan, parammap)
sol = solve(sprob, EM(), dt=.01, saveat = 2.0)
p2 = plot(sol, title = "SDE")

# solve as jump process
jumpsys = convert(JumpSystem, rs)
jumpsys = complete(jumpsys)
u₀map    = [S => 999, I => 1, R => 0]
dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap)
jprob = JumpProblem(jumpsys, dprob, Direct())
sol = solve(jprob)
sol = solve(jprob; seed = 1234) # hide
p3 = plot(sol, title = "jump")
plot(p1, p2, p3; layout = (3,1))
Catalyst.PNG(plot(p1, p2, p3; layout = (3,1), fmt = :png, dpi = 200)) # hide

@reaction_network
@network_component
make_empty_network
@reaction
Reaction
ReactionSystem

[Options for the @reaction_network DSL](@id api_dsl_options)

We have [previously described](@ref dsl_advanced_options) how options allow one to supply additional information to a ReactionSystem created with the DSL. Here follows a list of all options currently available.

  • [parameters](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system parameters.
  • [species](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system species.
  • [variables](@ref dsl_advanced_options_declaring_species_and_parameters): Allows the designation of a set of symbols as system non-species variables.
  • [ivs](@ref dsl_advanced_options_ivs): Allows the designation of a set of symbols as system independent variables.
  • [compounds](@ref chemistry_functionality_compounds): Allows the designation of compound species.
  • [observables](@ref dsl_advanced_options_observables): Allows the designation of compound observables.
  • [default_noise_scaling](@ref simulation_intro_SDEs_noise_saling): Enables the setting of a default noise scaling expression.
  • [differentials](@ref constraint_equations_coupling_constraints): Allows the designation of differentials.
  • [equations](@ref constraint_equations_coupling_constraints): Allows the creation of algebraic and/or differential equations.
  • [continuous_events](@ref constraint_equations_events): Allows the creation of continuous events.
  • [discrete_events](@ref constraint_equations_events): Allows the creation of discrete events.
  • [combinatoric_ratelaws](@ref faq_combinatoric_ratelaws): Takes a single option (true or false), which sets whether to use combinatorial rate laws.
  • [require_declaration](@ref dsl_advanced_options_require_dec): Turns off all inference of parameters, species, variables, the default differential, and observables (requiring these to be explicitly declared using e.g. @species).

[ModelingToolkit and Catalyst accessor functions](@id api_accessor_functions)

A ReactionSystem is an instance of a ModelingToolkit.AbstractTimeDependentSystem, and has a number of fields that can be accessed using the Catalyst API and the ModelingToolkit.jl Abstract System Interface. Below we overview these components.

There are three basic sets of convenience accessors that will return information either from a top-level system, the top-level system and all sub-systems that are also ReactionSystems (i.e. the full reaction-network), or the top-level system, all subs-systems, and all constraint systems (i.e. the full model). To retrieve info from just a base ReactionSystem rn, ignoring sub-systems of rn, one can use the ModelingToolkit accessors (these provide direct access to the corresponding internal fields of the ReactionSystem)

  • ModelingToolkit.get_unknowns(rn) is a vector that collects all the species defined within rn, ordered by species and then non-species variables.
  • Catalyst.get_species(rn) is a vector of all the species variables in the system. The entries in get_species(rn) correspond to the first length(get_species(rn)) components in get_unknowns(rn).
  • ModelingToolkit.get_ps(rn) is a vector that collects all the parameters defined within reactions in rn.
  • ModelingToolkit.get_eqs(rn) is a vector that collects all the Reactions and Symbolics.Equation defined within rn, ordering all Reactions before Equations.
  • Catalyst.get_rxs(rn) is a vector of all the Reactions in rn, and corresponds to the first length(get_rxs(rn)) entries in get_eqs(rn).
  • ModelingToolkit.get_iv(rn) is the independent variable used in the system (usually t to represent time).
  • ModelingToolkit.get_systems(rn) is a vector of all sub-systems of rn.
  • ModelingToolkit.get_defaults(rn) is a dictionary of all the default values for parameters and species in rn.

The preceding accessors do not allocate, directly accessing internal fields of the ReactionSystem.

To retrieve information from the full reaction network represented by a system rn, which corresponds to information within both rn and all sub-systems, one can call:

  • ModelingToolkit.unknowns(rn) returns all species and variables across the system, all sub-systems, and all constraint systems. Species are ordered before non-species variables in unknowns(rn), with the first numspecies(rn) entries in unknowns(rn) being the same as species(rn).
  • species(rn) is a vector collecting all the chemical species within the system and any sub-systems that are also ReactionSystems.
  • ModelingToolkit.parameters(rn) returns all parameters across the system, all sub-systems, and all constraint systems.
  • ModelingToolkit.equations(rn) returns all Reactions and all Symbolics.Equations defined across the system, all sub-systems, and all constraint systems. Reactions are ordered ahead of Equations with the first numreactions(rn) entries in equations(rn) being the same as reactions(rn).
  • reactions(rn) is a vector of all the Reactions within the system and any sub-systems that are also ReactionSystems.

These accessors will generally allocate new arrays to store their output unless there are no subsystems. In the latter case the usually return the same vector as the corresponding get_* function.

Below we list the remainder of the Catalyst API accessor functions mentioned above.

Basic system properties

See [Programmatic Construction of Symbolic Reaction Systems](@ref programmatic_CRN_construction) for examples and [ModelingToolkit and Catalyst Accessor Functions](@ref api_accessor_functions) for more details on the basic accessor functions.

species
Catalyst.get_species
nonspecies
reactions
Catalyst.get_rxs
nonreactions
numspecies
numparams
numreactions
speciesmap
paramsmap
isautonomous

Coupled reaction/equation system properties

The following system property accessor functions are primarily relevant to reaction system [coupled to differential and/or algebraic equations](@ref constraint_equations).

ModelingToolkit.has_alg_equations
ModelingToolkit.alg_equations
ModelingToolkit.has_diff_equations
ModelingToolkit.diff_equations

Basic species properties

The following functions permits the querying of species properties.

isspecies
Catalyst.isconstant
Catalyst.isbc
Catalyst.isvalidreactant

Basic reaction properties

ismassaction
dependents
dependants
substoichmat
prodstoichmat
netstoichmat
reactionrates

[Reaction metadata](@id api_rx_metadata)

The following functions permits the retrieval of [reaction metadata](@ref dsl_advanced_options_reaction_metadata).

Catalyst.hasnoisescaling
Catalyst.getnoisescaling
Catalyst.hasdescription
Catalyst.getdescription
Catalyst.hasmisc
Catalyst.getmisc

[Functions to extend or modify a network](@id api_network_extension_and_modification)

ReactionSystems can be programmatically extended using ModelingToolkit.extend and ModelingToolkit.compose.

setdefaults!
ModelingToolkit.extend
ModelingToolkit.compose
Catalyst.flatten

Network analysis and representations

Note, currently API functions for network analysis and conservation law analysis do not work with constant species (currently only generated by SBMLToolkit).

conservationlaws
conservedquantities
conservedequations
conservationlaw_constants
ReactionComplexElement
ReactionComplex
reactioncomplexmap
reactioncomplexes
incidencemat
complexstoichmat
complexoutgoingmat
incidencematgraph
linkageclasses
deficiency
subnetworks
linkagedeficiencies
isreversible
isweaklyreversible
reset_networkproperties!

Network comparison

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

[Network visualization](@id network_visualization)

Latexify can be used to convert networks to LaTeX equations by

using Latexify
latexify(rn)

An optional argument, form allows using latexify to display a reaction network's ODE (as generated by the reaction rate equation) or SDE (as generated by the chemical Langevin equation) form:

latexify(rn; form=:ode)
latexify(rn; form=:sde)

(As of writing this, an upstream bug causes the SDE form to be erroneously displayed as the ODE form)

Finally, another optional argument (expand_functions=true) automatically expands functions defined by Catalyst (such as mm). To disable this, set expand_functions=false.

Reaction networks can be plotted using the GraphMakie extension, which is loaded whenever all of Catalyst, GraphMakie, and NetworkLayout are loaded (note that a Makie backend, like CairoMakie, must be loaded as well). The two functions for plotting networks are plot_network and plot_complexes, which are two distinct representations.

plot_network(::ReactionSystem)
plot_complexes(::ReactionSystem)

[Rate laws](@id api_rate_laws)

As the underlying ReactionSystem is comprised of ModelingToolkit expressions, one can directly access the generated rate laws, and using ModelingToolkit tooling generate functions or Julia Exprs from them.

oderatelaw
jumpratelaw
mm
mmr
hill
hillr
hillar

Transformations

Base.convert
JumpInputs
ModelingToolkit.structural_simplify
set_default_noise_scaling

Chemistry-related functionalities

Various functionalities primarily relevant to modelling of chemical systems (but potentially also in biology).

@compound
@compounds
iscompound
components
coefficients
component_coefficients

Unit validation

validate(rx::Reaction; info::String = "")
validate(rs::ReactionSystem, info::String="")

Utility functions

symmap_to_varmap

[Spatial modelling](@id api_lattice_simulations)

The first step of spatial modelling is to create a so-called LatticeReactionSystem:

LatticeReactionSystem

The following functions can be used to querying the properties of LatticeReactionSystems:

reactionsystem
Catalyst.spatial_reactions
Catalyst.lattice
Catalyst.num_verts
Catalyst.num_edges
Catalyst.num_species
Catalyst.spatial_species
Catalyst.vertex_parameters
Catalyst.edge_parameters
Catalyst.edge_iterator
Catalyst.is_transport_system
has_cartesian_lattice
has_masked_lattice
has_grid_lattice
has_graph_lattice
grid_size
grid_dims

In addition, most accessor functions for normal ReactionSystems (such as species and parameters) works when applied to LatticeReactionSystems as well.

The following two helper functions can be used to create non-uniform parameter values.

make_edge_p_values
make_directed_edge_values

The following functions can be used to access, or change, species or parameter values stored in problems, integrators, and solutions that are based on LatticeReactionSystems.

lat_getu
lat_setu!
lat_getp
lat_setp!
rebuild_lat_internals!

Finally, we provide the following helper functions to plot and animate spatial lattice simulations.

lattice_plot
lattice_animation
lattice_kymograph