<![CDATA[The Large-Number Limit for Reaction Networks (Part 3)]]>

<![CDATA[The Large-Number Limit for Reaction Networks (Part 3)]]>joint with Arjun Jain

We used to talk about reaction networks quite a lot here. When Arjun Jain was visiting the CQT, we made a lot of progress understanding how the master equation reduces to the rate equation in the limit where there are very large numbers of things of each kind. But we never told you the end of the story, and by now it’s been such a long time that you may need a reminder of some basic facts!

So…

The rate equation treats the number of things of each kind as continuous—a nonnegative real number—and says how it changes in a deterministic way.

The master equation treats the number of things of each kind as discrete—a nonnegative integer—and says how it changes in a probabilistic way.

You can think of the master equation as the ‘true’ description, and the rate equation as an approximation that’s good in some limit where there are large numbers of molecules — or more precisely, where the probability distribution of having some number of molecules of each kind is sharply peaked near some large value.

You may remember that in the master equation, the state of a chemical system is described by a vector $\psi$ in a kind of ‘Fock space’, while time evolution is described with the help of an operator on this space, called the ‘Hamiltonian’ $H$ :

$\displaystyle{ \frac{d}{dt} \psi(t) = H \psi(t) }$

The Hamiltonian is built from annihilation and creation operators, so all of this looks very much like quantum field theory. The details are here, and we won’t try to review them all:

• John Baez and Jacob Biamonte, Quantum Techniques for Stochastic Mechanics.

The point is this: the ‘large-number limit’ where the master equation reduces to the rate equation smells a lot like the ‘classical limit’ of quantum field theory, where the description of light in terms of photons reduces to the good old Maxwell equations. So, maybe we can understand the large-number limit by borrowing techniques from quantum field theory!

How do we take the classical limit of quantum electromagnetism and get the classical Maxwell equations? For simplicity let’s ignore charged particles and consider the ‘free electromagnetic field’: just photons, described by the quantum version of Maxwell’s equations. When we take the classical limit we let Planck’s constant $\hbar$ go to zero: that much is obvious. However, that’s not all! The energy of each photon is proportional to $\hbar,$ so to take the classical limit and get a solution of the classical Maxwell’s equations with nonzero energy we also need to increase the number of photons. We cleverly do this in such a way that the total energy remains constant as $\hbar \to 0.$

So, in quantum electromagnetism the classical limit is also a large-number limit!

That’s a good sign. It suggests the same math will also apply to our reaction network problem.

But then we hit an apparent roadblock. What’s the analogue of Planck’s constant in chemical reaction networks? What should go to zero?

We told you the answer to that puzzle a while ago: it’s the reciprocal of Avogadro’s number!

You see, chemists measure numbers of molecules in ‘moles’. There’s a large number of molecules in each mole: Avogadro’s number. If we let the reciprocal of Avogadro’s number go to zero, we are taking a limit where chemistry becomes ‘continuous’ and the discreteness of molecules goes away. Of course this is just a mathematical trick, but it’s a very useful one.

So, we got around that roadblock. And then something nice happened.

When taking the classical limit of quantum electromagnetism, we focus attention on certain quantum states that are the ‘best approximations’ to classical states. These are called ‘coherent states’, and it’s very easy to study how the behave as we simultaneously let $\hbar \to 0$ and let the expected number of photons go to infinity.

And the nice thing is that these coherent states are also important in chemistry! But because chemistry involves probabilities rather than amplitudes, they have a different name: ‘Poisson distributions’. On this blog, Brendan Fong used them to give a new proof of a great result in mathematical chemistry, the Anderson–Craciun–Kurtz theorem.

So, we have most of the concepts and tools in place, and we can tackle the large-number limit using quantum techniques.

You can review the details here:

• The large-number limit for reaction networks (part 1).

• The large-number limit for reaction networks (part 2) .

So, after a quick refresher on the notation, we’ll plunge right in.

As you’ll see, we solve the problem except for one important technical detail: passing a derivative through a limit! This means our main result is not a theorem. Rather, it’s an idea for how to prove a theorem. Or if we act like physicists, we can call it a theorem.

Review of notation

The rate equation says

$\displaystyle{\frac{d}{dt}{x}(t) = \sum_{\tau \in T} r(\tau) (t(\tau)-s(\tau)) {x}(t)^{s(\tau)} }$

where:

• $x(t) \in \mathbb{R}^k$ is a vector describing concentrations of $k$ different species at time $t.$ In chemistry these species could be different kinds of molecules.

