The Azimuth Project
Experiments in multiple equilibrium states

We have so far always first considered the rate equation and then talked about the master equation and given heuristic arguments on how they are different. I wanted to start from the master equation and show how the rate equation arises from this. That’s exactly what we are going to do in this post.

The master equation

We can also arrive at a master equation. Now say we have a probability,

ψ i 1...i k\psi_{i_1...i_k} of having i 1i_1 things in state 1, i ki_k things in state kk etc. We write

ψ=ψ i 1...i kz 1 i 1z k i k\psi = \sum \psi_{i_1...i_k}z_1^{i_1}\cdots z_k^{i_k}

and a i ψ=z iψa_i^\dagger \psi = z_i \psi, a iψ=z iψa_i\psi = \frac{\partial}{\partial z_i} \psi then the master equation says

ddtψ=Hψ\frac{d}{dt}\psi = H\psi

Each transition corresponds to an operator. There is a term a 1 n 1a k n 1a 1 m 1a k m ka_1^{\dagger n_1}\cdots a_k^{\dagger n_1} a_1^{m_1}\cdots a_k^{m_k}

In general, the master equation is written as

H= τTr(τ)(a 1 n 1(τ)a k n k(τ)a 1 m 1(τ)a k m k(τ)N 1 m̲ 1N k m̲ k) H = \sum_{\tau \in T} r(\tau) (a_1^{\dagger n_1(\tau)}\cdots a_k^{\dagger n_k(\tau)}a_1^{m_1(\tau)}\cdots a_k^{m_k(\tau)}- N_1^{\underline m_1}\cdots N_k^{\underline m_k})

We will attempt to solve HΨ=0H\Psi = 0, by trying

Ψ= i=1 k n i=0 α i n in i!z i n i= i=1 ke α iz i \Psi = \prod_{i=1}^k \sum_{n_i = 0}^\infty \frac{\alpha_i^{n_i}}{n_i!} z_i^{n_i} = \prod_{i=1}^k e^{\alpha_i z_i}
(a 1 n 1(τ)a k n k(τ))Ψ=(z 1 n 1(τ)z k n k(τ))Ψ (a_1^{\dagger n_1(\tau)}\cdots a_k^{\dagger n_k(\tau)}) \Psi = (z_1^{n_1(\tau)}\cdots z_k^{n_k(\tau)})\Psi
(a 1 m 1(τ)a k m k(τ))Ψ=(α 1 m 1(τ)α k m k(τ))Ψ (a_1^{m_1(\tau)}\cdots a_k^{m_k(\tau)}) \Psi = (\alpha_1^{m_1(\tau)}\cdots \alpha_k^{m_k(\tau)}) \Psi
(N 1 m̲ 1N k m̲ k)Ψ=(a 1 m 1(τ)a 1 m 1(τ))(a k m k(τ)a k m k(τ))Ψ=(α 1 m 1(τ)α k m k(τ))(z 1 m 1(τ)z k m k(τ))Ψ (N_1^{\underline m_1}\cdots N_k^{\underline m_k})\Psi = (a_1^{\dagger m_1(\tau)} a_1^{m_1(\tau)})\cdots (a_k^{\dagger m_k(\tau)} a_k^{m_k(\tau)})\Psi = (\alpha_1^{m_1(\tau)}\cdots \alpha_k^{m_k(\tau)})(z_1^{m_1(\tau)}\cdots z_k^{m_k(\tau)})\Psi

So the master equation acting on Ψ\Psi becomes

HΨ= τTr(τ)(α 1 m 1(τ)α k m k(τ))(z 1 n 1(τ)z k n k(τ)z 1 m 1(τ)z k m k(τ))Ψ H\Psi = \sum_{\tau \in T} r(\tau)(\alpha_1^{m_1(\tau)}\cdots \alpha_k^{m_k(\tau)})\left(z_1^{n_1(\tau)}\cdots z_k^{n_k(\tau)} - z_1^{m_1(\tau)}\cdots z_k^{m_k(\tau)}\right) \Psi

We will assume that

HΨ=0 H\Psi = 0

From which it follows that

ddz iHΨ=0 \frac{d}{d z_i}H\Psi = 0

We will make use of the compact notation for the Hamiltonian as

H= τr(τ)α m(τ)(z n(τ)z m(τ)) H = \sum_{\tau} r(\tau)\alpha^{m(\tau)}(z^{n(\tau)}-z^{m(\tau)})
ddz iHΨ=dHdz iΨ+dΨdz iH \frac{d}{d z_i}H\Psi = \frac{d H}{d z_i}\Psi + \frac{d\Psi}{dz_i}H

This becomes

