Lecture notes. Written from the material of MIT 8.592J, licensed CC BY-NC-SA 4.0. These are notes, not a transcript: the material has been reorganised and rewritten. This adaptation carries the same licence, and the original is linked above.

1. Sequence Entropy and Evolving Probabilities¶

What this covers¶

This chapter answers two questions from the opening lecture of 8.592J. First: if a biological sequence — a stretch of DNA, a protein — were generated by drawing each letter independently at random, how much does knowing the letter frequencies actually save you, and what is a good way to put a number on that? Second: how do you describe a probability distribution that is changing in time, whether it is the composition of a mutating sequence, the size of a mutant sub-population, or the number of active molecules in a cell? It assumes elementary probability (independence, normalisation), Stirling's approximation, and enough linear algebra to read an eigenvector equation.

Sequences as independent random draws¶

Text is a sequence of characters: English uses an alphabet of \(m = 26\) letters (plus spaces and punctuation), DNA an alphabet of \(m = 4\) bases, protein an alphabet of \(m = 20\) amino acids. Real biological sequences are far from random — there are constraints and correlations at every scale that make them meaningful — but a model that assumes randomness is a useful baseline: comparing it to real sequences is one way to find where the structure is.

The simplest such model treats a sequence of length \(N\) as \(N\) independent, identically distributed (IID) draws from an alphabet \(\alpha = 1, \ldots, m\) with probabilities \(\{p_\alpha\}\), normalised as

\[\sum_{\alpha=1}^m p_\alpha = 1.\]

The probability of a specific sequence \(S = \{\alpha_1, \ldots, \alpha_N\}\) is the product of the probabilities of its letters,

\[p(S \mid \{p_\alpha\}) = \prod_{\ell=1}^N p_{\alpha_\ell}.\]

This probability depends only on how many times each letter occurs — the composition \(\{N_\alpha\}\) — not on the order the letters appear in. The number of ways to reorder a sequence with a given composition is the multinomial coefficient

\[\mathcal{N} = \frac{N!}{\prod_{\alpha=1}^m N_\alpha!},\]

which is exactly the coefficient that appears in the multinomial expansion

