Example: compensatory evolution at a salt bridge

We construct a pair landscape in which only two amino-acid configurations have high fitness: lysine–glutamate and glutamate–lysine, the classic charge-swap of an electrostatic interaction. The single-mismatch configurations (K,K and E,E) are deleterious but passable, and every other configuration is strongly deleterious. With this landscape the pair toggles between K,E and E,K through short-lived intermediates: a substitution at one site followed almost immediately by the compensating substitution at the other.

Setup

# Load phylogeny
mammals <- Phylogeny("../inst/extdata/mammals.newick")

# Flat site fitness: all structure comes from the pair coupling
f <- rep(1, 20)
Ne <- 1000

# Amino acid order is:
# A, C, D, E, F, G, H, I, K, L, M, N, P, Q, R, S, T, V, W, Y
K <- 9; E <- 4

# Population-scaled selection coefficients (2*Ne*s) of each pair
# configuration, relative to the compatible charge-swap pairs:
S <- matrix(-20, 20, 20)   # everything else:  2Ns = -20
S[K, K] <- -8              # mismatch K,K:     2Ns =  -8
S[E, E] <- -8              # mismatch E,E:     2Ns =  -8
S[K, E] <- 0               # compatible:       2Ns =   0
S[E, K] <- 0

# The scaled fitness of pair (i, j) is 2*Ne*(f[i] + f[j])*delta[i, j];
# with flat f this is 4*Ne*delta, so a coefficient of S corresponds to:
delta <- 1 + S / (4 * Ne)

The intermediate depth sets how tightly the two substitutions couple: the rate of escaping an intermediate grows roughly like \(e^{|2Ns|}\), while non-synonymous scaling (below) pins the overall rate of excursions. A deeper intermediate therefore gives the same number of toggles, each resolved much faster.

Model

Non-synonymous scaling makes one unit of branch length one expected amino-acid-changing substitution per codon, so the compensatory events are delivered at a rate the tree’s branch lengths control directly.

use_genetic_code("Standard nuclear")

hky <- HasegawaKishinoYano(equilibrium = rep(0.25, 4))

co <- CoEvolution(
    population_size = Ne,
    mutation_rate = 1e-6,
    nucleotide_model = hky,
    fitness_1 = f,
    fitness_2 = f,
    delta = delta,
    scaling_type = "non-synonymous")

At equilibrium the two compatible configurations hold essentially all the mass:

eq <- as.numeric(co$equilibrium)
sum(sort(eq, decreasing = TRUE)[1:8])   # the 8 codon pairs spelling K,E / E,K

## [1] 0.9996637

Simulation

s <- sample_sequence(model = co, length = 30)

sim <- simulate_over_phylogeny(
    phylogeny = mammals,
    model = co,
    sequence = s)

The interaction keeps the pair under strong purifying selection, so almost every event is synonymous churn within the K and E codon families; the compensatory dynamics live in the small non-synonymous remainder:

nrow(sim$substitutions)

## [1] 71618

ns <- subset(sim$substitutions, !(first_synonymous & second_synonymous))
nrow(ns)

## [1] 409

Substitutions into the two intermediates (excursion starts) are balanced by non-synonymous substitutions back into the compatible configurations (excursion ends), and no other amino-acid pair is ever reached:

sort(table(ns$pair_amino_acid_to), decreasing = TRUE)

## 
## E,K K,K E,E K,E 
## 115 105 101  88

Each excursion starts with a substitution into K,K or E,E and resolves with the next non-synonymous substitution at the same pair – either completing the toggle (landing on the opposite compatible configuration) or reverting. A small helper classifies each excursion and records how long the intermediate lasted:

count_toggles <- function(substitutions, compatible = c("K,E", "E,K")) {
  toggles <- 0; reversions <- 0; waits <- c()
  for (s0 in unique(substitutions$site)) {
    for (nd in unique(substitutions$node)) {
      g <- substitutions[substitutions$site == s0 & substitutions$node == nd, ]
      if (nrow(g) < 2) next
      g <- g[order(g$time), ]
      for (i in seq_len(nrow(g) - 1)) {
        entered  <- !(g$pair_amino_acid_to[i] %in% compatible)
        resolved <- g$pair_amino_acid_to[i + 1] %in% compatible
        if (entered && resolved) {
          waits <- c(waits, g$time[i + 1] - g$time[i])
          if (g$pair_amino_acid_from[i] != g$pair_amino_acid_to[i + 1]) {
            toggles <- toggles + 1
          } else {
            reversions <- reversions + 1
          }
        }
      }
    }
  }
  list(toggles = toggles, reversions = reversions, waits = waits)
}

excursions <- count_toggles(ns)
excursions$toggles

## [1] 104

excursions$reversions

## [1] 99

median(excursions$waits)

## [1] 0.0002205953

Completed toggles and reversions occur in the expected 1:1 ratio, and the median wait in the intermediate is well under one percent of a typical branch length – the compensating substitution follows essentially immediately.

Substitution history

Plotting every event would bury the signal in synonymous churn, so plot only the non-synonymous events:

sim_ns <- sim
sim_ns$substitutions <- ns
plot(sim_ns)

Each event draws one dot for the member of the pair that changed: a dot above the branch is a substitution at the first site, a dot below is a substitution at the second (green synonymous, red non-synonymous – here every plotted event is non-synonymous). The compensated toggles are the tight vertical pairs of red dots, one above and one below the branch at almost the same position: one site substitutes and the other follows within a fraction of a percent of the branch length.

The entry–exit structure is visible directly in the filtered substitution table:

head(ns)

node	time	from	to	pair_from	pair_to	pair_amino_acid_from	pair_amino_acid_to	first_synonymous	second_synonymous	first_color	second_color
51	0.0260742	2496	2480	AAG,GAG	AAG,AAG	K,E	K,K	TRUE	FALSE	#16A35E	#E84B2A
51	0.0261711	2480	3456	AAG,AAG	GAG,AAG	K,K	E,K	FALSE	TRUE	#E84B2A	#16A35E
63	0.1372808	2496	2480	AAG,GAG	AAG,AAG	K,E	K,K	TRUE	FALSE	#16A35E	#E84B2A
63	0.1386381	2419	3395	AAA,AAG	GAA,AAG	K,K	E,K	FALSE	TRUE	#E84B2A	#16A35E
64	0.0496161	2434	2418	AAA,GAA	AAA,AAA	K,E	K,K	TRUE	FALSE	#16A35E	#E84B2A
64	0.0497110	2418	3394	AAA,AAA	GAA,AAA	K,K	E,K	FALSE	TRUE	#E84B2A	#16A35E

Alignment

To export FASTA sequence:

as.fasta(sim$alignment, file = "palantir_sim.fa")

Show alignment in a scrollable viewer: