diff --git a/dev/.documenter-siteinfo.json b/dev/.documenter-siteinfo.json index eb15b7b..cac0c56 100644 --- a/dev/.documenter-siteinfo.json +++ b/dev/.documenter-siteinfo.json @@ -1 +1 @@ -{"documenter":{"julia_version":"1.10.0","generation_timestamp":"2024-02-02T10:56:03","documenter_version":"1.2.1"}} \ No newline at end of file +{"documenter":{"julia_version":"1.10.0","generation_timestamp":"2024-02-02T11:15:36","documenter_version":"1.2.1"}} \ No newline at end of file diff --git a/dev/predator_prey_tutorial/17593e40.svg b/dev/PredatorPreyTutorial/68aa82f6.svg similarity index 76% rename from dev/predator_prey_tutorial/17593e40.svg rename to dev/PredatorPreyTutorial/68aa82f6.svg index 128471e..d94635c 100644 --- a/dev/predator_prey_tutorial/17593e40.svg +++ b/dev/PredatorPreyTutorial/68aa82f6.svg @@ -1,176 +1,176 @@ - + - + - + - + - + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/dev/predator_prey_tutorial/index.html b/dev/PredatorPreyTutorial/index.html similarity index 55% rename from dev/predator_prey_tutorial/index.html rename to dev/PredatorPreyTutorial/index.html index dfe5c7c..906abde 100644 --- a/dev/predator_prey_tutorial/index.html +++ b/dev/PredatorPreyTutorial/index.html @@ -1,5 +1,5 @@ -Tutorial · LinearNoiseApproximation.jl

Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model

This example demonstrates how to use LinearNoiseApproximation.jl package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic trajectories obtained using the Gillespie algorithm.

using LinearNoiseApproximation
+Tutorial · LinearNoiseApproximation.jl
Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model
This example demonstrates how to use LinearNoiseApproximation.jl package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic trajectories obtained using the Gillespie algorithm.
using LinearNoiseApproximation
 using DifferentialEquations
 using Catalyst
 using JumpProcesses
@@ -78,4 +78,4 @@
         xlim=(0, duration + 1),
     )
 end
-plt
Example block output
save the plot Plots.savefig(plt, "predatorpreyfig.png")
This page was generated using Literate.jl.

+plt
Example block output

save the plot Plots.savefig(plt, "predatorpreyfig.png")

This page was generated using Literate.jl.

diff --git a/dev/api/index.html b/dev/api/index.html index 3c277e3..dff50f0 100644 --- a/dev/api/index.html +++ b/dev/api/index.html @@ -1,2 +1,2 @@ -API · LinearNoiseApproximation.jl

Main API

Types

LinearNoiseApproximation.LNASystem — Type

A type that contains the information of the LNA system.

Fields

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • odesys::ODESystem: the linear noise approximation system, it is a ODESystem that extends the ReactionSystem by adding the covariance terms
  • kwargs::Any: Optional keyword arguments
source

Functions

LinearNoiseApproximation.get_ode_propensity — Function

Function to obtain the ODEs according to the rate equations

Arguments

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • combinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false
source
LinearNoiseApproximation.get_symbolic_matrix — Function

Function to obtain the symbolic matrix of the LNA system

Arguments

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • combinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false
source
+API · LinearNoiseApproximation.jl

Main API

Types

LinearNoiseApproximation.LNASystem — Type

A type that contains the information of the LNA system.

Fields

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • odesys::ODESystem: the linear noise approximation system, it is a ODESystem that extends the ReactionSystem by adding the covariance terms
  • kwargs::Any: Optional keyword arguments
source

Functions

LinearNoiseApproximation.get_ode_propensity — Function

Function to obtain the ODEs according to the rate equations

Arguments

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • combinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false
source
LinearNoiseApproximation.get_symbolic_matrix — Function

Function to obtain the symbolic matrix of the LNA system

Arguments

  • rn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl
  • combinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false
source
diff --git a/dev/index.html b/dev/index.html index cdc6d5f..ac2d38b 100644 --- a/dev/index.html +++ b/dev/index.html @@ -1,5 +1,5 @@ -Home · LinearNoiseApproximation.jl

LinearNoiseApproximation

Installation

using Pkg
+Home · LinearNoiseApproximation.jl
LinearNoiseApproximation
Installation
using Pkg
 Pkg.add("LinearNoiseApproximation")
 using LinearNoiseApproximation
This package provides a numerical method of applying linear noise approximation (LNA) to a given reaction system using Catalyst. The derived expanded system can be solved using DifferentialEquations.jl.
Usage
Reaction network with system size Ω
Suppose we have a reaction network
rn = @reaction_network TwoStage begin
     @parameters  v0 v1 d0 d1 Ω
