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.
- Lotka-Volterra predator-prey equations inspired by Lotka parallelling physical chemistry with physical biology.

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.

Volterra and Lotka’s exchange on the theory of predator-prey interactions in Nature, 1926-7.

- What are some other examples of important interactions where two things have to come together, and when they do, things change?

Folk-knowledge regarding the natural laws of population change can be traced back to the beginnings of recorded history, but modern western thought on the topic is conventionally credit to Thomas Malthus’s “An Essay on the Principle of Population”

The simplest model of population growth is Thomas Malthus’s. Malthus proposed that the rate of proliferation of a population is proportional to the number of individuals in that population.

Take a look at Huffaker’s mite experiment, as see their observations for the numbers of prey mites and predator mites over time.

- Law of mass action, Guldberg, Waage, and others
Wilhelm Ostwald and the birth of physical chemistry (Lotka was his student)

- requirements (strong mixing, no spatial structure importance)
- ideal (reaction network sufficies for description)
consequence – reaction rate equation

Let us consider an irreversible reaction \(x + y \rightarrow z\).

\[\begin{align} \frac{dx}{d t} &= - f(x,y),\\ \frac{dy}{d t} &= - g(x,y),\\ \frac{dz}{d t} &= h(x,y). \end{align}\] The system has two conservation laws at play: \(x + z\) and \(y + z\). To see that the total amount of \(x + z\) never changes. We can differentiate and use the system above to simplify the result. \[\begin{gather} \frac{d}{d t} \left(x + z \right)= \frac{dx}{d t} + \frac{dz}{d t} = - f(x,y) + h(x,y) = 0 \end{gather}\] which implies \(f = h\). By conservation of mass of type y, we also have \(g=h\), so we need only consider \[\begin{align} \frac{dx}{d t} &= - h(x,y), \\ \frac{dy}{d t} &= - h(x,y), \\ \frac{dz}{d t} &= h(x,y). \end{align}\]

