Relaxation & Dynamics Analysis¶

This tutorial demonstrates how to use synth-nmr to probe the theoretical internal dynamics of a protein model using the Lipari-Szabo Model-Free formalism and analyze Nuclear Spin Relaxation limits. A new final section compares heuristic S² predictions with S² computed geometrically from the full 10-conformer NMR ensemble.

1. Setup Environment¶

In [ ]:

!pip install -q synth-nmr biotite matplotlib

!pip install -q synth-nmr biotite matplotlib

2. Load Protein Structure¶

In [ ]:

!wget -q https://files.rcsb.org/download/1D3Z.pdb -O ubiquitin.pdb

import biotite.structure.io as strucio
import matplotlib.pyplot as plt
import numpy as np

plt.rcParams.update(
    {"figure.dpi": 110, "font.size": 11, "axes.spines.top": False, "axes.spines.right": False}
)

# 1D3Z is an NMR ensemble with 10 conformers; take model 1 as the reference structure
ensemble_stack = strucio.load_structure("ubiquitin.pdb")
structure = ensemble_stack[0]
print(
    f"Loaded 1D3Z: {len(ensemble_stack)} conformers; using model 1 ({structure.array_length()} atoms)"
)

!wget -q https://files.rcsb.org/download/1D3Z.pdb -O ubiquitin.pdb import biotite.structure.io as strucio import matplotlib.pyplot as plt import numpy as np plt.rcParams.update( {"figure.dpi": 110, "font.size": 11, "axes.spines.top": False, "axes.spines.right": False} ) # 1D3Z is an NMR ensemble with 10 conformers; take model 1 as the reference structure ensemble_stack = strucio.load_structure("ubiquitin.pdb") structure = ensemble_stack[0] print( f"Loaded 1D3Z: {len(ensemble_stack)} conformers; using model 1 ({structure.array_length()} atoms)" )

3. Calculating Relaxation Rates¶

We need to specify the global tumbling time (tau_m_ns) and the strength of the external magnetic field (field_mhz). Higher fields increase the contribution of Chemical Shift Anisotropy (CSA) to the $R_1$ and $R_2$ rates.

In [ ]:

from synth_nmr import calculate_relaxation_rates

rates = calculate_relaxation_rates(
    structure,
    field_mhz=600.0,  # 600 MHz spectrometer
    tau_m_ns=4.1,  # ~4.1 ns tumbling time for Ubiquitin at 300 K
)

print(f"Calculated relaxation parameters for {len(rates)} backbone amide sites.")

