Question

- Why do some animal populations like lynx appear to evolve through boom and bust cycles?

Topics

- Law of mass action and the birth of physical chemistry at the end of the 1800’s.

Further reading

150 Years of Mass Action Law by E. O. Voit, H. A. Martens, S. W. Omholt, PLOS Computational Biology, 2015.

A Defence of Beanbag Genetics by J. B. S. Haldane in Perspectives in Biology and Medicine, Volume 7, Number 3, Spring 1964, pp. 343-360.

Wikipedia page on law of mass action.

The linear compartmental differential-equation models we have studied so far are very useful. We’ll see more of this whe we study consider telephone networks. But they can not represent the effects of interactions. For many natural systems, particularly in chemistry and biology, the most interesting features are driven by interactions where the system is changed when different moving pieces come together. If you take a moment to think, you can probably come up with a bunch of examples of important interactions.

The puzzle of find ways to describe rates of change in response particle interactions was confronted from the late 1700’s through the 1800’s during the maturation of the field of chemistry from its alchemical roots into a modern science. Following the successes of Newton’s laws and the force of gravity at explaining the motion of the planets, people began to search for the “forces” that would similarly explain the dynamics of chemical processes. The initial efforts along this line failed. The concept of “force”, *e.g.* the rate of change of momentum of the system, is of much less importance in chemistry. An answer was finally discovered and published in 1864 by mathematician Cato Guldberg and his chemist brother-in-law Peter Waage. But the value of their discovery was not recognized until Jacobus van’t Hoff’s rediscovery in 1877. Today, we call their discovery the “Law of Mass Action”.

Let’s imagine a well-mixed bucket of water into which we have added a few drops of blue and yellow dye particles. Initially, there are no green particles, but when a blue and a yellow particle bump into each other, they occasionally fuse to become a green particle. Let \(b(t)\), \(y(t)\), and \(g(t)\) represent the concentrations of these blue, yellow, and green particles, respectively, over time. The we can summarize the creation of green particles as an irreversible reaction \(b + y \rightarrow g\). We call the blue and yellow particles the reactants, and the green particles the product. Aside from this reaction, we will assume nothing creates or removes any of these particles.

Now, let’s hypothesize that, since change only occurs when blue and yellow particles come together, that the rate of change in the number each particle depends only on the current concentrations of the blue and yellow particles. Then the rates of change of the 3 concentrations have the form \[\begin{align} \frac{db}{d t} &= - U(b,y),\\ \frac{dy}{d t} &= - V(b,y),\\ \frac{dg}{d t} &= W(b,y), \end{align}\] for some unknown functions \(U\), \(V\), and \(W\). But we can simplify this. Since blue particles are not added or removed, and only change color to green when fused with a yellow particle, the sum of the numbers of blue and green particles at any time must always be the same as the initial number of blue particles. Thus, \(b(t) + g(t) = b(0)\). If we differentiate this “conservation law”, we find \[\begin{gather} \frac{db}{d t} = - \frac{dg}{d t}. \end{gather}\] The same kind of observation holds for yellow and green particles; \(y(t) + g(t) = y(0)\), from which we deduce that \[\begin{gather} \frac{dy}{d t} = - \frac{dg}{d t}. \end{gather}\] Thus, by substitution, \[\begin{align} \frac{db}{d t} &= - W(b,y), \\ \frac{dy}{d t} &= - W(b,y), \\ \frac{dg}{d t} &= W(b,y). \end{align}\]

