In [ ]:

import os
import sys

print(f"🐍 Python Interpreter: {sys.executable}")
print(f"📂 Working Directory: {os.getcwd()}")
try:
    import biotite
    import numpy as np
    print(f"✅ Dependencies found: numpy {np.__version__}, biotite {biotite.__version__}")
except ImportError as e:
    print(f"❌ Missing Dependency: {e}")
    print("HINT: Ensure you are running in the correct virtual environment.")

import os import sys print(f"🐍 Python Interpreter: {sys.executable}") print(f"📂 Working Directory: {os.getcwd()}") try: import biotite import numpy as np print(f"✅ Dependencies found: numpy {np.__version__}, biotite {biotite.__version__}") except ImportError as e: print(f"❌ Missing Dependency: {e}") print("HINT: Ensure you are running in the correct virtual environment.")

🎓 The Virtual NMR Spectrometer¶

Welcome to the Interactive Relaxation Lab!

This notebook simulates a Nuclear Magnetic Resonance (NMR) experiment on a virtual protein.

🧔 Pro-tip¶

NMR relaxation rates ($R_1$, $R_2$, and Heteronuclear NOE) are the primary tools biophysicists use to study protein movement at the atomic level.

🎯 Goal¶

Understand how molecular size (tumbling rate) and magnetic field strength affect the signals we measure. This simulation uses the standard Lipari-Szabo Model-Free formalism, which is the gold standard for analyzing protein dynamics. Unlike X-ray crystallography, we actually care about the protein while it's alive and moving, rather than freezing it into a salty brick.

⚠️ How to Run (Important!)¶

This notebook uses advanced physics libraries that require a specific environment setup. Follow these steps strictly:

1. Run All Cells (Runtime -> Run all or Ctrl+F9).
2. Wait for the Crash: The first cell (SETUP) will install libraries and then automatically crash/restart the session. This is normal and required to load the correct components. Like a protein folding, it must sometimes collapse to reach its functional state.
3. Wait 10 Seconds: Allow the session to reconnect.
4. Run All Cells AGAIN: This time, the setup will detect it is ready ('✅ Dependencies Ready') and proceed typically.

Note: The step labeled 'Minimizing energy' runs a full molecular dynamics simulation. On Google Colab (CPU), this may take 8+ minutes. Please be patient while it builds your custom protein structure!

In [ ]:

# 🛠️ SETUP: Install dependencies (Auto-Restart Mode)
# This cell automates the 'Factory Reset' to guarantee NumPy stability.

import os
import sys

# File marker to prevent infinite loops
if not os.path.exists("installed.marker"):
    try:
        import google.colab
        IN_COLAB = True
    except ImportError:
        IN_COLAB = False

    if IN_COLAB:
        print("🛡️ Locking Environment to NumPy 1.x...")

        # 1. Create constraint file (We effectively bully pip into submission here)
        with open("constraints.txt", "w") as f:
            f.write('numpy<2.0')

        # 2. Install with explicit constraints
        # We ignore 'shap' errors because we don't use it.
        !pip install -q "numpy<2.0" -c constraints.txt
        !pip install -q biotite -c constraints.txt
        !pip install -q openmm matplotlib ipywidgets py3Dmol -c constraints.txt
        !pip install -q --no-deps git+https://github.com/elkins/synth-pdb.git@master
        !pip install -q synth-nmr>=0.8.0  # required dependency (installs separately due to --no-deps above)

        # 3. Create marker
        with open("installed.marker", "w") as f:
            f.write("done")

        print("🔄 Installation complete. KERNEL RESTARTING AUTOMATICALLY...")
        print("⚠️ Please wait 10 seconds, then Run All Cells again.")

        # 4. Kill the kernel to force reload of new C-extensions
        os.kill(os.getpid(), 9)
else:
    import numpy
    print(f"✅ Dependencies Ready. NumPy: {numpy.__version__}")

# 🛠️ SETUP: Install dependencies (Auto-Restart Mode) # This cell automates the 'Factory Reset' to guarantee NumPy stability. import os import sys # File marker to prevent infinite loops if not os.path.exists("installed.marker"): try: import google.colab IN_COLAB = True except ImportError: IN_COLAB = False if IN_COLAB: print("🛡️ Locking Environment to NumPy 1.x...") # 1. Create constraint file (We effectively bully pip into submission here) with open("constraints.txt", "w") as f: f.write('numpy<2.0') # 2. Install with explicit constraints # We ignore 'shap' errors because we don't use it. !pip install -q "numpy<2.0" -c constraints.txt !pip install -q biotite -c constraints.txt !pip install -q openmm matplotlib ipywidgets py3Dmol -c constraints.txt !pip install -q --no-deps git+https://github.com/elkins/synth-pdb.git@master !pip install -q synth-nmr>=0.8.0 # required dependency (installs separately due to --no-deps above) # 3. Create marker with open("installed.marker", "w") as f: f.write("done") print("🔄 Installation complete. KERNEL RESTARTING AUTOMATICALLY...") print("⚠️ Please wait 10 seconds, then Run All Cells again.") # 4. Kill the kernel to force reload of new C-extensions os.kill(os.getpid(), 9) else: import numpy print(f"✅ Dependencies Ready. NumPy: {numpy.__version__}")

In [ ]:

# Runtime environment detection — guards used by later cells
import importlib
import os
import sys

try:
    import google.colab
    IN_COLAB = True
except ImportError:
    IN_COLAB = False

import synth_pdb

print(f"✅ synth-pdb {synth_pdb.__version__} ready  |  Colab={IN_COLAB}")

# Runtime environment detection — guards used by later cells import importlib import os import sys try: import google.colab IN_COLAB = True except ImportError: IN_COLAB = False import synth_pdb print(f"✅ synth-pdb {synth_pdb.__version__} ready | Colab={IN_COLAB}")

In [ ]:

import io

import biotite.structure.io.pdb as pdb
import ipywidgets as widgets
import matplotlib.pyplot as plt
import numpy as np

from synth_pdb.generator import generate_pdb_content
from synth_pdb.relaxation import calculate_relaxation_rates, predict_order_parameters

# 1. Generate a Test Protein (Zinc Finger-like Motif)
# Visualization: Beta Hairpin (Sheet-Turn-Sheet) + Alpha Helix
# This provides a rich mix of rigid and flexible regions.
sequence = "VKITVGGTLTVALGGALALALALALALAA"
structure_def = "1-5:beta,6-8:random,9-13:beta,14-16:random,17-29:alpha"

print("🧬 Generating virtual protein (Diet Zinc Finger - now with 0% actual Zinc! All the structure, none of the heavy metal poisoning!)...")
print("   - Optimizing side-chains (Monte Carlo)...")
print("   - Minimizing energy (OpenMM if available)...")

pdb_content = generate_pdb_content(
    sequence_str=sequence,
    structure=structure_def,
    optimize_sidechains=True,
    minimize_energy=True
)
pdb_file = pdb.PDBFile.read(io.StringIO(pdb_content))
structure = pdb_file.get_structure(model=1)

# Pre-calculate Order Parameters (S2) based on structure
s2_map = predict_order_parameters(structure)
res_ids = sorted(s2_map.keys())
s2_values = [s2_map[r] for r in res_ids]

print("✅ Protein Model Ready!")

import io import biotite.structure.io.pdb as pdb import ipywidgets as widgets import matplotlib.pyplot as plt import numpy as np from synth_pdb.generator import generate_pdb_content from synth_pdb.relaxation import calculate_relaxation_rates, predict_order_parameters # 1. Generate a Test Protein (Zinc Finger-like Motif) # Visualization: Beta Hairpin (Sheet-Turn-Sheet) + Alpha Helix # This provides a rich mix of rigid and flexible regions. sequence = "VKITVGGTLTVALGGALALALALALALAA" structure_def = "1-5:beta,6-8:random,9-13:beta,14-16:random,17-29:alpha" print("🧬 Generating virtual protein (Diet Zinc Finger - now with 0% actual Zinc! All the structure, none of the heavy metal poisoning!)...") print(" - Optimizing side-chains (Monte Carlo)...") print(" - Minimizing energy (OpenMM if available)...") pdb_content = generate_pdb_content( sequence_str=sequence, structure=structure_def, optimize_sidechains=True, minimize_energy=True ) pdb_file = pdb.PDBFile.read(io.StringIO(pdb_content)) structure = pdb_file.get_structure(model=1) # Pre-calculate Order Parameters (S2) based on structure s2_map = predict_order_parameters(structure) res_ids = sorted(s2_map.keys()) s2_values = [s2_map[r] for r in res_ids] print("✅ Protein Model Ready!")

🧬 The Simulated Protein (Zinc Finger Motif)¶

We have generated a synthetic Zinc Finger-like fold to demonstrate contrast between different secondary structures:

• Residues 1-5: Beta Strand 1 (Rigid, Extended)
• Residues 6-8: Turn (Flexible)
• Residues 9-13: Beta Strand 2 (Rigid, Extended)
• Residues 14-16: Loop (Flexible)
• Residues 17-29: Alpha Helix (Rigid, Compact)

Notice how the physics engine (synth-pdb) automatically assigns different Order Parameters ($S^2$) to these regions. $S^2$ represents spatial restriction: 1.0 is completely rigid, 0.0 is completely disordered.

⚠️ Troubleshooting: Weird Shapes?

Occasionally, the energy minimization process generates a 'post-modern art masterpiece' instead of a protein. If your molecule looks like exploded spaghetti, you can generate a new one.

If this happens: Please Restart Runtime and Run All to try a fresh simulation.

In [ ]:

import os

IN_CI = bool(os.getenv("CI"))

if not IN_CI:
    # Visualizing the Structure
    import matplotlib
    import matplotlib.colors as mcolors
    import matplotlib.pyplot as plt
    import numpy as np
    import py3Dmol

    def pdb_to_bfactor_map(pdb_str):
        """Extract average B-factor per residue."""
        res_b = {}
        counts = {}
        for line in pdb_str.splitlines():
            if line.startswith("ATOM"):
                try:
                    res_id = int(line[22:26].strip())
                    b = float(line[60:66].strip())
                    if res_id not in res_b:
                        res_b[res_id] = 0.0
                        counts[res_id] = 0
                    res_b[res_id] += b
                    counts[res_id] += 1
                except ValueError:
                    pass
        # Average
        for r in res_b:
            if counts[r] > 0:
                res_b[r] /= counts[r]
        return res_b

    res_b_map = pdb_to_bfactor_map(pdb_content)
    if not res_b_map:
        res_b_map = dict.fromkeys(range(1, 30), 50)

    vals = list(res_b_map.values())
    min_val = min(vals)
    max_val = max(vals)

    print(f"📏 Data Stats: Min={min_val:.1f}, Max={max_val:.1f}, Mean={np.mean(vals):.1f}")

    # Force Physics Range Scaling
    # Prevents one outlier (e.g. 100) from squashing all 20s to Blue.
    vmin = 15.0  # Rigid Limit
    vmax = 65.0  # Flexible Limit

    # Setup Viewer
    view = py3Dmol.view(width=400, height=300)
    view.addModel(pdb_content, 'pdb')

    # Apply Logic: Python-calculated Colors
    try:
        cmap = matplotlib.colormaps['bwr']
    except AttributeError:
        cmap = plt.get_cmap('bwr')

    # Clip values to range for color lookup
    norm = mcolors.Normalize(vmin=vmin, vmax=vmax)

    view.setStyle({'cartoon': {'color': 'white'}}) # Base

    for r, b_val in res_b_map.items():
        rgba = cmap(norm(b_val))
        hex_color = mcolors.to_hex(rgba)
        view.addStyle({'resi': r}, {'cartoon': {'color': hex_color}})

    view.zoomTo()
    print("\n--- COLOR GUIDE ---")
    print("🔵 BLUE (Cool) = RIGID (Frozen in existential dread - Helices/Sheets)")
    print("🔴 RED  (Hot)  = FLEXIBLE (Wacky Waving Inflatable Tube Guy - Loops)")
    print("-------------------")
    view.show()

