Algorithms and Learning for Protein Science

Contents

• Introduction and overview
• 0. Review: rigid motions, Boltzmann distributions, Markov chains, polytopes and attention
• 1. Proteins at the atomic level: structures, mechanisms, free energies, time scales
• 2. Structure prediction with AlphaFold2 and transformers
• 3. Molecular kinematics, inverse problems and loop sampling
• 4. Boltzmann samplers for sequences and structures
• 5. Computer practical: molecular manipulations
• 6. High-dimensional sampling: volumes of polytopes and densities of states
• 7. Spatial partitions and applications
• 8. Complements for self-study
• 9. Research directions and open problems
• Exercises
• Projects and assessment
• References

Introduction and overview

Audience: a student with an M1 in mathematics, computer science or physics who knows Markov chains, linear algebra and basic deep learning, and who wants to work on proteins without a degree in biochemistry.

Empirical first: AlphaFold2 predicts one structure per sequence, yet a membrane transporter only works by cycling through several conformations. A Metropolis chain on a double well with a barrier of 8\,kT crosses it six times in 2 \times 10^5 steps and misestimates the population of the second well by 30 % (Code 0.9). This course explains both facts and the algorithms that address them.

The course treats a protein as a point in a space of dimension about 15,000 that moves under thermal noise. It covers the prediction of structures from sequences, the kinematics of the backbone, generative models of sequences, and the sampling of high-dimensional spaces. It uses tools from Computational Statistics, Stochastic Calculus in Machine Learning, Geometric Data Analysis and Deep Learning.

What you will be able to do

• Describe a protein structure in Cartesian and internal coordinates, and compare two structures with the RMSD after optimal superposition, GDT, lDDT and the TM-score.
• Distinguish energies from free energies, and compute the population of a state from its energy and its width.
• Explain metastability and local ergodicity, and estimate an escape time with the Eyring–Kramers formula.
• Align sequences by dynamic programming, and explain how AlphaFold2 turns a multiple sequence alignment into a structure.
• Prove that invariant point attention is invariant under global rigid motions.
• Build a backbone from dihedral angles, count the degrees of freedom of a closed loop, and sample loop conformations.
• Fit mixtures of von Mises distributions to dihedral angles on the flat torus.
• Derive the Potts model from maximum entropy, fit it by pseudo-likelihood, and relate it to self-attention.
• Sample uniformly from a polytope by hit-and-run and estimate its volume by multiphase Monte Carlo.
• Compute a density of states with the Wang–Landau algorithm.
• Search nearest neighbours with random projection trees, and attribute a prediction to features with Shapley values, exactly for decision trees.

Official objectives and syllabus (2026–2027)

The course description on master-mva.com starts from a fact: proteins underlie all biological functions, yet their mechanisms at the atomic scale remain a fundamental open problem. With about 5000 atoms and three coordinates per atom, a polypeptide chain of median size lives in a configuration space of dimension 15,000.

AlphaFold2 provides plausible structures of well-folded regions, but it gives no insight into dynamics. The goal of the class is therefore twofold. First, it casts the main problems of protein dynamics into a rigorous mathematical and algorithmic framework. Second, it presents major ongoing developments at the interface of theoretical biophysics, geometry, topology and machine learning.

The organisation is six lectures of three hours plus one computer practical, which introduces standard molecular manipulations on real data and illustrates selected methods. Each announced item is matched below.

Beyond the syllabus, this page adds a review chapter (ch. 0), complements with longer proofs (ch. 8), a research chapter (ch. 9), graded exercises, project subjects and a written-exam proposal.

Where this course fits in the programme

Official prerequisites: training in algorithms or machine learning, and an interest in biophysics, biology or medicine (master-mva.com). The instructor listed on master-mva.com is Frédéric Cazals (Inria). The course carries the label of the Santé track and is taught in semester 2.

• Refreshers — Probability and Statistics for Markov chains and maximum likelihood; Linear Algebra, Analysis and Optimization for the SVD, rotations and linear programs; Scientific Python, NumPy and PyTorch for the code.
• Semester 1 courses that prepare it — Computational Statistics (MCMC), Stochastic Calculus in Machine Learning (Langevin dynamics), Geometric Data Analysis (distances, equivariance), Deep Learning and Foundations of Large Language Models (transformers), Introduction to Probabilistic Graphical Models (Markov random fields), Robotics (kinematic chains) and Topological Data Analysis.
• Related courses of semester 2 — Graphical Models: Discrete Inference and Learning (Potts models), Bayesian Machine Learning, Diffusion and Multiscale Models in Artificial Intelligence, Geometry and Shape Spaces, Interactions: Complexity, Stochasticity, Extrema and Rare Events (rare events), Multimodal Explainable AI (Shapley values) and Representation Learning for Computer Vision and Medical Imaging.
• Afterwards — the Research Internship; the Hospital Immersion of the Santé track gives a clinical context to the molecular mechanisms of ch. 1.

Course outline

• Review: rigid motions, Boltzmann distributions, Markov chains, polytopes and attention — the M1 tools used later, with short proofs.
• Proteins at the atomic level — structures, the RMSD and the Kabsch theorem, three molecular mechanisms, energies versus free energies, time scales and local ergodicity, open challenges.
• Structure prediction with AlphaFold2 and transformers — sequence alignment, multiple sequence alignments, the Evoformer and the structure module, invariant point attention, accuracy metrics.
• Molecular kinematics, inverse problems and loop sampling — internal coordinates, forward kinematics, loop closure, sampling on the closure variety, von Mises mixtures on the torus.
• Boltzmann samplers for sequences and structures — maximum entropy, Potts models, pseudo-likelihood, autoregressive models, factored attention, protein design.
• Computer practical — visualisation, structure prediction, distances and alignments, loop sampling.
• High-dimensional sampling — hardness of volume computation, hit-and-run, multiphase Monte Carlo, PDMPs, densities of states and Wang–Landau.
• Spatial partitions and applications — random projection trees, nearest neighbours, regression, dimension estimation, Shapley and Shapley–Shubik values, exact Shapley values of decision trees, mixing assessment.
• Complements for self-study — the 1D Eyring–Kramers formula, pseudo-likelihood consistency, invariance of PDMPs, Wang–Landau as stochastic approximation, tripeptide closure, conditioning on closure, the proofs on Shapley values and on von Mises fits.
• Research directions and open problems — conformational ensembles, generative design, enhanced sampling, high-dimensional volumes (self-study).

Session plan

21 h in class = 7 sessions of 3 h.

The volume is the official one of the Université Paris-Saclay programme table (18 h of lectures and 3 h of labs). master-mva.com announces six lectures of three hours and one computer practical, placed fifth; the plan keeps that order, one session per official item. The split of each session and the content of the practical are a proposal.

Level 1 = foundations, 2 = core methods, 3 = advanced material close to research; the level never decreases from one session to the next.

Self-study, outside the 21 h (about 7 h in total): ch. 0 Review — to read before session 1 (≈ 2 h, with Code 0.9) · ch. 8 Complements — the 1D Eyring–Kramers formula, pseudo-likelihood, PDMP invariance, Wang–Landau as stochastic approximation (≈ 3 h) · ch. 9 Research directions and open problems (≈ 2 h).

How to study this module

• Read ch. 0 before session 1. The trace inequality (Lemma 0.2) gives the Kabsch theorem, and the Gibbs variational principle (Proposition 0.4) gives free energies and Potts models.
• Follow one thread: a Boltzmann distribution \pi \propto e^{-\beta U} is the target of molecular dynamics (ch. 1), of sequence models (ch. 4) and of Wang–Landau (ch. 6). The obstacle is always a high dimension combined with rare transitions.
• After each session, run the code of the chapter and change one parameter (barrier height, number of dihedrals, coupling strength, dimension) to see what breaks.
• Do the exercises of the session plan before the next session. The ★★★ exercises and the projects prepare the tandem project, which carries 15 of the 20 points.
• Python environment: NumPy, SciPy, scikit-learn and PyTorch on CPU; every code block runs in under 30 seconds on synthetic data or on the digits data set of scikit-learn. PyMOL, VMD and AlphaFold2 are described in ch. 5 but are not needed to run the code.

0. Review: rigid motions, Boltzmann distributions, Markov chains, polytopes and attention

The course compares structures up to rigid motions, samples Boltzmann distributions with Markov chains, cuts polytopes with lines and reads transformers. This chapter states these M1 tools in the form used later, with short proofs. It is meant to be read before session 1.

0.1 Rigid motions and the SVD

Every real matrix H \in \mathbb{R}^{3 \times 3} has a singular value decomposition H = U \Sigma V^\top, \qquad U, V \in O(3), \qquad \Sigma = \mathrm{diag}(\sigma_1, \sigma_2, \sigma_3),\ \sigma_1 \geq \sigma_2 \geq \sigma_3 \geq 0 .

0.2 Boltzmann distributions and entropy

The same identity, read with constraints on the averages instead of a given energy, gives the maximum-entropy principle of Jaynes (1957). Chapter 4 uses it to derive the Potts model.

0.3 Markov chains and Monte Carlo

import numpy as np
from scipy.integrate import quad
rng = np.random.default_rng(10)

U = lambda x, h: h * (x**2 - 1) ** 2 + 0.5 * x          # tilted double well, barrier about h
def metropolis(h, n=200_000, step=0.25):
    x, out = -1.0, np.empty(n)
    for k in range(n):
        y = x + step * rng.standard_normal()             # symmetric proposal
        if np.log(rng.random()) < U(x, h) - U(y, h):     # accept with min(1, exp(-(U(y) - U(x))))
            x = y
        out[k] = x
    return out
for h in (1.0, 8.0):
    Z = quad(lambda x: np.exp(-U(x, h)), -3, 3)[0]
    p_right = quad(lambda x: np.exp(-U(x, h)), 0, 3)[0] / Z
    xs = metropolis(h)
    well = np.sign(xs[np.abs(xs) > 0.5])                  # which well, away from the barrier top
    hops = np.count_nonzero(np.diff(well))                # transitions between the two wells
    print(f"barrier {h:4.1f}: exact P(x > 0) = {p_right:.3f}, chain estimate = {np.mean(xs > 0):.3f}, transitions {hops}")

0.4 Polytopes and linear programs

0.5 Attention

1. Proteins at the atomic level: structures, mechanisms, free energies, time scales

A protein works by changing shape. This chapter defines a conformation, the distance between two conformations, the free energy of a state and the time needed to leave it, then lists the open problems of the course.

1.1 Proteins: structures

A protein is a chain of amino acids, of 20 kinds, joined by peptide bonds. The repeated atoms N, Cα and C form the backbone; a side chain hangs from each Cα. The order of the amino acids is the primary structure. Local patterns of hydrogen bonds give the secondary structure (α-helices, β-strands). The fold of the whole chain is the tertiary structure, and the assembly of several chains is the quaternary structure.

Experimental structures (X-ray crystallography, NMR, cryo-electron microscopy) are deposited in the Protein Data Bank (PDB), with atomic coordinates in ångströms (1 Å = 10^{-10} m). Branden and Tooze, and Kuriyan, Konforti and Wemmer, cover this material.

1.2 Comparing two structures

The superposition costs one 3 \times 3 SVD after an O(n) pass over the atoms. Code 2.12 implements it.

1.3 Three molecular mechanisms

In all three cases the function lies in a change of conformation or in the rigidity of one, not in a single static structure.

1.4 Energies versus free energies

1.5 Time scales and local ergodicity

Protein motions span more than fifteen decades of time. Molecular dynamics integrates Newton’s equations with steps of 1 to 2 fs, so a millisecond takes about 10^{12} steps (Adcock and McCammon, 2006).

import numpy as np
rng = np.random.default_rng(0)

def mfpt(h, beta=1.0, dt=2e-3, walkers=300, tmax=600.0):
    """Mean first-passage time from x=-1 to x=+1 for dX = -U'(X) dt + sqrt(2/beta) dW,
    with the double well U(x) = h (x^2 - 1)^2 (barrier height h)."""
    x = -np.ones(walkers)
    t_hit = np.full(walkers, np.nan)
    alive = np.ones(walkers, bool)
    for k in range(int(tmax / dt)):
        xa = x[alive]
        xa += -4 * h * xa * (xa**2 - 1) * dt + np.sqrt(2 * dt / beta) * rng.standard_normal(xa.size)
        x[alive] = xa
        hit = alive.copy(); hit[alive] = xa >= 1.0
        t_hit[hit] = (k + 1) * dt
        alive &= ~hit
        if not alive.any():
            break
    return np.nanmean(t_hit), alive.sum()