• Each $\tau \in T$ is a transition, or in chemistry, a reaction.

• $s(\tau) \in \mathbb{N}^k$ is a vector of natural numbers saying how many items of each species appear as inputs to the reaction $\tau.$ This is called the source of the reaction.

• $t(\tau) \in \mathbb{N}^k$ is a vector of natural numbers saying how many items of each species appear as outputs of the reaction $\tau.$ This is called the target of the reaction. So, $t(\tau) - s(\tau)$ says the net change of the number of items of each species in the reaction $\tau.$

• The rate at which the reaction $\tau$ occurs is proportional to the rate constant $r(\tau)$ times the number

$x(t)^{s(\tau)}$

Here we are raising a vector to a vector power and getting a number, using this rule:

$x^r = x_1^{r_1} \cdots x_k^{r_k}$

where $r$ is any vector of natural numbers $(r_1, \dots, r_k),$ and $x$ is any vector of nonnegative real numbers. From now on we’ll call a vector of natural numbers a multi-index.

In this paper:

• John Baez, Quantum techniques for reaction networks.

it was shown that the master equation implies

$\displaystyle{\frac{d}{dt}\langle N_i \psi(t)\rangle = \sum_{\tau \in T} r(\tau) (t_i(\tau)-s_i(\tau)) \; \langle N^{\, \underline{s(\tau)}} \, \psi(t)\rangle }$

Here:

• $\psi(t)$ is the stochastic state saying the probability of having any particular number of items of each species at each time $t.$ We won’t review the precise details of how this work; for that reread the relevant bit of Part 8.

• $N_i$ is the $i$ th number operator, defined using annihilation and creation operators as in quantum mechanics:

$N_i = a_i^\dagger a_i$

For the annihilation and creation operators, see Part 8.

• $\langle N_i \psi(t) \rangle$ is the expected number of items of the $i$ th species at time $t.$

• Similarly, $\langle N^{\underline{s(\tau)}} \psi(t)\rangle$ is the expected value of a certain product of operators. For any multi-index $r$ we define the falling power

$N_i^{\underline{r}_i} = N_i (N_i - 1) (N_i - 2) \cdots (N_i - r_i +1)$

and then we define

$N^{\underline{r}} = N_1^{\underline{r_1}} \cdots N_k^{\underline{r_k}}$

The large-number limit

Okay. Even if you don’t understand any of what we just said, you’ll see the master and rate equation look similar. The master equation implies this:

$\displaystyle{\frac{d}{dt}\langle N_i \psi(t)\rangle = \sum_{\tau \in T} r(\tau) (t_i(\tau)-s_i(\tau)) \; \langle N^{\, \underline{s(\tau)}} \, \psi(t)\rangle }$

while the rate equation says this:

$\displaystyle{\frac{d}{dt}{x}(t) = \sum_{\tau \in T} r(\tau) (t(\tau)-s(\tau)) \; {x}(t)^{s(\tau)} }$

So, we will try to get from the first to the second second with the help of a ‘large-number limit’.

We start with a few definitions. We introduce an adjustable dimensionless quantity which we call $\hbar.$ This is just a positive number, which has nothing to do with quantum theory except that we’re using a mathematical analogy to quantum mechanics to motivate everything we’re doing.

Definition. The rescaled number operators are defined as $\widetilde{N}_i = \hbar N_i.$ This can be thought of as a rescaling of the number of objects, so that instead of counting objects individually, we count them in bunches of size $1/\hbar.$

Definition. For any multi-index $r$ we define the rescaled falling power of the number operator $N_i$ by:

$\widetilde{N}_i^{\underline{r_i}} = \widetilde{N}_i (\widetilde{N}_i - \hbar)(\widetilde{N}_i-2\hbar)\cdots (\widetilde{N}_i-r_i\hbar+\hbar)$

and also define

$\widetilde{N}^{\underline{r}} = \widetilde{N}_1^{\underline{r_1}} \; \cdots \;\widetilde{N}_k^{\underline{r_k}}$

for any multi-index $r.$

Using these, we get the following equation:

$\displaystyle{\frac{1}{\hbar}\frac{d}{dt} \langle \widetilde{N}_i \psi(t)\rangle = \sum_{\tau \in T} r(\tau) (t_i(\tau)-s_i(\tau)) \; \langle \widetilde{N}^{\underline{s(\tau)}} \psi(t)\rangle \; \frac{1}{\hbar^{|s(\tau)|}} }$

