Reaction Networks, Part 2

Table of Contents

• Introduction
• From Master to Rate
• Poisson States
• The ACK Theorem
• Symmetry Operators
• Symmetries of Equilibrium Solutions
• Noether’s Theorem
• Symmetries of Poisson States

Introduction

This is Part 2/? of a series.

In part 1 we encountered the following “levels of abstraction” of Reaction Networks 0. Equilibria: the behavior in the limit.

1. Population Dynamics: described by the Rate Equation itself. (Or continuous dynamics?)
2. Stochastic dynamics: described by the Master Equation
3. Histories: we briefly mentioned a path-integral interpretation of time-evolution.
4. Micro Dynamics: we’ve assumed a lot about this in adopting the “Reaction network” model, but haven’t addressed it directly.

My chief interest in reaction networks, again, is in understanding how the different levels of abstraction of the same system relate to each other. We can start by asking:

• When do equilibria solutions of the Master equation imply equilibria of the Rate equation, and v.v.?
• What information is lost in the act of abstracting from the Master to Rate levels, and how important is it?
• Do the differences between their solution sets amount to characterizations of solutions, or of the networks generating the equations, or both?

From Master to Rate

Here again are the Rate and Master equations side by side; see part 1 for notation.

Each of these is derived from the same reaction network, with the same parameters, but they represent different “granularities” of description of the underlying physical system. The Rate Equation describes the time evolution of a single tuple of “populations” or “concentrations” with , while the Master equation describes the time evolution of a distribution over exact species numbers , with .

As a starting point, we can ask if the former—the coarser description—can be derived by some reductive operation from the latter. How can we convert a distribution over exact states into a single concentration tuple ?

The obvious guess is to postulate that , and look at how expectations of the particle numbers evolve. This will be a function of time through the time-dependence of the underlying distribution , so we’re looking for something like an Ehrenfest Theorem from quantum mechanics. Our expectation is analogous to a Q.M. expectation , but we will omit explicit mention of in the manner of a probabilistic expectation. If we need it, we’ll write or .

We can test this guess by plugging into the Master equation and seeing what we get. We’ll work in the g.f. representation where , where expectations like are calculated by setting . Then we have:

Not too bad—we almost have a Rate Equation, except that the r.h.s. expectation is wrong—if it were we’d have the Rate Equation exactly. Instead there are two discrepancies:

• It’s a falling power
• The power is inside the expectation rather than outside of it.

The first discrepancy is the kind of thing that arises from the discreteness of and would disappear in an limit; we’ll almost always have , so the differences between and ought will be minimal.

We could maybe work around the second discrepancy by thinking of the Rate Equation as an operator equation for the number operator, though this sounds like it would require a “Heisenberg picture”, and I can’t think of what that would mean right now.

So these aren’t fatal flaws, but still—maybe we can do better.

Poisson States

How to proceed? We can try to find a state for which happens to be equal to . If such a state obeys the Master Equation, it will automatically obey the Rate Equation in the variable . The condition is:

This is an extremely general condition on the distribution —there’s no hope of finding every possible solution. But maybe we can guess our way to some solution.

The problem will get easier if we can actually evaluate one side or the other of this equation. The L.H.S. expression suggests a direction: we consider states for which , because the two factors of will cancel and the resulting expression will have the same form we started with. This at least is a simplification. Writing , where is an unspecified function, the L.H.S. of the above becomes:

This now has the same denominator as the state will have under the same ansatz. If additionally is taken to be some specific tuple to the power , then we get even closer—we get a constant multiple of the exact state :

Note that if , in any tuple component, then ; we will fudge this to mean “the product of zero terms” and therefore equal zero, rather than representing a division by zero. Then we can evaluate the whole expectation:

So for this distribution, . We have actually fixed both of the discrepancies in our pseudo-Rate Equation: the falling-power has disappeared, and the power has moved outside the expectation.