for h in [2.0, 4.0, 6.0]:
    kramers = 2 * np.pi / (np.sqrt(8 * h * 4 * h)) * np.exp(h)   # Eyring-Kramers, beta = 1
    tau, left = mfpt(h)
    print(f"barrier {h:.0f} kT: simulated MFPT {tau:6.1f}, Eyring-Kramers {kramers:6.1f}, not crossed {left}")

1.6 Open challenges

• Sampling. Barriers break ergodicity on accessible time scales, so populations and free energies need biased or adaptive samplers (ch. 4, ch. 6).
• Structures to ensembles. A predicted structure is one point of \mathcal{C}; function needs the distribution over metastable states (ch. 2, ch. 3).
• Geometry of constraints. Bonds, angles and closed loops confine the chain to curved subsets of \mathcal{C} (ch. 3).
• High dimension. Volumes, densities of states and nearest neighbours are hard except through randomisation (ch. 6, ch. 7).
• Learning with physics. Models must respect rigid motions, and their predictions must be explained (ch. 2, ch. 7).

2. Structure prediction with AlphaFold2 and transformers

Can the structure be read from the sequence? This chapter follows AlphaFold2 from the alignment of related sequences to pair features, then to a structure built by an attention layer that respects rigid motions, and ends with the metrics of success.

2.1 Protein sequences and multiple sequence alignments

import numpy as np

def needleman_wunsch(a, b, match=2, mismatch=-1, gap=-2):
    """Global alignment by dynamic programming (Needleman-Wunsch, linear gap cost)."""
    n, m = len(a), len(b)
    S = np.zeros((n + 1, m + 1))
    S[:, 0] = gap * np.arange(n + 1)
    S[0, :] = gap * np.arange(m + 1)
    for i in range(1, n + 1):
        for j in range(1, m + 1):
            s = match if a[i - 1] == b[j - 1] else mismatch
            S[i, j] = max(S[i - 1, j - 1] + s, S[i - 1, j] + gap, S[i, j - 1] + gap)
    i, j, top, bot = n, m, "", ""                       # traceback
    while i > 0 or j > 0:
        s = match if i > 0 and j > 0 and a[i - 1] == b[j - 1] else mismatch
        if i > 0 and j > 0 and S[i, j] == S[i - 1, j - 1] + s:
            top, bot, i, j = a[i - 1] + top, b[j - 1] + bot, i - 1, j - 1
        elif i > 0 and S[i, j] == S[i - 1, j] + gap:
            top, bot, i = a[i - 1] + top, "-" + bot, i - 1
        else:
            top, bot, j = "-" + top, b[j - 1] + bot, j - 1
    return S[n, m], top, bot

for a, b in [("WHAT", "WAT"), ("HEAGAWGHEE", "PAWHEAE")]:
    score, top, bot = needleman_wunsch(a, b)
    print(f"score {score:g}\n  {top}\n  {bot}")

Pevsner’s textbook covers database search and alignment heuristics.

2.2 AlphaFold2: architecture

AlphaFold2 (Jumper and coauthors, 2021) takes an MSA, and optionally structural templates, and returns atomic coordinates with a confidence score. It has two parts: a trunk of 48 Evoformer blocks and a structure module of 8 layers with shared weights. The whole network runs three more times on its own outputs (recycling).

import numpy as np
from scipy.spatial.transform import Rotation
rng = np.random.default_rng(9)
N, P, c = 6, 4, 8                                    # residues, points per residue, channels

R = Rotation.random(N, random_state=1).as_matrix()   # residue frames T_i = (R_i, t_i)
t = rng.normal(0, 5, (N, 3))
s = rng.normal(size=(N, c))                          # single representation
Wq, Wk = rng.normal(size=(c, c)), rng.normal(size=(c, c))
Wqp, Wkp, Wvp = (rng.normal(size=(c, 3 * P)) for _ in range(3))
bias = rng.normal(size=(N, N))                       # pair bias b_ij from the pair representation

def ipa(R, t, gamma=0.1):
    """One head of invariant point attention: logits and the point output in local frames."""
    q, k = s @ Wq, s @ Wk
    to_global = lambda W: np.einsum("nab,npb->npa", R, (s @ W).reshape(N, P, 3)) + t[:, None]
    qg, kg, vg = to_global(Wqp), to_global(Wkp), to_global(Wvp)
    dist2 = ((qg[:, None] - kg[None]) ** 2).sum((-1, -2))          # sum_p |T_i q_ip - T_j k_jp|^2
    logits = q @ k.T / np.sqrt(c) + bias - gamma / 2 * dist2
    a = np.exp(logits - logits.max(1, keepdims=True)); a /= a.sum(1, keepdims=True)
    out = np.einsum("ij,jpa->ipa", a, vg)                          # attend in the global frame
    return logits, np.einsum("nba,npb->npa", R, out - t[:, None])  # back to the local frame of i

G = Rotation.random(random_state=2).as_matrix(); g = rng.normal(0, 10, 3)   # a global rigid motion
l1, o1 = ipa(R, t)
l2, o2 = ipa(G @ R, t @ G.T + g)
print("max change of the logits:", np.abs(l1 - l2).max(), " of the local outputs:", np.abs(o1 - o2).max())

Each layer of the structure module updates the frames by composing them with predicted rigid motions, and predicts the side-chain torsion angles. Training uses the frame-aligned point error (FAPE), with true frames \tilde T_i, true positions \tilde x_j, distances in ångströms and a clamp at 10 Å: \mathrm{FAPE} = \frac{1}{10\, N_{\mathrm{frames}} N_{\mathrm{atoms}}} \sum_{i, j} \min\Big(10,\ \big\| T_i^{-1} x_j - \tilde T_i^{-1} \tilde x_j \big\|\Big).

Positions are seen from each residue, so FAPE is invariant under rigid motions but penalises mirror images. The network also predicts a per-residue confidence (pLDDT, from 0 to 100) and a predicted aligned error (PAE) for each pair.

2.3 AlphaFold2: performance assessment

At CASP14 (2020), the blind assessment of structure prediction, AlphaFold2 reached a median backbone accuracy of 0.96 Å (Cα RMSD at 95 % coverage), against 2.8 Å for the next best method (Jumper and coauthors, 2021). Moussad, Roche and Bhattacharya (2023) review how transformers, including protein language models that need no MSA, spread through the field.

import numpy as np
rng = np.random.default_rng(1)

def kabsch(P, Q):
    """Rotation R and translation t minimising sum ||R p_i + t - q_i||^2 (Theorem 1.3)."""
    p0, q0 = P.mean(0), Q.mean(0)
    H = (P - p0).T @ (Q - q0)                  # 3 x 3 cross-covariance
    U, S, Vt = np.linalg.svd(H)
    D = np.diag([1.0, 1.0, np.sign(np.linalg.det(Vt.T @ U.T))])
    R = Vt.T @ D @ U.T
    return R, q0 - R @ p0

def rmsd(P, Q):
    R, t = kabsch(P, Q)
    return np.sqrt(((P @ R.T + t - Q) ** 2).sum(1).mean())

def best_score(P, Q, score, cut, seed=8, iters=10):
    """Heuristic max over superpositions: seed fragments, then refit on residues closer than cut."""
    best = 0.0
    for s in range(0, len(P) - seed + 1, 2):
        idx = np.arange(s, s + seed)
        for _ in range(iters):
            R, t = kabsch(P[idx], Q[idx])
            d = np.linalg.norm(P @ R.T + t - Q, axis=1)
            best = max(best, score(d))
            new = np.flatnonzero(d < cut)
            if len(new) < 3 or np.array_equal(new, idx):
                break
            idx = new
    return best

def helix(n):                                    # C-alpha trace of an ideal alpha-helix (Angstrom)
    k = np.arange(n)
    return np.c_[2.3 * np.cos(np.radians(100 * k)), 2.3 * np.sin(np.radians(100 * k)), 1.5 * k]

L = 60
target = helix(L)
model = target.copy()                            # hinge: rotate residues 40..59 by 40 degrees
a = np.radians(40); Rx = np.array([[1, 0, 0], [0, np.cos(a), -np.sin(a)], [0, np.sin(a), np.cos(a)]])
model[40:] = (model[40:] - model[40]) @ Rx.T + model[40]
model += 0.3 * rng.standard_normal(model.shape)  # thermal-like noise

d0 = 1.24 * (L - 15) ** (1 / 3) - 1.8
tm = best_score(model, target, lambda d: np.mean(1 / (1 + (d / d0) ** 2)), cut=d0)
gdt = np.mean([best_score(model, target, lambda d, c=c: np.mean(d < c), cut=c) for c in (1, 2, 4, 8)])
print(f"global RMSD {rmsd(model, target):.2f} A, RMSD of residues 0-39 {rmsd(model[:40], target[:40]):.2f} A")
print(f"GDT_TS {100 * gdt:.1f}, TM-score {tm:.3f} (d0 = {d0:.2f} A)")

The code prints a global RMSD of 3.98 Å, but 0.43 Å on residues 0–39, with GDT_TS = 79.2 and TM-score = 0.715.

A single hinge inflates the global RMSD, while GDT_TS and the TM-score see that two thirds of the chain are almost exact; CASP therefore ranks models with GDT and lDDT.

AlphaFold2 returns one structure per sequence, its accuracy drops with shallow MSAs, and a low pLDDT often marks disorder. It does not predict the populations of alternative conformations, the subject of ch. 3, ch. 4 and ch. 6.

3. Molecular kinematics, inverse problems and loop sampling

Bond lengths and angles barely move, so a backbone moves mostly by rotations about its bonds. This chapter describes a chain by its dihedral angles, like a robot arm by its joints, then closes and samples loops and fits densities on the torus of angles.

3.1 Internal coordinates

import numpy as np

BOND = {"N-CA": 1.458, "CA-C": 1.525, "C-N": 1.329}            # Angstrom
ANGLE = {"N-CA-C": 111.2, "CA-C-N": 116.2, "C-N-CA": 121.7}     # degrees

def place(a, b, c, bond, angle, torsion):
    """NeRF: position of d with |cd| = bond, angle(b,c,d) = angle, dihedral(a,b,c,d) = torsion."""
    bc = (c - b) / np.linalg.norm(c - b)
    n = np.cross(b - a, bc); n /= np.linalg.norm(n)
    M = np.c_[bc, np.cross(n, bc), n]
    th, to = np.radians(angle), np.radians(torsion)
    return c + M @ (bond * np.array([-np.cos(th), np.sin(th) * np.cos(to), np.sin(th) * np.sin(to)]))

def backbone(phi, psi, omega=None):
    """Backbone N, CA, C of len(phi) residues from internal coordinates (ideal bonds and angles)."""
    k = len(phi); omega = np.full(k, 180.0) if omega is None else omega
    X = [np.zeros(3), np.array([BOND["N-CA"], 0, 0])]            # N0, CA0
    t = np.radians(180 - ANGLE["N-CA-C"])
    X.append(X[1] + BOND["CA-C"] * np.array([np.cos(t), np.sin(t), 0]))   # C0
    for i in range(1, k):
        X.append(place(X[-3], X[-2], X[-1], BOND["C-N"], ANGLE["CA-C-N"], psi[i - 1]))      # N_i
        X.append(place(X[-3], X[-2], X[-1], BOND["N-CA"], ANGLE["C-N-CA"], omega[i - 1]))   # CA_i
        X.append(place(X[-3], X[-2], X[-1], BOND["CA-C"], ANGLE["N-CA-C"], phi[i]))         # C_i
    return np.array(X)

def dihedral(p0, p1, p2, p3):
    b0, b1, b2 = p0 - p1, p2 - p1, p3 - p2
    b1 = b1 / np.linalg.norm(b1)
    v, w = b0 - (b0 @ b1) * b1, b2 - (b2 @ b1) * b1
    return np.degrees(np.arctan2(np.cross(b1, v) @ w, v @ w))

k = 8
phi, psi = np.full(k, -60.0), np.full(k, -45.0)                 # alpha-helix
X = backbone(phi, psi)
phi_back = [dihedral(X[3*i - 1], X[3*i], X[3*i + 1], X[3*i + 2]) for i in range(1, k)]
psi_back = [dihedral(X[3*i], X[3*i + 1], X[3*i + 2], X[3*i + 3]) for i in range(k - 1)]
print("round trip error (deg):", np.abs(np.r_[phi_back - phi[1:], psi_back - psi[:-1]]).max().round(9))
print("CA(i)-CA(i+1) distances:", np.linalg.norm(np.diff(X[1::3], axis=0), axis=1).round(3))
print("rise per residue along the helix (A):", round(np.linalg.norm(X[1::3][-1] - X[1::3][0]) / (k - 1), 2))

def end_atoms(theta, k):                                        # last residue's N, CA, C
    phi = np.r_[0.0, theta[: k - 1]]; psi = np.r_[theta[k - 1:], 0.0]
    return backbone(phi, psi)[-3:].ravel()