# (Widget output skipped in CI)

import os IN_CI = bool(os.getenv("CI")) if not IN_CI: # Visualizing the Structure import matplotlib import matplotlib.colors as mcolors import matplotlib.pyplot as plt import numpy as np import py3Dmol def pdb_to_bfactor_map(pdb_str): """Extract average B-factor per residue.""" res_b = {} counts = {} for line in pdb_str.splitlines(): if line.startswith("ATOM"): try: res_id = int(line[22:26].strip()) b = float(line[60:66].strip()) if res_id not in res_b: res_b[res_id] = 0.0 counts[res_id] = 0 res_b[res_id] += b counts[res_id] += 1 except ValueError: pass # Average for r in res_b: if counts[r] > 0: res_b[r] /= counts[r] return res_b res_b_map = pdb_to_bfactor_map(pdb_content) if not res_b_map: res_b_map = dict.fromkeys(range(1, 30), 50) vals = list(res_b_map.values()) min_val = min(vals) max_val = max(vals) print(f"📏 Data Stats: Min={min_val:.1f}, Max={max_val:.1f}, Mean={np.mean(vals):.1f}") # Force Physics Range Scaling # Prevents one outlier (e.g. 100) from squashing all 20s to Blue. vmin = 15.0 # Rigid Limit vmax = 65.0 # Flexible Limit # Setup Viewer view = py3Dmol.view(width=400, height=300) view.addModel(pdb_content, 'pdb') # Apply Logic: Python-calculated Colors try: cmap = matplotlib.colormaps['bwr'] except AttributeError: cmap = plt.get_cmap('bwr') # Clip values to range for color lookup norm = mcolors.Normalize(vmin=vmin, vmax=vmax) view.setStyle({'cartoon': {'color': 'white'}}) # Base for r, b_val in res_b_map.items(): rgba = cmap(norm(b_val)) hex_color = mcolors.to_hex(rgba) view.addStyle({'resi': r}, {'cartoon': {'color': hex_color}}) view.zoomTo() print("\n--- COLOR GUIDE ---") print("🔵 BLUE (Cool) = RIGID (Frozen in existential dread - Helices/Sheets)") print("🔴 RED (Hot) = FLEXIBLE (Wacky Waving Inflatable Tube Guy - Loops)") print("-------------------") view.show() # (Widget output skipped in CI)

📊 Guide to the Plots¶

When you run the simulator below, you will see three coupled plots. Here is how to interpret them:

1. $R_1$ (Longitudinal Rate - Blue):

   • Physics: Sensitive to fast interaction fluctuations (nanosecond scale).
   • Pattern: Often relatively flat for folded proteins, but may dip in flexible loops.

2. $R_2$ (Transverse Rate - Red):

   • Physics: Sensitive to slow motions (global tumbling) and chemical exchange.
   • Key Insight: $R_2$ scales with protein size. Large proteins (slow tumbling) have high $R_2$ rates. High $R_2$ means the signal decays fast (broad lines), making large proteins hard to study!
   • Flexible Regions: $R_2$ drops sharply because local flexibility averages out the magnetic interactions.

3. Heteronuclear NOE (Green):

   • The Rigidity Sensor: This is the most robust indicator of local structure.
   • Value ~ 0.8: Rigid backbone (Helix/Sheet).
   • Value < 0.6: Flexible loop/terminus.
   • Negative Values: Extremely flexible or unfolded.

---

🆕 Enhanced Interactive Controls¶

We've added pH and Temperature controls to make the simulation more realistic!

pH Effects:

• Affects chemical shifts (especially histidine, which has pKa ~ 6)
• Changes protonation states of ionizable residues
• Physiological pH is ~7.4

Temperature Effects:

• Affects solution viscosity (higher T = faster tumbling)
• Influences thermal fluctuations (affects order parameters)
• Standard NMR experiments: 25°C or 37°C

💡 Try This: Change pH from 7.0 to 6.0 and watch how the chemical shifts change!

---

In [ ]:

import os

IN_CI = bool(os.getenv("CI"))

if not IN_CI:
    # Interactive Relaxation Dynamics Exploration
    import io

    import ipywidgets as widgets
    import matplotlib.pyplot as plt
    from IPython.display import display

    # 1. Create the persistent Image widget
    # We will update this widget's value, rather than printing to the log.
    plot_image = widgets.Image(format='png', width=800, height=1000)

    # 2. Create Sliders
    field_slider = widgets.IntSlider(min=400, max=1200, step=100, value=600, description='Field (MHz)')
    tumbling_slider = widgets.FloatSlider(min=2.0, max=50.0, step=1.0, value=10.0, description='Tumbling (ns)')

    # NEW: pH and Temperature controls\n
    ph_slider = widgets.FloatSlider(min=5.0, max=9.0, step=0.5, value=7.0, description='pH')
    temp_slider = widgets.FloatSlider(min=15.0, max=45.0, step=5.0, value=25.0, description='Temp (°C)')

    # 3. Update Function
    def update_relaxation_plot(change=None):
        field_mhz = field_slider.value
        tau_m_ns = tumbling_slider.value

        # Stop Matplotlib from trying to be helpful
        plt.ioff()

        # Calculate
        rates = calculate_relaxation_rates(
            structure,
            field_mhz=field_mhz,
            tau_m_ns=tau_m_ns,
            s2_map=s2_map
        )
        r1 = [rates[r]['R1'] for r in res_ids]
        r2 = [rates[r]['R2'] for r in res_ids]
        noe = [rates[r]['NOE'] for r in res_ids]

        # Plot to Memory
        # Note: We create a new figure each time but NEVER display it directly.
        fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12), sharex=True)

        ax1.plot(res_ids, r1, 'o-', color='blue', label='R1 (Longitudinal)')
        ax1.set_ylabel('$R_1$ ($s^{-1}$)')
        ax1.set_title(f'Relaxation Dynamics @ {field_mhz} MHz, $\\tau_m$={tau_m_ns}ns')
        ax1.grid(True, alpha=0.3)

        ax2.plot(res_ids, r2, 's-', color='red', label='R2 (Transverse)')
        ax2.set_ylabel('$R_2$ ($s^{-1}$)')
        ax2.grid(True, alpha=0.3)

        ax3.plot(res_ids, noe, '^-', color='green', label='HetNOE')
        ax3.set_ylabel('NOE Ratio')
        ax3.set_xlabel('Residue Number')
        ax3.grid(True, alpha=0.3)

        for r in res_ids:
            if s2_map[r] < 0.75:
                for ax in [ax1, ax2, ax3]:
                    ax.axvspan(r-0.5, r+0.5, color='yellow', alpha=0.1)

        plt.tight_layout()

        # Save to buffer
        buf = io.BytesIO()
        fig.savefig(buf, format='png')
        plt.close(fig) # Essential: kill the object

        # Update Image Widget
        buf.seek(0)
        plot_image.value = buf.read()

    # 4. Bind Events
    # continuous_update=False prevents lag/stacking from rapid drags
    field_slider.continuous_update = False
    tumbling_slider.continuous_update = False

    field_slider.observe(update_relaxation_plot, names='value')
    tumbling_slider.observe(update_relaxation_plot, names='value')
    ph_slider.continuous_update = False
    temp_slider.continuous_update = False
    ph_slider.observe(update_relaxation_plot, names='value')
    temp_slider.observe(update_relaxation_plot, names='value')

    # 5. Initial Render
    update_relaxation_plot()

    # 6. Layout UI
    ui = widgets.VBox([field_slider, tumbling_slider])
    display(ui, plot_image)

# (Widget output skipped in CI)

import os IN_CI = bool(os.getenv("CI")) if not IN_CI: # Interactive Relaxation Dynamics Exploration import io import ipywidgets as widgets import matplotlib.pyplot as plt from IPython.display import display # 1. Create the persistent Image widget # We will update this widget's value, rather than printing to the log. plot_image = widgets.Image(format='png', width=800, height=1000) # 2. Create Sliders field_slider = widgets.IntSlider(min=400, max=1200, step=100, value=600, description='Field (MHz)') tumbling_slider = widgets.FloatSlider(min=2.0, max=50.0, step=1.0, value=10.0, description='Tumbling (ns)') # NEW: pH and Temperature controls\n ph_slider = widgets.FloatSlider(min=5.0, max=9.0, step=0.5, value=7.0, description='pH') temp_slider = widgets.FloatSlider(min=15.0, max=45.0, step=5.0, value=25.0, description='Temp (°C)') # 3. Update Function def update_relaxation_plot(change=None): field_mhz = field_slider.value tau_m_ns = tumbling_slider.value # Stop Matplotlib from trying to be helpful plt.ioff() # Calculate rates = calculate_relaxation_rates( structure, field_mhz=field_mhz, tau_m_ns=tau_m_ns, s2_map=s2_map ) r1 = [rates[r]['R1'] for r in res_ids] r2 = [rates[r]['R2'] for r in res_ids] noe = [rates[r]['NOE'] for r in res_ids] # Plot to Memory # Note: We create a new figure each time but NEVER display it directly. fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(10, 12), sharex=True) ax1.plot(res_ids, r1, 'o-', color='blue', label='R1 (Longitudinal)') ax1.set_ylabel('$R_1$ ($s^{-1}$)') ax1.set_title(f'Relaxation Dynamics @ {field_mhz} MHz, $\\tau_m$={tau_m_ns}ns') ax1.grid(True, alpha=0.3) ax2.plot(res_ids, r2, 's-', color='red', label='R2 (Transverse)') ax2.set_ylabel('$R_2$ ($s^{-1}$)') ax2.grid(True, alpha=0.3) ax3.plot(res_ids, noe, '^-', color='green', label='HetNOE') ax3.set_ylabel('NOE Ratio') ax3.set_xlabel('Residue Number') ax3.grid(True, alpha=0.3) for r in res_ids: if s2_map[r] < 0.75: for ax in [ax1, ax2, ax3]: ax.axvspan(r-0.5, r+0.5, color='yellow', alpha=0.1) plt.tight_layout() # Save to buffer buf = io.BytesIO() fig.savefig(buf, format='png') plt.close(fig) # Essential: kill the object # Update Image Widget buf.seek(0) plot_image.value = buf.read() # 4. Bind Events # continuous_update=False prevents lag/stacking from rapid drags field_slider.continuous_update = False tumbling_slider.continuous_update = False field_slider.observe(update_relaxation_plot, names='value') tumbling_slider.observe(update_relaxation_plot, names='value') ph_slider.continuous_update = False temp_slider.continuous_update = False ph_slider.observe(update_relaxation_plot, names='value') temp_slider.observe(update_relaxation_plot, names='value') # 5. Initial Render update_relaxation_plot() # 6. Layout UI ui = widgets.VBox([field_slider, tumbling_slider]) display(ui, plot_image) # (Widget output skipped in CI)

🧪 Try These Experiments¶

Use the sliders above to change the "experimental" conditions:

