Demo: Simulation of Tyrosine NMR Spectrum

This notebook shows how the nmrsim library can be used to compose an entire 1H NMR spectrum from scratch.

The nmrsim.plt routines are convenient for quick plots, but for entire spectrums their small size and low resolution is noticeable (e.g. misleading signal intensities).

{TODO: provide ways to customize the plots (e.g. have ``plt.mplplot`` return the actual matplotlib object for customization, or use the peaklist data in another visualization library).}

This tutorial is adapted from the nmrmint tutorial.

(If you’re interested in an app for the simulation of a complete NMR spectrum, see the `nmrmit project <https://github.com/sametz/nmrmint>`__.)

import os
import sys
import numpy as np
import matplotlib as mpl
mpl.rcParams['figure.dpi']= 300
%config InlineBackend.figure_format = 'svg'  # makes inline plot look less blurry
%matplotlib inline
home_path = os.path.abspath(os.path.join('..', '..', '..'))
if home_path not in sys.path:
    sys.path.append(home_path)

tests_path = os.path.abspath(os.path.join('..', '..', '..', 'tests'))
if tests_path not in sys.path:
    sys.path.append(tests_path)

Here is the data for the spectrum of tyrosine in D2O:

1H NMR (500 MHz, Deuterium Oxide) δ 7.18 (d, J = 8.5 Hz, 1H), 6.89 (d, J = 8.5 Hz, 1H), 3.93 (dd, J = 7.7, 5.1 Hz, 1H),
3.19 (dd, J = 14.7, 5.1 Hz, 1H), 3.05 (dd, J = 14.7, 7.8 Hz, 1H).

Data is provided in ppm on a 500 MHz spectrometer. We’ll create a function to perform ppm-to-Hz conversions for us:

def ppm_to_hz(ppm, spec_freq):
    """Given a chemical shift in ppm and spectrometer frequency in MHz, return the corresponding chemical shift in Hz."""
    return [d * spec_freq for d in ppm]

The two “doublets” in the aromatic region actually comprise an AA’XX’ system. This 4-nuclei spin system can be modeled using the SpinSystem class:

from nmrsim import SpinSystem

Create a frequency list (in Hz) for the A, A’, X, and X’ nuclei:

v_aaxx = ppm_to_hz([7.18, 7.18, 6.89, 6.89], 500)
v_aaxx

[3590.0, 3590.0, 3445.0, 3445.0]

For the J values, as a first approximation we’ll assume JAX (an JA’X’) are close to the faux-doublet splitting of 8.5 Hz. We’ll estimate that JAA’ and JXX’ are about 2 Hz, and that the JAX’ and JA’X couplings are about 0 Hz.

j_aaxx = [[0, 2, 8.5, 0],
          [2, 0, 0, 8.5],
          [8.5, 0, 0, 2],
          [0, 8.5, 2, 0]]

aaxx = SpinSystem(v_aaxx, j_aaxx)

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

from nmrsim.plt import mplplot, mplplot_lineshape

mplplot(aaxx.peaklist());

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
[<matplotlib.lines.Line2D object at 0x7f3fe0ba3690>]

Next, we’ll create the ABX system for the aliphatic protons. For this exercise, we are assuming that the coupling constants that the first-order analysis provided are close enough.

(If accuracy is critical, there are methods for solving the ABX system. For example, see https://www.chem.wisc.edu/areas/reich/nmr/05-hmr-12-abx.htm#solving%20ABX )

v_abx = ppm_to_hz([3.93,3.19, 3.05], 500)
j_abx = [[0, 5.1, 7.75],
        [5.1, 0, -14.7],  # geminal Js should be negative
        [7.75, -14.7, 0]]
abx = SpinSystem(v_abx, j_abx)

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

mplplot(abx.peaklist(), y_max=0.2);

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
[<matplotlib.lines.Line2D object at 0x7f3fde429290>]

These spin systems can be combined into a spectrum:

tyr_spectrum = aaxx + abx
mplplot(tyr_spectrum.peaklist(), y_max=0.2)
type(tyr_spectrum)

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
[<matplotlib.lines.Line2D object at 0x7f3fde3ad9d0>]

nmrsim._classes.Spectrum

Addition of the two SpinSystem objects returned a Spectrum object.

If peak intensities look off, try using more data points for the lineshape. Here is the same example with ~ 10 data points per Hz:

points=int((tyr_spectrum.vmax - tyr_spectrum.vmin) * 10)
print(points)
mplplot(tyr_spectrum.peaklist(), y_max=0.5, points=points);

21838
[<matplotlib.lines.Line2D object at 0x7f3fde3bc910>]

The Spectrum class can also provide lineshape data for the spectrum:

mplplot_lineshape(*tyr_spectrum.lineshape(points=points));

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

The Spectrum.linewidth() method has an advantage over the .peaklist() method: it can take into account the linewidths specified by its component Multiplet/SpinSystem objects. The default value is 0.5 Hz, but this can be set to other values.

In D2O, the -OH and -NH protons are exchanged for D and are not seen in the spectrum. If we wanted to include these in the spectrum for pedagogical reasons, we could create broad singlets with the Multiplet class:

from nmrsim import Multiplet

# frequency in Hz, integration, [empty list for no coupling constants], peakwidth = 20 Hz
nh3 = Multiplet(8.3 * 500, 3, [], 20)
tyr_oh = Multiplet(9.8 * 500, 1, [], 10)
tyr_spectrum2 = tyr_spectrum + nh3 + tyr_oh

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

A Spectrum can have its .vmin and .vmax attributes reset to give a full spectral window (defaults are to provide a 50 Hz margin):

tyr_spectrum2.default_limits()  # resets limits, and returns vmin, vmax tuple

(1462.9870013439968, 4950.0)

points2 = int((tyr_spectrum2.vmax - tyr_spectrum2.vmin) * 10)
mplplot_lineshape(*tyr_spectrum2.lineshape(points=points2));

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>

What if you want the x axis to be in ppm?

# A future version of nmrsim should extend the API to facilitate using ppm in simulations.
# For now, simulations use Hz only, and ppm conversions need to be done manually.

tyr_spectrum2.vmin = -0.5 * 500
tyr_spectrum2.vmax = 10.5 * 500
x, y = tyr_spectrum2.lineshape(points=50000)
x_ppm = x / 500
mplplot_lineshape(x_ppm, y, limits=(-0.5, 10.5));

From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>
From hamiltonian_sparse:
Lz is type:  <class 'sparse._coo.core.COO'>
Lproduct is type:  <class 'sparse._coo.core.COO'>