We can now ponder what the rate function \(W(b,y)\) looks like. We expect the rate to be increasing – the more particles are in the bucket, the faster they should fuse. And we expect this increase to be smooth – small changes in concentration should not cause large changes in the rate. In this case, one way to represent the rate may be it’s Maclaurin series under our known constraints. The Maclaurin expansion up to third order in \(b\) and \(y\) is \[W(b,y) = c_{00} + c_{10} b + c_{01} y + c_{20} b^2 + c_{11} by + c_{02} y^2 + \ldots .\] Now, when either of the reactants is exhausted, the reaction can proceed no further. Put in terms of numbers and applied to our Taylor series, when \(b = 0\), we know \[0 = W(b,y) = c_{00} + c_{01} y + c_{02} y^2 + O(y^3).\] From this equation, we deduce that obvious solution \(0 = c_{00} = c_{01} = c_{02}\). In the same way, when \(y = 0\), we know \[0 = W(b,0) = c_{00} + c_{10} b + c_{20} b^2 + O(b^3),\] which implies \(0 = c_{00} = c_{10} = c_{20}\). Thus, \[W(b,y) = c _ {11} by + O(b^2,y^2)\] If we neglect higher order terms, we recover what is commonly known at the **law of mass action**.

**Law of mass action**: In a strongly-mixed system, the rate of a reaction with two reactants (\(b + y \rightarrow g\)) is proportional to the product of the concentrations of those two reactants. \[W(b,y) \propto b y.\]

Applying the law of mass-action to our simple system, \[\begin{align} \frac{db}{d t} &= - \alpha b y \\ \frac{dy}{d t} &= - \alpha b y \\ \frac{dg}{d t} &= \alpha b y \end{align}\] Using the conservation laws stated above, we can reduce our system of 3 equations to a single differential equation for \(g(t)\) and two scalar equations. \[\begin{align} b(t) &= b(0) - g(t), \\ y(t) &= y(0) - g(t), \\ \frac{dg}{d t} &= \alpha (b(0) - g) (y(0) - g). \end{align}\]

The use of the Maclaurin series is not well-justified, however. Maclaurin series are intended for functions that are smooth on both sides of zero. Our reaction rate function \(W(b,y)\) is not defined for negative values. More elaborate derivations can be given, but real reason the Law of Mass Action is so important is that it explains the data. The following two examples will show that it explains both dynamic and equilibrium observations of some chemical systems.