\[(p_1 + \cdots + p_m)^N = \sum_{\{N_\alpha\}}{}' \; \frac{N!}{\prod_\alpha N_\alpha!}\, p_1^{N_1}\cdots p_m^{N_m},\]

the sum running over all compositions with \(\sum_\alpha N_\alpha = N\). Both sides equal \(1\), by normalisation, and each term on the right is the multinomial probability of observing composition \(\{N_\alpha\}\),

\[p(N_1, \ldots, N_m) = \frac{N!}{\prod_\alpha N_\alpha!}\, p_1^{N_1} \cdots p_m^{N_m}.\]

Under this model, the most likely composition (the mode) coincides with the average composition,

\[N_\alpha^* = \langle N_\alpha \rangle = p_\alpha N.\]

Counting compositions: where entropy comes from¶

For large \(N\), the overwhelming majority of sequences generated by the model will have this typical composition \(N_\alpha = p_\alpha N\). The number of distinct sequences sharing it, \(\mathcal{N}\) from above, is a natural way to quantify how much "information" a sequence of length \(N\) carries: it counts how many equally-likely messages there are to choose among. A sensible measure of information should scale roughly linearly with \(N\) — a message twice as long should carry about twice the information, in the absence of context or syntax — which suggests taking a logarithm. Using Stirling's approximation \(\log N! \approx N\log N - N\) (valid once every \(N_\alpha \gg 1\)),

\[\log \mathcal{N} = \log N! - \sum_\alpha \log N_\alpha! \approx N\log N - N - \sum_\alpha (N_\alpha \log N_\alpha - N_\alpha) = -N \sum_\alpha \frac{N_\alpha}{N}\log\frac{N_\alpha}{N}.\]

This is exactly the entropy of mixing from thermodynamics. More generally, for any distribution \(\{p_\alpha\}\) define the mixing entropy

\[S(\{p_\alpha\}) = -\sum_\alpha p_\alpha \log p_\alpha.\]

Entropy is usually pictured as a measure of disorder among the possibilities; the information content of picking out one specific outcome from that disorder is tied to this same quantity, made concrete in the example below.

Worked example: how much do you learn from base composition?¶

To send a DNA sequence \(\text{ACTG}\cdots\) down a binary channel with no assumptions about it, you need \(2\) bits per base, since there are \(m=4\) equally likely letters and \(2^2 = 4\). Now suppose prior analysis of an organism's genome tells you the typical base composition is biased, i.e. \(p_\alpha = N_\alpha / N\) known and unequal. Only \(\mathcal{N} = N!/\prod_\alpha N_\alpha! \ll 4^N\) sequences are then "typical" ones, so

\[\log_2 \mathcal{N} = -N\sum_\alpha p_\alpha \log_2 p_\alpha < 2N.\]

Knowing the bias in advance saves you bits: instead of \(2\) bits per base you need only the entropy (in bits) of \(\{p_\alpha\}\). Define the information gained per letter as the shortfall from the naive \(2\)-bit code,

\[\mathcal{I}(\{p_\alpha\}) = 2 - \sum_\alpha p_\alpha \log_2\!\left(\frac{1}{p_\alpha}\right).\]

Two checks: if \(p_\alpha = 1/4\) for all four bases, \(\mathcal{I} = 0\) — you already assumed uniformity, so there is nothing to gain. If instead \(p_A = p_T = 0\) and \(p_C = p_G = 1/2\) (only two letters ever occur, equally often),

\[\mathcal{I} = 2 - \sum_{G,C} \tfrac{1}{2}\log_2 2 = 1 \text{ bit per base} —\]

knowing in advance that only two letters are possible halves the code you need, exactly as it should.

From a fixed distribution to an evolving one¶

Reproduction copies genetic material imperfectly, so a population's genetic composition diverges and evolves over generations. Because the changes are stochastic, the right object to track is not a fixed distribution \(\{p_\alpha\}\) but an evolving one, \(\{p_\alpha(t)\}\).

Mutations, and why some matter more than others. Take the DNA fragment \(\text{ATT CGC ATG}\), transcribed into the complementary mRNA \(\text{UAA GCG UAC}\) and translated into the peptide \(\text{Leu Ala Tyr}\). A single replication error changing the last codon from ATG to ATC turns it into a stop signal: \(\text{Leu Ala STOP}\), and translation halts wherever this happens in the sequence. Such a mutation is deleterious: it can prevent viable offspring. But the genetic code is redundant, so many mutations are synonymous — they change the DNA letter without changing the resulting amino acid, and so have no effect on the organism. Because synonymous variation carries no cost, it accumulates and differs between individuals at otherwise-identical genes: these single-nucleotide polymorphisms (SNPs) are what DNA fingerprinting compares. (Non-synonymous mutations are not automatically deleterious either — they may still produce viable offspring.)

Classical genetics, which predates the molecular picture, already had language for hereditary variation: a gene is a heritable unit, and diploid organisms like humans carry two copies of each, called alleles. Some alleles are dominant, others recessive: if \(A_1\) (brown eyes) is dominant over \(A_2\) (blue eyes), then genotypes \(A_1A_1\) and \(A_1A_2\) both give brown eyes, and only \(A_2A_2\) gives blue. This vocabulary — alleles, genotype counts — is what the population model below is built from.

The transition matrix and the discrete-time master equation¶

Return to the \(N\)-independent-site picture, but now let each site change state from generation to generation. Write \(\pi_{\alpha\beta}\) for the probability that a site in state \(\beta\) moves to state \(\alpha\) in one generation (first index destination, second index source); collectively these form the \(m\times m\) transition probability matrix \(\pi\). Probabilities evolve as

\[p_\alpha(\tau+1) = \sum_{\beta=1}^m \pi_{\alpha\beta}\, p_\beta(\tau), \qquad \text{or, in matrix form,} \qquad \vec p(\tau+1) = \pi\, \vec p(\tau) = \pi^\tau\, \vec p(1),\]

the last equality following by iterating, provided \(\pi\) itself does not change over time.

Normalisation forces every column of \(\pi\) to sum to one — a site ends up somewhere — so

\[\sum_\alpha \pi_{\alpha\beta} = 1, \qquad \text{equivalently} \qquad \pi_{\beta\beta} = 1 - \sum_{\alpha\neq\beta}\pi_{\alpha\beta}:\]

the probability of staying put is the complement of the probability of leaving. Substituting this back gives a flux-balance form of the same recursion,

\[p_\alpha(\tau+1) = p_\alpha(\tau) + \sum_{\beta\neq\alpha}\bigl[\pi_{\alpha\beta}p_\beta(\tau) - \pi_{\beta\alpha}p_\alpha(\tau)\bigr]:\]

the change at \(\alpha\) is what flows in from every other state minus what flows out to every other state. (Without the independence assumption across sites, the full transition matrix on the \(q^N\) joint states would be the tensor product \(\pi\otimes\pi\otimes\cdots\otimes\pi\) of \(N\) copies of the single-site matrix — independence is exactly what keeps the object \(m\times m\) instead of \(q^N\times q^N\).)

The continuous-time limit: the master equation¶

When generations are frequent and probabilities drift slowly, replace the recursion by a difference quotient over a small time step \(\Delta t\),

\[\frac{p_\alpha(\tau+1) - p_\alpha(\tau)}{\Delta t} = \sum_{\beta\neq\alpha}\left[\frac{\pi_{\alpha\beta}}{\Delta t}p_\beta(\tau) - \frac{\pi_{\beta\alpha}}{\Delta t}p_\alpha(\tau)\right],\]

and let \(\Delta t \to 0\). The left side becomes \(dp_\alpha/dt\); the off-diagonal ratios converge to the transition rates \(R_{\alpha\beta} = \lim_{\Delta t \to 0} \pi_{\alpha\beta}/\Delta t\) (\(\alpha \neq \beta\)), and conservation of probability again fixes the diagonal, \(\sum_\alpha R_{\alpha\beta} = 0\), i.e. \(R_{\beta\beta} = -\sum_{\alpha\neq\beta}R_{\alpha\beta}\) (the total rate out of \(\beta\)). The result is the master equation,

\[\frac{dp_\alpha(t)}{dt} = \sum_{\beta\neq\alpha}\bigl(R_{\alpha\beta}\,p_\beta(t) - R_{\beta\alpha}\,p_\alpha(t)\bigr):\]

the rate of change of the probability of being in state \(\alpha\) is the flux in from every other state minus the flux out to every other state. This single equation is the workhorse for every example in the rest of the chapter — a mutating base, a mutating allele, or a switching enzyme.

Steady states and how fast they are reached¶

Because both \(\pi\) and \(R\) conserve total probability column by column, the all-ones row vector \(\overleftarrow v^{\,*} = (1,\ldots,1)\) is a left eigenvector of both,

\[\overleftarrow v^{\,*}\pi = \overleftarrow v^{\,*} \quad (\text{eigenvalue } 1), \qquad \overleftarrow v^{\,*} R = 0 \quad (\text{eigenvalue } 0).\]

Every eigenvalue of a matrix has a matching right eigenvector, so \(\pi\) and \(R\) each also have a right eigenvector \(\vec p^{\,*}\) with

\[\pi\,\vec p^{\,*} = \vec p^{\,*}, \qquad R\,\vec p^{\,*} = 0.\]

Its entries are the steady-state probabilities: once reached, they no longer change with time. The remaining eigenvalues govern how quickly an arbitrary starting distribution relaxes toward this steady state.

Each site flips between alleles A1 and A2; μ2 is the A1→A2 rate, μ1 the A2→A1 rate. The steady state and the relaxation rate both come out of this pair of numbers alone.

Worked example: the two-state mutation model¶

Take \(m = 2\) states \(A_1, A_2\) with \(A_1 \to A_2\) at rate \(\mu_2\) and \(A_2 \to A_1\) at rate \(\mu_1\) (pictured above). The master equation is the linear system

\[\frac{d}{dt}\begin{pmatrix}p_1\\p_2\end{pmatrix} = \begin{pmatrix}-\mu_2 & \mu_1 \\ \mu_2 & -\mu_1\end{pmatrix}\begin{pmatrix}p_1\\p_2\end{pmatrix}.\]

The rate matrix has eigenvector \(\left(\frac{\mu_1}{\mu_1+\mu_2}, \frac{\mu_2}{\mu_1+\mu_2}\right)\) with eigenvalue \(0\) — already normalised, this is the steady state \(\vec p^{\,*}\) — and eigenvector \((1,-1)\) with eigenvalue \(-(\mu_1+\mu_2)\), which governs relaxation.

Starting from all sites in state \(A_1\), i.e. \(p_1(0)=1\), \(p_2(0)=0\), decompose the initial condition onto the two eigenvectors,

\[\begin{pmatrix}1\\0\end{pmatrix} = \begin{pmatrix}\frac{\mu_1}{\mu_1+\mu_2}\\ \frac{\mu_2}{\mu_1+\mu_2}\end{pmatrix} + \frac{\mu_2}{\mu_1+\mu_2}\begin{pmatrix}1\\-1\end{pmatrix},\]

and apply \(e^{tR}\) term by term — the steady-state eigenvector is unaffected (eigenvalue \(0\)), the other decays as \(e^{-(\mu_1+\mu_2)t}\) (eigenvalue \(-(\mu_1+\mu_2)\)) — giving the closed-form solution

\[p_1(t) = \frac{\mu_1}{\mu_1+\mu_2} + e^{-(\mu_1+\mu_2)t}\frac{\mu_2}{\mu_1+\mu_2}, \qquad p_2(t) = \frac{\mu_2}{\mu_1+\mu_2} - e^{-(\mu_1+\mu_2)t}\frac{\mu_2}{\mu_1+\mu_2}.\]

At long times \(p_1/p_2 \to \mu_1/\mu_2\), the steady-state ratio; note the relaxation rate is the sum \(\mu_1 + \mu_2\), not either rate alone — a fast escape route out of a state speeds up equilibration even if the return route is slow.

Many sites at once: the mutating-population master equation¶

The same two-state model can be reread as a population rather than a single site. Let \(A_1, A_2\) denote two alleles, and let each of \(N\) individuals be replaced each generation by an offspring that usually keeps its parent's allele but mutates to the other one at the same rates \(\mu_1, \mu_2\) as before. The population's state is fully described by the count \(n = N_1\) of \(A_1\)-alleles (\(N_2 = N-n\)), so there are \(N+1\) possible states, with probability \(p(n,t)\) of being in each.

Because a mutation event changes exactly one individual, \(n\) can only move to a neighbouring value in one step: this is a birth-death chain. Reading off which individuals are at risk of mutating in state \(n\): the \(n\) copies of \(A_1\) can each mutate away (total rate \(\mu_2 n\), taking \(n \to n-1\)), and the \(N-n\) copies of \(A_2\) can each mutate in (total rate \(\mu_1(N-n)\), taking \(n \to n+1\)).

Only neighbouring counts communicate — one mutation moves the A1-count n by exactly one, so the population's composition performs a random walk with state-dependent rates.

Collecting the two ways in and the two ways out of an interior state \(0<n<N\) gives the master equation

\[\frac{dp(n,t)}{dt} = \mu_2(n+1)\,p(n+1,t) + \mu_1(N-n+1)\,p(n-1,t) - \mu_2 n\,p(n,t) - \mu_1(N-n)\,p(n,t),\]

with boundary states having only one neighbour,

\[\frac{dp(0,t)}{dt} = \mu_2\,p(1,t) - \mu_1 N\,p(0,t), \qquad \frac{dp(N,t)}{dt} = \mu_1\,p(N-1,t) - \mu_2 N\,p(N,t).\]

The same equation describes a chemical reaction¶

The value of this formalism is that it does not care what the two states are called. Consider an enzyme \(E\) interconverting a substrate between forms \(A\) and \(B\),

\[A + E \underset{b'}{\overset{a'}{\rightleftharpoons}} B + E.\]