from synth_nmr import calculate_relaxation_rates rates = calculate_relaxation_rates( structure, field_mhz=600.0, # 600 MHz spectrometer tau_m_ns=4.1, # ~4.1 ns tumbling time for Ubiquitin at 300 K ) print(f"Calculated relaxation parameters for {len(rates)} backbone amide sites.")

4. Visualizing Internal Dynamics (Order Parameters)¶

The Generalized Order Parameter ($S^2$) ranges from 0 (completely flexible) to 1 (completely rigid). synth-nmr predicts $S^2$ based on the protein's secondary structure and structural elements.

In [ ]:

residues = list(rates.keys())
r1_data = [info["R1"] for info in rates.values()]
r2_data = [info["R2"] for info in rates.values()]
noe_data = [info["NOE"] for info in rates.values()]
s2_data = [info["S2"] for info in rates.values()]

fig, axes = plt.subplots(2, 2, figsize=(13, 8), sharex=True)

for ax, data, label, color in zip(
    axes.flat,
    [r1_data, r2_data, noe_data, s2_data],
    ["R₁ (s⁻¹)", "R₂ (s⁻¹)", "Steady-state NOE", "S² (order parameter)"],
    ["#2980b9", "#e74c3c", "#27ae60", "#8e44ad"],
):
    ax.bar(residues, data, color=color, alpha=0.80, width=0.9)
    ax.set_ylabel(label, fontsize=11)
    ax.grid(axis="y", ls="--", alpha=0.4)

for ax in axes[1]:
    ax.set_xlabel("Residue Number", fontsize=11)

fig.suptitle(
    "Backbone NMR Relaxation Parameters — Ubiquitin (1D3Z model 1)\n600 MHz, \u03c4_m = 4.1 ns",
    fontsize=13,
    fontweight="bold",
)
plt.tight_layout()
plt.show()

residues = list(rates.keys()) r1_data = [info["R1"] for info in rates.values()] r2_data = [info["R2"] for info in rates.values()] noe_data = [info["NOE"] for info in rates.values()] s2_data = [info["S2"] for info in rates.values()] fig, axes = plt.subplots(2, 2, figsize=(13, 8), sharex=True) for ax, data, label, color in zip( axes.flat, [r1_data, r2_data, noe_data, s2_data], ["R₁ (s⁻¹)", "R₂ (s⁻¹)", "Steady-state NOE", "S² (order parameter)"], ["#2980b9", "#e74c3c", "#27ae60", "#8e44ad"], ): ax.bar(residues, data, color=color, alpha=0.80, width=0.9) ax.set_ylabel(label, fontsize=11) ax.grid(axis="y", ls="--", alpha=0.4) for ax in axes[1]: ax.set_xlabel("Residue Number", fontsize=11) fig.suptitle( "Backbone NMR Relaxation Parameters — Ubiquitin (1D3Z model 1)\n600 MHz, \u03c4_m = 4.1 ns", fontsize=13, fontweight="bold", ) plt.tight_layout() plt.show()

Notice how the termini are highly flexible (low $S^2$, low $R_2$), while the core helices and sheets are rigid (high $S^2$, high $R_2$). The next section compares this heuristic S² to one computed directly from the NMR ensemble geometry.

5. Trajectory-Based S² from the Full NMR Ensemble¶

The $S^2$ values above were predicted from a single structure using a secondary-structure heuristic model. We can also measure $S^2$ geometrically from the 10-conformer ensemble directly, using the Lipari-Szabo vector formula:

$$S^2 = |\langle \hat{\mu} \rangle|^2$$

where $\hat{\mu}_i$ is the unit N–H bond vector in conformer $i$. If the N–H bond barely moves across the ensemble, all $\hat{\mu}_i$ point the same way, their mean has large magnitude, and $S^2 \approx 1$ (rigid). If the bond samples all directions, the vectors cancel and $S^2 \approx 0$ (floppy).

Comparing these two sources of $S^2$ reveals where the heuristic agrees with the actual conformational diversity in the NMR ensemble.

In [ ]:

from synth_nmr.trajectory import compute_s2_from_trajectory, load_trajectory

# Wrap all 10 conformers into a TrajectoryEnsemble
all_frames = [ensemble_stack[i] for i in range(len(ensemble_stack))]
traj_ensemble = load_trajectory(all_frames)
print(f"TrajectoryEnsemble: {len(traj_ensemble)} conformers")

# S² = |<μ̂>|² computed geometrically from all 10 NH vector orientations
s2_traj = compute_s2_from_trajectory(traj_ensemble)
print(f"Geometric S² computed for {len(s2_traj)} backbone amide sites.")

from synth_nmr.trajectory import compute_s2_from_trajectory, load_trajectory # Wrap all 10 conformers into a TrajectoryEnsemble all_frames = [ensemble_stack[i] for i in range(len(ensemble_stack))] traj_ensemble = load_trajectory(all_frames) print(f"TrajectoryEnsemble: {len(traj_ensemble)} conformers") # S² = |<μ̂>|² computed geometrically from all 10 NH vector orientations s2_traj = compute_s2_from_trajectory(traj_ensemble) print(f"Geometric S² computed for {len(s2_traj)} backbone amide sites.")

In [ ]:

# Residues present in both S² sources
common_res = sorted(set(rates.keys()) & set(s2_traj.keys()))
s2_rate = [rates[r]["S2"] for r in common_res]  # heuristic prediction
s2_tr = [s2_traj[r] for r in common_res]  # trajectory geometry

fig, axes = plt.subplots(1, 2, figsize=(13, 4.5))

# Left: side-by-side bar chart
ax = axes[0]
x = np.array(common_res)
w = 0.42
ax.bar(
    x - w / 2, s2_rate, w, label="Heuristic (single-structure model)", color="#16a085", alpha=0.82
)
ax.bar(x + w / 2, s2_tr, w, label="Geometric (10-conformer ensemble)", color="#8e44ad", alpha=0.82)
ax.set_ylim(0, 1.15)
ax.set_xlabel("Residue Number", fontsize=12)
ax.set_ylabel("S²", fontsize=12)
ax.set_title("S² Comparison: Heuristic vs. Ensemble Geometry", fontsize=12)
ax.legend(fontsize=10)
ax.grid(axis="y", ls="--", alpha=0.35)

# Right: correlation scatter
ax2 = axes[1]
ax2.scatter(s2_rate, s2_tr, s=50, alpha=0.78, color="#2c3e50", edgecolors="none")
ax2.plot([0, 1.05], [0, 1.05], "r--", lw=1.3, alpha=0.6, label="Perfect agreement")
corr = float(np.corrcoef(s2_rate, s2_tr)[0, 1])
ax2.set_xlabel("S² — heuristic prediction", fontsize=12)
ax2.set_ylabel("S² — 10-conformer ensemble", fontsize=12)
ax2.set_title(f"Pearson r = {corr:.3f}", fontsize=12)
ax2.annotate(
    "High r → heuristic captures\nreal conformational dynamics",
    xy=(0.05, 0.92),
    xycoords="axes fraction",
    fontsize=9,
    color="#555",
    ha="left",
)
ax2.legend(fontsize=10)

plt.tight_layout()
plt.show()

print(f"Pearson correlation (heuristic vs trajectory S²): {corr:.3f}")
print("r close to 1.0 means the heuristic correctly captures the true backbone dynamics.")

# Residues present in both S² sources common_res = sorted(set(rates.keys()) & set(s2_traj.keys())) s2_rate = [rates[r]["S2"] for r in common_res] # heuristic prediction s2_tr = [s2_traj[r] for r in common_res] # trajectory geometry fig, axes = plt.subplots(1, 2, figsize=(13, 4.5)) # Left: side-by-side bar chart ax = axes[0] x = np.array(common_res) w = 0.42 ax.bar( x - w / 2, s2_rate, w, label="Heuristic (single-structure model)", color="#16a085", alpha=0.82 ) ax.bar(x + w / 2, s2_tr, w, label="Geometric (10-conformer ensemble)", color="#8e44ad", alpha=0.82) ax.set_ylim(0, 1.15) ax.set_xlabel("Residue Number", fontsize=12) ax.set_ylabel("S²", fontsize=12) ax.set_title("S² Comparison: Heuristic vs. Ensemble Geometry", fontsize=12) ax.legend(fontsize=10) ax.grid(axis="y", ls="--", alpha=0.35) # Right: correlation scatter ax2 = axes[1] ax2.scatter(s2_rate, s2_tr, s=50, alpha=0.78, color="#2c3e50", edgecolors="none") ax2.plot([0, 1.05], [0, 1.05], "r--", lw=1.3, alpha=0.6, label="Perfect agreement") corr = float(np.corrcoef(s2_rate, s2_tr)[0, 1]) ax2.set_xlabel("S² — heuristic prediction", fontsize=12) ax2.set_ylabel("S² — 10-conformer ensemble", fontsize=12) ax2.set_title(f"Pearson r = {corr:.3f}", fontsize=12) ax2.annotate( "High r → heuristic captures\nreal conformational dynamics", xy=(0.05, 0.92), xycoords="axes fraction", fontsize=9, color="#555", ha="left", ) ax2.legend(fontsize=10) plt.tight_layout() plt.show() print(f"Pearson correlation (heuristic vs trajectory S²): {corr:.3f}") print("r close to 1.0 means the heuristic correctly captures the true backbone dynamics.")