Super-Resolved Nuclear Magnetic Resonance Spectroscopy

We present a novel method that breaks the resolution barrier in nuclear magnetic resonance (NMR) spectroscopy, allowing one to accurately estimate the chemical shift values of highly overlapping or broadened peaks. This problem is routinely encountered in NMR when peaks have large linewidths due to rapidly decaying signals, hindering its application. We address this problem based on the notion of finite-rate-of-innovation (FRI) sampling, which is based on the premise that signals such as the NMR signal, can be accurately reconstructed using fewer measurements than that required by existing approaches. The FRI approach leads to super-resolution, beyond the limits of contemporary NMR techniques. Using this method, we could measure for the first time small changes in chemical shifts during the formation of a Gold nanorod-protein complex, facilitating the quantification of the strength of such interactions. The method thus opens up new possibilities for the application and acceleration of multidimensional NMR spectroscopy across a wide range of systems.


S1.2 Parameter Estimation of FIDs
A typical one-dimensional free induction decay (FID) could be modeled in the form of a sum of weighted exponentials (SWE), as follows: where r denotes the largest integer less than or equal to r. For uniqueness of the estimated ω s from f (nT s ), the sampling interval T s should satisfy the inequality T s < 2π ω ,max , where ω ,max denotes the maximum frequency ω s could achieve.
Among the three unknown parameters of FID signals mentioned above, estimation of ω s is most important. A standard technique to determine ω s is to first estimate the Fourier spectrum from f (nT s ) using the discrete-time Fourier transform (DTFT). The frequencies are estimated as the locations of the peaks in the spectrum. The sampled DTFT is computed efficiently using the fast Fourier transform S3 algorithm. The Fourier approach has two major drawbacks. The first one stems from the fact that due to a finite observation window and non-zero values of α s, the signal f (t) is not bandlimited. Hence, sampling f (t) at any finite rate f s = 1 Ts will result in aliasing 1,2 , which affects the locations of the peaks in the spectra. The Fourier transform of f (t) is given as a e s Tobs−jω − 1 s − jΩ , (S. 3) and it can be shown that aliasing reduces as f s increases. For T s < 2π ω ,max , the aliasing error depends on T obs and the minimum value of α . The other limitation of the Fourier approach is its limited frequency resolution, which is determined by the uncertainty principle. For α = 0, the DTFT of f (nT s ) cannot resolve between two ω s separated apart by less than 2π Tobs . In order to resolve closely spaced frequencies, one has to observe the FID for a longer duration, which requires a longer experimentation time. Also, with longer durations, typically, the FIDs decay and noise begins to dominate making the estimates inaccurate and unreliable.
Fourier-based FID analysis does not take into account the SWE model of FIDs. By using this a priori signal model, several alternatives have been developed to obtain a resolution higher than that of the Fourier approach for a given number of measurements. A popular approach in the NMR community is linear prediction (LP) [3][4][5][6][7] , where each sample f (nT s ) is related to past L contiguous samples by means of a linear combination, that is, there exist L complex numbers c s such that The prediction coefficients {c } depend on {−α + jω } L =1 . In practice, measurement noise would make the linear combination on the right-hand side of (S.4) an estimate of the left-hand side. By using the linear predictability of f (nT s ), Prony 8 showed that {a , s } L =1 could be computed from 2L contiguous values of f (nT s ). Based on Prony's analysis, several LP-based methods have been proposed in both NMR and signal processing literature 3,4,[8][9][10][11][12][13][14][15] . There are several variants that use the SWE structure of FID to improve spectral resolution, for example, entropy-based methods [16][17][18] , Bayesian methods [19][20][21] , maximumlikelihood methods 22,23 , filter-diagonalization methods [24][25][26] , Padé-Laplace analysis 27,28 , etc. We briefly discuss the standard LP technique before proceeding with the development of our technique.