In a well-mixed, deterministic treatment that ignores fluctuations, the mean numbers obey

\[\frac{dN_A}{dt} = -a'N_E N_A + b'N_E N_B = -aN_A + b(N-N_A), \qquad a := N_E a',\ \ b := N_E b',\]

with steady state \(N_A^*/N_B^* = b/a\). But when the number of molecules of a species is small — a handful of copies of a particular protein in a single cell, say — the mean is not a good summary of the state, and the full distribution matters. That distribution obeys exactly the birth-death master equation derived above, with \(n = N_A\) in place of the allele count and \((a,b)\) playing the role of \((\mu_2, \mu_1)\) (reaction \(A\to B\) mirrors \(A_1 \to A_2\); \(B \to A\) mirrors \(A_2 \to A_1\)). The steady state is therefore inherited without further work: each of the \(N\) molecules is independently in state \(A\) with probability \(p_A^* = b/(a+b)\) and in state \(B\) with probability \(p_B^* = a/(a+b)\), so the count \(N_A\) is Binomially distributed,

\[p^*(n) = \binom{N}{n}\frac{b^n a^{N-n}}{(a+b)^N}.\]

The general lesson: once a birth-death master equation with a given rate structure has been solved once, it is solved for every process that reduces to that structure, whatever biological label is attached to the two states.