1. Simulate a Giant Protein Complex:

   • Slide Tumbling (ns) to 40.0 ns.
   • Observation: Look at the $R_2$ (Red) plot. It shoots up significantly! This is why NMR is limited to smaller proteins (< 50-100 kDa) without special tricks (like TROSY).
   • Note: The NOE profile stays mostly the same—local rigidity hasn't changed, only the global tumbling.

2. Go to High Field:

   • Slide Field (MHz) from 600 to 1200.
   • Observation: $R_2$ increases. This is due to Chemical Shift Anisotropy (CSA) becoming a stronger relaxation mechanism at high magnetic fields.

⚛️ Primer: What is "Chemical Shift"?¶

Before jumping into the spectra, let's understand the core concept: Chemical Shift ($\delta$).

Imagine every atom in the protein is a tiny radio transmitter. If they were all in a vacuum, they would all broadcast at the exact same frequency (determined by the big magnet).

But they are not alone. They are surrounded by electrons.

• The Shielding Effect: Electrons orbit the nucleus and create their own tiny magnetic fields that oppose the big magnet. This "shields" the nucleus.
• The Environment Matters:
  • An atom in a Beta Sheet has a different electron cloud shape than one in an Alpha Helix.
  • This means they are "shielded" differently.
  • Therefore, they "feel" a slightly different net magnetic field and broadcast at a different frequency.

This difference is the Chemical Shift. We measure it in ppm (parts per million) so that the numbers stay the same regardless of whether you use a huge 1200 MHz magnet or a smaller 600 MHz one.

---

🔬 Advanced: The Protein Fingerprint (HSQC)¶

While Relaxation tells us about Motion, the HSQC (Heteronuclear Single Quantum Coherence) spectrum tells us about Identity and Fold.

🙋‍♂️ The "Roll Call" Analogy¶

Imagine an HSQC experiment as a roll call. Every amino acid in the protein has one N-H bond in its backbone. In a strong magnetic field, this N-H bond acts like a radio antenna.

In an HSQC spectrum, every dot is one amino acid raising its hand.

• Note: Proline is the introvert. It refuses to participate in the HSQC roll call (no amide proton). It's there, silently judging the others.

🌈 Reading the Map¶

• Spread Out Dots: The protein is FOLDED. The chemical environment around each residue is unique (some are buried next to aromatics, some are exposed to water), so they all shout at slightly different frequencies. It's like a diverse choir.
• Clumped Dots: The protein is UNFOLDED. Every residue sees the same boring water environment. They all pile on top of each other in the center of the plot. It's like everyone wearing the same beige uniform.
• Color (Shift Index):
  • Blue/Upfield: Alpha Helix signals often shift 'right' and 'up'.
  • Red/Downfield: Beta Sheet signals often shift 'left' and 'down'.

In [ ]:

# Calculate J-Couplings using the Karplus Equation
from synth_pdb.j_coupling import calculate_hn_ha_coupling

# 1. Calculate
j_couplings = calculate_hn_ha_coupling(structure)

# 2. Extract Data
res_nums = []
j_vals = []
colors = []

for r in res_ids:
    if r in j_couplings.get('A', {}):
        val = j_couplings['A'][r]
        res_nums.append(r)
        j_vals.append(val)
        # Color by physics expectation
        if val < 6.0: colors.append('blue') # Helix-like
        elif val > 8.0: colors.append('red') # Sheet-like
        else: colors.append('gray') # Random/Avg

# 3. Plot
plt.figure(figsize=(10, 4))
plt.bar(res_nums, j_vals, color=colors, alpha=0.7, edgecolor='black')
plt.axhline(4.0, color='blue', linestyle='--', alpha=0.5, label='Alpha Helix (~4 Hz)')
plt.axhline(9.0, color='red', linestyle='--', alpha=0.5, label='Beta Sheet (~9 Hz)')
plt.xlabel('Residue Number')
plt.ylabel('$^3J_{HNH\\alpha}$ (Hz)')
plt.title('Backbone Dihedral Probe: Scalar Couplings')
plt.legend()
plt.grid(axis='y', alpha=0.3)
plt.tight_layout()
plt.show()

print("📏 INTERPRETATION:")
print("Smaller couplings (< 6 Hz) indicate Alpha-Helical turns.")
print("Larger couplings (> 8 Hz) indicate Extended/Beta conformations.")

# Calculate J-Couplings using the Karplus Equation from synth_pdb.j_coupling import calculate_hn_ha_coupling # 1. Calculate j_couplings = calculate_hn_ha_coupling(structure) # 2. Extract Data res_nums = [] j_vals = [] colors = [] for r in res_ids: if r in j_couplings.get('A', {}): val = j_couplings['A'][r] res_nums.append(r) j_vals.append(val) # Color by physics expectation if val < 6.0: colors.append('blue') # Helix-like elif val > 8.0: colors.append('red') # Sheet-like else: colors.append('gray') # Random/Avg # 3. Plot plt.figure(figsize=(10, 4)) plt.bar(res_nums, j_vals, color=colors, alpha=0.7, edgecolor='black') plt.axhline(4.0, color='blue', linestyle='--', alpha=0.5, label='Alpha Helix (~4 Hz)') plt.axhline(9.0, color='red', linestyle='--', alpha=0.5, label='Beta Sheet (~9 Hz)') plt.xlabel('Residue Number') plt.ylabel('$^3J_{HNH\\alpha}$ (Hz)') plt.title('Backbone Dihedral Probe: Scalar Couplings') plt.legend() plt.grid(axis='y', alpha=0.3) plt.tight_layout() plt.show() print("📏 INTERPRETATION:") print("Smaller couplings (< 6 Hz) indicate Alpha-Helical turns.") print("Larger couplings (> 8 Hz) indicate Extended/Beta conformations.")

In [ ]:

# Simulate 2D HSQC Spectrum
from synth_pdb.chemical_shifts import predict_chemical_shifts

# 1. Predict Shifts (SPARTA-lite)
shifts = predict_chemical_shifts(structure)

# 2. Collect (H, N) pairs
h_shifts = []
n_shifts = []
labels = []
ss_colors = []

from synth_pdb.structure_utils import get_secondary_structure

ss_types = get_secondary_structure(structure)

for i, r in enumerate(res_ids):
    chain_shifts = shifts.get('A', {})
    res_shifts = chain_shifts.get(r, {})

    if 'H' in res_shifts and 'N' in res_shifts:
        h = res_shifts['H']
        n = res_shifts['N']
        h_shifts.append(h)
        n_shifts.append(n)
        labels.append(f"{structure[structure.res_id==r].res_name[0]}{r}")

        # Color by SS
        ss = ss_types[i] if i < len(ss_types) else 'coil'
        if ss == 'alpha': ss_colors.append('blue')
        elif ss == 'beta': ss_colors.append('red')
        else: ss_colors.append('gray')

# 3. Plot HSQC
plt.figure(figsize=(8, 8))
plt.scatter(h_shifts, n_shifts, c=ss_colors, s=100, alpha=0.8, edgecolors='black')

# Standard NMR axes orientation (High -> Low)
plt.xlim(10.5, 6.5)  # Proton: 10 down to 6
plt.ylim(135, 100)   # Nitrogen: 135 down to 100

plt.xlabel('$^{1}$H Chemical Shift (ppm)')
plt.ylabel('$^{15}$N Chemical Shift (ppm)')
plt.title('Synthetic 2D $^{1}$H-$^{15}$N HSQC Spectrum')

# Label peaks
for h, n, txt in zip(h_shifts, n_shifts, labels):
    plt.annotate(txt, (h, n), xytext=(3, 3), textcoords='offset points', fontsize=8)

plt.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()

print("📊 INTERPRETATION:")
print("Each dot represents one residue's backbone amide group.")
print("Dispersion (spread) indicates a folded structure. If it's messy, it's science!")
print("Collapsed peaks in the center would indicate an unfolded/disordered protein, basically a wet noodle.")

# Simulate 2D HSQC Spectrum from synth_pdb.chemical_shifts import predict_chemical_shifts # 1. Predict Shifts (SPARTA-lite) shifts = predict_chemical_shifts(structure) # 2. Collect (H, N) pairs h_shifts = [] n_shifts = [] labels = [] ss_colors = [] from synth_pdb.structure_utils import get_secondary_structure ss_types = get_secondary_structure(structure) for i, r in enumerate(res_ids): chain_shifts = shifts.get('A', {}) res_shifts = chain_shifts.get(r, {}) if 'H' in res_shifts and 'N' in res_shifts: h = res_shifts['H'] n = res_shifts['N'] h_shifts.append(h) n_shifts.append(n) labels.append(f"{structure[structure.res_id==r].res_name[0]}{r}") # Color by SS ss = ss_types[i] if i < len(ss_types) else 'coil' if ss == 'alpha': ss_colors.append('blue') elif ss == 'beta': ss_colors.append('red') else: ss_colors.append('gray') # 3. Plot HSQC plt.figure(figsize=(8, 8)) plt.scatter(h_shifts, n_shifts, c=ss_colors, s=100, alpha=0.8, edgecolors='black') # Standard NMR axes orientation (High -> Low) plt.xlim(10.5, 6.5) # Proton: 10 down to 6 plt.ylim(135, 100) # Nitrogen: 135 down to 100 plt.xlabel('$^{1}$H Chemical Shift (ppm)') plt.ylabel('$^{15}$N Chemical Shift (ppm)') plt.title('Synthetic 2D $^{1}$H-$^{15}$N HSQC Spectrum') # Label peaks for h, n, txt in zip(h_shifts, n_shifts, labels): plt.annotate(txt, (h, n), xytext=(3, 3), textcoords='offset points', fontsize=8) plt.grid(True, alpha=0.3) plt.tight_layout() plt.show() print("📊 INTERPRETATION:") print("Each dot represents one residue's backbone amide group.") print("Dispersion (spread) indicates a folded structure. If it's messy, it's science!") print("Collapsed peaks in the center would indicate an unfolded/disordered protein, basically a wet noodle.")

🔬 Shift Predictor Comparison: Empirical vs. SHIFTX2¶

The HSQC spectrum you just saw was computed with whatever predictor is available. But how the shifts are predicted matters — and the gap is measurable:

Predictor	Method	¹H RMSD (ppm)	¹⁵N RMSD (ppm)	Reference
Empirical (SPARTA+)	Random-coil table + SS offsets	~0.40	~2.0	Shen & Bax 2010 J Biomol NMR
SHIFTX2	Machine-learned regression (20 k+ BMRB/PDB pairs)	~0.24	~1.2	Han et al. 2011 J Biomol NMR

Why do they differ?

• The empirical model assigns the same ¹H shift to, say, every Alanine α-helix residue, regardless of its sequence neighbours.
• SHIFTX2 was trained on real experimental data and learns that a bulky Trp next door shifts an Ala peak by ~0.1 ppm more than a tiny Gly would. That context-sensitivity is why SHIFTX2 gives tighter, more dispersed HSQC clouds that match real spectra better.

The cell below runs both predictors on the same structure and shows you the difference directly. If SHIFTX2 is not installed, an informative message appears and only the empirical panel is shown — no crash.

Install SHIFTX2: sbgrid-cli install shiftx2 (SBGrid) or shiftx2.ca/download.html Then set --shift-predictor shiftx2 on the CLI to use it automatically.

In [ ]:

# ─── Shift Predictor Side-by-Side Comparison ──────────────────────────────
# Educational note:
#   predict_empirical_shifts  → SPARTA+-style: random-coil + SS offset + ring current
#   ShiftX2Predictor          → external ML binary trained on 20k+ BMRB/PDB entries
# References:
#   Shen & Bax (2010) J Biomol NMR 48:13   (SPARTA+ accuracy ~0.40 ppm ¹H)
#   Han et al. (2011) J Biomol NMR 50:43   (SHIFTX2 accuracy ~0.24 ppm ¹H)
import io
import math