For example, in 1933, Grant and Hinshelwood published measurements of the rate of progress of a chemical reaction of potassium hydroxide and alkyl bromide, \[KOH + C_2 H_5 Br \rightarrow \text{Stuff}.\]

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Time (minutes), Acid titer
#
# @article{bib:Grant1933,
# title={71. The kinetics of reactions in solution. The interaction of potassium hydroxide and the alkyl halides in ethyl alcohol},
# ISSN={0368-1769},
# url={http://dx.doi.org/10.1039/JR9330000258},
# DOI={10.1039/jr9330000258},
# journal={Journal of the Chemical Society (Resumed)},
# publisher={Royal Society of Chemistry (RSC)},
# author={Grant, G. H. and Hinshelwood, C. N.},
# year={1933},
# pages={258},
# annote={Table II, b},
# }
#
# KOH + C2H5Br -> stuff
# Initial concentrations, 0.0125 M KOH, 0.050M C2H5Br
#
0, 10.
12, 8.59
24, 7.33
41, 6.09
63, 4.71
84, 3.74
108, 2.90
131, 2.25
148, 1.93
171, 1.45
218, 0.96
```

Suppose \(b(t)\) is concentration of potassium hydroxide(\(b(t) = [KOH]\)), \(y(t)\) is concentration of alkyl bromide(\(y(t) = [C_2 H_5 Br]\)), and \(g(t)\) is the concentration of one of the products. In order check the agreement between their observations and the predictions of the law of mass action, we must calculate \(g(t)\) when \(g(0) = 0\). While we can not completely solve for \(g(t)\) analytically, we can derive a testable formula. Starting from the equation \[\begin{align} \frac{dg}{d t} &= \alpha (b(0) - g) (y(0) - g), \end{align}\] we can separate variables and integrate, so \[\begin{align} \frac{1}{b(0) - y(0)} \ln\left(\frac{b(0) - g(t)}{y(0) - g(t)}\right) &= \alpha t + C \end{align}\] The left-hand side can be plotted based on Grant and Hinshelwood’s data, and shows a straight line. Thus, the law of mass action provides a good description of some chemical reactions.

[Show code]

```
#!/usr/bin/env python3
from pylab import *
from scipy import linalg
fname = 'GrantHinshelwood1933.csv'
data = loadtxt(fname, delimiter=',')
B0 = 0.050
Y0 = 0.0125
def proper():
t = data[:,0]*60 # convert minutes to seconds
g = Y0 - data[:,1]*Y0/10. # renormalized collected data
rhs = log( (B0-g)/(Y0-g) )/(B0-Y0)
A = vstack([t**0, t**1]).T
ans = linalg.solve(A.T.dot(A), A.T.dot(rhs))
plot(t, rhs, 'o',label='Data')
plot(t, ans[0]+ans[1]*t, 'r-',label='%.3f+%.4f t'%(ans[0],ans[1]))
title("Transformed concentration $g(t)$ over time", fontsize=18)
xlabel("Time (seconds)", fontsize=18)
ylabel(r"$\frac{1}{b(0)-y(0)}\left[ \ln(1-\frac{g(t)}{b(0)}) - \ln(1-\frac{g(t)}{y(0)})\right]$", fontsize=18)
legend(framealpha=0.0)
tight_layout()
savefig('GrantHinshelwood1933.png', transparent=True)
savefig('GrantHinshelwood1933.pdf')
proper()
```

A second context where we can check the validity of the law of mass action is in the establishment of chemical equilibrium. According to the law of mass action, a reaction is in equilibrium when its forward rate is equal to its backward rate.

\[b + y \; \underset{k _ -}{\overset{k _ +}{\rightleftharpoons}} \; g\]

The corresponding system of differential equations is given by \[\begin{align} \frac{db}{d t} &= - k_{+} b y + k_{-} g \\ \frac{dy}{d t} &= - k_{+} b y + k_{-} g \\ \frac{dg}{d t} &= k_{+} b y - k_{-} g \end{align}\] And so, at equilibrium, when none of the concentrations are changing, we expect \[\frac{b_{EQ} y_{EQ}}{ g_{EQ} } = \frac{k_{-}}{k_{+}}\] or, when expressed in terms of initial concentrations, \[k_{+} (b(0) + g(0) - g_{EQ}) (y(0) + g(0) - g_{EQ}) = k_{-} g_{EQ}.\]

A standard example of this is found measurements of the chemical equilibrium of iron thiocyanate (see Frank and Oswalt, 1947). In solution, iron thiocyanate creates a red coloring, the intensity of which can be measured using now-standard spectrometry techniques. The intensity of light \(T _ {\lambda}\) at wavelength \(\lambda\) passing through a liquid decays exponentially with the distance traveled \(\ell\) at a rate that depends on the concentration \(g _ {EQ}\) and the extinction coefficient \(\beta _ {\lambda}\) according to \[T _ {\lambda} = e^{-\beta _ {\lambda} g _ {EQ} \ell}\]

The extinction coefficient, equilibrium constant, and equilibrium concentration are not known a-priori, but can be fit from the data, and once done, we find very good agreement with the equilibrium concentration predicted by the law of mass action.

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
# Iron concentration (M), thiocyanate concentration (M), frac transmission at 4500 angstroms
#
# Data from Frank and Oswalt, 1947, Table 1
# http://dx.doi.org/10.1021/ja01198a026
#
0.001, 0.0003, 0.615
0.002, 0.0003, 0.399
0.002, 0.0003, 0.425
0.003, 0.0003, 0.306
0.003, 0.0003, 0.313
0.005, 0.0003, 0.190
0.008, 0.0003, 0.113
0.008, 0.0003, 0.119
```

[Show code]