Exercises¶

1. Open reading frames. Assume the four DNA bases A, G, T, C occur independently and with equal probability along a stretch of DNA.

(a) Using the genetic code, find the probability \(p_s\) that a randomly chosen codon (a triplet of bases) is a stop codon.
(b) Find the probability of an open reading frame (ORF) of length \(N\) — a run of \(N\) non-stop codons immediately followed by a stop codon.
(c) E. coli's genome is a closed loop of roughly \(5\times10^6\) bases per strand, readable in six possible frames. Under the random-sequence assumption, how many ORFs of length 600 (about one typical protein) would you expect to see by chance alone?

2 (optional). ORFs in the real genome. Download the complete genome of E. coli K-12.

(a) Scan all three reading frames on the given strand for stop codons and record the distance \(L\) between consecutive stop codons in each frame.
(b) Plot the distribution of \(L\) and compare it with the prediction for a random sequence from Exercise 1.
(c) Estimate a length cutoff \(L_{\mathrm{cut}}\) above which an observed ORF is unlikely to have arisen by chance.

3. Point mutations and base composition. Suppose transversions — mutations between the two purines A, G, or between the two pyrimidines T, C — occur with probability \(q\) per generation, and transitions — any purine-to-pyrimidine change or its reverse — occur with probability \(p\).

