T1 analysis

Import the necessary libraries¶

import sys

from scipy.optimize import curve_fit

sys.path.append("../../src")

from relaxometrynmr.core import T1Functions

import numpy as np

import matplotlib.pyplot as plt

import ipywidgets as widgets

from IPython.display import display

specify path to the data file and ensure that "\\" is appended to the end of the path¶

create an instance t1 of T1Functions

filepath = r"../../data/T1_data/1//"

t1 = T1Functions(filepath)

Read and convert Bruker NMR data to NMRPipe and CSDM formats: read_and_convert_bruker_data.¶

The function automatically detects and loads the variable delay list (vdlist, vplist, vclist) used in the experiment.

It returns a tuple containing three elements: a list of 1D NMR (spectra), the variable delay list (vd_list), and the complete dataset in CSDM format (csdm_ds)

spectra, vd_list, csdm_ds = t1.read_and_convert_bruker_data(filepath)

Process the returned 1D NMR spectra¶

apply the Gaussian apodisation (fwhm)
zero-filling for increased digital resolution (zero_fill_factor)
0^th order phase correction (ph0)
1^st order phase correction (ph1) -- this phase correction is a bit nuanced and so far, a value of 0 - 0.6 ° has worked quite well: see the "Understanding Phasing" example under User Guide
In applying the 1^st order correction, you would have to experiment with the mentioned values to obtain a pure absorption line-shape signal

Interactively tune the phase for each spectrum¶

Instead of hard-coding ph0/ph1 per spectrum index in an if/else block, pick a spectrum with the dropdown and drag the sliders until the lineshape is pure absorption. Each change is saved into phase_params, keyed by spectrum index. Revisiting an index restores the last value you set for it. Any index you never touch falls back to DEFAULT_PH0 / DEFAULT_PH1 when exp_spectra is built below.

FWHM = "500 Hz"
ZERO_FILL_FACTOR = 10

# fallback phase used for any spectrum index you don't explicitly tune below
DEFAULT_PH0 = 50
DEFAULT_PH1 = 0.5

phase_params = {}  # spectrum index -&gt; (ph0, ph1), filled in as you tune each spectrum

spectrum_selector = widgets.Dropdown(
    options=list(range(len(spectra))), value=0,
    description='Spectrum idx', style={'description_width': 'initial'}
)
ph0_slider = widgets.IntSlider(
    value=DEFAULT_PH0, min=0, max=360, step=1,
    description='PH0 (\u00b0)', continuous_update=False,
    style={'description_width': 'initial'}, layout=widgets.Layout(width='500px')
)
ph1_slider = widgets.FloatSlider(
    value=DEFAULT_PH1, min=0.0, max=1.0, step=0.0005,
    description='PH1', continuous_update=False, readout_format='.3f',
    style={'description_width': 'initial'}, layout=widgets.Layout(width='500px')
)
out = widgets.Output()

def update_preview():
    i = spectrum_selector.value
    ph0, ph1 = phase_params.get(i, (DEFAULT_PH0, DEFAULT_PH1))
    exp = t1.process_spectrum(
        spectra[i], fwhm=FWHM, zero_fill_factor=ZERO_FILL_FACTOR, ph0=ph0, ph1=ph1
    )
    ppm = exp.dimensions[0].coordinates.value
    y = exp.dependent_variables[0].components[0].real
    with out:
        out.clear_output(wait=True)
        fig, ax = plt.subplots(figsize=(7, 3))
        ax.plot(ppm, y)
        ax.invert_xaxis()
        ax.set_xlabel('$^{17}$O chemical shift (ppm)')
        ax.set_ylabel('Intensity (a.u.)')
        ax.set_title(f'Spectrum {i}  |  PH0 = {ph0}\u00b0   PH1 = {ph1:.3f}')
        ax.set_xlim(615, 530)
        plt.tight_layout()
        plt.show()