We’re nearly at the Rate-Equation-for-. All we need to do is take and normalize. Then this specific time-dependent distribution will have expectations which obey the Rate Equation.

To get the normalization constant, we need to evaluate:

This looks like the series expansion of , except that it’s written with tuples. Does it matter? Let’s check:

…and indeed, this is nothing but a multivariate exponential. (In the third line we combined the double sum into a single unrestricted sum over all the ).

Therefore we need a factor of to normalize the distribution, and the final distribution’s g.f. will be , which can be written in a few ways:

But look—the expression in the last line is nothing but a Poisson distribution:

So our full distribution is a product of independent Poissons with means: ; of course I knew this when I called them . In tuple notation we can write this whole combination as .

We could also call these product-of-Poisson distributions “coherent states”, borrowing from quantum mechanics, where the same distribution appears and acts in a similar way. I will call them “Poisson states”, however, to avoid introducing too much jargon.

Recapping, we have found that, if a Poisson state satisfies the Master Equation, then its expectation also will satisfy the Rate Equation. Apparently this works because the Poisson states are eigenstates of the operator and therefore of , but the physical meaning of this fact is obscure to me.

These results are found in a similar Q.M.-inspired notation in this paper of Baez, though the results themselves are older. I have tried, though, to present the derivation in the same way you might arrive at it yourself.

This result has two limitations:

• It requires that the Poisson state actually occurs, which is not guaranteed unless we take it as an initial condition
• Even if a Poisson state arises, there’s no guarantee it will actually remain “coherent” over time under stochastic time-evolution.

Still, it’s a starting point—we’ve learned something about the relationship between the two levels of abstraction. From here we can ask:

• Can we go the other way—does any solution to the Rate Equation always produce a Poisson state which is an equilibrium of the Master Equation?
• Under what conditions does a Poisson state remain “coherent” state over time? As a step towards this, could we characterize the first-order “divergence-from-coherence” under time evolution with ?

We’ll take on the first question next, filing the second away for later.

The ACK Theorem

The question: does any solution to the Rate Equation give us a Poisson-state solution to the Master Equation?

Here are the Rate and Master equations once again:

We can now assume we have solved the Rate with some . We construct the corresponding Poisson state and look for a condition that makes :

Well, there is our condition: the term in parentheses must vanish. The obvious way for this to happen is for each , but that would not be very interesting; it would mean all of our reactions do nothing.

Instead we can try to write this in another way. Recall that and were the output and input complexes of reaction ; a complex was just an integer combination of species, like . Then the above sum can be written as a double sum, first over all complexes, which I’ll call to distinguish from the known values. Then we write two separate sums over reactions, one for those where is equal to the input of that reaction, and one where it’s the output.

Now if the term in the brackets is identically zero, then . This is a more interesting condition! It is a condition on a solution rather than on the whole network. Physically, this means that our chosen Rate Equation solution has to produce exactly-cancelling inflows and outflows of each complex, rather than each species. This is a weak form of equilibrium where we are effectively reducing our network to only its complexes, without any regard for how these complexes are composed of the underlying species.

This condition is called “complex-balance”. This is the “Anderson-Craciun-Kurtz” or “ACK” theorem, which can found in full, in this Q.M.-derived notation, in this paper of Baez. (The proof in the paper seems needlessly complicated, though—that or mine is in error!)

Notice that it is a condition on the rate equation solution , and we can rewrite the Rate Equation in the same way:

In the last line I’ve defined to be the contribution to from inflows or outflows of complex . This is really an abuse of notation, because does not really correspond to a particular “component” of the overall concentration tuple . But if we reduced our network to only count complexes rather than species, so that the reactions converted , then this would be one of the distinct tuple-components of the reduced network.

From the above form of the Rate Equation we see that the “complex-balance” condition on automatically solves the Rate Equation as well, even if we hadn’t assumed that. If the complexes are balanced, the species must be too—this is rather obvious; the weaker condition follows from the stronger one.

