Toric Geometry in Reaction Networks

I want to figure out how to use toric geometry in chemistry. This is a good intro to toric geometry:

• William Fulton, Introduction to Toric Varieties, Princeton U. Press, 1993.

and this is a great explanation of how it shows up in chemistry:

• Mercedes Perez Millan, Alicia Dickenstein, Anne Shiu and Carsten Conradi, Chemical reaction systems with toric steady states, Bulletin of Mathematical Biology 74 (2012), 1027–1065.

You don’t need to read Fulton’s book to understand this paper! But you don’t need to read either to understand what I’m about to say. It’s very simple.

Suppose we have a bunch of chemical reactions. For example, just one:

\mathrm{A} \mathrel{\substack{\alpha_{\rightarrow} \\\longleftrightarrow\\ \alpha_{\leftarrow}}} \mathrm{B} + \mathrm{C}

or more precisely two: the forward reaction

\mathrm{A} \to \mathrm{B} + \mathrm{C}

with its rate constant \alpha_\to, and the reverse reaction

\mathrm{B} + \mathrm{C} \to \mathrm{A}

with its rate rate constant \alpha_{\rightarrow}. Then as I recently explained, these reactions are in a detailed balanced equilibrium when

\alpha_{\to} [\mathrm{A}] = \alpha_{\rightarrow} [\mathrm{B}] [\mathrm{C}]

This says the forward reaction is happening at the same rate as the reverse reaction.

Note: we have three variables, the concentrations [\mathrm{A}], [\mathrm{B}] and [\mathrm{C}], and they obey a polynomial equation. But it’s a special kind of polynomial equation! It just says that one monomial—a product of variables, times a constant—equals another monomial. That’s the kind of equation that’s allowed in toric geometry.

Let’s look at another example:

\mathrm{B} + \mathrm{C} \mathrel{\substack{\beta_{\rightarrow} \\\longleftrightarrow\\ \beta_{\leftarrow}}} \mathrm{D} + \mathrm{D} + \mathrm{A}

Now we have a detailed balance equilibrium when

\beta_{\to} [\mathrm{B}] [\mathrm{C}] = \beta_{\leftarrow} [\mathrm{D}]^2 [\mathrm{A}]

Again, one monomial equals another monomial.

Now let’s look at a bigger reaction network, formed by combining the two so far:

\mathrm{A} \mathrel{\substack{\alpha_{\rightarrow} \\\longleftrightarrow\\ \alpha_{\leftarrow}}} \mathrm{B} + \mathrm{C} \mathrel{\substack{\beta_{\rightarrow} \\\longleftrightarrow\\ \beta_{\leftarrow}}} \mathrm{D} + \mathrm{D} + \mathrm{A}

Detailed balance is a very strong condition: it says that each reaction is occurring at the same rate as its reverse. So, it happens when

\alpha_{\to} [\mathrm{A}] = \alpha_{\rightarrow} [\mathrm{B}] [\mathrm{C}]

and

\beta_{\to} [\mathrm{B}] [\mathrm{C}] = \beta_{\leftarrow} [\mathrm{D}]^2 [\mathrm{A}]

So, we can have more than one equation, but all of them simply equate two monomials. That’s how it always works in a detailed balanced equilibrium.

Definition. An affine toric variety is a subset of \mathbb{R}^n defined by a system of equations, each of which equates two monomials in the coordinates x_1, \dots, x_n.

So, if we ignore the restriction that our variables should be ≥ 0, the space of detailed balanced equilibria for a reaction network where every reaction is reversible is an affine toric variety. And the point is, there’s a lot one can say about such spaces!

A simple example of an affine toric variety is the twisted cubic, which is the subset

\{ (x,x^2,x^3) \} \subset \mathbb{R}^3

Here it is, as drawn by Claudio Rocchini:

I may say more about this, but today I just wanted to get the ball rolling.

Puzzle. What’s a reaction network whose detailed balanced equilibrium equations give the twisted cubic?

This entry was posted on Tuesday, July 3rd, 2018 at 6:20 pm and is filed under chemistry, mathematics. You can follow any responses to this entry through the RSS 2.0 feed. You can leave a response, or trackback from your own site.