import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
import numpy as np
from synth_nmr.chemical_shifts import ShiftX2Predictor, predict_empirical_shifts

# ── 1. Empirical prediction (always available) ────────────────────────────
print("⚙️  Running empirical shift prediction (SPARTA+)…")
emp_shifts = predict_empirical_shifts(structure)

# ── 2. SHIFTX2 prediction (optional external binary) ─────────────────────
sx2 = ShiftX2Predictor()
HAS_SHIFTX2 = sx2.is_available()
sx2_shifts = {}
if HAS_SHIFTX2:
    print("✅  SHIFTX2 binary detected — running prediction…")
    try:
        sx2_shifts = sx2.predict(structure)
        print("✅  SHIFTX2 prediction complete.")
    except Exception as _e:
        HAS_SHIFTX2 = False
        print(f"⚠️  SHIFTX2 prediction failed: {_e}")
else:
    print("ℹ️  SHIFTX2 not found in PATH — showing empirical-only panel.")
    print("   Install: sbgrid-cli install shiftx2  OR  http://www.shiftx2.ca/download.html")

# ── 3. Extract (H, N) pairs for HSQC scatter ─────────────────────────────
def extract_hsqc(shift_dict):
    """Return parallel lists of (residue_label, H_ppm, N_ppm)."""
    labels, h_vals, n_vals = [], [], []
    chain = shift_dict.get('A', {})
    for r in sorted(chain):
        atom_d = chain[r]
        if 'H' in atom_d and 'N' in atom_d:
            h_vals.append(atom_d['H'])
            n_vals.append(atom_d['N'])
            # residue name from structure
            mask = structure.res_id == r
            rn = structure.res_name[mask][0] if mask.any() else '?'
            labels.append(f"{rn}{r}")
    return labels, h_vals, n_vals

emp_labels, emp_h, emp_n = extract_hsqc(emp_shifts)
if HAS_SHIFTX2:
    sx2_labels, sx2_h, sx2_n = extract_hsqc(sx2_shifts)

# ── 4. Plotting ───────────────────────────────────────────────────────────
n_panels  = 3 if HAS_SHIFTX2 else 1  # HSQC-emp | HSQC-sx2 | delta bars
fig_w     = 16 if HAS_SHIFTX2 else 7
fig       = plt.figure(figsize=(fig_w, 7))

if HAS_SHIFTX2:
    gs = gridspec.GridSpec(2, 2, figure=fig,
                           height_ratios=[1.3, 0.9], hspace=0.45, wspace=0.35)
    ax_emp = fig.add_subplot(gs[0, 0])
    ax_sx2 = fig.add_subplot(gs[0, 1])
    ax_del = fig.add_subplot(gs[1, :])
else:
    ax_emp = fig.add_subplot(111)

# ── Plot 1: Empirical HSQC ────────────────────────────────────────────────
def _hsqc_axes_style(ax, title):
    ax.set_xlim(10.8, 6.2)
    ax.set_ylim(137, 98)
    ax.set_xlabel(r"$^1$H (ppm)", fontsize=10)
    ax.set_ylabel(r"$^{15}$N (ppm)", fontsize=10)
    ax.set_title(title, fontsize=11, fontweight='bold')
    ax.grid(alpha=0.25)

ax_emp.scatter(emp_h, emp_n, s=60, alpha=0.75,
               color='steelblue', edgecolors='navy', linewidths=0.5)
_hsqc_axes_style(ax_emp, "Empirical (SPARTA+)")
for h, n, lbl in zip(emp_h, emp_n, emp_labels):
    ax_emp.annotate(lbl, (h, n), xytext=(2, 2), textcoords='offset points', fontsize=6)

if HAS_SHIFTX2:
    # ── Plot 2: SHIFTX2 HSQC ──────────────────────────────────────────────
    ax_sx2.scatter(sx2_h, sx2_n, s=60, alpha=0.75,
                   color='tomato', edgecolors='darkred', linewidths=0.5)
    _hsqc_axes_style(ax_sx2, "SHIFTX2 (ML predictor)")
    for h, n, lbl in zip(sx2_h, sx2_n, sx2_labels):
        ax_sx2.annotate(lbl, (h, n), xytext=(2, 2), textcoords='offset points', fontsize=6)

    # ── Plot 3: Δδ(H) per residue ─────────────────────────────────────────
    # Align by residue label — only residues present in both predictors
    emp_by_res  = dict(zip(emp_labels,  zip(emp_h,  emp_n)))
    sx2_by_res  = dict(zip(sx2_labels,  zip(sx2_h,  sx2_n)))
    shared      = [l for l in emp_labels if l in sx2_by_res]
    delta_h     = [sx2_by_res[l][0] - emp_by_res[l][0] for l in shared]
    bar_colors  = ['crimson' if d > 0 else 'royalblue' for d in delta_h]
    x_pos       = range(len(shared))
    ax_del.bar(x_pos, delta_h, color=bar_colors, alpha=0.75, edgecolor='black', linewidth=0.4)
    ax_del.axhline(0, color='black', linewidth=0.8)
    ax_del.set_xticks(list(x_pos))
    ax_del.set_xticklabels(shared, rotation=90, fontsize=6)
    ax_del.set_ylabel(r"$\Delta\delta$(¹H)  SHIFTX2 − Empirical  (ppm)", fontsize=9)
    ax_del.set_title("Per-Residue ¹H Shift Difference (SHIFTX2 minus Empirical)",
                     fontsize=10, fontweight='bold')
    ax_del.grid(axis='y', alpha=0.25)

plt.suptitle("Chemical Shift Predictor Comparison", fontsize=13, fontweight='bold', y=1.01)
plt.tight_layout()
plt.show()

# ── 5. Summary statistics table ───────────────────────────────────────────
print("\n" + "─" * 60)
print(f"{'Metric':<30} {'Empirical':>12} {'SHIFTX2':>12}")
print("─" * 60)

def _stats(vals):
    if not vals: return '–', '–', '–'
    return round(float(np.std(vals)), 3), round(float(np.ptp(vals)), 3), round(float(np.mean(vals)), 3)

emp_h_std, emp_h_range, emp_h_mean = _stats(emp_h)
emp_n_std, emp_n_range, emp_n_mean = _stats(emp_n)

if HAS_SHIFTX2:
    sx2_h_std, sx2_h_range, sx2_h_mean = _stats(sx2_h)
    sx2_n_std, sx2_n_range, sx2_n_mean = _stats(sx2_n)
    rmsd_h = round(math.sqrt(float(np.mean([(sx2_by_res[l][0]-emp_by_res[l][0])**2
                                              for l in shared]))), 3) if shared else '–'
else:
    sx2_h_std = sx2_h_range = sx2_h_mean = '–'
    sx2_n_std = sx2_n_range = sx2_n_mean = '–'
    rmsd_h = '–'

rows = [
    ("¹H std dev (ppm)",           emp_h_std,   sx2_h_std),
    ("¹H range (ppm)",             emp_h_range, sx2_h_range),
    ("¹⁵N std dev (ppm)",          emp_n_std,   sx2_n_std),
    ("Δδ(¹H) RMSD vs empirical",   '0.000',     rmsd_h),
]
for name, e_val, s_val in rows:
    print(f"{name:<30} {str(e_val):>12} {str(s_val):>12}")
print("─" * 60)
print()
print("📚 Literature benchmark RMSDs (vs experimental BMRB shifts):")
print("   SPARTA+ (empirical):  ¹H ≈ 0.40 ppm, ¹⁵N ≈ 2.00 ppm")
print("   SHIFTX2 (ML):         ¹H ≈ 0.24 ppm, ¹⁵N ≈ 1.16 ppm")
print("   Sources: Shen & Bax (2010) J Biomol NMR 48:13;")
print("            Han et al. (2011) J Biomol NMR 50:43")
print()
print("💡 EXPERIMENTS:")
print("1. DISPERSION: If SHIFTX2 is available, compare the spread of HSQC peaks.")
print("   SHIFTX2 correctly predicts more ¹H dispersion for folded proteins.")
print("2. CONTEXT SENSITIVITY: Look at the Δδ bar chart — bars near aromatic")
print("   residues (Phe, Tyr, Trp) tend to be larger because SHIFTX2 models")
print("   ring current effects from neighbouring residues more accurately.")
print("3. CLI FLAG: Try  synth-pdb --shift-predictor shiftx2  vs")
print("             synth-pdb --shift-predictor empirical")
print("   and compare the output CSV files — you'll see the same effect.")

