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

The stationary-distribution and generic-transient tests use fitness profiles shipped with the package (psi_c50_1). The dS tests additionally construct a standard Zhang–Nielsen–Yang branch-site model and a mutation–selection fitness and population-size shift.

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 event-count scaling

Expectation. A model built with scaling_type = "synonymous-per-codon", "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-per-codon", "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-per-codon` = 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-per-codon	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: dS calibration under homogeneous and heterogeneous models

Homogeneous neutral mutation–selection model

The default scaling_type = "synonymous" is a dS gauge rather than an event-count gauge. With amino-acid fitness differences removed, synonymous histories divided by the model’s neutral synonymous opportunity count should recover dS = 1. The neutral total rate is also three substitutions per codon per dS unit.

ms_ds <- MutationSelection(1000, 1e-8, hky, rep(0, 20),
                           scaling_type = "synonymous")
set_palantir_seed(params$seed + 1L)
sim_ds <- simulate_over_phylogeny(p2, ms_ds, sample_sequence(ms_ds, n_sites))
per_site_total <- tabulate(sim_ds$substitutions$site + 1L,
                           nbins = n_sites) / 2
per_site_syn <- tabulate(sim_ds$substitutions$site[sim_ds$substitutions$synonymous] + 1L,
                         nbins = n_sites) / 2
ds_by_site <- per_site_syn / ms_ds$synonymous_opportunities
ds_checks <- data.frame(
    quantity = c("dS from synonymous histories", "neutral events per codon"),
    expected = c(1, 3),
    observed = c(mean(ds_by_site), mean(per_site_total)),
    se = c(sd(ds_by_site), sd(per_site_total)) / sqrt(n_sites))
ds_checks$z <- (ds_checks$observed - ds_checks$expected) / ds_checks$se
knitr::kable(ds_checks, digits = 4)

quantity	expected	observed	se	z
dS from synonymous histories	1	1.0032	0.0155	0.2033
neutral events per codon	3	3.0005	0.0223	0.0224

Heterogeneous GY94 and mutation–selection trees

Expectation. Under the dS gauge, branch lengths can be recovered from complete synonymous histories after dividing by the appropriate synonymous opportunity exposure. This should hold for (i) all branches and site classes of a standard GY94 branch-site model and (ii) a mutation–selection branch entered after a simultaneous fitness-profile and population-size shift. The second expectation is computed independently by Runge–Kutta integration of the selected master equation and the neutral dS-clock functional.

Criterion. Each input branch length is within four Monte-Carlo standard errors of its history estimate. Four standard errors are used because the GY94 check evaluates six branches simultaneously.

n_ds <- max(12000L, 4L * ceiling(n_sites / 4L))

ds_tree_file <- tempfile(fileext = ".newick")
ds_mode_file <- tempfile(fileext = ".newick")
writeLines("((A:0.15,B:0.30):0.20,(C:0.45,D:0.60):0.25);", ds_tree_file)
writeLines("((A:0,B:1):0,(C:0,D:0):0);", ds_mode_file)
ds_tree <- Phylogeny(ds_tree_file)
ds_modes <- Phylogeny(ds_mode_file, type = "mode")
ds_branches <- phylogeny_to_intervals(ds_tree, ds_modes)
ds_branches$input_dS <- ds_branches$to - ds_branches$from

set_palantir_seed(params$seed + 101L)
zny <- ZNYBranchSiteModel(
    mode_phylogeny = ds_modes,
    site_counts = rep(n_ds / 4L, 4L),
    omega0 = 0.2,
    omega2 = 3.0,
    kappa = 2.5,
    frequencies = c(T = .20, C = .30, A = .25, G = .25),
    frequency_model = "F1x4",
    scaling_type = "dS")
zny_models <- unlist(lapply(zny$models, `[[`, "models"), recursive = FALSE)
zny_NS <- vapply(zny_models, `[[`, numeric(1), "synonymous_opportunities")
stopifnot(max(zny_NS) - min(zny_NS) < 1e-12)
zny_sim <- simulate_gy94_site_model(ds_tree, zny, rescale_method = "exact")