for kk in (4, 7):
    theta0 = np.r_[np.full(kk - 1, -60.0), np.full(kk - 1, -45.0)]
    J = np.array([(end_atoms(theta0 + 1e-4 * e, kk) - end_atoms(theta0 - 1e-4 * e, kk)) / 2e-4
                  for e in np.eye(len(theta0))]).T
    print(f"{len(theta0)} dihedrals: Jacobian {J.shape}, rank {np.linalg.matrix_rank(J, tol=1e-6)}")

Steric clashes restrict (\phi, \psi) to a few regions of the Ramachandran plot (Ramachandran, Ramakrishnan and Sasisekharan, 1963): the right-handed helix near (-63^\circ, -43^\circ), the β region near (-120^\circ, 130^\circ) and a left-handed region near (60^\circ, 45^\circ).

3.2 Kinematics and loop closure

So a loop with fewer than 6 dihedrals generically cannot close, one with exactly 6 closes in isolated conformations, and one with more has continuous self-motions. Code 3.3 checks the rank 6 on a helix.

Longer loops are closed numerically: cyclic coordinate descent rotates one dihedral at a time towards closure (Canutescu and Dunbrack, 2003), and Noonan, O’Brien and Snoeyink (2005) adapt robotic inverse kinematics to protein backbones.

3.3 Inverse problems and protein loop sampling

Loop sampling is an inverse problem: the end constraints are the data, and the unknowns are the conformations that satisfy them. One wants a distribution over the closure set, not one solution. A standard scheme splits \theta = (\theta_F, \theta_C): draw the free dihedrals \theta_F, then solve for the six closure dihedrals \theta_C.

The matrix J_k follows from implicit differentiation of F(\theta_F, g_k(\theta_F)) = T^\star. A sampler that ignores this factor under-samples the steep regions of the closure set, where a small range of \theta_F covers a large area. Conditioning a density on closure adds a second factor (Section 8.6).

Work by O’Donnell, Agashe and Cazals (2023) characterises the geometric constraints within tripeptides that decide whether reconstructions exist. O’Donnell and Cazals (2023) use a global parameterisation of the backbone geometry to explore loop conformations more broadly.

3.4 Mixtures of densities on flat tori

import numpy as np
from scipy.special import i0e, i1e
rng = np.random.default_rng(2)

# Synthetic Ramachandran-like data on the flat torus [-pi, pi)^2: three components
mu_true = np.radians([[-63, -43], [-120, 130], [60, 45]])     # alpha, beta, left-handed
kap_true = np.array([[12, 10], [6, 4], [15, 15]])
w_true = np.array([0.55, 0.35, 0.10])
z = rng.choice(3, size=3000, p=w_true)
X = rng.vonmises(mu_true[z], kap_true[z])                      # angles in [-pi, pi)

def log_vm(x, mu, kap):                                        # log density of a von Mises on S^1
    return kap * np.cos(x - mu) - np.log(2 * np.pi * i0e(kap)) - kap

def A_inv(r):                                                  # approximate inverse of I1/I0 (Banerjee et al., 2005)
    return r * (2 - r**2) / (1 - r**2)

def em_torus(X, K, iters=200):
    mu = X[rng.choice(len(X), K, replace=False)]
    kap, w = np.full((K, 2), 1.0), np.full(K, 1 / K)
    for _ in range(iters):
        logp = np.log(w) + sum(log_vm(X[:, None, j], mu[None, :, j], kap[None, :, j]) for j in range(2))
        ll = np.logaddexp.reduce(logp, axis=1)
        R = np.exp(logp - ll[:, None])                         # E-step: responsibilities
        w = R.mean(0)                                          # M-step
        C, S = R.T @ np.cos(X), R.T @ np.sin(X)                # weighted resultant vectors, per angle
        mu = np.arctan2(S, C)
        rbar = np.sqrt(C**2 + S**2) / R.sum(0)[:, None]
        kap = A_inv(np.clip(rbar, 1e-6, 0.999))
    return w, mu, kap, ll.sum()

w, mu, kap, ll = em_torus(X, 3)
order = np.argsort(-w)
print("weights:", w[order].round(3))
print("means (deg):", np.degrees(mu[order]).round(1).tolist())
print("concentrations:", kap[order].round(1).tolist())
print("check of the approximate inverse at kappa = 10:", round(A_inv(i1e(10) / i0e(10)), 2))

The number of components is chosen by BIC or by minimum message length, which Amarasinghe and coauthors (2023) use for joint distributions of backbone and side-chain dihedrals.

4. Boltzmann samplers for sequences and structures

A family of homologous sequences is a sample from an unknown distribution. This chapter fits it with the least biased model compatible with the data, samples from it, compares it with self-attention, uses it for design, and ends with samplers for structures.

4.1 Maximum entropy and Potts models

Write a sequence of length L as a = (a_1, \dots, a_L), with each a_i in an alphabet of q = 21 symbols (20 amino acids and the gap). Let f_i and f_{ij} be the single-site and pair frequencies of the MSA (Definition 2.5).

Pseudo-likelihood (Besag, 1975) is consistent (Section 8.2), and its use for proteins is due to Ekeberg and coauthors (2013).

Unlike mutual information (Definition 2.5), couplings separate direct from indirect correlations. With a mean-field approximation J \approx -C^{-1}, Morcos and coauthors (2011) showed that top-ranked couplings capture native contacts across many families.

4.2 Sampling sequences

A Gibbs sampler is the natural Boltzmann sampler for a Potts model. It picks a site and redraws it from the conditional of Proposition 4.3. This is Metropolis–Hastings with acceptance probability 1, so it leaves P invariant (Algorithm 0.8). Its mixing slows down when strong couplings create several modes.

import numpy as np
from sklearn.linear_model import LogisticRegression
rng = np.random.default_rng(3)
L, q, M = 12, 4, 4000

# A synthetic "family": a Potts model with 8 coupled pairs (the "contacts")
pairs = [(0, 5), (1, 9), (2, 3), (3, 10), (4, 7), (6, 11), (8, 9), (0, 11)]
J = np.zeros((L, L, q, q))
for i, j in pairs:
    B = rng.normal(0, 1.0, (q, q)); J[i, j], J[j, i] = B, B.T
h = rng.normal(0, 0.5, (L, q))

# Boltzmann sampler: Gibbs sampling, M chains in parallel, 300 sweeps
S = rng.integers(q, size=(M, L))
for sweep in range(300):
    for i in range(L):
        field = h[i] + sum(J[i, j][:, S[:, j]].T for j in range(L) if j != i)   # (M, q)
        p = np.exp(field - field.max(1, keepdims=True)); p /= p.sum(1, keepdims=True)
        S[:, i] = (p.cumsum(1) > rng.random((M, 1))).argmax(1)