# ─── Shift Predictor Side-by-Side Comparison ────────────────────────────── # Educational note: # predict_empirical_shifts → SPARTA+-style: random-coil + SS offset + ring current # ShiftX2Predictor → external ML binary trained on 20k+ BMRB/PDB entries # References: # Shen & Bax (2010) J Biomol NMR 48:13 (SPARTA+ accuracy ~0.40 ppm ¹H) # Han et al. (2011) J Biomol NMR 50:43 (SHIFTX2 accuracy ~0.24 ppm ¹H) import io import math import matplotlib.gridspec as gridspec import matplotlib.pyplot as plt import numpy as np from synth_nmr.chemical_shifts import ShiftX2Predictor, predict_empirical_shifts # ── 1. Empirical prediction (always available) ──────────────────────────── print("⚙️ Running empirical shift prediction (SPARTA+)…") emp_shifts = predict_empirical_shifts(structure) # ── 2. SHIFTX2 prediction (optional external binary) ───────────────────── sx2 = ShiftX2Predictor() HAS_SHIFTX2 = sx2.is_available() sx2_shifts = {} if HAS_SHIFTX2: print("✅ SHIFTX2 binary detected — running prediction…") try: sx2_shifts = sx2.predict(structure) print("✅ SHIFTX2 prediction complete.") except Exception as _e: HAS_SHIFTX2 = False print(f"⚠️ SHIFTX2 prediction failed: {_e}") else: print("ℹ️ SHIFTX2 not found in PATH — showing empirical-only panel.") print(" Install: sbgrid-cli install shiftx2 OR http://www.shiftx2.ca/download.html") # ── 3. Extract (H, N) pairs for HSQC scatter ───────────────────────────── def extract_hsqc(shift_dict): """Return parallel lists of (residue_label, H_ppm, N_ppm).""" labels, h_vals, n_vals = [], [], [] chain = shift_dict.get('A', {}) for r in sorted(chain): atom_d = chain[r] if 'H' in atom_d and 'N' in atom_d: h_vals.append(atom_d['H']) n_vals.append(atom_d['N']) # residue name from structure mask = structure.res_id == r rn = structure.res_name[mask][0] if mask.any() else '?' labels.append(f"{rn}{r}") return labels, h_vals, n_vals emp_labels, emp_h, emp_n = extract_hsqc(emp_shifts) if HAS_SHIFTX2: sx2_labels, sx2_h, sx2_n = extract_hsqc(sx2_shifts) # ── 4. Plotting ─────────────────────────────────────────────────────────── n_panels = 3 if HAS_SHIFTX2 else 1 # HSQC-emp | HSQC-sx2 | delta bars fig_w = 16 if HAS_SHIFTX2 else 7 fig = plt.figure(figsize=(fig_w, 7)) if HAS_SHIFTX2: gs = gridspec.GridSpec(2, 2, figure=fig, height_ratios=[1.3, 0.9], hspace=0.45, wspace=0.35) ax_emp = fig.add_subplot(gs[0, 0]) ax_sx2 = fig.add_subplot(gs[0, 1]) ax_del = fig.add_subplot(gs[1, :]) else: ax_emp = fig.add_subplot(111) # ── Plot 1: Empirical HSQC ──────────────────────────────────────────────── def _hsqc_axes_style(ax, title): ax.set_xlim(10.8, 6.2) ax.set_ylim(137, 98) ax.set_xlabel(r"$^1$H (ppm)", fontsize=10) ax.set_ylabel(r"$^{15}$N (ppm)", fontsize=10) ax.set_title(title, fontsize=11, fontweight='bold') ax.grid(alpha=0.25) ax_emp.scatter(emp_h, emp_n, s=60, alpha=0.75, color='steelblue', edgecolors='navy', linewidths=0.5) _hsqc_axes_style(ax_emp, "Empirical (SPARTA+)") for h, n, lbl in zip(emp_h, emp_n, emp_labels): ax_emp.annotate(lbl, (h, n), xytext=(2, 2), textcoords='offset points', fontsize=6) if HAS_SHIFTX2: # ── Plot 2: SHIFTX2 HSQC ────────────────────────────────────────────── ax_sx2.scatter(sx2_h, sx2_n, s=60, alpha=0.75, color='tomato', edgecolors='darkred', linewidths=0.5) _hsqc_axes_style(ax_sx2, "SHIFTX2 (ML predictor)") for h, n, lbl in zip(sx2_h, sx2_n, sx2_labels): ax_sx2.annotate(lbl, (h, n), xytext=(2, 2), textcoords='offset points', fontsize=6) # ── Plot 3: Δδ(H) per residue ───────────────────────────────────────── # Align by residue label — only residues present in both predictors emp_by_res = dict(zip(emp_labels, zip(emp_h, emp_n))) sx2_by_res = dict(zip(sx2_labels, zip(sx2_h, sx2_n))) shared = [l for l in emp_labels if l in sx2_by_res] delta_h = [sx2_by_res[l][0] - emp_by_res[l][0] for l in shared] bar_colors = ['crimson' if d > 0 else 'royalblue' for d in delta_h] x_pos = range(len(shared)) ax_del.bar(x_pos, delta_h, color=bar_colors, alpha=0.75, edgecolor='black', linewidth=0.4) ax_del.axhline(0, color='black', linewidth=0.8) ax_del.set_xticks(list(x_pos)) ax_del.set_xticklabels(shared, rotation=90, fontsize=6) ax_del.set_ylabel(r"$\Delta\delta$(¹H) SHIFTX2 − Empirical (ppm)", fontsize=9) ax_del.set_title("Per-Residue ¹H Shift Difference (SHIFTX2 minus Empirical)", fontsize=10, fontweight='bold') ax_del.grid(axis='y', alpha=0.25) plt.suptitle("Chemical Shift Predictor Comparison", fontsize=13, fontweight='bold', y=1.01) plt.tight_layout() plt.show() # ── 5. Summary statistics table ─────────────────────────────────────────── print("\n" + "─" * 60) print(f"{'Metric':<30} {'Empirical':>12} {'SHIFTX2':>12}") print("─" * 60) def _stats(vals): if not vals: return '–', '–', '–' return round(float(np.std(vals)), 3), round(float(np.ptp(vals)), 3), round(float(np.mean(vals)), 3) emp_h_std, emp_h_range, emp_h_mean = _stats(emp_h) emp_n_std, emp_n_range, emp_n_mean = _stats(emp_n) if HAS_SHIFTX2: sx2_h_std, sx2_h_range, sx2_h_mean = _stats(sx2_h) sx2_n_std, sx2_n_range, sx2_n_mean = _stats(sx2_n) rmsd_h = round(math.sqrt(float(np.mean([(sx2_by_res[l][0]-emp_by_res[l][0])**2 for l in shared]))), 3) if shared else '–' else: sx2_h_std = sx2_h_range = sx2_h_mean = '–' sx2_n_std = sx2_n_range = sx2_n_mean = '–' rmsd_h = '–' rows = [ ("¹H std dev (ppm)", emp_h_std, sx2_h_std), ("¹H range (ppm)", emp_h_range, sx2_h_range), ("¹⁵N std dev (ppm)", emp_n_std, sx2_n_std), ("Δδ(¹H) RMSD vs empirical", '0.000', rmsd_h), ] for name, e_val, s_val in rows: print(f"{name:<30} {str(e_val):>12} {str(s_val):>12}") print("─" * 60) print() print("📚 Literature benchmark RMSDs (vs experimental BMRB shifts):") print(" SPARTA+ (empirical): ¹H ≈ 0.40 ppm, ¹⁵N ≈ 2.00 ppm") print(" SHIFTX2 (ML): ¹H ≈ 0.24 ppm, ¹⁵N ≈ 1.16 ppm") print(" Sources: Shen & Bax (2010) J Biomol NMR 48:13;") print(" Han et al. (2011) J Biomol NMR 50:43") print() print("💡 EXPERIMENTS:") print("1. DISPERSION: If SHIFTX2 is available, compare the spread of HSQC peaks.") print(" SHIFTX2 correctly predicts more ¹H dispersion for folded proteins.") print("2. CONTEXT SENSITIVITY: Look at the Δδ bar chart — bars near aromatic") print(" residues (Phe, Tyr, Trp) tend to be larger because SHIFTX2 models") print(" ring current effects from neighbouring residues more accurately.") print("3. CLI FLAG: Try synth-pdb --shift-predictor shiftx2 vs") print(" synth-pdb --shift-predictor empirical") print(" and compare the output CSV files — you'll see the same effect.")

🧲 Primer: What is "NOESY"?¶

NOESY stands for Nuclear Overhauser Effect Spectroscopy.

Think of it as "Molecular Stalking".

• Through-Space vs. Through-Bond: Most NMR experiments rely on bonds. The NOE is different—it measures how much atoms are whispering to each other through space (usually < 5 Angstroms).
• The Ruler of Structure: The strength of the NOE signal is proportional to $\frac{1}{r^6}$, where $r$ is the distance between atoms. This means it is incredibly sensitive to distance.
• Folding: If we see an NOE signal between Residue 1 and Residue 100, we know the protein chain must loop back around so those two residues are touching! This is the primary data used to calculate 3D structures.

---

🗺️ 3. The Contact Map (Virtual NOESY)¶

The Contact Map is the "Fingerprint" of a protein fold. It shows which residues are touching in 3D space.

• Diagonal: Residues are always close to their neighbors (i, i+1).
• Off-Diagonal: These are the interesting Long-Range Contacts. They tell us the protein has folded back on itself (e.g., a Beta Sheet or Alpha Helix packing).
• AI Relevance: This 2D matrix is exactly what tools like AlphaFold predict internally (the "Distogram").

In [ ]:

# Calculate Contact Map
from synth_pdb.contact import compute_contact_map

# Compute Alpha-Carbon distances (Standard for Fold Recognition)
contact_map = compute_contact_map(structure, method='ca', power=None)
binary_map = compute_contact_map(structure, method='ca', threshold=8.0, power=0)

# Plotting
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))

# 1. Distance Matrix (Continuous)
im1 = ax1.imshow(contact_map, cmap='viridis_r', interpolation='nearest')
ax1.set_title('Distance Matrix (Angstroms)')
ax1.set_xlabel('Residue Index')
ax1.set_ylabel('Residue Index')
plt.colorbar(im1, ax=ax1, label='Distance (Å)')

# 2. Binary Contact Map (AI/NMR Style)
# Threshold < 8 Angstroms is standard for CASP/AlphaFold evaluation
im2 = ax2.imshow(binary_map, cmap='Greys', interpolation='nearest')
ax2.set_title('Binary Contact Map (Threshold < 8Å)')
ax2.set_xlabel('Residue Index')
ax2.set_yticks([]) # Hide Y ticks for cleanliness

plt.tight_layout()
plt.show()

print("📊 INTERPRETATION:")
print("Dark patterns away from the diagonal indicate FOLDED structure.")
print("A Thick Diagonal with no off-diagonal spots = DISORDERED/UNFOLDED chain.")

# Calculate Contact Map from synth_pdb.contact import compute_contact_map # Compute Alpha-Carbon distances (Standard for Fold Recognition) contact_map = compute_contact_map(structure, method='ca', power=None) binary_map = compute_contact_map(structure, method='ca', threshold=8.0, power=0) # Plotting fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5)) # 1. Distance Matrix (Continuous) im1 = ax1.imshow(contact_map, cmap='viridis_r', interpolation='nearest') ax1.set_title('Distance Matrix (Angstroms)') ax1.set_xlabel('Residue Index') ax1.set_ylabel('Residue Index') plt.colorbar(im1, ax=ax1, label='Distance (Å)') # 2. Binary Contact Map (AI/NMR Style) # Threshold < 8 Angstroms is standard for CASP/AlphaFold evaluation im2 = ax2.imshow(binary_map, cmap='Greys', interpolation='nearest') ax2.set_title('Binary Contact Map (Threshold < 8Å)') ax2.set_xlabel('Residue Index') ax2.set_yticks([]) # Hide Y ticks for cleanliness plt.tight_layout() plt.show() print("📊 INTERPRETATION:") print("Dark patterns away from the diagonal indicate FOLDED structure.") print("A Thick Diagonal with no off-diagonal spots = DISORDERED/UNFOLDED chain.")

📊 3. Chemical Shift Index (CSI)¶

The Chemical Shift Index is a method to visualize secondary structure directly from the raw data, without calculating a full 3D structure.

We plot the deviation of the Alpha-Carbon ($C_\alpha$) shift from its random coil value: $$ \Delta \delta = \delta_{measured} - \delta_{random\_coil} $$

• Positive Clusters (> +0.7 ppm): Indicate Alpha Helices.
• Negative Clusters (< -0.7 ppm): Indicate Beta Sheets.
• Near Zero: Indicates Random Coil / Loops.

In [ ]:

# Calculate CSI
import synth_pdb.chemical_shifts

# Force reload to pick up the new function definition from disk
importlib.reload(synth_pdb.chemical_shifts)
from synth_pdb.chemical_shifts import calculate_csi

# 1. Get Deviations
csi_data = calculate_csi(shifts, structure)
chain_id = list(csi_data.keys())[0]
deltas = csi_data[chain_id]

res_nums = sorted(deltas.keys())
values = [deltas[r] for r in res_nums]

# 2. Plotting
plt.figure(figsize=(10, 4))
colors = ['red' if v > 0 else 'blue' for v in values]
plt.bar(res_nums, values, color=colors, alpha=0.7)

# 3. Add Threshold Lines (The "CSI" definition)
plt.axhline(0.7, color='black', linestyle='--', linewidth=1, label='Helix Threshold (+0.7)')
plt.axhline(-0.7, color='black', linestyle=':', linewidth=1, label='Sheet Threshold (-0.7)')
plt.axhline(0, color='black', linewidth=0.5)

plt.title('Chemical Shift Index (CSI) - $C_\\alpha$ Deviations')
plt.ylabel(r'$\Delta\delta_{C\alpha}$ (ppm)')
plt.xlabel('Residue Number')
plt.legend()
plt.grid(axis='y', alpha=0.3)
plt.show()