zny_dS <- do.call(rbind, lapply(seq_len(nrow(ds_branches)), function(i) {
    node_i <- ds_branches$node[i]
    keep <- zny_sim$substitutions$node == node_i &
            zny_sim$substitutions$synonymous
    by_site <- tabulate(zny_sim$substitutions$site[keep] + 1L, nbins = n_ds)
    values <- by_site / zny_NS[1]
    estimate <- mean(values)
    se <- sd(values) / sqrt(n_ds)
    data.frame(model = "GY94 branch-site", node = node_i,
               mode = ds_branches$state[i],
               input_dS = ds_branches$input_dS[i],
               observed_dS = estimate, se = se,
               z = (estimate - ds_branches$input_dS[i]) / se)
}))
knitr::kable(zny_dS, digits = 4)

model	node	mode	input_dS	observed_dS	se	z
GY94 branch-site	1	0	0.20	0.1972	0.0044	-0.6284
GY94 branch-site	2	0	0.15	0.1485	0.0038	-0.3962
GY94 branch-site	3	1	0.30	0.2942	0.0053	-1.0924
GY94 branch-site	4	0	0.25	0.2515	0.0050	0.2917
GY94 branch-site	5	0	0.45	0.4542	0.0068	0.6195
GY94 branch-site	6	0	0.60	0.5924	0.0078	-0.9735

The mutation–selection test uses a two-tip tree so that the shifted branch is entered directly from the ancestral equilibrium. Its independently integrated synonymous exposure is reported with the simulation estimate.

nuc <- c(T = .22, C = .28, A = .31, G = .19)
hky_ds <- HasegawaKishinoYano(nuc, transition_rate = 2.5,
                              transversion_rate = 1)
aa_order <- strsplit("ACDEFGHIKLMNPQRSTVWY", "")[[1]]
freq0 <- freq1 <- setNames(rep(.20 / 19, 20), aa_order)
freq0["W"] <- .80
freq1["A"] <- .80
Ne0 <- 500; Ne1 <- 20000
fit0 <- equilibrium_to_fitness(as.numeric(freq0), Ne0,
                               nucleotide_equilibrium = nuc)
fit1 <- equilibrium_to_fitness(as.numeric(freq1), Ne1,
                               nucleotide_equilibrium = nuc)
ms0 <- MutationSelection(Ne0, .01, hky_ds, fit0, scaling_type = "dS")
ms1 <- MutationSelection(Ne1, .01, hky_ds, fit1, scaling_type = "dS")

shift_tree_file <- tempfile(fileext = ".newick")
shift_mode_file <- tempfile(fileext = ".newick")
writeLines("(A:0.30,B:0.10);", shift_tree_file)
writeLines("(A:1,B:0);", shift_mode_file)
shift_tree <- Phylogeny(shift_tree_file)
shift_modes <- Phylogeny(shift_mode_file, type = "mode")
shift_length <- 0.30

p0_ds <- as.numeric(ms0$equilibrium)
h1_ds <- as.numeric(ms1$dS_outflux)
g1_syn <- syn_outflux(ms1)
target <- shift_length * sum(as.numeric(ms1$equilibrium) * h1_ds)
integrated <- rk4_to_budget(p0_ds, ms1$transition, h1_ds, g1_syn, target)
mean_syn_opportunity <- integrated$syn_exposure / shift_length

set_palantir_seed(params$seed + 201L)
shift_root <- sample_sequence(ms0, n_ds)
shift_sim <- simulate_over_interval_phylogeny(
    phylogeny = shift_tree,
    mode_phylogeny = shift_modes,
    models = list(ms0, ms1),
    sequence = shift_root,
    start_mode = 0,
    rescale_method = "exact")
shift_node <- shift_sim$intervals$node[shift_sim$intervals$state == 1]
shift_events <- shift_sim$substitutions[
    shift_sim$substitutions$node == shift_node &
    shift_sim$substitutions$synonymous, ]
shift_by_site <- tabulate(shift_events$site + 1L, nbins = n_ds)
shift_dS_by_site <- shift_by_site / mean_syn_opportunity
shift_estimate <- mean(shift_dS_by_site)
shift_se <- sd(shift_dS_by_site) / sqrt(n_ds)
shift_z <- (shift_estimate - shift_length) / shift_se
mutsel_dS <- data.frame(
    model = "Mutation-selection fitness + Ne shift",
    input_dS = shift_length,
    observed_dS = shift_estimate,
    se = shift_se,
    z = shift_z,
    intrinsic_time = integrated$tau,
    synonymous_opportunities_per_codon_per_dS = mean_syn_opportunity)
