Multiquantum Chemical Exchange Saturation Transfer NMR to Quantify Symmetrical Exchange: Application to Rotational Dynamics of the Guanidinium Group in Arginine Side Chains

Chemical exchange saturation transfer (CEST) NMR experiments have emerged as a powerful tool for characterizing dynamics in proteins. We show here that the CEST approach can be extended to systems with symmetrical exchange, where the NMR signals of all exchanging species are severely broadened. To achieve this, multiquantum CEST (MQ-CEST) is introduced, where the CEST pulse is applied to a longitudinal multispin order density element and the CEST profiles are encoded onto nonbroadened nuclei. The MQ-CEST approach is demonstrated on the restricted rotation of guanidinium groups in arginine residues within proteins. These groups and their dynamics are essential for many enzymes and for noncovalent interactions through the formation of hydrogen bonds, salt-bridges, and π-stacking interactions, and their rate of rotation is highly indicative of the extent of interactions formed. The MQ-CEST method is successfully applied to guanidinium groups in the 19 kDa L99A mutant of T4 lysozyme.


Sample Preparation.
A sample of L99A T4 Lysozyme was prepared according to previously published protocols 1 .

NMR Spectroscopy.
For all experiments, the spectrometer temperature was calibrated using a sample of d4-methanol and following standard protocols.
Longitudinal EXSY 3 zz-exchange experiments on free [ 13 C6, 15 N4]-L-arginine were carried out and analysed as described previously 1 . All zz-exchange NMR spectra were recorded at a static magnetic field of 14.1 T on a Bruker Avance II spectrometer and using a roomtemperature TXI HCN inverse probe. In all of these experiments, 20 mixing times between 2 ms and 300 ms were used, including two repeats for error estimation. The carrier frequencies were set to 4.7 ppm ( 1 H) and 70.3 ppm ( 15 N). The sweep width in the indirect dimension was set to 16.5 ppm with 40 complex points recorded. Four scans were recorded per transient and the recycle recovery delay (d1) was 1 s, giving a total recording time of ~2 hours per experiment.
The presented MQ-CEST profiles were all recorded using the pulse sequence in Fig 1D (with 13 C detection). The parameters employed in the MQ-CEST experiments are detailed in Table S1. The B1 field strengths and inhomogeneity were calibrated for each radio-frequency probe using the approach of Guenneugues et al 4 .