@@ -23,4 +23,4 @@
     \end{aligned}
 \end{equation}\]
where $P(\boldsymbol{n}, t)$ is the probability of the system being in state $\boldsymbol{n}$ at time $t$, 
\[f_r(\boldsymbol{n}, t) = k_r\Omega \prod_{i=1}^{N} \frac{n_i!}{(n_i-s_{ir})!\Omega^{s_{ir}}} \]
is the propensity function of reaction $r$ at state $\boldsymbol{n}$, and $\mathbf{S}_r$ is the stoichiometric vector of reaction $r$, $\Omega$ is the system size (or volume of the system).
The chemical master equation is written directly from the rate constants and stoichiometries of all the elementary reactions of a chemical system, but neither analytical nor numerical solutions are in general available. Fortunately, the chemical master equation can often be simplified in a linear noise approximation. Linear noise approximation is an expansion of the CME taking the inverse system size $1/\Omega$ as the perturbed variable, which is originally developed by [3]. The idea is to separate concentrations into a deterministic part, given by the solution of the deterministic rate equations, and a part describing the fluctuations about the deterministic mean
\[\frac{n_i}{\Omega} = \phi_i  + \Omega^{-1/2} \epsilon_i,\]
where $\phi_i$ is the solution of the deterministic rate equations (2), and $\epsilon_i$ represents fluctuations about the deterministic mean. Define $\boldsymbol{\Sigma}=\left(\epsilon_{ij}\right)_{N\times N}$ to be the covariance matrix of the fluctuations, the linear noise approximation is given by:
\[\begin{equation}
     \partial_t \boldsymbol{\Sigma} = \mathbf{A} \boldsymbol{\Sigma} + \boldsymbol{\Sigma} \mathbf{A}^T +  \Omega^{-1} \mathbf{B}, \tag{3}
-\end{equation}\]
where $\mathbf{A}=\mathbf{A}(\boldsymbol{\phi}),\, \mathbf{B}=\mathbf{B}(\boldsymbol{\phi})$ both depend on the solution $\boldsymbol{\phi}$ of the rate equation \eqref{eq:rate_equations}, which are defined by 
\[ \mathbf{A} = \left(A_{ij}\right)_{N\times N}, A_{ij}=\sum_{r=1}^{R} S_{ir}\partial_{\phi_j}g_r(\boldsymbol{\phi}),\]
as the Jacobian matrix of the deterministic rate equations, and 
\[\mathbf{B} =\left(B_{ij}\right)_{N\times N}, B_{ij} = \sum_{r=1}^{R} S_{ir}S_{jr}g_r(\boldsymbol{\phi}).\]
In this formulation, the LNA allows for analytical solutions that are locally valid close to macroscopic trajectories (solution of the rate equations) of the system. We refer to the review paper [2] for more details on the CME and LNA.
References
[1] Nicolaas Godfried van Kampen. The Expansion of the Master Equation. John Wiley & Sons, Inc., 2007. 
[2] David Schnoerr, Guido Sanguinetti, and Ramon Grima. Approximation and inference methods for stochastic biochemical kinetics — a tutorial review (2017). Journal of Physics A: Mathematical and Theoretical, 50(9):093001.
[3] Nicolaas Godfried van Kampen. Stochastic Processes in Physics and Chemistry, volume 1. Elsevier, 1992.