To say this another way—we have encountered s of Rate-Equation solutions:

1. (Trivial) All ; this is a network that does nothing.
2. (Complex-balance) All complexes have zero time-derivative, .
3. (Detailed-balance?) All species have zero time-derivative,

(1) implies all states are solutions of (2) and (3), but this isn’t interesting. Note that this is a condition on networks rather than on states.

The solutions (2) are the solutions of the “reduced” network where we forget about how complexes are comprised of species. It is also fairly trivial that (2) implies (3), but the result here is that (2) also implies a Poisson solution to the Master Equation.

And (3) is the full set of equilibria of the master equation—being simply the definition of equilibrium.

But clearly there will be other equilibria—perhaps other classes we can characterize like (2)—a worthy line of investigation, but one which I’ll file away with the others. For now, we’ll try to say more about the complex-balanced cases in class (2).

Symmetry Operators

We have so far found that:

• Poisson solutions to the Master Equation imply solutions to the Rate equation.
• Complex-balanced solutions to the Rate Equation imply Poisson solutions to the Master Equation.

Both of these are really results about the “reduced” network on complexes, but they’re a start. Two worthy questions at this point are:

• how many of these complex-balanced solutions are there?
• are they stable?

To answer these we’ll start to look at the “symmetries”—also called “conserved quantities” or “invariants”—of a reaction network. This will also offer an example of the Poisson-state methods from the previous sections.

We consider, following Baez, a simple reversible network with two reactions: creates two molecules from a , and creates a from two s:

The corresponding Rate and Master equations are:

Note that both equations can be factored into somewhat simpler forms, because the network consists of only two reactions between the same two species.

Note also that the Rate Equation will be satisfied for any i wth . At any such concentrations the rates of the two reactions are equal and the net change in concentration is zero.

Now we observe that time-evolution by the Master equation can never alter the quantity from its initial value; hence the following operator can be said to be a “symmetry” or “conserved quantity” of the system:

Here I’ve given the name to the conserved tuple , because this tuple acts like a “particle number” of something other than the species in our reaction. In a network representing chemical reactions there ought to be such a conserved operator for each underlying element; in particle physics there would be an operator for each conserved quantum number like “lepton number”. In this example it clearly corresponds to the total number of atoms, if we interpret as representing a diatomic molecule of two s. But it is not the same as , which counts the free (as opposed to bound) molecules, so we use a different name which is still suggestive of a “particle number”. (I could see renaming to too, but I won’t for now.)

In the generating-function representation we can write the action of this operator on a state as:

In the last line I’ve used to represent a particular eigenvalue of the operator.

Apparently, any monomial is an eigenstate of :

Any state for which will have the same eigenvalue. Usually there are a great many, but for very small there might be only a few, for example has only three:

We can visualize this on a grid of . Any line of gridpoints related by a multiple of will have the same value; a few examples are visualized below:

We can also visualize how divides up the overall space of all possible . (See my drawing dimensions post for more on this diagram style):

If there were two independent symmetries , they would carve out two independent degrees-of-freedom from the full state space of :

We can see this as describing the “coordinates” of a term :

• first pick a value for each independent symmetry
• then, choose any coordinate in the remaining subspace . In the present case, once is known, either or suffices to determine the exact state. Of course, all of this applies only to “pure” monomial states—a general will be composed of components along many . This decomposition is:

This is analogous to a decomposition of a quantum state into, first, energy levels, then some other quantum number like orbital angular momentum:

In the decomposition “basis”, the operator will be block-diagonal:

And we could represent the space of all as a sum of subspaces of fixed :

I’ve drawn these as double lines because the space for a fixed is still -dimensional. A matrix representation of itself will be block diagonal for the same reason. This diagram conveys that any state defined on the full space of can be factored into one independent (multi-dimensional) coordinate on each subspace. Note that this is describing the space of , which is much larger than the space of a single monomial. A monomial has a single value while potentially has a component along each value—how best to render this? These diagrams are a stab at relating my mental model, but I think they haven’t reached their final form yet.