where for any multi-index $r$ we set

$|r| = r_1 + \cdots + r_k$

This suggests a way to rescale the rate constants in the master equation:

Definition. The rescaled rate constants are

$\displaystyle{\widetilde{r}(\tau) = \frac{r(\tau)}{\hbar^{|s(\tau)|-1}}}$

From here onwards, we change our viewpoint. We consider the rescaled rate constants $\widetilde{r}(\tau)$ to be fixed, instead of the original rate constants $r(\tau).$ So, as we decrease $\hbar,$ we are studying situations where the original rate constants change to ensure that the rescaled rate constants stays fixed!

So, we switch to working with a rescaled master equation:

Definition. The rescaled master equation is:

$\displaystyle{\frac{d}{dt} \langle \widetilde{N}_i\widetilde{\psi}(t)\rangle = \sum_{\tau \in T} \widetilde{r}(\tau) (t_i(\tau)-s_i(\tau)) \; \langle \widetilde{N}^{\underline{s(\tau)}} \widetilde{\psi}(t)\rangle }$

This is really a one-parameter family of equations, depending on $\hbar\in (0,\infty).$ We write a solution of the rescaled master equation as $\widetilde{\psi}(t),$ but it is really one solution for each value of $\hbar.$

Following the same procedure as above, we can rescale the rate equation, using the same definition of the rescaled rate constants:

Definition. The rescaled number of objects of the $i\mathrm{th}$ species is defined as $\widetilde{x_i}=\hbar x_i,$ where $x_i$ is the original number of objects of the $i\mathrm{th}$ species. Here again, we are counting in bunches of $1/\hbar.$

Using this to rescale the rate equation, we get

Definition. The rescaled rate equation is

$\displaystyle{\frac{d}{dt}\widetilde{x}(t) = \sum_{\tau \in T} \widetilde{r}(\tau) (t(\tau)-s(\tau))\; \widetilde{x}(t)^{s(\tau)} }$

where

$\widetilde{x}(t)=(\widetilde{x_1}(t),\widetilde{x_2}(t),\dots, \widetilde{x_k}(t))$

Therefore, to go from the rescaled master equation to the rescaled rate equation, we require that

$\langle\widetilde{N}^{\underline{r}}\, \widetilde{\psi}(t)\rangle \to \langle\widetilde{N}\widetilde{\psi}(t)\rangle^r$

as $\hbar \to 0.$ If this holds, we can identify $\langle\widetilde{N}\widetilde{\psi}(t)\rangle$ with $\widetilde{x}(t)$ and get the rate equation from the master equation!

To this end, we introduce the following crucial idea:

Definition. A semiclassical family of states, $\widetilde{\psi},$ is defined as a one-parameter family of states depending on $\hbar \in (0,\infty)$ such that for some $\widetilde{c} \in [0,\infty)^k$ we have

$\langle\widetilde{N}^r\widetilde{\psi}\rangle \to \widetilde{c}^{\, r}$

for every $r\in \mathbb{N}^k$ as $\hbar \to 0.$

In particular, this implies

$\langle\widetilde{N}_i\widetilde{\psi}\rangle \to \widetilde{c}_i$

for every index $i.$

Intuitively, a semiclassical family is a family of probability distributions that becomes more and more sharply peaked with a larger and larger mean as $\hbar$ decreases. We would like to show that in this limit, the rescaled master equation gives the rescaled rate equation.

We make this precise in the following propositions.

Proposition 1. If $\widetilde{\psi}$ is a semiclassical family as defined above, then in the $\hbar \to 0$ limit, we have $\langle\widetilde{N}^{\underline{r}}\widetilde{\psi}\rangle \to \widetilde{c}^{\; r}$ as well.

Proof. For each index $i,$

$\displaystyle{ \langle\widetilde{N}_i^{\; \underline{r_i}}\, \widetilde{\psi}\rangle = \displaystyle{ \langle \widetilde{N}_i (\widetilde{N}_i - \hbar)(\widetilde{N}_i-2\hbar)\cdots(\widetilde{N}_i-r_i\hbar+\hbar)\,\widetilde{\psi}\rangle} }$

$\displaystyle{ = \Big\langle\Big(\widetilde{N}_i^r + \hbar\frac{r_i(r_i-1)}{2}\widetilde{N}_i^{r_i-1}+\cdots + \hbar^{r_i-1}(r_i-1)!\Big)\,\widetilde{\psi}\Big\rangle }$

