MomentClosure.jl API
Model definition
MomentClosure is fully compatible with reaction network models defined using Catalyst and stored as a ModelingToolkit.ReactionSystem. Note that previously we had implemented our own ReactionSystemMod, that allowed us to consider systems containing reactions which products are independent geometrically distributed random variables. However, this is now deprecated as Catalyst has added support for such parameteric stoichiometries offering a much more complete and efficient feature set.
Basic model properties
Moreover, we include a couple of tiny extensions to the Catalyst API:
MomentClosure.propensities — Functionpropensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaws=true)Return a vector of propensity functions of all reactions in the given ReactionSystem.
Notes:
combinatoric_ratelaws=trueuses binomials in calculating the propensity functions of aReactionSystem, see the notes forModelingToolkit.jumpratelaw.
MomentClosure.get_stoichiometry — Functionget_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.netstoichmatthat is used internally to deal with reactions involving symbolic stoichiometry coefficients. - The function also allows custom
smap, so it is not limited to the defaultCatalyst.speciesmapordering. - TODO: remove once this Catalyst issue is resolved.
Moment Equations
MomentClosure.generate_raw_moment_eqs — Functiongenerate_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_orderthroughout 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 howq_orderis 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=trueuses binomials in calculating the propensity functions of aReactionSystem, see the notes forModelingToolkit.jumpratelaw. Note that this field is irrelevant usingReactionSystemModas then the propensities are defined directly by the user.smapsets 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 withCatalyst.speciesmap.
generate_raw_moment_eqs(sys::SDESystem, m_order::Int)Given an SDESystem, return the RawMomentEquations of the system generated up to m_order.
MomentClosure.RawMomentEquations — Typestruct RawMomentEquations <: MomentClosure.MomentEquationsRaw moment equations generated for the given system plus a number of helper parameters (used internally).
Fields
odes:ModelingToolkit.ODESystemconsisting 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 toq_order.iter_m: Vector of all index combinations up tom_order.iter_q: Vector of all index combinations of order greater thanm_orderup toq_order.iter_1: Vector of index combinations of order 1.
MomentClosure.generate_central_moment_eqs — Functiongenerate_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_orderis not specified by the user, it is assumed that the reaction network contains only polynomial propensity functions and henceq_orderis determined automatically as ingenerate_raw_moment_eqs. However,q_ordermust 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=trueuses binomials in calculating the propensity functions of aReactionSystem, see the notes forModelingToolkit.jumpratelaw. Note that this field is irrelevant usingReactionSystemModas then the propensities are defined directly by the user.smapsets 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 withCatalyst.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.
MomentClosure.CentralMomentEquations — Typestruct CentralMomentEquations <: MomentClosure.MomentEquationsCentral moment equations generated for the given system plus a number of helper parameters (used internally).
Fields
odes:ModelingToolkit.ODESystemconsisting 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 toq_order.iter_m: Vector of all index combinations up tom_order.iter_q: Vector of all index combinations of order greater thanm_orderup toq_order.iter_1: Vector of index combinations of order 1.
MomentClosure.bernoulli_moment_eqs — Functionbernoulli_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
MomentClosure.moment_closure — Functionmoment_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_varsmust be specified for conditional closures as an array of indices of all species (as inCatalyst.speciesmap) which molecule number is a Bernoulli variable. Although not necessary for other closures, specifyingbinary_varsis 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.
MomentClosure.ClosedMomentEquations — Typestruct ClosedMomentEquations <: MomentClosure.MomentEquationsClosed moment equations and the corresponding closure functions.
Fields
odes:ModelingToolkit.ODESystemconsisting 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).
MomentClosure.deterministic_IC — Functiondeterministic_IC(u₀::Array{T, 1}, eqs::MomentEquations) where T<:RealGiven 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 theCatalyst.speciesmapfunction). - 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).
Displaying Equations and Closures
The generated moment equations can be converted into LaTeX expressions using Latexify as:
using Latexify
-latexify(moment_eqs)A ModelingToolkit.ODESystem (saved as moment_eqs.odes) can also be passed to latexify function directly but the output will be different as we apply additional formatting to the symbolic expressions.
Given ClosedMomentEquations, the closure functions can be visualised in the same way by adding a :closure argument:
latexify(moment_eqs, :closure)Note that this will print out only those higher order moments which are found in the given moment equations. It is possible to print the closure functions of all higher order moments using print_all=true argument:
latexify(moment_eqs, :closure, print_all=true)Linear Mapping Approximation
MomentClosure.linear_mapping_approximation — Function linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0;
- combinatoric_ratelaws = true) where T <: ReactionSystemGiven 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_nonlinearandrn_linearmust be identical in layout in order to be interpreted correctly, and the nonlinear reactions contained inrn_nonlinearmust all be linearised inrn_linearwith 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_varsmust be specified for conditional closures as an array of indices of all species (as inCatalyst.speciesmap) which molecule number is a Bernoulli variable. Note thatrn_nonlinearandrn_linearmay internally order the species differently:binary_varsmust 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_orderargument may be provided to increase the expansion order manually. combinatoric_ratelaws=trueuses binomials in calculating the propensity functions of aReactionSystem, see the notes forModelingToolkit.jumpratelaw.
Stochastic Simulation Utilities
We provide provides functions for higher-order moment extraction from SSA and FSP data:
MomentClosure.get_raw_moments — Functionget_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
dtused by the integrator or values were saved usingsaveat, 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 sizebcan 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.
MomentClosure.get_central_moments — Functionget_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.
MomentClosure.get_cumulants — Functionget_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.
MomentClosure.get_moments_FSP — Functionget_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
ODESolutionrepresents 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.
![](\"https://augustinas1.github.io/MomentClosure.jl/dev/assets/nonlinear_genetic_feedback_loop_scheme.png\")
⠀","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"Here we have a two-state gene promoter which can be in either state G or G^*. As in a previous example tutorial, we interpret the gene as a distinct species modelled as a Bernoulli variable g, associating the states G and G^* with values g=1 and g=0 respectively. Protein P is produced from both gene states G and G^* (with different rates ρ_u and ρ_b) and can subsequently decay. The switch between G and G^* (or the feedback) is introduced via protein binding to the gene in state G.","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"The first step in applying the LMA is to transform the nonlinear GRN into a linear GRN. This can be done by removing the second-order reaction between P and G, so that the reversible reaction G+P undersetσ_ustackrelσ_brightleftharpoons G^* is replaced by G undersetσ_ustackrelbarσ_brightleftharpoons G^*. This equivalent linear GRN can be visualised as (adapted from [1]):","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"![](\"https://augustinas1.github.io/MomentClosure.jl/dev/assets/linear_genetic_feedback_loop_scheme.png\")
⠀","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"Note that MomentClosure cannot automate this linearisation step as the choice of how the other reactions in the network are changed due to the removal of nonlinear reactions is arbitrary. For this reason, we define both nonlinear and linear GRNs using Catalyst as follows:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using Catalyst\r\n\r\n# NOTE: both models MUST preserve the same ordering of reactions in order to detect\r\n# how the nonlinear reactions are to be transformed using LMA\r\n\r\nrn_nonlinear = @reaction_network begin\r\n @parameters σ_b σ_u ρ_b ρ_u\r\n σ_b, g + p → 0\r\n σ_u*(1-g), 0 ⇒ g + p\r\n ρ_u, g → g + p\r\n ρ_b*(1-g), 0 ⇒ p\r\n 1, p → 0\r\nend \r\n\r\nrn_linear = @reaction_network begin\r\n @parameters σ_b_LMA σ_u ρ_b ρ_u\r\n σ_b_LMA, g → 0 # typing ̄σ_b is not allowed it seems\r\n σ_u*(1-g), 0 ⇒ g\r\n ρ_u, g → g+p\r\n (ρ_b*(1-g)), 0 ⇒ p\r\n 1, p → 0\r\nend ","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"We can now apply the LMA to find the effective parameter barσ_b and generate the corresponding moment equations of the linear GRN using MomentClosure's linear_mapping_approximation:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using MomentClosure\r\n\r\n# NOTE: we have to provide the indices of binary variables in the system as they are ordered in the *nonlinear* GRN.\r\n# The distinction here between linear and nonlinear GRNs is important as in some cases the internal ordering of variables of the two Catalyst models can differ\r\n@parameters t\r\n@variables g(t)\r\nbinary_vars = [speciesmap(rn_nonlinear)[g]]\r\n\r\nLMA_eqs, effective_params = linear_mapping_approximation(rn_nonlinear, rn_linear, binary_vars, combinatoric_ratelaws=false)\r\ndisplay(effective_params)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"OrderedDict{Any, Any} with 1 entry:\r\n σ_b_LMA => σ_b*μ₁₁(t)*(μ₁₀(t)^-1)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"We can also print out the moment equations:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using Latexify\r\nlatexify(LMA_eqs)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"beginalign*\r\nfracdmu_10dt = sigma_u - sigma_b mu_11 - sigma_u mu_10 \r\nfracdmu_01dt = rho_b + rho_u mu_10 - mu_01 - rho_b mu_10 \r\nfracdmu_11dt = rho_u mu_10 + sigma_u mu_01 - mu_11 - sigma_u mu_11 - sigma_b mu_10^-1 mu_11^2 \r\nfracdmu_02dt = rho_b + rho_u mu_10 + 2 rho_b mu_01 + 2 rho_u mu_11 + mu_01 - 2 mu_02 - rho_b mu_10 - 2 rho_b mu_11\r\nendalign*","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"Note that the results agree with Eqs. (1) and (2) (after a corresponding substitution) given in [1]. The moment equations are already closed, so we can solve them numerically and plot the mean protein number over time:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using OrdinaryDiffEq, Sundials, Plots\r\n\r\n# [g, p] as in `species(rn_nonlinear)`\r\nu₀ = [1.0, 0.001]\r\np = [0.004, 0.25, 25.0, 60.0]\r\ntspan = (0., 15.)\r\ndt = 0.1\r\n\r\nu₀map = deterministic_IC(u₀, LMA_eqs)\r\noprob_LMA = ODEProblem(LMA_eqs, u₀map, tspan, p)\r\nsol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)\r\n\r\nplot(sol_LMA, idxs=[2], label=\"LMA\", ylabel=\"⟨p⟩\", xlabel=\"time\", fmt=\"svg\")","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"(Image: LMA feedback loop mean protein number)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"To compare the LMA result to the true moment dynamics we use the Finite State Projection (FSP) algorithm implemented in FiniteStateProjection.jl. FSP can be much more efficient than SSA in accurately approximating the full time-dependent probability distribution of the given chemical system when the system's state space is small (hence it is particularly effective for the nonlinear GRN here). Having obtained the FSP solution, we can extract the time-evolution of moments using get_moments_FSP function and finally compare it to the LMA prediction. This can all be done as follows:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using FiniteStateProjection\r\n\r\nfsp_sys = FSPSystem(rn_nonlinear, combinatoric_ratelaw=false)\r\n# Truncate the state space of the system\r\n# The gene has two states (G or G*) whereas we consider protein number from 0 to 100\r\nstate_space = [2, 101]\r\n\r\n# The initial condition is the matrix of probabilities representing the state of the system\r\n# We assume zero protein and the gene to be in the state G, hence the probability of this\r\n# specific state should be set to 1 initially\r\nu0 = zeros(state_space...)\r\nu0[2, 1] = 1.0\r\n\r\n# construct an ODE problem from the FSPSystem and solve it\r\nfsp_prob = ODEProblem(fsp_sys, u0, tspan, p)\r\nsol_FSP = solve(fsp_prob, CVODE_BDF(), saveat=dt)\r\n\r\n# extract the 1st order raw moments from the FSP solution\r\nμ_FSP = get_moments_FSP(sol_FSP, 1, \"raw\")\r\nplot!(sol_FSP.t, μ_FSP[(0,1)], label=\"FSP\", legend=:bottomright)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"(Image: LMA and FSP feedback loop mean protein number)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"Using the LMA as implemented in MomentClosure, we can generate and solve the closed moment equations for any LMA-suitable nonlinear chemical reaction network, given that its linear equivalent is also provided. Proceeding further to obtain an approximate time-dependent probability distribution of the nonlinear network is more involved as only a handful of closed-form solutions of linear networks are known. Moreover, simply computing the solutions can be a challenge due to their complicated analytical form. For these reasons, MomentClosure does not provide an automated approach to computing the probability distributions of the relevant nonlinear systems using LMA—these steps have to performed manually on a case-by-case basis. Nevertheless, below we demonstrate how this could be done for the nonlinear GRN we have been considering so far.","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"The CME of the nonlinear feedback loop can be solved exactly as described in [1] and, for completeness, we write down the main result here. The probability distribution of the number of proteins n_p at time t is given by (assuming the initial conditions to be zero protein in state G):","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"P left( n_p t right) = frac1n_p fracd^n_pdw^n_pleft( G_0(w t) + G_1(w t) right) _w=-1","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"where","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"beginalign*\r\nG_0(w t) = exp(ρ_b w) left f(we^-t) (-rho_Delta w)^1-Sigma M(1-barsigma_b 2-Sigma -rho_Delta w) \r\n+ g(we^-t)M(1+sigma_uSigma-rho_Delta w) right\r\nG_1(w t) = sigma_u^-1 exp(rho_b w) left -sigma_u f(we^-t)(-rho_Delta w)^1-Sigma M(-barsigma_b 2-Sigma -rho_Delta w) \r\n+barsigma_b g(we^-t) M(sigma_u Sigma -rho_Delta w) right\r\nendalign*","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"where w=z-1 and the probability generating functions are given by G_i(z t) = sum_n_p=0^infty z^n_p P_i(n_p t). The function M(cdotcdotcdot) stands for the Kummer's (confluent hypergeometric) function and we also use the following definitions:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"beginalign*\r\n rho_Delta = rho_b - rho_u \r\n Sigma = sigma_u + barsigma_b + 1 \r\n f(w) = fracbarsigma_bSigma-1(-rho_Delta w)^Sigma-1e^-rho_u wM(sigma_u Sigma -rho_Delta w) \r\n g(w) = fracsigma_uSigma-1e^-rho_u w M(-barsigma_b 2-Sigma -rho_Delta w)\r\nendalign*","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"One can observe that having a closed-form solution does not make life that much easier as now we have to find a way how to numerically evaluate these complicated mathematical expressions in a reasonably efficient manner. The main problem here is that in order to construct the probability distribution we have to take arbitrarily high n^th order derivatives. Doing that symbolically is possible in principle but can be extremely computationally expensive as the symbolic expressions will be expanding exponentially in size with each higher derivative order: keeping them in memory is expensive and the same terms would have to be evaluated repeatedly many times. Similarly, naively nesting numerical derivative operations using an automatic differentiation framework is possible but still many redundant operations would need to be performed.","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"As pointed out by the authors in [1], a more efficient way is to expand the generating functions G_0(w t) + G_1(w t) around w=-1 as then the n^th Taylor series coefficient will be equivalent to the probability of having n proteins at time t. This approach can be particularly effective: recursion relations for the Taylor coefficients of various elementary functions can be defined and in turn utilised to minimise the amount of redundant work which has to be performed when Taylor expanding any given mathematical expression. Such higher-order automatic differentiation framework in Julia is provided by TaylorSeries.jl.","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"However, TaylorSeries only supports elementary function operations at the time and hence evaluating the Kummer's function M(cdotcdotcdot) requires some more work (these specialised numerics are readily available in more established scientific computing frameworks such as Mathematica but there's no fun in that). We can extend the TaylorSeries framework by constructing a function t_pFq that implements a recurrence relation between the Taylor coefficients for the generalized hypergeometric function pFq as defined in HypergeometricFunctions.jl. This can be done as follows (note that our construction is valid only for a single-variable Taylor series Taylor1):","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"using TaylorSeries, HypergeometricFunctions\r\nusing HypergeometricFunctions: pFqweniger\r\n\r\n# please let me know if a simpler and more efficient way to do this exists!\r\nfunction t_pFq(α::AbstractVector, β::AbstractVector, a::Taylor1)\r\n order = a.order\r\n aux = pFqweniger(α, β, constant_term(a))\r\n c = Taylor1(aux, order)\r\n\r\n iszero(order) && return c\r\n\r\n coeffs = t_pFq(α.+1, β.+1, Taylor1(a[0:end-1], a.order-1))\r\n factor = prod(α)/prod(β)\r\n for k in 1:length(a)-1\r\n c[k] = sum(i * a[i] * coeffs[k-i] for i in 1:k) * factor / k\r\n end\r\n\r\n return c\r\n\r\nend","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"We can now move on with the LMA procedure to obtain the probability distribution of the nonlinear GRN for the number of proteins at a chosen time point t=4. Step 4 of the LMA involves obtaining the time-independent effective parameter barsigma_b which is given by its time-averaged value:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"# calculate the raw moments up to time t at a fine temporal resolution\r\nT = 4.0\r\ntspan = (0., T)\r\ndt = 0.001\r\noprob_LMA = remake(oprob_LMA, tspan=tspan)\r\nsol_LMA = solve(oprob_LMA, CVODE_BDF(), saveat=dt)\r\n\r\n# rebuild the symbolic expression for the effective parameter as a function of raw moments\r\nμ_sym = [LMA_eqs.odes.states...]\r\np_sub = Dict.(Pair.(LMA_eqs.odes.ps, p)...)\r\navg_σ_b_sym = collect(values(effective_params))[1]\r\nfn = build_function(substitute(avg_σ_b_sym, p_sub), μ_sym)\r\navg_σ_b = eval(fn)\r\n# evaluate the time-averaged value of the effective parameter\r\n@time σ_b_avg = sum(avg_σ_b.(sol_LMA.u)) * dt / T","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"We proceed with the very last steps of the LMA to obtain the probability distribution:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"# need higher-precision numerics as Float64 can be unstable here due to very small numbers\r\n# DoubleFloats is sufficient for this example and much more efficient than BigFloat\r\nusing DoubleFloats\r\n\r\n# define the numerical values of the parameters\r\nσ_u = p[2]; ρ_b = p[3]; ρ_u = p[4]\r\nΣ = 1 + σ_b_avg + σ_u\r\nρ_Δ = ρ_b - ρ_u\r\n\r\nn = 100 # expansion order (or max protein number to evaluate)\r\nw₀ = -1 # value around which to expand\r\n\r\n# compute the Taylor expansion (note the use of Double64)\r\nw = w₀ + Taylor1(Double64, n)\r\n@time f = σ_b_avg/(Σ-1)*exp(-T*(Σ-1))*exp(-ρ_u*w*exp(-T))*t_pFq([σ_u], [Σ], -ρ_Δ*w*exp(-T))\r\n@time g = σ_u/(Σ-1)*exp(-ρ_u*w*exp(-T))*t_pFq([-σ_b_avg], [2-Σ], -ρ_Δ*w*exp(-T))\r\n\r\n@time G00 = exp(ρ_b*w)*(f * t_pFq([1-σ_b_avg], [2-Σ], -ρ_Δ*w) +\r\n g * t_pFq([1+σ_u], [Σ], -ρ_Δ*w) )\r\n\r\n@time G11 = σ_u^(-1) * exp(ρ_b*w) * (-σ_u*f*t_pFq([-σ_b_avg], [2-Σ], -ρ_Δ*w) +\r\n σ_b_avg*g*t_pFq([σ_u], [Σ], -ρ_Δ*w))\r\n\r\nprobs = (G00+G11).coeffs\r\n\r\n# check that the probability distribution is more or less normalised to 1\r\n# need higher numerical precision if not\r\nisapprox(sum(probs), 1.0, rtol=1e-2)","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"0.009289 seconds (41.03 k allocations: 1.603 MiB)\r\n0.009126 seconds (41.03 k allocations: 1.603 MiB)\r\n0.067889 seconds (85.06 k allocations: 3.613 MiB)\r\n0.066194 seconds (85.07 k allocations: 3.618 MiB)\r\ntrue","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"Finally, we can plot the distribution and compare it to the FSP result:","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"plot(0:n, probs, xlabel=\"n\", ylabel=\"P(n, t=4)\", label=\"LMA\", fmt=\"svg\")\r\n# plot the FSP probability of protein number by marginalising over the gene states\r\nplot!(0:n, sum(sol_FSP[:, 41], dims=1)'[1:n+1], label=\"FSP\")","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"(Image: LMA FSP feedback loop distribution)","category":"page"},{"location":"tutorials/LMA_example/#References","page":"LMA Example","title":"References","text":"","category":"section"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"[1]: Z. Cao and R. Grima, \"Linear mapping approximation of gene regulatory networks with stochastic dynamics\", Nature Communications 9, 3305 (2018). https://doi.org/10.1038/s41467-018-05822-0","category":"page"},{"location":"tutorials/time-dependent_propensities/#time-dependent-propensities","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"","category":"section"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"Here we consider an example of using time-dependent propensities in modelling chemical reaction networks using MomentClosure and DifferentialEquations. Following Schnoerr et al. (2015), we modify the Brusellator introduced earlier to include entrainment:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"beginalign*\r\n2X + Y stackrelc_1rightarrow 3X \r\nX stackrelc_2(t)rightarrow Y \r\n undersetc_4stackrelc_3rightleftharpoons X\r\nendalign*","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"The rate constant of the second reaction is now a sinusoidal function of time:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"beginalign*\r\n c_2(t) = c_2^0left(1+frac12sin(omega t) right) quad textif t leq tau \r\n c_2(t) = c_2^0 quad textif t tau\r\nendalign*","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"where c_2^0 is a fixed value and the sinusoidal modulation with frequency omega is switched off after a certain time tau. We can define the model and generate the corresponding moment equations as follows:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"using Catalyst, MomentClosure, Latexify\r\n\r\nrn = @reaction_network begin\r\n @parameters c₁ c₂ c₃ c₄ Ω ω τ\r\n (c₁/Ω^2), 2X + Y → 3X\r\n (c₂*(1+0.5*sin(ω*(t<τ)*t))), X → Y\r\n (c₃*Ω, c₄), 0 ↔ X\r\nend\r\n\r\nraw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)\r\nlatexify(raw_eqs)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"beginalign*\r\nfracdmu_10dt = c_3 Omega + c_1 mu_21 Omega^-2 - c_2 mu_10 - c_4 mu_10 - c_1 mu_11 Omega^-2 - 05 c_2 mu_10 sinleft( t omega left( t tau right) right) \r\nfracdmu_01dt = c_2 mu_10 + c_1 mu_11 Omega^-2 + 05 c_2 mu_10 sinleft( t omega left( t tau right) right) - c_1 mu_21 Omega^-2 \r\nfracdmu_20dt = c_2 mu_10 + c_3 Omega + c_4 mu_10 + 2 c_1 mu_31 Omega^-2 + 05 c_2 mu_10 sinleft( t omega left( t tau right) right) + 2 c_3 Omega mu_10 - 2 c_2 mu_20 - 2 c_4 mu_20 - c_1 mu_11 Omega^-2 - c_1 mu_21 Omega^-2 - c_2 mu_20 sinleft( t omega left( t tau right) right) \r\nfracdmu_11dt = c_2 mu_20 + c_1 mu_11 Omega^-2 + c_1 mu_22 Omega^-2 + c_3 Omega mu_01 + 05 c_2 mu_20 sinleft( t omega left( t tau right) right) - c_2 mu_10 - c_2 mu_11 - c_4 mu_11 - c_1 mu_12 Omega^-2 - c_1 mu_31 Omega^-2 - 05 c_2 mu_10 sinleft( t omega left( t tau right) right) - 05 c_2 mu_11 sinleft( t omega left( t tau right) right) \r\nfracdmu_02dt = c_2 mu_10 + c_1 mu_21 Omega^-2 + c_2 mu_11 sinleft( t omega left( t tau right) right) + 2 c_2 mu_11 + 2 c_1 mu_12 Omega^-2 + 05 c_2 mu_10 sinleft( t omega left( t tau right) right) - c_1 mu_11 Omega^-2 - 2 c_1 mu_22 Omega^-2\r\nendalign*","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"We can now easily close the moment equations using normal closure and solve the resulting system of ODEs:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"using OrdinaryDiffEq, Plots\r\n\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"normal\")\r\n\r\n# parameter values [c₁, c₂, c₃, c₄, Ω, ω, τ]\r\np = [0.9, 2., 1., 1., 5., 1., 40.]\r\n\r\n# initial molecule numbers of species [X, Y]\r\nu₀ = [1., 1.]\r\n\r\n# deterministic initial conditions\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\n\r\n# time interval to solve one on\r\ntspan = (0., 100.)\r\n\r\n# convert the closed raw moment equations into a DifferentialEquations ODEProblem\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\n\r\n# solve using Tsit5 solver\r\nsol = solve(oprob, Tsit5(), saveat=0.2)\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"(Image: Time-dependent Brusselator normal)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"It would be great to compare our results to the true dynamics. Using DifferentialEquations, we can run a modified SSA taking into account the time-dependent propensity functions (VariableRateJumps). This requires some care and is done by combining an ODEProblem with a JumpProblem:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"using JumpProcesses\r\n\r\n# convert ReactionSystem into JumpSystem\r\njsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)\r\n\r\n# Define an ODEProblem to integrate between reaction events\r\n# note that u₀ elements must be Floats (otherwise throws an error)\r\nf(du,u,p,t) = du .= 0\r\noprob = ODEProblem(f, u₀, tspan, p)\r\n\r\n# Create a modified SSA problem that now correctly incorporates VariableRateJumps\r\njprob = JumpProblem(jsys, oprob, Direct())","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"Finally, we can define a corresponding EnsembleProblem to simulate multiple SSA trajectories. However, the saveat argument does not work with VariableRateJumps (a known bug): without it, the solution is saved at each reaction event, in turn generating large data arrays that can get extremely memory-intensive when many trajectories are considered. Our workaround is to simply modify the output_func in EnsembleProblem so that each SSA trajectory is saved to the output array only at the specified timepoints (albeit a lot of garbage collection must be done):","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"# timestep at which the solution data is saved\r\ndt = 0.2\r\n# the corresponding time iterator (0:0.2:100 in our case)\r\nts = tspan[1]:dt:tspan[2]\r\n# save data for each trajectory only at the specified timepoints (interpolating the ODESolution)\r\nfout = (sol, i) -> (sol(ts), false)\r\nensembleprob = EnsembleProblem(jprob, output_func=fout)\r\n\r\n# simulate 10000 SSA trajectories (can get very slow...)\r\n@time sol_SSA = solve(ensembleprob, Tsit5(), trajectories=10000)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"2975.249760 seconds (1.07 G allocations: 105.883 GiB, 81.30% gc time, 0.02% compilation time)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"Finally, we can compute the mean SSA trajectories and compare to the moment closure estimates:","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"using SciMLBase.EnsembleAnalysis\r\n\r\nmeans_SSA = timeseries_steps_mean(sol_SSA)\r\nplot!(means_SSA.t, [means_SSA[1,:], means_SSA[2,:]], lw=1.5, labels=[\"SSA μ₁₀(t)\" \"SSA μ₀₁(t)\"], linestyle=:dash,\r\n linecolor=[1 2], background_color_legend=nothing, legend=:bottomright)","category":"page"},{"location":"tutorials/time-dependent_propensities/","page":"Time-dependent Propensity Functions","title":"Time-dependent Propensity Functions","text":"(Image: Time-dependent Brusselator SSA)","category":"page"},{"location":"tutorials/derivative_matching_example/#derivative-matching-example","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"","category":"section"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"In this example, we aim to reproduce a number of results from the original derivative matching paper by Singh and Hespanha [1]. We consider the bimolecular reaction system given by:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nX_1 stackrelc_1rightarrow 2X_1 + X_2 \r\nX_1 + X_2stackrelc_2rightarrow X_2\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"The reaction network and its parameters can be defined as follows:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"using Catalyst\r\n\r\nrn = @reaction_network begin\r\n @parameters c₁ c₂\r\n (c₁), x₁ → 2x₁+x₂\r\n (c₂), x₁+x₂ → x₂\r\nend\r\n\r\n# parameter values\r\np = [1.0, 1.0]\r\n# initial conditions\r\nu0 = [20, 10]\r\n# time interval to solve on\r\ntspan = (0., 0.5)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"We are interested in extracting the time-evolution of a specific third order cumulant, κ_03, using second and third order moment expansions with derivative matching, and comparing the obtained estimates to the SSA prediction.","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Let's start with a second order moment expansion and print out the third-order moment closure functions obtained with derivative matching:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"using MomentClosure, Latexify\r\n\r\neqs2 = generate_raw_moment_eqs(rn, 2)\r\ndm2_eqs = moment_closure(eqs2, \"derivative matching\")\r\nlatexify(dm2_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_30 = mu_10^-3 mu_20^3 \r\nmu_21 = mu_20 mu_01^-1 mu_10^-2 mu_11^2 \r\nmu_12 = mu_02 mu_01^-2 mu_10^-1 mu_11^2 \r\nmu_03 = mu_01^-3 mu_02^3\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Note that all closure functions are consistent with the ones shown in Table II of [1]. We can then move on to solving the generated system of moment ODEs:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"using OrdinaryDiffEq\r\n\r\nu0map = deterministic_IC(u0, dm2_eqs) # assuming deterministic initial conditions\r\noprob = ODEProblem(dm2_eqs, u0map, tspan, p)\r\ndm2_sol = solve(oprob, Tsit5(), saveat=0.01)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Now the question is how can we extract the time evolution of the cumulant kappa_03. Firstly, note that using the standard moment relationships it can be expressed in terms of raw moments as:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nkappa_03 = 2 mu_01^3 - 3mu_02mu_01 + mu_03\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"As we were solving for moments up to second order, we do not have any direct information on the third order moment mu_03. Nevertheless, we can manually approximate it using the corresponding closure function given above, i.e., mu_03 = mu_01^-3 mu_02^3. The time trajectories of mu_01 and mu_02 can be extracted from dm2_sol and their order in the array can be checked with:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"dm2_eqs.odes.states","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"5-element Array{Term{Real,Nothing},1}:\r\n μ₁₀(t)\r\n μ₀₁(t)\r\n μ₂₀(t)\r\n μ₁₁(t)\r\n μ₀₂(t)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Finally, we can combine all the steps to obtain the kappa_03 estimate:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"μ₀₁ = dm2_sol[2, :]\r\nμ₀₂ = dm2_sol[5, :]\r\nμ₀₃ = (μ₀₁ .^(-3)) .* (μ₀₂ .^3)\r\ndm2_κ₀₃ = 2 .* μ₀₁ .^3 .- 3 .* μ₀₂ .* μ₀₁ .+ μ₀₃","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Next we consider a third order moment expansion:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"eqs3 = generate_raw_moment_eqs(rn, 3)\r\ndm3_eqs = moment_closure(eqs3, \"derivative matching\")\r\nlatexify(dm3_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_40 = mu_10^4 mu_20^-6 mu_30^4 \r\nmu_31 = mu_01 mu_30 mu_10^3 mu_11^-3 mu_20^-3 mu_21^3 \r\nmu_22 = mu_01^2 mu_02^-1 mu_10^2 mu_11^-4 mu_12^2 mu_20^-1 mu_21^2 \r\nmu_13 = mu_03 mu_10 mu_01^3 mu_02^-3 mu_11^-3 mu_12^3 \r\nmu_04 = mu_01^4 mu_02^-6 mu_03^4\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"As expected, the closure functions agree with those given in Table III of [1]. We again check the order of variables","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"dm3_eqs.odes.states","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"9-element Array{Term{Real,Nothing},1}:\r\n μ₁₀(t)\r\n μ₀₁(t)\r\n μ₂₀(t)\r\n μ₁₁(t)\r\n μ₀₂(t)\r\n μ₃₀(t)\r\n μ₂₁(t)\r\n μ₁₂(t)\r\n μ₀₃(t)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"and solve the moment equations, computing the required cumulant:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"u0map = deterministic_IC(u0, dm3_eqs)\r\noprob = ODEProblem(dm3_eqs, u0map, tspan, p)\r\ndm3_sol = solve(oprob, Tsit5(), saveat=0.01, abstol=1e-8, reltol=1e-8)\r\n\r\nμ₀₁ = dm3_sol[2,:]\r\nμ₀₂ = dm3_sol[5,:]\r\nμ₀₃ = dm3_sol[9,:]\r\ndm3_κ₀₃ = 2 .* μ₀₁ .^ 3 - 3 .* μ₀₂ .* μ₀₁ .+ μ₀₃","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Note that we could have also obtained kappa_03 estimate in an easier way by using central moment equations, as third order central moments are equal to the corresponding third order cumulants:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"central_eqs3 = generate_central_moment_eqs(rn, 3)\r\ndm3_central_eqs = moment_closure(central_eqs3, \"derivative matching\")\r\n\r\nu0map = deterministic_IC(u0, dm3_central_eqs)\r\noprob = ODEProblem(dm3_central_eqs, u0map, tspan, p)\r\ndm3_central_sol = solve(oprob, Tsit5(), saveat=0.01, abstol=1e-8, reltol=1e-8)\r\n\r\n# check that the two estimates are equivalent\r\ndm3_κ₀₃ ≈ dm3_central_sol[9,:]","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"true","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"The last ingredient we need for a proper comparison between the second and third order moment expansions is a reference value predicted by the SSA. We can simulate 10^5 SSA trajectories as follows:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"using JumpProcesses\r\n\r\ndprob = DiscreteProblem(rn, u0, tspan, p)\r\njprob = JumpProblem(rn, dprob, Direct(), save_positions=(false, false))\r\n\r\nensembleprob = EnsembleProblem(jprob)\r\n@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.01, trajectories=100000)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"3.874558 seconds (7.45 M allocations: 714.845 MiB, 46.60% gc time)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"The time evolution of kappa_03 can be extracted from SSA data using the get_cumulants function:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"ssa_κ₀₃ = get_cumulants(sol_SSA, 3)[0, 3]","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"Finally, we plot the results:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"using Plots, LaTeXStrings\r\n\r\nplot(dm2_sol.t, dm2_κ₀₃, lw=2, label=\"2nd order DM\")\r\nplot!(dm2_sol.t, dm3_κ₀₃, lw=2, label=\"3rd order DM\")\r\nplot!(dm2_sol.t, ssa_κ₀₃, lw=2, label=\"SSA\")\r\nplot!(ylabel=L\"\\kappa_{03}\", xlabel=L\"t\", guidefontsize=14, legend=:bottomright)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"(Image: Derivative Matching cumulant)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"We observe that the third order moment truncation using derivative matching performs significantly better than the second order truncation, accurately matching the true SSA prediction (consistent with the figure in [1]).","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"It is also interesting to note that the closure functions of third order moments obtained using derivative matching and log-normal closures are equivalent. We can see that it is indeed the case by printing out the log-normal closure functions","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"ln2_eqs = moment_closure(eqs2, \"log-normal\")\r\nlatexify(ln2_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_30 = mu_10^-3 mu_20^3 \r\nmu_21 = mu_20 mu_01^-1 mu_10^-2 mu_11^2 \r\nmu_12 = mu_02 mu_01^-2 mu_10^-1 mu_11^2 \r\nmu_03 = mu_01^-3 mu_02^3\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"and the corresponding derivative matching functions obtained previously:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"latexify(dm2_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_30 = mu_10^-3 mu_20^3 \r\nmu_21 = mu_20 mu_01^-1 mu_10^-2 mu_11^2 \r\nmu_12 = mu_02 mu_01^-2 mu_10^-1 mu_11^2 \r\nmu_03 = mu_01^-3 mu_02^3\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"However, the equivalence holds only for third order moments. For example, the closure functions of fourth order moments differ:","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"ln3_eqs = moment_closure(eqs3, \"log-normal\")\r\nlatexify(ln3_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_40 = mu_10^-8 mu_20^6 \r\nmu_31 = mu_01^-2 mu_10^-6 mu_11^3 mu_20^3 \r\nmu_22 = mu_02 mu_20 mu_01^-4 mu_10^-4 mu_11^4 \r\nmu_13 = mu_01^-6 mu_02^3 mu_10^-2 mu_11^3 \r\nmu_04 = mu_01^-8 mu_02^6\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"latexify(dm3_eqs, :closure, print_all=true)","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"beginalign*\r\nmu_40 = mu_10^4 mu_20^-6 mu_30^4 \r\nmu_31 = mu_01 mu_30 mu_10^3 mu_11^-3 mu_20^-3 mu_21^3 \r\nmu_22 = mu_01^2 mu_02^-1 mu_10^2 mu_11^-4 mu_12^2 mu_20^-1 mu_21^2 \r\nmu_13 = mu_03 mu_10 mu_01^3 mu_02^-3 mu_11^-3 mu_12^3 \r\nmu_04 = mu_01^4 mu_02^-6 mu_03^4\r\nendalign*","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"This difference is expected as fourth order moments using log-normal closure are expressed exclusively in terms of first and second order moments, whereas derivative matching additionally incorporates third order moment information. Singh and Hespanha [1, 2] elaborate on this point: as there is no unique way to define higher order moment closure functions for log-normal distribution, both closures are consistent with the assumption that the population is jointly log-normally distributed. We urge the reader to consult the mentioned papers [1, 2] and the references therein for a more complete discussion comparing the two approaches.","category":"page"},{"location":"tutorials/derivative_matching_example/#References","page":"Derivative Matching Example","title":"References","text":"","category":"section"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"[1]: A. Singh and J. P. Hespanha, \"Lognormal Moment Closures for Biochemical Reactions\", in Proceedings of the 45th IEEE Conference on Decision and Control, ISSN:0191-2216 (Dec. 2006), pp. 2063–2068. https://doi.org/10.1109/CDC.2006.376994","category":"page"},{"location":"tutorials/derivative_matching_example/","page":"Derivative Matching Example","title":"Derivative Matching Example","text":"[2]: A. Singh and J. P. Hespanha, \"Approximate Moment Dynamics for Chemically Reacting Systems\", IEEE Transactions on Automatic Control 56, 414–418 (2011). https://doi.org/10.1109/TAC.2010.2088631","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/#geometric-and-conditional","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"","category":"section"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"In this tutorial, we demonstrate using Catalyst how to define chemical systems involving reactions which products are geometrically distributed random variables. As an example, we consider an autoregulatory (repressive) genetic feedback loop where proteins are expressed in bursts with a geometric burst size distribution, as described by Soltani et al. [1]. Moreover, the state of the gene is modelled as a binary variable—we demonstrate how conditional derivative matching and conditional gaussian closures can be used to approximate such systems, in turn reproducing a number of results from Ref. [1].","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We consider a negative feedback loop described by the following reactions:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nG^* stackrelk_onrightarrow G \r\nG + stackrelk_off*P^2rightarrow G^* \r\nG stackrelk_prightarrow G + mP \r\nP stackrelgamma_prightarrow \r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"A gene in the network switches between ON (G) and OFF (G^*) states: proteins are produced in the transcriptionally active ON state but the gene can be turned OFF by two protein molecules binding to the promoter region and thus blocking transcription (proteins decay at a constant rate irrespective of the gene state). Note that the gene state can be interpreted as a distinct species that have either zero or one copy number per cell. In other words, it is a Bernoulli random variable: 0 in the OFF state and 1 in the ON state.","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"The transcription (mRNA) dynamics are not modelled explicitly in this gene circuit. Instead, under the assumption of fast mRNA decay, proteins are taken to be produced in bursts of size m, where m is a random variable sampled from the geometric distribution phi(m) = b^m(1+b)^m+1 (with b denoting the mean burst size) [2]. Chemical reaction networks that include reactions which products are independent geometrically distributed random variables (or other symbolic variables) can be defined using Catalyst as demonstrated in this tutorial. Note that the referenced tutorial refers back to this tutorial because in the ancient times Catalyst did not provided such functionality and the previous iteration of MomentClosure implemented a ReactionSystemMod type that offered limited support for such systems (now deprecated as Catalyst does it better in a unified API).","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"Using Catalyst, our gene network model can be constructed as follows:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"# load all the packages we will need\r\nusing MomentClosure, Catalyst, Distributions, JumpProcesses, DiffEqBase, OrdinaryDiffEq, DiffEqBase.EnsembleAnalysis, Plots, Latexify\r\n\r\n# Proteins are produced in bursts of size m, where m is a geometric random variable with mean b.\r\n# Note that if b is the mean burst size, then p = 1/(1+b). \r\n# Implemented by first registering the distribution with Symbolics\r\n@register_symbolic Distributions.Geometric(b)\r\n@parameters b\r\nm = rand(Distributions.Geometric(1/(1+b)))\r\n\r\nrn = @reaction_network begin\r\n @parameters k_on k_off k_p γ_p\r\n k_on*(1-g), 0 --> g # G* -> G\r\n k_off*P^2, g --> 0 # G -> G*\r\n k_p, g --> g + $m*P # G -> G + mP, m ~ Geometric(p)\r\n γ_p, P --> 0 # P -> ∅\r\nend","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We can now generate the raw moment equations up to third order:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"eqs = generate_raw_moment_eqs(rn, 3)\r\nlatexify(eqs)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nfracdmu_1 0dt = k_on - k_off mu_1 2 - k_on mu_1 0 \r\nfracdmu_0 1dt = b k_p mu_1 0 - gamma_p mu_0 1 \r\nfracdmu_2 0dt = k_on + k_off mu_1 2 + k_on mu_1 0 - 2 k_on mu_2 0 - 2 k_off mu_2 2 \r\nfracdmu_1 1dt = k_on mu_0 1 + b k_p mu_2 0 - k_on mu_1 1 - gamma_p mu_1 1 - k_off mu_1 3 \r\nfracdmu_0 2dt = gamma_p mu_0 1 + b k_p mu_1 0 + 2 b k_p mu_1 1 + 2 k_p b^2 mu_1 0 - 2 gamma_p mu_0 2 \r\nfracdmu_3 0dt = k_on + 3 k_off mu_2 2 + 2 k_on mu_1 0 - k_off mu_1 2 - 3 k_off mu_3 2 - 3 k_on mu_3 0 \r\nfracdmu_2 1dt = k_on mu_0 1 + k_on mu_1 1 + k_off mu_1 3 + b k_p mu_3 0 - 2 k_on mu_2 1 - gamma_p mu_2 1 - 2 k_off mu_2 3 \r\nfracdmu_1 2dt = k_on mu_0 2 + gamma_p mu_1 1 + b k_p mu_2 0 + 2 k_p b^2 mu_2 0 + 2 b k_p mu_2 1 - k_off mu_1 4 - k_on mu_1 2 - 2 gamma_p mu_1 2 \r\nfracdmu_0 3dt = b k_p mu_1 0 + 3 gamma_p mu_0 2 + 3 b k_p mu_1 1 + 3 b k_p mu_1 2 + 6 k_p b^2 mu_1 0 + 6 k_p b^3 mu_1 0 + 6 k_p b^2 mu_1 1 - gamma_p mu_0 1 - 3 gamma_p mu_0 3\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"A lot of information in this system of ODEs is redundant as the gene state is a Bernoulli variable that (in our case) has the following properties:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nμ_j0 = μ_10 quad j in 2 3 dotsc \r\nμ_jk = μ_1k quad jk in 1 2 3 dotsc\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We can clean the moment equations accordingly by specifying which molecule numbers are Bernoulli variables and using the bernoulli_moment_eqs function:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"# specify the indices of species which molecule numbers are binary\r\nbinary_vars = [1]\r\n# simplify the moment equations using properties of Bernoulli variables\r\nclean_eqs = bernoulli_moment_eqs(eqs, binary_vars)\r\nlatexify(clean_eqs)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nfracdmu_1 0dt = k_on - k_off mu_1 2 - k_on mu_1 0 \r\nfracdmu_0 1dt = b k_p mu_1 0 - gamma_p mu_0 1 \r\nfracdmu_1 1dt = k_on mu_0 1 + b k_p mu_1 0 - k_on mu_1 1 - gamma_p mu_1 1 - k_off mu_1 3 \r\nfracdmu_0 2dt = gamma_p mu_0 1 + b k_p mu_1 0 + 2 b k_p mu_1 1 + 2 k_p b^2 mu_1 0 - 2 gamma_p mu_0 2 \r\nfracdmu_1 2dt = k_on mu_0 2 + gamma_p mu_1 1 + b k_p mu_1 0 + 2 b k_p mu_1 1 + 2 k_p b^2 mu_1 0 - k_off mu_1 4 - k_on mu_1 2 - 2 gamma_p mu_1 2 \r\nfracdmu_0 3dt = b k_p mu_1 0 + 3 gamma_p mu_0 2 + 3 b k_p mu_1 1 + 3 b k_p mu_1 2 + 6 k_p b^2 mu_1 0 + 6 k_p b^3 mu_1 0 + 6 k_p b^2 mu_1 1 - gamma_p mu_0 1 - 3 gamma_p mu_0 3\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"The system of ODEs is now much simpler and we can see that there are two higher-order moments we need to truncate: mu_13 and μ_14. We consider normal, derivative matching, conditional gaussian and conditional derivative matching closures to see how well they compare. First we apply different closures and print out the corresponding higher-order moment expressions in order to check that our results are consistent with those published in [1].","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"For normal closure:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"normal_eqs = moment_closure(eqs, \"normal\", binary_vars)\r\nlatexify(normal_eqs, :closure)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nmu_13 = mu_03 mu_10 + 3 mu_01 mu_12 + 3 mu_02 mu_11 + 6 mu_10 mu_01^3 - 6 mu_11 mu_01^2 - 6 mu_01 mu_02 mu_10 \r\nmu_14 = mu_04 mu_10 + 4 mu_01 mu_13 + 6 mu_02 mu_12 + 4 mu_03 mu_11 + 24 mu_11 mu_01^3 + 36 mu_02 mu_10 mu_01^2 - 24 mu_10 mu_01^4 - 6 mu_10 mu_02^2 - 12 mu_12 mu_01^2 - 24 mu_01 mu_02 mu_11 - 8 mu_01 mu_03 mu_10\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"Note that our expression for μ_14 is different from the equivalent langle gp^4 rangle in Eq. (8) of [1] as it seems that the authors introduced a mistake considering the central moment and cumulant relationship in Eq. (7).","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"Using derivative matching, we obtain:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"dm_eqs = moment_closure(eqs, \"derivative matching\", binary_vars)\r\nlatexify(dm_eqs, :closure)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nmu_1 3 = fracmu_0 1^3 mu_1 2^3 mu_0 3 mu_1 0mu_0 2^3 mu_1 1^3 \r\nmu_1 4 = fracmu_0 2^6 mu_1 1^4 mu_1 3^4 mu_0 4mu_0 1^4 mu_0 3^4 mu_1 2^6 mu_1 0\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"For conditional gaussian closure:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"cond_gaussian_eqs = moment_closure(eqs, \"conditional gaussian\", binary_vars)\r\nlatexify(cond_gaussian_eqs, :closure)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nmu_1 3 = frac-2 mu_1 1^3mu_1 0^2 + frac3 mu_1 1 mu_1 2mu_1 0 \r\nmu_1 4 = frac3 mu_1 2^2mu_1 0 + frac6 mu_1 1^4mu_1 0^3 + frac-12 mu_1 1^2 mu_1 2mu_1 0^2 + frac4 mu_1 1 mu_1 3mu_1 0\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"And, finally, for conditional derivative matching:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"cond_dm_eqs = moment_closure(eqs, \"conditional derivative matching\", binary_vars)\r\nlatexify(cond_dm_eqs, :closure)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"beginalign*\r\nmu_1 3 = fracmu_0 1^3 mu_1 2^3 mu_0 3 mu_1 0mu_0 2^3 mu_1 1^3 \r\nmu_1 4 = fracmu_0 2^6 mu_1 1^4 mu_1 3^4 mu_0 4mu_0 1^4 mu_0 3^4 mu_1 2^6 mu_1 0\r\nendalign*","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"All these results are consistent with [1], reassuring that the model and closures are implemented as intended. Finally, we can proceed to solve the resulting ODEs and compare the resulting means and standard deviations to the SSA. Following Soltani et al. [1], we define all model parameters as:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"mean_p = 200\r\nmean_b = 70\r\nγ_p_val = 1\r\nk_off_val = 0.001\r\nk_on_val = 0.05\r\n\r\nk_p_val = mean_p * γ_p_val * (k_off_val * mean_p^2 + k_on_val) / (k_on_val * mean_b)\r\n\r\nsymmap = [:k_on => k_on_val,\r\n :k_off => k_off_val,\r\n :k_p => k_p_val,\r\n :γ_p => γ_p_val,\r\n :b => mean_b]\r\npmap = symmap_to_varmap(rn, symmap)\r\n\r\n# initial gene state and protein number, order [g, p]\r\nu₀ = [1, 1]\r\n\r\n# time interval to solve on\r\ntspan = (0., 6.0)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"The reaction network with geometric bursts can be simulated using the SSA as usual:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"# convert the reaction network into a system of jump processes\r\njsys = convert(JumpSystem, rn; combinatoric_ratelaws=false)\r\n\r\n# create a discrete problem setting the simulation parameters\r\ndprob = DiscreteProblem(u₀, tspan, pmap)\r\n\r\n# create a JumpProblem compatible with ReactionSystemMod\r\njprob = JumpProblem(rn, dprob, Direct(), save_positions=(false, false))\r\n\r\n# simulate 2×10⁴ SSA trajectories\r\nensembleprob = EnsembleProblem(jprob)\r\n@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.1, trajectories=20000)\r\n# compute the means and variances\r\nmeans_ssa, vars_ssa = timeseries_steps_meanvar(sol_SSA)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"4.858357 seconds (14.30 M allocations: 901.032 MiB, 6.25% gc time, 63.45% compilation time)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We continue to solve the moment equations for each closure:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plt_m = plot() # plot mean protein number\r\nplt_std = plot() # plot ssd of protein number\r\n\r\n# construct the initial molecule number mapping\r\nu₀map = deterministic_IC(u₀, dm_eqs)\r\n\r\n# solve moment ODEs for each closure and plot the results\r\nfor closure in [\"normal\", \"derivative matching\",\r\n \"conditional gaussian\", \"conditional derivative matching\"]\r\n\r\n # it is very quick so we just apply all closures again\r\n closed_eqs = moment_closure(eqs, closure, binary_vars)\r\n\r\n # solve the system of moment ODEs\r\n oprob = ODEProblem(closed_eqs, u₀map, tspan, pmap)\r\n sol = solve(oprob, AutoTsit5(Rosenbrock23()), saveat=0.01)\r\n\r\n # μ₀₁ is 2nd and μ₀₂ is 4th element in sol\r\n # can check the order with `closed_eqs.odes.states`\r\n plt_m = plot!(plt_m, sol, idxs=[2], label=closure)\r\n plt_std = plot!(plt_std, sol.t, sqrt.(sol[4, :] .- sol[2, :].^2), label=closure)\r\n\r\nend\r\n\r\nplt_m = plot!(plt_m, xlabel=\"Time [hr]\", ylabel=\"Protein mean level\")\r\nplt_m = plot!(plt_m, means_ssa.t, means_ssa[2,:], label=\"SSA\", linestyle=:dash, color=\"gray\")\r\nplt_std = plot!(plt_std, xlabel=\"Time [hr]\", ylabel=\"Protein standard deviation\")\r\nplt_std = plot!(plt_std, vars_ssa.t, sqrt.(vars_ssa[2,:]), label=\"SSA\", linestyle=:dash, color=\"gray\")","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plot(plt_m, lw=2)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"(Image: Gene 1 means)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plot(plt_std, lw=2)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"(Image: Gene 2 stds)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We observe that the two conditional closures give the most accurate results. Therefore, closing moments conditioned on the gene being active is an effective approximation in this case. Note that our results for normal closure are different from [1], most likely due to the incorrect fifth-order moment closure function given in the paper.","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We now turn to a more complex network involving two genes (two Bernoulli variables), a repressor-activator circuit: an activator protein Y expressed by one gene can turn ON another gene which produces a repressor protein X, that subsequently can bind to the gene promoter region of the activator-producing gene turning it OFF (see [1] for more details). As before, X and Y proteins are expressed in geometrically-distributed bursts. The system can be specified as follows:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"@parameters b_x b_y\r\n# could redefine in terms of success probabilities p_x and p_y as above\r\nm = rand(Distributions.Geometric(1/(1+b_x)))\r\nl = rand(Distributions.Geometric(1/(1+b_y)))\r\n\r\n# g_x - gene state of X protein producing gene\r\n# g_y - gene state of Y protein producing gene\r\n# x, y - proteins X and Y\r\nrn = @reaction_network begin\r\n @parameters kx_on kx_off ky_on ky_off k_x γ_x k_y γ_y\r\n kx_on*(1-g_x)*y, 0 --> g_x # 0 -> g_x\r\n kx_off, g_x --> 0 # g_x -> 0\r\n ky_on*(1-g_y), 0 --> g_y # 0 -> g_y\r\n ky_off*x, g_y --> 0 # g_y -> 0\r\n k_x*g_x, 0 --> $m*x # 0 -> mx, m ~ Geometric(mean=b_x)\r\n γ_x, x --> 0 # x -> 0\r\n k_y*g_y, 0 --> $l*y # 0 -> ly, l ~ Geometric(mean_b_y)\r\n γ_y, y --> 0 # y -> 0\r\nend\r\n\r\n# both g_x and g_y are Bernoulli random variables\r\nbinary_vars = [1, 2]\r\n\r\n# Parameter initialisation\r\n\r\nmean_x = 100\r\nmean_y = 100\r\nmean_b_x = 5\r\nmean_b_y = 5\r\nγ_x_val = 1\r\nγ_y_val = 1\r\nkx_off_val = 4\r\nky_on_val = 0.3\r\nkx_on_val = 0.05\r\nky_off_val = 0.05\r\n\r\nk_x_val = mean_x * γ_x_val * (kx_off_val * mean_y^2 + kx_on_val) / (kx_on_val * mean_b_x)\r\nk_y_val = mean_y * γ_y_val * (ky_off_val * mean_x^2 + ky_on_val) / (ky_on_val * mean_b_y)\r\n\r\n# unclear if Soltani et al. (2015) actually used this parameter set as X numbers\r\n# jump to millions making SSA extremely slow...\r\n\r\n# introduce additional rescaling (otherwise rate coefficients are too high)\r\nk_x_val *= 0.00003\r\nk_y_val *= 0.01\r\n\r\n# parameter mapping\r\nsymmap = [:kx_on => kx_on_val,\r\n :kx_off => kx_off_val,\r\n :ky_on => ky_on_val,\r\n :ky_off => ky_off_val,\r\n :k_x => k_x_val,\r\n :k_y => k_y_val,\r\n :γ_x => γ_x_val,\r\n :γ_y => γ_y_val,\r\n :b_x => mean_b_x,\r\n :b_y => mean_b_y]\r\npmap = symmap_to_varmap(rn, symmap)\r\n\r\n# initial gene state and protein number, order [g_x, g_y, x, y]\r\nu₀ = [1, 1, 1, 1]\r\n\r\n# time interval to solve on\r\ntspan = (0., 12.0)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"Note that here we use a different parameter set from the one considered in [1] as otherwise the number of X proteins goes into millions, making the comparison between closure results difficult due to SSA being extremely slow—it appears that either we have included a mistake setting the parameter values or there is a typo in the paper.","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We can run SSA as follows:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"jsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)\r\ndprob = DiscreteProblem(jsys, u₀, tspan, pmap)\r\njprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false, false))\r\n\r\nensembleprob = EnsembleProblem(jprob)\r\n@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.1, trajectories=10000)\r\nmeans_ssa, vars_ssa = timeseries_steps_meanvar(sol_SSA)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"13.712331 seconds (43.27 M allocations: 1.697 GiB, 6.56% gc time, 10.01% compilation time)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"We use fourth order moment expansion and apply different closure methods as done for the negative feedback loop, considering the number of activator proteins Y and its standard deviation:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plt_m = plot() # plot mean activator protein number\r\nplt_std = plot() # plot ssd of activator protein number\r\n\r\nfor closure in [\"derivative matching\", \"conditional derivative matching\"]\r\n\r\n closed_eqs = moment_closure(eqs, closure, binary_vars)\r\n\r\n u₀map = deterministic_IC(u₀, closed_eqs)\r\n oprob = ODEProblem(closed_eqs, u₀map, tspan, pmap)\r\n sol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\n # μ₀₀₀₁ is the 4th and μ₀₀₀₂ is the 12th element in sol (can check with closed_eqs.odes.states)\r\n plt_m = plot!(plt_m, sol, idxs=[4], label=closure)\r\n plt_std = plot!(plt_std, sol.t, sqrt.(sol[12, :] .- sol[4, :].^2), label=closure)\r\nend\r\n\r\nplt_m = plot!(plt_m, xlabel=\"Time [hr]\", ylabel=\"Activator mean level\")\r\nplt_m = plot!(plt_m, means_ssa.t, means_ssa[4,:], label=\"SSA\", linestyle=:dash, color=\"gray\")\r\nplt_std = plot!(plt_std, xlabel=\"Time [hr]\", ylabel=\"Activator standard deviation\")\r\nplt_std = plot!(plt_std, vars_ssa.t, sqrt.(vars_ssa[4,:]), label=\"SSA\", linestyle=:dash, color=\"gray\")","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"Note that here we apply only derivative matching and conditional derivative matching, as both normal and conditional gaussian closures failed to be evaluated for the entire time course (irrespective of the solver used). Finally, we can plot the results:","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plot(plt_m, lw=2)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"(Image: Gene 2 means)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"plot(plt_std, lw=2, xlims=(0., 12.))","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"(Image: Gene 2 stds)","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"The trajectories obtained with conditional derivative matching quite accurately match the true dynamics, showing that it is indeed a powerful closure method for reaction networks involving species with very low (binary) copy numbers.","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/#References","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"References","text":"","category":"section"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"[1]: M. Soltani, C. A. Vargas-Garcia, and A. Singh, \"Conditional Moment Closure Schemes for Studying Stochastic Dynamics of Genetic Circuits\", IEEE Transactions on Biomedical Circuits and Systems 9, 518–526 (2015). https://doi.org/10.1109/tbcas.2015.2453158","category":"page"},{"location":"tutorials/geometric_reactions+conditional_closures/","page":"Geometrically Distributed Reaction Products and Conditional Closures","title":"Geometrically Distributed Reaction Products and Conditional Closures","text":"[2]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Approximation and inference methods for stochastic biochemical kinetics - a tutorial review\", Journal of Physics A: Mathematical and Theoretical 50, 093001 (2017). https://doi.org/10.1088/1751-8121/aa54d9","category":"page"},{"location":"tutorials/common_issues/#common_issues","page":"Common Issues","title":"Common Issues","text":"","category":"section"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Moment closure approximations are based on ad hoc assumptions and no rigorous and general predictions can be made on whether the results will be accurate or even physically meaningful [1-2]. Moreover, the truncated moment equations are prone to numerical instabilities and it may not be possible to solve them for the entire time course [3]. In this tutorial, we walk through a number of such issues encountered in the analysis of the Brusselator model introduced in the previous tutorial.","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"We first redefine the system and its parameters for completeness:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"using MomentClosure, Catalyst, OrdinaryDiffEq, Plots\r\n\r\nrn = @reaction_network begin\r\n @parameters c₁ c₂ c₃ c₄ Ω\r\n (c₁/Ω^2), 2X + Y → 3X\r\n (c₂), X → Y\r\n (c₃*Ω, c₄), 0 ↔ X\r\nend\r\n\r\np = [0.9, 2, 1, 1, 100]\r\nu₀ = [1, 1]\r\ntspan = (0., 100.)\r\n\r\nraw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"As we have seen earlier, second-order moment expansion using normal closure approximates the true system dynamics sufficiently accurately but it's interesting to see how other closures compare. Let's try applying zero closure:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"closed_raw_eqs = moment_closure(raw_eqs, \"zero\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 1)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"The trajectory of μ₀₁ becomes negative and so zero closure fails to provide physically meaningful results for this parameter set. Note that is important to correctly choose the ODE solver depending on the stiffness of the system and the accuracy required. We tried a number of recommended DifferentialEquations solvers here but none seemed to improve the results.","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Let's apply log-normal closure next:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"closed_raw_eqs = moment_closure(raw_eqs, \"log-normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2, legend=:bottomright)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"We observe sustained oscillatory behaviour instead of the expected damped oscillations. This result is unphysical: single SSA trajectories (that may display sustained oscillations) get dephased over time and hence the ensemble average should always show damped or overdamped oscillations [1].","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Normal closure is also quite fragile. This can be seen by simply including the combinatorial scaling of the mass-action propensity functions with combinatoric_ratelaw=true which leads to unphysical sustained oscillatory trajectories:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true)\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 3)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Nevertheless, this can be improved upon by increasing the order of moment expansion:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true)\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2, legend=:bottomright)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 4)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Some dampening in the system is now visible. Increasing the expansion order to 4 finally leads to physically sensible results:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"raw_eqs = generate_raw_moment_eqs(rn, 4, combinatoric_ratelaws=true)\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 5)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"For dessert, we consider unphysical divergent trajectories—a frequent problem with moment equations [3]. A good example is the second-order moment expansion including the combinatorial scaling of propensities with log-normal closure applied:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"raw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=true)\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"log-normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Rodas4P(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 6)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"In contrast to normal closure, increasing the expansion order makes the problem worse:","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"raw_eqs = generate_raw_moment_eqs(rn, 3, combinatoric_ratelaws=true)\r\nclosed_raw_eqs = moment_closure(raw_eqs, \"log-normal\")\r\n\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)\r\noprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)\r\nsol = solve(oprob, Rodas4P(), saveat=0.1)\r\n\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"┌ Warning: Interrupted. Larger maxiters is needed.\r\n└ @ SciMLBase C:\\Users\\asukys\\.julia\\packages\\SciMLBase\\Afx1r\\src\\integrator_interface.jl:331","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"(Image: Brusselator issue 7)","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"Note that the solver throws a warning being unable to evaluate the trajectories for the entire time course (other solvers perform similarly in this case). This usually implies that the moment ODE system is too stiff and cannot be solved: the time derivatives grow unboundedly and the solver timestep is being constantly reduced, requiring an ever-increasing number of the solver iterations (hence the maxiters warning).","category":"page"},{"location":"tutorials/common_issues/#References","page":"Common Issues","title":"References","text":"","category":"section"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"[1]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Comparison of different moment-closure approximations for stochastic chemical kinetics\", The Journal of Chemical Physics 143, 185101 (2015). https://doi.org/10.1063/1.4934990","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"[2]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Validity conditions for moment closure approximations in stochastic chemical kinetics\", The Journal of Chemical Physics 141, 084103 (2014). https://doi.org/10.1063/1.4892838","category":"page"},{"location":"tutorials/common_issues/","page":"Common Issues","title":"Common Issues","text":"[3]: E. Lakatos, A. Ale, P. D. W. Kirk, and M. P. H. Stumpf, \"Multivariate moment closure techniques for stochastic kinetic models\", The Journal of Chemical Physics 143, 094107 (2015). https://doi.org/10.1063/1.4929837","category":"page"},{"location":"tutorials/parameter_estimation_SDE/#parameter_estimation_SDE","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"","category":"section"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"In this tutorial, we look at the problem of model parameter identification of a diffusion process given timeseries data of its moments. Namely, we use moment closure approximations (MAs) to make the optimisation process more efficient by reducing the model evaluation time. All credit for this tutorial goes to Flemming Holtorf!","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"We consider a noisy variation of the Lotka-Volterra model, describing the interaction between a predator and prey species:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"beginalign*\n beginbmatrix dx dy endbmatrix = \n beginbmatrix gamma_1 x(t) - gamma_2 x(t) y(t) gamma_4 x(t) y(t) - gamma_3 y(t) - frac12 y(t) endbmatrix dt \n + beginbmatrix gamma_5 x(t) 0 endbmatrix dW_t\nendalign*","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"We can define this system of SDEs using ModelingToolkit as follows:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using ModelingToolkit\n\n@variables t, x(t), y(t)\n@parameters γ1, γ2, γ3, γ4, γ5\nγ = [γ1, γ2, γ3, γ4, γ5] \ndrift_eqs = [Differential(t)(x) ~ γ[1] * x - γ[2] * x * y ;\n Differential(t)(y) ~ γ[4] * x * y - γ[3] * y - y*0.5]\ndiff_eqs = [γ[5]*x; 0]\nLV = SDESystem(drift_eqs, diff_eqs, t, [x,y], γ, name = :LV)","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"Next, we generate some data used for the parameter estimation. Namely, we collect timeseries data of means and variances of both species:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using DifferentialEquations, DifferentialEquations.EnsembleAnalysis\n\nN_samples = 1000\nTf = 10\nt_data = 0:0.2:Tf\np_true = [γ[1] => 1, γ[2] => 2, γ[3] => 1, γ[4] => 2, γ[5] => 0.1]\nu0 = [1.0, 0.25]\nLV_data = solve(EnsembleProblem(SDEProblem(LV, u0, (0.0, Tf), p_true)), saveat = t_data, trajectories = N_samples)\nmeans, vars = timeseries_steps_meanvar(LV_data)","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"Now we are ready to test if we can estimate the model parameters gamma_1 dotsc gamma_5 solely from the data collected above. We first approach this parameter identification problem using the asymptotically exact approach of estimating the means and variances of the process with ensemble averages. Accordingly, we construct the following loss function:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using LinearAlgebra\n\nLV_sde = SDEProblem(LV, u0, (0.0, Tf), zeros(5))\nfunction obj(p)\n prob = remake(LV_sde, p = p)\n sol = solve(EnsembleProblem(prob), saveat = t_data, trajectories = 1000)\n sol_mean, sol_vars = timeseries_steps_meanvar(sol)\n obj = sum(norm(sol_mean[i] - means[i])^2 for i in 1:length(t_data))\n obj += 1e4*sum(norm(sol_vars[i] - vars[i])^2 for i in 1:length(t_data))\n return obj\nend","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"We can use this loss function with any suitable optimisation routine to identify a reasonable choice of the model parameters. For example, we can use a very simple derivative-free optimizer (Nelder-Mead method) implemented in the Optim.jl package. Since a single evaluation of the objective function requires sampling and hence is relatively expensive, we impose the constraint that the optimiser cannot run more than 2 minutes:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using Optim\n\npmap = Dict(i => findfirst(isequal(γ[i]), parameters(LV)) for i in 1:5)\np_init = Dict(γ[1] => 1.3, γ[2] => 1.5, γ[3] => 1.4, γ[4] => 2.2, γ[5] => 0.1)\np = [p_init[γ[pmap[i]]] for i in 1:5]\nopt_sampling = Optim.optimize(obj, p, Optim.Options(time_limit = 120))","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"* Status: failure (exceeded time limit of 120.0)\n\n * Candidate solution\n Final objective value: 3.550929e+00\n\n * Found with\n Algorithm: Nelder-Mead\n\n * Convergence measures\n √(Σ(yᵢ-ȳ)²)/n ≰ 1.0e-08\n\n * Work counters\n Seconds run: 121 (vs limit 120)\n Iterations: 139\n f(x) calls: 337","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"We can visualise the moment statistics obtained using the estimated parameters and compare to the true data:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using Plots\n\nt_detail = collect(0:0.01:Tf) \nopt_sol = solve(EnsembleProblem(SDEProblem(LV, u0, (0.0, Tf), opt_sampling.minimizer)), saveat = t_detail, trajectories = 1000)\nopt_means, opt_vars = timeseries_steps_meanvar(opt_sol)\n\nmean_comp = scatter(t_data, [m[1] for m in means], color = :blue,\n xlabel = \"time\", ylabel = \"population size\", \n grid = false, title = \"means\", label = \"⟨x⟩ data\")\nscatter!(mean_comp, t_data, [m[2] for m in means], color = :red, label = \"⟨y⟩ data\")\nplot!(mean_comp, t_detail, [m[1] for m in opt_means], linewidth = 2, color = :blue, label = \"⟨x⟩ SDE model\")\nplot!(mean_comp, t_detail, [m[2] for m in opt_means], linewidth = 2, color = :red, label = \"⟨y⟩ SDE model\")\n\nvar_comp = scatter(t_data, [v[1] for v in vars], color = :blue, grid = false,\n xlabel = \"time\", title = \"variances\", label = \"σ²(x) data\", legend = :topleft)\nscatter!(var_comp, t_data, [v[2] for v in vars], color = :red, label = \"σ²(y) data\")\nplot!(var_comp, t_detail, [v[1] for v in opt_vars], color = :blue, label = \"σ²(x) SDE model\")\nplot!(var_comp, t_detail, [v[2] for v in opt_vars], color = :red, label = \"σ²(y) SDE model\")\n\nplot(mean_comp, var_comp, size = (1200.0, 400.0))","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"(Image: LV SDE fig1)","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"We observe that the identified parameters match the data reasonably well. However, one may suspect that the optimiser is converging to a local minimum as the fit is not perfect.","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"Now we approach the same model identification problem via MAs in the hope of cutting down model evaluation cost, allowing us to identify better parameters in the same (or less) time. To that end, we construct an approximation of the moment dynamics of the process assuming that the distribution of the system state is approximately log-normal over the simulation horizon (using log-normal MA). Then we can implement a simple loss function by comparing the moments predicted by the approximate model with those obtained from data:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"using MomentClosure\n\nLV_moments = moment_closure(generate_raw_moment_eqs(LV, 2), \"log-normal\")\nu0map = deterministic_IC(u0, LV_moments)\nclosed_moment_prob = ODEProblem(LV_moments, u0map, (0.0, Tf), zeros(5))\nfunction obj_MCA(p)\n prob = remake(closed_moment_prob; p=p)\n sol = solve(prob, Tsit5(), saveat = t_data)\n if sol.retcode == ReturnCode.Success\n obj = sum(norm(sol.u[i][1:2] - means[i])^2 for i in 1:length(t_data))\n obj += 1e4*sum((sol.u[i][3] - sol.u[i][1]^2 - vars[i][1])^2 for i in 1:length(t_data))\n obj += 1e4*sum((sol.u[i][5] - sol.u[i][2]^2 - vars[i][2])^2 for i in 1:length(t_data))\n else\n obj = 1e6\n end\n return obj\nend","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"As before, any suitable optimisation routine can now be used to identify parameter values that result in a match between data and model prediction:","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"pmap = Dict(i => parameters(LV_moments.odes)[i] for i in 1:5)\np_init = Dict(γ[1] => 1.3, γ[2] => 1.5, γ[3] => 1.4, γ[4] => 2.2, γ[5] => 0.1)\np = [p_init[pmap[i]] for i in 1:5]\nopt_MCA = Optim.optimize(obj_MCA, p, Optim.Options(time_limit = min(120, opt_sampling.time_run)))","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"* Status: success\n\n * Candidate solution\n Final objective value: 1.160452e+00\n\n * Found with\n Algorithm: Nelder-Mead\n\n * Convergence measures\n √(Σ(yᵢ-ȳ)²)/n ≤ 1.0e-08\n\n * Work counters\n Seconds run: 0 (vs limit 120)\n Iterations: 362\n f(x) calls: 751","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"By visualising the results we see that indeed the identified parameters now provide a better match between data and the model, even if the original SDE model is being evaluated (note that the moment equations using log-normal MA provide a reasonably accurate approximation to the ensemble averages):","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"p_opt = [pmap[i] => opt_MCA.minimizer[i] for i in 1:5]\nt_detail = collect(0:0.01:Tf) \nopt_sol = solve(EnsembleProblem(SDEProblem(LV, u0, (0.0, Tf), p_opt)), saveat = t_detail, trajectories = 1000)\nopt_means = [timestep_mean(opt_sol, i) for i in 1:length(t_detail)]\nopt_vars = [timestep_meanvar(opt_sol, i)[2] for i in 1:length(t_detail)]\n\nopt_sol_approx = solve(ODEProblem(LV_moments, u0map, (0.0, Tf), p_opt), saveat = t_detail)\n\nmean_comp = scatter(t_data, [m[1] for m in means], color = :blue,\n xlabel = \"time\", ylabel = \"population size\", \n grid = false, title = \"means\", label = \"⟨x⟩ data\")\nscatter!(mean_comp, t_data, [m[2] for m in means], color = :red, label = \"⟨y⟩ data\")\nplot!(mean_comp, t_detail, [m[1] for m in opt_means], linewidth = 2, color = :blue, label = \"⟨x⟩ SDE model\")\nplot!(mean_comp, t_detail, [m[2] for m in opt_means], linewidth = 2, color = :red, label = \"⟨y⟩ SDE model\")\nplot!(mean_comp, t_detail, [m[1] for m in opt_sol_approx.u], linewidth = 2, color = :black, linestyle = :dash, label = \"log-normal MA\")\nplot!(mean_comp, t_detail, [m[2] for m in opt_sol_approx.u], linewidth = 2, color = :black, linestyle = :dash, label = nothing)\n\nvar_comp = scatter(t_data, [v[1] for v in vars], color = :blue,\n xlabel = \"time\", title = \"variances\", grid = false, label = \"σ²(x) data\", legend = :topleft)\nscatter!(var_comp, t_data, [v[2] for v in vars], color = :red, label = \"σ²(y) data\")\nplot!(var_comp, t_detail, [v[1] for v in opt_vars], color = :blue, label = \"σ²(x) SDE model\")\nplot!(var_comp, t_detail, [v[2] for v in opt_vars], color = :red, label = \"σ²(y) SDE model\")\nplot!(var_comp, t_detail, [m[3] - m[1]^2 for m in opt_sol_approx.u], linewidth = 2, color = :black, linestyle = :dash, label = \"log-normal MA\")\nplot!(var_comp, t_detail, [m[5] - m[2]^2 for m in opt_sol_approx.u], linewidth = 2, color = :black, linestyle = :dash, label = nothing)\n\nplot(mean_comp, var_comp, size = (1200.0, 400.0))","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"(Image: LV SDE fig2)","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"Finally, note that the parameter estimation is much faster using moment equations (compare to 120 s using the basic approach)","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"opt_MCA.time_run","category":"page"},{"location":"tutorials/parameter_estimation_SDE/","page":"Parameter Estimation of Diffusion Processes","title":"Parameter Estimation of Diffusion Processes","text":"0.10686397552490234","category":"page"},{"location":"tutorials/using_momentclosure_SDE/#tutorial_SDE","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"","category":"section"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"In the previous tutorial, we have shown how to generate and close moment equations for chemical reaction networks. Similarly, MomentClosure can be applied to systems of stochastic differential equations (SDEs) and here we demonstrate this functionality by working through a couple of practical examples. Namely, we consider two specific models from a relevant paper by Ghusinga et al. (2017) [1], 1) the Van der Pol oscillator and 2) a swinging pendulum, and apply moment closure approximations (MA) to reproduce some of the published results. For the theory behind the SDE moment expansion see here.","category":"page"},{"location":"tutorials/using_momentclosure_SDE/#Model-1","page":"Moment Equations from SDEs","title":"Model 1","text":"","category":"section"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"The deterministic Van der Pol oscillator is described by [1]:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\n fracd^2 xdt^2 - epsilon(1-x^2)fracdxdt + omega_n^2 x = Acos(omega_g t)\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"where epsilon is the bifurcation parameter, omega_n the natural frequency, omega_g the force frequency and A the force amplitude. We formulate a stochastic equivalent of this system by assuming that the force is noisy and taking x_1 = x and x_2 = fracdxdt, so that","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\n dx_1 = x_2 dt \r\n dx_2 = left( epsilon(1-x_1^2) x_2 - omega_n^2 x_1 right)dt + Acos(omega_g t)dt + A dw_t\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"In Julia, systems of SDEs can be conveniently described using ModelingToolkit.jl that provides an SDESystem type. Defining such models is straightforward (as in this example) and in our case is done as follows:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"using ModelingToolkit\r\n\r\n@variables t, x₁(t), x₂(t)\r\n@parameters ϵ, ω_n, ω_g, A\r\n\r\ndrift_eqs = [Differential(t)(x₁) ~ x₂;\r\n Differential(t)(x₂) ~ ϵ*(1-x₁^2)*x₂ - ω_n^2*x₁ + A*cos(ω_g*t)]\r\ndiff_eqs = [0; A]\r\n\r\nvdp_model = SDESystem(drift_eqs, diff_eqs, t, [x₁, x₂], [ϵ, ω_n, ω_g, A], name = :VdP)\r\n\r\nps = [ϵ => 0.1, ω_n => 120*pi, ω_g => 120*pi, A => 2.5] # parameter values\r\nu0 = [0.1, 0.1] # initial conditions\r\ntspan = (0., 0.1) # simulation time limit","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"We can now easily extract the raw moment equations up to second order using MomentClosure","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"using MomentClosure, Latexify\r\n\r\nmoment_eqs = generate_raw_moment_eqs(vdp_model, 2) \r\nlatexify(moment_eqs)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\nfracdmu_10dt = mu_01 \r\nfracdmu_01dt = A cosleft( t omega_g right) + epsilon mu_01 - epsilon mu_21 - omega_n^2 mu_10 \r\nfracdmu_20dt = 2 mu_11 \r\nfracdmu_11dt = epsilon mu_11 + A mu_10 cosleft( t omega_g right) + mu_02 - epsilon mu_31 - omega_n^2 mu_20 \r\nfracdmu_02dt = A^2 + 2 epsilon mu_02 + 2 A mu_01 cosleft( t omega_g right) - 2 omega_n^2 mu_11 - 2 epsilon mu_22\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"The obtained moment equations are equivalent to the published ones (except for a single term in mu_11 equation as there is a typo in the paper), giving some confidence that everything works as intended!","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"Note that the moment equations are not closed as there are two fourth order moments involved, mu_31 and mu_22, and hence we need to resort to MAs. We proceed to apply the derivative matching approximation and solve the system of ODEs as follows:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"using DifferentialEquations\r\n\r\nclosed_eqs = moment_closure(moment_eqs, \"derivative matching\")\r\n\r\nu0map = deterministic_IC(u0, closed_eqs)\r\noprob = ODEProblem(closed_eqs, u0map, tspan, ps)\r\n\r\nsol_MA = solve(oprob, Tsit5(), saveat=0.0001)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"Finally, we can also use DifferentialEquations to solve the SDEs for many trajectories and compare the ensemble statistics to the derivative matching approximation (getting a great match):","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"using DifferentialEquations.EnsembleAnalysis, Plots\r\n\r\nprob_SDE = SDEProblem(vdp_model, u0, tspan, ps)\r\n@time sol_SDE = solve(EnsembleProblem(prob_SDE), SRIW1(), saveat=0.0001, trajectories=100)\r\nmeans_SDE = timeseries_steps_mean(sol_SDE)\r\n\r\nplot(sol_MA.t, sol_MA[1, :], lw=2, label=\"MA\", ylabel=\"⟨x₁⟩\", xlabel=\"time\")\r\nplot!(sol_MA.t, means_SDE[1, :], lw=2, label=\"SDE\", linecolor=:red,\r\n linestyle=:dash, background_color_legend=nothing, legend=:topright, grid=false)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"(Image: Model 1 SDE)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/#Model-2","page":"Moment Equations from SDEs","title":"Model 2","text":"","category":"section"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"We next consider a simple pendulum, deterministically described as [1]","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\n fracd^2thetadt^2 + frackmfracdthetadt + fracgl sintheta = 0\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"where theta is the angular displacement, m the mass of the pendulum, g the gravitational acceleration constant, l the pendulum length and k the friction constant. The stochasticity is modelled as a white noise term arising due to pendulum randomly interacting with air molecules (assumed to be inversely proportional to m). Letting x_1 = theta and x_2 = fracdthetadt we write down the SDE as:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\n dx_1 = x_2 dt \r\n dx_2 = left( -frackmx_2 - fracglsinx_1 right)dt + frac1mdw_t\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"Using ModelingToolkit again we can define the model as","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"@variables t, x₁(t), x₂(t)\r\n@parameters k, l, m, g\r\n\r\ndrift_eqs = [Differential(t)(x₁) ~ x₂;\r\n Differential(t)(x₂) ~ -k/m*x₂ - g/l*sin(x₁)]\r\ndiff_eqs = [0; 1/m]\r\n\r\npendulum_model = SDESystem(drift_eqs, diff_eqs, t, [x₁, x₂], [k, l, m, g], name = :pendulum)\r\nps = [k => 10, m => 10, l => 10, g => 10]\r\nu0 = [3, 3]\r\ntspan = (0., 15.)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"As the drift term is non-polynomial, we can no longer write down the raw moment equations in a straightforward manner. Nevertheless, we can expand the expressions up to a certain Taylor expansion order q and generate the corresponding time-evolution equations for the means and higher order central moments—this situation is identical to moment expansion for chemical reaction networks with non-polynomial propensities. Hence we can obtain the equations for the means and variances with q=3:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"moment_eqs = generate_central_moment_eqs(pendulum_model, 2, 3) \r\nlatexify(moment_eqs) # the output here is maybe not the most visually pleasing","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"beginalign*\r\nfracdmu_10dt = mu_01 \r\nfracdmu_01dt = fracleft( - k right) mu_01m + fracleft( - g right) sinleft( mu_10 right)l + fracfrac16 g M_30 cosleft( mu_10 right)l + fracfrac12 g M_20 sinleft( mu_10 right)l \r\nfracdM_20dt = 2 M_11 \r\nfracdM_11dt = fracleft( - k right) M_11m + fracleft( - g right) M_20 cosleft( mu_10 right)l + fracfrac12 g M_30 sinleft( mu_10 right)l + M_02 \r\nfracdM_02dt = fracfrac11m^2 + frac-2 k M_02m + fracg M_21 sinleft( mu_10 right)l + frac-2 g M_11 cosleft( mu_10 right)l\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"We can now approximate the moment equations using gamma closure, solve the ODEs numerically and compare the predictions to true SDE solution, finding good agreement:","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"closed_eqs = moment_closure(moment_eqs, \"gamma\")\r\n\r\nu0map = deterministic_IC(u0, closed_eqs)\r\noprob = ODEProblem(closed_eqs, u0map, tspan, ps)\r\nsol_MA = solve(oprob, Tsit5(), saveat=0.01)\r\n\r\nprob_SDE = SDEProblem(pendulum_model, u0, tspan, ps)\r\nsol_SDE = solve(EnsembleProblem(prob_SDE), SRIW1(), saveat=0.01, trajectories=100)\r\nmeans_SDE = timeseries_steps_mean(sol_SDE)\r\n\r\nplot(sol_MA.t, sin.(sol_MA[1, :]), lw=2, label=\"MA\", ylabel=\"sin(⟨x₁⟩)\", xlabel=\"time\")\r\nplot!(sol_MA.t, sin.(means_SDE[1, :]), lw=2, label=\"SDE\", linecolor=:red,\r\n linestyle=:dash, background_color_legend=nothing, legend=:topright, grid=false)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"(Image: Model 2 SDE)","category":"page"},{"location":"tutorials/using_momentclosure_SDE/#References","page":"Moment Equations from SDEs","title":"References","text":"","category":"section"},{"location":"tutorials/using_momentclosure_SDE/","page":"Moment Equations from SDEs","title":"Moment Equations from SDEs","text":"[1]: K. R. Ghusinga, M. Soltani, A. Lamperski, S. V. Dhople, and A. Singh, \"Approximate moment dynamics for polynomial and trigonometric stochastic systems\", IEEE 56th Annual Conference on Decision and Control (2017). https://doi.org/10.1109/CDC.2017.8263922","category":"page"},{"location":"tutorials/using_momentclosure/#main_tutorial","page":"Using MomentClosure","title":"Using MomentClosure","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"This tutorial is an introduction to using MomentClosure to define chemical reaction network models, generate the corresponding moment equations, apply moment closure approximations and finally solve the resulting system of ODEs. To demonstrate this functionality, we will consider a specific case of an oscillatory chemical system known as the Brusselator, characterised by the reactions","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\n2X + Y stackrelc_1rightarrow 3X \r\nX stackrelc_2rightarrow Y \r\n undersetc_4stackrelc_3rightleftharpoons X\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"We have chosen this particular model to start with as it has been studied with different moment closures before by Schnoerr et al. (2015) and so it is useful as a reference point. The plots of moment trajectories we obtain in this tutorial fully reproduce some of the figures published in the paper, hence (partially) proving the validity of our implementation.","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The terminology and the notation used throughout is consistent with the Theory section of the docs and we advise giving it a skim-through.","category":"page"},{"location":"tutorials/using_momentclosure/#Model-Initialisation","page":"Using MomentClosure","title":"Model Initialisation","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Catalyst.jl provides a comprehensive interface to modelling chemical reaction networks in Julia and can be used to construct models fully-compatible with MomentClosure. For more details on how to do so we recommend reading Catalyst's tutorial. This way, the Brusselator can be defined as:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"using Catalyst\r\nrn = @reaction_network begin\r\n # including system-size parameter Ω\r\n @parameters c₁ c₂ c₃ c₄ Ω\r\n (c₁/Ω^2), 2X + Y → 3X\r\n (c₂), X → Y\r\n (c₃*Ω, c₄), 0 ↔ X\r\nend","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The returned rn is an instance of ModelingToolkit.ReactionSystem. The net stoichiometry matrix and an array of the corresponding propensities, if needed, can be extracted directly from the model using Catalyst.netstoichmat and MomentClosure function propensities respectively.","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Note that MomentClosure also supports systems containing geometrically distributed reaction products that can be defined using Catalyst—see this tutorial for more details.","category":"page"},{"location":"tutorials/using_momentclosure/#Generating-Moment-Equations","page":"Using MomentClosure","title":"Generating Moment Equations","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"We can now obtain the moment equations. The system follows the law of mass action, i.e., all propensity functions are polynomials in molecule numbers X(t) and Y(t), and so we can generate either raw or central moment equations, as described in the Theory section on moment expansion.","category":"page"},{"location":"tutorials/using_momentclosure/#Raw-Moment-Equations","page":"Using MomentClosure","title":"Raw Moment Equations","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Let's start with the raw moment equations which we choose to generate up to second order (m=2):","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"using MomentClosure\r\nraw_eqs = generate_raw_moment_eqs(rn, 2, combinatoric_ratelaws=false)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Note that we have set combinatoric_ratelaw=false in order to ignore the factorial scaling factors which Catalyst adds to mass-action reactions. The function generate_raw_moment_eqs returns an instance of RawMomentEquations that contains a ModelingToolkit.ODESystem composed of all the moment equations (accessed by raw_eqs.odes).","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"We can use Latexify to look at the generated moment equations:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"using Latexify\r\nlatexify(raw_eqs)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\nfracdmu_1 0dt = c_3 Omega + c_1 Omega^-2 mu_2 1 - c_2 mu_1 0 - c_4 mu_1 0 - c_1 Omega^-2 mu_1 1 \r\nfracdmu_0 1dt = c_2 mu_1 0 + c_1 Omega^-2 mu_1 1 - c_1 Omega^-2 mu_2 1 \r\nfracdmu_2 0dt = c_2 mu_1 0 + c_3 Omega + c_4 mu_1 0 + 2 c_3 Omega mu_1 0 + 2 c_1 Omega^-2 mu_3 1 - 2 c_2 mu_2 0 - 2 c_4 mu_2 0 - c_1 Omega^-2 mu_1 1 - c_1 Omega^-2 mu_2 1 \r\nfracdmu_1 1dt = c_2 mu_2 0 + c_3 Omega mu_0 1 + c_1 Omega^-2 mu_1 1 + c_1 Omega^-2 mu_2 2 - c_2 mu_1 0 - c_2 mu_1 1 - c_4 mu_1 1 - c_1 Omega^-2 mu_1 2 - c_1 Omega^-2 mu_3 1 \r\nfracdmu_0 2dt = c_2 mu_1 0 + 2 c_2 mu_1 1 + c_1 Omega^-2 mu_2 1 + 2 c_1 Omega^-2 mu_1 2 - c_1 Omega^-2 mu_1 1 - 2 c_1 Omega^-2 mu_2 2\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The raw moments are defined as","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"mu_ij(t) = langle X(t)^i Y(t)^j rangle","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"where langle rangle denote the expectation value and we have explicitly included the time-dependence for completeness (made implicit in the formatted moment equations). Note that the ordering of species (X first and Y second) is consistent with the order these variables appear within the Catalyst.@reaction_network macro. The ordering can also be checked using Catalyst.speciesmap function:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"speciesmap(rn)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Dict{Term{Real},Int64} with 2 entries:\r\n X(t) => 1\r\n Y(t) => 2","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Coming back to the generated moment equations, we observe that they depend on higher-order moments. For example, the ODE for mu_02 depends on third order moments μ_12 and μ_21 and the fourth order moment mu_22. Consider the general case of raw moment equations: if a network involves reactions that are polynomials (in molecule numbers) of at most order k, then its m^textth order moment equations will depend on moments up to order m+k-1. Hence the relationship seen above is expected as the Brusselator involves a trimolecular reaction whose corresponding propensity function is a third order polynomial in X(t) and Y(t). The number denoting the highest order of moments encountered in the generated RawMomentEquations can also be accessed as raw_eqs.q_order (returning 4 in this case).","category":"page"},{"location":"tutorials/using_momentclosure/#Central-Moment-Equations","page":"Using MomentClosure","title":"Central Moment Equations","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The corresponding central moment equations can also be easily generated:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"central_eqs = generate_central_moment_eqs(rn, 2, combinatoric_ratelaws=false)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Note that in case of non-polynomial propensity functions the Taylor expansion order q must also be specified, see the P53 system example for more details. Luckily, the Brusselator contains only mass-action reactions and hence q is automatically determined by the highest order (polynomial) propensity. The function generate_central_moment_eqs returns an instance of CentralMomentEquations. As before, we can visualise the central moment equations:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"latexify(central_eqs)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\nfracdmu_1 0dt = c_3 Omega + c_1 Omega^-2 M_2 1 + c_1 Omega^-2 mu_1 0^2 mu_0 1 + c_1 Omega^-2 M_2 0 mu_0 1 + 2 c_1 Omega^-2 M_1 1 mu_1 0 - c_2 mu_1 0 - c_4 mu_1 0 - c_1 Omega^-2 M_1 1 - c_1 Omega^-2 mu_0 1 mu_1 0 \r\nfracdmu_0 1dt = c_2 mu_1 0 + c_1 Omega^-2 M_1 1 + c_1 Omega^-2 mu_0 1 mu_1 0 - c_1 Omega^-2 M_2 1 - c_1 Omega^-2 mu_1 0^2 mu_0 1 - c_1 Omega^-2 M_2 0 mu_0 1 - frac21 c_1 Omega^-2 M_1 1 mu_1 0 \r\nfracdM_2 0dt = c_3 Omega + c_2 mu_1 0 + c_4 mu_1 0 + 2 c_1 Omega^-2 M_3 1 + c_1 Omega^-2 mu_1 0^2 mu_0 1 + 2 c_1 Omega^-2 mu_1 0^2 M_1 1 + 2 c_1 Omega^-2 M_3 0 mu_0 1 + 4 c_1 Omega^-2 M_2 1 mu_1 0 + 4 c_1 Omega^-2 M_2 0 mu_0 1 mu_1 0 - 2 c_2 M_2 0 - 2 c_4 M_2 0 - c_1 Omega^-2 M_1 1 - c_1 Omega^-2 M_2 1 - c_1 Omega^-2 M_2 0 mu_0 1 - c_1 Omega^-2 mu_0 1 mu_1 0 \r\nfracdM_1 1dt = c_2 M_2 0 + c_1 Omega^-2 M_1 1 + c_1 Omega^-2 M_2 2 + c_1 Omega^-2 mu_1 0^2 M_0 2 + c_1 Omega^-2 M_2 1 mu_0 1 + c_1 Omega^-2 mu_0 1 mu_1 0 + 2 c_1 Omega^-2 M_1 2 mu_1 0 + 2 c_1 Omega^-2 M_1 1 mu_0 1 mu_1 0 - c_2 M_1 1 - c_4 M_1 1 - c_2 mu_1 0 - c_1 Omega^-2 M_1 2 - c_1 Omega^-2 M_3 1 - c_1 Omega^-2 mu_1 0^2 M_1 1 - c_1 Omega^-2 mu_1 0^2 mu_0 1 - c_1 Omega^-2 M_0 2 mu_1 0 - c_1 Omega^-2 M_1 1 mu_0 1 - c_1 Omega^-2 M_3 0 mu_0 1 - c_1 Omega^-2 M_1 1 mu_1 0 - 2 c_1 Omega^-2 M_2 1 mu_1 0 - 2 c_1 Omega^-2 M_2 0 mu_0 1 mu_1 0 \r\nfracdM_0 2dt = c_2 mu_1 0 + 2 c_2 M_1 1 + c_1 Omega^-2 M_2 1 + 2 c_1 Omega^-2 M_1 2 + c_1 Omega^-2 mu_1 0^2 mu_0 1 + c_1 Omega^-2 M_2 0 mu_0 1 + 2 c_1 Omega^-2 M_0 2 mu_1 0 + 2 c_1 Omega^-2 M_1 1 mu_0 1 + 2 c_1 Omega^-2 M_1 1 mu_1 0 - c_1 Omega^-2 M_1 1 - 2 c_1 Omega^-2 M_2 2 - 2 c_1 Omega^-2 mu_1 0^2 M_0 2 - c_1 Omega^-2 mu_0 1 mu_1 0 - 2 c_1 Omega^-2 M_2 1 mu_0 1 - 4 c_1 Omega^-2 M_1 2 mu_1 0 - 4 c_1 Omega^-2 M_1 1 mu_0 1 mu_1 0\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Unfortunately, central moment equations often take a visually painful form. Note that the first two ODEs, as before, indicate the means, and the central moments are denoted as","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"M_ij(t) = langle (X(t)-mu_10(t))^i (Y(t)-mu_01(t))^j rangle","category":"page"},{"location":"tutorials/using_momentclosure/#Applying-Moment-Closure","page":"Using MomentClosure","title":"Applying Moment Closure","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"As observed above, the moment equations of the Brusselator are coupled and depend on higher order moments—we have an infinite hierarchy of ODEs in our hands which cannot be solved directly and requires approximate treatment. One way of approaching the problem is to apply moment closure approximations (MAs), in which higher order moments are expressed as functions of lower order moments, thus effectively truncating the hierarchy and enabling a numerical solution. A variety of MAs have been proposed in literature and are implemented in MomentClosure.jl, see the Theory section on MAs for more details.","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Let's apply normal closure to the raw moment equations raw_eqs we have generated earlier using moment_closure function:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"closed_raw_eqs = moment_closure(raw_eqs, \"normal\")","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The function returns ClosedMomentEquations that consists of both the ModelingToolkit.ODESystem containing all closed moment equations as well as the specific closure functions for each higher order raw moments. We can use Latexify again to look at the closed ODEs:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"latexify(closed_raw_eqs)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\nfracdmu_10dt = c_3 Omega + c_1 mu_01 mu_20 Omega^-2 + 2 c_1 mu_10 mu_11 Omega^-2 - c_2 mu_10 - c_4 mu_10 - c_1 mu_11 Omega^-2 - 2 c_1 mu_01 Omega^-2 mu_10^2 \r\nfracdmu_01dt = c_2 mu_10 + c_1 mu_11 Omega^-2 + 2 c_1 mu_01 Omega^-2 mu_10^2 - c_1 mu_01 mu_20 Omega^-2 - 2 c_1 mu_10 mu_11 Omega^-2 \r\nfracdmu_20dt = c_2 mu_10 + c_3 Omega + c_4 mu_10 + 2 c_3 Omega mu_10 + 2 c_1 mu_01 Omega^-2 mu_10^2 + 6 c_1 mu_11 mu_20 Omega^-2 - 2 c_2 mu_20 - 2 c_4 mu_20 - c_1 mu_11 Omega^-2 - c_1 mu_01 mu_20 Omega^-2 - 4 c_1 mu_01 Omega^-2 mu_10^3 - 2 c_1 mu_10 mu_11 Omega^-2 \r\nfracdmu_11dt = c_2 mu_20 + c_1 mu_11 Omega^-2 + c_3 Omega mu_01 + c_1 mu_02 mu_20 Omega^-2 + 2 c_1 Omega^-2 mu_11^2 + 2 c_1 mu_01 Omega^-2 mu_10^3 + 2 c_1 mu_10 Omega^-2 mu_01^2 - c_2 mu_10 - c_2 mu_11 - c_4 mu_11 - 2 c_1 mu_01 mu_11 Omega^-2 - c_1 mu_02 mu_10 Omega^-2 - 3 c_1 mu_11 mu_20 Omega^-2 - 2 c_1 Omega^-2 mu_01^2 mu_10^2 \r\nfracdmu_02dt = c_2 mu_10 + 2 c_2 mu_11 + c_1 mu_01 mu_20 Omega^-2 + 4 c_1 mu_01 mu_11 Omega^-2 + 2 c_1 mu_02 mu_10 Omega^-2 + 2 c_1 mu_10 mu_11 Omega^-2 + 4 c_1 Omega^-2 mu_01^2 mu_10^2 - c_1 mu_11 Omega^-2 - 4 c_1 Omega^-2 mu_11^2 - 2 c_1 mu_01 Omega^-2 mu_10^2 - 2 c_1 mu_02 mu_20 Omega^-2 - 4 c_1 mu_10 Omega^-2 mu_01^2\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The closure functions can also be displayed by adding :closure argument:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"latexify(closed_raw_eqs, :closure)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\nmu_21 = mu_01 mu_20 + 2 mu_10 mu_11 - 2 mu_01 mu_10^2 \r\nmu_12 = mu_02 mu_10 + 2 mu_01 mu_11 - 2 mu_10 mu_01^2 \r\nmu_31 = mu_01 mu_30 + 6 mu_01 mu_10^3 + 3 mu_10 mu_21 + 3 mu_11 mu_20 - 6 mu_11 mu_10^2 - 6 mu_01 mu_10 mu_20 \r\nmu_22 = mu_02 mu_20 + 2 mu_11^2 + 2 mu_01 mu_21 + 2 mu_10 mu_12 + 6 mu_01^2 mu_10^2 - 2 mu_02 mu_10^2 - 2 mu_20 mu_01^2 - 8 mu_01 mu_10 mu_11\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Similarly, we can close central moment equations using","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"closed_central_eqs = moment_closure(central_eqs, \"normal\")","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"and print out the corresponding closure functions:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"latexify(closed_central_eqs, :closure)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"beginalign*\r\nM_30 = 0 \r\nM_21 = 0 \r\nM_12 = 0 \r\nM_31 = 3 M_11 M_20 \r\nM_22 = M_02 M_20 + 2 M_11^2\r\nendalign*","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Higher order central moments under normal closure take a rather simple form compared to their raw moment equivalents. That can be expected due to the relationship between central moments and cumulants, on which the closure is based.","category":"page"},{"location":"tutorials/using_momentclosure/#Solving-Moment-Equations","page":"Using MomentClosure","title":"Solving Moment Equations","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The closed moment equations can be solved numerically using DifferentialEquations.jl that provides a variety of highly-efficient solvers and analysis tools. In order to do so, first we need to specify the values of all system parameters, the initial condition and the timespan to solve over. Then the ModelingToolkit.ODESystem corresponding to the moment equations can be directly converted into an ODEProblem which can finally be solved. Let's go through the procedure step-by-step for the closed raw moment equations (closed_raw_eqs). Most of what is covered below is closely based on this Catalyst tutorial.","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"We start with the parameters. Note that they are ordered as they appear after the end statement in the @reaction_network macro and this ordering can also be checked using the paramsmap function:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"paramsmap(rn)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Dict{Sym{ModelingToolkit.Parameter{Real}},Int64} with 5 entries:\r\n c₃ => 3\r\n Ω => 5\r\n c₄ => 4\r\n c₂ => 2\r\n c₁ => 1","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"We can now create a vector of parameter values as:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"# parameters [c₁, c₂, c₃, c₄, Ω]\r\np = [0.9, 2, 1, 1, 100]","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Next, we can specify the initial condition. Usually when working with moment equations we consider deterministic initial conditions so that the molecule numbers at initial time take the specified integer values with probability one. We can define the initial molecule numbers as X(t=0) = X_0 and Y(t=0) = Y_0. Probability one implies that initially the means will be equal to the molecule numbers, i.e., μ_10(t=0) = X_0 and μ_01(t=0) = Y_0, whereas all higher order raw moments will be products of the corresponding powers of the means, e.g., μ_21 = X_0^2 Y_0. Note that all central moments would be set to zero in this case. To make life easier we use deterministic_IC function which, given the initial molecule numbers, automatically constructs the variable mapping under deterministic initial conditions:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"# initial molecule numbers [X, Y] (order as in speciesmap(rn))\r\nu₀ = [1, 1]\r\nu₀map = deterministic_IC(u₀, closed_raw_eqs)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The next ingredient, the time interval to solve on, can be specified simply as:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"tspan = (0., 100.)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Now we are able to create the corresponding ODEProblem:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"oprob = ODEProblem(closed_raw_eqs, u₀map, tspan, p)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Finally, we have everything we need to solve the raw moment equations which can be done using any ODE solver implemented within DifferentialEquations.jl. We use the default Tsit5() solver and then plot the obtained mean molecule numbers:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"# using only ODE solvers from DifferentialEquations (faster to load)\r\nusing OrdinaryDiffEq\r\nsol = solve(oprob, Tsit5(), saveat=0.1)\r\n\r\nusing Plots\r\nplot(sol, idxs=[1,2], lw=2)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"(Image: Brusselator means 1)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The obtained moment dynamics show damped oscillations which is the expected averaged behaviour of the Brusselator in a vast swathe of parameter space [1]. However, to establish more clearly how well a second order moment expansion with normal moment closure performs for this system and this specific set of parameters, we can compare the resulting moment trajectories to the true moment estimates obtained using Gillespie's Stochastic Simulation Algorithm (SSA).","category":"page"},{"location":"tutorials/using_momentclosure/#Stochastic-Simulation","page":"Using MomentClosure","title":"Stochastic Simulation","text":"","category":"section"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"To run the SSA for a given reaction network, we build a JumpProcesses (formerly known as DiffEqJump) JumpProblem using Gillespie's Direct method. Note that other SSA variants are also available, see the documentation. Moreover, in order to run many realisations of the jump process, we define a corresponding EnsembleProblem. All of this can be done as follows:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"using JumpProcesses\r\n\r\n# convert ReactionSystem into JumpSystem\r\njsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)\r\n\r\n# create a DiscreteProblem encoding that the molecule numbers are integer-valued\r\ndprob = DiscreteProblem(jsys, u₀, tspan, p) # same parameters as defined earlier\r\n\r\n# create a JumpProblem: specify Gillespie's Direct Method as the solver\r\n# and SET save_positions to (false, false) as otherwise time of each\r\n# reaction occurence would be saved (complicating moment estimates)\r\njprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false, false))\r\n\r\n# define an EnsembleProblem to simulate multiple trajectories\r\nensembleprob = EnsembleProblem(jprob)\r\n\r\n# simulate 10000 SSA trajectories\r\n@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.1, trajectories=10000)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"Now we use the DifferentialEquations ensemble statistics tools to calculate the SSA mean values and plot them:","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"using DiffEqBase.EnsembleAnalysis\r\n\r\nmeans_SSA = timeseries_steps_mean(sol_SSA)\r\nplot!(means_SSA, lw=2, labels=[\"SSA μ₁₀(t)\" \"SSA μ₀₁(t)\"], linestyle=:dash,\r\n linecolor=[1 2], background_color_legend=nothing, legend=:bottomright)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"(Image: Brusselator SSA 1)","category":"page"},{"location":"tutorials/using_momentclosure/","page":"Using MomentClosure","title":"Using MomentClosure","text":"The comparison to the SSA reveals that the second-order moment expansion using normal closure captures the correct qualitative behaviour of the Brusselator and provides reasonably accurate moment estimates given this particular parameter set. Note, however, that moment closure approximations can lead to unphysical results and suffer from numerical instabilities, please see the Common Issues tutorial for more details.","category":"page"},{"location":"#MomentClosure","page":"Home","title":"MomentClosure","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"MomentClosure.jl is a tool to automatically obtain time-evolution equations of moments up to an arbitrary order for virtually any chemical reaction network or system of stochastic differential equations (SDEs), implementing a wide array of moment closure approximations commonly used in stochastic biochemical kinetics [1]. MomentClosure is (attempted to be) fairly well-integrated within the broader Julia ecosystem utilising a number of familiar packages:","category":"page"},{"location":"","page":"Home","title":"Home","text":"MomentClosure can be immediately applied to reaction network models defined using Catalyst and SDE systems built with ModelingToolkit.\nMoment equations are generated as a ModelingToolkit ODESystem (with some extra help from Symbolics and SymbolicUtils).\nThe resulting ODESystem can be solved using any DifferentialEquations ODE solvers, enabling further study of the system using parameter estimation, sensitivity analysis and bifurcation analysis tools.","category":"page"},{"location":"#Features","page":"Home","title":"Features","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Chemical reaction networks containing any number of molecular species and reactions with any type of smooth propensity functions are supported. Models can be defined using Catalyst as ModelingToolkit.ReactionSystem.\nAdded support for reaction networks involving reaction products that are geometrically distributed random variables. An example of such network is an autoregulatory gene network with bursty protein production where the burst size follows a geometric distribution, see this example for more details.\nEquations describing the time evolution of means and central moments of the number of molecules of each species in the system can be generated up to arbitrary order [2, 3]. Note that non-polynomial propensity functions are Taylor expanded to a specified order. Raw moment equations can also be generated for mass-action systems (where all propensity functions are polynomials). \nSDE systems defined as ModelingToolkit.SDESystem are supported (big thanks to @FHoltorf). Similarly to reaction networks with non-polynomial propensities, non-polynomial drift and diffusion coefficients are Taylor expanded.\nCurrently implemented moment closure approximations include:\nzero (central-moment neglect) closure [4, 5]\nnormal closure [4]\npoisson closure [4]\nlog-normal closure [5]\ngamma closure [5]\nderivative matching [6]\nconditional gaussian closure [7]\nconditional derivative matching [7]\nlinear mapping approximation [8]\nMoment equations are constructed as a ModelingToolkit.ODESystem that can be solved using any DifferentialEquations ODE solver. Moreover, parameter estimation, sensitivity analysis and bifurcation analysis tools can be applied to further study the resulting system of equations.\nLatexify can be used to generate LaTeX expressions of the corresponding moment equations.","category":"page"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"MomentClosure can be installed through the Julia package manager:","category":"page"},{"location":"","page":"Home","title":"Home","text":"]add MomentClosure\nusing MomentClosure","category":"page"},{"location":"#Citation","page":"Home","title":"Citation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"If you use MomentClosure in your work, please cite our paper:","category":"page"},{"location":"","page":"Home","title":"Home","text":"@article{MomentClosure2021,\n author = {Sukys, Augustinas and Grima, Ramon},\n title = \"{MomentClosure.jl: automated moment closure approximations in Julia}\",\n journal = {Bioinformatics},\n volume = {38},\n number = {1},\n pages = {289-290},\n year = {2021},\n month = {06},\n issn = {1367-4803},\n doi = {10.1093/bioinformatics/btab469},\n url = {https://doi.org/10.1093/bioinformatics/btab469},\n eprint = {https://academic.oup.com/bioinformatics/article-pdf/38/1/289/41891091/btab469.pdf},\n}","category":"page"},{"location":"#References","page":"Home","title":"References","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"[1]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Approximation and inference methods for stochastic biochemical kinetics - a tutorial review\", Journal of Physics A: Mathematical and Theoretical 50, 093001 (2017). https://doi.org/10.1088/1751-8121/aa54d9","category":"page"},{"location":"","page":"Home","title":"Home","text":"[2]: A. Ale, P. Kirk, and M. P. H. Stumpf, \"A general moment expansion method for stochastic kinetic models\", The Journal of Chemical Physics 138, 174101 (2013). https://doi.org/10.1063/1.4802475","category":"page"},{"location":"","page":"Home","title":"Home","text":"[3]: C. H. Lee, \"A Moment Closure Method for Stochastic Chemical Reaction Networks with General Kinetics\", MATCH Communications in Mathematical and in Computer Chemistry 70, 785-800 (2013). https://match.pmf.kg.ac.rs/electronicversions/Match70/n3/match70n3785-800.pdf","category":"page"},{"location":"","page":"Home","title":"Home","text":"[4]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Comparison of different moment-closure approximations for stochastic chemical kinetics\", The Journal of Chemical Physics 143, 185101 (2015). https://doi.org/10.1063/1.4934990","category":"page"},{"location":"","page":"Home","title":"Home","text":"[5]: E. Lakatos, A. Ale, P. D. W. Kirk, and M. P. H. Stumpf, \"Multivariate moment closure techniques for stochastic kinetic models\", The Journal of Chemical Physics 143, 094107 (2015). https://doi.org/10.1063/1.4929837","category":"page"},{"location":"","page":"Home","title":"Home","text":"[6]: A. Singh and J. P. Hespanha, \"Lognormal Moment Closures for Biochemical Reactions\", in Proceedings of the 45th IEEE Conference on Decision and Control, ISSN: 0191-2216 (Dec. 2006), pp. 2063-2068. https://doi.org/10.1109/CDC.2006.376994","category":"page"},{"location":"","page":"Home","title":"Home","text":"[7]: M. Soltani, C. A. Vargas-Garcia, and A. Singh, \"Conditional Moment Closure Schemes for Studying Stochastic Dynamics of Genetic Circuits\", IEEE Transactions on Biomedical Circuits and Systems 9, 518-526 (2015). https://doi.org/10.1109/TBCAS.2015.2453158","category":"page"},{"location":"","page":"Home","title":"Home","text":"[8]: Z. Cao and R. Grima, \"Linear mapping approximation of gene regulatory networks with stochastic dynamics\", Nature Communications 9, 3305 (2018). https://doi.org/10.1038/s41467-018-05822-0","category":"page"},{"location":"api/momentclosure_api/#api","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"CurrentModule = MomentClosure","category":"page"},{"location":"api/momentclosure_api/#Model-definition","page":"MomentClosure.jl API","title":"Model definition","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"MomentClosure is fully compatible with reaction network models defined using Catalyst and stored as a ModelingToolkit.ReactionSystem. Note that previously we had implemented our own ReactionSystemMod, that allowed us to consider systems containing reactions which products are independent geometrically distributed random variables. However, this is now deprecated as Catalyst has added support for such parameteric stoichiometries offering a much more complete and efficient feature set.","category":"page"},{"location":"api/momentclosure_api/#api-basic-network-properties","page":"MomentClosure.jl API","title":"Basic model properties","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"Moreover, we include a couple of tiny extensions to the Catalyst API:","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"propensities\r\nget_stoichiometry","category":"page"},{"location":"api/momentclosure_api/#MomentClosure.propensities","page":"MomentClosure.jl API","title":"MomentClosure.propensities","text":"propensities(rn::Union{ReactionSystem, ReactionSystemMod}; combinatoric_ratelaws=true)\n\nReturn a vector of propensity functions of all reactions in the given ReactionSystem.\n\nNotes:\n\ncombinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.get_stoichiometry","page":"MomentClosure.jl API","title":"MomentClosure.get_stoichiometry","text":"get_stoichiometry(rn::ReactionSystem, smap::AbstractDict)\n\nReturn the net stoichiometry matrix using the specified mapping of species to their indices.\n\nNotes:\n\nThis is a modification of Catalyst.netstoichmat that is used internally to deal with reactions involving symbolic stoichiometry coefficients.\nThe function also allows custom smap, so it is not limited to the default Catalyst.speciesmap ordering.\nTODO: remove once this Catalyst issue is resolved.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#Moment-Equations","page":"MomentClosure.jl API","title":"Moment Equations","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"generate_raw_moment_eqs\r\nRawMomentEquations\r\ngenerate_central_moment_eqs\r\nCentralMomentEquations\r\nbernoulli_moment_eqs","category":"page"},{"location":"api/momentclosure_api/#MomentClosure.generate_raw_moment_eqs","page":"MomentClosure.jl API","title":"MomentClosure.generate_raw_moment_eqs","text":"generate_raw_moment_eqs(rn::ReactionSystem, m_order::Int;\n langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn))\n\nGiven a ReactionSystem return the RawMomentEquations of the system generated up to m_order.\n\nNotes:\n\nThe 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.\nif 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.\ncombinatoric_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.\nsmap 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.\n\n\n\n\n\ngenerate_raw_moment_eqs(sys::SDESystem, m_order::Int)\n\nGiven an SDESystem, return the RawMomentEquations of the system generated up to m_order.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.RawMomentEquations","page":"MomentClosure.jl API","title":"MomentClosure.RawMomentEquations","text":"struct RawMomentEquations <: MomentClosure.MomentEquations\n\nRaw moment equations generated for the given system plus a number of helper parameters (used internally).\n\nFields\n\nodes: ModelingToolkit.ODESystem consisting of the time-evolution equations of raw moments.\nμ: Symbolic variables defining the raw moments.\nN: Number of species within the system.\nm_order: Order of moment equations.\nq_order: Expansion order.\niter_all: Vector of all index combinations up to q_order.\niter_m: Vector of all index combinations up to m_order.\niter_q: Vector of all index combinations of order greater than m_order up to q_order.\niter_1: Vector of index combinations of order 1.\n\n\n\n\n\n","category":"type"},{"location":"api/momentclosure_api/#MomentClosure.generate_central_moment_eqs","page":"MomentClosure.jl API","title":"MomentClosure.generate_central_moment_eqs","text":"generate_central_moment_eqs(rn::ReactionSystem, m_order::Int, q_order::Int=0;\n langevin::Bool=false, combinatoric_ratelaws::Bool=true, smap=speciesmap(rn))\n\nGiven a ReactionSystem return the CentralMomentEquations of the system generated up to m_order.\n\nNotes:\n\nif 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).\nif 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.\ncombinatoric_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.\nsmap 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.\n\n\n\n\n\ngenerate_central_moment_eqs(sys::SDESystem, m_order::Int, q_order::Int=0)\n\nGiven an SDESystem, return the CentralMomentEquations of the system generated up to m_order.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.CentralMomentEquations","page":"MomentClosure.jl API","title":"MomentClosure.CentralMomentEquations","text":"struct CentralMomentEquations <: MomentClosure.MomentEquations\n\nCentral moment equations generated for the given system plus a number of helper parameters (used internally).\n\nFields\n\nodes: ModelingToolkit.ODESystem consisting of the time-evolution equations of central moments.\nμ: Symbolic variables defining the means.\nM: Symbolic variables defining the central moments.\nN: Number of species within the system.\nm_order: Order of moment equations.\nq_order: Expansion order.\niter_all: Vector of all index combinations up to q_order.\niter_m: Vector of all index combinations up to m_order.\niter_q: Vector of all index combinations of order greater than m_order up to q_order.\niter_1: Vector of index combinations of order 1.\n\n\n\n\n\n","category":"type"},{"location":"api/momentclosure_api/#MomentClosure.bernoulli_moment_eqs","page":"MomentClosure.jl API","title":"MomentClosure.bernoulli_moment_eqs","text":"bernoulli_moment_eqs(sys::MomentEquations, binary_vars::Array{Int,1})\n\nGiven 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.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#Moment-Closure","page":"MomentClosure.jl API","title":"Moment Closure","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"moment_closure\r\nClosedMomentEquations\r\ndeterministic_IC","category":"page"},{"location":"api/momentclosure_api/#MomentClosure.moment_closure","page":"MomentClosure.jl API","title":"MomentClosure.moment_closure","text":"moment_closure(sys::MomentEquations, closure::String, binary_vars::Array{Int,1}=Int[])\n\nGiven MomentEquations, apply the specified moment closure approximation and return the ClosedMomentEquations.\n\nThe supported closure options are:\n\n\"zero\"\n\"normal\"\n\"log-normal\"\n\"poisson\"\n\"gamma\"\n\"derivative matching\"\n\"conditional gaussian\"\n\"conditional derivative matching\"\n\nNotes\n\nbinary_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.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.ClosedMomentEquations","page":"MomentClosure.jl API","title":"MomentClosure.ClosedMomentEquations","text":"struct ClosedMomentEquations <: MomentClosure.MomentEquations\n\nClosed moment equations and the corresponding closure functions.\n\nFields\n\nodes: ModelingToolkit.ODESystem consisting of the time-evolution equations of closed moments.\nclosure: Dictionary of moment closure functions for each higher order moment.\nopen_eqs: Original raw or central moment equations (before closure was applied).\n\n\n\n\n\n","category":"type"},{"location":"api/momentclosure_api/#MomentClosure.deterministic_IC","page":"MomentClosure.jl API","title":"MomentClosure.deterministic_IC","text":"deterministic_IC(u₀::Array{T, 1}, eqs::MomentEquations) where T<:Real\n\nGiven 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.\n\nNotes\n\nThe 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.\nThe 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).\nAs 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).\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#visualisation_api","page":"MomentClosure.jl API","title":"Displaying Equations and Closures","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"The generated moment equations can be converted into LaTeX expressions using Latexify as:","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"using Latexify\r\nlatexify(moment_eqs)","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"A ModelingToolkit.ODESystem (saved as moment_eqs.odes) can also be passed to latexify function directly but the output will be different as we apply additional formatting to the symbolic expressions.","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"Given ClosedMomentEquations, the closure functions can be visualised in the same way by adding a :closure argument:","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"latexify(moment_eqs, :closure)","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"Note that this will print out only those higher order moments which are found in the given moment equations. It is possible to print the closure functions of all higher order moments using print_all=true argument:","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"latexify(moment_eqs, :closure, print_all=true)","category":"page"},{"location":"api/momentclosure_api/#Linear-Mapping-Approximation","page":"MomentClosure.jl API","title":"Linear Mapping Approximation","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"linear_mapping_approximation","category":"page"},{"location":"api/momentclosure_api/#MomentClosure.linear_mapping_approximation","page":"MomentClosure.jl API","title":"MomentClosure.linear_mapping_approximation","text":" linear_mapping_approximation(rn_nonlinear::T, rn_linear::T, binary_vars::Array{Int,1}=Int[], m_order::Int=0;\n combinatoric_ratelaws = true) where T <: ReactionSystem\n\nGiven 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.\n\nNotes:\n\nrn_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.\nbinary_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.\nBy 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.\ncombinatoric_ratelaws=true uses binomials in calculating the propensity functions of a ReactionSystem, see the notes for ModelingToolkit.jumpratelaw.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#stochastic_simulation_utilities","page":"MomentClosure.jl API","title":"Stochastic Simulation Utilities","text":"","category":"section"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"We provide provides functions for higher-order moment extraction from SSA and FSP data:","category":"page"},{"location":"api/momentclosure_api/","page":"MomentClosure.jl API","title":"MomentClosure.jl API","text":"get_raw_moments\r\nget_central_moments\r\nget_cumulants\r\nget_moments_FSP","category":"page"},{"location":"api/momentclosure_api/#MomentClosure.get_raw_moments","page":"MomentClosure.jl API","title":"MomentClosure.get_raw_moments","text":"get_raw_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)\n\nGiven an EnsembleSolution of DifferentialEquations ensemble simulation, return a Dictionary of raw moments computed up to the specified order at each time step.\n\nNotes\n\nFor example, the dictionary key (2,0,1) maps to an array containing the values of the raw moment μ_201 at each time step.\nIt 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).\nMoments 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.\nOnly 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.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.get_central_moments","page":"MomentClosure.jl API","title":"MomentClosure.get_central_moments","text":"get_central_moments(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)\n\nGiven 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.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.get_cumulants","page":"MomentClosure.jl API","title":"MomentClosure.get_cumulants","text":"get_cumulants(sol::EnsembleSolution, order::Int; naive::Bool=true, b::Int=2)\n\nGiven 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.\n\n\n\n\n\n","category":"function"},{"location":"api/momentclosure_api/#MomentClosure.get_moments_FSP","page":"MomentClosure.jl API","title":"MomentClosure.get_moments_FSP","text":"get_moments_FSP(sol::ODESolution, order::Int, moment_type::String)\n\nGiven 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.\n\nNotes\n\nThe 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.\n\n\n\n\n\n","category":"function"}]
+[{"location":"theory/linear_mapping_approximation/#linear_mapping_approximation","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"","category":"section"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"The Linear Mapping Approximation (LMA) provides a novel way of approximating the solution of the CME and has been shown to be accurate for a variety of models of gene regulatory networks (GRNs) [1]. It is based on mapping a nonlinear GRN onto an equivalent linear GRN so that the exact solution of the linear system gives an approximate solution of the nonlinear system. The LMA is restricted in its applicability to chemical reaction networks where one of the substrates in each nonlinear reaction is a molecular species which copy number can be either zero or one (a Bernoulli/binary random variable). Note that a network can contain an arbitrary number of such species but more than one of them cannot be involved in any nonlinear reaction. Below we provide a short overview of the LMA and urge the reader to see the original paper for a more comprehensive description [1].","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"To explain the LMA procedure, we start by considering a simple example of a two-state GRN as done in a previous section on conditional closures, denoting the binary gene state by g and the protein number by p. We assume that there is a single nonlinear reversible reaction in the network given by","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"beginalign*\nG+P undersetσ_ustackrelσ_brightleftharpoons G^*\nendalign*","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"where P denotes the protein and the gene can be in either state G (g=1) or G^* (g=0). Our aim is to find an approximate time-dependent probability distribution of protein numbers p at time t.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"The steps of performing the LMA can then be described as follows:","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"1. Find the linear network by replacing any reversible nonlinear reaction (it must involve one binary species) in the nonlinear network by a reversible pseudo first-order reaction between the binary species' states.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"In our example, we replace the reaction above with G undersetσ_ustackrelbarσ_brightleftharpoons G^*, noting that the rate parameter is changed from σ_b to barσ_b. Moreover, this approach is just as applicable in case of cooperativity, e.g., G+nP rightleftharpoons G^* (where n is an integer indicating the cooperative order) would be similarly replaced with G rightleftharpoons G^*.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"2. Approximate the changed rate parameters of the linearised reactions by their expectation values.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"As noted in [1], the first-order reaction G stackrelbarσ_brightarrow G^* maps onto the second-order reaction G+P stackrelσ_brightarrow G^* if we choose barσ_b = σ_b left(p g=1 right), where p g=1 indicates the instantaneous protein number given the gene is in the state G. In LMA, we use the mean-field approximation taking the expectation value of the rate so that\nbeginalign*\nbarσ_b = σ_b left langle p g=1 right rangle = σ_b fracleftlangle pg rightrangleleftlangle g rightrangle\nendalign*\nThe same procedure can be extended to the general nonlinear reaction where n proteins bind cooperatively. The effective parameter is then given by\nbarσ_b = σ_b fracleftlangle prod_i=0^n-1 left( p-iright)g rightrangleleftlangle g rightrangle","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"3. Write down the moment equations for the linear network using the approximated stochastic rates.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"Note that the moment equations must be generated up the order given by the highest order nonlinear reaction in the network. If the only nonlinear reaction is the second-order reaction G+P stackrelσ_brightarrow G^*, we need to consider only moments up to the second order (as hinted by the functional form of barσ_b above). However, the moment hierarchy is otherwise closed, no additional moment closure approximations need to be performed, and therefore we can solve the moment equations in a straightforward manner.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"4. Solve the moment equations numerically up to time t and plug the resulting moment values into the equations for the effective parameters. Proceed to calculate the time-average of these parameters over the time-interval 0 t.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"In our example, plugging the solved-for moment values into the equation for barσ_b allows us to interpret the effective parameter as a time-dependent function barσ_b(t). However, as the time-dependent probability distribution solution of the CME for the nonlinear network with a general time-dependent barσ_b is most likely intractable, we transform barσ_b(t) into a time-independent constant by taking its time-average barσ_b^* = int_0^t barσ_b(t) dt t. This approach is justified in [1] by considering the Magnus expansion.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"5. Obtain the time-dependent probability distribution solution of the CME of the linear network assuming that the rate parameters of the linearised reactions are time-independent constants.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"Note that this step is the major limitation of the LMA as closed-form solutions are available only for a handful of systems (consult [1] for more details).","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"6. Finally, construct the approximate probability distribution of the nonlinear network at time t by replacing the respective rate parameters with their time-averaged equivalents obtained in the previous step.","category":"page"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"MomentClosure.jl provides automated generation of the closed moment equations using LMA given a nonlinear chemical reaction network and its linear equivalent. This encapsulates the first three steps of the LMA procedure outlined above which are general and can be seen as an original moment closure approximation. We apply the LMA on simple models of nonlinear GRNs and also discuss how the subsequent LMA steps can be performed in Julia on a case-by-case basis in this tutorial example.","category":"page"},{"location":"theory/linear_mapping_approximation/#References","page":"Linear Mapping Approximation","title":"References","text":"","category":"section"},{"location":"theory/linear_mapping_approximation/","page":"Linear Mapping Approximation","title":"Linear Mapping Approximation","text":"[1]: Z. Cao and R. Grima, \"Linear mapping approximation of gene regulatory networks with stochastic dynamics\", Nature Communications 9, 3305 (2018). https://doi.org/10.1038/s41467-018-05822-0","category":"page"},{"location":"tutorials/P53_system_example/#P53-System-Example","page":"P53 System Example","title":"P53 System Example","text":"","category":"section"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"In this example tutorial, we consider the oscillatory p53-Mdm2 system which has been studied by Lakatos et al. (2015) using central moment expansion up to 5th order with zero, normal, log-normal and gamma closures. We attempt to reproduce a number of results published in their paper.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"In contrast to the Brusselator we considered before, the p53-Mdm2 system contains non-polynomial propensities and hence is a good showcase for generating central moment equations using different values of q and comparing their accuracy. Note that this was not considered by Lakatos et al. who set q = m+1 throughout their study. We remind the reader that m denotes the moment expansion order and q is the order of Taylor expansion of the propensity functions: the moment equations are generated up to m^textth order and may include central moments up q^textth order so that all moments of orders m+1 dotsc q must be approximated using some moment closure method (see the Moment Expansion theory section for more details).","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"The p53-Mdm2 reaction network is described by the stoichiometric matrix","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"S = beginbmatrix\n 1 -1 -1 0 0 0\n 0 0 0 1 -1 0 \n 0 0 0 0 1 -1\n endbmatrix","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"and six reaction propensities","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"beginalign*\na_1 = k_1 \na_2 = k_2 x_ \na_3 = k_3 fracx yx + k_7 \na_4 = k_4 x \na_5 = k_5 y₀ \na_6 = k_6 y\nendalign*","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"where the variables are","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"number of p53 molecules x\nnumber of precursor of Mdm2 molecules y_0\nnumber of Mdm2 molecules y","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"with parameters","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"p53 production rate k_1\nMdm2-independent p53 degradation rate k_2\nsaturating p53 degradation rate k_3\np53-dependent Mdm2 production rate k_4\nMdm2 maturation rate k_5\nMdm2 degradation rate k_6\np53 threshold for degradation by Mdm2 k_7","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We begin by loading all the packages we will need","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"using Catalyst, MomentClosure, OrdinaryDiffEq, JumpProcesses,\n DiffEqBase.EnsembleAnalysis, Plots, Plots.PlotMeasures","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"and then build the model using Catalyst and set its parameters as follows:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"# → for mass-actions rate\n# ⇒ for non mass-actions rate\nrn = @reaction_network begin\n @parameters k₁ k₂ k₃ k₄ k₅ k₆ k₇\n (k₁), 0 → x\n (k₂), x → 0\n (k₃*x*y/(x+k₇)), x ⇒ 0\n (k₄*x), 0 ⇒ y₀\n (k₅), y₀ → y\n (k₆), y → 0\nend\n\n# parameters [k₁, k₂, k₃, k₄, k₅, k₆, k₇]\np = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01]\n\n# initial molecule numbers [x, y₀, y]\nu₀ = [70, 30, 60]","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Let's first simulate the reaction network using SSA in order to have a reference point of the real system dynamics. We choose a relatively long simulation time span in order to clearly see how the molecule numbers converge to their steady-state values and opt for 5 times 10^4 SSA realisations:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"# time interval to solve one on\ntspan = (0., 200.)\n\n# constructing the discrete jump problem using DifferentialEquations\njsys = convert(JumpSystem, rn, combinatoric_ratelaws=false)\ndprob = DiscreteProblem(jsys, u₀, tspan, p)\n\njprob = JumpProblem(jsys, dprob, Direct(), save_positions=(false, false))\nensembleprob = EnsembleProblem(jprob)\n\n# @time is only a rough estimate and not a proper benchmark\n@time sol_SSA = solve(ensembleprob, SSAStepper(), saveat=0.2, trajectories=50000)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"1332.443829 seconds (2.90 G allocations: 69.838 GiB, 42.81% gc time)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Single SSA trajectories show sustained oscillations in molecule numbers:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plot(sol_SSA[666], labels=[\"p53\" \"pre-Mdm2\" \"Mdm2\"], lw=2, tspan=(0., 100.),\n linecolor=[1 3 2], xlabel=\"Time [h]\", ylabel=\"Number of molecules\", size=(700, 400))","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 SSA one trajectory)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"The mean population dynamics show damped oscillations as the single trajectories get dephased over time and eventually reach steady-state values:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"means_SSA, vars_SSA = timeseries_steps_meanvar(sol_SSA)\nplot(means_SSA, labels=[\"p53\" \"pre-Mdm2\" \"Mdm2\"], lw=2, linecolor=[1 3 2],\n xlabel=\"Time [h]\", ylabel=\"Number of molecules\", size=(700, 400))","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 SSA mean trajectories)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Note that many more SSA runs are needed to completely denoise the mean trajectories (especially at later times).","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Finally, we can look at the marginal distributions noting a high level of asymmetry at the chosen time point t=25:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"data = componentwise_vectors_timepoint(sol_SSA, 25.0)\nh1 = histogram(data[1], normalize=true, xlabel=\"x\", ylabel=\"P(x)\")\nh2 = histogram(data[2], normalize=true, xlabel=\"y₀\", ylabel=\"P(y₀)\")\nh3 = histogram(data[3], normalize=true, xlabel=\"y\", ylabel=\"P(y)\")\nusing Plots.PlotMeasures\nplot(h1, h2, h3, legend=false, layout=(1,3), size = (1050, 250),\n left_margin = 5PlotMeasures.mm, bottom_margin = 7PlotMeasures.mm, guidefontsize=10)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 distribution)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Given the long runtime of the SSA above, it makes sense to turn to approximation methods such as moment closure which is usually much faster (suitable if only information about moments is needed). We start with the simplest second order, m =2, moment expansion and consider normal, log-normal and gamma closures. Our interest here is to see how the accuracy of each closure changes when the Taylor expansion order q is increased. Note that zero closure is equivalent to normal closure when q=3 (as third order central moments and cumulants are equivalent) and does not change when q is increased as all the higher order moments will be simply set to zero.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We can generate the corresponding plots of the mean p53 molecule number, x, as follows:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"closures = [\"normal\", \"log-normal\", \"gamma\"]\n\n# initialise separate plot for each closure\nplts = [plot() for i in 1:length(closures)]\n\nfor q in 3:6\n eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)\n for (closure, plt) in zip(closures, plts)\n closed_eqs = moment_closure(eqs, closure)\n\n u₀map = deterministic_IC(u₀, closed_eqs)\n oprob = ODEProblem(closed_eqs, u₀map, tspan, p)\n\n sol = solve(oprob, Tsit5(), saveat=0.1)\n plt = plot!(plt, sol, idxs=[1], lw=3, label = \"q = \"*string(q))\n end\nend\n\nfor plt in plts\n plt = plot!(plt, xlabel = \"Time [h]\", ylabel = \"Mean number of p53 molecules\")\n plt = plot!(plt, means_SSA.t, means_SSA[1,:], lw=2, linestyle=:dash, label = \"SSA\", color=\"gray\")\nend","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Normal closure:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plot(plts[1], size=(750, 450), leftmargin=2PlotMeasures.mm)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 normal means 2nd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Log-normal closure:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plot(plts[2], size=(750, 450), leftmargin=2PlotMeasures.mm)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 log-normal means 2nd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Zoomed-in log-normal closure (to better see the initial dampening in trajectories):","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plot(plts[2], xlims=(0., 50.), lw=3)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 log-normal means 2nd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Gamma closure:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plot(plts[3], size=(750, 450), leftmargin=2PlotMeasures.mm)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 gamma means 2nd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"When using normal closure, including central moments all the way up to 6th order (q=6) is needed in order to obtain qualitatively correct system behaviour. As the true distribution (as plotted above) is highly asymmetric, one may expect that performing normal closure, i.e., assuming that the system has a (symmetric) Gaussian distribution, is not the most optimal choice. For log-normal and gamma closures, quite accurate predictions are obtained with q=4 whereas q=6 tends to overestimate the dampening initially but better captures the steady-state behaviour. Note that for odd values of q log-normal and gamma closures introduced numerical instabilities and the moments could not be evaluated for the entire time course.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We can also plot the variance predictions:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"# rerunning the same calculations as they are reasonably fast\nplt = plot()\nfor q in [4,6]\n eqs = generate_central_moment_eqs(rn, 2, q, combinatoric_ratelaws=false)\n for closure in closures\n closed_eqs = moment_closure(eqs, closure)\n\n u₀map = deterministic_IC(u₀, closed_eqs)\n oprob = ODEProblem(closed_eqs, u₀map, tspan, p)\n sol = solve(oprob, Tsit5(), saveat=0.1)\n\n # index of M₂₀₀ can be checked with `u₀map` or `closed_eqs.odes.states`\n plt = plot!(plt, sol, idxs=[4], lw=3, label = closure*\" q = \"*string(q))\n end\nend\n\nplt = plot!(plt, xlabel = \"Time [h]\", ylabel = \"Variance of p53 molecule number\", legend=:topleft)\nplt = plot!(plt, means_SSA.t, vars_SSA[1,:], lw=2, linestyle=:dash, label = \"SSA\", color=\"gray\")\nplot(plt, size=(750, 450))","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 variances 2nd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Log-normal closure with q=6 clearly gives the best estimate of the variance. Therefore, including increasingly higher-order moment information in second-order moment expansion tends to leads to an overall more accurate statistical description of the P53-Mdm2 system.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We can extend the analysis further by considering higher order moment expansions (increasing m). Let's look at third-order expansion with q values 4 and 6 (odd q leads to heavy numerical instabilities as before):","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"closures = [\"zero\", \"normal\", \"log-normal\", \"gamma\"]\n\nplt_means = [plot() for i in 1:2]\nplt_vars = [plot() for i in 1:2]\n\nq_vals = [4, 6]\n\nfor (q, plt_m, plt_v) in zip(q_vals, plt_means, plt_vars)\n\n eqs = generate_central_moment_eqs(rn, 3, q, combinatoric_ratelaws=false)\n for closure in closures\n\n closed_eqs = moment_closure(eqs, closure)\n\n u₀map = deterministic_IC(u₀, closed_eqs)\n oprob = ODEProblem(closed_eqs, u₀map, tspan, p)\n\n sol = solve(oprob, Tsit5(), saveat=0.1)\n plt_m = plot!(plt_m, sol, idxs=[1], label = closure) \n plt_v = plot!(plt_v, sol, idxs=[4], label = closure)\n\n end\n\n plt_m = plot!(plt_m, means_SSA.t, means_SSA[1,:], title=\"m = 3, q = \"*string(q),\n linestyle=:dash, label = \"SSA\", color=\"gray\", legend=false)\n\n plt_v = plot!(plt_v, vars_SSA.t, vars_SSA[1,:], linestyle=:dash, label = \"SSA\", color=\"gray\", legend=false)\n\nend\n\nplt_means[1] = plot(plt_means[1], ylabel = \"Mean p53 molecule number\")\nplt_vars[1] = plot(plt_vars[1], ylabel = \"Variance of p53 molecule number\", legend=:topleft)\nplot(plt_means..., plt_vars..., size=(1250, 750), lw=1.5, xlabel=\"Time [h]\",\n guidefontsize=10, titlefontsize=12, legendfontsize=8, leftmargin=4mm, bottommargin=2mm)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 3rd order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"The results are rather interesting: the mean trajectories initially show the characteristic dampened oscillations, which, however, fail to converge to a steady-state value and instead start growing in amplitude as the time increases—this behaviour is also more pronounced for higher q. Surprisingly, log-normal and gamma closures still provide fairly accurate variance estimates. Note that zero closure ODEs fail to be evaluated for the entire time course at all, displaying strongly oscillatory mean trajectories and even exponentially growing variances. Therefore, third order moment expansion and all of the applied closure methods fail to give physically meaningful results.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"This outcome seems concerning as Lakatos et al. (2015) state that the same closures reached steady-state at extended simulation time (using identical parameter set). In order to investigate this discrepancy further, we turn to Python MEANS package, the implementation of which is fully based on the technical details presented in Lakatos et al. (2015). Note that the package has not been maintained for years, only works on Python 2 and setting up its dependencies may be fiddly. Nevertheless, having done that, we can generate and solve the moment equations up to third order running the following code in a separate Python 2 notebook:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"import means\nimport means.examples\nimport numpy as np\nimport matplotlib.pyplot as plt\n\nsys = means.examples.MODEL_P53\n\nps = [90, 0.002, 1.7, 1.1, 0.93, 0.96, 0.01]\nic_values = [70, 30, 60]\n\ntspan = np.arange(0, 200, 0.1)\n\nplt.figure(figsize=(10,6))\n\nfor closure in [\"normal\", \"log-normal\", \"gamma\"]:\n odes = means.mea_approximation(sys, 3, closure=closure, multivariate=True)\n simulation = means.Simulation(odes)\n trajectories = simulation.simulate_system(ps, ic_values, tspan)\n trajectories[0].plot(label=closure)\n\nplt.margins(x=0)\nplt.xlabel('Time [h]', fontsize=12)\nplt.ylabel('Mean p53 molecule number', fontsize=12)\nplt.tick_params(axis='both', which='major', labelsize=10)\nplt.legend(fontsize=12)\nplt.show()","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 3rd order expansion MEANS)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We observe that the mean trajectories fail to converge after all, hinting that some of the results published in Lakatos et al. (2015) may benefit from a more critical analysis. Note that MEANS does not allow arbitrary q values and automatically sets q = m+1, hence third-order moment expansion using MEANS is equivalent to our m=3 q=4 case.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"We stress that the obtained mean trajectories using the two packages are not identical due to different implementations of the closure methods (different solvers do not have a significant effect in this case). Lakatos et al. (2015) (and hence MEANS) formulate the normal, log-normal and gamma closure functions exclusively in terms of first and second order moments—higher order moments, whose corresponding ODEs are being solved for simultaneously, are not included in the construction of closure functions. One could say that such formulation makes a stronger assumption about the underlying distribution and does not utilise all the information available about the system—our implementation takes that into account. For example, in case of m=3 q=4 moment expansion, we approximate the fourth order central moments in terms of third and lower order moments, whereas MEANS would approximate it using only first and second order moments. Curiously, MEANS formulation appears to be more numerically stable (shows weaker oscillations) at longer times in this particular case, as can be seen from the plots above. Nevertheless, as both approaches failed to produce physically meaningful results, we do not investigate the numerical differences further—one would still generally expect that the inclusion of the computed higher-order moment values will improve the numerical stability of the ODEs and the accuracy of moment estimates.","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Finally, we can check whether better estimates can be obtained using even higher, fifth, order expansion with q=6:","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"plt = plot()\nclosures = [\"zero\", \"normal\", \"log-normal\", \"gamma\"]\n\neqs = generate_central_moment_eqs(rn, 5, 6, combinatoric_ratelaws=false)\n# faster to store than recompute in case we want to try different solvers/params\noprobs = Dict()\n\nfor closure in closures\n closed_eqs = moment_closure(eqs, closure)\n\n u₀map = deterministic_IC(u₀, closed_eqs)\n oprobs[closure] = ODEProblem(closed_eqs, u₀map, tspan, p)\n sol = solve(oprobs[closure], Tsit5(), saveat=0.1)\n\n plt = plot!(plt, sol, idxs=[1], label = closure) \nend\n\nplt = plot!(plt, xlabel = \"Time [h]\", ylabel = \"Mean p53 molecule number\")\nplt = plot!(plt, means_SSA.t, means_SSA[1, :], linestyle=:dash, label = \"SSA\", color=\"gray\")\nplot(plt, size=(750, 450), lw=2)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 5th order expansion)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"All closures are highly numerically unstable and fail to provide physically meaningful results. Performing similar fifth order moment expansion using the MEANS package we obtain (full code can be found in the example Jupyter notebook):","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"(Image: P53-Mdm2 5th order expansion MEANS)","category":"page"},{"location":"tutorials/P53_system_example/","page":"P53 System Example","title":"P53 System Example","text":"Note that normal closure is not included in the plot as it failed to converge using MEANS while log-normal and gamma closures seem to be much more stable in comparison to MomentClosure results (albeit still unphysical). In contrast, normal closure is the only one which succeeded in being evaluated for the entire time course using MomentClosure (full time not shown in the plot) which again is a result of the implementational differences between the two packages.","category":"page"},{"location":"tutorials/SIR_example/#SIR_example","page":"SIR Model Example","title":"SIR Model Example","text":"","category":"section"},{"location":"tutorials/SIR_example/","page":"SIR Model Example","title":"SIR Model Example","text":"See here for an example application of MomentClosure to the SIR model (thanks to @sdwfrost).","category":"page"},{"location":"theory/moment_expansion_SDE/#moment_expansion_SDE","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"","category":"section"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"Moment equations can be also written down for models described by N-dimensional Stochastic Differential Equations (SDEs) of the form ","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n dmathbfn = mathbff left(mathbfn t right) dt + mathbfGleft( mathbfn t right) dmathbfw_t\nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"where mathbfn=(n_1 dotsc n_N) is a N times 1 state vector, mathbff left(mathbfn t right) = left( f_1(mathbfn t) dotsc f_N(mathbfn t) right)^top is the deterministic (drift) term, mathbfG(mathbfn t) is a N times m (diffusion) matrix and mathbfw_t is a m times 1 vector of Wiener processes so that langle dmathbfw_t rangle = mathbf0 and dmathbfw_t dmathbfw_t^top = mathbfIdt, with mathbfI being an m times m identity matrix. Note that the notation used throughout is consistent with that of the previous section, and the state vector here is the same as for a chemical reaction network of N species (the only difference being that its elements are now continuous variables). The derivations below are based on Ghusinga et al. (2017) [1] and Bover (1978) [2].","category":"page"},{"location":"theory/moment_expansion_SDE/#raw_moment_eqs_SDE","page":"Moment Expansion (SDE)","title":"Raw Moment Equations","text":"","category":"section"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"To obtain the raw moment equations, we begin with Itô's lemma stating that for any smooth scalar function h(mathbfn) one has ","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n dh(mathbfn) = left sum_i=1^N fracpartial h(mathbfn) partial n_i f_i + frac12 sum_i=1^N sum_j=1^N fracpartial^2 h(mathbfn) partial n_i partial n_j (mathbfGmathbfG^top)_ij rightdt + sum_i=1^N sum_j=1^N fracpartial h(mathbfn) partial n_i mathbfG_ij dw_j \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"Taking the expectation of both sides we obtain","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracd langle h(mathbfn) rangledt = sum_i=1^N leftlangle fracpartial h(mathbfn)partial n_i f_i rightrangle + frac12 sum_i=1^N sum_j=1^N leftlangle fracpartial^2 h(mathbfn)partial n_i partial n_j (mathbfGmathbfG^top)_ijrightrangle \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"where we have used the fact that langle dmathbfw rangle = mathbf0. Finally, the general form of raw moment equations up to m^textth order is obtained by setting h(mathbfn) = mathbfn^mathbfi = n_1^i_1 dotsm n_N^i_N, so that langle h(mathbfn) rangle = mu_mathbfi and mathbfi = sum_j=1^N i_j leq m, leading to ","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracdmu_mathbfidt = sum_j=1^N leftlangle fracpartial mathbfn^mathbfi partial n_j f_j rightrangle + frac12 sum_j=1^N sum_k=1^N leftlangle fracpartial^2 mathbfn^mathbfi partial n_j partial n_k (mathbfGmathbfG^top)_jkrightrangle \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"The equation for the mean vector mathbfμ = langle mathbfn rangle = (mu_1 dotsc mu_N) take a particularly simple form (as the second derivatives are zero):","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracdmathbfμdt = langle mathbff(mathbfn t) rangle\nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/#central_moment_eqs_SDE","page":"Moment Expansion (SDE)","title":"Central Moment Equations","text":"","category":"section"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"Note that the raw moment expansion formulated above is valid only when the SDE terms mathbff left(mathbfn t right) and mathbfG(mathbfn t) are polynomial functions. If these functions are non-polynomial, we have to consider central moment equations and Taylor expand both mathbff left(mathbfn t right) and mathbfG(mathbfn t) around the mean mathbfμ up to the expansion order q (as done in case of the CME with non-polynomial propensities).","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"The corresponding Taylor expansions can be written down as:","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n f_j(mathbfn) = sum_mathbfk=0^q frac1mathbfkD^mathbfk f_j(mathbfμ) (mathbfn - mathbfμ)^mathbfk + dotsb \n left( mathbfG(mathbfn) mathbfG(mathbfn)^top right)_jk = sum_mathbfl=0^q frac1mathbfl D^mathbfl left( mathbfG(mathbfμ)mathbfG(mathbfμ)^top right)_jk (mathbfn - mathbfμ)^mathbfl + dotsb \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"and we remind the reader that we have previously introduced the notation","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\nmathbfk = k_1dotsm k_N \nD^mathbfk f = fracpartial^mathbfkpartial n_1^k_1 dotsm partial n_N^k_N f\n(mathbfn - mathbfμ)^mathbfk = (n_1-mu_1)^k_1 dotsm (n_N-mu_N)^k_N \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"It follows that the equations for the means can now be written down as:","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracdmu_idt = langle f_i rangle = sum_mathbfj=0^q frac1mathbfj D^mathbfj f_i(mathbfμ) M_mathbfj + dotsb \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"where M_mathbfj = langle (mathbfn-mathbfμ)^mathbfj rangle are central moments.","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"To obtain the central moment equations, we first perform a change of variables mathbfy = mathbfn - langle mathbfn rangle and note again that dlangle mathbfn rangledt = langle mathbff(mathbfn t) rangle, so that the SDE becomes","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n dmathbfy = left mathbffleft(mathbfy + langlemathbfnrangle t right) - leftlangle mathbff left(mathbfy + langlemathbfnrangle t right) rightrangle right dt + mathbfGleft(mathbfy + langlemathbfnrangle t right) dmathbfw_t \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"Then, using Itô's lemma for a smooth function h(mathbfy) and taking the expectation of both sides we find","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracd langle h(mathbfy) rangledt = sum_i=1^N left leftlangle fracpartial h(mathbfy)partial y_i f_i rightrangle - leftlangle fracpartial h(mathbfy)partial y_i rightrangle leftlangle f_i rightrangle right + frac12 sum_i=1^N sum_j=1^N leftlangle fracpartial^2 h(mathbfy)partial y_i partial y_j (mathbfGmathbfG^top)_ijrightrangle \nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"Letting h(mathbfy) = (mathbfn - mathbfμ)^mathbfi and substituting in the needed Taylor expansions we obtain the final form of central moment equations:","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"beginalign*\n fracd M_mathbfidt = sum_j=1^N Bigg leftlangle fracpartial (mathbfn - mathbfμ)^mathbfi partial (n_j - mu_j) sum_mathbfk=0^q-mathbfi+1 frac1mathbfk D^mathbfk f_j(mathbfμ) (mathbfn-mathbfμ)^mathbfk + dotsb rightrangle \n - leftlangle fracpartial (mathbfn-mathbfμ)^mathbfi partial(n_j - μ_j) rightrangle leftlangle sum_mathbfk=0^q frac1mathbfk D^mathbfk f_j(mathbfμ) (mathbfn-mathbfμ)^mathbfk + dotsb rightrangle Bigg \n + frac12 sum_j=1^N sum_k=1^N leftlangle fracpartial^2 (mathbfn-mathbfμ)^mathbfipartial(n_j-mu_j)partial(n_k-mu_k) sum_mathbfl=0^q-mathbfi+2 frac1mathbfl D^mathbfl left( mathbfG(mathbfμ) mathbfG(mathbfμ)^top right)_jk (mathbfn-mathbfμ)^mathbfl + dotsc rightrangle\nendalign*","category":"page"},{"location":"theory/moment_expansion_SDE/#References","page":"Moment Expansion (SDE)","title":"References","text":"","category":"section"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"[1]: K. R. Ghusinga, M. Soltani, A. Lamperski, S. V. Dhople, and A. Singh, \"Approximate moment dynamics for polynomial and trigonometric stochastic systems\", IEEE 56th Annual Conference on Decision and Control (2017). https://doi.org/10.1109/CDC.2017.8263922","category":"page"},{"location":"theory/moment_expansion_SDE/","page":"Moment Expansion (SDE)","title":"Moment Expansion (SDE)","text":"[2]: D. C. C. Bover, \"Moment Equation Methods for Nonlinear Stochastic Systems\", Journal of Mathematical Analysis and Applications 65, 306-320 (1978). https://doi.org/10.1016/0022-247X(78)90182-8","category":"page"},{"location":"theory/moment_expansion_CME/#moment_expansion_CME","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"","category":"section"},{"location":"theory/moment_expansion_CME/#chemical_master_equation","page":"Moment Expansion (CME)","title":"The Chemical Master Equation","text":"","category":"section"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Consider a chemical reaction network with N different molecular species X_i (i=1 dotsc N) and R chemical reactions, so that the system can be described by [1]:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n sum_i=1^N s_ij X_i xrightarrowk_j sum_i=1^N r_ij X_i quad j=1dotscR\nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"where s_ij and r_ij respectively denote the numbers of reactant and product molecules of species i in the chemical reaction j. The stoichiometric matrix is defined as S_ij = r_ij - s_ij (the net change in the number of molecules of species X_i when the r^textth reaction occurs). The state of the system is determined by the state vector mathbfn=(n_1 dotsc n_N), where n_i is the number of X_i molecules. The time evolution of the probability distribution of mathbfn is described by the Chemical Master Equation (CME)","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginequation\n fracd P(mathbfn t)dt = sum^R_r=1 Big a_r(mathbfn-S_r)P(mathbfn-S_r t) - a_r(mathbfn) P(mathbfn t) Big tag1\nendequation","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"where a_r(mathbfn) is the propensity function of the r^textth reaction and S_r is the r^textth column of the stoichiometric matrix S.","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Modelling using the CME framework is common in the study of biochemical and gene networks within cells. Although such master equations are rather simple in structure, for most systems there are no known analytical solutions and their stochastic simulations can be very computationally expensive. One possible approach to investigate the system at hand is to approximate the whole probability distribution solution of the CME in terms of its first few moments (e.g. mean and variance).","category":"page"},{"location":"theory/moment_expansion_CME/#raw_moment_eqs","page":"Moment Expansion (CME)","title":"Raw Moment Equations","text":"","category":"section"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"From Eq. (1) we can obtain a system of ordinary differential equations (ODEs) governing the time-evolution of the raw moments of the system up to specified moment expansion order m, given by mu_mathbfi = mu_i_1 dotsc i_N = langle n_1^i_1 dotsm n_N^i_N rangle, where the indices mathbfi=(i_1 dotsc i_N) are restricted such that mathbfi = sum_j=1^N i_j leq m. Note that first order raw moments are simply the means. For example, when N = 3, the mean molecule number of second species in the system is given by mu_2 = mu_0 1 0 = langle n_2 rangle (when N=3). We chose to relax the notation throughout so that the means can be indicated by a single index of the molecular species (used in this section) or by the corresponding one-hot vector (used in the code).","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"To find the first moment equations, we multiply the CME by n_i and sum over all possible states:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n sum_mathbfn n_i fracdP(mathbfnt)dt = sum_n_1^infty sum_n_2^infty dotsm sum_n_N^infty n_i fracdP(mathbfnt)dt \n = sum_r sum_mathbfn n_i a_r(mathbfn-S_r)P(mathbfn-S_r t) - n_i a_r(mathbfn) P(mathbfn t) \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Applying a transformation mathbfn-S_r rightarrow mathbfn on the first term in the sum leads to","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n fracdmu_idt = sum_r sum_mathbfn (n_i+S_ir) a_r(mathbfn)P(mathbfn t) - n_i a_r(mathbfn) P(mathbfn t) \n = sum_r sum_mathbfn S_ir a_r(mathbfn) P(mathbfn t) \n = sum_r S_ir langle a_r(mathbfn) rangle tag2\nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"The derivation can be extended to the multivariate case of higher order moments mu_mathbfi with mathbfi leq m:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n fracdmu_mathbfidt = sum_mathbfn n_1^i_1dotsm n_N^i_N fracdP(mathbfn t)dt \n = sum_r sum_mathbfn left (n_1+S_1r)^i_1dotsm (n_N +S_2r)^i_N - n_1^i_1dotsm n_N^i_N right a_r(mathbfn)P(mathbfn t) \n = sum_r sum_mathbfn left sum^i_1_j_1=0 binomi_1j_1 n_1^j_1S_1r^i_1-j_1 dotsm sum^i_N_j_N=0 binomi_Nj_N n_N^j_NS_Nr^ i_N-j_N - n_1^i_1dotsm n_N^i_N right a_r(mathbfn)P(mathbfn t) \n = sum_r sum_mathbfn sum^mathbfi-1_mathbfj=0 binomi_1j_1 dotsm binomi_Nj_N S_1r^i_1-j_1 dotsm S_Nr^ i_N-j_N n_1^j_1 dotsm n_N^j_N a_r(mathbfn)P(mathbfn t) \n = sum_r sum^mathbfi-1_mathbfj=0 binom mathbfi mathbfj S_r^mathbfi-mathbfj langle mathbfn^mathbfj a_r(mathbfn) rangle \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"where we have introduced multi-index notation:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign\n mathbfi = sum^N_j=1 i_j tag3 \n sum^mathbfi-1_mathbfj=0 = sum_substack0 leq j_1 leq i_1 dotsc 0 leq j_N leq i_N 0 leq j_1+dotsb+j_N leq mathbfi-1 tag4 \n binom mathbfi mathbfj = binomi_1j_1 dotsm binomi_Nj_N tag5 \n mathbfn^mathbfj = n_1^j_1dotsm n_N^j_N tag6 \n S_r^mathbfi-mathbfj = S_1r^ i_1-j_1 dotsm S_Nr^ i_N-j_N tag7\nendalign","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"It is crucial to stress three key points:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Throughout all derivations presented in this section we assume that the components of the net stoichiometry matrix S are constant values. However, in certain cases the reaction product may itself be a stochastic variable so that the corresponding expectation values, langle S_r^mathbfi-mathbfj rangle, must be taken into account. We have implemented a limited support for such systems where components of S can be independent geometrically distributed variables, see the specific gene network example.\nGeneration of raw moment equations is only possible if the kinetics of the system at hand are governed by the law of mass action where all propensity functions are polynomials in mathbfn, otherwise the expectation value terms langle a_r(mathbfn) rangle and langle n_i_1^j_1 dotsm n_i_N^j_N a_r(mathbfn) rangle are ill-defined. This key issue can be overcome by the general central moment expansion method presented below.\nThe order of the polynomials a_r(mathbfn) determines the order of moments encountered in the generated system of ODEs. If the system is linear (contains only zeroth and first order reactions), the m^textth order moments will depend only of moments of order m or lower, hence constituting a finite hierarchy of moment ODEs that can be readily solved numerically or otherwise without approximations. However, if the system is non-linear (involves second or higher order reactions), moment equations will depend on higher order moments. For example, if the reaction network contains bimolecular reactions, the corresponding propensity functions will be second order polynomials and hence m^textth order moment equations will now depend on (m+1)^th order moments. This leads to an infinite hierarchy of coupled moment equations where each moment will depend on higher order moments—it cannot be solved directly and needs to be truncated. This can be achieved using one of many moment closure approximations that express (m+1)^textth order moments in terms of m^textth and lower order moments using different distributional assumptions, effectively closing the hierarchy and enabling one to solve the moment equations up to m^textth order [1]. Details of all closure methods currently implemented within the package can be found in the next Theory section.","category":"page"},{"location":"theory/moment_expansion_CME/#central_moment_eqs","page":"Moment Expansion (CME)","title":"Central Moment Equations","text":"","category":"section"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"As we have seen, the moment equations of P(mathbfn t) can be obtained in a straightforward manner if the kinetics of the system are governed by the law of mass action where all propensity functions are polynomials in mathbfn [1]. Similarly, given all propensities are rational functions, a polynomial form can also be recovered [2]. However, problems arise when the propensities take more complicated non-polynomial functions. This can nevertheless be overcome by considering a more general method of moment expansion that enables us to obtain mean and central moment equations up to arbitrary order for virtually any chemical reaction network with any type of smooth (infinitely differentiable) propensity functions. Such framework was first independently formulated by Lee [3] and Ale et al. [4]—our derivation below closely follows these works.","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"We start by Taylor-expanding the propensity functions around the mean mathbfμ = langle mathbfn rangle = (mu_1 dotsc mu_N). This allows us to consider any general propensity function under the assumption that it is infinitely differentiable (smooth). The expansion leads to","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n a_r(mathbfn) = a_r(mathbfμ) + sum_j_1^N fracpartial a_r(mathbfμ)partial n_j_1(n_j_1-mu_j_1) \n + frac12sum_j_1 j_2 fracpartial^2 a_r(mathbfμ)partial n_j_1 partial n_j_2 (n_j_1-mu_j_1) (n_j_2-mu_j_2) + dotsb \n + frac1q sum_j_1 j_2 dotsc j_q fracpartial^q a_r(mathbfμ)partial n_j_1 partial n_j_2dotsmpartial n_j_q (n_j_1-mu_j_1) (n_j_2-mu_j_2) dotsm(n_j_q - mu_j_q) + dotsb \n = sum_mathbfj=0^q frac1mathbfj D^mathbfj a_r(mathbfμ) (mathbfn - mathbfμ)^mathbfj + dotsb \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"where q controls the expansion order and we have simplified the expression by using Eqs. (3-7) in addition to:","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\nmathbfj = j_1dotsm j_N \nD^mathbfj f = fracpartial^mathbfjpartial n_1^j_1 dotsm partial n_N^j_N f \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Now we can obtain equations governing the time-evolution of the means μ_i and the central moments M_mathbfi = M_i_1 dotsc i_N = langle (n_1 - mu_1)^i_1dotsm(n_N - mu_N)^i_N rangle, where again mathbfi leq m. We first consider the equations for the means immediately starting from Eq. (2):","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n fracdmu_idt = sum_r S_ir langle a_r(mathbfn) rangle \n = sum_r S_ir Big( sum_mathbfj=0^q frac1mathbfj D^ mathbfj a_r(mathbfμ) langle (mathbfn - mathbfμ)^mathbfj rangle + dotsb Big) \n = sum_r S_ir Big( sum_mathbfj=0^q frac1mathbfj D^ mathbfj a_r(mathbfμ) M_mathbfj + dotsb Big) \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"We can now derive the equations for central moments by taking the time derivative of M_mathbfi and using Eq. (1):","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n fracdM_mathbfidt = sum_mathbfn (n_1 - mu_1)^i_1dotsm(n_N - mu_N)^i_N fracdP(mathbfn t)dt \n - sum_mathbfn sum_j=1^N i_j fracdmu_jdt (n_1-mu_1)^i_1 dotsm(n_j-mu_j)^i_j-1dotsm(n_N-mu_N)^i_NP(mathbfn t) \n = sum_r sum_mathbfn Big (n_1 - mu_1 + S_1r)^i_1 dotsm (n_N-mu_N+S_Nr)^i_N \n - (n_1-mu_1)^i_1dotsm(n_N-mu_N)^i_N Biga_r(mathbfn)P(mathbfn t) - sum_j i_j fracdmu_jdtM_i_1dotsci_j-1dotsc i_N \n = sum_r sum_mathbfn Bigg sum_mathbfj=0^mathbfi-1 binomi_1j_1 dotsm binomi_Nj_N S_1r^ i_1-j_1 dotsm S_Nr^ i_N-j_N(n_1-mu_1)^j_1dotsm(n_N-mu_N)^j_NBigg \n times Bigg sum_mathbfj=0^q frac1mathbfj D^mathbfj a_r(mathbfμ) (mathbfn - mathbfμ)^mathbfj + dotsb Bigg-sum_j i_j fracdmu_jdt M_mathbfi_j- \n = sum_r sum_mathbfj=0^mathbfi-1 Bigg binomi_1j_1 dotsm binomi_Nj_N S_1r^ i_1-j_1 dotsm S_Nr^ i_N-j_N sum_mathbfk=0^q-mathbfj frac1mathbfkD^mathbfk a_r(mathbfμ) M_mathbfj+mathbfk + dotsb Bigg - sum_j i_j fracdmu_jdt M_mathbfi_j- \n = sum_r sum_mathbfj=0^mathbfi-1 Bigg binommathbfimathbfj S_r^mathbfi-mathbfj sum_mathbfk=0^q-mathbfj frac1mathbfkD^mathbfk a_r(mathbfμ) M_mathbfj+mathbfk + dotsb Bigg - sum_j i_j fracdmu_jdt M_mathbfi_j- \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"where we have also defined","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"beginalign*\n M_mathbfi_j-=M_i_1dotsci_j-1dotsc i_N \n M_mathbfj+mathbfk = M_j_1+k_1dotscj_N+k_N \nendalign*","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"Although raw moment equations for non-linear mass-action systems already require approximate treatment using moment closure, here we have an additional complication: if a system contains non-polynomial propensity functions, the equations for both means and central moments will in principle depend on an infinite number of higher order central moments. Hence the Taylor expansion order q is of utmost importance as it controls the degree of approximation: the propensity functions are expanded up to q^textth order so that m^textth central moment equations will depend on central moments of order qm and lower. Finally, moment closure approximations can be applied similarly as in the case of raw moment equations.","category":"page"},{"location":"theory/moment_expansion_CME/#References","page":"Moment Expansion (CME)","title":"References","text":"","category":"section"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"[1]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Approximation and inference methods for stochastic biochemical kinetics - a tutorial review\", Journal of Physics A: Mathematical and Theoretical 50, 093001 (2017). https://doi.org/10.1088/1751-8121/aa54d9","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"[2]: P. Milner, C. S. Gillespie, and D. J. Wilkinson, \"Moment closure approximations for stochastic kinetic models with rational rate laws\", Mathematical Biosciences 231, 99-104 (2011). https://doi.org/10.1016/j.mbs.2011.02.006","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"[3]: A. Ale, P. Kirk, and M. P. H. Stumpf, \"A general moment expansion method for stochastic kinetic models\", The Journal of Chemical Physics 138, 174101 (2013). https://doi.org/10.1063/1.4802475","category":"page"},{"location":"theory/moment_expansion_CME/","page":"Moment Expansion (CME)","title":"Moment Expansion (CME)","text":"[4]: C. H. Lee, \"A Moment Closure Method for Stochastic Chemical Reaction Networks with General Kinetics\", MATCH Communications in Mathematical and in Computer Chemistry 70, 785-800 (2013). https://match.pmf.kg.ac.rs/electronicversions/Match70/n3/match70n3785-800.pdf","category":"page"},{"location":"theory/moment_closure_approximations/#moment_closure_approximations","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"In the previous section, we have shown that for a non-linear system an infinite hierarchy of coupled moment equations is obtained that cannot be solved directly and, therefore, needs to be truncated in an approximate way. This can be achieved using moment closure approximations (MAs), in which all moments above a certain order m are expressed in terms of m^textth and lower order moments using various (usually distributional) assumptions [1]. Doing so enables us to effectively close the moment hierarchy, leading to a finite set of ODEs which can then be solved numerically. In this section, we present some the commonly used MA methods that are implemented in MomentClosure. Please see the Tutorials section for examples showing different MAs applied to a variety of systems.","category":"page"},{"location":"theory/moment_closure_approximations/#zero_closure","page":"Moment Closure Approximations","title":"Zero Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The simplest MA is the \"central-moment-neglect\" MA (CMN-MA) [2], also referred to as \"zero-closure\" [3] or \"low dispersion moment closure\" [4], where CMN-MA at m^textth order means that the moment equations are truncated by setting all central moments above order m to zero. For example, in the simple case of 2nd order truncation, the moment equations for the means mu_i and covariances C_ij become:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n fracd mu_idt = sum_r S_ir Big( a_r(mathbfμ) + frac12sum_i_1 i_2 fracpartial^2 a_r(mathbfμ)partial n_i_1 partial n_i_2 M_i_1 i_2 Big) \n fracdC_ijdt = fracdM_mathbf0_i+ j+dt \n = fracdlangle (n_i-mu_i) (n_j-mu_j) rangledt = \n = sum_r Big S_ir sum_k fracpartial a_r(mathbfμ)partial n_k C_jk + S_jr sum_k fracpartial a_r(mathbfμ)partial n_k C_ik \n + S_irS_jr Big( a_r(mathbfμ) + frac12 sum_k l fracpartial^2 a_r(mathbfμ)partial n_k partial n_l C_kl Big) Big \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/#normal_closure","page":"Moment Closure Approximations","title":"Normal Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Another popular MA is the \"normal moment-closure\", pioneered by Goodman [5] and Whittle [6], where all cumulants kappa_mathbfi above order m are set to zero, approximating the probability distribution of the system with the normal distribution [2]:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n kappa_mathbfi = 0 quad textfor mathbfi m\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"In order to truncate the higher order central or raw moments M_mathbfi using normal closure we express them in terms of cumulants kappa_mathbfi using the multivariate moment and cumulant relationships formalised by Balakrishnan et al. [7].","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Note that a different implementation of normal closure can be found in literature [3], where the higher order central moments are expressed in terms of a sum of product of covariances using Isserlis' theorem. However, one could argue that such formulation is not advisable as it assumes stronger \"Gaussianity\" of the underlying distribution than setting the higher order cumulants to zero which is less of an approximation on the form of the distribution and hence is preferable in the development of MAs. For example, in case we are truncating the moment equations at 4th order, the truncation-order central moments would be expressed only in terms of covariances whereas our formulation using cumulants would explicitly include information about the computed values of third central moments, which is expected to improve numerical stability and lead to more accurate moment estimates.","category":"page"},{"location":"theory/moment_closure_approximations/#poisson_closure","page":"Moment Closure Approximations","title":"Poisson Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Although the Poisson distribution lacks a general formulation for multiple variables [3], \"Poisson MA\" has been formulated [2, 8] assuming that the joint multivariate distribution is a product of univariate Poisson distributions, i.e., n_i sim textPoisson(mu_i). The cumulants of a univariate Poisson distribution are equal to the mean, hence in Poisson closure we set all higher order diagonal cumulants to the corresponding mean values and mixed cumulants to zero [2], which in our notation can be expressed as:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n kappa_mathbfi = mu_j quad textif i_j m textand i_jneq k = 0 quad textfor some j k in 1dotscN \n kappa_mathbfi = 0 quad textif mathbfi m textand i_j neq i_k (textwhere i_j neq 0 i_k neq 0) quad textfor some j k in 1dotscN\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Similarly to normal closure, the higher order central/raw moments can be expressed in terms of cumulants as described in [7].","category":"page"},{"location":"theory/moment_closure_approximations/#log-normal_closure","page":"Moment Closure Approximations","title":"Log-normal Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"\"Log-normal\" MA, first applied by Keeling [9], allows one to truncate the moment equations under the assumption that the distribution of the underlying stochastic process is log-normal. A positive multi-dimensional random variable mathbfn follows a log-normal distribution if its logarithm is normally distributed, so that mathbfy = ln mathbfn, and mathbfy sim mathcalN(mathbfnu Sigma), where mathbfnu and Sigma denote the vector of means and the covariance matrix respectively. By considering the moment generating function of the normal distribution, mathcalN(mathbfnu Sigma), one can show that the raw moments are given by [3, 10]:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n mu_mathbfi = expleft(mathbfi^topmathbfnu + frac12mathbfi^topSigmamathbfiright)\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"It follows that","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n nu_i = ln mu_i - frac12 Sigma_ii \n Sigma_ij = ln left( 1 + fraclangle (n_i - mu_i)(n_j - mu_j) rangleexp left( nu_i + nu_j + frac12 left( Sigma_ii + Sigma_jj right) right) right) \n = ln left( 1 + fracC_ijexp left( nu_i + nu_j + frac12 left( Sigma_ii + Sigma_jj right) right) right) \n Sigma_ii = ln left( 1 + fracC_iimu_i^2 right) \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Note that central moments can be obtained from raw moments by utilising their general multivariate relationship [11].","category":"page"},{"location":"theory/moment_closure_approximations/#gamma_closure","page":"Moment Closure Approximations","title":"Gamma Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The method of \"gamma closure\" was originally implemented by Lakatos et al. [3], where the authors acknowledged the ambiguity arising in defining multivariate gamma distributions, and, building upon previous definitions in the literature (e.g. [12] and [13]), proposed a new formulation of a multivariate gamma distribution. Here we reproduce the definition by closely following the description in [3] and elucidating some of the derivation steps.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"We denote a random variable drawn from gamma distribution with shape alpha and scale beta as n sim textGamma(alpha beta). The probability density function of the univariate gamma distribution with the corresponding shape-scale parameterisation is","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n f(n alpha beta) = frac1Gamma(alpha)beta^alpha n^alpha-1e^-fracnbeta\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where Gamma is the gamma function. The i^textth raw moment of n is given by","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n μ_i = fracGamma(alpha+i)beta^iGamma(alpha) = (alpha)_i beta^i taga\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where (alpha)_i = alpha (alpha+1) dotsm (alpha+i-1). Note that the moment generating function of X is","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n G_n(k) = langle e^kn rangle = left( 1 - beta k right)^-alpha\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"In order to construct a multivariate gamma distribution, we start by considering independent gamma variables Y_kl, k l = 1 dotsc N, with shape and scale parameters alpha_kl and beta_kl respectively. Here we define Y_kl to be symmetric in indices, i.e., Y_kl = Y_lk. Now consider an N-dimensional random vector mathbfn = left( n_1 n_2 dotsc n_N right), where n_i is a linear combination of independent gamma variables:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n n_i = sum_j=1^N fracbeta_iibeta_ij Y_ij\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The i^textth marginal moment generating function of the joint distribution of mathbfn is given by:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n G_n_i(k_i) = leftlangle e^k_i sum_j=1^N fracbeta_iibeta_ij rightrangle \n = G_Y_i1left(k_i fracbeta_iibeta_i1 right) G_Y_i2left(k_i fracbeta_iibeta_i2 right) dotsm G_Y_iNleft(k_i fracbeta_iibeta_iN right) \n = left(1-beta_ii k_i right)^sum_j=1^N alpha_ij\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"so that n_i sim textGamma(alpha_i beta_i), where alpha_i = sum_j=1^N alpha_ij and beta_i = beta_ii. Therefore, we have obtained an N-variate gamma distribution, which can be denoted as mathbfn sim MG(mathbfalpha mathbfbeta), where the vectors of shape and scale parameters are given by mathbfalpha = left( alpha_1 dotsc alpha_N right) and mathbfbeta = left( beta_1 beta_2 dotsc beta_N right) respectively.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"We can now readily obtain the raw m^textth order moment of n_i:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n langle n_i^m rangle = leftlangle left( sum_j=1^N fracbeta_iibeta_ij Y_ij right)^m rightrangle \n = leftlangle sum_k_1+k_2+dotsb+k_N = m fracmk_1k_2 dotsm k_N prod_j=1^N left( fracbeta_iibeta_ij Y_ij right)^k_j rightrangle \n = sum_mathbfk=m fracmmathbfk prod_j=1^N leftlangle left( fracbeta_iibeta_ij Y_ij right)^k_j rightrangle \n = beta_i^m sum_mathbfk=m fracmmathbfk left(\n prod_j=1^N left( alpha_ij right)_k_j right)\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Note that we have used Eq. (a) to get to the last line. The mixed raw moments are computed in a similar fashion:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n mu_mathbfi = leftlangle left( sum_j=1^N fracbeta_11beta_1jY_1j right)^i_1 dotsm left( sum_j=1^N fracbeta_NNbeta_NjY_Nj right)^i_N rightrangle \n = mathbfbeta^mathbfi leftlangle sum_mathbfk_1=i_1 dotsm sum_mathbfk_N=i_N fracmathbfimathbfk_1dotsmmathbfk_N prod_j=1^N left( fracY_1jbeta_1j right)^k_1_j dotsm prod_j=1^N left( fracY_Njbeta_Nj right)^k_N_j rightrangle \n = mathbfbeta^mathbfi leftlangle sum_mathbfk_1=i_1 dotsm sum_mathbfk_N=i_N fracmathbfimathbfk_1dotsmmathbfk_N prod_q=1^N left( fracY_qqbeta_qq right)^k_q_q prod_r=q+1^N left( fracY_qrbeta_qr right)^k_q_r + k_r_q rightrangle \n = mathbfbeta^mathbfi sum_mathbfk_1=i_1 dotsm sum_mathbfk_N=i_N fracmathbfimathbfk_1dotsmmathbfk_N prod_q=1^N left( alpha_qq right)_k_q_q prod_r=q+1^N left( alpha_qr right)_k_q_r + k_r_q\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where we have taken into account the symmetry in indices and defined each mathbfk_i as an N-dimensional vector mathbfk_i = left( k_i_1 dotsc k_i_N right) and mathbfbeta^mathbfi = beta_1^i_1 dotsm beta_N^i_N. Note that the mean and variance of n_i can be obtained from:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n mu_i = sum_j=1^N alpha_ij beta_i tagb \n C_ii = (alpha_i)_2 beta_i^2 - alpha_i^2beta_i^2 \n = sum_j=1^N alpha_ij beta_i^2\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Similarly, from Eq. (a) it follows that","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n langle n_i n_j rangle = beta_i beta_j left( sum_substackk l (k l) neq (j i)^N alpha_ik alpha_jl + (alpha_ij)_2 right)\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"which together with Eq. (b) allows us to express the covariance as:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n C_ij = langle n_i n_j rangle - beta_i beta_j left( sum_kl alpha_ik alpha_jl right) \n = alpha_ij beta_i beta_j\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Finally, from the equations above we can obtain all shape and scale parameters:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n beta_i = fracC_iimu_i \n alpha_ij = fracC_ijbeta_i beta_j \n alpha_ii = fracmu_ibeta_i - sum_substackk k neq i alpha_ik\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/#derivative_matching","page":"Moment Closure Approximations","title":"Derivative Matching","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The derivative matching MA [14, 15] is based on expressing moments above order m in terms of lower order moments in such a way that their time derivatives match those of the exact moments at some initial time and initial conditions. We outline the approach below, closely following the complete exposition found in the original papers of Singh and Hespanha [14, 15].","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The raw moment equations up to order m for any mass-action reaction network containing at most bimolecular (second order) reactions can be written down concisely in the matrix form:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n fracdmathbfμdt = hatmathbfa + Amathbfμ + B barmathbfμ \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where mathbfμ is a vector containing all raw moments of the system up to order m and barmathbfμ consists of all (m+1)^textth order raw moments which the equations depend on. The constant vector hatmathbfa and constant matrices A and B are chosen appropriately for the system at hand. In this case, an MA can be defined as a procedure where each moment in barmathbfμ, barμ_mathbfi, is approximated by a moment closure function varphi_mathbfi (mathbfμ) of moments up to order m. Then the moment equations can be rewritten as","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":" fracdmathbfνdt = hatmathbfa + Amathbfν + Bbarmathbfφ(mathbfν) ","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where the state of the system is now denoted by mathbfν instead of mathbfμ, stressing the fact that we are considering the approximation of the true moment dynamics, and barmathbfφ(mathbfν) is the corresponding vector of moment closure functions.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The idea behind derivative matching is to determine a map barmathbfφ so that the time derivatives between the exact moments, mathbfμ(t), and the approximate moments, mathbfν(t), would match at some initial time t_0 under the initial condition mathbfμ(t_0) = mathbfν(t_0):","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n left fracd^i mathbfμdt rightrvert_t=t_0 = left fracd^i mathbfνdt rightrvert_t=t_0\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"If these conditions hold, one can expect from a Taylor series approximation argument that mathbfμ(t) and mathbfν(t) will stay close at least locally in time and hence the MA will be sufficiently accurate.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"In order to move forward, Singh and Hespanha present what can be understood as essentially an ansatz. Firstly, moment closure functions for each mathbfi, where mathbfi m, are chosen to have a separable form given by","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n φ_mathbfi(mathbfμ) = prod_substack1 leq j_1+dotsb+j_N leq m left( μ_mathbfj right)^γ_mathbfj = prod_mathbfj=1^m left( μ_mathbfj right)^γ_mathbfj \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where γ_mathbfj are constants (unique for each vector mathbfi) that can be determined by solving the following linear equation system:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n C^mathbfi_mathbfj = sum_mathbfk=1^m gamma_mathbfk C^mathbfk_mathbfj quad textfor each mathbfj textwhere mathbfjleq m \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"were we have introduced multi-index scalars","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n C^mathbfu_mathbfv = C^u_1_v_1C^u_2_v_2 dotsm C^u_N_v_N\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"with each element defined as","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nC^a_b =\nbegincases\n fraca(a-b)b a geq b \n 0 a lt b\nendcases\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Using the specific construction of barmathbfφ described above, it can be shown [15] that for every deterministic initial condition, i.e., mathbfn(t_0) = mathbfn_0 with probability one, we will have","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n mathbfμ(t_0) = mathbfν(t_0) implies left fracd mathbfμdt rightrvert_t=t_0 = left fracd mathbfνdt rightrvert_t=t_0 \n implies left fracd^2mathbfμdt^2 rightrvert_t=t_0 = left fracd^2 mathbfνdt^2 rightrvert_t=t_0 + mathbfϵ(mathbfn_0) \nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"where all elements of mathbfϵ(mathbfn_0) are zero except the ones corresponding to m^textth order raw moments—these elements are second order polynomials in mathbfn_0. Note, however, that these results hold only for mass-action systems containing no higher than second order chemical reactions. While the derivative matching MA can be applied in the same way to systems containing higher order polynomial and non-polynomial propensity functions, it has not been rigorously analysed in such scenarios, where, naturally, we expect significantly larger approximation errors.","category":"page"},{"location":"theory/moment_closure_approximations/#conditional_closures","page":"Moment Closure Approximations","title":"Conditional Closures","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"As standard MAs often fail to provide sufficiently accurate approximations of chemical reaction networks involving both high- and low-abundance species, some novel approaches suggest using moments conditioned on the low-copy number species, which can lead to a more effective description of the system dynamics [1, 16]. Here we discuss the conditional moment closure proposed by Soltani et al. (2015) [16], applicable to networks containing molecular species which copy number can be either zero or one, i.e., a binary/Bernoulli random variable. Such conditional MA can be very useful in the study of gene networks where two-state gene systems are often considered—the gene state itself can be treated as a distinct species which molecule number is a Bernoulli variable.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The conditional MA is based on conditioning the higher order moments of high-abundance species on the binary species being in state 1 (instead of 0) and then applying standard MAs on the conditional moments. Closely following Soltani et al. [16], we use a two-state gene circuit to illustrate the conditional MAs, denoting the binary gene state by g and the protein number (high-abundance species) by p. Firstly, note that as g is a Bernoulli variable the following identities hold","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle g^j p rangle = langle g rangle quad j in 2 3 dotsc \nlangle g^j p^k rangle = langle g p^k rangle quad jk in 1 2 3 dotsc\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Therefore, we only need to concern ourselves with moments of the form langle p^k rangle and langle gp^k rangle. The former can be approximated using the standard MAs (no conditioning needed), whereas the latter can be written down as:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle gp^j rangle = langle p^j g=1 rangle langle g rangle quad jin 12 dotsc\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Now the conditional moment langle p^j g= 1 rangle can be expressed in terms of lower order conditional moments using one of the standard MAs, e.g., normal closure or derivative matching—the two methods (including the conditioning step) are respectively known as the conditional gaussian and conditional derivative matching MAs. For completeness, we show how the two methods are applied to approximate a specific higher order moment:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n langle gp^3 rangle = langle p^3 g=1 rangle langle g rangle\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/#conditional_gaussian_closure","page":"Moment Closure Approximations","title":"Conditional Gaussian Closure","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Under the conditional Gaussian MA, we assume that the number of protein molecules conditioned on the gene being active, p g = 1, follows a Gaussian distribution. In other words, we apply normal MA on the conditional moment langle p^j g=1 rangle. Hence we obtain (following the example above):","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle p^3 g= 1 rangle = 3 langle p^2 g = 1 rangle langle p g=1 rangle - 2 langle p g=1 rangle^3\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Using langle gp^j rangle = langle p^j g=1 rangle langle g rangle, we can rewrite the equation as:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle p^3 g= 1 rangle = 3 fraclangle gp^2 rangle langle gp ranglelangle g rangle^2 - 2 fraclangle gp rangle^3langle g rangle^3\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Plugging this into the previous expression of the higher-order moment langle gp^3 rangle we finally obtain:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\n langle gp^3 rangle = 3 fraclangle gp^2 rangle langle gp ranglelangle g rangle - 2 fraclangle gp rangle^3langle g rangle^2\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/#conditional_derivative_matching","page":"Moment Closure Approximations","title":"Conditional Derivative Matching","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"The conditional derivative matching boils down to approximating the higher order conditional moments in terms of lower order conditional moments using the standard derivative matching:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle p^3 g = 1 rangle = fraclangle p^2 g = 1 rangle^3langle p g= 1 rangle^3\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Using langle gp^j rangle = langle p^j g=1 rangle langle g rangle again, we find:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle gp^3 rangle = fraclangle gp^2 rangle^3 langle g ranglelangle gp rangle^3\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Note that conditional moment closure is fully applicable to systems containing multiple binary species. For example, given two two-state genes, g_1 and g_2, the same Bernoulli variable properties hold and the higher-order moments can be expressed as:","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"beginalign*\nlangle g_1 g_2 p^j rangle = langle p^j g_1 = g_2 = 1 rangle langle g_1 g_2 rangle quad j in 12dotsc\nendalign*","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Such higher order conditional moments can again be closed using normal closure or derivative matching.","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"Note that the description here is taken from [16] and we urge the reader to see the paper for more details. In addition, we have used MomentClosure to apply conditional closures on genetic feedback loops and reproduce some of the published results in this example.","category":"page"},{"location":"theory/moment_closure_approximations/#References","page":"Moment Closure Approximations","title":"References","text":"","category":"section"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[1]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Approximation and inference methods for stochastic biochemical kinetics - a tutorial review\", Journal of Physics A: Mathematical and Theoretical 50, 093001 (2017). https://doi.org/10.1088/1751-8121/aa54d9","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[2]: D. Schnoerr, G. Sanguinetti, and R. Grima, \"Comparison of different moment-closure approximations for stochastic chemical kinetics\", The Journal of Chemical Physics 143, 185101 (2015). https://doi.org/10.1063/1.4934990","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[3]: E. Lakatos, A. Ale, P. D. W. Kirk, and M. P. H. Stumpf, \"Multivariate moment closure techniques for stochastic kinetic models\", The Journal of Chemical Physics 143, 094107 (2015). https://doi.org/10.1063/1.4929837","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[4]: J. Hespanha, \"Moment closure for biochemical networks\", in 2008 3rd International Symposium on Communications, Control and Signal Processing (Mar. 2008), pp. 142–147. https://doi.org/10.1109/ISCCSP.2008.4537208","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[5]: L. A. Goodman, \"Population Growth of the Sexes\", Biometrics 9, Publisher: [Wiley, International Biometric Society], 212–225 (1953). https://doi.org/10.2307/3001852","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[6]: P. Whittle, \"On the use of the normal approximation in the treatment of stochastic processes\", Journal of the Royal Statistical Society: Series B (Methodological) 19, 268–281 (1957). https://doi.org/10.1111/j.2517-6161.1957.tb00263.x","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[7]: N. Balakrishnan, N. L. Johnson, and S. Kotz, “A note on relationships between moments, central moments and cumulants from multivariate distributions”, Statistics & Probability Letters 39, 49–54 (1998). https://doi.org/10.1016/S0167-7152(98)00027-3","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[8]: I. Nasell, \"An extension of the moment closure method\", Theoretical Population Biology 64, 233–239 (2003). https://doi.org/10.1016/S0040-5809(03)00074-1","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[9]: M. J. Keeling, \"Multiplicative Moments and Measures of Persistence in Ecology\", Journal of Theoretical Biology 205, 269–281 (2000). https://doi.org/10.1006/jtbi.2000.2066","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[10]: E. L. Crow and K. Shimizu, eds., Lognormal Distributions: Theory and Applications (Marcel Dekker, 1988).","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[11]: N. L. Johnson, S. Kotz, and N. Balakrishnan, Discrete Multivariate Distributions (Wiley, Feb. 1997).","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[12]: A. M. Mathal and P. G. Moschopoulos, \"A form of multivariate gamma distribution\", Annals of the Institute of Statistical Mathematics 44, 97–106 (1992). https://doi.org/10.1007/BF00048672","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[13]: E. Furman, \"On a multivariate gamma distribution\", Statistics & Probability Letters 78, 2353–2360 (2008). https://doi.org/10.1016/j.spl.2008.02.012","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[14]: A. Singh and J. P. Hespanha, \"Lognormal Moment Closures for Biochemical Reactions\", in Proceedings of the 45th IEEE Conference on Decision and Control, ISSN:0191-2216 (Dec. 2006), pp. 2063–2068. https://doi.org/10.1109/CDC.2006.376994","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[15]: A. Singh and J. P. Hespanha, \"Approximate Moment Dynamics for Chemically Reacting Systems\", IEEE Transactions on Automatic Control 56, 414–418 (2011). https://doi.org/10.1109/TAC.2010.2088631","category":"page"},{"location":"theory/moment_closure_approximations/","page":"Moment Closure Approximations","title":"Moment Closure Approximations","text":"[16]: M. Soltani, C. A. Vargas-Garcia, and A. Singh, \"Conditional Moment Closure Schemes for Studying Stochastic Dynamics of Genetic Circuits\", IEEE Transactions on Biomedical Circuits and Systems 9, 518–526 (2015). https://doi.org/10.1109/tbcas.2015.2453158","category":"page"},{"location":"tutorials/LMA_example/#linear_mapping_approximation_example","page":"LMA Example","title":"LMA Example","text":"","category":"section"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"In this example, we demonstrate how the Linear Mapping Approximation (LMA) can be applied on chemical reaction networks using MomentClosure. We illustrate the approach on models of a simple nonlinear feedback loop and a genetic toggle switch, in turn reproducing some of the results published in the original LMA paper by Cao and Grima [1].","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"We start by considering a simple nonlinear gene regulatory network (GRN) shown in the diagram below (adapted from [1]):","category":"page"},{"location":"tutorials/LMA_example/","page":"LMA Example","title":"LMA Example","text":"
⠀
⠀
