This notebook tests the core simulation engine against exact expectations. Each test states its expectation, runs a fresh simulation, compares the result statistically and visually, and prints a PASS or FAIL verdict against a criterion stated in advance. To rerun everything with a different seed or simulation size:

rmarkdown::render("engine_validation.Rmd", params = list(seed = 1, n_sites = 3000))

All tests use the fitness profiles shipped with the package (psi_c50_1, amino-acid equilibrium frequencies converted to fitness with equilibrium_to_fitness) and the mutation–selection codon model.

data(psi_c50_1)
hky <- HasegawaKishinoYano(equilibrium = rep(0.25, 4))
psi_1 <- equilibrium_to_fitness(as.numeric(psi_c50_1[1, ]), 1000)
psi_2 <- equilibrium_to_fitness(as.numeric(psi_c50_1[2, ]), 1000)
n_sites <- params$n_sites

Test 1: branch-length calibration under the three scaling modes

Expectation. A model built with scaling_type = "synonymous", "non-synonymous", or "substitution" is normalized so that one unit of branch length delivers, at stationarity, one expected substitution of that class per codon site: synonymous, non-synonymous, or of any kind, respectively. Simulating 3000 sites over a two-taxon tree with branch lengths 1, the mean number of class events per site per unit branch length should equal 1.

Criterion. For each scaling mode, the observed mean is within three standard errors of 1.

tree2 <- tempfile(fileext = ".newick"); writeLines("(A:1.0,B:1.0);", tree2)
p2 <- Phylogeny(tree2)

t1 <- do.call(rbind, lapply(c("synonymous", "non-synonymous", "substitution"),
  function(sc) {
    ms <- MutationSelection(1000, 1e-8, hky, psi_1, scaling_type = sc)
    sim <- simulate_over_phylogeny(p2, ms, sample_sequence(ms, n_sites))
    sub <- sim$substitutions
    in_class <- switch(sc,
        synonymous       = sub$synonymous,
        "non-synonymous" = !sub$synonymous,
        substitution     = rep(TRUE, nrow(sub)))
    per_site <- tabulate(sub$site[in_class] + 1L, nbins = n_sites) / 2
    data.frame(scaling = sc, mean = mean(per_site),
               se = sd(per_site) / sqrt(n_sites),
               z = (mean(per_site) - 1) / (sd(per_site) / sqrt(n_sites)))
}))
knitr::kable(t1, digits = 4)

scaling	mean	se	z
synonymous	1.0042	0.0140	0.2975
non-synonymous	1.0155	0.0179	0.8672
substitution	1.0115	0.0141	0.8143

verdict("Test 1", all(abs(t1$z) < 3),
        sprintf("class events per site per unit length: %s (|z| max %.2f, criterion < 3)",
                paste(sprintf("%.3f", t1$mean), collapse = ", "), max(abs(t1$z))))

Test 1: PASS — class events per site per unit length: 1.004, 1.015, 1.012 (|z| max 0.87, criterion < 3)

Test 2: stationary distribution

Expectation. The model’s stationary distribution (model$equilibrium) is both the distribution the root sequence is sampled from and the invariant distribution of the substitution process, so the states observed at the end of a branch of any length are draws from it. Comparing the codon frequencies at a tip across 3000 independent sites against the computed stationary distribution tests the equilibrium computation, the root sampler, and the propagation together.

Criterion. Monte-Carlo chi-square goodness of fit, observed tip codon counts against the stationary distribution, p > 0.01 at both the codon and the amino-acid level.

ms <- MutationSelection(1000, 1e-8, hky, psi_1)
sim <- simulate_over_phylogeny(p2, ms, sample_sequence(ms, n_sites))
codons <- names(GeneticCode())
obs_codon <- as.numeric(table(factor(sim$alignment["A", ], levels = codons)))
pi_codon <- as.numeric(ms$equilibrium)

chi_codon <- chisq.test(obs_codon, p = pi_codon, rescale.p = TRUE,
                        simulate.p.value = TRUE, B = 5000)
obs_aa <- codon_to_aa(obs_codon); th_aa <- codon_to_aa(pi_codon)
chi_aa <- chisq.test(obs_aa, p = th_aa, rescale.p = TRUE,
                     simulate.p.value = TRUE, B = 5000)

The stationary and observed amino-acid distributions as sequence logos (letter height = frequency):

plot_freqs(data.frame(stationary = th_aa, observed = obs_aa / sum(obs_aa)),
           c("stationary", "observed"))

verdict("Test 2", chi_codon$p.value > 0.01 && chi_aa$p.value > 0.01,
        sprintf("chi-square p = %.3f (codon), %.3f (amino acid); TV(observed, stationary) = %.4f at the amino-acid level; criterion p > 0.01",
                chi_codon$p.value, chi_aa$p.value, tv(obs_aa / sum(obs_aa), th_aa)))