One way to determine \(h(x,y)\) is to look at it’s McLaurin series under our known constraints. The McLaurin expansion up to third order in \(x\) and \(y\) is \[h(x,y) = c_{00} + c_{10} x + c_{01} y + c_{20} x^2 + c_{11} xy + 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 \(x = 0\), we know \[0 = h(0,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 = h(x,0) = c_{00} + c_{10} x + c_{20} x^2 + O(x^3),\] which implies \(0 = c_{00} = c_{10} = c_{20}\). Thus, \[h(x,y) = c _ {11} xy + O(x^2,y^2)\] If we neglect higher order terms, we recover what is commonly known at the **law of mass action**.–

**Basic Law of mass action**: The rate of a reaction with two reactants (\(x + y \rightarrow z\)) is proportional to the product of the concentrations of those two reactants. (\(h(x,y) \propto x y.\))

Applying the basic law of mass-action to our simple system, \[\begin{align} \frac{dx}{d t} &= - k x y \\ \frac{dy}{d t} &= - k x y \\ \frac{dz}{d t} &= k x y \end{align}\] Since the the sum \(x + z\) is conserved, it must always be the same as whatever it started: \(x(t) + z(t) = x(0) + z(0)\). And since the sum \(y + z\) is also conserved, \(y(t) + z(t) = y(0) + z(0)\). Using these conservation laws, we can reduce our system of 3 equations to a single differential equation for \(z(t)\) and two scalar equations. \[\begin{align} x(t) &= x(0) + z(0) - z(t), \\ y(t) &= y(0) + z(0) - z(t), \\ \frac{dz}{d t} &= k (x(0) + z(0) - z) (y(0) + z(0) - z). \end{align}\]

In 1933, Grant and Hinshelwood published measurements of the rate of progress of a chemical reaction of potassium hydroxide and alkyl halides.

\[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
```

In order check the agreement between their observations and the predictions of the law of mass action, we must calculate \(z(t)\) when \(z(0) = 0\).

\[\begin{align} \frac{dz}{d t} &= k (x(0) - z) (y(0) - z), \\ \frac{dz}{(x(0) - z) (y(0) - z)} &= k dt \\ \frac{1}{x(0) - y(0)} \ln\left(\frac{x(0) - z(t)}{y(0) - z(t)}\right) &= k t + C \end{align}\] Using the initial condition \(z(0)= 0\) to determine the integration constant, we find \[\begin{align} \frac{1}{x(0) - y(0)} \ln\left(\frac{1 - z(t)/x(0)}{1 - z(t)/y(0)}\right) &= k t . \end{align}\]

This can be plotted and compared against Grant and Hinshelwood’s data, and shows the law of mass action can provide 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=',')
A0 = 0.0125
B0 = 0.050
def proper():
t = data[:,0]*60
u = data[:,1]*A0/10. # renormalized collected data
b = log((A0/B0)*((B0 - (A0-u))/u) )/(B0-A0)
A = vstack([t**0, t**1]).T
ans = linalg.solve(A.T.dot(A), A.T.dot(b))
plot(t, b, 'o',label='Data')
plot(t, ans[0]+ans[1]*t, 'r-',label='%.3f+%.4f t'%(ans[0],ans[1]))
title("Transformed concentration $z(t)$ over time", fontsize=18)
xlabel("Time (seconds)", fontsize=18)
ylabel(r"$\frac{1}{y(0)-x(0)}\left[ \ln(1-\frac{z(t)}{y(0)}) - \ln(1-\frac{z(t)}{x(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.

\[x + y \; \underset{k _ -}{\overset{k _ +}{\rightleftharpoons}} \; z\]

The corresponding compartmental system of differential equations is given by \[\begin{align} \frac{dx}{d t} &= - k_{+} x y + k_{-} z \\ \frac{dy}{d t} &= - k_{+} x y + k_{-} z \\ \frac{dz}{d t} &= k_{+} x y - k_{-} z \end{align}\] And so, at equilibrium, when none of the concentrations are changing, we expect \[\frac{x_{EQ} y_{EQ}}{ z_{EQ} } = \frac{k_{-}}{k_{+}}\] or, when expressed in terms of initial concentrations, \[0 = (x(0) + z(0) - z_{EQ}) (y(0) + z(0) - z_{EQ}) - \frac{k_{-}}{k_{+}} z_{EQ}.\]

A standard exammple 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 \(z _ {EQ}\) and the extinction coefficient \(\beta _ {\lambda}\) according to \[T _ {\lambda} = e^{-\beta _ {\lambda} z _ {EQ} \ell}\]

The extinction coefficent, 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)
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')
main()
```

Thus, the simple law of mass action is atleast 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 simple 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 constrainst on our rate function and have been used for modelling. Other rate functions that have been used include \[ x^m y^n \] and the Michaelis-Menton form \[\frac{k x y}{1 + j_x x + j_y y}\] Another common form uses the Hill function, \[\frac{k x y^n}{h^n + y^n}\] which allow’s for step-like switch of a reaction rate when the exponent \(n\) is large.

The law of mass-action has a more general form as well that sometimes appears in chemistry and mathematics and can be applied to reactions with 3 or more reactants.

**General Law of mass action**: The rate of a reaction with any number of reactants should be proportional to the product of the concentrations of all the reactants. If \(a\) is a reactant species, \([a]\) is the concentration of species \(a\), and \(n _ a\) is the stochiometry coefficient of species \(a\), then, given a reaction of the form \[ \left( \sum _ a n _ a a \right) \rightarrow \left( \sum _ z m _ z z \right),\] the rate of change in the concentration of the products \(z\) is given by \[\frac{d[z]}{dt} = k m _ z \prod _ a [a]^{n _ a} \] while the rate of reduction of the concentrations of the individual reactants obeys \[\frac{d[a]}{dt} = - k n _ a \prod _ a [a]^{n _ a}.\]

However, while mathematically convenient, this general form of the law of mass action rarely provides an accurate description of a real reaction. If you think about different kinds of molecules bouncing around in a bottle of gas, the strong-mixing principle implies 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 even 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.

- Lotka-Volterra predator-prey
- equations
- phase-plane analysis matching Huffaker’s
- Model flaws

Let \(N(t)\) be the number of prey at time \(t\) while \(P(t)\) is the number of predators at time \(t\). The simplest model that captures the interactions needs only 3 types of events. Prey have to be able to reproduce (or else predators will surely run out of food), predators have to be able to consume prey to make more predators, and predators have to die (or else they are sure to drive the prey to extinction). Thus, the simpliest predator-prey reaction network is \[\begin{gather} N \xrightarrow{r} 2 N , \quad P \xrightarrow{d} \emptyset , \quad N + P \xrightarrow{c} 2 P \end{gather}\] When we convert these reactions into a differential equation system for the rates of change in the numbers of prey and predators, we find \[\begin{align*} \frac{dN}{d t} &= r N - c P N, \\ \frac{dP}{d t} &= c N P - d P. \end{align*}\]

Here is a python script showing how to calculate time-series solutions based on the Lotka–Volterra predator-prey differential equations we derived in class.

[Show code]

```
#!/usr/bin/env python3
from scipy import array, linspace
from scipy import integrate
from matplotlib.pyplot import *
def vector_field(X, t, r, c, d):
# The differenti equations are
#
# dN
# -- = -c N P + r N
# dt
#
# dP
# -- = c N P - d P
# dt
N = X[0] # prey density
P = X[1] # predator density
return array([r*N - c*N*P, c*N*P - d*P])
def main():
# Parameters for our predator-prey differential equations
r = 1. # net density-independent growth rate of prey population
c = 2.
d = 1.
# set up our initial conditions
N0 = 1.
P0 = 1.
X0 = array([N0, P0])
# choose the time's we'd like to know the approximate solution
t = linspace(0., 20., 100)
X = integrate.odeint(vector_field, X0, t, args=(r,c,d))
# now, plot the solution curves
figure(1)
subplot(2,1,1) # a time-series plot
plot(t, X[:,0], 'bx-', linewidth=2, label='Prey')
plot(t, X[:,1], 'g+-', linewidth=2, label='Predator')
xlabel('Time (days)')
ylabel('Number')
legend(framealpha=0)
subplot(2,1,2) # a phase-space plot
plot(X[:,0], X[:,1], 'k-')
xlabel('Number of Prey')
ylabel('Number of Predators')
tight_layout()
savefig('predatorprey.pdf')
savefig('predatorprey.png', transparent=True)
show()
main()
```

- Alternative rate laws
- handling times
- power-law
- Gompertzian
- Hill functions …

- Volterra and DAcona’s original predator-prey fish data from Trieste has been extended as far as 1961. Plot this data, and attempt to fit it was the Lotka-Volterra model.

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

```
# year, selachian tons, total tons
#
# Yearly data on Trieste market selachians and totals
# collected by DAncona and given to Volterra.
# 1961. selachians, a group at the top of the food web
# or very close ot it. The lower part gives the total catch of selachians, bony fishes, mollusks,
# and crustaceans in 1Oaton units.
#
# @article{bib:Scudo1971,
# title={Vito Volterra and theoretical ecology},
# volume={2},
# ISSN={0040-5809},
# url={http://dx.doi.org/10.1016/0040-5809(71)90002-5},
# DOI={10.1016/0040-5809(71)90002-5},
# number={1},
# journal={Theoretical Population Biology},
# publisher={Elsevier BV},
# author={Scudo, Francesco M.},
# year={1971},
# month={Mar},
# pages={1-23},
# }
1904,505,5249
1905,524,4935
1906,470,4823
1907,470,5469
1908,800,9834
1909,457,6472
1910,369,6663
1911,624,7224
1912,749,7925
1913,1003,6357
1914,760,5162
1915,69,914
1916,138,846
1917,116,744
1918,NaN,NAN
1919,484,2499
1920,896,5747
1921,853,6449
1922,867,8216
1923,770,7593
1924,704,7987
1925,901,7547
1926,1000,7699
1927,913,8052
1928,990,8049
1929,965,6073
1930,723,8400
1931,716,10018
1932,790,7631
1933,519,7530
1934,582,8642
1935,699,7073
1936,809,9413
1937,764,9940
1938,913,10078
1939,1172,9684,partial
1940,715,7843
1941,196,5331
1942,326,4418
1943,140,2312
1944,49,1361
1945,149,1414
1946,1062,7846
1947,988,9979
1948,1159,11050
1949,934,8468
1950,863,10477
1951,995,9279
1952,866,8648
1953,631,8195
1954,580,8195
1955,754,9452
1956,754,7280
1957,833,7369
1958,826,8487
1959,806,7899
1960,853,7509
1961,844,7513
```

- 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
```