By the definition of a semiclassical family,

$\displaystyle{ \lim_{\hbar\to 0} \langle\Big(\widetilde{N}_i^{r_i} + \hbar\frac{r_i(r_i-1)}{2}\widetilde{N}_i^{r_i-1}+ \cdots + \hbar^{r_i-1}(r_i-1)!\Big)\;\widetilde{\psi}\rangle} = \widetilde{c}_i^{\; r_i}$

since every term but the first approaches zero. Thus, we have

$\displaystyle{ \lim_{\hbar \to 0} \langle\widetilde{N}_i^{\; \underline{r_i}}\, \widetilde{\psi}\rangle = \widetilde{c}_i^{\; r_i} }$

A similar but more elaborate calculation shows that

$\displaystyle{ \lim_{\hbar \to 0} \langle\widetilde{N}_1^{\, \underline{r_1}} \cdots\widetilde{N}_k^{\, \underline{r_k}} \, \widetilde{\psi}\rangle = \lim_{\hbar \to 0} \langle\widetilde{N}_1^{\, r_1}\cdots \widetilde{N}_k^{\, r_k} \widetilde{\psi}\rangle= \lim_{\hbar \to 0}\langle\widetilde{N}^{\, r} \, \widetilde{\psi}\rangle }$

or in other words

$\langle\widetilde{N}^{\, \underline{r}}\,\widetilde{\psi}\rangle \to \widetilde{c}^{\, r}$

as $\hbar \to 0.$ █

Proposition 2. If $\widetilde{\psi}$ is a semiclassical family of states, then

$\displaystyle{ \langle (\widetilde{N}-\widetilde{c})^{r}\, \widetilde{\psi}\rangle \to 0 }$

for any multi-index $r.$

Proof. Consider the $r_i\mathrm{th}$ centered moment of the $i\mathrm{th}$ number operator:

$\displaystyle{\langle(\widetilde{N}_i-\widetilde{c}_i)^{r_i}\widetilde{\psi}\rangle = \sum_{p =0}^{r_i} {r_i \choose p}\langle\widetilde{N}_i^p\widetilde{\psi}\rangle(-\widetilde{c}_i)^{r_i-p} }$

Taking the limit as $\hbar$ goes to zero, this becomes

$\begin{array}{ccl} \displaystyle{ \lim_{\hbar \to 0}\sum_{p =0}^{r_i} {r_i \choose p}\langle\widetilde{N}_i^p\widetilde{\psi}\rangle(-\widetilde{c}_i)^{r_i-p} } &=& \displaystyle{ \sum_{p =0}^{r_i} {r_i \choose p}(\widetilde{c}_i)^p(-\widetilde{c}_i)^{r_i-p} } \\ \\ &=& (\widetilde{c}_i-\widetilde{c}_i)^{r_i} \\ \\ &=& 0 \end{array}$

For a general multi-index $r$ we can prove the same sort of thing with a more elaborate calculation. First note that

$\langle (\widetilde{N}-\widetilde{c})^{r}\widetilde{\psi}\rangle=\langle(\widetilde{N_1}-\widetilde{c_1})^{r_1} \cdots (\widetilde{N_k}-\widetilde{c_k})^{r_k})\widetilde{\psi}\rangle$

The right-hand side can be expanded as

$\displaystyle{ \langle(\sum_{p_1 =0}^{r_1} {r_1 \choose p_1}\widetilde{N}_1^{p_1}(-\widetilde{c}_1)^{r_1-p_1} ) \cdots (\sum_{p_k =0}^{r_k} {r_k \choose p_k}\widetilde{N}_k^{p_k}(-\widetilde{c}_k)^{r_k-p_k} )\widetilde{\psi} }\rangle$

We can write this more tersely as

$\displaystyle{ \sum_{p=0}^r} {r\choose p} \langle \widetilde{N}^p\widetilde{\psi}\rangle (-\widetilde{c})^{r-p}$

where for any multi-index $r$ we define

$\displaystyle{{\sum_{p=0}^{r}}= \sum_{p_1 =0}^{r_1} \cdots \sum_{p_k =0}^{r_k} }$

and for any multi-indices $r, p$ we define

$\displaystyle{ {r \choose p}={r_1 \choose p_1}{r_2 \choose p_2}\cdots {r_k \choose p_k}}$

Now using the definition of a semiclassical state, we see