Data analysis.
Simulation and Fitting of MQ-CEST Data. For the simulation and fitting of MQ-CEST data in arginine guanidinium groups the three-spin system, C !" # N !$ % ! N !$ % " , is considered. The composite decoupling scheme ensures that the H ! % nuclei have no significant effect on this spin system during the CEST period, although these protons are an important source of autorelaxation. All combinations of direct products of the individual spin product operators in C !" & N !$ % ! N !$ % " and the identity operator give rise to 64 product operators 5 required to fully describe this spin system. To account for the rotational exchange about the C ζ -N ε bond S3 (symmetrical exchange of N η1 and N η2 ) a further 63 product operators are required: for convenience (in common with literature conventions) we can refer to the two states as the ground (G) and excited (E) state although the symmetric exchange here means that these states have equal energy and population. This gives a total of 127 product operators, describing the full spin system. In order to fit or simulate MQ-CEST data the Bloch-McConnell equation for the system 6,7 needs to be solved, ' '( Here, m(t) is an (n×1) column vector containing the n product operator terms for the spin system and L is the (n×n) Liouvillian matrix that describes its dynamics and time-dependence.
To calculate the Liouvillian matrix we first write down the Hamiltonian for each spin interaction: in the present case the chemical shift of the two 15 N η nuclei, the CEST pulse and Jcouplings between the 15 N η nuclei and the 13 C ζ nucleus are included. The small scalar-coupling between the two 15 N η is ignored 8 . For each Hamiltonian, -, the commutator superoperator, --, is calculated as: In which 2 is the identity matrix and the tilde denotes the complex conjugate. The elements of L can then be found by transforming -to the product operator basis 9 .
Having added these interactions to the Liouvillian, relaxation terms are added to the matrix phenomenologically, following the approach of Helgstrand et al 10 . Due to the small gyromagnetic ratios of 13 C and 15 N dipolar interactions between them will be a relatively small source of relaxation 11 , which are ignored. Instead, the auto-relaxations will be dominated by dipolar relaxation to proximal 1 H nuclei and the chemical shift anisotropy of the individual spins. Consequently, to a good approximation, we can write the spin relaxation of higher order product operators as the sum of the individual longitudinal and transverse single order spin operators: Where , , ∈ ( , , ). It is also assumed that * % & ! = * % & " . This approach considerably reduces the number of parameters to be optimised in the model and, as shown below, the derived kex and Δω values are relatively uncorrelated from the transverse relaxation rates ( Fig   S8). It should however be noted that the MQ-CEST approach should not be treated as a method of attaining relaxation rates of the various product operator terms. It should also be noted that the 15 N η nuclei will undergo multiexponential relaxation, which we approximate here with a single exponential, as a result of cross-correlated dipolar relaxation with each of the two directly bonded 1 H η nuclei 12 .
To complete the Liouvillian matrix, L, chemical exchange between the two 15 N η nuclei needs to be added. The 127×127 exchange matrix k is defined as: where kGE and kEG are the forward and backwards rates of rotation, respectively, kex = kGE + kEG and -" is the 63×63 identity matrix. Since the exchange between 15 N η1 and 15 N η2 is symmetric, kEG = kGE. Finally, to obtain the full Liouvillian, L, k is added to Eq 2.
To reduce the size of the Liouvillian matrix used to propagate the spin-system, it is noted that, in the MQ-CEST experiment, transverse 13 C ζ magnetisation cannot be created and so any terms containing these are ignored. Removing these terms, we arrive at a 63×63 matrix, L, that describes the evolution of the three-spin . / ! / " spin system during the MQ-CEST experiment.
Once L is known and since L is time-independent, Eq. 1 can be solved according to, where Eq. 5 in turn can be solved by directly taking the matrix exponential of L, using the Padé approximation (as implemented in the python module scipy 13 ) or similar. For analysis of the MQ-CEST data, the matrix exponential can be calculated more rapidly through a Taylor

S5
When simulating MQ-CEST data it is important to take B1 field inhomogeneity into account. To do this, for each CEST offset frequency, the CEST intensity is calculated for ten B1 fields evenly spaced between ±2σ of the average CEST pulse strength 15 (in simulations σ is set to 10% of the CEST pulse strength while for fitting experimentally determined values are used) and the final intensity result is calculated as the sum of the results weighted by Gaussian coefficients.
For fitting experimental data, the results of Eq. 5 need to be passed to an optimiser to minimise the standard χ 2 equation 16 : Here the sum extends over all experimental data points per residue, including data at multiple CEST field strengths and/or multiple static magnetic field strengths, and x refers to the parameters to be optimised, including kex, Δω, average 15 N η chemical shifts and relaxation rates.
The uncertainties, @AB,> are estimated from the scatter in the CEST profiles 17 . A python script, which makes extensive use of the numpy and scipy modules 13 has been written for the simulation and fitting of MQ-CEST data. The least-squares optimiser is implemented using the LMFIT module 18 . This script is available upon request.
Analysis of longitudinal zz-exchange spectra. For each temperature the intensities of the cross and diagonal peaks for the N η nuclei are quantified using FuDA 19 . They are then fit to the following equation with kex as a fitting parameter 1 : Where Iab and Iba are the cross-peak intensities and Iaa and Ibb are the intensities of the two diagonal peaks.

S6
Supporting Table   Table S1: Acquisition parameters used in MQ-CEST NMR Experiments (all experiments were carried out using the 13 C detect pulse sequence in Fig 1D) a) TXI: HCN inverse room-temperature probe. TCI: Cryogenic HNC inverse probe with cooled 1 H and 13 C preamplifiers, 1 H inner coil and 13 C outer coil. TXO: 13 C-optimised cryogenic CNH probe with cooled 1 H, 13   a) In all cases errors are taken from the covariance matrix of the fits 16 . These errors are cross validated by performing fits using grid searched kex values (see Fig S8).  All gradients are 1 ms in length and applied with the following strengths: g1 = 12.3 G/cm, g2 = 5.9 G/cm, g3 = 16.6 G/cm, g4 = 10.2 G/cm, g5 = 7.0 G/cm, g6 = 15.5 G/cm, g7 = 9.1 G/cm and g8 = 19.8 G/cm.