(a) Write the resulting \(4\times4\) transition matrix \(\Pi_1\) for the frequencies \((p_A, p_G, p_T, p_C)\), checking that it preserves \(p_A+p_G+p_T+p_C=1\).
(b) Find the eigenvalues of \(\Pi_1\) (its symmetry lets you guess the eigenvectors).
(c) Find \(\Pi_t = \Pi_1^{\,t}\), the transition matrix after \(t\) generations.
(d) Show that, after many generations, all four bases become equally frequent, and estimate — in terms of \(p\) and \(q\) — how many generations that takes.
(e) The result in (d) relies on forward and backward rates between any pair of bases being equal. In the human genome C and G are rarer than A and T, partly because methylated CG pairs are more mutable. Model this by letting A→C and T→G occur at a higher rate \(p_+\) and the reverse (C→A, G→T) at a lower rate \(p_-\), with the other transitions still at \(p\) and transversions still at \(q\). Write the modified transition matrix and find the resulting steady-state base frequencies.

4 (optional). Correlations along the genome. The independent-site models above all predict \(p_{XY} = p_X p_Y\) for the joint probability of finding bases \(X, Y\) at two different positions. Test this against the E. coli genome.

(a) Compute the four single-base frequencies.
(b) Count all 16 pairs of neighbouring bases to get the joint probabilities \(p_{XY}\), and form the correlation matrix \(c^{(1)}_{XY} = p_{XY}/(p_Xp_Y)\).
(c) Repeat for bases separated by \(n = 2, 3, \ldots\) positions, and describe how \(c^{(n)}_{XY}\) decays with the separation \(n\).