18 Responses to Toric Geometry in Reaction Networks

  1. tomate says:
    3 July, 2018 at 6:43 pm

    I read this one by Carciun Shiu et al. https://arxiv.org/abs/0708.3431.

    The way I understood it is as follows.

    Let W be the rate matrix of a continuous-time Markov process on a finite state space (equivalently: of a linear chemical network where each complex is constituted by exactly one specie).

    Let c be its steady state Wc = 0 (assumed unique). We want to find an expression for p in terms of the rates. One-line derivation: The determinant of W is zero. Consider the Laplace expansion

    0 = det W = sum_j (-1)^{i+j} W_{ij} det W(i,j) (1)

    where W(i,j) is the matrix where you remove the i-th row and the j-th column (or the other way around…).

    Then c_j = (-1)^{i+j} det W(i,j) (up to signs…), independently of i. Turns up that this has a representation in terms of spanning trees by the matrix-tree theorem, but this is not relevant here.

    Now, consider a chemical network with mean-field equations. If deficiency = 0, we can still formally solve the above equation for the complexes, but we want to know the steady state in terms of the species. So for example if the first complex is C1= X1+2 X2 then c1 = x1 x2^2. Now Eq. (1) becomes an algebraic equation for the concentrations, which (if I understood well) describes a toric algebraic variety.

    If deficiency /= 0, I have no idea.

    Reply
  2. tomate says:
    3 July, 2018 at 6:44 pm

    an expression for p

    Reply
  3. tomate says:
    3 July, 2018 at 6:44 pm

    an expression for c !!!

    Reply
  4. leporcinet11 says:
    4 July, 2018 at 10:30 am

    Pretty interesting stuff. But I wonder, in reality, is the rate constant constant ? does the exponent on the concentration always matches the stoichiometric coefficient ?

    Reply
    • John Baez says:
      4 July, 2018 at 6:36 pm

      You’re asking if the law of mass action is exactly true. No: no laws are exactly true. The law of mass action can be derived from certain assumptions: for example, one should assume a well-mixed solution with molecules moving around randomly, where each reaction occurs with a given probability when the reactants get close enough to each other, moving in the right range of speeds, with enough molecules of each kind that it’s good to approximate the resulting stochastic differential equation by a deterministic time evolution for the mean values of the number of each kind of molecule.

      In reality these assumptions are often good approximations, but not always. They are probably not very good for complex reactions like this one in our description of the citric acid cycle:

      \alpha-\mathrm{ketoglutarate} + \mathrm{NAD}^+ + \text{CoA-SH} \longleftrightarrow \mathrm{succinyl-CoA} + \mathrm{NADH} + \mathrm{H}^+ + \textrm{CO}_2

      I doubt that α-ketoglutarate, NAD+ and CoA-SH all meet simultaneously and react; more likely two of them meet, react (probably with the help of an enzyme), and then meet the third.

      The Wikipedia article Law of mass action says:

      In biochemistry, there has been significant interest in the appropriate mathematical model for chemical reactions occurring in the intracellular medium. This is in contrast to the initial work done on chemical kinetics, which was in simplified systems where reactants were in a relatively dilute, pH-buffered, aqueous solution. In more complex environments, where bound particles may be prevented from disassociation by their surroundings, or diffusion is slow or anomalous, the model of mass action does not always describe the behavior of the reaction kinetics accurately. Several attempts have been made to modify the mass action model, but consensus has yet to be reached. Popular modifications replace the rate constants with functions of time and concentration. As an alternative to these mathematical constructs, one school of thought is that the mass action model can be valid in intracellular environments under certain conditions, but with different rates than would be found in a dilute, simple environment.

      Luckily, as a mathematician, I’m comfortable taking some equations that are only approximately true in the real world and spending a lot of time thinking about them. Reality is always more complicated than we understand, but working with simple approximations tends to be useful—and it’s certainly fun.

      Reply
  5. Scott Hotton says:
    5 July, 2018 at 1:21 am

    This has been a stimulating series for me. I found the rate equation

    \left( \begin{array}{c} \dot{[X]} \\ \dot{[Y]} \\ \dot{[XY]} \end{array} \right) = (\alpha_{\leftarrow}[XY] - \alpha_{\rightarrow} [X][Y]) \left( \begin{array}{r} 1 \\ 1 \\ -1 \end{array} \right)

    for reaction (3)

    X + Y \rightleftharpoons XY

    in part 2 to be interesting because it has aspects like a conservative system and like a dissipative system.

    Here are two independently conserved quantities:

    M = [X] + [Y] + 2 [XY]
    N = [X] - [Y]

    The level surfaces of M form one collection of parallel planes and the level surfaces of N form another collection of parallel planes. The intersection of any level surface of M with any level surface of N is a line parallel to the vector (1,1,-1). I will simply call these the “invariant lines” of the rate equation. Every orbit is contained within an invariant line. In this respect the rate equation is like a conservative system.

    The set of equilibria for the rate equation is a hyperboloid of one sheet. I will designate it by E. It has equation:

    \alpha_{\leftarrow} [XY] - \alpha_{\rightarrow}[X][Y] = 0

    The union of the invariant lines that intersects E tangentially is a parabolic cylinder. I will designate it by P. It has equation:

    (\alpha_{\leftarrow} + \alpha_{\rightarrow} ([X]+[Y]))^2 \; + \; 4 \, \alpha_{\rightarrow}(\alpha_{\leftarrow} [XY] - \alpha_{\rightarrow}[X][Y]) = 0

    The invariant lines “below” P do not intersect E.

    The invariant lines “above” P intersect E in two
    points.

    E \cap P is a parabola contained in the vertical plane with
    equation:

    \alpha_{\leftarrow} + \alpha_{\rightarrow} ([X]+[Y]) = 0

    I will call the intersection of E with

    \{\;([X],[Y],[XY]) \; \colon \; \alpha_{\leftarrow} + \alpha_{\rightarrow} ([X]+[Y]) \leq 0 \; \}

    E_-. It is the portion of E that is on or “behind” the vertical intersection plane. The complement of E_- in E is E_+. It is the portion of E that is “in front” of the vertical plane. Similarly with P_- and P_+.

    Every invariant line above P intersects E_-.
    Let ([X]_{E_-},[Y]_{E_-},[XY])_{E_-} \in E_- and set

    \lambda_{E_-} = -(\alpha_{\leftarrow} + \alpha_{\rightarrow} ([X]_{E_-}+[Y]_{E_-})) > 0

    If the function \lambda(t) satisfies the logistic ODE:

    \dot{\lambda} = \alpha_{\rightarrow} \lambda(\lambda_{E_-} - \lambda)

    then

    \left( \begin{array}{c} {[X]} \\ {[Y]} \\ {[XY]} \end{array} \right) = \left( \begin{array}{c} {[X]_{E_-}} \\ {[Y]_{E_-}} \\ {[XY]_{E_-}} \end{array} \right) + \lambda(t) \left( \begin{array}{r} 1 \\ 1 \\ -1 \end{array} \right)

    satisfies the rate equation above. \lambda=0 corresponds to the point in E_-. \lambda=\lambda_{E_-} corresponds to a point in E_+. For \lambda \in [0,\lambda_{E_-}] the corresponding point is above E_+ and converges towards it. For \lambda > \lambda_{E_-} the corresponding point is below E_+ and converges towards it.

    However the logistic ODE only has solutions for all future time for non-negative initial conditions. Non-negative values of \lambda correspond to points on or above E_- \cup P_+.

    The restriction of the rate equation to an invariant line below P produces a quadratic ODE that does not have any solutions that do not blow up in finite time. Solutions to the rate equation exist for all future time only for initial conditions on or above E_- \cup P_+. Call this S.

    E_+ separates the interior of S into two connected components. Every point in the interior of S above E_+ converges down to E_+. Every point in the interior of S below E_+ converges up to E_+. E_+ is globally attracting on the interior of S (which contains the positive octant). In this respect the rate equation is like a dissipative system.

    It sounds like more complicated chemical reactions can also have aspects like a conservative system and like a dissipative system. Perhaps this example could have been better expressed in terms of reaction vectors and “stoichiometric subspaces” introduced in Part 6.

    Reply
    • John Baez says:
      5 July, 2018 at 1:47 am

      Thanks for your long and interesting comment! It will take me a while to absorb it. For now a trivial remark. You mention two conserved quantities:

      M = [X] + [Y] + 2 [XY]
      N = [X] - [Y]

      For a second I found these surprising, but they’re linear combinations of conserved quantities I could have more easily guessed from the chemistry. If we imagine that an XY molecule is made of an X atom and a Y atom stuck together, “conservation of the total number of X atoms” gives the conserved quantity

      [X] + [XY]

      while “conservation of the total number of Y atoms” gives the conserved quantity

      [Y] + [XY]

      Your quantities are the sum and difference of these. I’m hoping to discover why yours are better in some way.

      Reply
      • Scott Hotton says:
        5 July, 2018 at 5:14 am

        I choose the letter ‘M’ for “mass”. I was thinking of Lavoisier’s law on the conservation of mass. The letter ‘N’ is just the next letter in the alphabet. The level sets of M are inclined planes and I imagined the state of the system moving a little bit like a ball rolling down an inclined plane and remaining within a vertical plane.

        Reply
      • John Baez says:
        5 July, 2018 at 6:14 am

        Hi! I worked through your analysis. It’s nice. In the end this system behaves quite similarly to the reaction network

        \mathrm{B} \mathrel{\substack{\alpha \\\longleftrightarrow\\ \beta}} \mathrm{A} + \mathrm{A}

        which I discussed here:

        • Network theory (part 18).

        Here there are only 2 species, not 3, and only 1 conservation law, not 2, but when you take advantage of the conservation law you’re led to the logistic equation again. There should be some slick way to say that when you take advantage of the conservation laws these systems become isomorphic.

        Reply
        • Scott Hotton says:
          5 July, 2018 at 9:06 am

          I am glad you liked it. There were
          some new ideas for me in that link. If
          i understand correctly the reaction
          vectors for the reaction:

          X + Y \rightleftharpoons XY

          are (1,1,-1) and (-1,-1,1) which
          span a 1 dimensional subspace which is the
          stoichiometric subspace for this
          reaction. What i called the invariant
          lines are parallel to the stoichiometric
          subspace.

          I think you may be right that there
          is a slick way to relate the two systems
          with conservation laws. Since the plane
          N=0 is invariant under the
          dynamical system we can restrict the
          dynamical system to this plane. In this
          plane [X] = [Y]. Even though
          X and Y are two distinct
          species within this plane their
          concentrations could vary the same as
          for the two atoms and diatomic molecule
          reaction, assuming they have the same
          values for their rate constants.

  6. John Baez says:
    5 July, 2018 at 2:16 am

    Any takers on my puzzle? Is it too hard or too easy?

    Reply
    • christopher upshaw says:
      9 July, 2018 at 5:40 am

      Would A <=> B + B <=> C + C + C work?

      Reply
      • John Baez says:
        9 July, 2018 at 3:48 pm

        Yes! If the rate constants are all 1, this gives equilibria with

        [\mathrm{A}] = [\mathrm{B}]^2 = [\mathrm{C}]^3

        The twisted cubic!

        Now some algebraic geometer should tell me what’s cool about the twisted cubic… some facts that might inspire new thoughts on chemistry.

        Reply
  7. Ammar Husain says:
    5 July, 2018 at 7:58 pm

    If you regarding the concentrations as complex variables, I wonder what would the class in K_0 (Var_\mathbb{C}) look like with all its cohomological information contained within. Then take norm on both sides of all the equations. My guess would be that the result is often just a disjoint union of very few points since it seems like it’s very hard to be in detailed balance with so many reactions in the network. So computing or bounding h^0 might be useful.

    Reply
    • John Baez says:
      5 July, 2018 at 10:13 pm

      I want to do something cool like this. But aren’t affine toric varieties too boring to have interesting cohomological invariants? I don’t know, this is just a guess based on the fact that over \mathbb{R} the examples I’ve seen all look contractible… and the fact that people who talk about toric varieties instantly move to more general ones.

      I don’t agree that one only gets a few points. The dimension of your toric variety will be (at least) the number of conserved quantities, since there will be (at least, but it turns out exactly) as many detailed balanced equilibria as there are choices for the values of all those conserved quantities.

      So, for example, in the urea cycle discussed earlier, the space of detailed balanced equilibria will be an 11-dimensional toric variety. As far as chemistry goes this variety is defined over the rig of nonnegative reals, but we can work over the complex numbers if we’re looking for interesting invariants to attach to this variety.

      But because of everything I’ve said, I have a fear that this toric varieity, like the space of possible choices of conserved quantities, will be contractible!

      That doesn’t rule out the possibility of interesting invariants, but it means they’ll be more like algebraic geomety than topology.

      Reply
      • Ammar Husain says:
        5 July, 2018 at 10:55 pm

        Ah. It is more typical to have some conserved quantities.

        Reply
        • Ammar Husain says:
          18 July, 2018 at 8:35 pm

          Which basic operations do you want on the reaction networks. Adding a reaction, but keeping the possible reactants the same and adding a reactant that doesn’t contribute to any reactions are the ones that come to mind. But would you want anything different? To determine the types of morphisms for the resulting varieties.

        • John Baez says:
          20 July, 2018 at 2:56 pm

          I haven’t thought about this! It’s an interesting question. There is a category of reaction networks, or more precisely stochastic Petri nets, where a map is a function sending species to species and reactions to reactions, preserving rate constants. There should be a similar category of detailed balanced reaction networks, and a functor from this to a category of affine toric varieties. Working out the details should tell us what kind of morphisms of toric varieties we want!

          First I should finish my paper with Jade Master, which among other things discusses the category of Petri nets.

You can use Markdown or HTML in your comments. You can also use LaTeX, like this: $latex E = m c^2 $. The word 'latex' comes right after the first dollar sign, with a space after it.

Fill in your details below or click an icon to log in:

Gravatar
WordPress.com Logo

You are commenting using your WordPress.com account. ( Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. ( Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. ( Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. ( Log Out /  Change )

Cancel

Connecting to %s

This site uses Akismet to reduce spam. Learn how your comment data is processed.