T2 cpmg analysis
Import the necessary libraries¶
specify path to the data file and ensure that "\\" is appended to the end of the path¶
- create an instance t1 of T1Functions
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. In this case, the vclist is loaded, but vdlist terminology will be used.
- Note: the vclist contains the variable count and needs to be converted to time delay -- this will be explained in a later stage.
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)
Process the returned 1D NMR spectra¶
- apply the Gaussian apodisation (fwhm)
- zero-filling for increased digital resolution (zero_fill_factor)
- 0th order phase correction (ph0)
- 1st 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 1st 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 = "50 Hz"
ZERO_FILL_FACTOR = 10
# fallback phase used for any spectrum index you don't explicitly tune below
DEFAULT_PH0 = 292
DEFAULT_PH1 = 0.1948
phase_params = {} # spectrum index -> (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=-360, 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('$^{7}$Li chemical shift (ppm)')
ax.set_ylabel('Intensity (a.u.)')
ax.set_title(f'Spectrum {i} | PH0 = {ph0}\u00b0 PH1 = {ph1:.3f}')
#ax.set_xlim(20, -20)
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
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).
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=-15, ppm_end=15)
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
Convert vd list to numpy array and ensure that the list and extracted intensities and areas are of the same length¶
- The variable count list (n) imported from the file_path is converted to a numpy array
- The conversion factor 2 x D2 is applied to convert counts to time (seconds)
- where [D2 - π - D2]n is used for this conversion
- 4.86 μs delay i.e. D2 is used in the experiment and n corresponds to the variable count list containerised in the vdlist variable
- To ensure consistency, the simps_ints and abs_ints are sliced to match the length of vd_list, since sometimes the experiment is stopped when the NMR user observes that the system has completely relaxed
Note, D2 is labelled D20 in Bruker.
# Enter your D2 or D20 value in the variable D2
D2 = 4.86e-6
#vd_list imported from the file_path and converted into a numpy array
vd_list = np.array(vd_list)
vd_list = vd_list * 2 * D2
# slicing the vd_list if some data points are missing
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.errorbar(vd_list, simps_ints, yerr=int_uncs, fmt='o', color='red', label='error')
# 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()
Single exponential fitting¶
#T1rho fitting
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, 'Absolute Intensities', 'Guess Abs Int Curve', 'Fitted Curve Absolute Intensity', 'b')
]
# Initial guess parameters
T1_guess = 10.6 * 10**-3
for i, (ints, label, guess_label, fitted_label, color) in enumerate(intensity_sets):
if i ==0:
# A_guess = np.max(ints)
A_guess = 747197159
B_guess = 43
C_guess = np.min(ints)
# Scatter plot
# ax.scatter(vd_list, ints, color=color, label=label)
#scatter plot with marker having no face color
ax.scatter(vd_list, ints, color=color, marker='o', facecolors='none', label=label)
# Guess curve
guess_integrated_int = t1.mono_expdec(vd_list, T1_guess, A_guess, B_guess, C_guess)
# guess_integrated_int = t1.expdec(vd_list, T1_guess, T2_guess, T3_guess, A_guess, B_guess)
# ax.plot(vd_list, guess_integrated_int, color='brown', linestyle='--', label=guess_label, alpha=0.9)
# Fit the data
popt, pcov = curve_fit(t1.mono_expdec, vd_list, ints, p0=[T1_guess, A_guess, B_guess, C_guess])
# Save the fitted params and uncertainties
T1_fitted, A_fitted, B_fitted, C_fitted = popt
T1_unc, A_unc, B_unc, C_unc = np.sqrt(np.diag(pcov))
#define T1 and T2
component_1 = A_fitted * (B_fitted)*np.exp(-vd_list/T1_fitted) + C_fitted
# Extract the fitted curve
fitted_curve = t1.mono_expdec(vd_list,T1_fitted, A_fitted, B_fitted, C_fitted)
ax.plot(vd_list, fitted_curve, linestyle='-', color=color, label=fitted_label)
# ax.scatter(vd_list, fitted_curve, color='black', marker='o', facecolors='none', label=fitted_label)
ax.plot(vd_list, component_1, linestyle='--', color='black', alpha=0.5, label='component_1')
# print the fitted parameters and uncertainties
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}')
print(f'C_{label.lower().replace(" ", "_")}: {C_fitted} ± {C_unc}')
# #Format the string and append fitted parameters
output_lines.append(f'M0_{label.lower().replace(" ", "_")}: {A_fitted} ± {A_unc}\n')
output_lines.append(f'T1_{label.lower().replace(" ", "_")}: {T1_fitted} ± {T1_unc}\n')
output_lines.append(f'B_{label.lower().replace(" ", "_")}: {B_fitted} ± {B_unc}\n')
output_lines.append(f'C_{label.lower().replace(" ", "_")}: {C_fitted} ± {C_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)')
#plot the covariance matrix in another figure and label the axes with the fitted parameters
plt.savefig(filepath+'mono_exp_T1_fitting.svg', bbox_inches='tight', transparent=True)
plt.tight_layout()
fig, ax = plt.subplots()
im = ax.imshow(np.log(np.abs(pcov)))
ax.set_xticks(np.arange(len(popt)))
ax.set_yticks(np.arange(len(popt)))
ax.set_xticklabels(['T1', 'A', 'B', 'C'])
ax.set_yticklabels(['T1', 'A', 'B', 'C'])
plt.colorbar(im)
plt.show()
plt.clf()
plt.close()
From the fit above, it is evident that there is more than 1 T2 spin-spin relaxation time constant, so the relaxation curve will be fitted with two T2 components.
Multiple exponential fitting¶
#T1rho fitting
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, 'Absolute Intensities', 'Guess Abs Int Curve', 'Fitted Curve Absolute Intensity', 'b')
]
# Initial guess parameters
T1_guess = 10.6 * 10**-3
T2_guess = 0.1 * 10**-3
for i, (ints, label, guess_label, fitted_label, color) in enumerate(intensity_sets):
if i ==0:
# A_guess = np.max(ints)
A_guess = 747197159
C_guess = 0.8
D_guess = 0.2
# Scatter plot
# ax.scatter(vd_list, ints, color=color, label=label)
#scatter plot with marker having no face color
ax.scatter(vd_list, ints, color=color, marker='o', facecolors='none', label=label)
# Guess curve
guess_integrated_int = t1.di_expdec(vd_list, T1_guess, T2_guess, A_guess, C_guess, D_guess)
# guess_integrated_int = t1.expdec(vd_list, T1_guess, T2_guess, T3_guess, A_guess, B_guess)
# ax.plot(vd_list, guess_integrated_int, color='brown', linestyle='--', label=guess_label, alpha=0.9)
# Fit the data
popt, pcov = curve_fit(t1.di_expdec, vd_list, ints, p0=[T1_guess,T2_guess, A_guess, C_guess, D_guess])
# Save the fitted params and uncertainties
T1_fitted, T2_fitted, A_fitted, C_fitted, D_fitted = popt
T1_unc, T2_unc, A_unc, C_unc, D_unc = np.sqrt(np.diag(pcov))
#define T1 and T2
component_1 = A_fitted * (C_fitted)*np.exp(-vd_list/T1_fitted)
component_2 = A_fitted * (D_fitted)*np.exp(-vd_list/T2_fitted)
# Extract the fitted curve
fitted_curve = t1.di_expdec(vd_list,T1_fitted, T2_fitted, A_fitted, C_fitted, D_fitted)
ax.plot(vd_list, fitted_curve, linestyle='-', color=color, label=fitted_label)
# ax.scatter(vd_list, fitted_curve, color='black', marker='o', facecolors='none', label=fitted_label)
ax.plot(vd_list, component_1, linestyle='--', color='black', alpha=0.5, label='component_1')
ax.plot(vd_list, component_2, linestyle='--', color='blue', alpha=0.5, label='component_2')
# print the fitted parameters and uncertainties
print(f'T1_{label.lower().replace(" ", "_")}: {T1_fitted} ± {T1_unc}')
print(f'T2_{label.lower().replace(" ", "_")}: {T2_fitted} ± {T2_unc}')
print(f'A_{label.lower().replace(" ", "_")}: {A_fitted} ± {A_unc}')
print(f'D_{label.lower().replace(" ", "_")}: {D_fitted} ± {D_unc}')
print(f'C_{label.lower().replace(" ", "_")}: {C_fitted} ± {C_unc}')
# #Format the string and append fitted parameters
output_lines.append(f'M0_{label.lower().replace(" ", "_")}: {A_fitted} ± {A_unc}\n')
output_lines.append(f'T1_{label.lower().replace(" ", "_")}: {T1_fitted} ± {T1_unc}\n')
output_lines.append(f'T2_{label.lower().replace(" ", "_")}: {T2_fitted} ± {T2_unc}\n')
output_lines.append(f'D_{label.lower().replace(" ", "_")}: {D_fitted} ± {D_unc}\n')
output_lines.append(f'C_{label.lower().replace(" ", "_")}: {C_fitted} ± {C_unc}\n')
#save the fitted params and uncertainties in a text file
with open(filepath+'di_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)')
#plot the covariance matrix in another figure and label the axes with the fitted parameters
plt.savefig(filepath+'di_exp_T1_fitting.svg', bbox_inches='tight', transparent=True)
plt.tight_layout()
fig, ax = plt.subplots()
im = ax.imshow(np.log(np.abs(pcov)))
ax.set_xticks(np.arange(len(popt)))
ax.set_yticks(np.arange(len(popt)))
ax.set_xticklabels(['T1', 'T2', 'A', 'C', 'D'])
ax.set_yticklabels(['T1', 'T2', 'A', 'C', 'D'])
plt.colorbar(im)
plt.show()
plt.clf()
plt.close()
The nucleus studied here is 7Li and the diffusing Li species have almost similar relaxation times, which is evident from the T2 relaxation time constants and the covariance matrix plot.