Nuclear Singlet Relaxation by Chemical Exchange
The population imbalance between nuclear singlet states and triplet states of strongly coupled spin-1/2 pairs, also known as nuclear singlet order, is well protected against several common relaxation mechanisms. We study the nuclear singlet relaxation of 13 C pairs in aqueous solutions of 1,2-13 C 2 squarate, over a range of pH values. The 13 C singlet order is accessed by introducing 18 O nuclei in order to break the chemical equivalence. The squarate dianion is in chemical equilibrium with hydrogen-squarate (SqH − ) and squaric acid (SqH 2 ) characterised by the dissociation constants pKa 1 = 1.5 and pKa 2 = 3.4. Surprisingly, we observe a striking increase in the singlet decay time constants T S when the pH of the solution exceeds ∼ 10, which is far above the acid-base equilibrium points. We derive general rate expressions for chemical-exchange-induced nuclear singlet relaxation and provide a qualitative explanation of the T S behaviour of the squarate dianion. We identify a kinetic contribution to the singlet relaxation rate constant which depends explicitly on kinetic rate constants. Qualitative agreement is achieved between the theory and the experimental data. This study shows that infrequent chemical events may have a strong effect on the relaxation of nuclear singlet order.
I. INTRODUCTION<!>A. Sample preparation<!>B. Instrumental details<!>T S measurements<!>tion.<!>Isotope shifts<!>B. Spin-lattice relaxation<!>Singlet decay time constants<!>B. Site-specific Relaxation Superoperator<!>C. Analytical relaxation rate constants<!>D. Application to 1,2-13 C-squarate<!>E. Physical Interpretation<!>V. DISCUSSION<!>CONFLICTS OF INTEREST<!>Combined matrix representations<!>(B2)<!>Perturbative treatment
<p>Long-lived nuclear spin states (LLS) are configurations of nuclear spins which are protected against common relaxation or dissipation mechanisms and which display extended decay time constants [1][2][3][4][5][6][7][8][9][10][11][12][13][14][15][16][17][18][19] . A seminal example is the singlet order of a spin-1/2 pair ensemble, corresponding to a population imbalance between the singlet and triplet states of the spin-1/2 pairs 10,16 . The decay time constants of such states are often many times the spinlattice relaxation time constant T 1 , with lifetimes exceeding 1 hour being observed in favourable cases 12,19 . Longlived states have been applied to the study of slow chemical and transport processes [20][21][22][23][24][25][26][27][28][29] , to the characterisation of biomolecular ligand binding [30][31][32][33][34] , and for the transport and storage of nuclear hyperpolarization 5,7,[35][36][37][38][39][40][41][42][43] . Longlived nuclear singlet order plays a central role in the generation of nuclear hyperpolarization from hydrogen gas enriched in the para spin isomer [44][45][46][47][48][49][50][51][52][53][54][55][56][57][58][59] .</p><p>The relaxation mechanisms of LLS have been extensively studied 3,6,15,[60][61][62][63][64][65][66][67][68][69][70][71] .</p><p>Most attention has been given to the dipole-dipole (DD), chemical shift anisotropy (CSA) and random field relaxation mechanisms 3,6,60,61,64,68 . More unusual relaxation mechanisms such as nuclear spin-rotation and scalar relaxation of the second kind have also been identified 4,15,[65][66][67]69 . Generally speaking, the efficiency of a relaxation mechanism is roughly proportional to the rotational correlation time τ c for the random tumbling of the molecules in solution. This suggests that small molecules, which tumble more rapidly in solution and have small values of τ c , tend to have the longest LLS lifetimes. In most cases, that is also the observed behaviour.</p><p>a) mhl@soton.ac.uk Figure 1. The 1,2-13 C2-squarate dianion. In the absence of 13 C labelling, the ion belongs to the point group D 4h and has 4-fold rotational symmetry. In the current work, random 18 O enrichment is used to slightly break the chemical equivalence of the 13 C sites by means of secondary isotope shifts.</p><p>When labelled with two 13 C nuclei, the small, highly symmetrical squarate molecule is an attractive candidate for exhibiting long-lived 13 C 2 singlet order in aqueous solution (see figure 1). The squarate dianion is stabilised by electron delocalisation and exhibits four-fold rotational symmetry [72][73][74][75] . Since the dominant oxygen isotope 16 O has no nuclear spin, the 13 C pair of 1,2-13 C 2 -squarate is magnetically isolated, eliminating all intramolecular dipole-dipole contributions to singlet relaxation. Furthermore, it is technically feasible to introduce non-magnetic 18 O nuclei into the squarate dianion, as has been achieved previously for oxalate 76 . The increased mass of the 18 O nuclei modifies the vibrational wavefunctions, inducing small secondary isotope shifts of the 13 C resonances 77 . These shifts may break the exact magnetic equivalence of the 13 C nuclei, providing access to long-lived 13 C 2 singlet order. Large values of the singlet relaxation time constant T S are therefore anticipated for the 13 C pair in 18 O-enriched 1,2-13 C 2 -squarate, especially at low magnetic field, when the chemical shift anisotropy mechanism of singlet relaxation is removed.</p><p>As described below, our initial measurements of the singlet decay time constant T S for the 13 C 2 pairs in 18 O-enriched 1,2-13 C 2 squarate were disappointing. In aqueous solution and neutral pH we observed a singlet decay time constant of around T S ∼ 30 s in a field of 9.4 T, which is slightly shorter than the value of T 1 under the same conditions. Conventional mechanisms such as CSA-induced relaxation are insufficient to explain this low value of T S . Furthermore, the observed value of T S for 18 O-enriched 1,2-13 C 2 squarate was found to be strongly dependent on pH, with a rather sharp transition from short T S values at low pH to long T S values at high pH, approaching 300 s in low magnetic field (see figure 2).</p><p>The strong pH-dependence of T S suggests that the reversible protonation of squarate in aqueous solution plays a role in the singlet relaxation. In aqueous solution, the symmetrical squarate dianion is in equilibrium with hydrogen squarate and squaric acid (see figure 3). The possibility arises that this dynamic chemical process could act as a relaxation mechanism for nuclear singlet order. However, at first sight, this appears unlikely. The pKa values for the two deprotonation stages of squaric acid are 1.5 and 3.4 respectively 78 . The equilibrium fractions of the monoprotonated and diprotonated species at a pH of 7.0 are only ∼250 ppm and ∼0.8 ppb respectively. Such minority species could only have a significant effect if a mechanism exists that leads to an extreme sensitivity of nuclear singlet order to rare protonation events.</p><p>A chemical exchange mechanism for singlet relaxation has been proposed before, in the context of ligand binding [30][31][32][33][34] . Those reports assume that the singlet relaxation rate constant T −1 S for an exchanging system is given by a population-weighted average of the relaxation rate constants for the participating chemical species 30 . This theory is used to show that singlet relaxation may be used to assess the binding of ligands to proteins, with significant biochemical applications [30][31][32][33][34] .</p><p>However, the theoretical framework used for ligand binding studies [30][31][32][33][34] cannot explain the data in figure 2. If T −1 S were given by a population-weighted average of the T −1 S values for the participating species, the strongest pHdependence would be expected in the region of the pKa values, i.e. below pH ∼ 4, not around pH ∼ 11, as is observed.</p><p>The data in figure 2 therefore indicate a strong chemical-exchange mechanism for singlet relaxation mechanism that has a much stronger effect than the mere averaging of the rate constant T −1 S over the participating species. We propose a mechanism of this kind below. The protonation of a squarate oxygen atom causes a substantial change in the electronic structure which leads in turn to a large isotropic chemical shift difference between the two 13 C sites. In high magnetic field the isotropic shift difference of the protonated species, measured in frequency units, greatly exceeds the scalar spin-spin coupling between the two 13 C nuclei. Hence, in the protonated form, the singlet state of the 13 C 2 pair is not even an approximate eigenstate of the coherent spin Hamiltonian. This leads to a rapid loss of singlet order upon each protonation event. As discussed below, within a certain kinetic regime, the singlet relaxation rate constant T −1 S becomes directly proportional to the rate of protonation of the squarate dianion. The protonation rate is proportional to hydronium concentration and is therefore strongly pH-dependent.</p><p>We provide approximate theoretical expressions for the singlet relaxation rate constant of a spin-1/2 pair ensemble, using a simplified model of 2-site exchange. Three terms are shown to contribute to the singlet relaxation rate constant T −1 S in the presence of chemical exchange. The first two terms are the population-weighted averages of the T −1 S values for the participating species. These terms correspond closely to the theory developed for ligand binding applications of singlet NMR [30][31][32][33][34] . Since these terms depend on equilibrium populations, they do not depend on the individual values of the kinetic rate constants, but only on their ratio. The third term is the most important one, in the current context. It depends on the individual kinetic rate constants of the chemical transformations and cannot be described as a populationweighted average of some quantity. We therefore call the third term the kinetic contribution to the singlet relaxation rate constant. Since proton exchange is usually a very fast process, the third term (kinetic term) can easily dominate, even when some of the participating species have negligible equilibrium concentrations. That is the key to understanding the remarkable pH-dependent T S data in figure 2.</p><p>Apart from leading to an understanding of singlet relaxation in systems undergoing chemical exchange, the results described here indicate that nuclear singlet relaxation may be an exquisitely sensitive reporter of rare chemical events. For example, an individual molecule of 18 O-enriched 1,2-13 C 2 squarate, in a solution with pH ∼ 11, only picks up a proton about once a minute on average, and that proton is resident for only a fraction of a millisecond, on average, before returning to the solvent. The 13 C 2 singlet relaxation measurement detects the rare protonation events by means of a significant fall in the T S value. It is hard to find other methods for the reliable characterisation of such infrequent chemical events under equilibrium conditions.</p><!><p>The synthesis of 1,2-13 C 2 squaric acid (3,4dihydroxycyclobut-3-ene-1,2-dione-1,2-13 C 2 ) will be described in a following publication. For all experiments 5.8 mg of 1,2-13 C 2 squaric acid[ 16 O 4 ] were dissolved in 0.5 mL of a 1:1 mixture of H 2 18 O:D 2 16 O leading to a concentration of 0.1 M. The isotopic purity of the H 2 18 O was 97.1%. Solutions were incubated at 80 ○ C for two hours allowing the 16 O/ 18 O exchange process to reach its dynamical equilibrium. The pH value of each sample was adjusted by step-wise addition of aqueous solutions of 0.2 M NaOH or 0.2 M HCl. In each instance the total volume added did not exceed 50 µL leading to a negligible effect on the sample concentration. The pH values were estimated by using a pH meter equipped with a Hamilton SpinTrode probe (www.hamiltoncompany.com), suitable for pH measurements inside NMR tubes. For the sake of simplicity, we ignore the potential complications associated with solvent deuteration, such as: (i) the influence of solvent deuteration on the acid-base equilibrium; (ii) the distinction between the concentration of the hydronium ion, indicated by pH, and that of its deuterated analogues; (iii) the modified response of the pH-meter, induced by solvent deuteration.</p><!><p>Experiments were performed on 400 MHz and 700 MHz Bruker Avance Neo systems. For all samples the recycling delays between transients were set to 60 s corresponding to at least three T 1 periods. 13 C pulse nutation frequencies were given by ∼13.9 kHz at 9.4 T, and ∼10.7 kHz at 16.5 T. T 1 measurements used a single transient for each data point. T S measurements (see below) used averaging over 1 to 8 transients in order to achieve a signal-to-noise of ∼16 for the first data point. The NMR signals were sampled with 131k data points and a spectral width of 146 ppm.</p><p>Apart from experiments at fixed magnetic fields strengths of 16.5 T (700 MHz) and 9.4 T (400 MHz), field cycling experiments were performed inside the fringe field of the 9.4 T magnet. An overview of the fringe field as a function of the detection coil center to center of the sample tube is given in figure 4. Samples were prepolarised at 9.4 T before being shuttled to a fringe field region of 1.9 T as indicated in figure 4. The sample transport was performed using a stepper motor driven sample carriage operated at a fixed velocity of 1 m s −1 , as described in reference 79. Typical shuttle times were on the order of ∼ 300 ms. Signal detection was performed by shuttling the sample back into the detection region of the high field NMR magnet and applying a π 2 pulse. Longitudinal relaxation time constants at 16.5, 9.4 and 1.9 T were determined with a standard inversion recovery experiment. For T 1 measurements at 1.9 T the sample was first prepolarised at 9.4 T, the magnetisation then inverted and shuttled to 1.9 T as described above. After a variable recovery period the sample was shuttled back to the detection region of the high field magnet to indirectly observe its recovery to thermal equilibrium at 1.9 T.</p><!><p>Figure 5. Pulse sequence protocol for T S measurements at variable magnetic strength. Singlet order is excited via an M2S block on the 13 C channel (a). After singlet order excitation the sample is shuttled (optional) to a desired magnetic field strength (b) and allowed to evolve for a time period τ evol (c). In case of a field cycling experiment the sample is shuttled back into the sweet spot of the detection magnet (d). After the free evolution period unwanted spin operator components are filtered out using a gradient based singlet filtration step (e). Signal detection is initialised by reconverting unobservable singlet order into observable magnetisation via the S2M sequence (f).</p><p>An overview of the experimental procedure for T S estimation is given in figure 5. (a) The initial M2S block (see figure 6a) converts thermally polarised 13 Cmagnetisation into homonuclear 13 C-singlet order 7,8 . Experimentally optimised M2S parameters for 18 O-enriched 1,2-13 C 2 squaric acid are given as a function of pH in table I. (b) After singlet order excitation an optional field cycling step is introduced to transport the sample to a desirable magnetic field strength. (c) Singlet order is allowed to freely evolve for a variable duration τ evol . (d) For the case of field cycling experiments the sample is shuttled back into the active region of the magnet. (e) Application of a gradient-based singlet order filter (see figure 6b) isolates NMR signals passing through singlet order 80 . The field gradient parameters for the singlet filter are summarised in table II. (f) The S2M sequence (the chronological reverse of the M2S sequence) reconverts the singlet order back into observable magnetisa-</p><!><p>In order to compensate for field homogeneity effects, the individual pulses for the M2S and S2M sequence were replaced by composite pulse counterparts. In detail, 180 φ pulses were replaced by 90 90+φ 180 φ 90 90+φ sequences; 90 φ pulses were replaced by 180 97.2+φ 360 291.5+φ 180 97.2+φ 90 φ sequences 81,82 . The free evolution delays were adjusted to account for the duration of the elongated pulse elements. For example, echo delays were shortened by half the duration of the refocusing element. The 13 C shifts may be interpreted qualitatively by assuming rapid exchange between the protonated and deprotonated forms in figure 3. If the exchange is assumed to be rapid relative to the 13 C chemical shift changes, the observed 13 C chemical shift corresponds to the mean value for the exchanging species, weighted by their populations. The pH-dependence is weak above pH ∼ 7 since the deprotonated Sq 2− species dominates in this regime.</p><p>The observed broadening of the 13 C spectrum around pH ∼ 4 indicates that in this pH regime, the exchange of the protonated forms occurs on a timescale comparable to the changes in 13 C chemical shift frequencies. A quantitative analysis of the pH-dependent 13 C chemical shifts and linewidths would require consideration of a complex exchange network with numerous exchanging species and has not been attempted.</p><!><p>At most pH values, the 13 C spectra display a triplet fine structure due to secondary isotope shifts induced by the randomly distributed 18 O nuclei (see figure 7b). A representative 13 C equilibrium spectrum of 18 O-enriched 1,2-13 C 2 squarate at 16.5 T and pH = 14.0 is shown in figure 8. The spectrum displays three sharp peaks. The two outer peaks are separated by ≃ −15.0 ppb and ≃ +15.0 ppb from the central peak. The integrated intensities of the three peaks are in the approximate ratios of 1.2 : 2.0 : 1.0.</p><p>In this work we follow the convention of an isotope shift being defined as the change in the chemical shift upon substitution of the lighter nucleus (l) by the heavier nucleus (h)</p><p>where X represents the observed nucleus and n the number of chemical bonds between X and h. This convention was used in a study of 18 O-enriched oxalate 76 , but differs in sign from the most common convention for secondary isotope shifts (see, for example, Hansen 77 ). At high pH, the Sq 2− species dominates. Under these conditions the spectrum may be interpreted in similar fashion to 18 O-enriched oxalate 76 . The following assumptions are made:</p><p>1. Two-bond and three-bond 18 O-induced isotope shifts are ignored.</p><p>2. The magnitude of one-bond 18 O-induced isotope shifts is small compared to the J-coupling of the 13 C pair.</p><p>These assumptions allow the 13 C spin system of each 18 O-isotopologues to be classified as either a magnetically equivalent A 2 spin system, or a near-equivalent AB spin system. Each A 2 spin systems gives rise to a single 13 C peak, while in general, each near-equivalent AB system gives rise to a four-peak subspectrum with two strong central peaks and two weak outer peaks. However, under the conditions that apply here, the outer peaks have small intensity, and the two inner peaks are unresolved and coalesce to a single peak at the mean of the two 13 C chemical shifts. The experimentally observed spectral pattern may therefore be understood as a superposition of these single 13 C resonances (see figure 9). The observed mean chemical shifts of the different peaks are thus given by</p><p>where 1 ∆ 13 C( 18 O) represents the 18 O-induced singlebond isotope shift. The spectral separation between the high-δ and low-δ peak is therefore equal to the one-bond 18 O-induced isotope shift experienced by a 13 C nucleus. At high pH values the experimental data (figure 8) give the following estimate of the one-bond 18 O-induced isotope shift:</p><p>The intensities of the different peaks are proportional to the abundances of the different isotopologues. From simple statistical arguments, assuming a random distribution of 18 O nuclei amongst the oxygen atoms, the abundances of the species giving rise to the three peaks are proportional to {(1 − x) 2 , 2x(1 − x), x 2 }, where x is the abundance of 18 O nuclei, and 1 − x is the abundance of 16 O nuclei. The observed spectral amplitudes are consistent with an 18 O enrichment level of x ≃ 46%. This corresponds well to the known composition of the solvent (see section II A).</p><p>The 18 O isotope shifts become significantly larger at low pH, as shown in figure 7(b). This effect has not been analyzed in detail, but is provisionally attributed to small changes of the equilibrium fractions of the different protonated species upon 18 O substitution. 18 O substitution changes the zero-point vibrational energy of -OH bonds and therefore influences the equilibrium fractions of the various protonated species, which leads in turn to a perturbation of the 13 C chemical shifts when averaged over all exchanging species. This mechanism has been identified as the source of relatively large and remote 18 O-induced isotope shifts in the 13 C spectra of some carboxylic acids 83,84 .  In contrast to the observed T S values, only a weak pHdependence of T 1 is observed. For example, at 9.4 T we observe T 1 values of approximately 30 s for all pH values apart from pH = 4.1, where T 1 decreases to ≈ 15 s. In particular, T 1 does not display any notable features either at the pKa values of squaric acid, nor at the inflection point pH ≃ 10 at which the observed T S values start to drastically increase.</p><!><p>Similar observations are made for the pH-dependent T 1 values at 16.5 T and 2 T with an overall decrease in T 1 as the magnetic field strength increases. This is consistent with a dominant CSA relaxation mechanism, proportional to the square of the magnetic field strength.  O:D2 16 O squarate at pH = 14.0 (black line). The spectrum was acquired at 16.5 T using the pulse sequence protocol shown in figure 5 with no evolution delay. The blue line shows the spectrum after application of a single π 2 pulse. Both spectra are referenced to a carrier frequency of 203.05 ppm.</p><p>As shown in figure 9, the central peak of the 13 C triplet is associated with 1,2-13 C 2 -squarate species with an 18 O atom attached to one 13 C atom and a 16 O atom attached to the other. These species contain strongly coupled 13 C pairs with a small chemical shift difference between the two 13 C nuclei, and are suitable candidates for singlet order excitation techniques. Figure 11 shows a singlet-filtered 18 O-enriched 1,2-13 C 2 squarate spectrum at pH = 14.0. The spectrum was acquired at 16.5 T using the pulse sequence protocol in figure 5 with the M2S parameters given in table I. The signals from the outer peaks of the triplet are suppressed in the singlet-filtered 13 C spectrum since these signals are generated by isotopologues with no 18 O atoms directly bonded to both 13 C atoms, or with no 16 O atoms directly bonded to both 13 C atoms. Although these species include isotopologues with asymmetric distributions of 18 O atoms on the other side of the molecule to the 13 C 2 pair, these remote heavy atom substitutions are too far away to allow appreciable excitation of 13 C singlet order, for the chosen pulse sequence parameters.</p><p>Depending on the pH of the solution the intensity of the singlet filtered peak displays approximately 10% to 30% of the intensity of the unfiltered spectrum. The small frequency shift of the singlet-filtered spectrum with respect to the pulse-acquire spectrum is attributed to a slight change in sample temperature between experiments.</p><p>The experimentally optimised M2S pulse sequence parameters (see table I) serve as an independent estimator for the 1 ∆ 13 C( 18 O) isotope shift and the scalar coupling constant J 12 . In the absence of relaxation losses the relationship between the optimal M2S parameters and the spin interaction parameters are given by 8,85 :</p><p>where ∆Ω is the resonance frequency difference between the two spins in rad s −1 . Inversion of these relations leads to an estimate of J 12 ≃ 54.86 ±0.5 Hz for the scalar coupling constant. The resonance frequency difference ∆Ω is consistent with the 18 O-induced isotope shift estimated from the spectral splittings (see section III A).</p><p>The experimentally optimised M2S parameters summarised in table I indicate that similar estimates are obtained for the majority of the other pH values. The only exception is the pH region 2 < pH < 7, where we were unable to observe any singlet-filtered signal components, presumably due to short T 2 values associated with the exchange broadening at intermediate pH which is readily visible in figure 7. The short T 2 values cause strong signal losses during the M2S/S2M pulse schemes.</p><!><p>The pulse sequence protocol shown in figure 5 was used to study nuclear singlet decay time constants of 18 Oenriched 1,2-13 C 2 squaric acid at a set of pH values and magnetic fields. Figure 2 summarises the nuclear singlet decay time constants as a function of pH at 16.5, 9.4 and 2 T.</p><p>At all magnetic fields, the T S values of 18 O-enriched 1,2-13 C 2 squaric acid vary strongly with the pH value of the solution. The overall pH dependence resembles a sigmoidal curve. For example, at 9.4 T the nuclear singlet decay time constant of 18 O-enriched 1,2-13 C 2 squaric acid at pH = 1.3 is given by T S = 18.3 s. Initially, increasing the pH does not lead to any significant changes in T S . However at a pH value of ∼ 11 a sharp increase in T S is observed reaching values up to T S = 53 s. The maximum T S value of 275 s is reached for pH = 14.0. Similar behaviour is observed at 16.5 T and 2 T with an overall increase in relaxation times as the magnetic field strength decreases. For example, at 16.5 T, T S reaches a maximum of T S = 121 s, whereas at 2 T the observer singlet decay time constant reaches a maximum of T S = 380 s. The increase in T S with a decrease in the magnetic field strength is consistent with a dominant chemical shift anisotropy relaxation mechanism in this pH range.</p><p>As pointed out above, one of the most surprising features of the pH-dependence of T S is that the transition from fast to slow relaxation occurs around the pH-value of 11, which is far above the pKa values of the acid-base equilibria for this system. The two spins supporting the nuclear singlet state may find themselves in two different molecular environments A (red) and B (blue). This change may be induced upon binding of an otherwise NMR inactive partner (grey ball). The exchange dynamics are described by an exchange matrix</p><p>The spin Hamiltonians are in general different for configuration A and B. For simplicity we assume that the scalar coupling constant J 12 = ω J (2π) is identical for both configurations.</p><p>Note that the subscripts in symbols such as k BA should be read from right to left, i.e. k BA denotes the rate constant for the chemical transformation of species A into B, and vice versa for k AB .</p><p>In the following discussion, we use the symbol X to denote a general chemical species, where X ∈ {A, B} for 2-site exchange.</p><p>The coherent dynamics of a spin system in chemical species X are governed by the following Hamiltonian:</p><p>(</p><p>The mean (rotating frame) resonance frequency of the spin pair is denoted Ω</p><p>2 ). The resonance frequency difference between the two spins in configuration X is given by ∆Ω X = Ω X 1 − Ω X 2 . At a magnetic field B, this is given by ∆Ω X = −γB∆δ X , (6)   where ∆δ X is the difference in chemical shifts of the two 13 C sites in configuration X.</p><p>The Hamiltonians H A coh and H B coh may be used to define the site-specific parts of the coherent Liouvillian</p><p>where ĤX is the commutation superoperator for Hamiltonian H X . In general, between exchange events the subensembles A and B evolve under a site-specific Liouvillian LX including coherent and incoherent parts</p><p>where ΓX represents a site-specific relaxation superoperator, discussed below. The evolution of the sub-ensembles A and B including exchange events may be described by coupling their site-specific evolution via the forward and backward rates k BA and k AB [86][87][88] . The resulting equations of motion are given by</p><p>The evolution of ρ A and ρ B are conveniently combined by introducing a joint density operator</p><p>and a joint generator of motion Ĝ</p><p>Equation 9 then reduces to a single system</p><p>and the evolution of the combined ensemble is completely characterised by Ĝ.</p><!><p>The site-specific relaxation superoperator ΓX is typically derived from a microscopic model and captures details about the fluctuating part of the spin Hamiltonian 64 . For isolated spin-1/2 pairs stochastic modulations of the dipolar and chemical shift interactions are the dominant contributions to the fluctuating Hamiltonian. Any additional relaxation features are approximated by coupling the nuclear spins to an external random field 89-91</p><p>Within the fast motion and isotropic diffusion limit the site-specific relaxation superoperator for dipolar driven spin relaxation takes the form</p><p>The rotational correlation time and dipole-dipole spinspin coupling constant (in rad s −1 ) for a molecule in configuration X are denoted by τ X c and b X</p><p>12 respectively. The operators T 12 2m (X) are site-specific second-rank spherical tensor operators resulting from the coupling of the two spin angular momentum operators 91 .</p><p>In general, the equilibration of the spin system with a finite-temperature environment should be taken into account by replacing the double commutators in equation 14, and similar equations, by thermally weighted Lindbladian dissipators 92 . This complication is not relevant to the current discussion and has been avoided.</p><p>The chemical shift anisotropy tensors are in general different for spin systems in configurations A and B since the electronic environment of the spins is altered upon a binding event. To simplify the discussion we assume axially symmetric chemical shift tensors for both sites, η i (X) = 0. X ∈ {A, B}. The corresponding relaxation superoperator may be expressed as follows</p><p>The frequency ω i csa (X) = −γB 0 δ i csa (X) represents the chemical shift anisotropy parameter in rad s −1 . The cross terms are weighted by the second rank reduced Wigner matrix element</p><p>which in turn depends on the angle between the unique principal axes of the two CSA tensors in configuration X, denoted β X P P , and the geometric mean of the anisotropy parameters ξ i csa (X). The operators X i 2m (X) are defined by</p><p>These operators are components of the irreducible spherical tensor resulting from coupling the spherical tensor spin operators with the spherical components of the static magnetic field 91 . The symbols C JM jmj ′ m ′ denote Clebsch-Gordan coefficients.</p><p>An external random-field contribution to the relaxation superoperator is often included to take into account unspecified mechanisms, such as interactions with paramagnetic species in solution, or internuclear dipole-dipole couplings with other magnetic nuclei within the same molecule or in different molecules. The corresponding site-specific relaxation superoperator is given by [89][90][91]</p><p>where τ X erf describes the correlation time of the random field fluctuations, ω X erf denotes the root-mean-square amplitude in rad s −1 , and −1 ≤ κ X ≤ 1 the correlation coefficient of the fluctuations experienced by spins 1 and 2, in configuration X.</p><!><p>An analytic approximation to the singlet order decay rate constant under chemical exchange may be constructed by considering the eigenvalues of the join generator of motion Ĝ. For a two-site exchange model the matrix representation of Ĝ is of dimension dim( Ĝ) = 2 × N 2 H , where N H is the dimension of the underlying Hilbert space. But as discussed in appendix B the singlet order dynamics for the combined ensemble (A+B) are well approximated by the following 4 × 4 matrix sub-block</p><p>The sub-block g is spanned by the operators</p><p>The operators SO X represent singlet order for spin systems in configuration A or B, whereas the operators O X ∆ represent zero-quantum operators that are anti-symmetric under spin permutation.</p><p>The terms R 0 S (X) and R 0 ∆ (X) in equation 19 are given by: R 0</p><p>We introduce the following symbol for a weighted average over the equilibrium populations of the exchange process:</p><p>Appendix B gives a perturbative treatment of two-site exchange under the assumptions k BA ≫ ω J and pH > pKa, and shows that the decay constant R S for singlet order may be approximated by</p><p>The first term is given by the population-weighted average of the R 0 S terms for the individual species, as given by equation 21:</p><p>The second term is given by</p><p>This term is due to the "leakage" of the singlet spin order operators into other terms, in particular O ∆ , through the influence of the chemical shift frequency differences ∆Ω X , partially suppressed by ω J 3,6 . This behaviour is captured by the inclusion of the term R leak S .</p><p>The third term is the most important one in the current context. It is given by</p><p>.</p><p>(26) This is called here the kinetic term since it depends on the kinetic rate constants k AB and k BA , and cannot be expressed as a population-weighted average. Unlike the R 0 S and R leak S terms, the kinetic term R kin S has no counterpart in the absence of chemical exchange. The term R kin S expresses the direct influence of the chemical exchange process on the singlet relaxation, through the stochastic modulation of the chemical shift differences by the chemical process.</p><p>The treatment of singlet relaxation under ligand binding assumed that the singlet relaxation rate constant may be treated as a population-weighted average of the singlet relaxation rate constants of the individual species [30][31][32][33][34] , omitting the kinetic term. Although this assumption might have some merit under the conditions of ligand binding, it is certainly not true in general.</p><!><p>The proton exchange of squarate in solution is a complex problem since three chemical species, Sq 2− , SqH − and SqH 2 are involved (see figure 3). Furthermore, each of the protonated 1,2-13 C 2 species displays several isotopomers, depending on the position of the 13 C 2 pair with respect to the protonation sites. For the sake of simplicity, we consider only the high-pH regime, ignoring the neutral SqH 2 species, and neglecting the distinction between the multiple 1,2-13 C 2 isotopomers of SqH − . Protonation (or deuteration) of squarate is taken into account solely through its influence on the 13 C chemical shifts. We ignore any direct effects of the proton or deuteron nuclear spins, through heteronuclear dipoledipole couplings or scalar couplings, and also ignore the effect of protonation on other interactions such as the 13 C-13 C J-coupling. We also ignore the effect of 18 O or 2 H on the chemical equilibria. Although this model is highly simplified, we show below that it reproduces the main features of the observed data. We consider the following chemical reaction</p><p>characterised by the protonation second-order rate constant k p (in units of M −1 s −1 ), and the pseudo-first-order rate constant k d for the deprotonation (in units of s −1 ).</p><p>The symbol M = mol dm −3 indicates molarity. In line with common practice, the solvent concentration is assumed to be effectively constant and factored into the deprotonation rate constant k d . The kinetic rate constant k AB is related to the protonation rate constant k p as follows</p><p>The concentration of the hydronium ions depends on the pH value of the solution through</p><p>The exchange rate constant k AB may therefore be expressed as follows:</p><p>The exchange rate constant k BA is simply equal to the deprotonation rate constant</p><p>At the equivalence point pH = pKa the equilibrium concentrations of Sq 2− and SqH − are equal. We thus have</p><p>The required exchange rate constants are then given by</p><p>Hence the kinetic rate constants k AB and k BA are solely characterised by the protonation rate constant k p , the pK a value of the acid-base equilibrium, and the sample pH. Since there is no available independent data on the protonation rate constant k p , we treat this quantity as an adjustable parameter. Equations 23 to 26 may be combined with equation 33 to provide an estimate for the decay rate constant R S of 13 C 2 singlet order in 18 O-enriched 1,2-13 C 2 squarate, in the high-pH regime. The singlet relaxation time constant T S is given by the inverse of the rate constant, T S = R −1 S . Evaluation of T S requires estimates of the isotropic shift differences ∆Ω A and ∆Ω B , the 13 C chemical shift anisotropies and dipole-dipole coupling parameters, and the correlation time τ c of the fluctuating interactions. The parameters used for the numerical modelling are given in table III. The pKa value was chosen to equal the known pKa 2 = 3.4 of squaric acid. The chemical shift differences and the CSA parameters are based on the computational chemistry estimates given in appendix A. The dipole-dipole coupling constant is based on the C-C bond length of ∼ 150 pm for the squarate dianion 93 . The rotational correlation time τ C = 100 ps agrees reasonably well with molecular dynamics simulations of the squarate dianion in aqueous solutions 93 . The minor influence of intermolecular interactions, including interactions with dissolved oxygen, is taken into account by a random-field contribution quantified by the parameter ξ erf (1 − κ), recognizing that it is difficult to disentangle the correlation coefficient κ from the mean-square fluctuation amplitude ξ erf . The parameter ξ erf (1 − κ) was chosen to match the experimentally observed T S values at pH = 14.</p><p>The solid curves in figure 13 shows analytical curves for T S = R −1 S as a function of pH and magnetic field, using equation 23 for the rate constant, with parameters given in table III. The protonation rate constant k p = 10 7 M −1 s −1 was chosen to obtain a rough qualitative match with experiment. This value is reasonably consistent with known protonation rate constants for carboxylic acids in solution, determined by NMR spectroscopy 94,95 and fluorescence spectroscopy 96,97 .</p><p>Figure 13 also shows experimental values for T S at the magnetic field of 9.4 T. Further values of T S at different magnetic fields are shown in figure 2. Although the correspondence between theory and experiment is qualitative at best, the salient features of the experimental data are reproduced. In particular the two-site model does reproduce the rapid transition from short T S to long T S at a pH value which greatly exceeds the pK a value of the acidbased equilibrium. The analytical curves also predict the decrease in the high-pH limit of T S with increasing magnetic field, which is due to the dominant chemical shift anisotropy relaxation in this regime.</p><p>The match between experiment and theory cannot be significantly improved within the framework of the twosite exchange model. In particular, adjustment of the protonation rate constant k p within reasonable limits does not lead to a significant improvement in the quality of the match. We attribute the lack of a quantitative correspondence between experiment and theory to the obvious shortcomings of the highly simplified two-site exchange model.</p><p>A quantitative description of the 13 C 2 singlet relaxation in 18 O-enriched 1,2-13 C 2 squarate would require a much more sophisticated kinetic model, taking into account the multiple SqH − isotopomers. Such a study would be a considerable undertaking and we have not attempted it.</p><!><p>Consider an individual Sq 2− molecule in a high-pH solution. Protonation events are rare. For example, at pH = 10, the protonation rate constant is estimated through equation 28 to be k AB ≃ 1 × 10 −3 s −1 . Hence, at pH=10, a given Sq 2− molecule has to wait for about 1000 s before being protonated, on average.</p><p>From the best-fit protonation rate constant k p = 10 7 M −1 s −1 and the pKa value of pKa 2 = 3.4, the deprotonation rate constant of SqH − is estimated through equation 32 to be given by k d ≈ 4.0 × 10 3 s −1 . This implies that the mean lifetime of a SqH − molecule, before losing its proton, is</p><p>at pH = 10. Hence, once a Sq 2− anion is protonated, the species SqH − persists in solution for around 250 µs on average, before the proton is abstracted by the solvent again.</p><p>The lifetime of SqH − molecules in solution is statistically distributed. The probability δp d for an individual SqH − molecule to lose its proton during the small interval between time points t = τ and t = τ + δτ is given by 98</p><p>Hence the probability of a SqH − molecule keeping its proton during the same small time interval is given by</p><p>where k d is the deprotonation rate constant. A finite interval τ may be subdivided into n = τ δτ small time intervals δτ . Assuming statistical independence of events in successive time intervals (Markov process) 99,100 , the probability p H (τ ) for a molecule to remain protonated during the entire interval τ may be calculated as the product of the individual probabilities. In the limit of large n the probability p H (τ ) becomes</p><p>It follows that the probability of a molecule remaining protonated for the finite interval τ , and then losing its proton during the following small time interval [τ, τ +δτ ], is given by</p><p>Consider now a 18 O-enriched 1,2-13 C 2 squarate molecule in the singlet state:</p><p>Since the isotropic shift difference between the 13 C sites (generated by the 18 O isotope shift) is very small, the singlet state is a good eigenstate of the coherent spin Hamiltonian. Hence a 18 O-enriched 1,2-13 C 2 squarate molecule in the nuclear singlet state remains in that state indefinitely -or more precisely, for times short compared to the relaxation time constants. But now suppose that the singlet Sq 2− molecule is protonated at time t = 0, giving a SqH − molecule. The change in electronic structure induces a chemical shift difference of ∆δ between the 13 C sites, corresponding to a precession frequency difference of ∆Ω = −γB∆δ in the magnetic field B. If the molecule remains in the protonated state for an interval τ , the spin state of the 13 C 2 pair state at time t = τ is given by</p><p>where the central triplet state is given by</p><p>Since the state probabilities are given by the square magnitudes of the superposition coefficients, the probability of the 13 C 2 pair being in the singlet state falls from unity to the value of cos 2 ∆Ωτ at time t = τ . Simultaneously the probability of the 13 C 2 pair being in the central triplet state increases from zero to sin 2 ∆Ωτ at time t = τ . Now consider a sub-ensemble of Sq 2− molecules, all protonated simultaneously at time t = 0, and all in the nuclear singlet state. We may ask what the state of this sub-ensemble is at a later time t = T , where T is long compared to the lifetime τ SqH − of a protonated squarate molecule, but sufficiently short that is unlikely that a squarate would become protonated during this time, and also short compared to the background relaxation rate constant R 0 S . The probability that one of these SqH − molecules remains protonated for an interval τ , and is then deprotonated during the following small interval δτ , is given by k d exp(−k d τ )δτ , through equation 38. The probability that an individual SqH − molecule has that particular history, and is still in the singlet state after deprotonation, is obtained by multiplying that factor by cos 2 ∆Ωτ . The probability of any member of the sub-ensemble still being in the singlet state at the end of the interval T is obtained by integrating over all possible values of the protonation duration τ . This gives</p><p>The upper limit of the integral has been extended from T to +∞ by assuming that T greatly exceeds the SqH − lifetime.</p><p>Similarly the probability that the same sub-ensemble of squarate molecules, initially in the S 0 ⟩ state and protonated at time point t = 0, are found in the central triplet state T 0 ⟩ after the interval T , is given by</p><p>The probabilities of transitions to the outer triplet states T ±1 ⟩ are negligible within the assumptions of this model. Singlet order corresponds to a difference between the population of the singlet state and the mean population of the three triplet states. The change in ensemble singlet order due to the chemical and spin dynamical processes during a finite interval T is therefore given by</p><p>The calculation above assumes that all molecules in the sub-ensemble are protonated at time point t = 0. To calculate the contribution of the protonation-deprotonation events to the relaxation rate constant for singlet order, this quantity must be multiplied by the protonation rate constant k AB , given in equation 30. This leads to the following estimate of the kinetic contribution to the relaxation of Sq 2− nuclear singlet order</p><p>If the leakage term is ignored for simplicity, an approximate expression for the singlet relaxation rate constant is given by</p><p>where the "background" singlet relaxation rate constant is given by equation 24, and the approximate form of the kinetic term is given by equation 45.</p><p>Equation 45 has a straightforward physical interpretation. The factor k p 10 −pH represents the rate constant for protonation of the Sq 2− anions. If the resulting SqH − species has a small frequency difference between the two 13 C sites (∆Ω → 0), the singlet state remains a good energy eigenstate in the SqH − species and there is no significant evolution of the singlet order before the molecule loses its proton again. The kinetic contribution to R S becomes very small in this limit, so that R S is dominated by the term R 0 S , which is the population-weighted average of the singlet relaxation rate constants in the individual species. This is the approach used previously for the studies of ligand binding by singlet NMR [30][31][32][33][34] .</p><p>If the chemical shift frequency difference in the SqH − species is large ( 2∆Ω ≫ k d ), the first term in equation 45 approaches 2 3, and the kinetic contribution to R S acquires a very simple form:</p><p>In this limit, every protonation event destroys 2/3 of the singlet order, so that the kinetic contribution to R S is simply equal to the protonation rate constant of Sq 2− anions in solution, multiplied by 2/3. Note that equation 47 is completely free of spin system parametersexcept for their role in setting the limits of the relevant approximations.</p><p>For the system considered here, the factor ∆ SO is given by ∼ 0.64 in a magnetic field of 9.4 T, i.e. slightly less than the 2 3 limit. Under these conditions, every protonation event destroys on average about 64% of the nuclear singlet order of a Sq 2− anion, and the kinetic contribution to R S is equal to about 64% of the protonation rate constant of squarate ions in solution.</p><p>Figure 13 compares pH-dependent T S values derived by the perturbative expression for R S (equation 23) with those given by the simplified approach (equations 45 and 46), as a function of magnetic field. There is reasonable agreement between these two forms of analysis.</p><!><p>The study of 13 C nuclear singlet relaxation in 18 Oenriched 1,2-13 C 2 squarate reveals a surprising dependence of the singlet relaxation time constant T S on the pH of the solution. A rather sharp transition is observed, from a short T S at pH values below ∼ 10, to long T S values above pH ∼ 12. The position of this transition is more than 6 pH points above the known pK a values of squaric acid.</p><p>A previously described model of nuclear singlet relaxation in exchanging systems, developed for ligand binding studies, assumed that the singlet relaxation rate constant is given by a population-weighted average of the relaxation rate constants for the exchanging species [30][31][32][33][34] . That model cannot explain the squaric acid observations, since the equilibrium populations of the protonated squarate species are negligible in the pH range of interest.</p><p>In this paper we show theoretically that there is an additional contribution to the singlet relaxation rate constant in chemically exchanging systems, which we call the kinetic term. This term cannot be described as a thermodynamic average over equilibrium concentrations. The kinetic term depends on kinetic rate constants and can be very significant for fast reactions such as protonation in solution.</p><p>The kinetic contribution to R S is derived for a two-site exchange model, under various approximations, and is shown to reproduce the main features of the experimental data -especially the position of the transition in singlet relaxation behaviour, which is at a pH value which is far displaced from the acid/base equilibrium point.</p><p>These observations and analysis indicate that nuclear singlet relaxation may be highly sensitive to chemical exchange, even in cases where the equilibrium overwhelmingly favours one of the exchanging species. In the current case, the concentration of the Sq 2− anion exceeds that of SqH − anion by a factor of more than 10 6 at pH = 10. Nevertheless, the rare and transient SqH − species has a very significant effect on 13 C 2 singlet relaxation. This effect presents both difficulties and opportunities. The sensitivity of nuclear singlet relaxation to rare protonation events is likely to make it even more difficult to utilize long-lived states for the storage and transport of hyperpolarized spin order, especially in protic solvents. It is probably no coincidence that the longest observed singlet lifetimes have been observed in aprotic solvents 12 , including liquid carbon dioxide 19 . It is likely that this effect has already contributed to some observations in the literature, one example being the relatively rapid decay of 13 C 2 singlet order in 18 O-labelled oxalate 76 .</p><p>On the other hand, the extreme sensitivity of nuclear singlet relaxation to rare chemical events suggest applications of this phenomenon to the determination of chemical rate constants in dynamic equilibria, under conditions beyond the range of other techniques. As shown by equation 47, nuclear singlet relaxation may provide a direct measurement of small kinetic rate constants, in suitable circumstances.</p><p>The theory given here may also be of value for other phenomena involving a combination of chemical and spin dynamics, such as Signal Amplification by Reversible Exchange (SABRE) [101][102][103] , Chemically Induced Dynamic Nuclear Polarization (CIDNP) 104,105 , and spin-isomer conversion of dissolved H 2 gas 106 .</p><p>One additional observation should be mentioned: The theory given above relies on the existence of a Larmor frequency difference between the 13 C nuclei in the protonated squarate species. This effect is field-dependent and should disappear in a sufficiently small magnetic field. However, experimental results do not fully verify this prediction. There appears to be a small pH-dependent contribution to R S which does not disappear at very low magnetic field. At the time of writing, it is not known whether this observation is due to a relatively trivial effect, such as relaxation by paramagnetic impurities, or whether it indicates a further relaxation mechanism of unknown origin.</p><!><p>There are no conflicts to declare.</p><p>Estimates of the 13 C nuclear shielding tensors for SqH − and Sq 2− , were obtained by the SCF-GIAO algorithm as implemented in Gaussian 110 . All calculations assumed a water solvation model and neglect secondary isotope shifts.</p><p>Since the Sq 2− molecule has a high symmetry, the principal axes of the symmetric parts of each shielding tensor are all precisely aligned with a molecular symmetry axis. The principal values of all four 13 C shielding tensors are the same and are calculated to be σ XX = −68.6 ppm, σ Y Y = −97.4 ppm, σ ZZ = +74.7 ppm. For each 13 C site, the X principal axis is parallel to the line joining that site with the centre of the square, the Z principal axis is perpendicular to the plane of the square, and the Y axis is perpendicular to both X and Y. Figure 14 depicts the calculated shielded tensors as blue and orange ovaloid surfaces 107,108 superposed on the Sq 2− molecular structure.  107,108 . Note the strong perturbation of the shielding tensors for the three 13 C sites closest to the protonation site, relative to Sq 2− (see figure 14).</p><p>The energy-minimized structure of the SqH − molecule has low symmetry, as shown in figure 15. Assume that the proton is attached to oxygen O1, which is bound to carbon C1. Define an axis system such that the x-axis is parallel to the line joining C1 to the centre of the square, the z-axis is perpendicular to the square, and the y-axis is perpendicular to the other two axes. The shielding tensors for the four carbon sites, expressed in this axis system, are as follows: (A1)</p><p>The symmetric parts of these tensors are depicted by the ovaloid surfaces in figure 15. The asymmetric disposition of the hydrogen atom breaks the symmetry; the numbering of the atoms is such that the H atom is closer to O4 than to O2. The calculated values of the isotropic parts of the shielding tensors are given in table A. The quantity ⟨σ iso ⟩ is the average of the isotropic 13 C shielding factors over all four positions. Since the shielding σ and chemical shift δ have opposite signs, the increase in mean isotropic shielding when Sq 2− is protonated corresponds roughly to the decrease in the isotropic chemical shift δ iso when the pH is decreased (see figure 7).</p><p>The quantity ⟨ ∆σ iso ⟩ in table A is the mean difference in isotropic shielding factor between adjacent 13 C sites, neglecting the sign. This quantity corresponds roughly to the parameter ∆δ used in the two-site calculations (table III). Table IV. Quantum chemistry estimates for the isotropic shielding parameters for the 13 C sites in Sq 2− and SqH − , assuming a water solvation model. The sites are numbered so that the hydrogen is attached to O1 and is closer to O4 than to O2. The row labelled ⟨σ iso ⟩ gives the average isotopic shielding of all four 13 C positions. The row labelled ⟨ ∆σ iso ⟩ gives the average unsigned difference in isotropic shielding factors between adjacent 13 C positions.</p><p>Appendix B: Detailed Relaxation Theory</p><!><p>Any site-specific density operators may be expanded as a superposition of orthonormal spin operators</p><p>(B1)</p><p>The coefficients q n (t) may be derived from the trace inner product q X n (t) = Tr{(Q X n ) † ρ X (t)}.</p><!><p>In general there are N L = N 2 H coefficients for each site, where N H is the dimension of the corresponding Hilbert space. The site-specific density operator is then fully characterised by the vector q X (t) = [q X 1 (t), . . . , q X N L (t)],</p><p>where boldfaced letters are used to emphasize an explicit matrix representation of any object. The matrix representation for a site-specific Liouvillian is defined via</p><p>Following equation 10 the state of the combined ensemble may be obtained by "stacking" the vectors q A (t) and q B (t) on top of each other q(t) = [q A 1 (t), . . . , q A N L (t), q B 1 (t), . . . , q B N L (t)].</p><p>(B5)</p><p>The equation of motion for q(t) is then given by d dt q(t) = Gq(t) = (L + K)q(t), (B6) with</p><p>Here, (⊕) represents a matrix direct sum, (⊗) a matrix direct product and 1 N L a unity matrix of dimension N L .</p><!><p>The generator Ĝ or its representation G may be characterised in terms of its eigenvalues and eigenvectors Gv j = Λ j v j .</p><p>(B8)</p><p>The eigenvalues Λ j are in general complex Λ j = +iω j − λ j (B9)</p><p>with ω j characterising the oscillatory behaviour and λ j ≥ 0 characterising the decaying behaviour of the eigenvector v j . The evolution of any initial density operator q(0) may then be expressed as follows q(t) = exp{Gt}q(0) = j c j (0) exp{Λ j t}v j .</p><p>(B10)</p><p>The coefficients c j (0) represent the projections of q(0) along the corresponding eigenvector v j at t = 0. In general, obtaining all eigenvalues and eigenoperators of G is difficult, and one typically has to resort to an approximate treatment of the problem. In general however the exchange rates k BA and k AB may be highly asymmetric, which poses a problem for typical approximation techniques. But we may symmetrise the exchange matrix K via a similarity transformation</p><p>where w represents the geometric mean of the forward and backward exchange rates 111 . The transformed generator of motion then takes the form</p><p>It is easy to see that G and G have the same eigenvalues. We therefore omit the tilde from now on, to keep the notation somewhat compact.</p><p>A SpinDynamica 109 calculation shows that the relevant subspace of G is spanned by the following ten operators</p>