Single-Enzyme MWC Port: Laplace vs HMC on Real CBS Allostery¶

Notebook 13 used a cheap Gauss–Newton Laplace approximation for the posterior, with the caveat that the real tool (Maud) uses HMC, and that Laplace is only a local Gaussian. This notebook closes that gap on a faithful, real-units, single-enzyme problem and asks directly: when does the cheap Laplace posterior agree with HMC, and when does it fail?

We port one enzyme — cystathionine β-synthase (CBS), the SAM-activated allosteric switch — using its actual rate law and priors from the Biosustain methionine model.

Real Maud rate law (modular rate law with MWC allostery). With activator SAM (amet), n subunits, transfer constant L, dissociation constant Kd:

$$v_\text{CBS} = V_\max \cdot \underbrace{\frac{1}{1 + L\,\big(1/(1+[\text{SAM}]/K_d)\big)^{n}}}_{\text{fraction in active (R) state}}$$

At low SAM the enzyme sits in the inactive state (v = Vmax/(1+L)); SAM shifts it to the active state (v → Vmax). Real Maud priors (priors.toml), lognormal:

parameter	Maud id	prior
dissociation const Kd	CBS1/amet	median 9.22e-5 mM, scale 2
transfer const L	CBS1	median 1, scale 2

This is a genuine in-vitro-style identifiability problem: measure CBS rate vs [SAM] and infer (Kd, L). Because it is 2-D we can compute the exact posterior on a grid as an unimpeachable referee for both Laplace and HMC. (Substrate is taken saturating and Vmax calibrated, so the two allosteric constants are the unknowns; the same machinery extends to Vmax as a nuisance.)

In [1]:

import time, torch
import numpy as np
import matplotlib.pyplot as plt
import matplotlib
matplotlib.rcParams['figure.dpi']=120
torch.set_default_dtype(torch.float64)

N_SUB = 2                 # CBS subunits (MWC)
Vmax  = 1.0               # normalized kcat*[E] at saturating substrate
def cbs_rate(amet, Kd, L):
    return Vmax / (1.0 + L*(1.0/(1.0 + amet/Kd))**N_SUB)

# real Maud lognormal priors (mM units)
Kd_med, Kd_logsd = 9.22e-5, 2.0
L_med,  L_logsd  = 1.0, 2.0
prior_mean = torch.tensor([np.log(Kd_med), np.log(L_med)])
prior_sd   = torch.tensor([Kd_logsd, L_logsd])
Kd_true, L_true = 9.22e-5, 3.0          # truth (within prior)
NAMES=['Kd (mM)','L']
print('CBS MWC model + real Maud priors ready;  true Kd=%.2e mM, L=%.1f' % (Kd_true, L_true))

import time, torch import numpy as np import matplotlib.pyplot as plt import matplotlib matplotlib.rcParams['figure.dpi']=120 torch.set_default_dtype(torch.float64) N_SUB = 2 # CBS subunits (MWC) Vmax = 1.0 # normalized kcat*[E] at saturating substrate def cbs_rate(amet, Kd, L): return Vmax / (1.0 + L*(1.0/(1.0 + amet/Kd))**N_SUB) # real Maud lognormal priors (mM units) Kd_med, Kd_logsd = 9.22e-5, 2.0 L_med, L_logsd = 1.0, 2.0 prior_mean = torch.tensor([np.log(Kd_med), np.log(L_med)]) prior_sd = torch.tensor([Kd_logsd, L_logsd]) Kd_true, L_true = 9.22e-5, 3.0 # truth (within prior) NAMES=['Kd (mM)','L'] print('CBS MWC model + real Maud priors ready; true Kd=%.2e mM, L=%.1f' % (Kd_true, L_true))

CBS MWC model + real Maud priors ready;  true Kd=9.22e-05 mM, L=3.0

Inference machinery: MAP, Laplace, exact grid, and preconditioned HMC¶

phi = [log Kd, log L]. The model is closed-form (no steady-state solve here — this is the nonlinear rate law that sits on top of the linear stoichiometric scaffold in the network notebooks), so gradients are direct. HMC is preconditioned with the Laplace Hessian as its mass matrix so it mixes well on the stiff posterior.

