MomentClosure.jl API

propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaws=true)

Return a vector of propensity functions of all reactions in the given ReactionSystem.

Notes:

• combinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw.

get_stoichiometry(rn::ReactionSystem, smap::AbstractDict)

Return the net stoichiometry matrix using the specified mapping of species to their indices.

Notes:

• This is a modification of Catalyst.netstoichmat that is used internally to deal with reactions involving symbolic stoichiometry coefficients.
• The function also allows custom smap, so it is not limited to the default Catalyst.speciesmap ordering.
• TODO: remove once this Catalyst issue is resolved.

generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int;
                        langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn))

Given a ReactionSystem return the RawMomentEquations of the system generated up to m_order.

Notes:

• The expansion order $q$, denoted by q_order throughout the docs, is automatically determined from the given polynomial form of the propensity functions, see the tutorial and the theory section for more details on how q_order is obtained.
• if langevin=true, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation.
• combinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw. Note that this field is irrelevant using ReactionSystemMod as then the propensities are defined directly by the user.
• smap sets the variable ordering in the moment equations (which index corresponds to which species in the reaction network). By default, this is consistent with the internal system ordering accessible with Catalyst.speciesmap.

generate_raw_moment_eqs(sys::SDESystem, m_order::Int)

Given an SDESystem, return the RawMomentEquations of the system generated up to m_order.

struct RawMomentEquations <: MomentClosure.MomentEquations

Raw moment equations generated for the given system plus a number of helper parameters (used internally).

Fields

• odes: ModelingToolkit.ODESystem consisting of the time-evolution equations of raw moments.

• μ: Symbolic variables defining the raw moments.

• N: Number of species within the system.

• m_order: Order of moment equations.

• q_order: Expansion order.

• iter_all: Vector of all index combinations up to q_order.

• iter_m: Vector of all index combinations up to m_order.

• iter_q: Vector of all index combinations of order greater than m_order up to q_order.

• iter_1: Vector of index combinations of order 1.

generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0;
                            langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn))

Given a ReactionSystem return the CentralMomentEquations of the system generated up to m_order.

Notes:

• if q_order is not specified by the user, it is assumed that the reaction network contains only polynomial propensity functions and hence q_order is determined automatically as in generate_raw_moment_eqs. However, q_order must be specified if non-polynomial propensities are included. Note that the expansion order $q$ denotes the highest order of central moments which will be included in the ODEs (due to the Taylor expansion of propensity functions).
• if langevin=true, instead of the Chemical Master Equation the Chemical Langevin Equation (diffusion approximation) is considered, and the moment equations are constructed from the corresponding SDE formulation.
• combinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw. Note that this field is irrelevant using ReactionSystemMod as then the propensities are defined directly by the user.
• smap sets the variable ordering in the moment equations (which index corresponds to which species in the reaction network). By default, this is consistent with the internal system ordering accessible with Catalyst.speciesmap.

generate_central_moment_eqs(sys::SDESystem, m_order::Int, q_order::Int=0)

Given an SDESystem, return the CentralMomentEquations of the system generated up to m_order.

struct CentralMomentEquations <: MomentClosure.MomentEquations

Central moment equations generated for the given system plus a number of helper parameters (used internally).

Fields

• odes: ModelingToolkit.ODESystem consisting of the time-evolution equations of central moments.

• μ: Symbolic variables defining the means.

• M: Symbolic variables defining the central moments.

• N: Number of species within the system.

• m_order: Order of moment equations.

• q_order: Expansion order.

• iter_all: Vector of all index combinations up to q_order.

• iter_m: Vector of all index combinations up to m_order.

• iter_q: Vector of all index combinations of order greater than m_order up to q_order.

• iter_1: Vector of index combinations of order 1.

bernoulli_moment_eqs(sys::MomentEquations, binary_vars::Array{Int,1})

Given MomentEquations and an array of indices specifying the species which molecule numbers are binary variables (either 0 or 1), apply identities of Bernoulli variables to remove the redundant ODEs and return the cleaned up MomentEquations. See here for example usage.

moment_closure(sys::MomentEquations, closure::String, binary_vars::Array{Int,1}=Int[])

Given MomentEquations, apply the specified moment closure approximation and return the ClosedMomentEquations.

The supported closure options are:

• "zero"
• "normal"
• "log-normal"
• "poisson"
• "gamma"
• "derivative matching"
• "conditional gaussian"
• "conditional derivative matching"