Test 2: PASS — chi-square p = 0.469 (codon), 0.340 (amino acid); TV(observed, stationary) = 0.0197 at the amino-acid level; criterion p > 0.01

Test 3: the transient between two equilibria

Expectation. In a time-heterogeneous simulation the sequence starts at the stationary distribution of the first model and relaxes to the stationary distribution of the second. The state distribution along the way follows the master equation of the second model’s generator, expressed in the branch’s own units: branch position is measured in expected substitutions of the scaled class, so the predicted trajectory is obtained by integrating $dp/du = pQ / g(p)$, where $g(p)$ is the class flux under the current state distribution. This is a demanding construction: the two fitness profiles used here are nearly disjoint, so the class flux starts at several hundred times its equilibrium value and the branch begins with a burst of substitutions. (This test found the defect fixed in version 1.2.1: the exact rescaler’s knot table could not resolve that burst and misplaced early event positions.)

Criterion. Under both rescale methods ("segments" and "exact"): cumulative class events per site track branch position to within $0.02 + 4\sqrt{u/n}$ at every checkpoint (the count’s sampling standard deviation grows with $u$, so the bound scales with the checkpoint; the 1.2.1 defect exceeded the early-checkpoint bounds five-fold), and Monte-Carlo chi-square of the observed amino-acid distribution against the predicted trajectory gives p ≥ 0.001 at every checkpoint with at most two of the twelve p-values below 0.05.

mA <- MutationSelection(1000, 1e-8, hky, psi_1, scaling_type = "non-synonymous")
mB <- MutationSelection(1000, 1e-8, hky, psi_2, scaling_type = "non-synonymous")
Q <- mB$transition
p0 <- as.numeric(mA$equilibrium)
piB <- as.numeric(mB$equilibrium)
W <- Q; W[outer(aa_of, aa_of, "==")] <- 0; diag(W) <- 0
w <- rowSums(W)                       # non-synonymous outflux per state

# Predicted state distribution as a function of branch position u:
# fourth-order Runge-Kutta on dp/du = pQ / g(p).
predict_traj <- function(p0, Q, w, L, du = 0.002) {
    f <- function(p) as.numeric(p %*% Q) / sum(p * w)
    n <- ceiling(L / du)
    P <- matrix(NA, n + 1, length(p0)); P[1, ] <- p0
    p <- p0
    for (k in 1:n) {
        k1 <- f(p); k2 <- f(p + du/2 * k1); k3 <- f(p + du/2 * k2); k4 <- f(p + du * k3)
        p <- p + du/6 * (k1 + 2*k2 + 2*k3 + k4)
        p[p < 0] <- 0; p <- p / sum(p)
        P[k + 1, ] <- p
    }
    list(u = seq(0, n * du, by = du), P = P)
}
L <- 8
traj <- predict_traj(p0, Q, w, L)
pred_at <- function(u) traj$P[which.min(abs(traj$u - u)), ]

treeT <- tempfile(fileext = ".newick"); writeLines(sprintf("(A:%f,B:0.01);", L), treeT)
modeT <- tempfile(fileext = ".newick"); writeLines("(A:1,B:0);", modeT)
pT <- Phylogeny(treeT); mT <- Phylogeny(modeT, type = "mode")

checkpoints <- c(0.25, 0.5, 1, 2, 4, 8)
run_method <- function(method) {
    root <- sample_sequence(mA, n_sites)
    sim <- simulate_over_interval_phylogeny(pT, mT, list(mA, mB), root,
                                            start_mode = 0, rescale_method = method)
    nodeA <- sim$intervals$node[sim$intervals$state == 1]
    subA <- sim$substitutions[sim$substitutions$node %in% nodeA, ]
    subA <- subA[order(subA$site, subA$time), ]
    ns_time <- subA$time[aa_of[subA$from + 1] != aa_of[subA$to + 1]]
    state_at <- function(u) {
        st <- root$index[, 1]
        sel <- subA[subA$time <= u, ]
        if (nrow(sel)) {
            last <- sel[!duplicated(sel$site, fromLast = TRUE), ]
            st[last$site + 1] <- last$to
        }
        st
    }
    rows <- lapply(checkpoints, function(u) {
        emp <- tabulate(state_at(u) + 1L, nbins = length(p0))
        pr <- pred_at(u)
        emp_aa <- codon_to_aa(emp); pr_aa <- codon_to_aa(pr)
        keep <- pr_aa * n_sites >= 5
        chi <- chisq.test(c(emp_aa[keep], sum(emp_aa[!keep])),
                          p = c(pr_aa[keep], sum(pr_aa[!keep])),
                          rescale.p = TRUE, simulate.p.value = TRUE, B = 5000)
        data.frame(method = method, u = u,
                   cum_class = sum(ns_time <= u) / n_sites,
                   p = chi$p.value,
                   tv_obs_pred = tv(emp_aa / n_sites, pr_aa),
                   tv_pred_eqB = tv(pr_aa, codon_to_aa(piB)))
    })
    list(table = do.call(rbind, rows), state_at = state_at)
}
res_seg <- run_method("segments")
res_ex  <- run_method("exact")
t3 <- rbind(res_seg$table, res_ex$table)
knitr::kable(t3, digits = 4)

