Paracetamol: which torsions move together¶
md_tools.analysis.t_mi measures the mutual information between torsion pairs -- and measures
the ESTIMATOR'S OWN BIAS against a shuffled null, which is the part that makes the number usable.
The plug-in histogram estimator is biased upward. An independent pair does not return zero; it returns a small positive number that reads as a weak coupling. On this system the pair with the LARGEST raw MI turns out to be the one that is entirely bias, which is the whole argument for measuring the null rather than eyeballing the raw value.
numpy only -- this notebook needs no scikit-learn and no rdkit.
import sys, pathlib
sys.path.insert(0, str(pathlib.Path.cwd().parents[1] / "_analysis")) # docs/tutorial/_analysis
import numpy as np
import matplotlib.pyplot as plt
import mdtraj as md
import tutorial_data as td
# Where the runs are. NOT defaulted: see docs/tutorial/_analysis/tutorial_data.py.
# The resolved path is deliberately NOT printed. These notebooks are committed WITH
# their outputs, so an echoed absolute path would put one machine's layout into a
# published file -- which this project forbids, and which executing a notebook turns
# from a parameter into a stored string.
print("tutorial runs:", "found" if td.runs_root().is_dir() else "MISSING")
tutorial runs: found
from md_tools.analysis import t_mi
PARA_TORSIONS = {"omega": (0, 1, 3, 4), "aryl": (1, 3, 4, 5), "phenol": (6, 7, 8, 17)}
names, quads = list(PARA_TORSIONS), list(PARA_TORSIONS.values())
cmd = td.solute("PARA-1us", stride=25)
tors = md.compute_dihedrals(cmd, quads)
report = t_mi(tors, names=names, bins=24, n_null=8)
for pair in report["pairs"]:
print(f"{pair['a']:7s}-{pair['b']:7s} raw {pair['mi_nats']:.4f} "
f"debiased {pair['mi_nats_debiased']:+.4f} "
f"excess {pair['excess_over_null_in_sd']:+6.1f} sd")
omega -aryl raw 0.0334 debiased +0.0206 excess +9.6 sd aryl -phenol raw 0.0462 debiased +0.0013 excess +0.4 sd omega -phenol raw 0.0117 debiased -0.0027 excess -1.2 sd
Reading the three rows¶
Sort them by RAW MI and the ranking is wrong. The honest column is the excess over the null, measured per pair by shuffling one axis -- which destroys the dependence while keeping both marginals and the binning exactly as they were.
A large excess is evidence of coupling. A small one is not evidence of independence: an underpowered comparison finds nothing, and so does a genuinely independent pair. The module deliberately has no verdict field for this reason -- a boolean would collapse +10 sd and +1 sd into the same answer, and +1 sd is exactly where a reader needs the number.
fig, ax = plt.subplots(figsize=(7.6, 4.0), layout="constrained")
labels = [f"{p['a']}-{p['b']}" for p in report["pairs"]]
raw = [p["mi_nats"] for p in report["pairs"]]
deb = [p["mi_nats_debiased"] for p in report["pairs"]]
pos = np.arange(len(labels))
ax.bar(pos - 0.2, raw, 0.4, label="raw (plug-in)", color="0.65")
ax.bar(pos + 0.2, deb, 0.4, label="debiased", color="#C1501A")
for i, p in enumerate(report["pairs"]):
ax.annotate(f"{p['excess_over_null_in_sd']:+.1f} sd", (i + 0.2, deb[i]),
ha="center", va="bottom", fontsize=9)
ax.axhline(0, color="k", lw=0.8)
ax.set_xticks(pos); ax.set_xticklabels(labels); ax.set_ylabel("MI (nats)")
ax.set_title("Raw MI overstates every pair; the null says by how much")
ax.legend()
<matplotlib.legend.Legend at 0x7f97383a2120>
Doing it over the whole set¶
With more torsions the question becomes WHICH pairs are coupled, and that is a different problem statistically: d torsions means d(d-1)/2 simultaneous comparisons, so a handful of "significant" pairs is the expected yield of pure noise.
analyse handles that: it sweeps the binning and the origin, bootstraps over BLOCKS rather than
frames, and applies a Benjamini-Hochberg false-discovery-rate correction across the pairs. Its
verdicts are three-valued and none of them says "independent" -- a pair where no dependence
was detected at this resolution and this sample size has not been shown independent.
dependence_graph then gives the connected components over supported pairs. A component is a
CANDIDATE analysis group, not proof of irreducible many-body coupling, and a missing edge rules
nothing out: pairwise MI is blind to dependence that appears only in three or more variables
jointly.
from md_tools.analysis._t_mi import analyse, dependence_graph
# `analyse` is the orchestrator: it sweeps the binning and the origin, bootstraps over BLOCKS
# (a correlated trajectory has far fewer independent samples than frames), and applies a
# Benjamini-Hochberg correction across the pairs before assigning a verdict.
result = analyse(tors, names=names)
for pair in result.pairs:
print(f"{pair['name_i']:7s}-{pair['name_j']:7s} {pair['verdict']}")
print("\ndependence graph:", dependence_graph(result))
omega -aryl supported_dependence omega -phenol no_dependence_detected_at_this_resolution aryl -phenol no_dependence_detected_at_this_resolution dependence graph: [['aryl', 'omega']]