Symmetries of Equilibrium Solutions

Now we’ll consider what a symmetry of the “Hamiltonian” can tell us about its solutions, in two senses:

1. What can a -symmetry tell us about general equilibrium solutions () to the Master equation?
2. If we start with a Rate Equation equilibrium solution, and construct a Poisson solution to the Master Equation using the ACK Theorem, what can a symmetry tell us about this solution in particular?

First question first. Let us continue to work on our simple example network:

And let us conjure again the “Hamiltonian” for this network:

What is the action of this on a general state? Even the action of this on a simple monomial will be a fairly large expression; to keep it under control we write it with tuple , representing the change in species induced by reaction . For this network , so we only need the one; I’ll write it simply as . Note this is orthogonal to the conserved tuple :

.

Then the action of is:

Evidently the monomial makes contributions to three states, all with the same eigenvalue of our symmetry as the original:

1. Adds to state with
2. Adds to state with
3. Subtracts from the original state

If we think of our as representing a flow , then any final state will receive contributions from itself and the two neighboring states; and we can reindex the above expression to reflect this:

Now if a state is in equilibrium with , then as we’ve just seen, cannot change the value of , so itself must decompose cleanly into a sum over :

And apparently, any equilibrium solution must have on every -component! And, if we find some set of equilibria on each component, then any arbitrary linear combination of these will, after normalization, be an equilibrium solution to the full .

The above expression has the form:

Which gives us the condition for equilibrium on a eigenstate in terms of the same constants:

where I’ve replaced the powers above with the “components” of as in the sum , and dropped the primes. Note that the coefficients are still functions of .

The above is a 3-term recursion relation for the components of the equilibrium distribution . We have come up with a discrete boundary-value-problem for any equilibrium solution, with boundary conditions arising at either or because neither particle-number is allowed to be negative.

Consequently, there should only be at most one equilibrium solution to the Master Equation for each value of .

That’s still a lot of equilibrium states, but at least we know something!

(This is the claim of Baez et al, but I’m a bit skeptical of this. For one thing, the recursion relation skips the states of immediately-adjacent , and it seems possible to contrive a system where both interlaced recursions can be satisfied at once, or neither can. We might get something like sine/cosine modes. Filing it away.)

The following depicts how the state space is reduced in equilibrium. A linear combination over the independent equilibria will suffice to describe any equilibrium solution.

(Note to self: find a coherent way to depict “infinitely many axes”)

If our has multiple symmetries , then any equilibria will clearly decompose along each symmetry, considered independently. But in general the equilibrium state of one symmetry will be some linear combination of those of the other; we don’t expect them to be “simultaneously diagonalizeable” unless one of the symmetries implies the other.

Note that this problem has turned out to be quite similar to a discrete Laplacian operator:

This is easier to see if we reindex our states by and :

Our three coefficients, however, are complicated functions of and , so it is very unlikely this expression would be the discrete form of any reasonable continuous equation.

I wonder at this point:

• Can we say anything about the relationship between the unique equilibria of neighboring values?

But we’ll put this aside for now, along with the question of stability posed above.

Noether’s Theorem

We can summarize the results of the previous section in a more general way.

We have seen, for the Master Equation:

• If can’t change the value of an operator on a single state, then time-evolution with can’t take us out of the subspace of constant .
• Furthermore, time-evolution will operate separately on components of distinct , and any equilibria will be a linear combination of independent equilibria on each component.

These facts amount to Noether’s Theorem for reaction networks. The intuitive form of this theorem is:

• if, for some operator , every “transition” defined by the operator links states with the same eigenvalue of
• …then no amount of time-evolution can ever take a state of one into another; eigenvalues of are therefore conserved by time-evolution.

This relationship between and turns out to have a number of equivalent mathematical expressions. A few of these are:

