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)

import plotly
import plotly.graph_objs as go
from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot
init_notebook_mode(connected=True)

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 nmrsim.plt import mplplot, mplplot_lineshape

mplplot(aaxx.peaklist());

[<matplotlib.lines.Line2D object at 0x623870650>]

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)

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

[<matplotlib.lines.Line2D object at 0x11cfb50d0>]

These spin systems can be combined into a spectrum:

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

[<matplotlib.lines.Line2D object at 0x11cf62890>]

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 0x62327a710>]

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

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

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

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));

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));