Notes

• binary_vars must be specified for conditional closures as an array of indices of all species (as in Catalyst.speciesmap) which molecule number is a Bernoulli variable. Although not necessary for other closures, specifying binary_vars is recommended as the properties of Bernoulli variables will be used to remove the redundant moment equations and simplify the expressions, which can significantly improve numerical stability.

struct ClosedMomentEquations <: MomentClosure.MomentEquations

Closed moment equations and the corresponding closure functions.

Fields

• odes: ModelingToolkit.ODESystem consisting of the time-evolution equations of closed moments.

• closure: Dictionary of moment closure functions for each higher order moment.

• open_eqs: Original raw or central moment equations (before closure was applied).

deterministic_IC(u₀::Array{T, 1}, eqs::MomentEquations) where T<:Real

Given an array of initial molecule numbers and the corresponding moment equations, return a mapping of each moment to its initial value under deterministic initial conditions.

Notes

• The means are set to initial molecule numbers (as they take the values specified in u₀ with probability one). The higher order raw moments are products of the corresponding powers of the means whereas the higher order central moments are simply zero.
• The ordering of u₀ elements must be consistent with the ordering of species in the corresponding reaction system (can be checked with the Catalyst.speciesmap function).
• As higher-order moment functions under log-normal, gamma, derivative matching and the conditional closures involve moments raised to negative powers, setting initial molecule numbers of certain species to zeros will result in NaN errors when solving the ODEs (the specifics depend on the system at hand).

  linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0;
                               combinatoric_ratelaws = true) where T <: ReactionSystem

Given a nonlinear ReactionSystem and an equivalent linear ReactionSystem, perform the Linear Mapping Approximation (LMA) and return the corresponding linear RawMomentEquations of the system as well as a Dictionary of reaction parameter substitutions obtained using LMA that are used to generate the moment equations. See the LMA theory section for more details.

Notes:

• rn_nonlinear and rn_linear must be identical in layout in order to be interpreted correctly, and the nonlinear reactions contained in rn_nonlinear must all be linearised in rn_linear with rate coefficients updated accordingly. Although this requires a lot of manual input, automating the linearisation further is difficult due to arbitrary choices that may be mane in constructing the reaction networks.
• binary_vars must be specified for conditional closures as an array of indices of all species (as in Catalyst.speciesmap) which molecule number is a Bernoulli variable. Note that rn_nonlinear and rn_linear may internally order the species differently: binary_vars must be consistent with the ordering in the nonlinear network.
• By default the moment equations will be generated up to the order determined by the degree of nonlinearity of the nonlinear system's reactions. However, if higher order moment information is required, the optional m_order argument may be provided to increase the expansion order manually.
• combinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw.

get_raw_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an EnsembleSolution of DifferentialEquations ensemble simulation, return a Dictionary of raw moments computed up to the specified order at each time step.

Notes

• For example, the dictionary key (2,0,1) maps to an array containing the values of the raw moment $μ_{201}$ at each time step.
• It is assumed that the time steps are all at the same time point for all trajectories (i.e., fixed dt used by the integrator or values were saved using saveat, as discussed here).
• Moments are computed using Cumulants.jl internally. The naive algorithm of moment tensor calculations (naive=true) is usually faster for small systems but the proposed novel algorithm (naive=false) should be more efficient in case of many marginal variables. The block size b can also be specified for the novel algorithm and may have a significant effect on its performance.
• Only useful if higher order moments are needed: DifferentialEquations has a number of far more efficient and flexible ensemble statistics functions for means, variances and correlations, see this tutorial for more details.

get_central_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an EnsembleSolution of DifferentialEquations ensemble simulation, return a Dictionary of central moment estimates computed up to the specified order at each time step. See the notes of get_raw_moments function for more information.

get_cumulants(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)

Given an EnsembleSolution of DifferentialEquations ensemble simulation, return a Dictionary of cumulant estimates computed up to the specified order at each time step. See the notes of get_raw_moments function for more information.

get_moments_FSP(sol::ODESolution, order::Int, moment_type::String)

Given an ODESolution obtained using FiniteStateProjection.jl, return a Dictionary of moments computed up to the specified order at each time step. Here, moment_type specifies the type of moments to be computed: available options are raw, central or cumulant.

Notes

• The ODESolution represents the time-evolution of the probability density function that is the solution of the Chemical Master Equation approximated using Finite State Projection algorithms. See the documentation of FiniteStateProjection.jl for more information.