Fig S5.
Comparison of signal/noise ratios for 13 C ( Fig S4A) and 1 H-detected ( Fig S4B) reference planes of the MQ-CEST experiment. The signal/noise are calculated for R52 in L99A T4 Lysozyme. The comparison here ignores the role of hydrogen-exchange with the solvent (very slow under the conditions employed here for R52 25 ). The spectra were acquired with the same number of scans, complex points (and acquisition time) in t1 and recycle delay on a 700 MHz spectrometer equipped with a TCI cryogenic probe and processed identically. The temperature of the sample was altered in steps of 5 K between 283 K and 298 K to change the correlation time. The global correlation time, ! , is estimated at each temperature using an empirical relationship 26 . The local correlation time is then estimated from this as " ! with the value for the 15 N-1 H order parameters S 2 taken from Werbeck et al 27 . The measured signal to noise ratios were plotted against the estimated local correlation times above. The dotted lines indicate fits to a decaying exponential function. This fit is extrapolated between 1 and 25 ns. Overall, 13 C detection is preferable when the local correlation time of a residue exceeds approximately 10 ns. On a 13 C optimised TXO cryogenic probe, 13 C detection is preferable for residues with even smaller local correlation times. Of note is that flexible arginine side chains, with small local correlation times, are often severely affected by hydrogen-exchange.   The colors indicate different confidence intervals calculated by taking the difference between the minimum " value and the " value at each point in the grid and then drawing contour lines for Δ " at the indicated confidence intervals for two degrees of freedom 16 . Overall, the fitted exchange rate, kex, is only very weakly correlated with the transverse relaxation rates, R2. This means that accurate kex are obtained although the obtained R2 values may be uncertain. The horizontal dashed lines represent the confidence interval for the experimentally measured R2 values ( R2 ± σR2 ). refocus and evolve respective J-couplings ; -Apply 80 us WALTZ64 decoupling on 1H to suppress J-coupling revolution ; -At the end of this period we will have 2Cz(zeta)Nz(eps) magnetisation -note that sweep widths may need ; to be adjusted slightly to avoid negative delays. ; ; Selective 13Czeta Excitation ; -1.5 ms Eburp2 pulse at 16.4 T ; ; CEST period ; -The saturation frequencies are taken from fq1list. If the absolute frequency value is set to >10000Hz a ; references plane is recorded ; -The CEST pulse length is set by d17. ;preprocessor-flags-start ;HALFDWELL: for initial sampling delay of half a dwell-time with refocus and evolve respective J-couplings ; -Apply 80 us WALTZ64 decoupling on 1H to suppress J-coupling revolution ; -At the end of this period we will have 2Cz(zeta)Nz(eps) magnetisation -note that sweep widths may need ; to be adjusted slightly to avoid negative delays. ; ; Selective 13Czeta Excitation ; -1.5 ms Eburp2 pulse at 16.4 T ; ; CEST period ; -The saturation frequencies are taken from fq1list. If the absolute frequency value is set to >10000Hz a ; references plane is recorded ; -Set cnst12 = 8.0 (move 1H carrier to middle of amide region) ; -The CEST pulse length is set by d17. ;preprocessor-flags-start ;HALFDWELL: for initial sampling delay of half a dwell-time with ; option -DHALFDWELL (eda: ZGOPTNS) ;preprocessor-flags-end