The system size expansion, also known as van Kampen's expansion or the Ω-expansion, is a technique pioneered by Nico van Kampen used in the analysis of stochastic processes. Specifically, it allows one to find an approximation to the solution of a master equation with nonlinear transition rates. The leading order term of the expansion is given by the linear noise approximation, in which the master equation is approximated by a Fokker–Planck equation with linear coefficients determined by the transition rates and stoichiometry of the system. Less formally, it is normally straightforward to write down a mathematical description of a system where processes happen randomly (for example, radioactive atoms randomly decay in a physical system, or genes that are expressed stochastically in a cell). However, these mathematical descriptions are often too difficult to solve for the study of the systems statistics (for example, the mean and variance of the number of atoms or proteins as a function of time). The system size expansion allows one to obtain an approximate statistical description that can be solved much more easily than the master equation.
Preliminaries Systems that admit a treatment with the system size expansion may be described by a probability distribution P ( X , t ) {\displaystyle P(X,t)} , giving the probability of observing the system in state X {\displaystyle X} at time t {\displaystyle t} . X {\displaystyle X} may be, for example, a vector with elements corresponding to the number of molecules of different chemical species in a system. In a system of size Ω {\displaystyle \Omega } (intuitively interpreted as the volume), we will adopt the following nomenclature: X {\displaystyle \mathbf {X} } is a vector of macroscopic copy numbers, x = X / Ω {\displaystyle \mathbf {x} =\mathbf {X} /\Omega } is a vector of concentrations, and ϕ {\displaystyle \mathbf {\phi } } is a vector of deterministic concentrations, as they would appear according to the rate equation in an infinite system. x {\displaystyle \mathbf {x} } and X {\displaystyle \mathbf {X} } are thus quantities subject to stochastic effects. A master equation describes the time evolution of this probability. Henceforth, a system of chemical reactions will be discussed to provide a concrete example, although the nomenclature of "species" and "reactions" is generalisable. A system involving N {\displaystyle N} species and R {\displaystyle R} reactions can be described with the master equation:
∂ P ( X , t ) ∂ t = Ω ∑ j = 1 R ( ∏ i = 1 N E − S i j − 1 ) f j ( x , Ω ) P ( X , t ) . {\displaystyle {\frac {\partial P(\mathbf {X} ,t)}{\partial t}}=\Omega \sum _{j=1}^{R}\left(\prod _{i=1}^{N}\mathbb {E} ^{-S_{ij}}-1\right)f_{j}(\mathbf {x} ,\Omega )P(\mathbf {X} ,t).}
Here, Ω {\displaystyle \Omega } is the system size, E {\displaystyle \mathbb {E} } is an operator which will be addressed later, S i j {\displaystyle S_{ij}} is the stoichiometric matrix for the system (in which element S i j {\displaystyle S_{ij}} gives the stoichiometric coefficient for species i {\displaystyle i} in reaction j {\displaystyle j} ), and f j {\displaystyle f_{j}} is the rate of reaction j {\displaystyle j} given a state x {\displaystyle \mathbf {x} } and system size Ω {\displaystyle \Omega } .
E − S i j {\displaystyle \mathbb {E} ^{-S_{ij}}} is a step operator, removing S i j {\displaystyle S_{ij}} from the i {\displaystyle i} th element of its argument. For example, E − S 23 f ( x 1 , x 2 , x 3 ) = f ( x 1 , x 2 − S 23 , x 3 ) {\displaystyle \mathbb {E} ^{-S_{23}}f(x_{1},x_{2},x_{3})=f(x_{1},x_{2}-S_{23},x_{3})} . This formalism will be useful later. The above equation can be interpreted as follows. The initial sum on the RHS is over all reactions. For each reaction j {\displaystyle j} , the brackets immediately following the sum give two terms. The term with the simple coefficient −1 gives the probability flux away from a given state X {\displaystyle \mathbf {X} } due to reaction j {\displaystyle j} changing the state. The term preceded by the product of step operators gives the probability flux due to reaction j {\displaystyle j} changing a different state X ′ {\displaystyle \mathbf {X'} } into state X {\displaystyle \mathbf {X} } . The product of step operators constructs this state X ′ {\displaystyle \mathbf {X'} } .
Example For example, consider the (linear) chemical system involving two chemical species X 1 {\displaystyle X_{1}} and X 2 {\displaystyle X_{2}} and the reaction X 1 → X 2 {\displaystyle X_{1}\rightarrow X_{2}} . In this system, N = 2 {\displaystyle N=2} (species), R = 1 {\displaystyle R=1} (reactions). A state of the system is a vector X = { n 1 , n 2 } {\displaystyle \mathbf {X} =\{n_{1},n_{2}\}} , where n 1 , n 2 {\displaystyle n_{1},n_{2}} are the number of molecules of X 1 {\displaystyle X_{1}} and X 2 {\displaystyle X_{2}} respectively. Let f 1 ( x , Ω ) = n 1 Ω = x 1 {\displaystyle f_{1}(\mathbf {x} ,\Omega )={\frac {n_{1}}{\Omega }}=x_{1}} , so that the rate of reaction 1 (the only reaction) depends on the concentration of X 1 {\displaystyle X_{1}} . The stoichiometry matrix is ( − 1 , 1 ) T {\displaystyle (-1,1)^{T}} . Then the master equation reads:
∂ P ( X , t ) ∂ t = Ω ( E − S 11 E − S 21 − 1 ) f 1 ( X Ω ) P ( X , t ) = Ω ( f 1 ( X + Δ X Ω ) P ( X + Δ X , t ) − f 1 ( X Ω ) P ( X , t ) ) , {\displaystyle {\begin{aligned}{\frac {\partial P(\mathbf {X} ,t)}{\partial t}}&=\Omega \left(\mathbb {E} ^{-S_{11}}\mathbb {E} ^{-S_{21}}-1\right)f_{1}\left({\frac {\mathbf {X} }{\Omega }}\right)P(\mathbf {X} ,t)\\&=\Omega \left(f_{1}\left({\frac {\mathbf {X} +\mathbf {\Delta X} }{\Omega }}\right)P\left(\mathbf {X} +\mathbf {\Delta X} ,t\right)-f_{1}\left({\frac {\mathbf {X} }{\Omega }}\right)P\left(\mathbf {X} ,t\right)\right),\end{aligned}}}
where Δ X = { 1 , − 1 } {\displaystyle \mathbf {\Delta X} =\{1,-1\}} is the shift caused by the action of the product of step operators, required to change state X {\displaystyle \mathbf {X} } to a precursor state X ′ {\displaystyle \mathbf {X} '} .
Linear noise approximation If the master equation possesses nonlinear transition rates, it may be impossible to solve it analytically. The system size expansion utilises the ansatz that the variance of the steady-state probability distribution of constituent numbers in a population scales like the system size. This ansatz is used to expand the master equation in terms of a small parameter given by the inverse system size. Specifically, let us write the X i {\displaystyle X_{i}} , the copy number of component i {\displaystyle i} , as a sum of its "deterministic" value (a scaled-up concentration) and a random variable ξ {\displaystyle \xi } , scaled by Ω 1 / 2 {\displaystyle \Omega ^{1/2}} :
X i = Ω ϕ i + Ω 1 / 2 ξ i . {\displaystyle X_{i}=\Omega \phi _{i}+\Omega ^{1/2}\xi _{i}.}
The probability distribution of X {\displaystyle \mathbf {X} } can then be rewritten in the vector of random variables ξ {\displaystyle \xi } :
P ( X , t ) = P ( Ω ϕ + Ω 1 / 2 ξ ) = Π ( ξ , t ) . {\displaystyle P(\mathbf {X} ,t)=P(\Omega \mathbf {\phi } +\Omega ^{1/2}\mathbf {\xi } )=\Pi (\mathbf {\xi } ,t).}
Consider how to write reaction rates f {\displaystyle f} and the step operator E {\displaystyle \mathbb {E} } in terms of this new random variable. Taylor expansion of the transition rates gives:
f j ( x ) = f j ( ϕ + Ω − 1 / 2 ξ ) = f j ( ϕ ) + Ω − 1 / 2 ∑ i = 1 N ∂ f j ( ϕ ) ∂ ϕ i ξ i + O ( Ω − 1 ) . {\displaystyle f_{j}(\mathbf {x} )=f_{j}(\mathbf {\phi } +\Omega ^{-1/2}\mathbf {\xi } )=f_{j}(\mathbf {\phi } )+\Omega ^{-1/2}\sum _{i=1}^{N}{\frac {\partial f_{j}(\mathbf {\phi } )}{\partial \phi _{i}}}\xi _{i}+O(\Omega ^{-1}).}
The step operator has the effect E f ( n ) → f ( n + 1 ) {\displaystyle \mathbb {E} f(n)\rightarrow f(n+1)} and hence E f ( ξ ) → f ( ξ + Ω − 1 / 2 ) {\displaystyle \mathbb {E} f(\xi )\rightarrow f(\xi +\Omega ^{-1/2})} :
∏ i = 1 N E − S i j ≃ 1 − Ω − 1 / 2 ∑ i S i j ∂ ∂ ξ i + Ω − 1 2 ∑ i ∑ k S i j S k j ∂ 2 ∂ ξ i ∂ ξ k + O ( Ω − 3 / 2 ) . {\displaystyle \prod _{i=1}^{N}\mathbb {E} ^{-S_{ij}}\simeq 1-\Omega ^{-1/2}\sum _{i}S_{ij}{\frac {\partial }{\partial \xi _{i}}}+{\frac {\Omega ^{-1}}{2}}\sum _{i}\sum _{k}S_{ij}S_{kj}{\frac {\partial ^{2}}{\partial \xi _{i}\,\partial \xi _{k}}}+O(\Omega ^{-3/2}).}
We are now in a position to recast the master equation.
∂ Π ( ξ , t ) ∂ t − Ω 1 / 2 ∑ i = 1 N ∂ ϕ i ∂ t ∂ Π ( ξ , t ) ∂ ξ i = Ω ∑ j = 1 R ( − Ω − 1 / 2 ∑ i S i j ∂ ∂ ξ i + Ω − 1 2 ∑ i ∑ k S i j S k j ∂ 2 ∂ ξ i ∂ ξ k + O ( Ω − 3 / 2 ) )
× ( f j ( ϕ ) + Ω − 1 / 2 ∑ i ∂ f j ( ϕ ) ∂ ϕ i ξ i + O ( Ω − 1 ) ) Π ( ξ , t ) . {\displaystyle {\begin{aligned}&{}\quad {\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial t}}-\Omega ^{1/2}\sum _{i=1}^{N}{\frac {\partial \phi _{i}}{\partial t}}{\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial \xi _{i}}}\\&=\Omega \sum _{j=1}^{R}\left(-\Omega ^{-1/2}\sum _{i}S_{ij}{\frac {\partial }{\partial \xi _{i}}}+{\frac {\Omega ^{-1}}{2}}\sum _{i}\sum _{k}S_{ij}S_{kj}{\frac {\partial ^{2}}{\partial \xi _{i}\,\partial \xi _{k}}}+O(\Omega ^{-3/2})\right)\\&{}\qquad \times \left(f_{j}(\mathbf {\phi } )+\Omega ^{-1/2}\sum _{i}{\frac {\partial f_{j}(\mathbf {\phi } )}{\partial \phi _{i}}}\xi _{i}+O(\Omega ^{-1})\right)\Pi (\mathbf {\xi } ,t).\end{aligned}}}
This rather frightening expression makes a bit more sense when we gather terms in different powers of Ω {\displaystyle \Omega } . First, terms of order Ω 1 / 2 {\displaystyle \Omega ^{1/2}} give
∑ i = 1 N ∂ ϕ i ∂ t ∂ Π ( ξ , t ) ∂ ξ i = ∑ i = 1 N ∑ j = 1 R S i j f j ( ϕ ) ∂ Π ( ξ , t ) ∂ ξ i . {\displaystyle \sum _{i=1}^{N}{\frac {\partial \phi _{i}}{\partial t}}{\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial \xi _{i}}}=\sum _{i=1}^{N}\sum _{j=1}^{R}S_{ij}f_{j}(\mathbf {\phi } ){\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial \xi _{i}}}.}
These terms cancel, due to the macroscopic reaction equation
∂ ϕ i ∂ t = ∑ j = 1 R S i j f j ( ϕ ) . {\displaystyle {\frac {\partial \phi _{i}}{\partial t}}=\sum _{j=1}^{R}S_{ij}f_{j}(\mathbf {\phi } ).}
The terms of order Ω 0 {\displaystyle \Omega ^{0}} are more interesting:
∂ Π ( ξ , t ) ∂ t = ∑ j ( ∑ i k − S i j ∂ f j ∂ ϕ k ∂ ( ξ k Π ( ξ , t ) ) ∂ ξ i + 1 2 f j ∑ i k S i j S k j ∂ 2 Π ( ξ , t ) ∂ ξ i ∂ ξ k ) , {\displaystyle {\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial t}}=\sum _{j}\left(\sum _{ik}-S_{ij}{\frac {\partial f_{j}}{\partial \phi _{k}}}{\frac {\partial (\xi _{k}\Pi (\mathbf {\xi } ,t))}{\partial \xi _{i}}}+{\frac {1}{2}}f_{j}\sum _{ik}S_{ij}S_{kj}{\frac {\partial ^{2}\Pi (\mathbf {\xi } ,t)}{\partial \xi _{i}\,\partial \xi _{k}}}\right),}
which can be written as
∂ Π ( ξ , t ) ∂ t = − ∑ i k A i k ∂ ( ξ k Π ) ∂ ξ i + 1 2 ∑ i k [ B B T ] i k ∂ 2 Π ∂ ξ i ∂ ξ k , {\displaystyle {\frac {\partial \Pi (\mathbf {\xi } ,t)}{\partial t}}=-\sum _{ik}A_{ik}{\frac {\partial (\xi _{k}\Pi )}{\partial \xi _{i}}}+{\frac {1}{2}}\sum _{ik}[\mathbf {BB} ^{T}]_{ik}{\frac {\partial ^{2}\Pi }{\partial \xi _{i}\,\partial