# Calculate CSI import synth_pdb.chemical_shifts # Force reload to pick up the new function definition from disk importlib.reload(synth_pdb.chemical_shifts) from synth_pdb.chemical_shifts import calculate_csi # 1. Get Deviations csi_data = calculate_csi(shifts, structure) chain_id = list(csi_data.keys())[0] deltas = csi_data[chain_id] res_nums = sorted(deltas.keys()) values = [deltas[r] for r in res_nums] # 2. Plotting plt.figure(figsize=(10, 4)) colors = ['red' if v > 0 else 'blue' for v in values] plt.bar(res_nums, values, color=colors, alpha=0.7) # 3. Add Threshold Lines (The "CSI" definition) plt.axhline(0.7, color='black', linestyle='--', linewidth=1, label='Helix Threshold (+0.7)') plt.axhline(-0.7, color='black', linestyle=':', linewidth=1, label='Sheet Threshold (-0.7)') plt.axhline(0, color='black', linewidth=0.5) plt.title('Chemical Shift Index (CSI) - $C_\\alpha$ Deviations') plt.ylabel(r'$\Delta\delta_{C\alpha}$ (ppm)') plt.xlabel('Residue Number') plt.legend() plt.grid(axis='y', alpha=0.3) plt.show()

🛡️ 4. Structure Quality & Validation¶

A structure is only useful if it obeys physics. We run a Full Validation Suite covering:

1. Bond Geometries: Lengths and angles.
2. Ramachandran Plot: Ensuring backbone torsion angles ($\phi, \psi$) are in allowed regions (alpha/beta/proline-restricted).
3. Steric Clashes: Atoms overlapping in space.
4. Rotamers: Side-chains in low-energy conformations.

In [ ]:

# Run Validation
import logging

from synth_pdb.validator import PDBValidator

# Configure logging to show only critical info in notebook
logging.basicConfig(level=logging.INFO, stream=sys.stdout, format='%(levelname)s: %(message)s', force=True)

print("🔍 Running Biophysical Validation...")
validator = PDBValidator(pdb_content)
validator.validate_all()

violations = validator.get_violations()
if not violations:
    print("✅ No Violations Found! Structure is physically sound.")
else:
    print(f"⚠️ Found {len(violations)} Violations:")
    for v in violations[:10]: # Show first 10
        print(f"  - {v}")
    if len(violations) > 10:
        print(f"  ... and {len(violations)-10} more.")

# Run Validation import logging from synth_pdb.validator import PDBValidator # Configure logging to show only critical info in notebook logging.basicConfig(level=logging.INFO, stream=sys.stdout, format='%(levelname)s: %(message)s', force=True) print("🔍 Running Biophysical Validation...") validator = PDBValidator(pdb_content) validator.validate_all() violations = validator.get_violations() if not violations: print("✅ No Violations Found! Structure is physically sound.") else: print(f"⚠️ Found {len(violations)} Violations:") for v in violations[:10]: # Show first 10 print(f" - {v}") if len(violations) > 10: print(f" ... and {len(violations)-10} more.")

In [ ]:

# Ramachandran Plot Visualization
# We extract Phi/Psi angles and plot them against the "Allowed" regions.

phi_angles = []
psi_angles = []
res_colors = []

# Extract Dihedrals (using Biotite)
from biotite.structure import dihedral_backbone

phi, psi, omega = dihedral_backbone(structure)

# Convert radians to degrees and filter non-existing (termini)
for i in range(len(phi)):
    if not np.isnan(phi[i]) and not np.isnan(psi[i]):
        p = np.degrees(phi[i])
        s = np.degrees(psi[i])
        phi_angles.append(p)
        psi_angles.append(s)

        # Color by Secondary Structure
        # Note: 'structure' array index matches phi/psi index
        # We'll use a simple heuristic for color
        if -100 < p < -30 and -80 < s < -10:
            res_colors.append('blue') # Alpha
        elif -180 < p < -40 and (90 < s < 180 or -180 < s < -160):
             res_colors.append('red') # Beta
        elif p > 0:
             res_colors.append('green') # Left-handed (Gly check)
        else:
             res_colors.append('gray')

plt.figure(figsize=(6, 6))
plt.scatter(phi_angles, psi_angles, c=res_colors, alpha=0.7, edgecolors='black')

# Underlying Regions (Simplified Boxes for reference)
plt.gca().add_patch(plt.Rectangle((-100, -70), 70, 60, color='blue', alpha=0.1, label='Alpha'))
plt.gca().add_patch(plt.Rectangle((-180, 90), 140, 90, color='red', alpha=0.1, label='Beta'))

plt.xlim(-180, 180)
plt.ylim(-180, 180)
plt.axhline(0, color='black', linewidth=0.5)
plt.axvline(0, color='black', linewidth=0.5)
plt.xlabel(r'Phi ($\phi$)')
plt.ylabel(r'Psi ($\psi$)')
plt.title('Ramachandran Plot (Backbone Geometry)')
plt.grid(True, alpha=0.3)
plt.legend(loc='upper right')
plt.show()

print("📊 INTERPRETATION:")
print("Blue Dots: Alpha Helices (Bottom Left)")
print("Red Dots: Beta Sheets (Top Left)")
print("Green Dots: Left-handed Helices (Top Right - Usually only Glycine)")

# Ramachandran Plot Visualization # We extract Phi/Psi angles and plot them against the "Allowed" regions. phi_angles = [] psi_angles = [] res_colors = [] # Extract Dihedrals (using Biotite) from biotite.structure import dihedral_backbone phi, psi, omega = dihedral_backbone(structure) # Convert radians to degrees and filter non-existing (termini) for i in range(len(phi)): if not np.isnan(phi[i]) and not np.isnan(psi[i]): p = np.degrees(phi[i]) s = np.degrees(psi[i]) phi_angles.append(p) psi_angles.append(s) # Color by Secondary Structure # Note: 'structure' array index matches phi/psi index # We'll use a simple heuristic for color if -100 < p < -30 and -80 < s < -10: res_colors.append('blue') # Alpha elif -180 < p < -40 and (90 < s < 180 or -180 < s < -160): res_colors.append('red') # Beta elif p > 0: res_colors.append('green') # Left-handed (Gly check) else: res_colors.append('gray') plt.figure(figsize=(6, 6)) plt.scatter(phi_angles, psi_angles, c=res_colors, alpha=0.7, edgecolors='black') # Underlying Regions (Simplified Boxes for reference) plt.gca().add_patch(plt.Rectangle((-100, -70), 70, 60, color='blue', alpha=0.1, label='Alpha')) plt.gca().add_patch(plt.Rectangle((-180, 90), 140, 90, color='red', alpha=0.1, label='Beta')) plt.xlim(-180, 180) plt.ylim(-180, 180) plt.axhline(0, color='black', linewidth=0.5) plt.axvline(0, color='black', linewidth=0.5) plt.xlabel(r'Phi ($\phi$)') plt.ylabel(r'Psi ($\psi$)') plt.title('Ramachandran Plot (Backbone Geometry)') plt.grid(True, alpha=0.3) plt.legend(loc='upper right') plt.show() print("📊 INTERPRETATION:") print("Blue Dots: Alpha Helices (Bottom Left)") print("Red Dots: Beta Sheets (Top Left)") print("Green Dots: Left-handed Helices (Top Right - Usually only Glycine)")

---

🔭 The Structural Frontier: Advanced Observables & Machine Learning¶

Welcome to the cutting edge of NMR biophysics. In this section, we will explore two advanced topics:

1. Residual Dipolar Couplings (RDCs): A global structural restraint.
2. Neural Network Chemical Shifts: Using Machine Learning to predict spectra.

📐 1. Residual Dipolar Couplings (RDCs)¶

Normally, molecules wildly tumble in solution, averaging out the magnetic interactions between nuclei (Dipolar Couplings) to exactly zero.

However, if we add an Alignment Medium (like a dilute liquid crystal, or Pf1 bacteriophage) to our sample tube, the protein bumps into the medium and aligns just slightly (maybe 1 in 1000 molecules is aligned at any given time).

This incomplete averaging leaves a small Residual Dipolar Coupling (RDC).

Why is this amazing? While NOEs measure distance, RDCs measure angle relative to a global coordinate system (the Alignment Tensor). This gives us long-range information about how different parts of the protein are oriented relative to each other, even if they are far apart in space!

The RDC value depends on:

1. The angle of the bond vector (e.g., N-H) relative to the magnetic field.
2. The Axial Magnitude ($D_a$): How strongly the protein is aligned.
3. The Rhombicity ($R$): How asymmetrical the alignment is.

Below you will see a bar chart of per-residue RDC values and a 3D color-by-RDC view of the structure (red = positive, blue = negative).

In [ ]:

import os

IN_CI = bool(os.getenv("CI"))