To read these, recall that our are p.d.f.s, so an expectation is taken and a time derivative is

Thus the last line is

Which is clearly equivalent to the previous line. It is also interesting to note that

for any and any , since is “infinitesimal-stochastic”—its columns sum to zero—leading to an expectation of zero in any state. This means that the expectations like can be trivially written in terms of the commutator

Beyond these observations I won’t prove these statements—they follow straightforwardly from the definitions. For more see Baez’ Noether paper and book.

We will use the first form when we want to say “operator is a symmetry of ”.

This Noether theorem is fairly trivial: is constant on connected components? Big deal. It is most interesting as an elementary example for reference when considering more complicated cases—of how symmetry operators of time-evolution lead to conserved quantities, and, more interestingly to me, how the eigenvalues of symmetry operators supply us with coordinates by which we identify states, and dynamics interact with these coordinates—a given dynamical law’s effect on the state can only be to either:

• move “along” coordinates of a symmetry (if )
• move “between” coordinates of a symmetry (like an operator—what is the general rule here?)
• or to “mix” the coordinate of a symmetry (the general case)

Symmetries of Poisson States

Now we’ll move on to the second question from above. This was: 2. If we start with a Rate Equation equilibrium solution, and construct a Poisson solution to the Master Equation using the ACK Theorem, what can a symmetry tell us about this solution in particular?

Any Rate-Equation equilibrium solution for our example network will have .

Then the ACK theorem implies that the following Poisson state is an equilibrium solution to the Master equation:

We can verify this:

By the discussion right before this, we know this solution is a particular linear combination of solutions with specific values. However there is no reason to think this solution is an eigenstate of and indeed it is not:

(Here we’ve used the notation for a “scalar product of tuples”; it looks like a dot product over any number of arguments.)

So a Poisson state is not an eigenstate of —but of course it isn’t, because the eigenstates of were the monomials. Instead acts on Poisson states to transform them in some obscure way. But whatever this thing is, we expect it to also be an equilibrium, because our original state was a superposition of equilibria for each , and won’t mix these. Let’s verify:

Gross, but it works. In the last line we used , which was the condition that made a solution of the Rate Equation.

’s action is to only skew all the subspaces by , so it can’t take us out of equilibrium. This means we can create many equilibria from our original Poisson equilibria by applying different powers of the operator. Any linear combination of single- equilibria would work, but these won’t be Poisson states in general.

We might now take some inspiration from the form of the time-evolution operator , and ask what it means to write a one-dimensional “shift-by- operator”:

Here is a continuous parameter which determines how much to shift by. On a monomial (a eigenstate) this operator is simply a rescaling:

What does do to a coherent state? We can find this a couple ways. The simplest is to write as a series in . But it will be slightly more interesting to note that each of the terms in can be written as , and therefore:

So our “shift-by-” operator is actually a scale by operator. (This might have been obvious if you recognize as a generator of scale transformations, but I didn’t remember this. See wiki.)

Therefore the action of this operator on a coherent state is (switching to now):

Which to be clear, if my notation has lost you, means:

But the second exponential is another Poisson state, with a transformed :

The transformed Poisson state is not normalized, though. Instead we have:

Observe also that this new concentration tuple still obeys the Rate equation “complex-balanced” condition!

We expected this, because no power of can take any state off the space of constant ; just as no power of could take any state off a plane of constant .

So we’ve found that this operator defines a group flow through a one-dimensional space of complex-balanced Rate Equation solutions. And furthermore, every complex-balanced solution should be reachable with some (for this simple network, at least).

Some reasonable questions at this point are:

• What happens with more reactions and more symmetries? Can we reach every Rate-equation equilibria from any starting point via some combination of scaling operators? Can we leave equilibria this way?
• What about non-complex-balanced equilibria, the weaker condition? Can we say anything about them at all? Perhaps the symmetry operator will help?
• Still, we keep punting this: what about stability?

We’ll leave it there for now, though.