method	u	cum_class	p	tv_obs_pred	tv_pred_eqB
segments	0.25	0.2407	0.8724	0.0172	0.8524
segments	0.50	0.4907	0.3697	0.0246	0.7805
segments	1.00	0.9920	0.1578	0.0136	0.6777
segments	2.00	2.0033	0.1490	0.0202	0.4595
segments	4.00	3.9630	0.2448	0.0182	0.1908
segments	8.00	7.9927	0.0070	0.0233	0.0648
exact	0.25	0.2440	0.3127	0.0256	0.8524
exact	0.50	0.4913	0.1602	0.0326	0.7805
exact	1.00	1.0010	0.4619	0.0123	0.6777
exact	2.00	2.0157	0.8138	0.0077	0.4595
exact	4.00	4.0160	0.3585	0.0189	0.1908
exact	8.00	7.9980	0.3677	0.0102	0.0648

The distance from the new equilibrium decays along the branch, and the observed distributions (points) follow the predicted trajectory (line) under both methods:

u_dense <- seq(0, L, by = 0.1)
tv_pred <- sapply(u_dense, function(u) tv(codon_to_aa(pred_at(u)), codon_to_aa(piB)))
plot(u_dense, tv_pred, type = "l", lwd = 2,
     xlab = "Branch position (expected non-synonymous substitutions per site)",
     ylab = "Total variation distance from the new equilibrium")
u_pts <- seq(0.1, L, by = 0.26)
for (r in list(list(res_seg, "#E84B2A", 19), list(res_ex, "#2166AC", 1))) {
    tv_obs <- sapply(u_pts, function(u)
        tv(codon_to_aa(tabulate(r[[1]]$state_at(u) + 1L, nbins = length(p0))) / n_sites,
           codon_to_aa(piB)))
    points(u_pts, tv_obs, col = r[[2]], pch = r[[3]])
}
legend("topright", bty = "n",
       c("predicted trajectory", "observed (segments)", "observed (exact)"),
       lwd = c(2, NA, NA), pch = c(NA, 19, 1), col = c("black", "#E84B2A", "#2166AC"))

Amino-acid distributions at the start, in mid-transient, and at the end of the branch — predicted against observed (segments method):

plot_freqs(data.frame(
    start_pred = codon_to_aa(p0),
    u0.5_pred  = codon_to_aa(pred_at(0.5)),
    u0.5_obs   = codon_to_aa(tabulate(res_seg$state_at(0.5) + 1L, nbins = length(p0))) / n_sites,
    u8_pred    = codon_to_aa(pred_at(8)),
    u8_obs     = codon_to_aa(tabulate(res_seg$state_at(8) + 1L, nbins = length(p0))) / n_sites,
    eq_B       = codon_to_aa(piB)),
    c("start", "u=0.5 pred", "u=0.5 obs", "u=8 pred", "u=8 obs", "equilibrium B"))

calib_tol <- 0.02 + 4 * sqrt(t3$u / n_sites)
calib_ok <- all(abs(t3$cum_class - t3$u) < calib_tol)
p_ok <- min(t3$p) >= 0.001 && sum(t3$p < 0.05) <= 2
verdict("Test 3", calib_ok && p_ok,
        sprintf("max |cumulative class events - u| / tolerance = %.2f (criterion < 1); chi-square p range %.3f-%.3f with %d of %d below 0.05 (criterion: min >= 0.001, at most 2)",
                max(abs(t3$cum_class - t3$u) / calib_tol), min(t3$p), max(t3$p), sum(t3$p < 0.05), nrow(t3)))

Test 3: PASS — max |cumulative class events - u| / tolerance = 0.22 (criterion < 1); chi-square p range 0.007-0.872 with 1 of 12 below 0.05 (criterion: min >= 0.001, at most 2)

Summary

for (id in names(verdicts)) {
    cat(sprintf("- **%s: %s**\n", id, ifelse(verdicts[[id]], "PASS", "FAIL")))
}

Test 1: PASS
Test 2: PASS
Test 3: PASS

cat(sprintf("\n**Overall: %s** (PalantiR %s, seed %d, %d sites per test)\n",
    ifelse(all(unlist(verdicts)), "PASS", "FAIL"),
    as.character(packageVersion("PalantiR")), params$seed, n_sites))

Overall: PASS (PalantiR 1.2.1, seed 20260821, 3000 sites per test)

PalantiR engine validation

Test 1: branch-length calibration under the three scaling modes

Test 2: stationary distribution

Test 3: the transient between two equilibria

Summary