In [2]:

def make_nlp(amet, data, sigma):
    def nlp(phi):
        Kd, L = torch.exp(phi[0]), torch.exp(phi[1])
        pred = cbs_rate(amet, Kd, L)
        return 0.5*(((pred-data)/sigma)**2).sum() + 0.5*(((phi-prior_mean)/prior_sd)**2).sum()
    return nlp

def fit_map_laplace(nlp):
    phi = prior_mean.clone().requires_grad_(True)
    opt = torch.optim.Adam([phi], lr=0.05)
    for _ in range(3000):
        opt.zero_grad(); v=nlp(phi); v.backward(); opt.step()
    pmap = phi.detach()
    H = torch.autograd.functional.hessian(nlp, pmap)
    cov = torch.linalg.inv(H); sd = torch.sqrt(torch.diag(cov))
    return pmap, cov, sd, H

def grid_posterior(nlp, pmap, half=(3.5,5.0), nx=240):
    gk = torch.linspace(pmap[0]-half[0], pmap[0]+half[0], nx)
    gl = torch.linspace(pmap[1]-half[1], pmap[1]+half[1], nx)
    GK, GL = torch.meshgrid(gk, gl, indexing='ij'); lp=torch.empty_like(GK)
    for i in range(nx):
        for j in range(nx):
            lp[i,j] = -nlp(torch.stack([GK[i,j], GL[i,j]]))
    P = torch.exp(lp-lp.max()); P/=P.sum()
    return gk, gl, P

