¹³C 2D MAT NMR of L-Histidine

The following is an illustration for fitting 2D MAT/PASS datasets. The example dataset is a \(^{13}\text{C}\) 2D MAT spectrum of L-Histidine from Walder et al. 1

import numpy as np
import csdmpy as cp
import matplotlib.pyplot as plt
from lmfit import Minimizer, report_fit

from mrsimulator import Simulator
from mrsimulator.methods import SSB2D
from mrsimulator import signal_processing as sp
from mrsimulator.utils import spectral_fitting as sf
from mrsimulator.utils import get_spectral_dimensions
from mrsimulator.utils.collection import single_site_system_generator

Import the dataset

filename = "https://sandbox.zenodo.org/record/814455/files/1H13C_CPPASS_LHistidine.csdf"
mat_data = cp.load(filename)

# standard deviation of noise from the dataset
sigma = 0.4192854

# For the spectral fitting, we only focus on the real part of the complex dataset.
mat_data = mat_data.real

# Convert the coordinates along each dimension from Hz to ppm.
_ = [item.to("ppm", "nmr_frequency_ratio") for item in mat_data.dimensions]

When using the SSB2D method, ensure the horizontal dimension of the dataset is the isotropic dimension. Here, we apply an appropriate transpose operation to the dataset.

mat_data = mat_data.T  # transpose

# plot of the dataset.
max_amp = mat_data.max()
levels = (np.arange(24) + 1) * max_amp / 25  # contours are drawn at these levels.
options = dict(levels=levels, alpha=0.75, linewidths=0.5)  # plot options

plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_data, colors="k", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()

Create a fitting model

Guess model

Create a guess list of spin systems.

shifts = [120, 128, 135, 175, 55, 25]  # in ppm
zeta = [-70, -65, -60, -60, -10, -10]  # in ppm
eta = [0.8, 0.4, 0.9, 0.3, 0.0, 0.0]

spin_systems = single_site_system_generator(
    isotopes="13C",
    isotropic_chemical_shifts=shifts,
    shielding_symmetric={"zeta": zeta, "eta": eta},
    abundance=100 / 6,
)

Method

Create the SSB2D method.

# Get the spectral dimension parameters from the experiment.
spectral_dims = get_spectral_dimensions(mat_data)

PASS = SSB2D(
    channels=["13C"],
    magnetic_flux_density=9.395,  # in T
    rotor_frequency=1500,  # in Hz
    spectral_dimensions=spectral_dims,
    experiment=mat_data,  # add the measurement to the method.
)

# Optimize the script by pre-setting the transition pathways for each spin system from
# the method.
for sys in spin_systems:
    sys.transition_pathways = PASS.get_transition_pathways(sys)

Guess Spectrum

# Simulation
# ----------
sim = Simulator(spin_systems=spin_systems, methods=[PASS])
sim.run()

# Post Simulation Processing
# --------------------------
processor = sp.SignalProcessor(
    operations=[
        # Lorentzian convolution along the isotropic dimensions.
        sp.FFT(axis=0),
        sp.apodization.Exponential(FWHM="50 Hz"),
        sp.IFFT(axis=0),
        sp.Scale(factor=60),
    ]
)
processed_data = processor.apply_operations(data=sim.methods[0].simulation).real