+\end{equation}\]
where $\mathbf{A}=\mathbf{A}(\boldsymbol{\phi}),\, \mathbf{B}=\mathbf{B}(\boldsymbol{\phi})$ both depend on the solution $\boldsymbol{\phi}$ of the rate equation \eqref{eq:rate_equations}, which are defined by 
\[ \mathbf{A} = \left(A_{ij}\right)_{N\times N}, A_{ij}=\sum_{r=1}^{R} S_{ir}\partial_{\phi_j}g_r(\boldsymbol{\phi}),\]
as the Jacobian matrix of the deterministic rate equations, and 
\[\mathbf{B} =\left(B_{ij}\right)_{N\times N}, B_{ij} = \sum_{r=1}^{R} S_{ir}S_{jr}g_r(\boldsymbol{\phi}).\]
In this formulation, the LNA allows for analytical solutions that are locally valid close to macroscopic trajectories (solution of the rate equations) of the system. We refer to the review paper [2] for more details on the CME and LNA.
References
[1] Nicolaas Godfried van Kampen. The Expansion of the Master Equation. John Wiley & Sons, Inc., 2007. 
[2] David Schnoerr, Guido Sanguinetti, and Ramon Grima. Approximation and inference methods for stochastic biochemical kinetics — a tutorial review (2017). Journal of Physics A: Mathematical and Theoretical, 50(9):093001.
[3] Nicolaas Godfried van Kampen. Stochastic Processes in Physics and Chemistry, volume 1. Elsevier, 1992.
diff --git a/dev/search_index.js b/dev/search_index.js index f7d503f..2dcc267 100644 --- a/dev/search_index.js +++ b/dev/search_index.js @@ -1,3 +1,3 @@ var documenterSearchIndex = {"docs": -[{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"EditURL = \"../../examples/predator_prey_tutorial.jl\"","category":"page"},{"location":"predator_prey_tutorial/#Tutorial:-Using-LinearNoiseApproximation.jl-to-solve-a-Lotka-Volterra-model","page":"Tutorial","title":"Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model","text":"","category":"section"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"This example demonstrates how to use LinearNoiseApproximation.jl package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic trajectories obtained using the Gillespie algorithm.","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"using LinearNoiseApproximation\nusing DifferentialEquations\nusing Catalyst\nusing JumpProcesses\nusing Plots","category":"page"},{"location":"predator_prey_tutorial/#Define-the-Lotka-Volterra-model","page":"Tutorial","title":"Define the Lotka-Volterra model","text":"","category":"section"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population U grows at a rate represented by the parameter alpha in the absence of predators, but decreases as they are consumed by the predator population V at a rate determined by the interaction strength parameter, beta. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, delta. The corresponding ODE system reads","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"beginequation\n\tbeginaligned\n\t\tfracmathrmdUmathrmd t = alpha U - beta U V\n\t\tfracmathrmdVmathrmd t = beta U V - delta V\n\tendalignedquad tin (0 t_textend)\nendequation","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"rn = @reaction_network begin\n @parameters α β δ\n @species U(t) V(t)\n α, U --> 2*U\n β, U + V --> 2*V\n δ, V --> 0\nend","category":"page"},{"location":"predator_prey_tutorial/#Convert-the-model-to-a-JumpSystem-for-the-Gillespie-algorithm","page":"Tutorial","title":"Convert the model to a JumpSystem for the Gillespie algorithm","text":"","category":"section"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"Using Catalyst.jl we can convert ReactionSystem to a JumpSystem which can be used to simulate the stochastic trajectories using the Gillespie algorithm.","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"jumpsys = convert(JumpSystem, rn)\n\n#define the initial conditions, parameters, and time span\nu0 = [120.0, 140.0]\nps = [0.8, 0.005, 0.4]\nduration = 30.0\ntspan = (0.0, duration)\n\n#solve the model using the Gillespie algorithm\ndprob = DiscreteProblem(jumpsys, u0, tspan, ps)\njprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))\njsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)\n\n# solve the model using the Linear Noise Approximation (LNA)\nlna = LNASystem(rn)\nalg = Vern7()\nprob = ODEProblem(lna, u0, tspan, ps)\nsol = solve(prob, alg, saveat=0.5, seed=2)\n\n#find the indices of the mean and variance states\nmean_idxs, var_idxs = find_states_cov_number([1, 2], lna)\n\n\n#plot the results\nplt = Plots.plot(; size=(800, 350))\ncolor_palette=[\"#1f78b4\" \"#ff7f00\"]\nlabel = [\"U\" \"V\"]\n\nfor (i, (mean_idx, var_idx)) in enumerate(zip(mean_idxs, var_idxs))\n mean = sol[mean_idx, :]\n var = sol[var_idx, :]\n\n Plots.plot!(\n jsol.t,\n jsol[mean_idx, :];\n label=string(\"SSA: \", label[i]),\n xlabel=\"Time\",\n ylabel=\"Numbers\",\n color=color_palette[i],\n legend=:topleft,\n st=:scatter,\n alpha=0.8,\n left_margin=5Plots.mm,\n bottom_margin=5Plots.mm,\n markerstrokewidth=0,\n xlim=(0, duration + 1),\n xticks=0:5:duration + 1,\n )\n Plots.plot!(\n sol.t,\n mean;\n label=string(\"LNA: \", label[i]),\n ribbon=sqrt.(var),\n ribbonalpha=0.3,\n color=color_palette[i],\n lw=5,\n xlim=(0, duration + 1),\n )\nend\nplt","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"save the plot Plots.savefig(plt, \"predatorpreyfig.png\")","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"","category":"page"},{"location":"predator_prey_tutorial/","page":"Tutorial","title":"Tutorial","text":"This page was generated using Literate.jl.","category":"page"},{"location":"api/#Main-API","page":"API","title":"Main API","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"CurrentModule = LinearNoiseApproximation","category":"page"},{"location":"api/#Types","page":"API","title":"Types","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"LNASystem","category":"page"},{"location":"api/#LinearNoiseApproximation.LNASystem","page":"API","title":"LinearNoiseApproximation.LNASystem","text":"A type that contains the information of the LNA system.\n\nFields\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\nodesys::ODESystem: the linear noise approximation system, it is a ODESystem that extends the ReactionSystem by adding the covariance terms\nkwargs::Any: Optional keyword arguments\n\n\n\n\n\n","category":"type"},{"location":"api/#Functions","page":"API","title":"Functions","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"get_ode_propensity\nget_symbolic_matrix\nfind_states_cov_number\nexpand_initial_conditions","category":"page"},{"location":"api/#LinearNoiseApproximation.get_ode_propensity","page":"API","title":"LinearNoiseApproximation.get_ode_propensity","text":"Function to obtain the ODEs according to the rate equations\n\nArguments\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\ncombinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.get_symbolic_matrix","page":"API","title":"LinearNoiseApproximation.get_symbolic_matrix","text":"Function to obtain the symbolic matrix of the LNA system\n\nArguments\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\ncombinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.find_states_cov_number","page":"API","title":"LinearNoiseApproximation.find_states_cov_number","text":"Function to find the index of the covariance matrix in the LNA system.\n\nArguments\n\nnum::Int: the index of the species in the reaction network\nlna_sys::LNASystem: the LNA system\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.expand_initial_conditions","page":"API","title":"LinearNoiseApproximation.expand_initial_conditions","text":"Function to expand the initial conditions of the LNA system.\n\nArguments\n\nsys::LNASystem: the LNA system\nu0::Vector: the initial conditions of the orignal system (rate equations)\n\n\n\n\n\n","category":"function"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = LinearNoiseApproximation","category":"page"},{"location":"#LinearNoiseApproximation","page":"Home","title":"LinearNoiseApproximation","text":"","category":"section"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"LinearNoiseApproximation\")\nusing LinearNoiseApproximation","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package provides a numerical method of applying linear noise approximation (LNA) to a given reaction system using Catalyst. The derived expanded system can be solved using DifferentialEquations.jl.","category":"page"},{"location":"#Usage","page":"Home","title":"Usage","text":"","category":"section"},{"location":"#Reaction-network-with-system-size-Ω","page":"Home","title":"Reaction network with system size Ω","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Suppose we have a reaction network","category":"page"},{"location":"","page":"Home","title":"Home","text":"rn = @reaction_network TwoStage begin\n @parameters v0 v1 d0 d1 Ω\n @species M(t) P(t)\n v0*Ω, ∅ --> M\n v1, M --> M + P\n d0, M --> ∅\n d1, P --> ∅\nend","category":"page"},{"location":"","page":"Home","title":"Home","text":"Here the system size Ω is integrated in the system (propensities), you can simply obtain the LNA system by calling","category":"page"},{"location":"","page":"Home","title":"Home","text":"LNAsys = LNASystem(rn)","category":"page"},{"location":"","page":"Home","title":"Home","text":"then solving the LNA by","category":"page"},{"location":"","page":"Home","title":"Home","text":"@variables v0, v1, d0, d1, Ω\nrates = [v0=>4.0,v1=>10.0,d0=>1.0,d1=>1.0,Ω=>5.0]\ntspan = (0.0, 20.0)\nu0 =[0, 0] # initial condition for M and P\nprob = ODEProblem(LNAsys, u0, tspan, rates)\nsol = solve(prob,Vern7(),abstol=1e-7, saveat =1.0)","category":"page"},{"location":"#Brief-introduction-to-linear-noise-approximation","page":"Home","title":"Brief introduction to linear noise approximation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Linear noise approximation (LNA) is a moment-based approximation method for stochastic chemical kinetics [1]. A general chemical kinetics reaction can be described as follows: given a set of chemical species X_i i = 1 ldots N, we define R chemical reactions by the notation","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n k_rsum_i=1^Ns_irX_i rightarrow sum_i=1^Ns_irX_i tag1\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where the stoichiometric coefficients s_ir and s_ir are non-negative integer numbers denoting the numbers of reactant and product molecules, respectively. The quantity k_r in Equation (1) is called the reaction rate constant of the r-th reaction. Classically, the dynamics of a chemical reaction system as in Equation (1) is modelled by the law of mass action. The law of mass action states that the rate of a reaction is proportional to the product of the concentrations of reactant molecules, which lead to the following rate equations as:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n fracmathrmdmathrmdt phi_i = sum_r=1^R S_ir g_r(boldsymbolphi) tag2\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where phi_i is the concentration of species X_i, ","category":"page"},{"location":"","page":"Home","title":"Home","text":"boldsymbolS = S_ir_Ntimes R S_ir=s_ir - s_ir i=1ldotsN r=1ldotsR notag","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the stoichiometric matrix, and ","category":"page"},{"location":"","page":"Home","title":"Home","text":"g_r(boldsymbolphi) = k_r prod_i=1^N phi_i^s_ir ","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the rate of the r-th reaction. ","category":"page"},{"location":"","page":"Home","title":"Home","text":"However, the law of mass action is only valid when the number of molecules is large. When the number of molecules is small, System (1) can instead be modelled by a continuous-time Markov jump process to study the probability of the system being in a particular state at a given time. The dynamics of such a system can be described by the chemical master equation (CME) [2]:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n beginaligned\n fracmathrmdmathrmdt P(boldsymboln t) = sum_r=1^R f_r(boldsymboln - mathbfS_r t) P(boldsymboln - mathbfS_r t) - sum_r=1^R f_r(boldsymboln t) P(boldsymboln t) notag\n endaligned\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where P(boldsymboln t) is the probability of the system being in state boldsymboln at time t, ","category":"page"},{"location":"","page":"Home","title":"Home","text":"f_r(boldsymboln t) = k_rOmega prod_i=1^N fracn_i(n_i-s_ir)Omega^s_ir ","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the propensity function of reaction r at state boldsymboln, and mathbfS_r is the stoichiometric vector of reaction r, Omega is the system size (or volume of the system).","category":"page"},{"location":"","page":"Home","title":"Home","text":"The chemical master equation is written directly from the rate constants and stoichiometries of all the elementary reactions of a chemical system, but neither analytical nor numerical solutions are in general available. Fortunately, the chemical master equation can often be simplified in a linear noise approximation. Linear noise approximation is an expansion of the CME taking the inverse system size 1Omega as the perturbed variable, which is originally developed by [3]. The idea is to separate concentrations into a deterministic part, given by the solution of the deterministic rate equations, and a part describing the fluctuations about the deterministic mean","category":"page"},{"location":"","page":"Home","title":"Home","text":"fracn_iOmega = phi_i + Omega^-12 epsilon_i","category":"page"},{"location":"","page":"Home","title":"Home","text":"where phi_i is the solution of the deterministic rate equations (2), and epsilon_i represents fluctuations about the deterministic mean. Define boldsymbolSigma=left(epsilon_ijright)_Ntimes N to be the covariance matrix of the fluctuations, the linear noise approximation is given by:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n partial_t boldsymbolSigma = mathbfA boldsymbolSigma + boldsymbolSigma mathbfA^T + Omega^-1 mathbfB tag3\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where mathbfA=mathbfA(boldsymbolphi) mathbfB=mathbfB(boldsymbolphi) both depend on the solution boldsymbolphi of the rate equation \\eqref{eq:rate_equations}, which are defined by ","category":"page"},{"location":"","page":"Home","title":"Home","text":" mathbfA = left(A_ijright)_Ntimes N A_ij=sum_r=1^R S_irpartial_phi_jg_r(boldsymbolphi)","category":"page"},{"location":"","page":"Home","title":"Home","text":"as the Jacobian matrix of the deterministic rate equations, and ","category":"page"},{"location":"","page":"Home","title":"Home","text":"mathbfB =left(B_ijright)_Ntimes N B_ij = sum_r=1^R S_irS_jrg_r(boldsymbolphi)","category":"page"},{"location":"","page":"Home","title":"Home","text":"In this formulation, the LNA allows for analytical solutions that are locally valid close to macroscopic trajectories (solution of the rate equations) of the system. We refer to the review paper [2] for more details on the CME and LNA.","category":"page"},{"location":"#References","page":"Home","title":"References","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"[1] Nicolaas Godfried van Kampen. The Expansion of the Master Equation. John Wiley & Sons, Inc., 2007. ","category":"page"},{"location":"","page":"Home","title":"Home","text":"[2] David Schnoerr, Guido Sanguinetti, and Ramon Grima. Approximation and inference methods for stochastic biochemical kinetics — a tutorial review (2017). Journal of Physics A: Mathematical and Theoretical, 50(9):093001.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[3] Nicolaas Godfried van Kampen. Stochastic Processes in Physics and Chemistry, volume 1. Elsevier, 1992.","category":"page"}] +[{"location":"api/#Main-API","page":"API","title":"Main API","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"CurrentModule = LinearNoiseApproximation","category":"page"},{"location":"api/#Types","page":"API","title":"Types","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"LNASystem","category":"page"},{"location":"api/#LinearNoiseApproximation.LNASystem","page":"API","title":"LinearNoiseApproximation.LNASystem","text":"A type that contains the information of the LNA system.\n\nFields\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\nodesys::ODESystem: the linear noise approximation system, it is a ODESystem that extends the ReactionSystem by adding the covariance terms\nkwargs::Any: Optional keyword arguments\n\n\n\n\n\n","category":"type"},{"location":"api/#Functions","page":"API","title":"Functions","text":"","category":"section"},{"location":"api/","page":"API","title":"API","text":"get_ode_propensity\nget_symbolic_matrix\nfind_states_cov_number\nexpand_initial_conditions","category":"page"},{"location":"api/#LinearNoiseApproximation.get_ode_propensity","page":"API","title":"LinearNoiseApproximation.get_ode_propensity","text":"Function to obtain the ODEs according to the rate equations\n\nArguments\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\ncombinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.get_symbolic_matrix","page":"API","title":"LinearNoiseApproximation.get_symbolic_matrix","text":"Function to obtain the symbolic matrix of the LNA system\n\nArguments\n\nrn::ReactionSystem: a reaction network that is used to construct the LNA system, which is from Catalyst.jl\ncombinatoric_ratelaws::Bool: (Optional) whether to use the combinatoric rate laws, default is false\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.find_states_cov_number","page":"API","title":"LinearNoiseApproximation.find_states_cov_number","text":"Function to find the index of the covariance matrix in the LNA system.\n\nArguments\n\nnum::Int: the index of the species in the reaction network\nlna_sys::LNASystem: the LNA system\n\n\n\n\n\n","category":"function"},{"location":"api/#LinearNoiseApproximation.expand_initial_conditions","page":"API","title":"LinearNoiseApproximation.expand_initial_conditions","text":"Function to expand the initial conditions of the LNA system.\n\nArguments\n\nsys::LNASystem: the LNA system\nu0::Vector: the initial conditions of the orignal system (rate equations)\n\n\n\n\n\n","category":"function"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"EditURL = \"../../examples/PredatorPreyTutorial.jl\"","category":"page"},{"location":"PredatorPreyTutorial/#Tutorial:-Using-LinearNoiseApproximation.jl-to-solve-a-Lotka-Volterra-model","page":"Tutorial","title":"Tutorial: Using LinearNoiseApproximation.jl to solve a Lotka-Volterra model","text":"","category":"section"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"This example demonstrates how to use LinearNoiseApproximation.jl package to solve a Lotka-Volterra model using the Linear Noise Approximation (LNA) and compare it to the stochastic trajectories obtained using the Gillespie algorithm.","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"using LinearNoiseApproximation\nusing DifferentialEquations\nusing Catalyst\nusing JumpProcesses\nusing Plots","category":"page"},{"location":"PredatorPreyTutorial/#Define-the-Lotka-Volterra-model","page":"Tutorial","title":"Define the Lotka-Volterra model","text":"","category":"section"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"The Lotka-Volterra model describes the dynamics of predator-prey interactions in an ecosystem. It assumes that the prey population U grows at a rate represented by the parameter alpha in the absence of predators, but decreases as they are consumed by the predator population V at a rate determined by the interaction strength parameter, beta. The predator population decreases in size if they cannot find enough prey to consume, which is represented by the mortality rate parameter, delta. The corresponding ODE system reads","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"beginequation\n\tbeginaligned\n\t\tfracmathrmdUmathrmd t = alpha U - beta U V\n\t\tfracmathrmdVmathrmd t = beta U V - delta V\n\tendalignedquad tin (0 t_textend)\nendequation","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"rn = @reaction_network begin\n @parameters α β δ\n @species U(t) V(t)\n α, U --> 2*U\n β, U + V --> 2*V\n δ, V --> 0\nend","category":"page"},{"location":"PredatorPreyTutorial/#Convert-the-model-to-a-JumpSystem-for-the-Gillespie-algorithm","page":"Tutorial","title":"Convert the model to a JumpSystem for the Gillespie algorithm","text":"","category":"section"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"Using Catalyst.jl we can convert ReactionSystem to a JumpSystem which can be used to simulate the stochastic trajectories using the Gillespie algorithm.","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"jumpsys = convert(JumpSystem, rn)\n\n#define the initial conditions, parameters, and time span\nu0 = [120.0, 140.0]\nps = [0.8, 0.005, 0.4]\nduration = 30.0\ntspan = (0.0, duration)\n\n#solve the model using the Gillespie algorithm\ndprob = DiscreteProblem(jumpsys, u0, tspan, ps)\njprob = JumpProblem(jumpsys, dprob, Direct(), save_positions=(false, false))\njsol = solve(jprob, SSAStepper(), saveat=0.5, seed=2)\n\n# solve the model using the Linear Noise Approximation (LNA)\nlna = LNASystem(rn)\nalg = Vern7()\nprob = ODEProblem(lna, u0, tspan, ps)\nsol = solve(prob, alg, saveat=0.5, seed=2)\n\n#find the indices of the mean and variance states\nmean_idxs, var_idxs = find_states_cov_number([1, 2], lna)\n\n\n#plot the results\nplt = Plots.plot(; size=(800, 350))\ncolor_palette=[\"#1f78b4\" \"#ff7f00\"]\nlabel = [\"U\" \"V\"]\n\nfor (i, (mean_idx, var_idx)) in enumerate(zip(mean_idxs, var_idxs))\n mean = sol[mean_idx, :]\n var = sol[var_idx, :]\n\n Plots.plot!(\n jsol.t,\n jsol[mean_idx, :];\n label=string(\"SSA: \", label[i]),\n xlabel=\"Time\",\n ylabel=\"Numbers\",\n color=color_palette[i],\n legend=:topleft,\n st=:scatter,\n alpha=0.8,\n left_margin=5Plots.mm,\n bottom_margin=5Plots.mm,\n markerstrokewidth=0,\n xlim=(0, duration + 1),\n xticks=0:5:duration + 1,\n )\n Plots.plot!(\n sol.t,\n mean;\n label=string(\"LNA: \", label[i]),\n ribbon=sqrt.(var),\n ribbonalpha=0.3,\n color=color_palette[i],\n lw=5,\n xlim=(0, duration + 1),\n )\nend\nplt","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"save the plot Plots.savefig(plt, \"predatorpreyfig.png\")","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"","category":"page"},{"location":"PredatorPreyTutorial/","page":"Tutorial","title":"Tutorial","text":"This page was generated using Literate.jl.","category":"page"},{"location":"","page":"Home","title":"Home","text":"CurrentModule = LinearNoiseApproximation","category":"page"},{"location":"#LinearNoiseApproximation","page":"Home","title":"LinearNoiseApproximation","text":"","category":"section"},{"location":"#Installation","page":"Home","title":"Installation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"using Pkg\nPkg.add(\"LinearNoiseApproximation\")\nusing LinearNoiseApproximation","category":"page"},{"location":"","page":"Home","title":"Home","text":"This package provides a numerical method of applying linear noise approximation (LNA) to a given reaction system using Catalyst. The derived expanded system can be solved using DifferentialEquations.jl.","category":"page"},{"location":"#Usage","page":"Home","title":"Usage","text":"","category":"section"},{"location":"#Reaction-network-with-system-size-Ω","page":"Home","title":"Reaction network with system size Ω","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Suppose we have a reaction network","category":"page"},{"location":"","page":"Home","title":"Home","text":"rn = @reaction_network TwoStage begin\n @parameters v0 v1 d0 d1 Ω\n @species M(t) P(t)\n v0*Ω, ∅ --> M\n v1, M --> M + P\n d0, M --> ∅\n d1, P --> ∅\nend","category":"page"},{"location":"","page":"Home","title":"Home","text":"Here the system size Ω is integrated in the system (propensities), you can simply obtain the LNA system by calling","category":"page"},{"location":"","page":"Home","title":"Home","text":"LNAsys = LNASystem(rn)","category":"page"},{"location":"","page":"Home","title":"Home","text":"then solving the LNA by","category":"page"},{"location":"","page":"Home","title":"Home","text":"@variables v0, v1, d0, d1, Ω\nrates = [v0=>4.0,v1=>10.0,d0=>1.0,d1=>1.0,Ω=>5.0]\ntspan = (0.0, 20.0)\nu0 =[0, 0] # initial condition for M and P\nprob = ODEProblem(LNAsys, u0, tspan, rates)\nsol = solve(prob,Vern7(),abstol=1e-7, saveat =1.0)","category":"page"},{"location":"#Brief-introduction-to-linear-noise-approximation","page":"Home","title":"Brief introduction to linear noise approximation","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"Linear noise approximation (LNA) is a moment-based approximation method for stochastic chemical kinetics [1]. A general chemical kinetics reaction can be described as follows: given a set of chemical species X_i i = 1 ldots N, we define R chemical reactions by the notation","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n k_rsum_i=1^Ns_irX_i rightarrow sum_i=1^Ns_irX_i tag1\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where the stoichiometric coefficients s_ir and s_ir are non-negative integer numbers denoting the numbers of reactant and product molecules, respectively. The quantity k_r in Equation (1) is called the reaction rate constant of the r-th reaction. Classically, the dynamics of a chemical reaction system as in Equation (1) is modelled by the law of mass action. The law of mass action states that the rate of a reaction is proportional to the product of the concentrations of reactant molecules, which lead to the following rate equations as:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n fracmathrmdmathrmdt phi_i = sum_r=1^R S_ir g_r(boldsymbolphi) tag2\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where phi_i is the concentration of species X_i, ","category":"page"},{"location":"","page":"Home","title":"Home","text":"boldsymbolS = S_ir_Ntimes R S_ir=s_ir - s_ir i=1ldotsN r=1ldotsR notag","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the stoichiometric matrix, and ","category":"page"},{"location":"","page":"Home","title":"Home","text":"g_r(boldsymbolphi) = k_r prod_i=1^N phi_i^s_ir ","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the rate of the r-th reaction. ","category":"page"},{"location":"","page":"Home","title":"Home","text":"However, the law of mass action is only valid when the number of molecules is large. When the number of molecules is small, System (1) can instead be modelled by a continuous-time Markov jump process to study the probability of the system being in a particular state at a given time. The dynamics of such a system can be described by the chemical master equation (CME) [2]:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n beginaligned\n fracmathrmdmathrmdt P(boldsymboln t) = sum_r=1^R f_r(boldsymboln - mathbfS_r t) P(boldsymboln - mathbfS_r t) - sum_r=1^R f_r(boldsymboln t) P(boldsymboln t) notag\n endaligned\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where P(boldsymboln t) is the probability of the system being in state boldsymboln at time t, ","category":"page"},{"location":"","page":"Home","title":"Home","text":"f_r(boldsymboln t) = k_rOmega prod_i=1^N fracn_i(n_i-s_ir)Omega^s_ir ","category":"page"},{"location":"","page":"Home","title":"Home","text":"is the propensity function of reaction r at state boldsymboln, and mathbfS_r is the stoichiometric vector of reaction r, Omega is the system size (or volume of the system).","category":"page"},{"location":"","page":"Home","title":"Home","text":"The chemical master equation is written directly from the rate constants and stoichiometries of all the elementary reactions of a chemical system, but neither analytical nor numerical solutions are in general available. Fortunately, the chemical master equation can often be simplified in a linear noise approximation. Linear noise approximation is an expansion of the CME taking the inverse system size 1Omega as the perturbed variable, which is originally developed by [3]. The idea is to separate concentrations into a deterministic part, given by the solution of the deterministic rate equations, and a part describing the fluctuations about the deterministic mean","category":"page"},{"location":"","page":"Home","title":"Home","text":"fracn_iOmega = phi_i + Omega^-12 epsilon_i","category":"page"},{"location":"","page":"Home","title":"Home","text":"where phi_i is the solution of the deterministic rate equations (2), and epsilon_i represents fluctuations about the deterministic mean. Define boldsymbolSigma=left(epsilon_ijright)_Ntimes N to be the covariance matrix of the fluctuations, the linear noise approximation is given by:","category":"page"},{"location":"","page":"Home","title":"Home","text":"beginequation\n partial_t boldsymbolSigma = mathbfA boldsymbolSigma + boldsymbolSigma mathbfA^T + Omega^-1 mathbfB tag3\nendequation","category":"page"},{"location":"","page":"Home","title":"Home","text":"where mathbfA=mathbfA(boldsymbolphi) mathbfB=mathbfB(boldsymbolphi) both depend on the solution boldsymbolphi of the rate equation \\eqref{eq:rate_equations}, which are defined by ","category":"page"},{"location":"","page":"Home","title":"Home","text":" mathbfA = left(A_ijright)_Ntimes N A_ij=sum_r=1^R S_irpartial_phi_jg_r(boldsymbolphi)","category":"page"},{"location":"","page":"Home","title":"Home","text":"as the Jacobian matrix of the deterministic rate equations, and ","category":"page"},{"location":"","page":"Home","title":"Home","text":"mathbfB =left(B_ijright)_Ntimes N B_ij = sum_r=1^R S_irS_jrg_r(boldsymbolphi)","category":"page"},{"location":"","page":"Home","title":"Home","text":"In this formulation, the LNA allows for analytical solutions that are locally valid close to macroscopic trajectories (solution of the rate equations) of the system. We refer to the review paper [2] for more details on the CME and LNA.","category":"page"},{"location":"#References","page":"Home","title":"References","text":"","category":"section"},{"location":"","page":"Home","title":"Home","text":"[1] Nicolaas Godfried van Kampen. The Expansion of the Master Equation. John Wiley & Sons, Inc., 2007. ","category":"page"},{"location":"","page":"Home","title":"Home","text":"[2] David Schnoerr, Guido Sanguinetti, and Ramon Grima. Approximation and inference methods for stochastic biochemical kinetics — a tutorial review (2017). Journal of Physics A: Mathematical and Theoretical, 50(9):093001.","category":"page"},{"location":"","page":"Home","title":"Home","text":"[3] Nicolaas Godfried van Kampen. Stochastic Processes in Physics and Chemistry, volume 1. Elsevier, 1992.","category":"page"}] }