Skip to content

Latest commit

 

History

History
383 lines (343 loc) · 14.6 KB

core_api.md

File metadata and controls

383 lines (343 loc) · 14.6 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 ModelingToolkitBase.AbstractSystems, specifically a ModelingToolkitBase.System representing an ODE, SDE, or jump model.

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 ModelingToolkitBase.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
oprob = ODEProblem(rs, u₀map, tspan, parammap)
sol = solve(oprob, Tsit5())
p1 = plot(sol, title = "ODE")

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

# solve as jump process
u₀map    = [S => 999, I => 1, R => 0]
jprob = JumpProblem(rs, u₀map, tspan, parammap)
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 permit the user to supply non-reaction information to ReactionSystem created through 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 coupled_models_dsl): Allows the designation of differentials.
  • [equations](@ref coupled_models_dsl): Allows the creation of algebraic and/or differential equations.
  • [continuous_events](@ref events): Allows the creation of continuous events.
  • [discrete_events](@ref events): Allows the creation of discrete events.
  • brownians: Allows the creation of brownian processes that can be used to add noise to non-species variables.
  • poissonians: Allows the creation of poissonian processes that can be used to add jump events to non-species variables.
  • discretes: Creates discrete parameters, i.e. time-dependent parameters.
  • [combinatoric_ratelaws](@ref faq_combinatoric_ratelaws): Takes a single option (true or false), which sets whether to use combinatorial rate laws.
  • unit_checks: Takes a single option (true or false) controlling whether unit validation runs during DSL construction (false by default).

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

A ReactionSystem is an instance of a ModelingToolkitBase.AbstractSystem, and has a number of fields that can be accessed using the Catalyst API and the ModelingToolkitBase.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)

  • ModelingToolkitBase.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).
  • ModelingToolkitBase.get_ps(rn) is a vector that collects all the parameters defined within reactions in rn.
  • ModelingToolkitBase.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).
  • ModelingToolkitBase.get_iv(rn) is the independent variable used in the system (usually t to represent time).
  • ModelingToolkitBase.get_systems(rn) is a vector of all sub-systems of rn.
  • ModelingToolkitBase.get_initial_conditions(rn) is a dictionary of all the initial 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:

  • ModelingToolkitBase.unknowns(rn) returns all species and variables across the system, all sub-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.
  • ModelingToolkitBase.parameters(rn) returns all parameters across the system, all sub-systems.
  • ModelingToolkitBase.equations(rn) returns all Reactions and all Symbolics.Equations defined across the system, all sub-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 [ModelingToolkitBase 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 coupled_models).

ModelingToolkitBase.has_alg_equations
ModelingToolkitBase.alg_equations
ModelingToolkitBase.has_diff_equations
ModelingToolkitBase.diff_equations

Basic species properties

The following functions permits the querying of species properties.

isspecies
Catalyst.isconstant
Catalyst.isbc
Catalyst.isvalidreactant

Symbolic variable properties

The following function from SymbolicIndexingInterface.jl is useful for getting the name of individual symbolic variables (species, parameters, or non-species variables).

SymbolicIndexingInterface.getname

Basic reaction properties

ismassaction
dependents
dependants
substoichmat
prodstoichmat
netstoichmat
reactionrates
PhysicalScale

[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

[System-level metadata](@id api_system_metadata)

The following types and functions allow storing and retrieving system-level metadata on a ReactionSystem. These are primarily intended for use by file parsers that need to attach extra mapping data (e.g., initial condition or parameter value maps) to a system.

Catalyst.U0Map
Catalyst.ParameterMap
Catalyst.has_u0_map
Catalyst.get_u0_map
Catalyst.set_u0_map
Catalyst.has_parameter_map
Catalyst.get_parameter_map
Catalyst.set_parameter_map

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

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

ModelingToolkitBase.extend
ModelingToolkitBase.compose
Catalyst.flatten

[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 = :ode, math_delimiters = true) # hide
latexify(rn; form=:sde)
latexify(rn; form = :sde, math_delimiters = true) # hide

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 ModelingToolkitBase expressions, one can directly access the generated rate laws, and using ModelingToolkitBase tooling generate functions or Julia Exprs from them.

oderatelaw
jumpratelaw
mm
mmr
hill
hillr
hillar

Transformations and ReactionSystem Conversions

ode_model
ss_ode_model
sde_model
jump_model
hybrid_model
ModelingToolkitBase.mtkcompile
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

Catalyst.validate_units(rx::Reaction; info::String = "", warn::Bool = true)
Catalyst.validate_units(rs::ReactionSystem; info::String = "", warn::Bool = true)
Catalyst.assert_valid_units(rx::Reaction; info::String = "")
Catalyst.assert_valid_units(rs::ReactionSystem; info::String = "")
Catalyst.unit_validation_report(rx::Reaction; info::String = "")
Catalyst.unit_validation_report(rs::ReactionSystem; info::String = "")
Catalyst.UnitValidationIssue
Catalyst.UnitValidationReport
Catalyst.UnitValidationError

[Spatial modelling](@id api_dspace_simulations)

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

DiscreteSpaceReactionSystem

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

reactionsystem
Catalyst.spatial_reactions
Catalyst.dspace
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_dspace
has_masked_dspace
has_grid_dspace
has_graph_dspace
grid_size
grid_dims

In addition, most accessor functions for normal ReactionSystems (such as species and parameters) works when applied to DiscreteSpaceReactionSystems 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 DiscreteSpaceReactionSystems.

spat_getu
spat_setu!
spat_getp
spat_setp!
rebuild_spat_internals!

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

dspace_plot
dspace_animation
dspace_kymograph