τr(τ)α m(τ)(n iz n(τ)m iz m(τ))Ψ+α iHΨ=0\sum_{\tau} r(\tau)\alpha^{m(\tau)}(n_i z^{n'(\tau)}-m_i z^{m'(\tau)}) \Psi + \alpha_i H \Psi=0

but from the assumption, the second term vanishes. Hence,

τr(τ)α m(τ)(n iz n(τ)m iz m(τ))Ψ=0\sum_{\tau} r(\tau)\alpha^{m(\tau)}(n_i z^{n'(\tau)}-m_i z^{m'(\tau)}) \Psi=0

and we let

z 1=z 2=z 3=...=z k=1 z_1=z_2=z_3=...=z_k=1

which yields a solution for non-zero Ψ\Psi as

τr(τ)α m(τ)(n im i)=0\sum_{\tau} r(\tau)\alpha^{m(\tau)}(n_i -m_i)=0

which gives a solution to the rate equation, as asserted. Creation ex nihilo, something from nothing.

Some notation

Consider the nth derivative of the product of Ψ\Psi and HH.

(ddz i) nHΨ= l+m=nH (l)Ψ (m) \left(\frac{d}{d z_i}\right)^n H\Psi = \sum_{l+m=n}H^{(l)}\Psi^{(m)}

and so H (l)H^{(l)} and Ψ (m)\Psi^{(m)} become short hand for the lth and mth derivative respectively. In the case we consider here, this takes the simpler form as

l+m=nH (l)Ψ (m)=( l+m=nα i (m)H (l))Ψ\sum_{l+m=n}H^{(l)}\Psi^{(m)} = \left(\sum_{l+m=n} \alpha_i^{(m)} H^{(l)}\right)\Psi

And so we are left with the scalar α i\alpha_i to the mth power times the lth derivative of the Hamiltonian operator. We will evaluate this at

z 1=z 2=z 3=...=z k=1 z_1=z_2=z_3=...=z_k=1


H (l)= τr(τ)α m(τ)(Π q=0 l(n iq)Π q=0 l(m iq)) H^{(l)} = \sum_{\tau}r(\tau)\alpha^{m(\tau)}\left(\Pi_{q=0}^l(n_i-q) - \Pi_{q=0}^l(m_i-q)\right)

Note that the term Π q=0 l(n iq)\Pi_{q=0}^l(n_i-q) will vanish iff we count over a value where n i=qn_i=q. The same is true for the term Π q=0 l(m iq)\Pi_{q=0}^l(m_i-q).

This is the moment

We will consider the number operator NN acting ll times on HΨH\Psi

N lHΨ=N^l H \Psi=

From above we know that

HΨ= τr(τ)α m(τ)(z n(τ)z m(τ))Ψ H\Psi = \sum_{\tau} r(\tau)\alpha^{m(\tau)}(z^{n(\tau)}-z^{m(\tau)})\Psi

and also that

N q(HΨ)= l+m=q(N lH)(N mΨ)= l+m=qα i m(N lH)Ψ N^q (H \Psi) = \sum_{l+m=q}(N^l H)(N^m \Psi)= \sum_{l+m=q}\alpha_i^m (N^l H)\Psi

For now we will consider each term (N lH)Ψ(N^l H)\Psi separately.

(N lH)Ψ= τr(τ)α m(τ)(n i lz n(τ)m i lz m(τ))Ψ (N^l H)\Psi = \sum_{\tau}r(\tau)\alpha^{m(\tau)} (n_i^l z^{n(\tau)} - m_i^l z^{m(\tau)})\Psi

When we set z 1=...=z k=1z_1=...=z_k=1 we arrive at

τr(τ)α m(τ)(n i l(τ)m i l(τ)) \sum_{\tau}r(\tau)\alpha^{m(\tau)} (n_i^l(\tau) - m_i^l(\tau))
  • (Result) For l=1l=1 we recover that the mean vanishes for appropriate choice of α m(τ)\alpha^{m(\tau)}, given that the rate equation vanishes with appropriate choice of X m(τ)X^{m(\tau)}.

It is a straight forward application of the factor theorem to show that (n im i)(n_i-m_i) is a factor of (n i l(τ)m i l(τ))(n_i^l(\tau) - m_i^l(\tau)) so we have