def on_index_change(change):
    i = spectrum_selector.value
    ph0, ph1 = phase_params.get(i, (DEFAULT_PH0, DEFAULT_PH1))
    ph0_slider.unobserve(on_phase_change, names='value')
    ph1_slider.unobserve(on_phase_change, names='value')
    ph0_slider.value = ph0
    ph1_slider.value = ph1
    ph0_slider.observe(on_phase_change, names='value')
    ph1_slider.observe(on_phase_change, names='value')
    update_preview()

def on_phase_change(change):
    phase_params[spectrum_selector.value] = (ph0_slider.value, ph1_slider.value)
    update_preview()


spectrum_selector.observe(on_index_change, names='value')
ph0_slider.observe(on_phase_change, names='value')
ph1_slider.observe(on_phase_change, names='value')

display(widgets.VBox([spectrum_selector, ph0_slider, ph1_slider, out]))
on_index_change(None)  # draw initial preview

VBox(children=(Dropdown(description='Spectrum idx', options=(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14)…

Build `exp_spectra` from the tuned phases¶

Applies the (ph0, ph1) you set per index above (falling back to the default for any index left untouched).

exp_spectra = []
for i, spectrum in enumerate(spectra):
    ph0, ph1 = phase_params.get(i, (DEFAULT_PH0, DEFAULT_PH1))
    exp_spectrum = t1.process_spectrum(
        spectrum, fwhm=FWHM, zero_fill_factor=ZERO_FILL_FACTOR, ph0=ph0, ph1=ph1
    )
    exp_spectra.append(exp_spectrum)

Find the area under the peak of interest using the `integrate_spectrum_region()` function¶

The integration function employed here integrate each spectrum using trapezoid and simpson function, respectively. ppm_start and ppm_end need to be defined as the starting and ending ppm region needed to be integrated. The integrated area of each spectrum is appended to trapz_ints and simps_ints, respectively. x_ and y_regions are regions of integration in the spectra -- needed for visuals.

There is no difference between trapz and simps, so you would have to use either of the two in a later stage

trapz_ints = []
simps_ints = []
x_regions = []
y_regions = []
int_uncs = []

for i, exp_spectrum in enumerate(exp_spectra):
    trapz_int, simps_int, x_region, y_region, int_unc = t1.integrate_spectrum_region(exp_spectrum, ppm_start=500, ppm_end=650)
    trapz_ints.append(trapz_int)
    simps_ints.append(simps_int)
    x_regions.append(x_region)
    y_regions.append(y_region)
    int_uncs.append(int_unc)

plot_spectra_and_zoom() function¶

creates plots of NMR spectra with both full view and zoomed regions (max and min x zoom)
highlights the integrated x_ and y_regions on the zoomed plot
returns maximum intensities from each spectrum (abs_ints): relevant for relaxometry just like integrated areas contained in trapz_ints and simps_ints

max_x_zoom = 1000

min_x_zoom = 200

abs_ints = t1.plot_spectra_and_zoomed_regions(exp_spectra, x_regions, y_regions, max_x_zoom, min_x_zoom)

No description has been provided for this image

Convert vd list to numpy array and ensure that the list and extracted intensities and areas are of the same length¶

# vd_list imported from the file_path is converted into a numpy array

vd_list = np.array(vd_list)

# slicing the vd_list if some data points are missing

vd_list = vd_list[:len(abs_ints)]


simps_ints = simps_ints[:len(vd_list)]

abs_ints = abs_ints[:len(vd_list)]

Viusalise the list and the extracted intensities and areas¶

the extracted areas from either trapz or simps integration and extracted max intensities of each spectrum are plotted against corresponding time in vd_list

fig, ax = plt.subplots()

ax.scatter(vd_list, abs_ints, color='blue', label='intensity extraction')


ax.scatter(vd_list, simps_ints, color='red', label='area extraction')


ax.semilogy()

ax.semilogx()

ax.legend(loc='best', frameon=True)

ax.set_xlabel(r'$\tau$ (s)')
ax.set_ylabel('Intensity (arbitrary unit)')



plt.tight_layout()

plt.show()

Let us fit one of the data with a simple exponential satrec function: here we are using the intensities, the areas can also be utilised¶

fig, ax = plt.subplots()

output_lines = []

# Define a list of tuples for the two sets of intensities
intensity_sets = [
    (simps_ints, 'Simps Area', 'Guess Curve', 'Fitted Curve', 'r'),
    (abs_ints, 'Intensities', 'Guess int. Curve', 'Fitted int. curve', 'b')
]


# Initial guess parameters
M0_guess = 0.9
T1_guess = 2.8e-3

for i, (ints, label, guess_label, fitted_label, color) in enumerate(intensity_sets):
    # here we are using the maximum intensity of each spectrum to 
    # extract the spin-lattice relaxation time, you can use the integrated area of each spectum as well by setting i == 0
    if i == 1:
        A_guess = np.max(ints)
        B_guess = np.min(ints)

        # Scatter plot
        ax.scatter(vd_list, ints, color=color, label=label)

        # Guess curve
        guess_integrated_int = t1.mono_satrec_func(vd_list, M0_guess, T1_guess, A_guess, B_guess)
        ax.plot(vd_list, guess_integrated_int, color='brown', linestyle='--', label=guess_label, alpha=0.5)

        # Fit the data
        popt, pcov = curve_fit(t1.mono_satrec_func, vd_list, ints, p0=[M0_guess, T1_guess, A_guess, B_guess])

        # Save the fitted params and uncertainties
        M0_fitted, T1_fitted, A_fitted, B_fitted = popt
        M0_unc, T1_unc, A_unc, B_unc = np.sqrt(np.diag(pcov))

        # Extract the fitted curve
        fitted_curve = t1.mono_satrec_func(vd_list, M0_fitted, T1_fitted, A_fitted, B_fitted)
        ax.plot(vd_list, fitted_curve, linestyle='-', color='grey', label=fitted_label)

        # Print the fitted parameters and uncertainties
        print(f'M0_{label.lower().replace(" ", "_")}: {M0_fitted} ± {M0_unc}')
        print(f'T1_{label.lower().replace(" ", "_")}: {T1_fitted} ± {T1_unc}')
        print(f'A_{label.lower().replace(" ", "_")}: {A_fitted} ± {A_unc}')
        print(f'B_{label.lower().replace(" ", "_")}: {B_fitted} ± {B_unc}')

        #Format the string and append fitted parameters

        output_lines.append(f'M0_{label.lower().replace(" ", "_")}: {M0_fitted} ± {M0_unc}\n')
        output_lines.append(f'T1_{label.lower().replace(" ", "_")}: {T1_fitted} ± {T1_unc}\n')
        output_lines.append(f'A_{label.lower().replace(" ", "_")}: {A_fitted} ± {A_unc}\n')
        output_lines.append(f'B_{label.lower().replace(" ", "_")}: {B_fitted} ± {B_unc}\n')
        #save the fitted params and uncertainties in a text file
        with open(filepath+'mono_exp_fitted_params.txt', 'w') as f:
            f.writelines(output_lines)

ax.semilogy()
ax.semilogx()
ax.legend(loc='best', frameon=False)
ax.set_xlabel(r'$\tau$ (s)')
ax.set_ylabel('Intensity (arbitrary unit)')
plt.tight_layout()
plt.savefig(filepath+'mono_exp_T1_fitting.svg', bbox_inches='tight', transparent=True)
plt.show()
plt.clf()
plt.close()

M0_intensities: 0.8505179934087478 ± 53790196376120.46
T1_intensities: 0.004649182967112021 ± 0.00015951628755413104
A_intensities: 172358414460.14325 ± 1.0900642940815778e+25
B_intensities: 32979920734.106155 ± 933572502.0359585

Understanding the optimised parameters¶

M0_intensities == equilibrium magnetisation
T1_intensities == spin-lattice relaxation time in seconds
B_intensities == Baseline offset
A_intensities == scaling factor for the overall amplitude