5. Activation/deactivation reaction. A finite number \(N\) of copies of a molecule interconvert between an active form \(A\) and an inactive form \(B\) (a phosphate group added or removed), at rates \(a\) and \(b\): \(A \rightleftharpoons B\), with the relevant enzyme concentrations already folded into \(a\) and \(b\).

(a) Write the master equation governing \(p(N_A = n, N_B = N-n, t)\).
(b) Starting from every molecule in state \(A\), i.e. \(p(n, 0) = \delta_{n,N}\), find \(p(n,t)\) at all later times. Guessing the form of the solution and checking it against the master equation from (a) is a reasonable way in.

6 (optional). Human polymorphisms. The Seattle SNP database records, for many polymorphisms, the frequency \(x\) of the minor allele in populations of European (ED) and African (AD) descent.

(a) Histogram the frequencies \(x\) (folding any \(x > 0.5\) to \(1-x\)) and compare against the steady-state frequency distribution \(f(x)\) derived in lecture. What does the comparison suggest about the value of \(N\mu\)?
(b) Compare \(f(x)\) between the two populations and, taking \(\mu = 10^{-8}\), fit \(f(x)\) to estimate an effective population size \(N\) for each.
(c) Compare the frequency spectra of synonymous and non-synonymous mutations against \(f(x)\): which are more common, and why would you expect that?

Sources¶

01-1-1-probability-information.md — MIT OCW 8.592J / HST.452J, Statistical Physics in Biology (Spring 2011), Lecture 1, part 1.1 "Probability & Information": the sequence-as-IID model, the multinomial coefficient and multinomial probability, the Stirling's-approximation derivation of the mixing entropy, and the DNA binary-channel information-content example.
02-1-2-evolving-probabilities.md — same course, Lecture 1, part 1.2 "Evolving Probabilities": the mutation/translation example, the classical-genetics vocabulary of dominant and recessive alleles, the discrete transition matrix, the continuous-time master equation, the steady-state eigenvector argument with its two-state worked example, the mutating-population birth-death master equation, and the enzymatic-reaction correspondence.
psets/01-questions.md — Assignment 1, Questions 1–6 (Questions 2, 4 and 6 are marked optional in the original).
No transcript or written notes were supplied for this session; the chapter follows the slide markdown only. Both slide files carry a front-matter notice that they are a model's reconstruction of a PDF with no text layer, with "every equation... unverified" — the mathematics above has been checked for internal consistency (for instance, the boundary terms of the mutating-population master equation were re-derived from the birth-death rates rather than copied) but not checked against the original PDF pages.
Exercise 6 refers to "the steady-state frequency distribution \(f(x)\) derived in lecture" for allele frequencies under mutation and genetic drift. That derivation (a diffusion-equation result) is not contained in either slide file supplied for this chapter, so it is named here but not covered above.

Contents · 2. Diffusion Limit of Allele Frequencies →