τr(τ)α m(τ)(n i(τ)m i(τ))Poly((n i(τ),m i(τ)) \sum_{\tau}r(\tau)\alpha^{m(\tau)} (n_i(\tau) - m_i(\tau))\text{Poly}((n_i(\tau), m_i(\tau))

We know from the rate equation that

τr(τ)α m(τ)(n i(τ)m i(τ))=0 \sum_{\tau}r(\tau)\alpha^{m(\tau)} (n_i(\tau) - m_i(\tau)) = 0

Example of multiple equilibrium states

For the general terms and definitions used in this project see

Here we consider examples from Table 1 in the following paper.

P. M. Schlosser and M. Feinberg, A theory of multiple steady states in isothermal homogeneous CFSTRs with many reactions, Chemical Engineering Science 49 (1994), 1749–1767.

  • By ‘just one equilibrium state’, I presume they mean just a 1-parameter family of equilibrium states. After all, if x A,x Bx_A, x_B is an equilibrium solution of the rate equation, so is cx A,cx Bc x_A, c x_B for any c0c \ge 0.

The chemical reaction network

The example we consider here is from Table 1 (7), giving the chemical reaction network defined by

2A+B3A 2 A + B \to 3 A
3A2A+B 3 A \to 2 A + B

The input and output functions then give.

m(τ 1)={2,1} m(\tau_1)=\{2,1\}
n(τ 1)={3,0} n(\tau_1)=\{3,0\}

and for τ 2\tau_2

m(τ 2)={3,0} m(\tau_2)=\{3,0\}
n(τ 2)={2,1} n(\tau_2)=\{2,1\}

The Stochastic Petri Net

multiple equilibrium states

The chemical rate equation

The general form of the rate equation is

ddtX i= τTr(τ)X m(τ)(n i(τ)m i(τ)) \frac{d}{d t}X_i=\sum_{\tau\in T}r(\tau)X^{m(\tau)}(n_i(\tau)-m_i(\tau))

The population of species AA (BB) will be denoted X 1(t)X_1(t) (X 2(t)X_2(t)) and the rate of the reaction τ 1\tau_1 (τ 2\tau_2) as r 1r_1 (r 2r_2).

ddtX 1=r 1X 1 2X 2r 2X 1 3 \frac{d}{d t}X_1 = r_1 X_1^2 X_2 - r_2 X_1^3


ddtX 2=r 1X 1 2X 2+r 2X 1 3 \frac{d}{d t} X_2 = r_1 X_1^2 X_2 + r_2 X_1^3

The master equation

The general form of the master equation is

H= τTr(τ)(a n(τ)a m(τ)a m(τ)a m(τ)) H = \sum_{\tau\in T} r(\tau)(a^\dagger^{n(\tau)}a^{m(\tau)} - a^\dagger^{m(\tau)} a^{m(\tau)})

In our case, this becomes

H=r 1(a 1 a 1 a 1 a 1a 1a 2a 1 a 1 a 2 a 1a 1a 2)+r 2(a 1 a 1 a 2 a 1a 1a 1a 1 a 1 a 1 a 1a 1a 1) H = r_1 (a_1^\dagger a_1^\dagger a_1^\dagger a_1 a_1 a_2 - a_1^\dagger a_1^\dagger a_2^\dagger a_1 a_1 a_2) + r_2 (a_1^\dagger a_1^\dagger a_2^\dagger a_1 a_1 a_1 - a_1^\dagger a_1^\dagger a_1^\dagger a_1 a_1 a_1 )

Proof that if the rate equation vanishes then so does the mater equation

Since the rate equation is known to vanish from the Deficiency Zero Theorem, we have that

ddtX 1=0=r 1X 1 2X 2r 2X 1 3 \frac{d}{d t}X_1 = 0 = r_1 X_1^2 X_2 - r_2 X_1^3


ddtX 2=0=r 1X 1 2X 2+r 2X 1 3 \frac{d}{d t} X_2 = 0 = r_1 X_1^2 X_2 + r_2 X_1^3
  • (Algebraic independence) Let
f 1,...,f m f_1,...,f_m

be polynomials in

F[x 1,...,x n] F[x_1,...,x_n]

are called algebraically independent if there is no non-zero polynomial

AF[y 1,...,y m] A \in F[y_1,...,y_m]

such that

A(f 1,...,f m)=0 A(f_1,...,f_m)=0

then its called the annihilating polynomial.

We consider a wave function of the form

Ψ:=e bz 1e cz 2 \Psi := e^{b z_1} e^{c z_2}

We then calculate HΨH\Psi. We find solutions where this vanishes, for non-zero Ψ\Psi as

(r 1b 2cr 2b 3)z 1 3+(r 2b 3r 1b 2c)z 1 2z 2 (r_1 b^2 c - r_2 b^3)z_1^3 + (r_2 b^3 - r_1 b^2 c)z_1^2 z_2

For this to vanish, each of the terms must vanish separately. Hence we arrive at a system of equations

r 1b 2cr 2b 3=0 r_1 b^2 c - r_2 b^3 = 0
r 2b 3r 1b 2c=0 r_2 b^3 - r_1 b^2 c = 0

In this case, and in every case I’ve considered so far, a solution of one of these equations, gives a solution of the second.

category: experiments