$\displaystyle{ \lim_{\hbar \to 0} \sum_{p=0}^r} {r\choose p} \langle \widetilde{N}^p\widetilde{\psi}\rangle (-\widetilde{c})^{r-p}= \displaystyle{ \sum_{p=0}^r} {r\choose p} (\widetilde{c})^{p} (-\widetilde{c})^{r-p}$

But this equals zero, as the last expression, expanded, is

$\displaystyle{ (\widetilde{c})^r \left( \sum_{p_1=0}^{r_1} {r_1\choose p_1} (-1)^{r_1-p_1}\right) \cdots \left( \sum_{p_k=0}^{r_k} {r_k\choose p_k} (-1)^{r_k-p_k} \right) }$

where each individual sum is zero. █

Here is the theorem that would finish the job if we could give a fully rigorous proof:

“Theorem.” If $\widetilde{\psi}(t)$ is a solution of the rescaled master equation and also a semiclassical family for the time interval $[t_0,t_1],$ then $\widetilde{x}(t) = \langle \widetilde{N} \widetilde{\psi}(t) \rangle$ is a solution of the rescaled rate equation for $t \in [t_0,t_1].$

Proof sketch. We sketch a proof that relies on the assumption that we can pass the $\hbar \to 0$ limit through a time derivative. Of course, to make this rigorous, we would need to justify this. Perhaps it is true only in certain cases.

Assuming that we can pass the limit through the derivative:

$\displaystyle{\lim_{\hbar \to 0}\frac{d}{dt} \langle \widetilde{N}\widetilde{\psi}(t)\rangle = \lim_{\hbar \to 0} \sum_{\tau \in T} \widetilde{r}(\tau) (t(\tau)-s(\tau))\langle \widetilde{N}^{\, \underline{s(\tau)}} \, \widetilde{\psi}(t)\rangle }$

and thus

$\displaystyle{\frac{d}{dt}\lim_{\hbar \to 0} \langle \widetilde{N}\widetilde{\psi}(t)\rangle = \sum_{\tau \in T} \widetilde{r}(\tau) (t(\tau)-s(\tau)) \lim_{\hbar \to 0}\langle \widetilde{N}^{\, \underline{s(\tau)}} \, \widetilde{\psi}(t)\rangle }$

and thus

$\displaystyle{\frac{d}{dt}\widetilde{x}(t) = \sum_{\tau \in T} \widetilde{r}(\tau) (t(\tau)-s(\tau))\widetilde{x}^{\, s(\tau)} }.$

As expected, we obtain the rescaled rate equation. █

Another question is this: if we start with a semiclassical family of states as our initial data, does it remain semiclassical as we evolve it in time? This will probably be true only in certain cases.

An example: rescaled coherent states

The best-behaved semiclassical states are the coherent states.
Consider the family of coherent states

$\displaystyle{\widetilde{\psi}_{\widetilde{c}} = \frac{e^{(\widetilde{c}/\hbar) z}}{e^{\widetilde{c}/\hbar}}}$

using the notation developed in the earlier mentioned paper. In that paper it was shown that for any multi-index $m$ and any coherent state $\Psi$ we have

$\langle N^{\underline{m}}\Psi\rangle = \langle N\Psi \rangle^m$

Using this result for $\widetilde{\psi}_{\widetilde{c}}$ we get

$\displaystyle{\langle \widetilde{N}^{\underline{m}}\widetilde{\psi}_{\widetilde{c}}\rangle~=~\hbar^{|m|}\langle N^{\underline{m}}\widetilde{\psi}_{\widetilde{c}}\rangle~=~\hbar^{|m|}\langle N\widetilde{\psi}_{\widetilde{c}}\rangle^m~=~\hbar^{|m|}\frac{\widetilde{c}^m}{\hbar^{|m|}}~=~\widetilde{c}^m}$

Since $\langle\widetilde{N}^{\underline{m}}\widetilde{\psi}_{\widetilde{c}}\rangle$ equals $\langle \widetilde{N}^{m} \widetilde{\psi}_{\widetilde{c}}\rangle$ plus terms of order $\hbar,$ as $\hbar \to 0$ we have

$\langle\widetilde{N}^{\underline{m}}\widetilde{\psi}_{\widetilde{c}}\rangle~\to~\langle\widetilde{N}^{m}\widetilde{\psi}_{\widetilde{c}}\rangle=\widetilde{c}^{m}$

showing that our chosen $\widetilde{\psi}_{\widetilde{c}}$ is indeed a semiclassical family.

]]>