S1.3 Least-Squares Approach with Polynomial Root-Finding
The linear prediction equation in (S.4) could be expressed in matrix form as . . .
For N ≥ 2L, the matrix F has rank L and one can uniquely determine c by solving (S.5) using least-squares regression, and the least-squares solution isĉ LS = (F H F) −1 F H f . Alternatively, structured decompositions of F such as Cholesky decomposition 29 , Householder QR decomposition 15 , and singularvalue decomposition (SVD) 12 , could also be deployed to efficiently solve (S.5). Once the c s are obtained, one could compute the roots of the polynomial p(z) = L =0 c z , where c 0 = 1, and the roots will be . This idea forms the basis for several techniques 8,10 .

S1.4 LP in the Presence of Noise and Model-Order Selection
In practice, the measurements {f (nT s )} are corrupted by noise. The noisy samples are modeled as a e s nTs + w(n), (S. 6) where w(n)s are assumed to be identically distributed, zero-mean, circular white Gaussian random variables 30 with variance σ 2 w . It is also assumed that w(n)s are independent of the signal f (nT s ). The noisy samples {f (nT s )} are not linearly predictable, and the prediction coefficients are estimated as the solution to the constrained minimization problem: where c a ∈ C L+1 and 2 denotes the Euclidean norm. The solution to the above minimization problem is computed by first performing the SVD of [f F] and then choosing the right singular vector corresponding to the minimum singular value 31 . This method is effectively linear prediction using SVD (LPSVD) 12 .
Several variants of this basic idea are available in the NMR literature.
Another important problem in both NMR and signal processing communities is to estimate the unknown , τ 2 is termed as the rate of innovation of y(t). When R y is finite, the signal y(t) is referred to as an FRI signal. Vetterli et al. showed that FRI signals, such as stream of Dirac impulses or differential Dirac impulses, piecewise splines, etc., could be reconstructed from samples taken at the rate of innovation using Gaussian or sinc sampling kernels. The reconstruction mechanism is based on high-resolution spectral estimation techniques 30 Hence, the rate of innovation is finite and equal to 2L T obs . Due to the FRI property, FID signals could be reconstructed using a minimum of 2L samples.
Since the signal structure is known, super-resolution is possible even with a finite number of samples.

S1.6 DEESPRIT -A New Technique for Handling Decaying Exponentials
In the section, we develop a new autocorrelation-based technique for parameter estimation from the FID samples shown in (S.2). In its standard manifestation, ESPRIT is capable of handling only complex sinusoids. On the other hand, FIDs are often decaying and hence, the standard ESPRIT is not directly applicable. We need a counterpart of ESPRIT that is capable of handling decaying exponentials -this is precisely one of our contributions in this paper. The difference between our technique and the other techniques in the signal processing literature that can handle damped exponentials 10,14 is that ours is S6 based on the autocorrelation. Before going into the details of the proposed method, we briefly discuss standard ESPRIT.