if not IN_CI:
    # Interactive RDC Explorer
    # Educational note: We use synth_pdb.rdc (the shim module) instead of
    # importing synth_nmr.rdc directly, to stay consistent with the rest of
    # this notebook and to benefit from synth_pdb's input validation.
    import io
    import io as _sio
    import os

    # Serialize the structure to a PDB string so py3Dmol can display it.
    # Deriving it here keeps this cell self-contained -- no reliance on a
    # variable that may or may not have been set in an earlier cell.
    import biotite.structure.io.pdb as _pdbio
    import ipywidgets as widgets
    import numpy as np
    import py3Dmol
    from IPython.display import display

    from synth_pdb.rdc import calculate_rdcs
    _pdb_buf = _sio.StringIO()
    _pdb_out  = _pdbio.PDBFile()
    _pdbio.set_structure(_pdb_out, structure)
    _pdb_out.write(_pdb_buf)
    pdb_string = _pdb_buf.getvalue()

    RDCS_TWO_THIRDS = 2.0 / 3.0  # Physical upper limit for rhombicity R

    # ─── Widgets ──────────────────────────────────────────────────────────────
    # Da slider: positive only. Convention is Da > 0; sign of individual RDCs
    # comes from the bond-vector angle, not from Da itself.
    da_slider = widgets.FloatSlider(
        min=1.0, max=20.0, step=0.5, value=10.0,
        description='Da (Hz)',
        tooltip='Axial magnitude of alignment tensor (Hz). Controls RDC scale.',
        style={'description_width': '70px'},
        layout=widgets.Layout(width='370px'),
        continuous_update=False,
    )
    r_slider = widgets.FloatSlider(
        min=0.0, max=0.665, step=0.025, value=0.15,
        description='R (rhombicity)',
        tooltip='Rhombicity: 0 = axially symmetric, 0.667 = maximum asymmetry.',
        style={'description_width': '100px'},
        layout=widgets.Layout(width='400px'),
        continuous_update=False,
    )

    # Image widget -- bar chart renders here without flickering
    rdc_plot_image = widgets.Image(format='png', width=820, height=370)

    # ─── Bar Chart ────────────────────────────────────────────────────────────
    def update_rdc_plot(change=None):
        da = da_slider.value
        # Clamp R defensively: slider steps can accumulate infinitesimally above 2/3,
        # which the synth-nmr validator rejects.
        r  = float(np.clip(r_slider.value, 0.0, RDCS_TWO_THIRDS - 1e-9))

        plt.ioff()

        rdcs = calculate_rdcs(structure, Da=da, R=r)
        if not rdcs:
            return

        r_nums = sorted(rdcs.keys())
        r_vals = [rdcs[i] for i in r_nums]

        # Physical limits for this tensor (Tjandra & Bax 1997, Science 278:1111)
        #   Maximum:  +2*Da        (bond vector parallel to Z-axis, theta=0 deg)
        #   Minimum:  -Da*(1+1.5R) (bond vector in XY-plane at phi=90, theta=90 deg)
        rdc_max =  2.0 * da
        rdc_min = -da * (1.0 + 1.5 * r)

        # Color bars by sign (red = positive, blue = negative)
        bar_colors = ['#e74c3c' if v >= 0 else '#2980b9' for v in r_vals]

        fig, ax = plt.subplots(figsize=(10, 4))

        # Secondary structure background shading
        for ind, ss in enumerate(ss_types):
            if ind < len(r_nums):
                if ss == 'alpha':
                    ax.axvspan(r_nums[ind] - 0.5, r_nums[ind] + 0.5, color='blue', alpha=0.08)
                elif ss == 'beta':
                    ax.axvspan(r_nums[ind] - 0.5, r_nums[ind] + 0.5, color='red', alpha=0.08)

        ax.bar(r_nums, r_vals, color=bar_colors, alpha=0.82, edgecolor='white', linewidth=0.5)

        # Physical limit dashed lines -- no measured value can exceed these bounds
        ax.axhline(rdc_max, color='#c0392b', linestyle='--', linewidth=1.4, alpha=0.75,
                   label=f'+2*Da = {rdc_max:.1f} Hz  (theta=0 deg, Z-axis)')
        ax.axhline(rdc_min, color='#1a5276', linestyle='--', linewidth=1.4, alpha=0.75,
                   label=f'Min = {rdc_min:.1f} Hz  (theta=90 deg, phi=90 deg)')
        ax.axhline(0, color='black', linewidth=0.7, alpha=0.45)

        ax.set_ylabel('1D_NH (Hz)', fontsize=11)
        ax.set_xlabel('Residue Number', fontsize=11)
        ax.set_title(
            f'N-H Residual Dipolar Couplings  |  Da = {da:.1f} Hz,  R = {r:.3f}\n'
            f'Red = positive (N-H near Z-axis)   Blue = negative (N-H near XY-plane)',
            fontsize=10,
        )
        ax.set_xlim(min(r_nums) - 0.8, max(r_nums) + 0.8)
        ax.set_ylim(rdc_min * 1.18, rdc_max * 1.18)
        ax.legend(fontsize=8, loc='upper right', framealpha=0.8)
        ax.grid(axis='y', linestyle='--', alpha=0.3)

        plt.tight_layout()
        buf = io.BytesIO()
        fig.savefig(buf, format='png', bbox_inches='tight')
        plt.close(fig)
        buf.seek(0)
        rdc_plot_image.value = buf.read()


    # ─── py3Dmol Color-by-RDC View ────────────────────────────────────────────
    # Colors the protein cartoon by the sign and magnitude of the N-H RDC at
    # each residue, connecting the bar chart to the 3D structure visually.
    # Red = positive RDC (N-H near tensor Z-axis)
    # Blue = negative RDC (N-H in XY-plane)
    # White = near-zero RDC (near the magic angle, theta ~ 54.7 deg)
    rdc_3d_output = widgets.Output()

    def update_rdc_3d(change=None):
        da = da_slider.value
        r  = float(np.clip(r_slider.value, 0.0, RDCS_TWO_THIRDS - 1e-9))
        rdcs = calculate_rdcs(structure, Da=da, R=r)
        if not rdcs:
            return

        # Normalize to the SYMMETRIC max so both red and blue extremes are used.
        # Using the data-max (not physical limit) keeps the color range tight.
        scale = max(abs(v) for v in rdcs.values()) or 1.0

        def rdc_to_hex(val):
            # Map [-scale, +scale] to blue-white-red.
            # t in [-1, 1]: positive -> red, negative -> blue, 0 -> white.
            t = max(-1.0, min(1.0, val / scale))
            if t >= 0:
                # White -> Red: keep r=255, fade g and b
                r_c = 255
                g_c = int(255 * (1.0 - t))
                b_c = int(255 * (1.0 - t))
            else:
                # White -> Blue: keep b=255, fade r and g
                t2  = -t
                r_c = int(255 * (1.0 - t2))
                g_c = int(255 * (1.0 - t2))
                b_c = 255
            return f'#{r_c:02x}{g_c:02x}{b_c:02x}'

        # Build ONE setStyle call per residue using setStyle (not addStyle).
        # Using addStyle accumulates across slider updates -- previous color calls
        # stack up in the py3Dmol command queue and older colors bleed through.
        # setStyle with a specific selector REPLACES the style for that residue.
        view = py3Dmol.view(width=820, height=380)
        view.addModel(pdb_string, 'pdb')

        # First pass: set all residues to white (fallback for missing/PRO residues)
        view.setStyle({}, {'cartoon': {'color': 'white'}})

        # Second pass: override each residue with its RDC-derived color.
        # Each setStyle call is scoped to a single residue, so they don't interfere.
        for res_id, val in rdcs.items():
            color = rdc_to_hex(val)
            view.setStyle({'resi': str(res_id)}, {'cartoon': {'color': color}})

        view.zoomTo()

        with rdc_3d_output:
            rdc_3d_output.clear_output(wait=True)
            view.show()


    # ─── Bind events and render ────────────────────────────────────────────────
    da_slider.observe(update_rdc_plot, names='value')
    r_slider.observe(update_rdc_plot,  names='value')
    da_slider.observe(update_rdc_3d,   names='value')
    r_slider.observe(update_rdc_3d,    names='value')

    update_rdc_plot()

    # The 3Dmol.js library initializes asynchronously in the browser.
    # Calling view.show() immediately at kernel execution time races with
    # the library bootstrap and shows the 'failed to load' error.
    # Instead, show a placeholder; the first slider drag fires update_rdc_3d()
    # after 3Dmol.js is ready.
    with rdc_3d_output:
        from IPython.display import HTML as _HTML
        from IPython.display import display as _disp
        _disp(_HTML(
            '<div style="text-align:center;padding:40px;color:#555;'
            'border:1px dashed #bbb;border-radius:8px;font-style:italic">'
            '&#128260; Move a slider above to load the 3D color-by-RDC view'
            '</div>'
        ))

    display(widgets.VBox([
        widgets.HTML(
            '<b>RDC Fingerprint -- drag sliders to explore</b><br>'
            '<i>Background shading: blue = alpha-helix regions, red = beta-strand regions</i>'
        ),
        widgets.HBox([da_slider, r_slider]),
        rdc_plot_image,
        widgets.HTML('<b>3D Color-by-RDC View</b> (cartoon colored by RDC sign and magnitude)'),
        rdc_3d_output,
    ]))

    print()
    print('EXPERIMENTS:')
    print('1. HELIX SIGNATURE: In the alpha-helix region (blue background), notice the')
    print('   RDC values vary smoothly. Consecutive N-H bonds rotate ~100 deg per residue')
    print('   in a helix, producing a characteristic quasi-sinusoidal RDC pattern.')
    print('2. SCALE WITH Da: Set R=0.0 and sweep Da. All bars scale proportionally.')
    print('   Same structural information, different magnitude -- set by alignment strength.')
    print('3. ASYMMETRY WITH R: Set Da=10, then sweep R from 0 to 0.665.')
    print('   The rhombicity term breaks X-Y symmetry. Residues that were mild become')
    print('   extreme -- R reveals orientation information that an axial tensor cannot.')
    print('4. 3D VIEW: Red residues have N-H bonds near the tensor Z-axis;')
    print('   blue residues have N-H bonds in the XY-plane. Secondary structure')
    print('   blocks should show coherent color patches.')

# (Widget output skipped in CI)

import os IN_CI = bool(os.getenv("CI")) if not IN_CI: # Interactive RDC Explorer # Educational note: We use synth_pdb.rdc (the shim module) instead of # importing synth_nmr.rdc directly, to stay consistent with the rest of # this notebook and to benefit from synth_pdb's input validation. import io import io as _sio import os # Serialize the structure to a PDB string so py3Dmol can display it. # Deriving it here keeps this cell self-contained -- no reliance on a # variable that may or may not have been set in an earlier cell. import biotite.structure.io.pdb as _pdbio import ipywidgets as widgets import numpy as np import py3Dmol from IPython.display import display from synth_pdb.rdc import calculate_rdcs _pdb_buf = _sio.StringIO() _pdb_out = _pdbio.PDBFile() _pdbio.set_structure(_pdb_out, structure) _pdb_out.write(_pdb_buf) pdb_string = _pdb_buf.getvalue() RDCS_TWO_THIRDS = 2.0 / 3.0 # Physical upper limit for rhombicity R # ─── Widgets ────────────────────────────────────────────────────────────── # Da slider: positive only. Convention is Da > 0; sign of individual RDCs # comes from the bond-vector angle, not from Da itself. da_slider = widgets.FloatSlider( min=1.0, max=20.0, step=0.5, value=10.0, description='Da (Hz)', tooltip='Axial magnitude of alignment tensor (Hz). Controls RDC scale.', style={'description_width': '70px'}, layout=widgets.Layout(width='370px'), continuous_update=False, ) r_slider = widgets.FloatSlider( min=0.0, max=0.665, step=0.025, value=0.15, description='R (rhombicity)', tooltip='Rhombicity: 0 = axially symmetric, 0.667 = maximum asymmetry.', style={'description_width': '100px'}, layout=widgets.Layout(width='400px'), continuous_update=False, ) # Image widget -- bar chart renders here without flickering rdc_plot_image = widgets.Image(format='png', width=820, height=370) # ─── Bar Chart ──────────────────────────────────────────────────────────── def update_rdc_plot(change=None): da = da_slider.value # Clamp R defensively: slider steps can accumulate infinitesimally above 2/3, # which the synth-nmr validator rejects. r = float(np.clip(r_slider.value, 0.0, RDCS_TWO_THIRDS - 1e-9)) plt.ioff() rdcs = calculate_rdcs(structure, Da=da, R=r) if not rdcs: return r_nums = sorted(rdcs.keys()) r_vals = [rdcs[i] for i in r_nums] # Physical limits for this tensor (Tjandra & Bax 1997, Science 278:1111) # Maximum: +2*Da (bond vector parallel to Z-axis, theta=0 deg) # Minimum: -Da*(1+1.5R) (bond vector in XY-plane at phi=90, theta=90 deg) rdc_max = 2.0 * da rdc_min = -da * (1.0 + 1.5 * r) # Color bars by sign (red = positive, blue = negative) bar_colors = ['#e74c3c' if v >= 0 else '#2980b9' for v in r_vals] fig, ax = plt.subplots(figsize=(10, 4)) # Secondary structure background shading for ind, ss in enumerate(ss_types): if ind < len(r_nums): if ss == 'alpha': ax.axvspan(r_nums[ind] - 0.5, r_nums[ind] + 0.5, color='blue', alpha=0.08) elif ss == 'beta': ax.axvspan(r_nums[ind] - 0.5, r_nums[ind] + 0.5, color='red', alpha=0.08) ax.bar(r_nums, r_vals, color=bar_colors, alpha=0.82, edgecolor='white', linewidth=0.5) # Physical limit dashed lines -- no measured value can exceed these bounds ax.axhline(rdc_max, color='#c0392b', linestyle='--', linewidth=1.4, alpha=0.75, label=f'+2*Da = {rdc_max:.1f} Hz (theta=0 deg, Z-axis)') ax.axhline(rdc_min, color='#1a5276', linestyle='--', linewidth=1.4, alpha=0.75, label=f'Min = {rdc_min:.1f} Hz (theta=90 deg, phi=90 deg)') ax.axhline(0, color='black', linewidth=0.7, alpha=0.45) ax.set_ylabel('1D_NH (Hz)', fontsize=11) ax.set_xlabel('Residue Number', fontsize=11) ax.set_title( f'N-H Residual Dipolar Couplings | Da = {da:.1f} Hz, R = {r:.3f}\n' f'Red = positive (N-H near Z-axis) Blue = negative (N-H near XY-plane)', fontsize=10, ) ax.set_xlim(min(r_nums) - 0.8, max(r_nums) + 0.8) ax.set_ylim(rdc_min * 1.18, rdc_max * 1.18) ax.legend(fontsize=8, loc='upper right', framealpha=0.8) ax.grid(axis='y', linestyle='--', alpha=0.3) plt.tight_layout() buf = io.BytesIO() fig.savefig(buf, format='png', bbox_inches='tight') plt.close(fig) buf.seek(0) rdc_plot_image.value = buf.read() # ─── py3Dmol Color-by-RDC View ──────────────────────────────────────────── # Colors the protein cartoon by the sign and magnitude of the N-H RDC at # each residue, connecting the bar chart to the 3D structure visually. # Red = positive RDC (N-H near tensor Z-axis) # Blue = negative RDC (N-H in XY-plane) # White = near-zero RDC (near the magic angle, theta ~ 54.7 deg) rdc_3d_output = widgets.Output() def update_rdc_3d(change=None): da = da_slider.value r = float(np.clip(r_slider.value, 0.0, RDCS_TWO_THIRDS - 1e-9)) rdcs = calculate_rdcs(structure, Da=da, R=r) if not rdcs: return # Normalize to the SYMMETRIC max so both red and blue extremes are used. # Using the data-max (not physical limit) keeps the color range tight. scale = max(abs(v) for v in rdcs.values()) or 1.0 def rdc_to_hex(val): # Map [-scale, +scale] to blue-white-red. # t in [-1, 1]: positive -> red, negative -> blue, 0 -> white. t = max(-1.0, min(1.0, val / scale)) if t >= 0: # White -> Red: keep r=255, fade g and b r_c = 255 g_c = int(255 * (1.0 - t)) b_c = int(255 * (1.0 - t)) else: # White -> Blue: keep b=255, fade r and g t2 = -t r_c = int(255 * (1.0 - t2)) g_c = int(255 * (1.0 - t2)) b_c = 255 return f'#{r_c:02x}{g_c:02x}{b_c:02x}' # Build ONE setStyle call per residue using setStyle (not addStyle). # Using addStyle accumulates across slider updates -- previous color calls # stack up in the py3Dmol command queue and older colors bleed through. # setStyle with a specific selector REPLACES the style for that residue. view = py3Dmol.view(width=820, height=380) view.addModel(pdb_string, 'pdb') # First pass: set all residues to white (fallback for missing/PRO residues) view.setStyle({}, {'cartoon': {'color': 'white'}}) # Second pass: override each residue with its RDC-derived color. # Each setStyle call is scoped to a single residue, so they don't interfere. for res_id, val in rdcs.items(): color = rdc_to_hex(val) view.setStyle({'resi': str(res_id)}, {'cartoon': {'color': color}}) view.zoomTo() with rdc_3d_output: rdc_3d_output.clear_output(wait=True) view.show() # ─── Bind events and render ──────────────────────────────────────────────── da_slider.observe(update_rdc_plot, names='value') r_slider.observe(update_rdc_plot, names='value') da_slider.observe(update_rdc_3d, names='value') r_slider.observe(update_rdc_3d, names='value') update_rdc_plot() # The 3Dmol.js library initializes asynchronously in the browser. # Calling view.show() immediately at kernel execution time races with # the library bootstrap and shows the 'failed to load' error. # Instead, show a placeholder; the first slider drag fires update_rdc_3d() # after 3Dmol.js is ready. with rdc_3d_output: from IPython.display import HTML as _HTML from IPython.display import display as _disp _disp(_HTML( '

