19 Drift is the null hypothesis
(ns meme-constructor.adaptation.post-iii2-drift-null
(:require [scicloj.kindly.v4.kind :as kind]
[clojure.string :as str]
[clojure.set :as set]
[meme-constructor.adaptation.proteins :as p]
[meme-constructor.population :as pop]
[meme-constructor.design :as d]
[meme-constructor.figures :as fig]))Post III.2 of Part III, “A Constructor Theory of Adaptation”.
A site that no selection watches changes at the rate it mutates. How far a site falls below that rate, Kimura’s constraint \(1 - k/\mu\), measures how fine-tuned it is. On real proteins the measure separates the hormone from the connecting peptide and the haem pocket from the surface, and it comes out graded: the C-peptide of insulin, at 0.56, is weakly designed. Design being made shows up not as an excess \(k > \mu\) but as a constraint that switches on at sites that gain a task and off at sites that lose one.
The measure is Kimura’s. The reading through the design conditions, the predictions and all computations are ours.
- 1 − k/μ: how far a site falls below the neutral rate
- graded: the insulin C-peptide, at 0.56, is weakly designed
- not an excess k > μ
- constraint switching on at sites that gain a task, off at sites that lose one
19.1 1. Drift as the null hypothesis
Kimura (1968, 1983) showed that a neutral mutation is substituted at the rate at which it arises. If \(\mu\) is the mutation rate at a site and \(f_0\) the fraction of its mutations that are neutral, the substitution rate is \(k = \mu f_0\). Post 7 proves the toy version (Proposition 7.7).
Definition III.2.1 (Drift as null hypothesis (AD3)). For a site of a gene, the null hypothesis is that no selection acts on it, so that \(k = \mu\). The constraint of the site is
\[1 - \frac{k}{\mu}.\]
It is the share of mutations at the site that selection removes, and so the degree to which the site is fine-tuned, in the sense of the fine-tuning condition (D3). Design being kept shows as constraint. Design being made shows as constraint that switches on at sites that gain a task and off at sites that lose one. The mutation rate \(\mu\) is estimated from the synonymous sites of the same gene, so that \(k/\mu\) is the familiar \(dN/dS\).
The other two design conditions are read off the same data:
Since design is graded (post I.5), a site with constraint 0.5 is weakly designed, not a borderline case.
- positive selection fixes variants faster than drift
- after a single change of the task it does not: selection works on variation already there (section 4)
- the constraint switches on at sites that gain a task
- and off at sites that lose one
- needs a task that keeps moving (post III.3, Proposition III.3.12)
19.2 2. The real-protein experiment (Y1)
Conjecture III.2.2 (Constraint separates design within a gene). Within one gene, the constraint is high where the protein does its task and low where it does not. Predictions:
- Insulin. The A and B chains, which form the hormone, are conserved. The C-peptide, which is cut out and discarded, changes at nearly the rate of the fibrinopeptides, the fastest proteins Kimura knew.
- Haemoglobin β. Residues in contact with the haem are more conserved than surface residues.
- Histone H4 is nearly invariant.
Refutation: the C-peptide is as conserved as the A and B chains, or the haem contacts are not more conserved than the surface.
19.2.1 Data
A library function downloads everything once into results/postIII2/ and aligns it; this post reads the cache. On a fresh checkout the first render downloads (about 1.5 minutes).
(def results (p/ensure-results! "results/postIII2"))(defn r3 [x] (if (and (number? x) (not (Double/isNaN (double x)))) (/ (Math/round (* 1000.0 x)) 1000.0) x))(kind/table
{:row-maps
(vec (for [[g label] [[:ins "preproinsulin"] [:hbb "haemoglobin β"]
[:h4 "histone H4"] [:fga "fibrinopeptide A"]]
:let [{:keys [alignment]} (get-in results [:genes g])]]
{:protein label
:sequences (count alignment)
:query (get-in results [:queries g :query])}))})| protein | sequences | query |
|---|---|---|
| preproinsulin | 18 | gene_exact:INS AND taxonomy_id:40674 AND reviewed:true AND length:[100 TO 120] |
| haemoglobin β | 122 | gene_exact:HBB AND taxonomy_id:40674 AND reviewed:true AND length:[145 TO 147] |
| histone H4 | 118 | protein_name:"histone H4" AND taxonomy_id:40674 AND length:[103 TO 103] |
| fibrinopeptide A | 40 | gene_exact:FGA AND taxonomy_id:40674 AND reviewed:true |
The queries went to https://rest.uniprot.org/uniprotkb/stream (UniProt Consortium 2025) on the date below; taxonomy_id:40674 is Mammalia. Rules fixed before any rate was computed:
- One entry per genus (smallest accession, human always kept), to thin out primates and bovids. This leaves 122 genera for haemoglobin β.
- Insulin: reviewed entries of length 100 to 120, the full precursor; most reviewed insulins are only the 51 residues of the mature chains.
INS_RODSP(no species) is dropped. - Histone H4: only two mammalian entries are reviewed, so we take all entries of length 103 and per species the most frequent sequence.
- Fibrinopeptide A: the annotated peptide, or the whole entry where the entry is the peptide.
- Coding DNA: the first ENA coding sequence of an entry whose translation equals the UniProt sequence.
- Structure: PDB 2HHB, human deoxyhaemoglobin (Fermi et al. 1984).
(select-keys results [:downloaded :mafft :pdb-url :ena-url]){:downloaded "2026-09-29",
:mafft
{:version "v7.526 (2024/Apr/26)",
:command
"tools/mafft-linux64/mafft.bat --auto --quiet in.fa > out.fa"},
:pdb-url "https://files.rcsb.org/download/2HHB.pdb",
:ena-url "https://www.ebi.ac.uk/ena/browser/api/fasta/<protein id>"}Alignment is MAFFT (Katoh & Standley 2013), version and call as above, with default parameters under --auto.
The measure on protein data. Without DNA we cannot estimate \(\mu\), so, as Kimura did, we state the constraint relative to the fastest region of the same gene: the C-peptide for insulin, the surface for haemoglobin β. Per-site diversity, the probability that two species differ at the site (gaps ignored), is used for the plots. The region rate, the fraction \(p\) of differing sites for each species pair, corrected as \(-\ln(1-p)\) and averaged over pairs, is used for the tables, with a bootstrap interval over sites (200 resamples, seed 1). This assumes \(\mu\) is the same along the gene; section 3 checks it.
19.2.2 Insulin
Region boundaries come from the UniProt features of human preproinsulin (P01308), not from us:
(def ins (get-in results [:genes :ins]))(p/feature-regions (:ref-features ins) ["SIGNAL" "PEPTIDE" "PROPEP"]){"C peptide" [[57 87]],
"Insulin A chain" [[90 110]],
"Insulin B chain" [[25 54]],
"SIGNAL" [[1 24]]}We add the two pairs of basic residues where the prohormone is cut (positions 55–56 and 88–89).
(def ins-ranges
{"signal peptide" [[1 24]] "B chain" [[25 54]] "C-peptide" [[57 87]]
"A chain" [[90 110]] "cleavage sites" [[55 56] [88 89]]})(def ins-region-of
(into {} (for [[nm rs] ins-ranges [a b] rs q (range a (inc b))] [q nm])))(def ins-profile (p/site-profile (:alignment ins) "P01308"))(def region-colour
{"signal peptide" "#9467bd" "B chain" "#d62728" "C-peptide" "#7f7f7f"
"A chain" "#1f77b4" "cleavage sites" "#2ca02c"})(kind/plotly
{:data (vec (for [[nm colour] region-colour
:let [xs (filter #(= nm (ins-region-of (:pos %))) ins-profile)]]
{:type "bar" :name nm :x (mapv :pos xs) :y (mapv :diversity xs)
:text (mapv :ref-aa xs) :marker {:color colour}}))
:layout {:title (str "Preproinsulin, per-site diversity (" (count (:alignment ins)) " mammals)")
:xaxis {:title "position in human preproinsulin"}
:yaxis {:title "P(two species differ)" :range [0 1]}
:height 400}})(def ins-summary
(p/region-summary (:alignment ins)
(update-vals ins-ranges #(p/region-columns (:alignment ins) "P01308" %))
"C-peptide" {:reps 200 :seed 1}))(defn table-rows [rows]
(mapv #(-> (update-vals % r3) (update :interval (fn [iv] (mapv r3 iv)))) rows))(kind/table {:row-maps (table-rows ins-summary)})| region | sites | mean-diversity | invariant | rate | relative-rate | interval | relative-constraint |
|---|---|---|---|---|---|---|---|
| signal peptide | 24.0 | 0.252 | 0.25 | 0.303 | 0.767 | |
0.233 |
| B chain | 30.0 | 0.101 | 0.5 | 0.122 | 0.309 | |
0.691 |
| C-peptide | 31.0 | 0.307 | 0.226 | 0.395 | 1.0 | |
0.0 |
| A chain | 21.0 | 0.133 | 0.524 | 0.156 | 0.394 | |
0.606 |
| cleavage sites | 4.0 | 0.0 | 1.0 | 0.0 | 0.0 | |
1.0 |
rate is the mean corrected distance over the 153 species pairs, relative-rate the rate divided by that of the C-peptide with its 95 % bootstrap interval, relative-constraint is one minus it.
(defn diversities [profile pred] (keep #(when (pred %) (:diversity %)) profile))(def ins-p
(p/permutation-p (diversities ins-profile #(#{"A chain" "B chain"} (ins-region-of (:pos %))))
(diversities ins-profile #(= "C-peptide" (ins-region-of (:pos %))))
5000 1))ins-p1.999600079984003E-4Observation III.2.3 (Insulin). 18 mammals, bootstrap and permutation seed 1. The B chain changes at 0.31 and the A chain at 0.39 of the C-peptide rate; both intervals lie well below 1. Half of their sites are invariant, against a fifth of the C-peptide. The per-site diversities of the A and B chains are lower than those of the C-peptide with one-sided permutation p = 1/5001, the smallest value 5000 permutations can give (sites are not independent, so this is a rough guide). The four cleavage residues never change. The first refutation condition did not occur.
The signal peptide was not in the prediction. It changes at 0.77 of the C-peptide rate, and its interval includes 1, although it has a task: it sends the precursor into the secretory pathway. Section 5 comes back to this.
19.2.3 Haemoglobin β
From 2HHB we compute, for the β chains B and D, the residues with any atom within 5 Å of a haem atom, and each residue’s relative solvent accessibility in the tetramer (Shrake & Rupley 1973, probe 1.4 Å, 200 points per atom, maxima of Tien et al. 2013). PDB residue \(r\) is UniProt position \(r+1\) of P68871, since the initiator methionine is removed. Haem contact: touches a haem in B or D. Surface: otherwise, with relative accessibility at least 0.25. Buried: the rest.
(def structure (:structure results))(def hbb-class (p/hbb-site-classes structure 0.25))(frequencies (vals hbb-class)){:surface 54, :buried 71, :haem 21}(def hbb (get-in results [:genes :hbb]))(def hbb-profile (p/site-profile (:alignment hbb) "P68871"))(def haem-positions (sort (keep (fn [[q c]] (when (= :haem c) q)) hbb-class)))(str/join " " (map #(str (:ref-aa (nth hbb-profile (dec %))) %) haem-positions))"L32 T39 F42 F43 F46 H64 K67 V68 A71 F72 F86 L89 L92 H93 L97 V99 N103 F104 L107 V138 L142"These are UniProt positions. In the usual mature numbering subtract 1: they include the distal histidine E7 (His64 here, His63 mature) and the proximal histidine F8 (His93 here, His92 mature).
(def class-colour {:haem "#d62728" :buried "#1f77b4" :surface "#7f7f7f"})(kind/plotly
{:data (vec (for [[c colour] class-colour
:let [xs (filter #(= c (hbb-class (:pos %))) hbb-profile)]]
{:type "bar" :name (name c) :x (mapv :pos xs) :y (mapv :diversity xs)
:text (mapv :ref-aa xs) :marker {:color colour}}))
:layout {:title (str "Haemoglobin β, per-site diversity (" (count (:alignment hbb)) " mammalian genera)")
:xaxis {:title "position in human HBB (UniProt, Met = 1)"}
:yaxis {:title "P(two species differ)" :range [0 1]}
:height 400}})(def hbb-ranges
(into {} (for [[c qs] (group-by hbb-class (keys hbb-class))] [c (mapv (fn [q] [q q]) qs)])))(def hbb-summary
(p/region-summary (:alignment hbb)
(update-vals hbb-ranges #(p/region-columns (:alignment hbb) "P68871" %))
:surface {:reps 200 :seed 1}))(kind/table {:row-maps (table-rows hbb-summary)})| region | sites | mean-diversity | invariant | rate | relative-rate | interval | relative-constraint |
|---|---|---|---|---|---|---|---|
| surface | 54.0 | 0.289 | 0.093 | 0.35 | 1.0 | |
0.0 |
| buried | 71.0 | 0.13 | 0.31 | 0.145 | 0.414 | |
0.586 |
| haem | 21.0 | 0.041 | 0.667 | 0.043 | 0.123 | |
0.877 |
(def hbb-p
(p/permutation-p (diversities hbb-profile #(= :haem (hbb-class (:pos %))))
(diversities hbb-profile #(= :surface (hbb-class (:pos %))))
5000 1))hbb-p1.999600079984003E-4Observation III.2.4 (Haemoglobin β). 122 genera, seed 1. The 21 haem contacts change at 0.12 of the surface rate (interval 0.02 to 0.31), buried non-contact residues at 0.41. Two thirds of the haem contacts are invariant, against one in eleven surface sites. Permutation p for “haem contacts less diverse than surface” is again 1/5001. The second refutation condition did not occur. The order haem < buried < surface shows a second, weaker task besides holding the haem: packing the fold.
19.2.4 Histone H4 and fibrinopeptide A
(def h4 (get-in results [:genes :h4 :alignment]))(def h4-human (:seq (first (filter #(= "P62805" (:id %)) h4)))){:genera (count h4)
:identical-to-human (count (filter #(= h4-human (:seq %)) h4))}{:genera 118, :identical-to-human 109}Observation III.2.5 (Histone H4). 109 of 118 genera carry exactly the human sequence. The nine others differ from it at one to five sites and all come from unreviewed, mostly genome-predicted entries, so some of these differences may be sequencing or prediction errors rather than real variants. Either way H4 is nearly invariant, as predicted.
For the comparison between genes the species must be the same, since the rate depends on how far apart the species are. The genera present in all four data sets form a small common panel:
(defn genus [e] (first (str/split (:species e) #"\s+")))(def panel
(apply set/intersection
(for [g [:ins :hbb :h4 :fga]] (set (map genus (get-in results [:genes g :alignment]))))))(sort panel)("Bos" "Canis" "Cavia" "Chlorocebus" "Homo" "Macaca" "Ovis" "Sus")(defn panel-rate [g ref ranges]
(let [aln (get-in results [:genes g :alignment])]
(p/region-rate (p/restrict aln panel)
(if ranges
(p/region-columns aln ref ranges)
(range (count (:seq (first aln))))))))(def panel-rows
(let [fpa (panel-rate :fga nil nil)]
(vec (for [[label g ref ranges] (concat
[["fibrinopeptide A" :fga nil nil]]
(for [nm ["C-peptide" "signal peptide" "A chain" "B chain"]]
[(str "insulin " nm) :ins "P01308" (ins-ranges nm)])
(for [c [:surface :buried :haem]]
[(str "haemoglobin β " (name c)) :hbb "P68871" (hbb-ranges c)])
[["histone H4" :h4 nil nil]])
:let [rate (panel-rate g ref ranges)]]
{:region label :rate (r3 rate) :relative-to-fibrinopeptide (r3 (/ rate fpa))}))))(kind/table {:row-maps panel-rows})| region | rate | relative-to-fibrinopeptide |
|---|---|---|
| fibrinopeptide A | 0.446 | 1.0 |
| insulin C-peptide | 0.332 | 0.745 |
| insulin signal peptide | 0.241 | 0.541 |
| insulin A chain | 0.17 | 0.381 |
| insulin B chain | 0.116 | 0.26 |
| haemoglobin β surface | 0.335 | 0.752 |
| haemoglobin β buried | 0.134 | 0.3 |
| haemoglobin β haem | 0.026 | 0.059 |
| histone H4 | 0.005 | 0.011 |
Observation III.2.6 (Common panel). Eight genera. On the same species the C-peptide changes at 0.75 of the fibrinopeptide rate, about as fast as the haemoglobin surface, and histone H4 at about 1 %. The ordering fibrinopeptide > C-peptide > A, B chains, and surface > buried > haem > H4, is the one Kimura’s constraint predicts. With eight genera the numbers are rough.
19.3 3. Synonymous sites: from relative to absolute (Y2)
Is the C-peptide itself free of constraint, and is \(\mu\) the same along the gene? Both need coding DNA. We threaded the ENA coding sequences onto a MAFFT alignment of their translations and counted, per species pair, synonymous and nonsynonymous differences by Nei & Gojobori (1986) with the Jukes–Cantor correction, dropping pairs with gene-wide \(p_S > 0.5\). This is the conservation half of the dN/dS experiment (Y2), whose other half is about MHC. Per region we give \(dN/dS\) with the region’s own \(dS\), and \(dN/dS_\text{gene}\) with \(dS\) of the whole gene, which assumes a shared \(\mu\) but is far less noisy: a region of 20 codons has about five synonymous sites.
(def ins-dnds (p/dnds-summary (get-in results [:cds :ins]) (assoc ins-ranges "whole gene" [[1 110]]) {:max-ps 0.5}))(def hbb-dnds (p/dnds-summary (get-in results [:cds :hbb]) (assoc hbb-ranges "whole gene" [[1 147]]) {:max-ps 0.5}))(defn dnds-rows [label {:keys [regions]}]
(vec (for [[nm {:keys [n-pairs dn ds ds-gene omega omega-gene]}] regions]
{:gene label :region (name nm) :pairs n-pairs
:dN (r3 dn) :dS (r3 ds) :dS-gene (r3 ds-gene)
:dN-over-dS (r3 omega) :dN-over-dS-gene (r3 omega-gene)
:constraint (r3 (- 1.0 omega-gene))}))){:insulin-cds (count (get-in results [:cds :ins :alignment]))
:haemoglobin-cds (count (get-in results [:cds :hbb :alignment]))}{:insulin-cds 18, :haemoglobin-cds 26}(kind/table {:row-maps (into (dnds-rows "insulin" ins-dnds) (dnds-rows "haemoglobin β" hbb-dnds))})| constraint | dN-over-dS | dS-gene | region | pairs | dS | gene | dN | dN-over-dS-gene |
|---|---|---|---|---|---|---|---|---|
| 0.646 | 0.489 | 0.475 | signal peptide | 144 | 0.344 | insulin | 0.168 | 0.354 |
| 0.863 | 0.095 | 0.475 | B chain | 144 | 0.685 | insulin | 0.065 | 0.137 |
| 0.557 | 0.416 | 0.475 | C-peptide | 144 | 0.507 | insulin | 0.211 | 0.443 |
| 0.88 | 0.11 | 0.475 | A chain | 144 | 0.519 | insulin | 0.057 | 0.12 |
| 0.996 | 0.003 | 0.46 | cleavage sites | 132 | 0.543 | insulin | 0.002 | 0.004 |
| 0.758 | 0.242 | 0.475 | whole gene | 144 | 0.475 | insulin | 0.115 | 0.242 |
| 0.569 | 0.28 | 0.401 | surface | 262 | 0.617 | haemoglobin β | 0.173 | 0.431 |
| 0.774 | 0.252 | 0.401 | buried | 262 | 0.36 | haemoglobin β | 0.091 | 0.226 |
| 0.943 | 0.085 | 0.401 | haem | 262 | 0.269 | haemoglobin β | 0.023 | 0.057 |
| 0.726 | 0.274 | 0.401 | whole gene | 262 | 0.401 | haemoglobin β | 0.11 | 0.274 |
Observation III.2.7 (dN/dS). Insulin: 18 coding sequences, 144 of 153 pairs kept; haemoglobin β: 26, 262 pairs. With \(\mu\) from the whole gene, the constraint \(1 - dN/dS_\text{gene}\) is 0.86 (B chain), 0.88 (A chain), 0.56 (C-peptide), 0.65 (signal peptide); for haemoglobin β 0.94 (haem contacts), 0.77 (buried), 0.57 (surface). Regional \(dS\) scatter from 0.27 to 0.69 around gene-wide 0.48 and 0.40. No region has \(dN/dS > 1\).
So the C-peptide is not near neutral: more than half of its nonsynonymous mutations are removed by selection. The protein-level comparison hid this, because the fibrinopeptide is not neutral either. Our prediction had it near neutral; it is weakly designed, which fits a role in folding the precursor. That contradicts nothing, since design is graded. Within the gene the measure still separates hormone from connecting peptide clearly.
The other half of the dN/dS experiment (Y2), at the antigen-binding sites of MHC class I, is not done. There Hughes & Nei (1988) found nonsynonymous exceeding synonymous substitution. It needs allele data, a different download.
19.4 4. Positive selection in the toy (Y5)
The positive-selection experiment (Y5) tested the obvious idea.
Conjecture III.2.8 (Making design shows as an excess). After the task T changes, the sites affected by the change show \(k > \mu\) for a short time and then \(k < \mu\). Refutation: \(k\) stays below \(\mu\) also at the change.
We reuse the world, copy noise and source variant of post 7. For 100 generations the population lives in the old world; then T changes and it runs on. Rates are counted along one line of descent per seed, traced back from host 0 of the last generation, in windows of equal length and pooled over the seeds.
(def landscape (pop/make-landscape {:block-size 5 :seed 1 :size-cost 0.0}))(def copy-noise {:p-loss 0.02 :p-gain 0.01 :p-jump 0.05})(def source
;; the source variant of Observation 7.8: all points found, clouds in their plateau, half the levers
(let [t (pop/tradition landscape {:n-points 4 :n-clouds 2 :psy-bits 5})]
(reduce (fn [v c] (bit-set v (first (remove #(bit-test v %) (pop/variant->tasks (:mask c))))))
t (pop/blocks-of-type landscape :cloud))))(defn lineage-across
"Line of descent over two runs: `before` generations in world w1 with
params p1, then `after` generations in w2 with params p2."
[w1 p1 w2 p2 before after seed]
(let [h1 (pop/run (merge p1 {:landscape w1 :init source :generations before :n-agents 60 :seed seed}))
h2 (pop/run (merge p2 {:landscape w2 :init (mapv :variant (peek (:generations h1)))
:generations after :n-agents 60 :seed (+ 1000 seed)}))
l2 (d/lineage h2 0)]
(into (d/lineage h1 (:parent (first l2))) l2)))(defn window-rates
"k/mu per site class and window of generations, pooled over lineages."
[lines classes width]
(let [n (count (first lines))]
(vec (for [a (range 0 n width)
:let [per (mapv #(d/site-counts copy-noise (subvec % a (min n (+ a width))) (:n-tasks landscape)) lines)]]
(into {:from (inc a)}
(for [[nm ss] classes]
(let [k (reduce + (for [pc per s ss] (:k (pc s))))
mu (reduce + (for [pc per s ss] (:mu (pc s))))]
[nm (/ k mu)])))))))(defn sites-of [b] (range (* 5 (:index b)) (* 5 (inc (:index b)))))(defn rate-plot [title rows classes colours change]
(kind/plotly
{:data (vec (for [[nm _] classes]
{:type "scatter" :mode "lines+markers" :name nm :line {:color (colours nm)}
:x (mapv :from rows) :y (mapv #(get % nm) rows)}))
:layout {:title title :height 400
:xaxis {:title "generation (start of window)"}
:yaxis {:title "k / μ"}
:shapes [{:type "line" :x0 change :x1 change :y0 0 :y1 1 :yref "paper"
:line {:dash "dot" :color "#444"}}
{:type "line" :x0 0 :x1 1 :xref "paper" :y0 1 :y1 1
:line {:dash "dash" :color "#aaa"}}]}}))19.4.1 First test: T changes where the small points are
The new world redraws the hidden pattern of every small point (pop/counterfactual-reality). We take the first redraw seed for which all four points change, seed 1. The criterion is \(c_\text{rat}\) before and after. Affected sites are the point sites whose bit differs between old and new pattern.
(def world-b (pop/counterfactual-reality landscape 1))(def point-pairs (filter #(= :point (:type (first %))) (map vector (:blocks landscape) (:blocks world-b))))(def must-change
(set (for [[a b] point-pairs s (sites-of a)
:when (not= (bit-test (:pattern a) s) (bit-test (:pattern b) s))] s)))(def points-classes
{"affected point sites" must-change
"other point sites" (set/difference (set (mapcat (comp sites-of first) point-pairs)) must-change)
"neutral sites" (set (mapcat sites-of (pop/blocks-of-type landscape :neutral)))})(def rat-lines
(mapv #(lineage-across landscape {:copy copy-noise :criterion {:kind :rat}}
world-b {:copy copy-noise :criterion {:kind :rat}} 100 500 %)
(range 1 11)))(def rat-rows (window-rates rat-lines points-classes 25))(def y5-colours {"affected point sites" "#d62728" "other point sites" "#ff7f0e" "neutral sites" "#7f7f7f"
"new lever sites" "#d62728" "old lever sites" "#1f77b4" "all other sites" "#7f7f7f"})(rate-plot "c_rat, new small points at generation 100 (seeds 1–10, windows of 25)"
rat-rows points-classes y5-colours 100)(mapv #(update-vals % r3) (take 8 rat-rows))[{:from 1.0,
"affected point sites" 0.0,
"other point sites" 0.0,
"neutral sites" 0.843}
{:from 26.0,
"affected point sites" 0.0,
"other point sites" 0.0,
"neutral sites" 1.128}
{:from 51.0,
"affected point sites" 0.0,
"other point sites" 0.0,
"neutral sites" 0.887}
{:from 76.0,
"affected point sites" 0.0,
"other point sites" 0.0,
"neutral sites" 0.972}
{:from 101.0,
"affected point sites" 0.973,
"other point sites" 0.792,
"neutral sites" 0.775}
{:from 126.0,
"affected point sites" 0.871,
"other point sites" 0.792,
"neutral sites" 1.072}
{:from 151.0,
"affected point sites" 0.59,
"other point sites" 0.783,
"neutral sites" 0.833}
{:from 176.0,
"affected point sites" 0.714,
"other point sites" 0.719,
"neutral sites" 0.736}]Observation III.2.9 (First test failed). Seeds 1 to 10. Before the change the point sites have \(k = 0\). In the first window after it the affected sites rise to \(k/\mu = 0.97\), the other sites of the changed points to 0.79, and then both decline slowly towards 0 over some 400 generations. The affected sites never exceed \(\mu\) in any window. The refutation condition occurred.
The reason lies in the model. Under \(c_\text{rat}\) selection acts by discarding: a host rejects a worse copy and copies the same parent again, and parents are drawn uniformly. A host that found the new point has no more offspring than any other, so the new pattern spreads by drift and is only protected from loss. An excess \(k > \mu\) needs a fixation probability above the neutral \(1/N\), and that needs carriers that leave more offspring. We conjecture that selection by mere discarding never raises the fixation probability, whatever T does. The argument is the one just given; a proof is pending.
19.4.2 Second test: differential reproduction
So we use that channel: criterion \(c_\text{none}\) (every copy kept), parents drawn with weight \(e^{6 F_\text{psy}}\). The task is the lever count. At generation 100 the two neutral blocks become the lever blocks and the old lever blocks become neutral. Seeds 1 to 20, windows of 10.
(def world-c (assoc landscape :psy-mask (reduce bit-or (map :mask (pop/blocks-of-type landscape :neutral)))))(defn mask-sites [m] (set (filter #(bit-test m %) (range (:n-tasks landscape)))))(def lever-classes
(let [new-l (mask-sites (:psy-mask world-c)) old-l (mask-sites (:psy-mask landscape))]
{"new lever sites" new-l
"old lever sites" old-l
"all other sites" (set/difference (set (range (:n-tasks landscape))) new-l old-l)}))(def repro {:copy copy-noise :criterion {:kind :none} :psy-bias 6.0})(def repro-lines (mapv #(lineage-across landscape repro world-c repro 100 200 %) (range 1 21)))(def repro-rows (window-rates repro-lines lever-classes 10))(rate-plot "c_none, parents weighted by exp(6 F_psy), lever set moved at generation 100 (seeds 1–20)"
repro-rows lever-classes y5-colours 100)(mapv #(update-vals % r3) (take 16 (drop 6 repro-rows)))[{:from 61.0,
"new lever sites" 1.151,
"old lever sites" 0.034,
"all other sites" 1.05}
{:from 71.0,
"new lever sites" 0.835,
"old lever sites" 0.057,
"all other sites" 1.012}
{:from 81.0,
"new lever sites" 1.081,
"old lever sites" 0.045,
"all other sites" 1.062}
{:from 91.0,
"new lever sites" 1.382,
"old lever sites" 0.296,
"all other sites" 1.12}
{:from 101.0,
"new lever sites" 0.646,
"old lever sites" 1.029,
"all other sites" 0.952}
{:from 111.0,
"new lever sites" 0.16,
"old lever sites" 1.086,
"all other sites" 1.09}
{:from 121.0,
"new lever sites" 0.0,
"old lever sites" 1.025,
"all other sites" 1.087}
{:from 131.0,
"new lever sites" 0.0,
"old lever sites" 0.98,
"all other sites" 0.912}
{:from 141.0,
"new lever sites" 0.0,
"old lever sites" 1.004,
"all other sites" 0.953}
{:from 151.0,
"new lever sites" 0.023,
"old lever sites" 1.013,
"all other sites" 1.133}
{:from 161.0,
"new lever sites" 0.057,
"old lever sites" 0.986,
"all other sites" 1.094}
{:from 171.0,
"new lever sites" 0.057,
"old lever sites" 0.977,
"all other sites" 1.063}
{:from 181.0,
"new lever sites" 0.034,
"old lever sites" 1.093,
"all other sites" 0.994}
{:from 191.0,
"new lever sites" 0.011,
"old lever sites" 0.891,
"all other sites" 0.702}
{:from 201.0,
"new lever sites" 0.0,
"old lever sites" 1.238,
"all other sites" 0.969}
{:from 211.0,
"new lever sites" 0.0,
"old lever sites" 0.757,
"all other sites" 1.205}]Observation III.2.10 (Second test). Seeds 1 to 20. Before the change the old lever sites are held (\(k/\mu\) 0.03 to 0.06) and the future ones drift (0.84 to 1.15). After it the new lever sites fall to 0.65, then 0.16, then about 0; the old ones rise to about 1, released like a pseudogene (Li, Gojobori & Nei 1981); the other sites stay near 1. No window after the change has \(k > \mu\) at the new lever sites: the refutation condition occurred here too. The one high window is just before the change (1.38 in generations 91 to 100), but the other sites reach 1.12 there and range from 0.70 to 1.13, so it is weak evidence.
Again the model explains it. When T changes, the future lever sites have drifted for 100 generations and about half their bits are set in any host. Selection does not wait for new mutations; it picks the hosts that already carry most of the new levers. The substitutions on the winning line of descent happened before T changed, while they were neutral. Selection on standing variation leaves no excess when it acts.
What a change of T does produce is the signature of Definition III.2.1: at the new lever sites \(k/\mu\) falls from about 1 to about 0 within twenty generations, so the constraint switches on, and at the old ones it switches off.
19.5 5. Limits of the measure
- a duplicated and repurposed gene is one variation step
- but shows only as many conserved sites
- Post I.4 measures variation algorithmically (C6) and shows how the count misstates it
- the insulin signal peptide has a task but low constraint
- it needs hydrophobic residues, and many will do
- its design lies in a property of many sites
- codon bias, splicing signals, mRNA structure
- then dS underestimates μ and 1 − dN/dS underestimates constraint
- diversity counts differing species, not substitutions on a tree
- primates and bovids overrepresented; sites not independent: p-values rough
19.6 Summary
- Drift is the null hypothesis (AD3): without selection \(k = \mu\). The constraint \(1 - k/\mu\) is the degree of fine-tuning of a site.
- On real proteins the measure separates design within genes. Insulin’s A and B chains change at a third of the C-peptide rate, haem contacts of haemoglobin β at an eighth of the surface rate, histone H4 is invariant in 109 of 118 genera. Neither refutation condition of the real-protein experiment (Y1) occurred.
- With synonymous sites the constraint is about 0.87 for the insulin chains, 0.94 for the haem contacts and 0.56 for the C-peptide, which is weakly designed, not neutral as we predicted.
- In the toy a single change of T gave no excess \(k > \mu\), neither under selection by discarding nor with differential reproduction (Y5). It switched the constraint on at the new sites and off at the old ones. That switch is the signature of design being made.
19.7 Open problems
- Selection by discarding. Prove that selection which only rejects worse copies, without giving carriers more offspring, cannot raise a fixation probability. This matters for the rational criterion \(c_\text{rat}\) of Part II.
- MHC. The positive-selection half of the dN/dS experiment (Y2), at MHC antigen-binding sites, is not done. Without it the excess \(k > \mu\) has no test on real data in this book.
- A measure of design in a property, for example conservation of hydrophobicity rather than residue identity, for the signal peptide.
- A tree. All rates are pairwise averages. Per-site rates on a phylogeny would remove most of the sampling bias.
19.8 References
- Fermi, G., Perutz, M. F., Shaanan, B. & Fourme, R. (1984). The crystal structure of human deoxyhaemoglobin at 1.74 Å resolution. Journal of Molecular Biology 175, 159–174. doi:10.1016/0022-2836(84)90472-8 (verified; the source of PDB 2HHB)
- Hughes, A. L. & Nei, M. (1988). Pattern of nucleotide substitution at major histocompatibility complex class I loci reveals overdominant selection. Nature 335, 167–170. doi:10.1038/335167a0 (verified)
- Katoh, K. & Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Molecular Biology and Evolution 30, 772–780. doi:10.1093/molbev/mst010 (verified)
- Kimura, M. (1968). Evolutionary rate at the molecular level. Nature 217, 624–626. (verified in post 7)
- Kimura, M. (1983). The Neutral Theory of Molecular Evolution. Cambridge University Press. (verified in post 7)
- Li, W.-H., Gojobori, T. & Nei, M. (1981). Pseudogenes as a paradigm of neutral evolution. Nature 292, 237–239. (verified in post 7)
- Nei, M. & Gojobori, T. (1986). Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Molecular Biology and Evolution 3, 418–426. doi:10.1093/oxfordjournals.molbev.a040410 (verified)
- Shrake, A. & Rupley, J. A. (1973). Environment and exposure to solvent of protein atoms. Lysozyme and insulin. Journal of Molecular Biology 79, 351–371. doi:10.1016/0022-2836(73)90011-9 (verified)
- Tien, M. Z., Meyer, A. G., Sydykova, D. K., Spielman, S. J. & Wilke, C. O. (2013). Maximum allowed solvent accessibilites of residues in proteins. PLoS ONE 8, e80635. doi:10.1371/journal.pone.0080635 (verified, including the theoretical values used in the code)
- UniProt Consortium (2025). UniProt: the Universal Protein Knowledgebase in 2025. Nucleic Acids Research 53, D609–D617. doi:10.1093/nar/gkae1010 (verified)
source: book2/notebooks/meme_constructor/adaptation/post_iii2_drift_null.clj