# Plot of the guess Spectrum
# --------------------------
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_data, colors="k", **options)
ax.contour(processed_data, colors="r", linestyles="--", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()

Least-squares minimization with LMFIT

Use the make_LMFIT_params() for a quick setup of the fitting parameters.

params = sf.make_LMFIT_params(sim, processor)
print(params.pretty_print(columns=["value", "min", "max", "vary", "expr"]))

Out:

Name                                      Value      Min      Max     Vary     Expr
SP_0_operation_1_Exponential_FWHM            50     -inf      inf     True     None
SP_0_operation_3_Scale_factor                60     -inf      inf     True     None
sys_0_abundance                           16.67        0      100     True     None
sys_0_site_0_isotropic_chemical_shift       120     -inf      inf     True     None
sys_0_site_0_shielding_symmetric_eta        0.8        0        1     True     None
sys_0_site_0_shielding_symmetric_zeta       -70     -inf      inf     True     None
sys_1_abundance                           16.67        0      100     True     None
sys_1_site_0_isotropic_chemical_shift       128     -inf      inf     True     None
sys_1_site_0_shielding_symmetric_eta        0.4        0        1     True     None
sys_1_site_0_shielding_symmetric_zeta       -65     -inf      inf     True     None
sys_2_abundance                           16.67        0      100     True     None
sys_2_site_0_isotropic_chemical_shift       135     -inf      inf     True     None
sys_2_site_0_shielding_symmetric_eta        0.9        0        1     True     None
sys_2_site_0_shielding_symmetric_zeta       -60     -inf      inf     True     None
sys_3_abundance                           16.67        0      100     True     None
sys_3_site_0_isotropic_chemical_shift       175     -inf      inf     True     None
sys_3_site_0_shielding_symmetric_eta        0.3        0        1     True     None
sys_3_site_0_shielding_symmetric_zeta       -60     -inf      inf     True     None
sys_4_abundance                           16.67        0      100     True     None
sys_4_site_0_isotropic_chemical_shift        55     -inf      inf     True     None
sys_4_site_0_shielding_symmetric_eta          0        0        1     True     None
sys_4_site_0_shielding_symmetric_zeta       -10     -inf      inf     True     None
sys_5_abundance                           16.67        0      100    False 100-sys_0_abundance-sys_1_abundance-sys_2_abundance-sys_3_abundance-sys_4_abundance
sys_5_site_0_isotropic_chemical_shift        25     -inf      inf     True     None
sys_5_site_0_shielding_symmetric_eta          0        0        1     True     None
sys_5_site_0_shielding_symmetric_zeta       -10     -inf      inf     True     None
None

Solve the minimizer using LMFIT

minner = Minimizer(sf.LMFIT_min_function, params, fcn_args=(sim, processor, sigma))
result = minner.minimize()
report_fit(result)

Out:

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 183
    # data points      = 32768
    # variables        = 25
    chi-square         = 54771.5124
    reduced chi-square = 1.67277013
    Akaike info crit   = 16883.5043
    Bayesian info crit = 17093.4345
[[Variables]]
    sys_0_site_0_isotropic_chemical_shift:  119.168697 +/- 0.00372790 (0.00%) (init = 120)
    sys_0_site_0_shielding_symmetric_zeta: -72.1291854 +/- 0.32647311 (0.45%) (init = -70)
    sys_0_site_0_shielding_symmetric_eta:   0.98544438 +/- 0.00762029 (0.77%) (init = 0.8)
    sys_0_abundance:                        16.2139175 +/- 0.07758057 (0.48%) (init = 16.66667)
    sys_1_site_0_isotropic_chemical_shift:  128.195863 +/- 0.00312967 (0.00%) (init = 128)
    sys_1_site_0_shielding_symmetric_zeta: -75.6251701 +/- 0.27375960 (0.36%) (init = -65)
    sys_1_site_0_shielding_symmetric_eta:   0.94619261 +/- 0.00580552 (0.61%) (init = 0.4)
    sys_1_abundance:                        20.4664816 +/- 0.07809427 (0.38%) (init = 16.66667)
    sys_2_site_0_isotropic_chemical_shift:  136.194877 +/- 0.00476933 (0.00%) (init = 135)
    sys_2_site_0_shielding_symmetric_zeta: -86.3262707 +/- 0.38698036 (0.45%) (init = -60)
    sys_2_site_0_shielding_symmetric_eta:   0.42642729 +/- 0.00815945 (1.91%) (init = 0.9)
    sys_2_abundance:                        12.3222669 +/- 0.07828852 (0.64%) (init = 16.66667)
    sys_3_site_0_isotropic_chemical_shift:  172.997746 +/- 0.00307179 (0.00%) (init = 175)
    sys_3_site_0_shielding_symmetric_zeta: -69.1130708 +/- 0.25253811 (0.37%) (init = -60)
    sys_3_site_0_shielding_symmetric_eta:   0.99996947 +/- 0.00633576 (0.63%) (init = 0.3)
    sys_3_abundance:                        19.3490339 +/- 0.07587813 (0.39%) (init = 16.66667)
    sys_4_site_0_isotropic_chemical_shift:  54.5170718 +/- 0.00145926 (0.00%) (init = 55)
    sys_4_site_0_shielding_symmetric_zeta: -20.0877990 +/- 0.12450012 (0.62%) (init = -10)
    sys_4_site_0_shielding_symmetric_eta:   0.16387265 +/- 0.13564086 (82.77%) (init = 0)
    sys_4_abundance:                        18.1300502 +/- 0.05430808 (0.30%) (init = 16.66667)
    sys_5_site_0_isotropic_chemical_shift:  26.9919519 +/- 0.00163183 (0.01%) (init = 25)
    sys_5_site_0_shielding_symmetric_zeta: -10.3631453 +/- 0.54821268 (5.29%) (init = -10)
    sys_5_site_0_shielding_symmetric_eta:   0.84480483 +/- 0.21592244 (25.56%) (init = 0)
    sys_5_abundance:                        13.5182498 +/- 0.05061883 (0.37%) == '100-sys_0_abundance-sys_1_abundance-sys_2_abundance-sys_3_abundance-sys_4_abundance'
    SP_0_operation_1_Exponential_FWHM:      99.1258201 +/- 0.31273906 (0.32%) (init = 50)
    SP_0_operation_3_Scale_factor:          101.312086 +/- 0.22835689 (0.23%) (init = 60)
[[Correlations]] (unreported correlations are < 0.100)
    C(sys_5_site_0_shielding_symmetric_zeta, sys_5_site_0_shielding_symmetric_eta) =  0.940
    C(sys_4_site_0_shielding_symmetric_zeta, sys_4_site_0_shielding_symmetric_eta) =  0.682
    C(SP_0_operation_1_Exponential_FWHM, SP_0_operation_3_Scale_factor)            =  0.562
    C(sys_1_site_0_shielding_symmetric_zeta, sys_1_site_0_shielding_symmetric_eta) =  0.438
    C(sys_3_site_0_shielding_symmetric_zeta, sys_3_site_0_shielding_symmetric_eta) =  0.434
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_site_0_shielding_symmetric_eta) =  0.430
    C(sys_2_site_0_shielding_symmetric_zeta, sys_2_site_0_shielding_symmetric_eta) =  0.340
    C(sys_4_site_0_shielding_symmetric_zeta, sys_4_abundance)                      = -0.291
    C(sys_0_site_0_shielding_symmetric_zeta, sys_0_abundance)                      = -0.291
    C(sys_3_site_0_shielding_symmetric_zeta, sys_3_abundance)                      = -0.285
    C(sys_4_abundance, SP_0_operation_3_Scale_factor)                              = -0.279
    C(sys_0_abundance, sys_1_abundance)                               = -0.274
    C(sys_1_site_0_shielding_symmetric_zeta, sys_1_abundance)                      = -0.270
    C(sys_1_abundance, sys_2_abundance)                               = -0.270
    C(sys_1_abundance, sys_3_abundance)                               = -0.263
    C(sys_0_abundance, sys_3_abundance)                               = -0.257
    C(sys_2_abundance, sys_3_abundance)                               = -0.247
    C(sys_0_abundance, sys_2_abundance)                               = -0.232
    C(sys_2_site_0_shielding_symmetric_eta, sys_2_abundance)                       =  0.223
    C(sys_2_site_0_shielding_symmetric_zeta, sys_2_abundance)                      = -0.218
    C(sys_2_abundance, sys_4_abundance)                               = -0.207
    C(sys_1_site_0_isotropic_chemical_shift, sys_1_abundance)                      =  0.200
    C(sys_4_site_0_shielding_symmetric_eta, sys_4_abundance)                       =  0.191
    C(sys_0_abundance, sys_4_abundance)                               = -0.183
    C(sys_4_site_0_isotropic_chemical_shift, SP_0_operation_1_Exponential_FWHM)    = -0.163
    C(sys_2_abundance, SP_0_operation_3_Scale_factor)                              =  0.162
    C(sys_1_site_0_isotropic_chemical_shift, SP_0_operation_1_Exponential_FWHM)    = -0.160
    C(sys_3_abundance, sys_4_abundance)                               = -0.155
    C(sys_1_abundance, sys_4_abundance)                               = -0.152
    C(sys_3_site_0_shielding_symmetric_zeta, SP_0_operation_3_Scale_factor)        = -0.119
    C(sys_0_site_0_shielding_symmetric_zeta, SP_0_operation_3_Scale_factor)        = -0.117
    C(sys_1_site_0_shielding_symmetric_zeta, SP_0_operation_3_Scale_factor)        = -0.114

The best fit solution

best_fit = sf.bestfit(sim, processor)[0]

# Plot of the best fit solution
plt.figure(figsize=(8, 3.5))
ax = plt.subplot(projection="csdm")
ax.contour(mat_data, colors="k", **options)
ax.contour(best_fit, colors="r", linestyles="--", **options)
ax.set_xlim(180, 15)
plt.grid()
plt.tight_layout()
plt.show()

1

B. J. Walder, K. K. Dey, D. C. Kaseman, J. H. Baltisberger, and P. J. Grandinetti, Sideband separation experiments in NMR with phase incremented echo train acquisition, J. Phys. Chem. 2013, 138, 174203-1-12. DOI: 10.1063/1.4803142