def hmc(nlp, pmap, H, n=6000, eps=0.6, Lf=12, seed=1):
    torch.manual_seed(seed)
    cholH = torch.linalg.cholesky(H); Minv = torch.linalg.inv(H)
    def gradU(q):
        q=q.clone().requires_grad_(True); g,=torch.autograd.grad(nlp(q), q); return g
    Kin = lambda p: 0.5*(p@(Minv@p))
    cur=pmap.clone(); curU=nlp(cur); samples=[]; acc=0
    for t in range(n):
        p0 = cholH@torch.randn(2); p=p0.clone(); q=cur.clone()
        p = p - 0.5*eps*gradU(q)
        for i in range(Lf):
            q = q + eps*(Minv@p); p = p - eps*gradU(q)
        p = p + 0.5*eps*gradU(q)
        dH = (nlp(q)+Kin(p)) - (curU+Kin(p0))
        if torch.log(torch.rand(1)) < -dH:
            cur=q.detach(); curU=nlp(cur); acc+=1
        samples.append(cur.clone())
    return torch.stack(samples)[n//6:], acc/n

print('inference machinery ready')

def make_nlp(amet, data, sigma): def nlp(phi): Kd, L = torch.exp(phi[0]), torch.exp(phi[1]) pred = cbs_rate(amet, Kd, L) return 0.5*(((pred-data)/sigma)**2).sum() + 0.5*(((phi-prior_mean)/prior_sd)**2).sum() return nlp def fit_map_laplace(nlp): phi = prior_mean.clone().requires_grad_(True) opt = torch.optim.Adam([phi], lr=0.05) for _ in range(3000): opt.zero_grad(); v=nlp(phi); v.backward(); opt.step() pmap = phi.detach() H = torch.autograd.functional.hessian(nlp, pmap) cov = torch.linalg.inv(H); sd = torch.sqrt(torch.diag(cov)) return pmap, cov, sd, H def grid_posterior(nlp, pmap, half=(3.5,5.0), nx=240): gk = torch.linspace(pmap[0]-half[0], pmap[0]+half[0], nx) gl = torch.linspace(pmap[1]-half[1], pmap[1]+half[1], nx) GK, GL = torch.meshgrid(gk, gl, indexing='ij'); lp=torch.empty_like(GK) for i in range(nx): for j in range(nx): lp[i,j] = -nlp(torch.stack([GK[i,j], GL[i,j]])) P = torch.exp(lp-lp.max()); P/=P.sum() return gk, gl, P def hmc(nlp, pmap, H, n=6000, eps=0.6, Lf=12, seed=1): torch.manual_seed(seed) cholH = torch.linalg.cholesky(H); Minv = torch.linalg.inv(H) def gradU(q): q=q.clone().requires_grad_(True); g,=torch.autograd.grad(nlp(q), q); return g Kin = lambda p: 0.5*(p@(Minv@p)) cur=pmap.clone(); curU=nlp(cur); samples=[]; acc=0 for t in range(n): p0 = cholH@torch.randn(2); p=p0.clone(); q=cur.clone() p = p - 0.5*eps*gradU(q) for i in range(Lf): q = q + eps*(Minv@p); p = p - eps*gradU(q) p = p + 0.5*eps*gradU(q) dH = (nlp(q)+Kin(p)) - (curU+Kin(p0)) if torch.log(torch.rand(1)) < -dH: cur=q.detach(); curU=nlp(cur); acc+=1 samples.append(cur.clone()) return torch.stack(samples)[n//6:], acc/n print('inference machinery ready')

inference machinery ready

Scenario A — informative dose-response (SAM straddles Kd)¶

12 measurements log-spaced across Kd, 4% noise. This is the well-designed experiment.

In [3]:

amet_A = torch.logspace(-5.3, -2.7, 12)        # ~5e-6 .. 2e-3 mM
torch.manual_seed(2)
cleanA = cbs_rate(amet_A, Kd_true, L_true)
dataA  = cleanA*(1 + 0.04*torch.randn(12)); sigA = 0.04*cleanA

nlpA = make_nlp(amet_A, dataA, sigA)
pmapA, covA, sdA, HA = fit_map_laplace(nlpA)
gkA, glA, PA = grid_posterior(nlpA, pmapA)
t0=time.perf_counter(); chainA, accA = hmc(nlpA, pmapA, HA); thmc=time.perf_counter()-t0

# grid marginal sd
mkA=(PA.sum(1)*gkA).sum(); skA=torch.sqrt((PA.sum(1)*(gkA-mkA)**2).sum())
mlA=(PA.sum(0)*glA).sum(); slA=torch.sqrt((PA.sum(0)*(glA-mlA)**2).sum())
print('SCENARIO A  (12 doses straddling Kd)   HMC accept %.2f, %.0fs' % (accA, thmc))
print('  MAP     : Kd=%.2e  L=%.2f   (true %.2e, %.1f)' % (torch.exp(pmapA[0]),torch.exp(pmapA[1]),Kd_true,L_true))
print('  log-sd   Laplace [%.3f, %.3f]' % (sdA[0], sdA[1]))
print('           grid    [%.3f, %.3f]' % (skA, slA))
print('           HMC     [%.3f, %.3f]' % (chainA[:,0].std(), chainA[:,1].std()))
print('  logL skew (grid) = %.2f   ->  near-Gaussian, Laplace is accurate' %
      ((PA.sum(0)*((glA-mlA)/slA)**3).sum()))

amet_A = torch.logspace(-5.3, -2.7, 12) # ~5e-6 .. 2e-3 mM torch.manual_seed(2) cleanA = cbs_rate(amet_A, Kd_true, L_true) dataA = cleanA*(1 + 0.04*torch.randn(12)); sigA = 0.04*cleanA nlpA = make_nlp(amet_A, dataA, sigA) pmapA, covA, sdA, HA = fit_map_laplace(nlpA) gkA, glA, PA = grid_posterior(nlpA, pmapA) t0=time.perf_counter(); chainA, accA = hmc(nlpA, pmapA, HA); thmc=time.perf_counter()-t0 # grid marginal sd mkA=(PA.sum(1)*gkA).sum(); skA=torch.sqrt((PA.sum(1)*(gkA-mkA)**2).sum()) mlA=(PA.sum(0)*glA).sum(); slA=torch.sqrt((PA.sum(0)*(glA-mlA)**2).sum()) print('SCENARIO A (12 doses straddling Kd) HMC accept %.2f, %.0fs' % (accA, thmc)) print(' MAP : Kd=%.2e L=%.2f (true %.2e, %.1f)' % (torch.exp(pmapA[0]),torch.exp(pmapA[1]),Kd_true,L_true)) print(' log-sd Laplace [%.3f, %.3f]' % (sdA[0], sdA[1])) print(' grid [%.3f, %.3f]' % (skA, slA)) print(' HMC [%.3f, %.3f]' % (chainA[:,0].std(), chainA[:,1].std())) print(' logL skew (grid) = %.2f -> near-Gaussian, Laplace is accurate' % ((PA.sum(0)*((glA-mlA)/slA)**3).sum()))

SCENARIO A  (12 doses straddling Kd)   HMC accept 0.96, 47s
  MAP     : Kd=9.04e-05  L=3.04   (true 9.22e-05, 3.0)
  log-sd   Laplace [0.074, 0.039]
           grid    [0.074, 0.039]
           HMC     [0.074, 0.040]
  logL skew (grid) = 0.07   ->  near-Gaussian, Laplace is accurate

Scenario B — data-poor (only high SAM, few points)¶

6 measurements, all above Kd, so the low-SAM plateau that pins L is never observed. This is the realistic information-poor case.

In [4]:

amet_B = torch.logspace(-4.0, -2.7, 6)
torch.manual_seed(2)
cleanB = cbs_rate(amet_B, Kd_true, L_true)
dataB  = cleanB*(1 + 0.04*torch.randn(6)); sigB = 0.04*cleanB

nlpB = make_nlp(amet_B, dataB, sigB)
pmapB, covB, sdB, HB = fit_map_laplace(nlpB)
gkB, glB, PB = grid_posterior(nlpB, pmapB, half=(3.0,6.0))
t0=time.perf_counter(); chainB, accB = hmc(nlpB, pmapB, HB, eps=0.5); thmcB=time.perf_counter()-t0

PlB = PB.sum(0); mlB=(PlB*glB).sum(); slB=torch.sqrt((PlB*(glB-mlB)**2).sum())
skewB=(PlB*((glB-mlB)/slB)**3).sum()
cdf=torch.cumsum(PlB,0)
loL=torch.exp(glB[(cdf>0.025).nonzero()[0,0]]); hiL=torch.exp(glB[(cdf>0.975).nonzero()[0,0]])
lapLo=torch.exp(pmapB[1]-1.96*sdB[1]); lapHi=torch.exp(pmapB[1]+1.96*sdB[1])
print('SCENARIO B  (6 doses, all above Kd)    HMC accept %.2f, %.0fs' % (accB, thmcB))
print('  L posterior skew (grid) = %.2f   ->  heavy upper tail, non-Gaussian' % skewB)
print('  L 95%% CI  grid/HMC : [%.2f, %.2f]   (true L=%.1f)' % (loL, hiL, L_true))
print('  L 95%% CI  Laplace  : [%.2f, %.2f]   <- underestimates the upper tail' % (lapLo, lapHi))

amet_B = torch.logspace(-4.0, -2.7, 6) torch.manual_seed(2) cleanB = cbs_rate(amet_B, Kd_true, L_true) dataB = cleanB*(1 + 0.04*torch.randn(6)); sigB = 0.04*cleanB nlpB = make_nlp(amet_B, dataB, sigB) pmapB, covB, sdB, HB = fit_map_laplace(nlpB) gkB, glB, PB = grid_posterior(nlpB, pmapB, half=(3.0,6.0)) t0=time.perf_counter(); chainB, accB = hmc(nlpB, pmapB, HB, eps=0.5); thmcB=time.perf_counter()-t0 PlB = PB.sum(0); mlB=(PlB*glB).sum(); slB=torch.sqrt((PlB*(glB-mlB)**2).sum()) skewB=(PlB*((glB-mlB)/slB)**3).sum() cdf=torch.cumsum(PlB,0) loL=torch.exp(glB[(cdf>0.025).nonzero()[0,0]]); hiL=torch.exp(glB[(cdf>0.975).nonzero()[0,0]]) lapLo=torch.exp(pmapB[1]-1.96*sdB[1]); lapHi=torch.exp(pmapB[1]+1.96*sdB[1]) print('SCENARIO B (6 doses, all above Kd) HMC accept %.2f, %.0fs' % (accB, thmcB)) print(' L posterior skew (grid) = %.2f -> heavy upper tail, non-Gaussian' % skewB) print(' L 95%% CI grid/HMC : [%.2f, %.2f] (true L=%.1f)' % (loL, hiL, L_true)) print(' L 95%% CI Laplace : [%.2f, %.2f] <- underestimates the upper tail' % (lapLo, lapHi))

SCENARIO B  (6 doses, all above Kd)    HMC accept 0.97, 48s
  L posterior skew (grid) = 1.07   ->  heavy upper tail, non-Gaussian
  L 95% CI  grid/HMC : [1.14, 6.95]   (true L=3.0)
  L 95% CI  Laplace  : [0.95, 4.35]   <- underestimates the upper tail

In [5]:

def ellipse(pmap, cov, nsd=2.45):
    vals,vecs=torch.linalg.eigh(cov); th=torch.linspace(0,2*np.pi,120)
    pts=vecs@torch.diag(nsd*torch.sqrt(vals))@torch.stack([torch.cos(th),torch.sin(th)])
    return torch.exp(pmap[0]+pts[0]), torch.exp(pmap[1]+pts[1])

fig, ax = plt.subplots(2, 2, figsize=(13, 9))

# (1) dose-response fit (A) with posterior predictive band
ad = torch.logspace(-5.6,-2.5,200)
samp = chainA[torch.randint(0,len(chainA),(200,))]
pp = torch.stack([cbs_rate(ad, torch.exp(s[0]), torch.exp(s[1])) for s in samp])
ax[0,0].fill_between(ad.numpy(), pp.quantile(0.025,0).numpy(), pp.quantile(0.975,0).numpy(),
                     alpha=0.25, color='tab:blue', label='95% predictive')
ax[0,0].plot(ad, cbs_rate(ad,Kd_true,L_true), 'r-', lw=2, label='true')
ax[0,0].errorbar(amet_A.numpy(), dataA.numpy(), yerr=sigA.numpy(), fmt='ko', ms=4, label='data (A)')
ax[0,0].set_xscale('log'); ax[0,0].axvline(Kd_true, ls=':', color='gray')
ax[0,0].set_xlabel('[SAM] (mM)'); ax[0,0].set_ylabel('CBS rate / Vmax')
ax[0,0].set_title('CBS allosteric activation (scenario A fit)'); ax[0,0].legend(fontsize=8)

# (2) scenario A posterior: grid + HMC + Laplace all agree
ax[0,1].contour(torch.exp(gkA), torch.exp(glA), PA.T, levels=6, cmap='Greys')
ax[0,1].scatter(torch.exp(chainA[:,0]), torch.exp(chainA[:,1]), s=2, alpha=0.15, color='tab:blue', label='HMC')
ex,ey=ellipse(pmapA,covA); ax[0,1].plot(ex,ey,'r-',lw=2,label='Laplace 95%')
ax[0,1].plot(Kd_true,L_true,'r*',ms=15)
ax[0,1].set_xlim(torch.exp(chainA[:,0]).quantile(.003), torch.exp(chainA[:,0]).quantile(.997))
ax[0,1].set_ylim(torch.exp(chainA[:,1]).quantile(.003), torch.exp(chainA[:,1]).quantile(.997))
ax[0,1].set_xlabel('Kd (mM)'); ax[0,1].set_ylabel('L')
ax[0,1].set_title('Scenario A: grid = HMC = Laplace'); ax[0,1].legend(fontsize=8)

# (3) scenario B posterior: skewed; Laplace ellipse misses the tail
ax[1,0].contour(torch.exp(gkB), torch.exp(glB), PB.T, levels=6, cmap='Greys')
ax[1,0].scatter(torch.exp(chainB[:,0]), torch.exp(chainB[:,1]), s=2, alpha=0.15, color='tab:orange', label='HMC')
ex,ey=ellipse(pmapB,covB); ax[1,0].plot(ex,ey,'r-',lw=2,label='Laplace 95%')
ax[1,0].plot(Kd_true,L_true,'r*',ms=15)
ax[1,0].set_xlim(torch.exp(chainB[:,0]).quantile(.005), torch.exp(chainB[:,0]).quantile(.99))
ax[1,0].set_ylim(0, torch.exp(chainB[:,1]).quantile(.99))
ax[1,0].set_xlabel('Kd (mM)'); ax[1,0].set_ylabel('L')
ax[1,0].set_title('Scenario B: HMC banana, Laplace ellipse misses it'); ax[1,0].legend(fontsize=8)

# (4) scenario B marginal of L: grid/HMC skewed vs Laplace Gaussian (densities in L-space)
Lgrid = torch.exp(glB)
dens_logL = PlB/(PlB.sum()*(glB[1]-glB[0]))          # p(log L)
ax[1,1].plot(Lgrid, dens_logL/Lgrid, 'k-', lw=2, label='exact grid')   # /L Jacobian -> p(L)
ax[1,1].hist(torch.exp(chainB[:,1]).numpy(), bins=80, density=True, alpha=0.4, color='tab:orange', label='HMC')
lap = torch.exp(-0.5*((glB-pmapB[1])/sdB[1])**2); lap=lap/(lap.sum()*(glB[1]-glB[0]))
ax[1,1].plot(Lgrid, lap/Lgrid, 'r--', lw=2, label='Laplace (Gaussian-in-log)')
ax[1,1].axvline(L_true, color='g', ls=':', lw=2, label='true L')
ax[1,1].set_xlim(0,12); ax[1,1].set_xlabel('L'); ax[1,1].set_yticks([])
ax[1,1].set_title('Scenario B: Laplace misses L upper tail'); ax[1,1].legend(fontsize=8)

plt.suptitle('CBS MWC allostery: Laplace vs HMC vs exact posterior (real Maud parameterization)',
             fontsize=13, fontweight='bold')
plt.tight_layout(); plt.savefig('cbs_mwc_laplace_vs_hmc.png', dpi=120, bbox_inches='tight'); plt.show()

def ellipse(pmap, cov, nsd=2.45): vals,vecs=torch.linalg.eigh(cov); th=torch.linspace(0,2*np.pi,120) pts=vecs@torch.diag(nsd*torch.sqrt(vals))@torch.stack([torch.cos(th),torch.sin(th)]) return torch.exp(pmap[0]+pts[0]), torch.exp(pmap[1]+pts[1]) fig, ax = plt.subplots(2, 2, figsize=(13, 9)) # (1) dose-response fit (A) with posterior predictive band ad = torch.logspace(-5.6,-2.5,200) samp = chainA[torch.randint(0,len(chainA),(200,))] pp = torch.stack([cbs_rate(ad, torch.exp(s[0]), torch.exp(s[1])) for s in samp]) ax[0,0].fill_between(ad.numpy(), pp.quantile(0.025,0).numpy(), pp.quantile(0.975,0).numpy(), alpha=0.25, color='tab:blue', label='95% predictive') ax[0,0].plot(ad, cbs_rate(ad,Kd_true,L_true), 'r-', lw=2, label='true') ax[0,0].errorbar(amet_A.numpy(), dataA.numpy(), yerr=sigA.numpy(), fmt='ko', ms=4, label='data (A)') ax[0,0].set_xscale('log'); ax[0,0].axvline(Kd_true, ls=':', color='gray') ax[0,0].set_xlabel('[SAM] (mM)'); ax[0,0].set_ylabel('CBS rate / Vmax') ax[0,0].set_title('CBS allosteric activation (scenario A fit)'); ax[0,0].legend(fontsize=8) # (2) scenario A posterior: grid + HMC + Laplace all agree ax[0,1].contour(torch.exp(gkA), torch.exp(glA), PA.T, levels=6, cmap='Greys') ax[0,1].scatter(torch.exp(chainA[:,0]), torch.exp(chainA[:,1]), s=2, alpha=0.15, color='tab:blue', label='HMC') ex,ey=ellipse(pmapA,covA); ax[0,1].plot(ex,ey,'r-',lw=2,label='Laplace 95%') ax[0,1].plot(Kd_true,L_true,'r*',ms=15) ax[0,1].set_xlim(torch.exp(chainA[:,0]).quantile(.003), torch.exp(chainA[:,0]).quantile(.997)) ax[0,1].set_ylim(torch.exp(chainA[:,1]).quantile(.003), torch.exp(chainA[:,1]).quantile(.997)) ax[0,1].set_xlabel('Kd (mM)'); ax[0,1].set_ylabel('L') ax[0,1].set_title('Scenario A: grid = HMC = Laplace'); ax[0,1].legend(fontsize=8) # (3) scenario B posterior: skewed; Laplace ellipse misses the tail ax[1,0].contour(torch.exp(gkB), torch.exp(glB), PB.T, levels=6, cmap='Greys') ax[1,0].scatter(torch.exp(chainB[:,0]), torch.exp(chainB[:,1]), s=2, alpha=0.15, color='tab:orange', label='HMC') ex,ey=ellipse(pmapB,covB); ax[1,0].plot(ex,ey,'r-',lw=2,label='Laplace 95%') ax[1,0].plot(Kd_true,L_true,'r*',ms=15) ax[1,0].set_xlim(torch.exp(chainB[:,0]).quantile(.005), torch.exp(chainB[:,0]).quantile(.99)) ax[1,0].set_ylim(0, torch.exp(chainB[:,1]).quantile(.99)) ax[1,0].set_xlabel('Kd (mM)'); ax[1,0].set_ylabel('L') ax[1,0].set_title('Scenario B: HMC banana, Laplace ellipse misses it'); ax[1,0].legend(fontsize=8) # (4) scenario B marginal of L: grid/HMC skewed vs Laplace Gaussian (densities in L-space) Lgrid = torch.exp(glB) dens_logL = PlB/(PlB.sum()*(glB[1]-glB[0])) # p(log L) ax[1,1].plot(Lgrid, dens_logL/Lgrid, 'k-', lw=2, label='exact grid') # /L Jacobian -> p(L) ax[1,1].hist(torch.exp(chainB[:,1]).numpy(), bins=80, density=True, alpha=0.4, color='tab:orange', label='HMC') lap = torch.exp(-0.5*((glB-pmapB[1])/sdB[1])**2); lap=lap/(lap.sum()*(glB[1]-glB[0])) ax[1,1].plot(Lgrid, lap/Lgrid, 'r--', lw=2, label='Laplace (Gaussian-in-log)') ax[1,1].axvline(L_true, color='g', ls=':', lw=2, label='true L') ax[1,1].set_xlim(0,12); ax[1,1].set_xlabel('L'); ax[1,1].set_yticks([]) ax[1,1].set_title('Scenario B: Laplace misses L upper tail'); ax[1,1].legend(fontsize=8) plt.suptitle('CBS MWC allostery: Laplace vs HMC vs exact posterior (real Maud parameterization)', fontsize=13, fontweight='bold') plt.tight_layout(); plt.savefig('cbs_mwc_laplace_vs_hmc.png', dpi=120, bbox_inches='tight'); plt.show()

Takeaways and caveats¶

The gap is closed, both ways.

• Laplace is validated where it should hold. On the well-designed dose-response (scenario A), the cheap Gauss-Newton/Laplace posterior used in notebook 13 matches both the exact grid and HMC to three decimals. So the fast method is trustworthy when the data are informative and the posterior is near-Gaussian-in-log.
• Laplace fails where Maud's HMC earns its keep. When the data don't excite the regulator across its active range (scenario B), the posterior on the transfer constant L is genuinely skewed with a heavy upper tail; Laplace's symmetric ellipse underestimates the 95% upper bound (≈4.4 vs the true-covering ≈6.8 from grid/HMC). This is exactly the regime where you must use full posterior sampling — and it is the same identifiability story as notebooks 12–13, now seen as posterior shape, not just width.

This is a faithful port — with honest limits.

• Real Maud rate-law form and real lognormal priors, in real mM units; only the measurements are synthetic (but a SAM dose-response is exactly the kind of in-vitro data this models). The subunit count n=2 and saturating-substrate assumption are modeling choices stated up front; Vmax, Km, and the other enzymes are fixed rather than co-fit. A complete fit of the whole methionine model against its real metabolomics is Maud's job — this notebook validates the approximation Maud would otherwise force you to pay full HMC for, and shows precisely when that shortcut is safe.
• Our HMC is a hand-rolled, mass-matrix-preconditioned sampler, cross-checked against the exact grid; for higher dimensions you would use NUTS (Maud uses Stan).