1) ESPRIT for frequency estimation:
Consider the problem of estimating the frequencies {ω } L =1 of the complex exponential signalf (n) in noise, given as where ω ∈ (−π, π] and a s are independent Gaussian random variables with respective variances σ 2 a , and w(n)s are independent of f (n) and i.i.d. circular Gaussian random variables with zero mean and variance , for n, m ∈ 0, N − 1 , and , p ∈ 1, L , where δ K denotes the Kronecker delta. The expectation is computed jointly over the amplitudes and noise, which are also mutually independent.
The autocorrelation sequence of the random processf (n) in (S.7) is given as the autocorrelation sequence rff (n, m) is written as By using (S.8), the M th -order autocorrelation matrix off (n) is decomposed as A ∈ R L×L is a diagonal matrix with entries A[ , ] = σ 2 a , and I M is the M × M identity matrix. Consider the matrix U ∈ C M ×L whose columns consist of eigenvectors corresponding to the L largest eigenvalues of R yy . It can be shown that, for M ≥ L, the matrices U and R yy have the same range S7 space 30 . Hence, for M ≥ L, there exists a nonsingular matrix P ∈ C L×L such that U = V P. By using the Vandermonde structure of V, one can write where Z ∈ C L×L is a diagonal matrix with the elements Z[ , ] = e −jω . By using the relation U = V P and (S.11), one can write Since P −1 ZP has same eigenvalues as Z, which are {e −jω } L =1 , ω s can be computed from the eigenvalues of U † F U L , which, in turn, can be computed from the eigenvalues of the autocorrelation matrix Rff .
2) DEESPRIT for parameter estimation of damped exponentials: Consider the samplesf (n) in (S. 6) for n ∈ 0, N − 1 . We assume that a s are independent Gaussian random variables with respective variances σ 2 a . The correlation between the random variables f (nT s ) = a z n and f p (mT s ) = a p z m p is given as where z = e −(α −jω )Ts , and the expectation is computed over the amplitudes. By using (S.13), the autocorrelation sequence of the random processf (nT s ) in (S.6) is written as where the expectation is taken jointly over the amplitudes and noise. Since |z | is not necessarily equal to unity, rff (nT s , mT s ) cannot be written as a function of only n − m. This is a consequence of the nonstationarity of the random process f (nT s ) in (S.6) unlike the process in (S.7), which is stationary.
Due to the nonstationarity, we cannot directly construct a matrix of the form given in (S.10).
Consider f (nT s ) in (S.6), and the M th -order autocorrelation matrix Rff (nT s ) = E{f n,Mf To remove the dependency of the matrix Rff (nT s )[k 1 , k 2 ] on n, we construct the matrix R ã ff , which is Rff (nT s ). The elements of R ã ff are given as The matrix R ã ff has a form similar to that of Rff in (S.10). Now, one could apply the matrix decompositions developed for classical ESPRIT to R ã ff . In practice, the autocorrelation matrix R ã ff is estimated from the noisy samples {f (n)} N −1 n=0 as where

S1.8 Parameter Estimation Using FRI-NMR -A Simulation
The first simulation is designed to demonstrate the parameter estimation capability of the FRI-NMR

S1.9 Accuracy of Frequency Estimation in FRI-NMR as Function of Damping and Amplitudes
In this section, we analyze experimentally the accuracy of frequency estimation (measured in Hz) as a function of the damping factor α and amplitude a. We consider 256 samples of a single-component FID with frequency and sampling rate of 500 Hz and 1500 Hz, respectively. The signal is corrupted by zeromean, additive white Gaussian noise (AWGN) with a fixed standard deviation of 2.5. The frequency is estimated by FRI-NMR as the amplitude is increased from 1 to 5 in steps of 1 and for α = 1, 2, · · · , 9, 10, and 20. For a given amplitude and damping, we computed average S/N and root mean-squared error in frequency estimation over 500 independent noise realizations.
The computed average S/N is shown in Figure S2(a) and RMSE in Figure S2

S2. COMPARISON OF FRI-NMR WITH LPSVD AND FOURIER METHODS
The goal of this experiment is to asses the resolution capability of the FRI-NMR method and compare it with that of the LPSVD (cf. Section S1) and Fourier methods. We considered N = 600 noisy FID samples (cf. (S.6)) for L = 2 with a 1 = a 2 = 1, α 1 = α 2 = 3. The frequencies of the two components are F 1 = 5000 Hz and F 2 = 5000 + ∆F Hz, where ∆F represents the frequency separation. The FID was sampled at 12 kHz and corrupted by zero mean additive white Gaussian noise (AWGN) to simulate a noisy measurement, with a variance σ 2 w chosen to achieve a given S/N. The samples are apodized by using cos(30) window function. To compute the Fourier spectra, the samples are zero-padded to have a length of 2 17 . We chose similar amplitudes and damping factors for the two FID components such that they have equal S/N. The performance of the technique in resolving the two frequencies is assessed by varying ∆F below the resolution limit given by 1 Tobs ≤ 1 (N −1)Ts = 20.034 Hz. We observed that for a given S/N, ∆F , and noise realization, the LPSVD and FRI-NMR methods may not always estimate two distinct frequencies. Based on the estimatesF 1 andF 2 , we classify the outcome of a trial for a given method as belonging to one of the following three categories: (i) Resolved  Figure S3 shows bar plots that indicate the number of times a given method resolves, misses, or does not resolve the frequencies, out of a total of 1000 trials. The classification metric shown in these bar plots show the resolution capability but does not quantify the error in the estimation. To quantify the accuracy of the methods, for a given S/N and ∆F , we compute normalized bias, root-mean variance, and root-mean-squared error (RMSE) in Hz. These three metrics are computed as respectively, whereF 1 andF 2 denote the estimated values of F 1 and F 2 , respectively, and E s denotes the sample mean measured over 1000 independent realizations of noise. In order to obtain meaningful BIAS, VAR, and RMSE when the frequencies are unresolved, that is when both F 1 and F 2 have only one estimate that lies within the interval [F 1 − 50, F 2 + 50] Hz, we assign the same value for bothF 1 andF 2 and compute these metrics as shown in Figure S4. Any conclusion on the results in this analysis should be made by simultaneously analyzing Figures S3 and S4.
To illustrate further, consider the case where S/N = 10. By applying apodization and by using dissipative spectra (real part of the Fourier spectra), in Figure S3, we observe that FRI-NMR is able to resolve frequencies at 5 Hz, whereas the Fourier method is able to distinguish frequencies separated by 17.5 Hz. Comparing LPSVD and FRI-NMR, we find that both techniques are able to estimate frequencies separated by 5 Hz but LPSVD is able to resolve frequencies only 10% of the time compared with FRI-NMR, which is able to resolve frequencies about 70% of the time.
Further, comparing the plots in Figure S4, the flat region in the bias and MSE curves of the Fourier plots signify that the Fourier method is unable to resolve the frequencies and gives only a single peak about the average (F 1 +F 2 )/2. Hence, the error is constant. As the frequency separation increases beyond 17.5 Hz, the Fourier method is able to resolve the frequencies, which is also shown in Figure   The dissociation constant for the Ubiquitin Gold-nanorod (AuNR) interaction was determined assuming a fast exchange. First, the number of Ubiquitin molecules binding one nanorod was estimated based on the surface area of the nanorods, which was calculated based on its dimensions observed in the transmission electron microscopy (TEM) (Fig. 5b of the main document) and depicted below.
Assuming the Ubiquitin molecules are uniformly and tightly packed on the surface of the nanorod, the total surface area was divided by the area occupied by one Ubiquitin molecule. The area occupied by one Ubiquitin molecule on the surface was conservatively estimated as a circle having a radius corresponding to the hydrodynamic radius of Ubiquitin (1.5 nm). Using this approach, the total number of Ubiquitin molecules adsorbed on the surface of the gold nanorod was approximately estimated as 600.
Next, we assumed that all the binding sites on the nanorods are independent of each other and the binding of the peptide to one site does not affect its binding to another site on the AuNR surface.
The dissociation constant (K D ) was then obtained by fitting Equation (S.22) to the chemical shift data obtained from the titration experiment, which is valid for fast exchange and for a system with multiple binding sites for the ligand 43 . Here, the protein represents the ligand and the nanorod represents the macromolecule on which binding takes place: where δ obs and δ free denote the observed chemical shifts of the bound peptide and free peptide, respectively.
The symbols [AuNR] total and [Ubq] total represent the total concentration of AuNR and Ubiquitin at a given point during the titration, respectively. In (S.22), we have δ max = δ bound −δ free , where δ bound is the chemical shift in the completely bound form. The δ max is obtained as part of the fitting procedure 43 . Using this equation, the K D was estimated as 27 ± 3 µM. The number of ubiquitin molecules adsorbed on the gold-nanorod at any given time was taken as 600 (see Section S5) and the dissociation constant (KD) indicated for each residue was obtained by fitting the shifts to the equation shown in Section S5 corresponding to fast exchange. The error bars correspond to the estimated uncertainty in the measurement of 1 H shifts (0.0025 ppm) as described in the main text.