' '🔄 Move a slider above to load the 3D color-by-RDC view' '

' )) display(widgets.VBox([ widgets.HTML( 'RDC Fingerprint -- drag sliders to explore
' 'Background shading: blue = alpha-helix regions, red = beta-strand regions' ), widgets.HBox([da_slider, r_slider]), rdc_plot_image, widgets.HTML('3D Color-by-RDC View (cartoon colored by RDC sign and magnitude)'), rdc_3d_output, ])) print() print('EXPERIMENTS:') print('1. HELIX SIGNATURE: In the alpha-helix region (blue background), notice the') print(' RDC values vary smoothly. Consecutive N-H bonds rotate ~100 deg per residue') print(' in a helix, producing a characteristic quasi-sinusoidal RDC pattern.') print('2. SCALE WITH Da: Set R=0.0 and sweep Da. All bars scale proportionally.') print(' Same structural information, different magnitude -- set by alignment strength.') print('3. ASYMMETRY WITH R: Set Da=10, then sweep R from 0 to 0.665.') print(' The rhombicity term breaks X-Y symmetry. Residues that were mild become') print(' extreme -- R reveals orientation information that an axial tensor cannot.') print('4. 3D VIEW: Red residues have N-H bonds near the tensor Z-axis;') print(' blue residues have N-H bonds in the XY-plane. Secondary structure') print(' blocks should show coherent color patches.') # (Widget output skipped in CI)

🤖 2. Neural Network Shift Predictor (The Machine Learning Frontier)¶

The Chemical Shift Index (CSI) we calculated earlier relies on a simple physics-based "Lookup Table" (e.g., SPARTA-lite): If amino acid X is in an alpha helix, add Y ppm to the baseline.

The Problem: This ignores almost all sequence context. The chemical shift of an Alanine in an alpha helix is slightly different depending on whether the neighboring residue is a bulky Tryptophan or a tiny Glycine.

The ML Solution: We use a Multilayer Perceptron (MLP) Neural Network! The network takes the explicit atomic geometry (dihedral angles $\phi, \psi$), the exact amino acid type, AND the sequence-neighbor identities (using one-hot encoding). It learns a non-linear continuous surface rather than a discrete step-function, predicting a $\Delta$CS correction term to add to the baseline.

In [ ]:

# Neural Array Predictor
try:
    import torch
    HAS_TORCH = True
except ImportError:
    HAS_TORCH = False
    print("⚠️ PyTorch is not installed. To run the Neural Simulator, you need to install it.")
    print("   Run: !pip install torch")

HAS_TORCH_GEO = False
if HAS_TORCH:
    try:
        import torch_geometric  # noqa: F401
        HAS_TORCH_GEO = True
    except ImportError:
        print("⚠️  torch_geometric is not installed — the GNN-based NeuralShiftPredictor")
        print("   requires it. Install the full ML extras with:")
        print("   !pip install synth-nmr[ml]")
        print("   (Note: torch_geometric must match your torch version — see")
        print("    https://pytorch-geometric.readthedocs.io/en/latest/install/installation.html )")
        print("   Skipping Neural Shift section.")

if HAS_TORCH and HAS_TORCH_GEO:
    from synth_nmr.neural_shifts import NeuralShiftPredictor

    print("🧠 Loading Neural Network Estimator...")
    # We initialize the predictor. If it has a pre-trained checkpoint, it will load it.
    # If not, it instantiates with random weights (predicting unstructured noise).
    neural_predictor = NeuralShiftPredictor()

    print("⚙️ Running Neural Inference (MLP forward pass)...")
    neural_shifts = neural_predictor.predict(structure)

    # Let's compare Empirical (Physics-only) vs Neural (Physics + ML) for the C-beta atom.
    # C-beta (CB) is highly sensitive to the rotamer state (side-chain geometry), making it a prime target for ML.

    emp_cb = []
    neu_cb = []
    plot_res = []

    for r in res_ids:
        # Check Chain A
        if 'A' in shifts and r in shifts['A'] and 'CB' in shifts['A'][r]:
            if 'A' in neural_shifts and r in neural_shifts['A'] and 'CB' in neural_shifts['A'][r]:
                emp_cb.append(shifts['A'][r]['CB'])
                neu_cb.append(neural_shifts['A'][r]['CB'])
                plot_res.append(r)

    if plot_res:
        # Calculate differences
        deltas = [neu - emp for neu, emp in zip(neu_cb, emp_cb)]

        plt.figure(figsize=(10, 4))
        colors = ['red' if d > 0 else 'blue' for d in deltas]
        plt.bar(plot_res, deltas, color=colors, alpha=0.7)
        plt.axhline(0, color='black', linewidth=1)
        plt.ylabel(r"$\Delta$CS (Neural - Empirical) (ppm)")
        plt.xlabel('Residue Number')
        plt.title('Neural Network Learned Corrections for $C_\\beta$ Shifts')
        plt.grid(axis='y', alpha=0.3)
        plt.tight_layout()
        plt.show()

        print("📊 INTERPRETATION:")
        print(r"The bars represent the $\Delta$CS 'correction' that the Neural Network applied over the basic empirical model.")
        print("Even if the model weights are untrained (random), you can see it extracting feature variance!")

# Neural Array Predictor try: import torch HAS_TORCH = True except ImportError: HAS_TORCH = False print("⚠️ PyTorch is not installed. To run the Neural Simulator, you need to install it.") print(" Run: !pip install torch") HAS_TORCH_GEO = False if HAS_TORCH: try: import torch_geometric # noqa: F401 HAS_TORCH_GEO = True except ImportError: print("⚠️ torch_geometric is not installed — the GNN-based NeuralShiftPredictor") print(" requires it. Install the full ML extras with:") print(" !pip install synth-nmr[ml]") print(" (Note: torch_geometric must match your torch version — see") print(" https://pytorch-geometric.readthedocs.io/en/latest/install/installation.html )") print(" Skipping Neural Shift section.") if HAS_TORCH and HAS_TORCH_GEO: from synth_nmr.neural_shifts import NeuralShiftPredictor print("🧠 Loading Neural Network Estimator...") # We initialize the predictor. If it has a pre-trained checkpoint, it will load it. # If not, it instantiates with random weights (predicting unstructured noise). neural_predictor = NeuralShiftPredictor() print("⚙️ Running Neural Inference (MLP forward pass)...") neural_shifts = neural_predictor.predict(structure) # Let's compare Empirical (Physics-only) vs Neural (Physics + ML) for the C-beta atom. # C-beta (CB) is highly sensitive to the rotamer state (side-chain geometry), making it a prime target for ML. emp_cb = [] neu_cb = [] plot_res = [] for r in res_ids: # Check Chain A if 'A' in shifts and r in shifts['A'] and 'CB' in shifts['A'][r]: if 'A' in neural_shifts and r in neural_shifts['A'] and 'CB' in neural_shifts['A'][r]: emp_cb.append(shifts['A'][r]['CB']) neu_cb.append(neural_shifts['A'][r]['CB']) plot_res.append(r) if plot_res: # Calculate differences deltas = [neu - emp for neu, emp in zip(neu_cb, emp_cb)] plt.figure(figsize=(10, 4)) colors = ['red' if d > 0 else 'blue' for d in deltas] plt.bar(plot_res, deltas, color=colors, alpha=0.7) plt.axhline(0, color='black', linewidth=1) plt.ylabel(r"$\Delta$CS (Neural - Empirical) (ppm)") plt.xlabel('Residue Number') plt.title('Neural Network Learned Corrections for $C_\\beta$ Shifts') plt.grid(axis='y', alpha=0.3) plt.tight_layout() plt.show() print("📊 INTERPRETATION:") print(r"The bars represent the $\Delta$CS 'correction' that the Neural Network applied over the basic empirical model.") print("Even if the model weights are untrained (random), you can see it extracting feature variance!")