knitr::kable(mutsel_dS, digits = 5)

model	input_dS	observed_dS	se	z	intrinsic_time	synonymous_opportunities_per_codon_per_dS
Mutation-selection fitness + Ne shift	0.3	0.28654	0.00931	-1.44582	0.40184	0.32631

all_ds_z_max <- max(abs(c(ds_checks$z, zny_dS$z, shift_z)))
verdict("Test 2", all_ds_z_max < 4,
        sprintf("homogeneous checks and all history estimates recover the dS scale (|z| max %.2f, criterion < 4)",
                all_ds_z_max))

Test 2: PASS — homogeneous checks and all history estimates recover the dS scale (|z| max 1.45, criterion < 4)

Test 3: 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, with both Holm-adjusted p-values greater than 0.01. The adjustment controls the family-wise false-positive rate across the codon- and amino-acid-level comparisons.

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)
stationary_p <- p.adjust(c(codon = chi_codon$p.value,
                           amino_acid = chi_aa$p.value), method = "holm")

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 3", all(stationary_p > 0.01),
        sprintf("raw chi-square p = %.3f (codon), %.3f (amino acid); minimum Holm-adjusted p = %.3f; TV(observed, stationary) = %.4f at the amino-acid level; criterion adjusted p > 0.01",
                chi_codon$p.value, chi_aa$p.value, min(stationary_p),
                tv(obs_aa / sum(obs_aa), th_aa)))

Test 3: PASS — raw chi-square p = 0.408 (codon), 0.302 (amino acid); minimum Holm-adjusted p = 0.603; TV(observed, stationary) = 0.0189 at the amino-acid level; criterion adjusted p > 0.01

Test 4: 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 no Monte-Carlo chi-square comparison of the observed amino-acid distribution against the predicted trajectory is rejected at family-wise alpha 0.01 after Holm adjustment across the twelve checkpoints.

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")
t4 <- rbind(res_seg$table, res_ex$table)
t4$p_holm <- p.adjust(t4$p, method = "holm")
knitr::kable(t4, digits = 4)

method	u	cum_class	p	tv_obs_pred	tv_pred_eqB	p_holm
segments	0.25	0.2393	0.8044	0.0198	0.8524	1.0000
segments	0.50	0.4873	0.8540	0.0198	0.7805	1.0000
segments	1.00	0.9997	0.8340	0.0072	0.6777	1.0000
segments	2.00	2.0400	0.2460	0.0165	0.4595	1.0000
segments	4.00	4.0340	0.2667	0.0214	0.1908	1.0000
segments	8.00	8.0170	0.0016	0.0293	0.0648	0.0192
exact	0.25	0.2313	0.1298	0.0244	0.8524	1.0000
exact	0.50	0.4817	0.0920	0.0236	0.7805	0.8398
exact	1.00	0.9947	0.0420	0.0181	0.6777	0.4619
exact	2.00	2.0360	0.8254	0.0066	0.4595	1.0000
exact	4.00	4.0560	0.5859	0.0091	0.1908	1.0000
exact	8.00	8.0117	0.0840	0.0185	0.0648	0.8398

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(t4$u / n_sites)
calib_ok <- all(abs(t4$cum_class - t4$u) < calib_tol)
p_ok <- all(t4$p_holm > 0.01)
verdict("Test 4", calib_ok && p_ok,
        sprintf("max |cumulative class events - u| / tolerance = %.2f (criterion < 1); raw chi-square p range %.3f-%.3f, minimum Holm-adjusted p %.3f (criterion > 0.01)",
                max(abs(t4$cum_class - t4$u) / calib_tol), min(t4$p),
                max(t4$p), min(t4$p_holm)))

Test 4: PASS — max |cumulative class events - u| / tolerance = 0.34 (criterion < 1); raw chi-square p range 0.002-0.854, minimum Holm-adjusted p 0.019 (criterion > 0.01)

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
Test 4: PASS

cat(sprintf("\n**Overall: %s** (PalantiR %s, seed %d, %d sites in standard tests; %d in tree-level dS tests)\n",
    ifelse(all(unlist(verdicts)), "PASS", "FAIL"),
    as.character(packageVersion("PalantiR")), params$seed, n_sites, n_ds))

Overall: PASS (PalantiR 1.3.0, seed 20260821, 3000 sites in standard tests; 12000 in tree-level dS tests)