```
from pylab import *
data = loadtxt('FrankOswalt.csv',delimiter=',')
a, b, T = data[:,0], data[:,1], data[:,2]
l = 1.305 # estimates transmission path length
@vectorize
def xeq(a,b,K):
vls = roots([1, -(a + b + 1/K), a*b])
return min(vls)
def main():
figure(1)
X = a+b
D = -log10(T)/l
Y = a*b/D
pcs = polyfit(X,Y,1)
beta = 1/pcs[0]
K = 1/(pcs[1]*beta)
plot(X,Y,'x')
plot(X,pcs[0]*X + pcs[1],'r-')
xlabel('a+b',fontsize=18)
ylabel('ab/D',fontsize=18)
title('Recreation of Frank Fig 1',fontsize=18)
tight_layout()
print("Lst Sqrs Fit: ", pcs, [beta, K, 1/K])
A = linspace(min(a),max(a),50)
x_prd = xeq(A, b[0], K)
x_obs = D/beta
#figure(2,figsize=(6,9))
#subplot(2,1,1)
#semilogy(a,T,'o')
#ylabel('Transmission at 4500 Ang', fontsize=18)
#subplot(2,1,2)
figure(2)
plot(A,x_prd,'r-',label='Mass Action Prediction')
plot(a,x_obs,'go',label='Observations')
ylabel('Equilibrium FeSCN++ (M)', fontsize=18)
xlabel('Initial Iron Concentration (M)', fontsize=18)
legend(framealpha=0.0, fontsize=12)
text(0.005,0.00008, 'K = %.2f'%K, fontsize=12)
text(0.005,0.00006, r'$\beta$ = %.2f'%beta, fontsize=12)
tight_layout()
savefig('massActionEqtest.png', transparent=True)
savefig('massActionEqtest.pdf')
show()
main()
```

Thus, the simple law of mass action is at least sometimes a good explanation of both the rates and equilibrium behavior of some chemical reactions and perhaps other interactions that can be modelled as a binary reaction.

The law of mass action is the most commonly used rate law for binary reactions, but it not the only one. There are a variety of functions that satisfy our basic constraints on our rate function and have been used for modelling. Other rate functions that have been used include \[ b^m y^n \] and the Michaelis-Menton form \[\frac{k b y}{1 + j_b b + j_y y}\] Another common form uses the Hill function, \[\frac{k b y^n}{h^n + y^n}.\] The Hill function allows step-like switching of a reaction rate when the exponent \(n\) is large.

All of these laws assume strong-mixing, so that the physical positions of particles are not important. In weakly mixed systems, the concentrations of reactants and products may very significantly between spatial locations, possibly leading to the creation of gradients and waves.

There are also rate laws that allow for 3 or more reactants at once. However, while mathematically convenient, these seldom provide good descriptions of a real reactions. If you think about different kinds of molecules bouncing around in a bottle of gas, the strong-mixing principle assumes that binary reactions, ones that occur when two molecules collide, are rare. For a reaction with 3 reactants to occur, all 3 kinds of molecule would have to collide at the same time, which would be much more rare than two molecules colliding. So, usually, reactions described with 3 or more reactants are actually short-hand for a series of binary molecular reactions, and one of these binary reactions controls the rate of the over-all reaction.

- Below is a data set on bromine’s reaction with xylene. If we hypothesize that the reaction rate has the form \[\frac{dy}{dt} = - k y^m\] where \(y\) is the concentration of bromine, \(k\) is the rate constant, and \(m\) is the reaction order, what do you estimate the reaction’s order to be?

[ * Data : hide , shown as table , shown as CSV shown as QR * ]

```
#Time (minutes), Bromine concentration (moles/liter)
#
# At 17 degrees Celcius, a small amount of iodine and bromine
# are added to pure liquid xylene. Iodine and xylene levels
# are approximately constant, so changes in the rate must
# be due to changes in bromine concentration.
#
# Data collected by A. Neyens, published in Hill, 1977, page 45
#
0 , 0.3335
2.25 , 0.2965
4.50 , 0.2660
6.33 , 0.2450
8.00 , 0.2255
10.25, 0.2050
12.00, 0.1910
13.50, 0.1794
15.60, 0.1632
17.85, 0.1500
19.60, 0.1429
27.00, 0.1160
30.00, 0.1053
38.00, 0.0830
41.00, 0.0767
45.00, 0.0705
47.00, 0.0678
57.00, 0.0553
63.00, 0.0482
```