# Pseudo-likelihood = one multinomial logistic regression per site (Proposition 4.3)
onehot = np.eye(q)[S].reshape(M, L * q)
Jhat = np.zeros((L, L, q, q))
for i in range(L):
    cols = [c for c in range(L * q) if c // q != i]
    clf = LogisticRegression(C=0.05, max_iter=2000).fit(onehot[:, cols], S[:, i])
    W = clf.coef_.reshape(q, L - 1, q)                     # W[a, j', b]
    for jj, j in enumerate([j for j in range(L) if j != i]):
        Jhat[i, j] = W[:, jj, :]

Jz = Jhat - Jhat.mean(2, keepdims=True) - Jhat.mean(3, keepdims=True) + Jhat.mean((2, 3), keepdims=True)
F = np.linalg.norm(Jz, axis=(2, 3)); F = (F + F.T) / 2      # zero-sum gauge, symmetrised score
iu = np.triu_indices(L, 1)
top = np.argsort(-F[iu])[:8]
found = {(int(iu[0][t]), int(iu[1][t])) for t in top}
is_contact = np.array([(i, j) in pairs for i, j in zip(*iu)])
print("precision of the top 8 pairs:", len(found & set(pairs)) / 8)
print("mean score, contacts vs non-contacts:", F[iu][is_contact].mean().round(3), F[iu][~is_contact].mean().round(3))

4.3 Potts models versus transformers

The couplings are “generalised”: nothing forces J_{ij}(a, b) = J_{ji}(b, a). Rende, Gerace, Laio and Goldt (2024) develop this mapping, while standard attention, whose weights depend on the symbols, represents higher-order interactions. Contacts can be read from the attention maps of protein language models (Rives and coauthors, 2021; Rao and coauthors, 2021).

4.4 Application in protein design, and samplers for structures

Russ and coauthors (2020) designed chorismate mutase enzymes by sampling a Potts model of the natural family and testing the sequences in bacteria lacking the enzyme. Many low-energy designs were functional, including sequences far from natural ones, whereas samples of a model without couplings were not: the couplings carry the information needed for function.

The same weights appear in Algorithm 3.8: a learned or geometric proposal, followed by a correction, is the common form of Boltzmann samplers for structures.

5. Computer practical: molecular manipulations

The practical connects the previous chapters to the tools of structural biology. Viewers and AlphaFold2 need installed software or network access, so this page describes them; the computations run in NumPy on synthetic backbones.

5.1 Molecular visualisation with PyMOL and VMD

A PDB file is a text file with one fixed-width ATOM line per atom: atom name, residue name, chain, residue number, coordinates in Å, occupancy and a B-factor column. AlphaFold2 writes its pLDDT in the B-factor column, so any viewer can colour a prediction by confidence.

import numpy as np
k = np.arange(20)
xyz = np.c_[2.3 * np.cos(np.radians(100 * k)), 2.3 * np.sin(np.radians(100 * k)), 1.5 * k]
plddt = 90 - 40 * np.exp(-((k - 10) / 3.0) ** 2)            # mock confidence, low in the middle
lines = [f"ATOM  {i + 1:5d}  CA  ALA A{i + 1:4d}    {x:8.3f}{y:8.3f}{z:8.3f}{1.0:6.2f}{b:6.2f}           C"
         for i, ((x, y, z), b) in enumerate(zip(xyz, plddt))]
pdb = "\n".join(lines + ["END"])
print("\n".join(pdb.splitlines()[9:12]))
print("columns of the first line:", len(lines[0]))                # 78 = up to the element symbol
# open("helix.pdb", "w").write(pdb)   # then, in PyMOL: load helix.pdb; spectrum b

The practical uses standard commands.

• PyMOL: fetch 1ubq downloads a PDB entry; show cartoon; spectrum b colours by B-factor; align mobile, target aligns sequences, then superposes with outlier rejection; super and cealign superpose without sequence identity; rms_cur gives the RMSD without fitting.
• VMD: in the Tcl console, measure fit returns the optimal transformation between two selections, move applies it, and measure rmsd gives the RMSD.

5.2 Structure predictions with AlphaFold2

ColabFold (Mirdita and coauthors, 2022) runs AlphaFold2 in a notebook and builds the MSA with the fast MMseqs2 search. The practical follows four steps.

1. Submit the sequence of a small protein whose structure is in the PDB.
2. Read the ranked models: the pLDDT per residue and the PAE matrix (Section 2.2).
3. Superpose the best model on the experimental structure (Theorem 1.3) and compute the RMSD, the GDT_TS and the TM-score (Code 2.12).
4. Check whether the residues of low pLDDT are the ones with large errors or are disordered in the experiment.

5.3 Molecular distances and structural alignments

So distance-based scores (dRMSD, contact maps, lDDT) cannot tell a structure from its mirror image, while the RMSD over SO(3) can.

5.4 Conformational loop sampling

import numpy as np
from scipy.optimize import least_squares
rng = np.random.default_rng(8)

def place(a, b, c, bond, angle, torsion):                    # NeRF, as in Code 3.3
    bc = (c - b) / np.linalg.norm(c - b)
    n = np.cross(b - a, bc); n /= np.linalg.norm(n)
    th, to = np.radians(angle), np.radians(torsion)
    d = bond * np.array([-np.cos(th), np.sin(th) * np.cos(to), np.sin(th) * np.sin(to)])
    return c + d[0] * bc + d[1] * np.cross(n, bc) + d[2] * n

def extend(X, phi, psi):
    """Append residues to the atoms X = [..., N, CA, C] with dihedrals (psi of the previous residue, phi)."""
    X = list(X)
    for ps, ph in zip(psi, phi):
        X.append(place(X[-3], X[-2], X[-1], 1.329, 116.2, ps))   # N   (psi of previous residue)
        X.append(place(X[-3], X[-2], X[-1], 1.458, 121.7, 180))  # CA  (omega = 180)
        X.append(place(X[-3], X[-2], X[-1], 1.525, 111.2, ph))   # C   (phi)
    return np.array(X)

# A 12-residue native: helix, a loop (residues 4-7), helix
phi = np.array([-60, -60, -60, -60, -80, 60, -90, -120, -60, -60, -60, -60.0])
psi = np.array([-45, -45, -45, 140, 160, 30, 0, 130, -45, -45, -45, -45.0])
t = np.radians(180 - 111.2)
start = np.array([[0, 0, 0], [1.458, 0, 0], [1.458 + 1.525 * np.cos(t), 1.525 * np.sin(t), 0]])
native = extend(start, phi[1:], psi[:-1])                     # atoms N0, CA0, C0, N1, ...

anchor, target = native[9:12], native[23:26]                 # N3, CA3, C3 fixed; C7, N8, CA8 must be reached

def loop(theta):                                             # theta = (psi3..psi7, phi4..phi7)
    return extend(anchor, np.r_[theta[5:], 0.0], theta[:5])[3:-1]   # N4 ... C7, N8, CA8

def closure(theta):
    return (loop(theta)[-3:] - target).ravel()

# Loop sampling: random starts from a Ramachandran-like prior, then closure by least squares
loops = []
for trial in range(40):
    a = rng.random(9) < 0.5                                  # alpha or beta basin, per dihedral
    theta0 = np.r_[np.where(a[:5], rng.normal(-43, 15, 5), rng.normal(130, 20, 5)),
                   np.where(a[5:], rng.normal(-63, 15, 4), rng.normal(-120, 20, 4))]
    sol = least_squares(closure, theta0)
    if np.abs(sol.fun).max() < 1e-3:
        loops.append(loop(sol.x)[:-2])                       # loop atoms N4 ... C7
loops = np.array(loops)
dev = np.sqrt(((loops - native[12:24]) ** 2).sum(-1).mean(-1))   # RMSD in place: anchors are fixed
print(f"closed loops: {len(loops)} / 40; RMSD to native: min {dev.min():.2f}, median {np.median(dev):.2f} A")

# Distances: the C-alpha distance matrix, the contact map and the superposition-free dRMSD
CA = native[1::3]
D = np.linalg.norm(CA[:, None] - CA[None], axis=-1)
sep = np.abs(np.subtract.outer(np.arange(12), np.arange(12)))
print("contacts (CA-CA < 8 A, |i - j| >= 3):", int(((D < 8) & (sep >= 3)).sum() // 2))
model = native.copy(); model[12:24] = loops[np.argmin(dev)]
Dm = np.linalg.norm(model[1::3][:, None] - model[1::3][None], axis=-1)
iu = np.triu_indices(12, 1)
print(f"dRMSD of the closest sampled loop: {np.sqrt(((D - Dm)[iu] ** 2).mean()):.2f} A")

The least-squares solver does not sample the closure set with a known density. Weighting the solutions as in Proposition 3.7 and Algorithm 3.8 is the subject of Exercise 17.

6. High-dimensional sampling: volumes of polytopes and densities of states

How large is a region of a high-dimensional space, and how can it be sampled? On polytopes, this chapter shows that exact computation is hopeless, builds the random walks and estimators that work, and carries them over to densities of states.

6.1 Polytopes and their volumes: hardness and approximability

6.2 Linear programs and hit-and-run

Berbee and coauthors (1987) use hit-and-run to find non-redundant constraints: almost surely the line leaves K through the interior of a facet, so the constraint attaining t_+ is not redundant.

6.3 Estimating volumes with hit-and-run and PDMPs

import numpy as np
from math import factorial, gamma, pi
rng = np.random.default_rng(4)

n = 5                                            # K = {x : x_i >= -1, sum x_i <= 1}, a simplex
A = np.vstack([-np.eye(n), np.ones((1, n))])     # A x <= b
b = np.r_[np.ones(n), 1.0]
exact = (n + 1) ** n / factorial(n)              # translate y = x + 1: y >= 0, sum y <= n + 1

def har_chain(x, r, steps):
    """Hit-and-run in K ∩ B(0, r): uniform direction, then uniform point on the chord."""
    out = []
    for _ in range(steps):
        u = rng.standard_normal(n); u /= np.linalg.norm(u)
        Au, slack = A @ u, b - A @ x                       # chord in K: A(x + t u) <= b
        tmax = np.min(slack[Au > 0] / Au[Au > 0]); tmin = np.max(slack[Au < 0] / Au[Au < 0])
        xu, xx = x @ u, x @ x                              # chord in the ball: |x + t u| <= r
        disc = np.sqrt(xu**2 - xx + r**2)
        lo, hi = max(tmin, -xu - disc), min(tmax, -xu + disc)
        x = x + rng.uniform(lo, hi) * u
        out.append(x)
    return x, np.array(out)

r0 = 1 / np.sqrt(n)                              # B(0, r0) ⊂ K (distance to the facet sum x = 1)
R = np.sqrt(n**2 + n - 1)                        # K ⊂ B(0, R) (farthest vertex)
m = int(np.ceil(n * np.log2(R / r0)))            # radii r_i = r0 * 2^(i/n), i = 0..m
radii = np.minimum(r0 * 2.0 ** (np.arange(m + 1) / n), R)
vol = pi ** (n / 2) / gamma(n / 2 + 1) * r0 ** n  # volume of the inner ball
x = np.zeros(n)
for i in range(1, m + 1):
    x, pts = har_chain(x, radii[i], 4000)
    pts = pts[::4]                                # thinning
    ratio = np.mean(np.linalg.norm(pts, axis=1) <= radii[i - 1])   # vol(K_{i-1}) / vol(K_i)
    vol /= ratio
print(f"n = {n}, phases m = {m}, estimate {vol:.1f}, exact {exact:.1f}")

Chevallier, Cazals and Fearnhead (2022) use such PDMPs instead of hit-and-run inside the multiphase scheme for polytopes in high dimension. Cousins and Vempala (2016) made the scheme practical with Gaussian cooling, which integrates Gaussians of increasing variance instead of indicators of balls.

6.4 Densities of states: the Wang–Landau algorithm

import numpy as np
rng = np.random.default_rng(5)
Lside = 4; N = Lside * Lside
nbr = np.array([[((r + dr) % Lside) * Lside + (c + dc) % Lside for dr, dc in ((0, 1), (1, 0), (0, -1), (-1, 0))]
                for r in range(Lside) for c in range(Lside)])

# Exact density of states by enumeration of the 2^16 configurations (periodic 4 x 4 Ising)
conf = 1 - 2 * ((np.arange(2**N)[:, None] >> np.arange(N)) & 1)
E_all = -(conf * conf[:, nbr[:, 0]]).sum(1) - (conf * conf[:, nbr[:, 1]]).sum(1)
levels, g_exact = np.unique(E_all, return_counts=True)
idx = {e: k for k, e in enumerate(levels)}

# Wang-Landau with the 1/t refinement of the modification factor
s = rng.choice([-1, 1], N); E = -(s * s[nbr[:, 0]]).sum() - (s * s[nbr[:, 1]]).sum()
lng, H, lnf = np.zeros(len(levels)), np.zeros(len(levels)), 1.0
for step in range(1, 1_500_001):
    i = rng.integers(N)
    E_new = E + 2 * s[i] * s[nbr[i]].sum()
    if np.log(rng.random()) < lng[idx[E]] - lng[idx[E_new]]:   # accept with min(1, g(E)/g(E'))
        s[i] = -s[i]; E = E_new
    k = idx[E]; lng[k] += lnf; H[k] += 1
    if step % 1000 == 0:
        if lnf > len(levels) / step and H.min() > 0.8 * H.mean():
            lnf, H = lnf / 2, np.zeros(len(levels))       # flat histogram: halve ln f
        if lnf <= len(levels) / step:
            lnf = len(levels) / step                        # 1/t regime
lng += np.log(2**N) - np.logaddexp.reduce(lng)              # normalise: sum_E g(E) = 2^N
err = np.abs(lng - np.log(g_exact))
print("energy levels:", levels.tolist())
print("exact g(E):", g_exact.tolist())
print(f"max |ln g_WL - ln g_exact| = {err.max():.3f}, final ln f = {lnf:.1e}")

Chevallier and Cazals (2020) design an adapted random walk to speed up the convergence of Wang–Landau in continuous spaces, where the levels are bins of energy.

7. Spatial partitions and applications

Sampled conformations are clouds of points in high dimension. This chapter partitions them with random hyperplanes, which adapt to the intrinsic dimension, explains predictions with Shapley values, and uses both tools to check whether samplers have mixed.

7.1 Random projection trees

With a = q - x_{(1)} and b = q - x_{(i)}, the lemma says that a random split rarely separates a query q from its nearest neighbour x_{(1)} while keeping a farther point x_{(i)} with it.

7.2 Nearest neighbours, regression and dimension estimation

Several independent trees reduce the failure probability, with the union of the leaves of q as candidates. The same partitions give a regression estimator, the mean response in the leaf of x, whose convergence rate depends on the intrinsic dimension rather than on D (Kpotufe and Dasgupta, 2012).

import numpy as np
from sklearn.datasets import load_digits
rng = np.random.default_rng(6)

def rp_tree(idx, X, leaf=20):
    """Random projection tree: split at the median of a projection on a random direction."""
    if len(idx) <= leaf:
        return idx
    u = rng.standard_normal(X.shape[1]); u /= np.linalg.norm(u)
    proj = X[idx] @ u
    t = np.median(proj)
    return (u, t, rp_tree(idx[proj <= t], X, leaf), rp_tree(idx[proj > t], X, leaf))

def leaf_of(node, x):
    while isinstance(node, tuple):
        u, t, left, right = node
        node = left if x @ u <= t else right
    return node

X = load_digits().data.astype(float)                 # 1797 images of 8 x 8 pixels, D = 64
queries, base = X[:200], X[200:]
true_nn = np.argmin(((queries[:, None] - base[None]) ** 2).sum(-1), axis=1)
for n_trees in (1, 5, 20):
    trees = [rp_tree(np.arange(len(base)), base) for _ in range(n_trees)]
    hits = 0
    for qi, x in enumerate(queries):
        cand = np.unique(np.concatenate([leaf_of(T, x) for T in trees]))
        hits += cand[np.argmin(((base[cand] - x) ** 2).sum(1))] == true_nn[qi]
    print(f"{n_trees:2d} trees: exact nearest neighbour found for {hits / len(queries):.0%} of queries")

# Cell radii shrink at a rate set by the intrinsic dimension, not by D = 100
def diam_profile(X, depth=8):
    cells, out = [np.arange(len(X))], []
    for _ in range(depth):
        new = []
        for c in cells:
            u = rng.standard_normal(X.shape[1]); proj = X[c] @ u; t = np.median(proj)
            new += [c[proj <= t], c[proj > t]]
        cells = new
        out.append(np.mean([np.max(np.linalg.norm(X[c] - X[c].mean(0), axis=1)) for c in cells]))
    return np.array(out) / out[0]
for d in (2, 20):
    Z = rng.uniform(-1, 1, (4096, d)) @ rng.standard_normal((d, 100)) / np.sqrt(d)
    prof = diam_profile(Z)
    slope = np.polyfit(np.arange(1, 9), np.log2(prof), 1)[0]          # log2 radius lost per level
    print(f"intrinsic dim {d:2d}: radius ratio after 8 levels {prof[-1]:.2f}, levels per halving {-1 / slope:.1f}")

7.3 Shapley values and Shapley–Shubik indices

If f has 0/1 output and the background is one point b with f(b) = 0, then v_x takes values in \{0, 1\}; when it is monotone, it is a simple game, and the Shapley values are Shapley–Shubik indices of the features.

A pass over the leaves costs O(\text{leaves} \times \text{depth}) per background point, instead of 2^d coalitions, and a forest averages its trees. Lundberg and coauthors (2020) develop exact algorithms for tree ensembles. Section 8.7 gives the proof.

7.4 Feature selection and mixing assessment

Features with a small mean |\phi_j| can be removed: this is Shapley-based feature selection. For mixing, train a classifier to tell apart the samples of two chains. A cross-validated accuracy near 0.5 means that they are indistinguishable; otherwise Shapley values point to the coordinates that are not mixed.

import numpy as np
from itertools import permutations, combinations
from math import factorial
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
rng = np.random.default_rng(7)

# 1. Shapley-Shubik index of the weighted vote [quota 51; weights 50, 49, 1]
w, quota = [50, 49, 1], 51
piv = np.zeros(3)
for order in permutations(range(3)):
    tot = 0
    for p in order:
        tot += w[p]
        if tot >= quota:
            piv[p] += 1; break
print("Shapley-Shubik indices:", (piv / factorial(3)).round(3).tolist())

# 2. Mixing assessment: can a classifier tell two MCMC chains apart, and on which coordinates?
def chain(mode, n=2000):
    X = rng.normal(0, 1, (n, 4))            # coordinates 1-3 well mixed in both chains
    X[:, 0] = rng.normal(mode, 0.4, n)       # coordinate 0: each chain stuck in one well (-1 or +1)
    return X
X = np.vstack([chain(-1.0), chain(+1.0)]); y = np.r_[np.zeros(2000), np.ones(2000)]
clf = RandomForestClassifier(n_estimators=50, max_depth=5, random_state=0)
print("cross-validated accuracy:", cross_val_score(clf, X, y, cv=5).mean().round(3), "(0.5 = indistinguishable)")
clf.fit(X, y)

def shapley(f, x, background):
    """Exact interventional Shapley values: v(S) = mean_b f(x_S, b_rest), 2^d coalitions."""
    d = len(x); phi = np.zeros(d)
    def v(S):
        Z = background.copy(); Z[:, list(S)] = x[list(S)]
        return f(Z).mean()
    for j in range(d):
        others = [k for k in range(d) if k != j]
        for size in range(d):
            for S in combinations(others, size):
                wS = factorial(size) * factorial(d - size - 1) / factorial(d)
                phi[j] += wS * (v(S + (j,)) - v(S))
    return phi
f = lambda Z: clf.predict_proba(Z)[:, 1]
bg = X[rng.choice(len(X), 100, replace=False)]
rows = rng.choice(len(X), 30, replace=False)
Phi = np.array([shapley(f, X[i], bg) for i in rows])
print("mean |Shapley value| per coordinate:", np.abs(Phi).mean(0).round(3).tolist())
print("efficiency, sum_j phi_j = f(x) - mean_b f(b):", np.allclose(Phi.sum(1), f(X[rows]) - f(bg).mean()))

# 3. Proposition 7.12: exact Shapley values of a tree for one background point, in one pass over the leaves
def tree_shap(est, x, b):
    t, phi = est.tree_, np.zeros(len(x))
    x, b = x.astype(np.float32), b.astype(np.float32)          # scikit-learn trees compare in float32
    def walk(node, A, B):
        if t.children_left[node] < 0:                          # leaf of value f_l = P(class 1)
            f_l, s, r = t.value[node, 0, 1] / t.value[node, 0].sum(), len(A), len(B)
            for j in A: phi[j] += f_l * factorial(s - 1) * factorial(r) / factorial(s + r)
            for j in B: phi[j] -= f_l * factorial(s) * factorial(r - 1) / factorial(s + r)
            return
        j, thr = t.feature[node], t.threshold[node]
        for child, left in ((t.children_left[node], True), (t.children_right[node], False)):
            px, pb = (x[j] <= thr) == left, (b[j] <= thr) == left   # does each point pass this test?
            if px and pb: walk(child, A, B)
            elif px and j not in B: walk(child, A | {j}, B)    # passed by x only: j must be in S
            elif pb and j not in A: walk(child, A, B | {j})    # passed by b only: j must not be in S
    walk(0, frozenset(), frozenset())
    return phi
Phi_tree = np.array([np.mean([tree_shap(est, X[i], bb) for est in clf.estimators_ for bb in bg], axis=0)
                     for i in rows[:3]])                       # a forest averages its trees
print("leaf formula of Proposition 7.12 = enumeration:", np.allclose(Phi_tree, Phi[:3]))

8. Complements for self-study

This chapter collects the proofs and side results that the sessions only state. Each section can be read right after the session that uses it.

8.1 The Eyring–Kramers formula in one dimension

Proof of Theorem 1.11 in one dimension. Let m < s < b, with m the lowest point left of the saddle s and b beyond s.

As \beta \to \infty, the outer integrand concentrates near y = s. There the inner integral has already collected the whole well around m. Laplace’s method gives two Gaussian integrals: \int_{-\infty}^{s} e^{-\beta U} \simeq e^{-\beta U(m)} \sqrt{\frac{2\pi}{\beta\, U''(m)}}, \qquad \int e^{\beta U(y)} dy \simeq e^{\beta U(s)} \sqrt{\frac{2\pi}{\beta\, |U''(s)|}} .

Their product, times \beta, is the formula of Theorem 1.11; the factors \beta cancel.

8.2 Pseudo-likelihood is consistent

If the conditionals determine \theta in the chosen gauge, the minimiser is unique. The empirical loss is convex (Proposition 4.3) and converges pointwise, so its minimiser converges to \theta^\star, by the standard argument for convex M-estimators.

8.3 Invariance of the billiard PDMP

The generator of the process of Definition 6.10 acts on functions f(x, v) as \mathcal{L} f = v \cdot \nabla_x f + \lambda \Big( \int f(x, w)\, d\sigma(w) - f(x, v) \Big), on functions that satisfy f(x, v) = f(x, R_\nu v) at the boundary, with R_\nu the reflection.

Let \mu be uniform on K \times S^{n-1}. The refreshment term integrates to zero against \mu, because \sigma is the refreshment law. By the divergence theorem, the transport term gives \int_{S^{n-1}} \int_{\partial K} (v \cdot \nu)\, f(x, v)\, dA(x)\, d\sigma(v).

The reflection maps outgoing velocities (v \cdot \nu > 0) to incoming ones, preserves \sigma and flips the sign of v \cdot \nu. With the boundary condition, the two halves cancel, so \int \mathcal{L} f\, d\mu = 0 and \mu is invariant.

8.4 Wang–Landau as stochastic approximation

Write \theta_t(E) = \log \tilde g_t(E). The update of Algorithm 6.12 with step \gamma_t reads \theta_{t+1}(E) = \theta_t(E) + \gamma_{t+1}\, \mathbf{1}\big[U(X_{t+1}) = E\big].

Its mean field is the vector of level probabilities under the biased target \pi_\theta(x) \propto e^{-\theta(U(x))}. With \varphi = \theta - \log g, these probabilities are p_\theta(E) = \frac{e^{-\varphi(E)}}{\sum_{E'} e^{-\varphi(E')}} .

The mean ODE \dot\theta = p_\theta raises \varphi fastest where it is smallest, so \max \varphi - \min \varphi decreases. Its limit points are the \theta = \log g + c, which give the flat histogram of Proposition 6.13. Fort and coauthors (2015) turn this picture into an almost-sure convergence theorem for decreasing steps with \sum_t \gamma_t = \infty and \sum_t \gamma_t^2 < \infty.

8.5 Tripeptide loop closure: the elimination

The rigid peptide units fix the distances between consecutive Cα atoms, and the anchors fix the first and last Cα. The three Cα of the tripeptide therefore form a triangle of known shape. The two peptide units can each turn only about their side of the triangle, and the triangle can turn about its fixed side: three rotation angles \tau_1, \tau_2, \tau_3 replace the six dihedrals.

At each Cα, the fixed bond angle N–Cα–C couples the angles of the two sides that meet there. Each such condition is an equation of degree 2 in each of the two half-angle tangents u_i = \tan(\tau_i/2) involved. Eliminating u_1 and u_3 from the three equations with resultants leaves a univariate polynomial of degree 16 in u_2 (Coutsias and coauthors, 2004). Its real roots give the closures; each root determines the dihedrals back.

8.6 Conditioning on closure: the co-area formula

The result depends on the coordinates chosen on SE(3) to express closure, an instance of the Borel–Kolmogorov paradox: conditioning on an event of probability zero needs a choice of limit. Combined with Proposition 3.7, the importance weight of a closed loop is \pi / (J_F\, q) for one branch chosen uniformly, and \pi / (J_F\, K q) when all K closures are kept (Algorithm 3.8).

8.7 Shapley values: uniqueness and decision trees

Proof of Theorem 7.8. The formula of Definition 7.7 satisfies the four axioms by inspection.

Conversely, the unanimity games u_T(S) = \mathbf{1}[T \subseteq S], for non-empty T, form a basis of the space of games.

In c\, u_T, players outside T are null and players in T are symmetric, so efficiency forces \phi_j(c\, u_T) = c / |T| for j \in T and 0 otherwise. Additivity extends the values to every game, which proves uniqueness. ∎

Proof of Proposition 7.12.

• Write z_S for the point equal to x on S and to b elsewhere, so that v_x(S) = f(z_S) - f(b).
• At a test on feature j, the point z_S behaves like x if j \in S and like b otherwise.
• So a test passed by both points is passed by z_S, and a test failed by both is failed. A test passed by x only requires j \in S; one passed by b only requires j \notin S.
• Hence a leaf with a test failed by both points, or with A_\ell \cap B_\ell \neq \emptyset, is never reached. Any other leaf is reached exactly when A_\ell \subseteq S \quad \text{and} \quad B_\ell \cap S = \emptyset .
• Let u_\ell(S) be 1 when z_S reaches \ell, and 0 otherwise. Then v_x(S) = \sum_\ell f_\ell\, u_\ell(S) - f(b). The formula of Definition 7.7 is linear in v and ignores additive constants, so \phi_j(v_x) = \sum_\ell f_\ell\, \phi_j(u_\ell).
• Features outside A_\ell \cup B_\ell are null for u_\ell. So \phi_j(u_\ell) is the average marginal contribution of j over uniform random orders of the s + r features of A_\ell \cup B_\ell.
• For j \in A_\ell, the contribution is 1 exactly when j comes after the other features of A_\ell and before those of B_\ell. Among the (s + r)! orders, this happens in (s-1)!\, r! of them.
• For j \in B_\ell, the contribution is -1 exactly when j comes after all of A_\ell and before the rest of B_\ell: this happens in s!\, (r-1)! orders.

Dividing both counts by (s + r)! gives the formulas of Proposition 7.12. ∎

8.8 Maximum likelihood for the von Mises density

Proof of Proposition 3.10. The weighted log-likelihood is \kappa \sum_k r_k \cos(x_k - \mu) - \Big(\sum_k r_k\Big) \log I_0(\kappa) + \text{const}.

The sum equals \sqrt{C^2 + S^2}\, \cos(\mu - \hat\mu), maximal at \mu = \hat\mu. The derivative in \kappa then vanishes when A(\kappa) = \bar R, since I_0' = I_1.

Now \log I_0(\kappa) is the log-partition function of an exponential family with statistic \cos(x - \mu), so its second derivative A'(\kappa) is a variance, hence positive. A increases strictly from 0 to 1, and the root is unique. ∎

9. Research directions and open problems

This chapter surveys active questions related to each part of the course. It is meant for self-study and for choosing a project; the cited works are starting points, not a complete bibliography.

From structures to conformational ensembles

Structure predictors return one model, while function needs populations of states (ch. 1). Del Alamo and coauthors (2022) obtained alternative conformations of transporters and receptors from AlphaFold2 by subsampling the MSA. Wayment-Steele and coauthors (2024) cluster the MSA to predict multiple conformations of fold-switching proteins. AlphaFold3 (Abramson and coauthors, 2024) extends prediction to complexes with nucleic acids, ligands and ions.

Open problem: whether such predictions carry thermodynamic information, that is, populations and free-energy differences, rather than a list of plausible states.

Long time scales and rare events

Special-purpose hardware made millisecond simulations of small proteins possible (Shaw and coauthors, 2010). Markov state models summarise many short trajectories into a kinetic network (Prinz and coauthors, 2011). Metadynamics biases the dynamics along collective variables to escape free-energy minima (Laio and Parrinello, 2002).

Open problem: finding good reaction coordinates, the \xi of Definition 1.7, in dimension 15,000, possibly by learning them.

Generative models for structures and design

Boltzmann generators (Noé and coauthors, 2019) learn proposals for equilibrium sampling. Diffusion models generate new backbones (Watson and coauthors, 2023), and inverse-folding networks design sequences for them (Dauparas and coauthors, 2022).

Open problems: unbiased Boltzmann sampling with learned proposals in large systems, where importance weights degenerate; and the experimental success rate of designs.

Language models and Potts models

Protein language models predict structure without MSAs (Lin and coauthors, 2023). The mapping between attention and Potts models (Rende and coauthors, 2024) suggests asking which higher-order interactions deeper models learn, and whether they improve design over pairwise models (Russ and coauthors, 2020).

Geometry and topology of conformation spaces

Closure sets of loops are real algebraic varieties (ch. 3), whose dimension, singularities and connected components govern sampling. Persistent homology has been used to describe protein structure, flexibility and folding (Xia and Wei, 2014). The study of tripeptide reconstructions by O’Donnell, Agashe and Cazals (2023) shows how geometric constraints decide whether closures exist.

High-dimensional volumes and densities of states

Practical volume algorithms (Cousins and Vempala, 2016; Chevallier, Cazals and Fearnhead, 2022) handle polytopes far beyond the reach of exact methods, and pushing them to larger dimensions and to non-convex bodies is active work. Convergence rates of Wang–Landau and of its variants remain partly open (Fort and coauthors, 2015). Energy landscapes are not convex, and sublevel sets of the energy can have many components, so the guarantees of ch. 6 do not apply.

Explaining and validating samplers

Classifier-based two-sample tests and Shapley values (ch. 7) give diagnostics for mixing. Their statistical calibration on correlated MCMC output, and the choice of the background distribution for Shapley values, are open questions.

Exercises

The exercises are sorted by level: ★ applies the course, ★★ combines results or proves a short statement, ★★★ is open-ended. The session plan says which ones to do after each session. Exercises 13, 15, 17, 18, 19 and 20 are computational.

1. ★ Counting degrees of freedom. (a) What is the dimension of the space of shapes of a molecule of 5000 atoms? (b) A chain of 300 residues has ideal bond lengths and angles and trans peptides. How many backbone dihedrals move its backbone atoms, and what fraction is this of the 3 \times 900 - 6 coordinates of its backbone shape?

   Solution: (a) 3 \times 5000 - 6 = 14{,}994 (Definition 1.1).

   2. Each residue has \phi and \psi, but \phi_1 and \psi_{300} move no backbone atom of the chain, which leaves 2 \times 300 - 2 = 598 dihedrals. The fraction is 598 / 2694 \approx 0.22: fixing bonds and angles removes almost four fifths of the coordinates.

2. ★ Energy against entropy. State B lies \Delta U = 2 kcal/mol above state A, with an entropy larger by \Delta S = 0.01 kcal/(mol·K). Compute \mathbb{P}(B)/\mathbb{P}(A) at 300 K, then the temperature at which both states are equally populated.

   Solution: by Proposition 1.5, \Delta F = \Delta U - T \Delta S = 2 - 3 = -1 kcal/mol at 300 K, so the ratio is e^{1/0.596} \approx 5.35. The populations are equal when \Delta F = 0, that is, at T = \Delta U / \Delta S = 200 K; below it A dominates.

3. ★ Needleman–Wunsch by hand. Align HEAG with HAG (match 2, mismatch -1, gap -2). Fill the table of Proposition 2.2 and trace back an optimal alignment.

   Solution: with rows H, E, A, G and columns H, A, G, after the borders 0, -2, -4, \dots, the rows of the table are (2, 0, -2), (0, 1, -1), (-2, 2, 0) and (-4, 0, 4). The score is 4. The traceback gives HEAG over H-AG: a match, a gap opposite E, then two matches.

4. ★ The TM-score scale. (a) Compute d_0 for L = 100 and L = 300. (b) For L = 100, a model has 80 % of its residues at 1 Å and 20 % at 10 Å after superposition. Compute its TM-score for this superposition.

   Solution: (a) d_0 = 1.24 \times 85^{1/3} - 1.8 = 3.65 Å and 6.36 Å; a longer chain tolerates larger deviations. (b) The score is 0.8 / (1 + (1/3.65)^2) + 0.2 / (1 + (10/3.65)^2) = 0.744 + 0.024 = 0.768 . The badly placed fifth contributes almost nothing, and does not pull the score down as it would pull the RMSD up.

5. ★ A Shapley–Shubik index. Compute the indices of the weighted vote with quota 4 and weights 3, 2, 1 and 1.

   Solution: the winning coalitions are those that contain player 1 and someone else, and \{2, 3, 4\}. Player 1 is pivotal exactly when it arrives second or third, in 12 of the 24 orders, so its index is 1/2. Players 2, 3 and 4 play symmetric roles in this game, so each gets 1/6: weight 2 brings no more power than weight 1.

6. ★ Rejection in high dimension. How many uniform proposals in [-1, 1]^n are needed, on average, to obtain 100 points of the unit ball, for n = 10 and n = 20?

   Solution: by Proposition 6.3 the acceptance rates are 2.49 \times 10^{-3} and 2.46 \times 10^{-8}, so about 4.0 \times 10^4 and 4.1 \times 10^9 proposals.

7. ★★ Superposition in the plane. For centred points p'_i, q'_i \in \mathbb{R}^2, show that the optimal rotation angle is \theta^\star = \operatorname{atan2}\Big( \sum_i p'_i \times q'_i,\ \sum_i p'_i \cdot q'_i \Big), with p \times q = p_x q_y - p_y q_x. Why can no reflection problem arise in this formula?

   Solution: with A the sum of dot products and B the sum of cross products, the cost is a constant minus 2 \sum_i q'_i \cdot R_\theta p'_i = 2 (A \cos\theta + B \sin\theta).

   The maximum is at \theta = \operatorname{atan2}(B, A). The search runs over rotations only, so the determinant is 1 by construction; Theorem 1.3 needs the sign s because the SVD searches over O(3).

8. ★★ Breaking the invariance. In Definition 2.8, replace the distance term by an inner product in the global frame: -\gamma \sum_p \big\langle T_i q_{ip},\ T_j k_{jp} \big\rangle . Is the attention still invariant under global rotations? Under global translations?

   Solution: rotations preserve inner products, so rotations about the origin leave the logits unchanged. A translation by t adds \langle a_i + b_j, t \rangle + P \|t\|^2 to the sum, where a_i and b_j are the sums over p of the transformed points and P is the number of points.

   The term P \|t\|^2 and the part \langle a_i, t \rangle do not depend on j, so they cancel in the softmax over j; the part \langle b_j, t \rangle does not cancel. The squared distance of IPA avoids this, because it depends on differences only.

9. ★★ Closing a planar arm. A planar arm has two links of lengths \ell_1, \ell_2 from the origin, with joint angles \theta_1, \theta_2. (a) For which end distances r can it reach a given point, and with how many solutions? (b) Solve for \ell_1 = \ell_2 = r = 1. (c) What is the dimension of the solution set for three links, and how does this match Proposition 3.5?

   Solution: (a) The law of cosines gives \cos\theta_2 = \frac{r^2 - \ell_1^2 - \ell_2^2}{2 \ell_1 \ell_2} . There are two solutions (elbow up and down) when |\ell_1 - \ell_2| < r < \ell_1 + \ell_2, one at the boundary and none outside. (b) \cos\theta_2 = -1/2, so \theta_2 = \pm 120^\circ. (c) The end position imposes 2 conditions on 3 angles, which leaves a curve: the planar analogue of n - 6 is n - 2 for a position, or n - 3 for a position and an orientation.

10. ★★ The von Mises concentration. (a) Show that A'(\kappa) = 1 - A(\kappa)/\kappa - A(\kappa)^2, using I_1' = I_0 - I_1/\kappa. (b) For \bar R = 0.9, compare the approximation of Code 3.11 with the exact root.

    Solution: (a) Differentiate A = I_1/I_0 with I_0' = I_1: A' = \frac{I_1'}{I_0} - \frac{I_1^2}{I_0^2} = 1 - \frac{A}{\kappa} - A^2 . It equals the variance of \cos(x - \mu) (Section 8.8), hence it is positive. (b) The approximation gives 0.9 \times 1.19 / 0.19 = 5.64, the exact root is 5.30: a 6 % bias, removed by one Newton step with the formula of (a).

11. ★★ Gauge and parameter count. (a) Show that replacing J_{ij}(a, b) by J_{ij}(a, b) + u(a) and h_i(a) by h_i(a) - u(a) leaves the Potts distribution unchanged. (b) Count the parameters for L = 100 and q = 21, before and after fixing the gauge.

    Solution: (a) The energy of every sequence changes by u(a_i) - u(a_i) = 0.

    2. Before fixing the gauge, the count is Lq + \binom{L}{2} q^2 = 2100 + 4950 \times 441 = 2{,}185{,}050 .

    In the zero-sum gauge, each field has q - 1 free entries and each coupling (q-1)^2: 100 \times 20 + 4950 \times 400 = 1{,}982{,}000.

12. ★★ Factored attention is economical. In Proposition 4.7, show that scaling A^h by c > 0 and W^h V^h by 1/c gives the same model. Then compare the number of coupling parameters of a Potts model and of 4 heads of factored attention, for 50 sites and 21 symbols, with one square symbol matrix per head.

    Solution: J only involves the products A^h_{ij} (W^h V^h)_{ab}, which are unchanged. The Potts model has \binom{50}{2} \times 441 = 540{,}225 coupling parameters. Four heads have 4 \times 2500 positional entries and 4 \times 441 symbol entries, 11,764 in total: the factorisation is a strong low-rank prior on the couplings.

13. ★★ Hit-and-run against a ball walk (computational). In the rectangle [0, 100] \times [0, 1], run 2 \times 10^5 steps of hit-and-run and of a random-walk Metropolis (Gaussian steps of 0.5, rejected outside), both from (50, 0.5). Compare the means and the autocorrelation times of x_1.

    Hints and expected result: estimate the autocorrelation time by batch means, with batches longer than the mixing time. Hit-and-run has an autocorrelation time of order 10^2 steps and a mean near 50. The ball walk moves diffusively, with steps of about 0.5, so crossing 100 units takes of order (100/0.5)^2 = 4 \times 10^4 steps. Its mean can still be off by ten or more after 2 \times 10^5 steps. Long chords make hit-and-run insensitive to the elongation, in line with the R^2/r^2 factor of Theorem 6.7 being a worst case.

14. ★★ Phases of the multiphase scheme. For the cube [-1, 1]^n, take r = 1 and R = \sqrt n. Compute the number of phases in dimensions 10 and 100, and the total number of points with N = 8m / \varepsilon^2 per phase for \varepsilon = 0.1. Why can the bound of Proposition 6.8 on the ratios fail for a non-convex body?

    Solution: m = \lceil (n/2) \log_2 n \rceil gives 17 and 333 phases. With N = 800\,m, the totals are 800\, m^2 = 2.3 \times 10^5 and 8.9 \times 10^7 points, each produced by a hit-and-run chain.

    The proof of Proposition 6.8 uses K \subset 2^{1/n} K, which needs convexity and 0 \in K; a thin non-convex shell can have \mathrm{vol}(K_{i-1}) much smaller than \mathrm{vol}(K_i).

15. ★★ Heat capacity from a density of states (computational). From the exact g(E) printed by Code 6.15, compute the heat capacity per spin C(T) = \mathrm{Var}(E) / (N T^2), with k_B = J = 1, on a grid of temperatures, and locate its maximum.

    Hints and expected result: use p_T(E) \propto g(E) e^{-(E - E_{\min})/T} to avoid overflow. The maximum is near T \approx 2.44, with C \approx 0.82 per spin. It lies above the infinite-lattice critical temperature 2/\log(1 + \sqrt 2) \approx 2.27 (Onsager), a finite-size effect. One density of states gives every temperature at once (Definition 6.11).

16. ★★ An exact projection probability. For orthogonal a, b \in \mathbb{R}^D and u uniform on the sphere, show that \mathbb{P}\big( |\langle u, a\rangle| \geq |\langle u, b \rangle| \big) = \frac{2}{\pi} \arctan \frac{\|a\|}{\|b\|}, and evaluate it for the ratios 0.5 and 0.1. Compare with Lemma 7.2.

    Solution: in the plane of a and b the direction angle \omega is uniform, and the event reads |\tan\omega| \geq \|b\|/\|a\| when \omega is measured from b. Its probability is 1 - \frac{2}{\pi}\arctan\frac{\|b\|}{\|a\|} = \frac{2}{\pi}\arctan\frac{\|a\|}{\|b\|} .

    The values are 0.295 and 0.063, below the bounds 0.5 and 0.1 of the lemma.

17. ★★★ Loop sampling with the Jacobian correction (computational). Start from Code 5.4. Split the nine dihedrals into three free and six closure dihedrals, compute J_k by implicit differentiation, and weight each closure by Proposition 3.7 and by the factor of Section 8.6. Compare the weighted and unweighted distributions of one dihedral of the loop.

    Hints and expected result: the closure Jacobian is the 9 \times 6 derivative of the three end atoms with respect to the closure dihedrals; use a least-squares solve for implicit differentiation. Expect weights that vary over at least one order of magnitude, and visible differences between weighted and unweighted histograms. Report the effective sample size (\sum w)^2 / \sum w^2.

18. ★★★ A mixing diagnostic (computational). Use the chains of Code 0.9, for h = 1 and h = 8, and test whether the first and second halves of each chain are distinguishable. Build a classifier test as in Code 7.13, on lagged features such as (x_t, x_{t-10}), and compare with the Gelman–Rubin statistic on four chains.

    Hints and expected result: for h = 1 the cross-validated accuracy should stay near 0.5. For h = 8 it depends on whether a transition happened in one half only. Discuss how the correlation of MCMC samples biases the cross-validated accuracy, and how thinning or block splits correct it.

19. ★★★ Billiard PDMP for volumes (computational). Implement the PDMP of Definition 6.10 in the polytope of Code 6.9, using it in place of hit-and-run, and compare the relative errors of both for n = 5, 10, 20 at equal numbers of oracle calls.

    Hints and expected result: between events, move to the first boundary hit (a ratio test), reflect, and refresh at rate \lambda of order 1/\mathrm{diam}(K). Sample the position at regular times, not at events. Chevallier, Cazals and Fearnhead (2022) report gains of PDMPs in high dimension; check whether they appear in your runs.

20. ★★★ Does factored attention learn the couplings? (computational) On the sequences of Code 4.6, train one layer of factored attention by masked prediction, in PyTorch, and compare its couplings (Proposition 4.7) with the pseudo-likelihood estimate.

    Hints and expected result: parametrise A^h by a softmax over j \neq i of free logits, and use 1 to 8 heads. With enough heads, the contact precision should approach that of pseudo-likelihood, as in Rende and coauthors (2024). With one head, the rank constraint should show in the coupling matrices.

Projects and assessment

master-mva.com announces projects for students working in tandem (15 points) and an individual quiz (5 points). A project reproduces or expands recently published results; students return a report and a git repository, notebook or code archive, and have one month. The catch-up is an oral exam on the lectures, typically one or two questions per lecture. The Université Paris-Saclay course detail describes the assessment as “Evaluation Continue non Intégrale” (partial continuous assessment).

The subjects below are proposals in that format. The group subjects are sized for a tandem, as in the official modality; the solo subjects suit a student who works alone or prepares the quiz. The written exam is a self-assessment proposal in the spirit of the quiz, not an official modality.

Solo projects

1. [Solo · ≈ 20 h] Paper study: Wang–Landau and its 1/t variant. Reproduce the convergence behaviour discussed by Belardinelli and Pereyra (2007) and analysed by Fort and coauthors (2015).
   • Work: extend Code 6.15 from the 4 \times 4 to the 5 \times 5 periodic Ising model, with exact g(E) by enumeration (2^{25} states, vectorised in chunks); compare the flat-histogram rule, the 1/t rule and deterministic steps \gamma_t = 1/t.
   • Success criterion: a log–log plot of \max_E |\log \tilde g - \log g| against the number of steps, with the saturation of the plain rule and the decay of the 1/t rule visible.
   • Deliverables: code, a 4–6 page report, a 10-minute talk. Chapters 6 and 8.
2. [Solo · ≈ 20 h] Paper study: random projection trees and intrinsic dimension. Reproduce the main claim of Dasgupta and Freund (2008): cell radii shrink at a rate set by the intrinsic dimension.
   • Work: synthetic manifolds (subspaces, spheres, swiss rolls) of dimension 1 to 10 in \mathbb{R}^{100}, with noise; RP trees with median and with max-type splits; the digits data set; the nearest-neighbour bound of Theorem 7.4 checked empirically against \Phi.
   • Success criterion: levels-per-halving estimates within 20 % of the true dimension up to d = 6, and failure rates that grow with \Phi.
   • Deliverables: code, report, talk. Chapter 7.
3. [Solo · ≈ 20 h] Mixtures on the torus with model selection. Fit mixtures of von Mises products and of bivariate sine models to synthetic Ramachandran-like data, and choose the number of components.
   • Work: EM for both families (Code 3.11, with Newton steps for \kappa); BIC and a held-out log-likelihood; data with 2 to 6 components, including correlated ones; the minimum-message-length approach of Amarasinghe and coauthors (2023) read and summarised.
   • Success criterion: the true number of components recovered in at least 8 of 10 synthetic data sets with 3000 points, and an explanation of the failures.
   • Deliverables: code, report, talk. Chapter 3.

Group projects

1. [Group of 2 · ≈ 30 h each] Paper study: arDCA against pseudo-likelihood. Reproduce, on synthetic families, the comparison of Trinquier and coauthors (2021) between autoregressive models and Potts models.
   • Data: Potts models with planted contacts on L = 30 to 60 sites and q = 4 to 21 (as in Code 4.6), sampled by Gibbs chains; MSAs of 500 to 5000 sequences.
   • Split: one member implements pseudo-likelihood and the contact score with APC; the other implements arDCA, its exact sampling and its likelihood; both evaluate the generated sequences.
   • Success criterion: contact precision against MSA size for both models, and the match of single and pair frequencies of generated sequences (correlation of connected pair frequencies above 0.9 for the larger MSAs).
   • Deliverables: code, a 6–8 page report, a 15-minute talk. Chapters 4 and 8.
2. [Group of 2 · ≈ 30 h each] Paper study: volumes of polytopes with hit-and-run and PDMPs. Reproduce the comparison of Chevallier, Cazals and Fearnhead (2022) on polytopes with known volumes.
   • Data: simplices, cross-polytopes and cubes cut by the half-space \sum_i x_i \leq t, whose exact volumes are known (for the cut cube, 2^n times the Irwin–Hall distribution function at (t + n)/2), in dimensions 5 to 40.
   • Split: one member implements multiphase Monte Carlo with hit-and-run (Code 6.9); the other implements the billiard PDMP (Definition 6.10, Exercise 19); both run the benchmark at equal numbers of oracle calls.
   • Success criterion: relative error against cost curves in each dimension, and a reasoned conclusion on which walk wins and from which dimension.
   • Deliverables: code, report, talk. Chapter 6.
3. [Group of 2 · ≈ 30 h each] Loop sampling with exact closure and reweighting. Implement tripeptide loop closure following Coutsias and coauthors (2004) and use it as a proposal for loop sampling.
   • Data: synthetic backbones built with Code 3.3, or a public PDB entry downloaded by the students, with one loop of 6 to 10 residues.
   • Split: one member implements the closure (numerical root finding of the closure equations is acceptable before the polynomial), and checks the at most 16 solutions; the other implements the weights of Proposition 3.7 and Section 8.6, the prior of Section 3.4 and a Metropolis–Hastings chain.
   • Success criterion: all solutions found for random tripeptides (checked against Code 5.4), and a weighted ensemble whose dihedral histograms match the target on a test case where it is known.
   • Deliverables: code, report, talk. Chapters 3, 5 and 8.

Written exam (proposal)

Format: 2 hours, one handwritten A4 sheet allowed, no electronic devices, four exercises.

1. Question (5 points). Explain the difference between an energy and a free energy, and why a protein state of higher energy can be more populated. Define local ergodicity, and explain what Code 0.9 shows about Monte Carlo estimates across a barrier.

   Expected: Definition 1.4 and Proposition 1.5, with F_A = \langle U \rangle_A - T S_A; a harmonic example such as Example 1.6. Definition 1.10 with the two inequalities, and the escape time of order e^{\beta \Delta U} (Theorem 1.11). For Code 0.9: the ergodic theorem holds, but with 6 transitions in 2 \times 10^5 steps the estimate of the second well is off by 30 %.

2. Question (5 points). (a) Align WHAT with WAT by Needleman–Wunsch (match 2, mismatch -1, gap -2). (b) Compute the Shapley–Shubik indices of the vote with quota 51 and weights 50, 49, 1. (c) Compute d_0 of the TM-score for L = 100.

   Expected: (a) score 4, WHAT over W-AT (Example 2.3). (b) (2/3, 1/6, 1/6) by enumerating the 6 orders (Example 7.10). (c) d_0 = 3.65 Å.

3. Question (5 points). Prove that the maximum-entropy distribution with given single and pair marginals is a Potts model, when such a Potts model exists. Then explain why pseudo-likelihood avoids the partition function, and why its loss is convex.

   Expected: Theorem 4.1, using that \mathbb{E}_P[\log P^\star] depends only on the marginals, so \mathrm{KL}(P \| P^\star) = S(P^\star) - S(P). Proposition 4.3: Z cancels in the conditional, and the loss is a sum of log-sum-exp terms of linear functions.

4. Question (5 points). Prove that hit-and-run leaves the uniform distribution of a convex body invariant. Then show that, in the multiphase scheme, each ratio \rho_i is at least 1/2, and deduce the order of the number of phases.

   Expected: Proposition 6.6, the symmetric kernel 2 / (|S^{n-1}|\, \ell(x, y)\, \|x - y\|^{n-1}) and detailed balance. Proposition 6.8: K_i \subset 2^{1/n} K_{i-1} by convexity, so \mathrm{vol}(K_i) \leq 2\,\mathrm{vol}(K_{i-1}), and m = \lceil n \log_2(R/r) \rceil phases.

Grading rubric (proposal)

References

Official references (master-mva.com, in the order of the course page)

General references listed on the course page: bioinformatics [10]; biophysics and theoretical biophysics [2], [8], [35]; algorithms and machine learning [9], [24], [36].

• [1] C. I. Branden and J. Tooze, Introduction to Protein Structure, Garland Science, 2012.
• [2] J. Kuriyan, B. Konforti and D. Wemmer, The Molecules of Life: Physical and Chemical Principles, Garland Science, 2012.
• [3] R. Peng, L.-A. Wu, Q. Wang, J. Qi and G. F. Gao, “Cell entry by SARS-CoV-2”, Trends in Biochemical Sciences 46(10):848–860, 2021.
• [4] S. Murakami, “Multidrug efflux transporter, AcrB — the pumping mechanism”, Current Opinion in Structural Biology 18(4):459–465, 2008.
• [5] A. Schmidt, H. Xu, A. Khan, T. O’Donnell, S. Khurana, L. King, J. Manischewitz, H. Golding, P. Suphaphiphat, A. Carfi, E. Settembre, P. Dormitzer, T. Kepler, R. Zhang, A. Moody, B. Haynes, H.-X. Liao, D. Shaw and S. Harrison, “Preconfiguration of the antigen-binding site during affinity maturation of a broadly neutralizing influenza virus antibody”, PNAS 110(1):264–269, 2013.
• [6] J. C. Schön and M. Jansen, “Prediction, determination and validation of phase diagrams via the global study of energy landscapes”, International Journal of Materials Research 100(2):135, 2009.
• [7] S. A. Adcock and J. A. McCammon, “Molecular dynamics: survey of methods for simulating the activity of proteins”, Chemical Reviews 106(5):1589–1615, 2006.
• [8] K. Dill and S. Bromberg, Molecular Driving Forces: Statistical Thermodynamics in Biology, Chemistry, Physics, and Nanoscience, Garland Science, 2010.
• [9] T. Lelièvre, M. Rousset and G. Stoltz, Free Energy Computations: A Mathematical Perspective, Imperial College Press / World Scientific, 2010.
• [10] J. Pevsner, Bioinformatics and Functional Genomics, John Wiley & Sons, 2015.
• [11] J. Jumper, R. Evans, A. Pritzel, T. Green, M. Figurnov, O. Ronneberger, K. Tunyasuvunakool, R. Bates, A. Žídek, A. Potapenko et al., “Highly accurate protein structure prediction with AlphaFold”, Nature 596(7873):583–589, 2021.
• [12] B. Moussad, R. Roche and D. Bhattacharya, “The transformative power of transformers in protein structure prediction”, PNAS 120(32):e2303499120, 2023.
• [13] T. O’Donnell and F. Cazals, “Modeling the dynamics of proteins: techniques from geometry and kinematics”, in preparation, 2024 (as listed on master-mva.com).
• [14] E. Coutsias, C. Seok, M. Jacobson and K. Dill, “A kinematic view of loop closure”, Journal of Computational Chemistry 25(4):510–528, 2004.
• [15] K. Noonan, D. O’Brien and J. Snoeyink, “Probik: protein backbone motion by inverse kinematics”, The International Journal of Robotics Research 24(11):971–982, 2005.
• [16] T. O’Donnell, V. Agashe and F. Cazals, “Geometric constraints within tripeptides and the existence of tripeptide reconstructions”, Journal of Computational Chemistry 44(13):1236–1249, 2023.
• [17] T. O’Donnell and F. Cazals, “Enhanced conformational exploration of protein loops using a global parameterization of the backbone geometry”, Journal of Computational Chemistry 44(11):1094–1104, 2023.
• [18] P. R. Amarasinghe, L. Allison, P. J. Stuckey, M. Garcia de la Banda, A. M. Lesk and A. S. Konagurthu, “Getting ‘ϕψχal’ with proteins: minimum message length inference of joint distributions of backbone and sidechain dihedral angles”, Bioinformatics 39(Supplement 1):i357–i367, 2023.
• [19] F. Morcos, A. Pagnani, B. Lunt, A. Bertolino, D. Marks, C. Sander, R. Zecchina, J. Onuchic, T. Hwa and M. Weigt, “Direct-coupling analysis of residue coevolution captures native contacts across many protein families”, PNAS 108(49):E1293–E1301, 2011.
• [20] J. Trinquier, G. Uguzzoni, A. Pagnani, F. Zamponi and M. Weigt, “Efficient generative modeling of protein sequences using simple autoregressive models”, Nature Communications 12(1):5800, 2021.
• [21] R. Rende, F. Gerace, A. Laio and S. Goldt, “Mapping of attention mechanisms to a generalized Potts model”, Physical Review Research 6(2):023057, 2024.
• [22] W. P. Russ, M. Figliuzzi, C. Stocker, P. Barrat-Charlaix, M. Socolich, P. Kast, D. Hilvert, R. Monasson, S. Cocco, M. Weigt et al., “An evolution-based model for designing chorismate mutase enzymes”, Science 369(6502):440–445, 2020.
• [23] S. Levy (ed.), Flavors of Geometry, Cambridge University Press, 1997.
• [24] A. Blum, J. Hopcroft and R. Kannan, Foundations of Data Science, Cambridge University Press, 2020.
• [25] H. Berbee, C. Boender, A. Ran, C. Scheffer, R. Smith and J. Telgen, “Hit-and-run algorithms for the identification of nonredundant linear inequalities”, Mathematical Programming 37(2):184–207, 1987.
• [26] L. Lovász, “Hit-and-run mixes fast”, Mathematical Programming, Series B 86:443–461, 1999.
• [27] B. Cousins and S. Vempala, “A practical volume algorithm”, Mathematical Programming Computation 8(2):133–160, 2016.
• [28] A. Chevallier, F. Cazals and P. Fearnhead, “Efficient computation of the volume of a polytope in high-dimensions using piecewise deterministic Markov processes”, AISTATS, 2022.
• [29] G. Fort, B. Jourdain, E. Kuhn, T. Lelièvre and G. Stoltz, “Convergence of the Wang–Landau algorithm”, Mathematics of Computation 84(295):2297–2327, 2015.
• [30] A. Chevallier and F. Cazals, “Wang–Landau algorithm: an adapted random walk to boost convergence”, Journal of Computational Physics 410:1–19, 2020.
• [31] S. Dasgupta and K. Sinha, “Randomized partition trees for exact nearest neighbor search”, JMLR Workshop and Conference Proceedings (COLT) 30:1–21, 2013.
• [32] S. Dasgupta and Y. Freund, “Random projection trees and low dimensional manifolds”, Proceedings of the 40th ACM Symposium on Theory of Computing, 537–546, 2008.
• [33] S. M. Lundberg and S.-I. Lee, “A unified approach to interpreting model predictions”, Advances in Neural Information Processing Systems 30, 2017.
• [34] S. M. Lundberg, G. Erion, H. Chen, A. DeGrave, J. M. Prutkin, B. Nair, R. Katz, J. Himmelfarb, N. Bansal and S.-I. Lee, “From local explanations to global understanding with explainable AI for trees”, Nature Machine Intelligence 2(1):56–67, 2020.
• [35] D. M. Zuckerman, Statistical Physics of Biomolecules: An Introduction, CRC Press, 2010.
• [36] T. Hastie, R. Tibshirani and J. Friedman, The Elements of Statistical Learning: Data Mining, Inference and Prediction, Springer, 2001.

Structures, alignment and structure prediction (ch. 1–2, 5)

• W. Kabsch, “A solution for the best rotation to relate two sets of vectors”, Acta Crystallographica A 32(5):922–923, 1976.
• S. B. Needleman and C. D. Wunsch, “A general method applicable to the search for similarities in the amino acid sequence of two proteins”, Journal of Molecular Biology 48(3):443–453, 1970.
• S. Henikoff and J. G. Henikoff, “Amino acid substitution matrices from protein blocks”, PNAS 89(22):10915–10919, 1992.
• A. Zemla, “LGA: a method for finding 3D similarities in protein structures”, Nucleic Acids Research 31(13):3370–3374, 2003.
• Y. Zhang and J. Skolnick, “Scoring function for automated assessment of protein structure template quality”, Proteins 57(4):702–710, 2004.
• J. Xu and Y. Zhang, “How significant is a protein structure similarity with TM-score = 0.5?”, Bioinformatics 26(7):889–895, 2010.
• V. Mariani, M. Biasini, A. Barbato and T. Schwede, “lDDT: a local superposition-free score for comparing protein structures and models using distance difference tests”, Bioinformatics 29(21):2722–2728, 2013.
• M. Mirdita, K. Schütze, Y. Moriwaki, L. Heo, S. Ovchinnikov and M. Steinegger, “ColabFold: making protein folding accessible to all”, Nature Methods 19:679–682, 2022.
• A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N. Gomez, Ł. Kaiser and I. Polosukhin, “Attention is all you need”, Advances in Neural Information Processing Systems 30, 2017.

Dynamics, kinematics and directional statistics (ch. 1, 3, 8)

• H. A. Kramers, “Brownian motion in a field of force and the diffusion model of chemical reactions”, Physica 7(4):284–304, 1940.
• A. Bovier, M. Eckhoff, V. Gayrard and M. Klein, “Metastability in reversible diffusion processes I: sharp asymptotics for capacities and exit times”, Journal of the European Mathematical Society 6(4):399–424, 2004.
• E. T. Jaynes, “Information theory and statistical mechanics”, Physical Review 106(4):620–630, 1957.
• G. N. Ramachandran, C. Ramakrishnan and V. Sasisekharan, “Stereochemistry of polypeptide chain configurations”, Journal of Molecular Biology 7(1):95–99, 1963.
• R. A. Engh and R. Huber, “Accurate bond and angle parameters for X-ray protein structure refinement”, Acta Crystallographica A 47(4):392–400, 1991.
• J. Parsons, J. B. Holmes, J. M. Rojas, J. Tsai and C. E. M. Strauss, “Practical conversion from torsion space to Cartesian space for in silico protein synthesis”, Journal of Computational Chemistry 26(10):1063–1068, 2005.
• N. Go and H. A. Scheraga, “Ring closure and local conformational deformations of chain molecules”, Macromolecules 3(2):178–187, 1970.
• A. A. Canutescu and R. L. Dunbrack, “Cyclic coordinate descent: a robotics algorithm for protein loop closure”, Protein Science 12(5):963–972, 2003.
• H. Singh, V. Hnizdo and E. Demchuk, “Probabilistic model for two dependent circular variables”, Biometrika 89(3):719–723, 2002.
• A. Banerjee, I. S. Dhillon, J. Ghosh and S. Sra, “Clustering on the unit hypersphere using von Mises–Fisher distributions”, Journal of Machine Learning Research 6:1345–1382, 2005.
• K. V. Mardia and P. E. Jupp, Directional Statistics, Wiley, 2000.

Sequence models and samplers (ch. 4)

• J. Besag, “Statistical analysis of non-lattice data”, Journal of the Royal Statistical Society, Series D (The Statistician) 24(3):179–195, 1975.
• M. Weigt, R. A. White, H. Szurmant, J. A. Hoch and T. Hwa, “Identification of direct residue contacts in protein–protein interaction by message passing”, PNAS 106(1):67–72, 2009.
• S. D. Dunn, L. M. Wahl and G. B. Gloor, “Mutual information without the influence of phylogeny or entropy dramatically improves residue contact prediction”, Bioinformatics 24(3):333–340, 2008.
• M. Ekeberg, C. Lövkvist, Y. Lan, M. Weigt and E. Aurell, “Improved contact prediction in proteins: using pseudolikelihoods to infer Potts models”, Physical Review E 87(1):012707, 2013.
• A. Rives, J. Meier, T. Sercu, S. Goyal, Z. Lin, J. Liu, D. Guo, M. Ott, C. L. Zitnick, J. Ma and R. Fergus, “Biological structure and function emerge from scaling unsupervised learning to 250 million protein sequences”, PNAS 118(15):e2016239118, 2021.
• R. Rao, J. Meier, T. Sercu, S. Ovchinnikov and A. Rives, “Transformer protein language models are unsupervised structure learners”, International Conference on Learning Representations, 2021.
• F. Noé, S. Olsson, J. Köhler and H. Wu, “Boltzmann generators: sampling equilibrium states of many-body systems with deep learning”, Science 365(6457):eaaw1147, 2019.

High-dimensional sampling and partitions (ch. 6–7)

• G. Elekes, “A geometric inequality and the complexity of computing volume”, Discrete & Computational Geometry 1(4):289–292, 1986.
• M. E. Dyer and A. M. Frieze, “On the complexity of computing the volume of a polyhedron”, SIAM Journal on Computing 17(5):967–974, 1988.
• M. Dyer, A. Frieze and R. Kannan, “A random polynomial-time algorithm for approximating the volume of convex bodies”, Journal of the ACM 38(1):1–17, 1991.
• L. Lovász and S. Vempala, “Simulated annealing in convex bodies and an O^*(n^4) volume algorithm”, Journal of Computer and System Sciences 72(2):392–417, 2006.
• R. L. Smith, “Efficient Monte Carlo procedures for generating points uniformly distributed over bounded regions”, Operations Research 32(6):1296–1308, 1984.
• P. Fearnhead, J. Bierkens, M. Pollock and G. O. Roberts, “Piecewise deterministic Markov processes for continuous-time Monte Carlo”, Statistical Science 33(3):386–412, 2018.
• F. Wang and D. P. Landau, “Efficient, multiple-range random walk algorithm to calculate the density of states”, Physical Review Letters 86(10):2050–2053, 2001.
• R. E. Belardinelli and V. D. Pereyra, “Fast algorithm to calculate density of states”, Physical Review E 75(4):046701, 2007.
• S. Kpotufe and S. Dasgupta, “A tree-based regressor that adapts to intrinsic dimension”, Journal of Computer and System Sciences 78(5):1496–1515, 2012.
• L. S. Shapley, “A value for n-person games”, in Contributions to the Theory of Games II, Annals of Mathematics Studies 28, Princeton University Press, 307–317, 1953.
• L. S. Shapley and M. Shubik, “A method for evaluating the distribution of power in a committee system”, American Political Science Review 48(3):787–792, 1954.

Research directions (ch. 9)

• D. del Alamo, D. Sala, H. S. Mchaourab and J. Meiler, “Sampling alternative conformational states of transporters and receptors with AlphaFold2”, eLife 11:e75751, 2022.
• H. K. Wayment-Steele, A. Ojoawo, R. Otten, J. M. Apitz, W. Pitsawong, M. Hömberger, S. Ovchinnikov, L. Colwell and D. Kern, “Predicting multiple conformations via sequence clustering and AlphaFold2”, Nature 625:832–839, 2024.
• J. Abramson, J. Adler, J. Dunger, R. Evans, T. Green, A. Pritzel, O. Ronneberger, L. Willmore, A. J. Ballard, J. Bambrick et al., “Accurate structure prediction of biomolecular interactions with AlphaFold 3”, Nature 630:493–500, 2024.
• D. E. Shaw, P. Maragakis, K. Lindorff-Larsen, S. Piana, R. O. Dror, M. P. Eastwood, J. A. Bank, J. M. Jumper, J. K. Salmon, Y. Shan and W. Wriggers, “Atomic-level characterization of the structural dynamics of proteins”, Science 330(6002):341–346, 2010.
• J.-H. Prinz, H. Wu, M. Sarich, B. Keller, M. Senne, M. Held, J. D. Chodera, C. Schütte and F. Noé, “Markov models of molecular kinetics: generation and validation”, Journal of Chemical Physics 134(17):174105, 2011.
• A. Laio and M. Parrinello, “Escaping free-energy minima”, PNAS 99(20):12562–12566, 2002.
• J. L. Watson, D. Juergens, N. R. Bennett, B. L. Trippe, J. Yim, H. E. Eisenach, W. Ahern, A. J. Borst, R. J. Ragotte, L. F. Milles et al., “De novo design of protein structure and function with RFdiffusion”, Nature 620:1089–1100, 2023.
• J. Dauparas, I. Anishchenko, N. Bennett, H. Bai, R. J. Ragotte, L. F. Milles, B. I. M. Wicky, A. Courbet, R. J. de Haas, N. Bethel et al., “Robust deep learning-based protein sequence design using ProteinMPNN”, Science 378(6615):49–56, 2022.
• Z. Lin, H. Akin, R. Rao, B. Hie, Z. Zhu, W. Lu, N. Smetanin, R. Verkuil, O. Kabeli, Y. Shmueli et al., “Evolutionary-scale prediction of atomic-level protein structure with a language model”, Science 379(6637):1123–1130, 2023.
• K. Xia and G.-W. Wei, “Persistent homology analysis of protein structure, flexibility, and folding”, International Journal for Numerical Methods in Biomedical Engineering 30(8):814–844, 2014.