Molecular dynamics simulations of ^1H NMR relaxation in Gd^{3+}--aqua
Philip M. Singer, Yunke Liu, Lawrence B. Alemany, George J. Hirasaki, Walter G. Chapman, Dilip Asthagiri
aa r X i v : . [ phy s i c s . c h e m - ph ] F e b Molecular dynamics simulations of H NMR relaxation in Gd –aqua Philip M. Singer, ∗ Yunke Liu, Lawrence B. Alemany, GeorgeJ. Hirasaki, Walter G. Chapman, and Dilip Asthagiri † Department of Chemical and Biomolecular Engineering,Rice University, 6100 Main St., Houston, TX 77005, USA Shared Equipment Authority and Department of Chemistry,Rice University, 6100 Main St., Houston, TX 77005, USA
Atomistic molecular dynamics simulations are used to investigate H NMR T relaxation of wa-ter from paramagnetic Gd ions in solution at 25 ◦ C. Simulations of the T relaxivity dispersionfunction r computed from the Gd – H dipole–dipole autocorrelation function agree within ≃ f & r above f & B & r in chelated Gd contrast-agents used for clinical MRI, withoutany adjustable parameters or models. Below f . r comparedto measurements, which is used to estimate the zero-field electron-spin relaxation time. The moststrongly bound water molecule to the Gd complex evinces an autocorrelation function consistentwith the mono-exponential decay used in the Solomon-Bloembergen-Morgan (SBM) inner-spheremodel. I. INTRODUCTION
The traditional theory of enhanced H NMR (nuclearmagnetic resonance) relaxation of water due to paramag-netic transition-metal ions and lanthanide ions in aque-ous solutions originates from Solomon [1], Bloembergenand Morgan [2, 3], a.k.a. the Solomon-Bloembergen-Morgan (SBM) model. The extended SBM model ac-counts for paramagnetic relaxation of both inner-sphere[1–3] and outer-sphere [4, 5] water of the paramagnetic-ion complex, the contact term [2, 3, 6, 7], the Curie term[5, 8], and the electron-spin relaxation [3, 5, 8–12]. TheSBM model is most widely used in the interpretation ofparamagnetic enhanced H relaxation of water due toGd -based contrast agents [13–34] used in clinical MRI(magnetic resonance imaging). The SBM model is alsoused in the interpretation of paramagnetic relaxation inwater-saturated sedimentary rocks [35–41].The inner-sphere water constitutes the ligands of theGd complex. The SBM inner-sphere model assumes arigid Gd – H dipole-dipole pair undergoing rotationaldiffusion, which according to the Debye model resultsin mono-exponential decay of the autocorrelation func-tion. The Debye model was previously used in theBloembergen-Purcell-Pound (BPP) model [42] for likespins (e.g. H– H pairs), and then adopted in the SBMmodel for unlike spins (e.g. Gd – H pairs). The SBMinner-sphere model also takes the electron-spin relaxationtime into account, resulting in 6 free parameters for the H NMR relaxivity r . The inner-sphere relaxation isgenerally considered the largest contribution to the totalrelaxivity.The second-sphere and outer-sphere water are lesstightly bound than the inner-sphere water. The extended ∗ [email protected] † [email protected] SBM outer-sphere model assumes that the Gd ion andH O are two force-free hard-spheres undergoing transla-tional diffusion, which according to the Torrey [43] andHwang-Freed [44] models results in stretched-exponentialdecay of the autocorrelation function. The SBM outer-sphere model adds an additional 2 free parameters, bring-ing the total to 8 free parameters for the extended SBMmodel. Relaxation from the contact term is negligiblefor H NMR in Gd –aqua [20, 27], and is therefore ne-glected. Likewise, the Curie term is negligible in thepresent case [5, 8].The application of the extended SBM model to Gd –aqua therefore requires fitting 8 free parameters overa large frequency range in measured r dispersion(NMRD). In chelated Gd complexes, an order param-eter plus a shorter correlation time [45, 46] are addedto the rotational motion of the complex [32], or to theelectron-spin relaxation [19, 26, 29], taking the total to 10free parameters. Not only is this an over-parameterizedinversion problem, but some of the underlying assump-tions are questionable, especially regarding electron-spinrelaxation time [4].Atomistic molecular dynamics (MD) simulations canhelp elucidate some of these issues. MD simulation werepreviously used for like-spins, e.g. H– H dipole-dipolepairs, such as liquid-state alkanes, aromatics, and water[47–50], as well as methane over a large range of densities[51]. In all these cases, good agreement was found be-tween simulated and measured H NMR relaxation anddiffusion, without any adjustable parameters in the inter-pretation of the simulations. With the simulations thusvalidated against measurements, simulations can then beused to separate the intramolecular (i.e. rotational) fromintermolecular (i.e. translational) contributions to relax-ation, and to explore the corresponding H– H dipole-dipole autocorrelation functions in detail. For instance,MD simulations revealed for the first time ever that wa-ter and alkanes do not conform to the BPP model of amono-exponential decay in the rotational autocorrelationfunction, except for highly symmetric molecules such asneopentane. More complex systems such as viscous poly-mers [52] and heptane confined in a polymer matrix [53]have also been simulated, which again saw good agree-ment with measurements, and which lead to insights intothe distribution in dynamic molecular modes.In this report, we extend these MD simulation tech-niques to Gd – H dipole-dipole pairs, i.e. unlike spins,in a Gd –aqua complex. The simulations show goodagreement with measurements above f & r in chelated Gd contrast agents used for clinical MRI,without relying on the over-parameterized extended SBMmodel. At the very least, these results could reduce thenumber of free parameters in the SBM model, and helpput constraints on its inherent assumptions for clinicalMRI applications. II. METHODOLOGYA. Molecular simulation
To model the Gd –aqua system, we use theAMOEBA polarizable force field to describe solvent wa-ter [54] and the ion [55]. The simulations are conductedusing the Tinker-8 package [56]. The Gd ion has anincomplete set of 4 f orbitals, but this incomplete shelllies under the filled 5 s and 5 p orbitals. Thus ligand-fieldeffects [57] are not expected to play a part in hydrationand a spherical model of the ion ought to be adequate.However, polarization effects would be important giventhe large charge on the cation, and hence our choice ofthe AMOEBA model.Experimental NMR studies on Gd –aqua use concen-trations of about 0.3 mM. This amounts to having oneGdCl molecule in a solvent bath in a cubic cell about176 ˚A in size. Such large simulations are computationallyrather demanding with the AMOEBA polarizable force-field. Further, at such dilutions, the ions essentially “see”only the water around them, not the other ions, and sinceunderstanding the behavior of water around the ion is offirst interest, here we study a single Gd ion in a cu-bic cell of size about 24 . − Debyes was maintained in the polarizationcalculations. The equation of motion was integrated us-ing the RESPA integrator [58, 59] with a time step of1 fs.In our exploratory calculations, we had first equili- brated the system for about 250 ps under
N pT con-ditions, where the pressure (of 1 atm.) and temper-ature (of 298 K) was maintained using a Nos´e-Hoover[60, 61] barostat and thermostat, respectively. These cal-culations predict an ion partial molar volume of about −
126 cc/mol, which can be compared with the experi-mental value of − . − . .
805 ˚A. This system was again equili-brated under
N V T conditions for over 200 ps now usingthe Andersen thermostat [63], which we have found to bea better first step in preparing the system for
N V E sim-ulations [47]. Subsequently, we propagated the systemunder
N V E conditions for 1.81 ns, saving frames every0.1 ps for analysis. We used the last 1.6384 ns of simula-tion (equal to 16384= 2 frames) for analysis. The meantemperature in the production phase was 297 . N V E phase is excellent, with the standard error relative to themean value being 10 − .We also performed exploratory calculations using onemolecule of GdCl in 512 water molecules. The broadertrends using the GdCl simulations were similar, howeverscatter in the G ( t ) data was larger, suggesting potentialinfluence of the Cl − ions.
1. Mean residence time analysis
Several strategies [64–66] have been proposed in theliterature to estimate the mean residence time of wa-ter molecules in the inner-hydration sphere of ions andbiomolecules. All these methods rely on keeping track ofthe particle as it moves from the defined inner sphere tothe second sphere or the bulk.The procedure due to Impey et al. [64] is perhapsthe most used. They define a survival probability (seebelow), but base this on an ad hoc choice of a coarse-graining time ( t ⋆ = 2 ps); transient escapes from theinner sphere that are shorter than t ⋆ are ignored. Theso-called “direct” procedure by Rode et al. [65] scansthe entire trajectory and counts all exchange events thatpersist for longer than t ⋆ = 0 . t = 0 is found in the reactant stateat time t . By construction this method is insensitive tobarrier re-crossings. But for the case of cation-hydration,the barrier separating the inner-sphere and the secondsphere will be high (as suggested by the pair-correlation)and accounting for re-crossings is likely not critical.Since our interest is in a qualitative estimate of theresidence time, following Garc´ıa and Stiller [67] we adoptthe Impey et al. [64] procedure, but with the coarseningtime set equal to the time separation between consecutiveframes, δt (=100 fs here).We follow the residence time of some defined watermolecule w using an indicator function χ w that equals1 if the water molecule is in the inner sphere and zerootherwise. For each saved configuration, we keep trackof χ w and then construct the autocorrelation to half thetotal length of time. The i th element of the autocorrela-tion vector is just the number of times the water moleculesurvives in the inner sphere for i consecutive time inter-vals. From this information, we construct the probability P ( t ) that w survives in the inner sphere for i = 1 , , . . . time intervals and t = iδt .Following Impey et al. , we model the initial decay ofthe survival probability by P ( t ) = P o exp (cid:18) − tτ m (cid:19) (1)where P o is a constant and τ m is defined as the meanresidence time. The results below show that the expo-nential model is good for times up to 100 ps. Both thelimitation of the total simulation time and the possibilityof water molecules escaping to the bulk without ever re-turning limit the applicability of this model to very longtimes. B. Structure and dynamics
Fig. 1 shows the ion-water radial distribution function.The location and magnitude of the peak is in good agree-ment with earlier studies founded on either ab initio [27]or empirical force field-based [29] simulations. Consistentwith these studies, we find that the first sphere (i.e. in-ner sphere) contains between q = 8 ↔ q ≃ r = 3 ˚A, Fig. 1),second shell (from 3 ˚A to 5.5 ˚A), and bulk domains (be-yond 5.5 ˚A). Fig. 2 shows the 32 most bound moleculesin order of decreasing time spent in the inner sphere. Wefind that 3 water molecules spend more than 50% of thetime in the inner sphere, called “group A” molecules, andthese 3 water molecules contribute a total of 2.5 to theaverage coordination number. The remaining 41 − q − . FIG. 1. Radial distribution function g ( r ) of water oxygenaround the Gd ion; the function n ( r ) gives the coordinationnumber, for both total and group A partitioned molecules.FIG. 2. Percent of time spent by water molecules that visitthe inner-sphere of Gd during the 1.64 ns production runof the simulation. Of those water molecules, N A = 3 (groupA) spend a large fraction (68% ↔ . N A = 3 most strongly bound molecules are partitionedinto group A. characterized by fewer water molecules than indicatedby the average coordination number. Our finding of 3( < q = 8 . ion is consistent with this ex-pectation. These molecules are also the ones that oughtto most influence the NMR relaxation.Fig. 3 shows the behavior of the survival probability FIG. 3. Survival probability for various water molecules. Thedata is fit to Eq. 1 between time 0 and the 100 ps boundaryindicated by the dashed vertical line. The residence times τ m are indicated in the legend. for water molecules, including group A molecules. Wecan immediately notice that the residence times for theindividual inner-sphere water molecules spans time scalesfrom the O (100 ns) for the most well bound water to O (100 ps) for the least well bound water. For watermolecule τ m ≃ O NMR data, where τ m ≃ . ↔ . C. H NMR relaxivity
The enhanced H NMR relaxation rate 1 /T for wateris given by the average over the N = 512 water moleculesin the L = 24 .
805 ˚A box containing one paramagneticGd ion: 1 T = 1 N N X i =1 T i , (2)1 T = 1 N N X i =1 T i , (3)where T i is the average T relaxation for the i ’th watermolecule. In other words, T i is averaged over the Hpair on the i ’th water molecule. The gadolinium molarconcentration is given by [ Gd ] = [ W ] /N , where [ W ] =55,705 mM is the molar concentration of water in thesimulation box, which leads to the following expression for the relaxivity: r = 1[ Gd ] 1 T = 1[ W ] N X i =1 T i , (4) r = 1[ Gd ] 1 T = 1[ W ] N X i =1 T i , (5)in units of (mM − s − ), which is independent of N (or L equivalently). The “fast-exchange” regime ( T i, i ≫ τ m ) is assumed [9], see below. The H– H dipole-dipolerelaxation [47] is not considered in these simulations.The computation of T , originates from the Gd – Hdipole-dipole autocorrelation function G ( t ) [42, 43, 71–74] shown in Fig. 4(a), where t is the lag time. Thisautocorrelation is well suited for computation using MDsimulations [75, 76]. Using the convention in the text byMcConnell [72], G i ( t ) in units of (s − ) for the i ’th watermolecule is determined by: G i ( t ) = 14 (cid:16) µ π (cid:17) ~ γ I γ S S ( S + 1) × X j =1 * (3 cos θ ij ( t + t ′ ) − r ij ( t + t ′ ) (3 cos θ ij ( t ′ ) − r ij ( t ′ ) + t ′ , (6)where µ is the vacuum permeability and ~ is the reducedPlanck constant. γ I / π = 42 .
58 MHz/T is the nucleargyro-magnetic ratio for H with spin I = 1 /
2, and γ S =658 γ I is the electron gyro-magnetic ratio for Gd withspin S = 7 /
2. The summation j = 1,2 is over the H pairon the i ’th water molecule, and the factor 1/2 outside thesummation emphasizes that G i ( t ) is the average value forthe i ’th water molecule. This convention is used since1 /T i in Eq. 2 refers to the average relaxation of the i ’th water molecule, and [ W ] = 55,705 mM is the molarconcentration of water in the simulation box. θ ij is theangle between the Gd – H vector r ij and the appliedmagnetic field B , where { ij } refers to the j ’th H onthe i ’th water molecule. For convenience, the i index isordered by the amount of time the i ’th water moleculespends in the inner sphere, as shown in Fig. 2.We assume a spherically symmetric (i.e. isotropic)system, and therefore G m ( t ) = G ( t ) is independent ofthe order m , which amounts to saying that the direc-tion of the applied magnetic field B = B z is arbitrary.This assumption was verified in Ref. [29]. For simplic-ity, we therefore use the m = 0 harmonic Y ( θ, φ ) = p / π (3 cos θ −
1) for the MD simulations [47].The second-moment ∆ ω i (i.e. strength) of the dipole-dipole interaction is defined as such [73]: G i (0) = 13 ∆ ω i . (7)Assuming the angular term in Eq. 6 is uncorre-lated with the distance term at t = 0, the relation (cid:10) (3 cos θ ij ( τ ) − (cid:11) τ = 4 / FIG. 4. (a) MD simulations of the autocorrelation functions G ( t ), where 1 in 10 data points are shown for clarity, for to-tal and inner-sphere (scaled by a factor 1/2 on the y -axis forclarity). Also shown are ILT (inverse Laplace transform) fits(Eq. 23) to the simulated data, and SBM model (Eq. 32).(b) Spectral density functions J ( ω ) from FFT (fast Fouriertransform) (Eq. 17), including the f = 0 data point repre-sented as a horizontal line placed at low frequency, for totaland inner-sphere (scaled by a factor 1/3 on the y -axis forclarity). Also shown is the ILT result (Eq. 28), and the f − behavior expected in the SBM model at high frequencies. indices ij ) reduces the second moment to:∆ ω i = 35 (cid:16) µ π (cid:17) ~ γ I γ S S ( S + 1) 12 X j =1 * r ij ( t ′ ) + t ′ . (8)The time-averaged Gd – H distances r i for the i ’th wa-ter molecule can also be determined as such:1 r i = 12 X j =1 * r ij ( t ′ ) + t ′ , (9) the results of which are shown in Fig. 5 for the first 32molecules.In accordance with Eq. 2, the summation over N =512 water molecules is then taken to define the totalquantities: G ( t ) = N X i =1 G i ( t ) , (10) h τ i = 1 G (0) Z ∞ G ( t ) dt (11)∆ ω = N X i =1 ∆ ω i . (12)In practice, only the 41 water molecules which visit theinner sphere are required in the summation since the re-maining 471 have negligible contributions to G ( t ). Theaverage correlation time h τ i is defined as the normalizedintegral [73].Of additional interest is the contribution from inner-sphere water alone. According to Fig. 2, molecule ≃ q = 8 . G in ( t ) = q G ( t ) . (13) τ in = 1 G in (0) Z ∞ G in ( t ) dt (14)∆ ω in = q ∆ ω (15) r in = r , (16)where G in ( t ) represents the full inner-sphere contribu-tion. Also defined is the average Gd – H distance r in for molecule in impliesinner-sphere (Eq. 13).The next step is to take the (two-sided even) fastFourier transform (FFT) of the G ( t ) to obtain the spec-tral density function: J ( ω ) = 2 Z ∞ G ( t ) cos( ωt ) dt, (17)the results of which are shown as the data symbols inFig. 4(b). The relaxation rates are then determined forunlike spins [72]:1 T = J ( ω ) + 73 J ( ω e ) , (18)1 T = 23 J (0) + 12 J ( ω ) + 136 J ( ω e ) , (19)assuming ω e ≫ ω , where ω = γ I B = 2 πf is the H NMR resonance frequency, and ω e = 658 ω is the FIG. 5. Time-averaged Gd – H distances r i for the i ’thwater molecule according to Eq. 9, where only the first N =32 molecules are shown, in the order presented in Fig. 2,where r = r in . electron resonance frequency. Finally, the relaxivity dis-persion functions are given by: r = 1[ W ] 1 T , (20) r = 1[ W ] 1 T , (21)where the chemical shift term in r [5] is neglected forsimplicity. We also assume the “fast-exchange” regime,namely that the condition T , ≫ τ m holds [9]. This isclearly the case here given that T , & µ s accordingto Eqs. 18, 19, Fig. 4(b), and given the residence time τ m ≃ O NMR T mea-surements [14, 17, 18, 24, 70]. The fast-exchange regimecan also be inferred directly from measurements since r increases with decreasing temperature [13]. Investi-gations are underway to extend the simulations to theslow-exchange regime.The low frequency (i.e. extreme narrowing) limit r (0) = r (0) = r , (0) is given by: r , (0) = 1[ W ] 209 ∆ ω h τ i . (22)Note how r , (0) (for unlike Gd – H spin pairs) is afactor 2/3 less than the equivalent expression for like H– H spin pairs [47].The above analysis also assumes the electron spin is apoint-dipole centered at the Gd ion [4]. Given that thesimulation agrees with measurements above f & H [5].
D. Inverse Laplace Transform
The FFT result for J ( ω ) in Fig. 4(b) is sparse. Be-sides the f = 0 data point, the lowest frequency FFTdata point is given by the resolution ∆ f = 1 / t max =2440 MHz, where t max = 204 ps is the longest lag timecomputed in G ( t ). Zero padding to increase t max andtherefore reduce ∆ f adds oscillations to J ( ω ), and there-fore does not improve the situation.In order to overcome these limitations in the FFT, weuse the inverse Laplace transform (ILT) of G ( t ) instead,the results of which are shown in Fig. 4(b). The followingadvantages are recognized: (1) the ILT filters out thenoise while still honoring the FFT data (including the f = 0 data point), (2) the ILT provides J ( ω ) for anydesired f value, and (3) the ILT leads to physical insightinto the distribution P ( τ ) of molecular correlation times τ . FIG. 6. Probability density function P ( τ ) obtained from theILT (Eq. 23), for total and inner-sphere (scaled by a factor1/6 on the y -axis for clarity). Correlation times h τ i and τ in (Table I), and τ T (Table II) also listed. The distribution P ( τ ) is determined from the ILT of G ( t ) [77, 78]: G ( t ) = Z ∞ P ( τ ) exp (cid:18) − tτ (cid:19) dτ. (23)The goal is to recover the distribution in molecular cor-relation times P ( τ ) from Eq. 23, which is a Fredholmintegral equation of the first kind. For concision, andbecause this terminology is in wide-spread use, we termthe procedure of recovering P an “inverse Laplace trans-form,”(ILT) but we emphasize that strictly speaking in-verting Eq. 23 to recover P is not formally a Laplaceinversion [79]. With this understanding, we note thatdetails of the ILT procedure can be found in Refs. [50–53] and the supplementary material in Refs. [49, 80]. Webriefly highlight the essential details here. G ( t ) is available only at discrete time intervals, and theinversion is a rather ill-posed problem. In our approach,we use Tikhonov regularization [49, 80], with the vector P being one for which P = arg min P ≥ || G − K P || + α || P || (24)is a minimum. Here G is the column vector representa-tion of the autocorrelation function G ( t ), P is the columnvector representation of the distribution function P ( τ ), α is the regularization parameter, and K is the kernelmatrix: K = K ij = exp (cid:18) − t i τ j (cid:19) . (25)The results for P ( τ ) are shown in Fig. 6, from which thefollowing are determined: h τ i = 1 G (0) Z ∞ P ( τ ) τ dτ, (26) G (0) = Z ∞ P ( τ ) dτ = 13 ∆ ω . (27)The spectral density J ( ω ) is then determined from theFourier transform (Eq. 17) of G ( t ) (Eq. 23): J ( ω ) = Z ∞ τ ωτ ) P ( τ ) dτ , (28)from which T , at any desired f can be determined (Eqs.18 and 19). Some alternatives to the above ILT formula-tion are discussed in Ref. [50].The P ( τ ) is binned from τ min = 10 − to τ max = 10 ps using 200 logarithmically spaced bins. τ max is chosensuch that τ max = 5 t max , consistent with previous studies[50, 52, 53]. In the present case of low-viscosity fluids( η ≃ τ max does not impact J ( ω ),and is therefore not a free parameter in the ILT analysis.The constant “div” in Fig. 6 is a “division” on a log-scale. More specifically, div = log ( τ i +1 ) − log ( τ i ) isindependent of the bin index i , and ensures unit areawhen P ( τ ) is of unit height and a decade wide [80].As shown in Fig. 4(a), the residual between the G ( t ) data and the ILT fit is not dominated by Gaus-sian noise. As such, the regularization parameter is fixedto α = 10 − in accordance with previous studies [49–53]. As shown in Fig. 4(b), we find that α = 10 − givesgood agreement with the parameter-free J ( ω ) from FFT,which validates the ILT analysis. E. Diffusion
An independent computation of translational diffusion D T was performed from MD simulations. We calculatethe mean square displacement (cid:10) ∆ r (cid:11) of the water oxygenand Gd ion as a function of the diffusion evolution time FIG. 7. MD simulations of mean-square displacement (cid:10) ∆ r (cid:11) versus time t for bulk water, water in Gd –aqua, and Gd in Gd –aqua. Solid lines show fitting region used to obtaintranslational diffusion coefficient D sim from Eq. 29 for t > D T values which include the correctionterm (Eq. 30). t ( <
10 ps), and additionally average over a sample of 50molecules to ensure adequate statistical convergence.As shown in Fig. 7, at long-times ( t ) the slope ofthe linear diffusive regime gives the translational self-diffusion coefficient D sim according to Einstein’s relation: D sim = 16 δ (cid:10) ∆ r (cid:11) δt (29)where D sim is the diffusion coefficient obtained in thesimulation using periodic boundary conditions in a cu-bic cell of length L . In the linear regression procedure,the early ballistic regime and part of the linear regime isexcluded to obtain a robust estimate of D sim .Following Yeh and Hummer [81] (see also D¨unweg andKremer [82]), we obtain the diffusion coefficient for aninfinite system D T from D sim using D T = D sim + k B T πη ξL (30)where η is the shear viscosity and ξ = 2 . ξ/L is the potentialat the charge site in a Wigner lattice.) For the systemsizes considered in this study, L ≃
25 ˚A, the correctionfactor constituted ≃
13% of D , and ≃
16% of D W . Thecorrection factor was not applied to D Gd . F. Measurements
We prepared a Gd –aqua solution in de-ionized waterat [ Gd ] = 0.3 mM and measured T ,meas at a controlledtemperature of 25 ◦ C, using static fields at f = 2 . f = 500MHz with a Bruker 500 MHz Spectrometer. The mea-sured relaxivity r meas was determined as such: r meas = 1[ Gd ] (cid:18) T ,meas − T ,bulk (cid:19) . (31)where T ,bulk = 3.13 s was found for bulk water (notde-oxygenated [83]) at 2.3 MHz and 500 MHz. Fieldcycling r meas data of a Gd –aqua solution at 25 ◦ C weretaken from Luchinat et al. [30] (supplementary material),which agreed with our measurements at f = 2.3 MHz(within ≃ f = 500 MHz significantlyextends the frequency range of the measurements. III. RESULTS AND DISCUSSIONS
In this section we compare the simulated total relax-ivity r with measurements. The simulated relaxivity r in for inner-sphere water is then compared with theSolomon-Bloembergen-Morgan (SBM) model. The zero-field electron-spin relaxation time is then incorporated tomatch the low-frequency relaxivity measurements. A. Comparison with measurements
The results of total simulated relaxivity r (Eq. 20) areshown in Fig. 8, alongside corresponding measurementsof Gd –aqua solution at 25 ◦ C by Luchinat et al. [30](supplementary material) and this work. A cross-plot ofsimulated r (total) versus measured r meas results arealso shown in Fig. 9. The simulated r dispersion iswithin ≃
5% of measurements r meas above f & r and r meas above f & ≃ h τ i (Eq. 26) and the square-root of the secondmoment (i.e. strength) of ∆ ω (Eq. 27) for the total.Note that the product ∆ ω h τ i ≃ .
006 indicates that theRedfield-Bloch-Wagness condition (∆ ω h τ i ≪
1) is sat-isfied [71–73], which justifies the relaxivity analysis usedhere.
B. Comparison with SBM model
The SBM inner-sphere model assumes a rigid Gd – H pair undergoing rotational diffusion, leading to the
FIG. 8. H NMR relaxivity measurements r meas of Gd –aqua solution at 25 ◦ C (from Ref. [30] supplementary mate-rial, and this work), compared with simulated total relaxivity r , and inner-sphere r in . Also shown is the simulated r e pre-dicted at low frequency ( f < . T e ≃
153 ps taken into account(Eq. 36).FIG. 9. Cross-plot of simulated r (total) versus measured r meas at 25 ◦ C for data taken from Fig. 8, both above f > • ) and below f < ◦ ). Also shown is thesimulated r e predicted at low frequency ( f < . T e ≃
153 ps takeninto account (Eq. 36) following mono-exponential decay in the autocorrelationfunction: G SBM ( t ) = G SBM (0) exp (cid:18) − tτ R (cid:19) . (32) h τ i ∆ ω/ π q τ in ∆ ω in / π r in T e (ps) (MHz) (ps) (MHz) (˚A) (ps)36 25.4 8.3 26 23.6 3.09 153TABLE I. Analysis of the simulation results including; (left)total correlation time h τ i (Eq. 26) and square-root of secondmoment ∆ ω (Eq. 27); (middle) inner-sphere coordinationnumber q (Fig. 1), correlation time τ in (Eq. 26), square-rootof second moment ∆ ω in (Eq. 27), and Gd – H distance r in (Fig. 5); (right) zero-field electron-spin relaxation time T e (Eq. 36). This functional form is identical to the BPP model [42]which is based on the Debye model, where the rotational-diffusion correlation time τ R is defined as the averagetime it takes the rigid pair to rotate by 1 radian.As shown in Fig. 4(a), G in ( t ) for molecule G SBM ( t ) (plottedwith τ R = τ in = 26 ps) over the first decade of decay, con-sistent with previous findings [29]. This result is expectedgiven that molecule ion. The SBM inner-sphere model also predicts a f − decrease in the spectraldensity J ( ω ) at high frequencies, which is consistent with J in ( ω ) shown in Fig. 4(b). This is also reflected in Fig.6 where P in ( τ ) for inner-sphere is dominated by a sharp(almost single-valued) peak at τ ≃
26 ps.Note that as expected, the scattering in J in ( ω ) (simu-lation of molecule J ( ω ) (simulation averaged over allmolecules) from FFT. However the ILT analysis filtersout the scattering while honoring the FFT data, includ-ing the f = 0 data point.The inner-sphere correlation time τ in = 26 ps (Table I)is consistent with published values from the SBM inner-sphere model, where a range of τ R ≃ ↔
45 ps is found[13, 15, 19, 20, 22, 29, 30]. Note also that τ in = 26 ps isa factor ≃
10 larger than that of bulk water ( τ R ≃ . –aqua complex. Also listed in TableI is the average Gd – H distance r in (Eq. 16) for inner-sphere water (Fig. 5). The average r in ≃ .
09 ˚A isconsistent with published values from the inner-sphereSBM model, where a range of r ≃ . ↔ . r in is shown in Fig. 8,which accounts for ≃
62% of the total r at all f .This is a smaller contribution than previously published[13, 15, 19, 20, 22, 29, 30], however note that we deter-mine r in and r without any adjustable parameters ormodels. The remaining ≃
38% contribution to r comesfrom second/outer-sphere water, which is a larger contri-bution than previously published.Fig. 4(a) shows that the total G ( t ) deviates from Eq.32 (plotted with τ R = h τ i = 36 ps). The deviation of G ( t )from Eq. 32 is expected since G ( t ) includes the inner-sphere plus the second/outer-sphere contributions. The second/outer-sphere contributions are expected to followthe Hwang-Freed (HF) model for the relative transla-tional diffusion between Gd and H O assuming twoforce-free hard-spheres [44]: G HF ( t ) = G HF (0) 54 π ∞ Z x
81 + 9 x − x + x × exp (cid:18) − x tτ D (cid:19) dx. (33)The translational-diffusion correlation time τ D is definedas the average time it takes the molecule to diffuse by onehard-core diameter d . G HF ( t ) is a multi-exponential de-cay by nature, and therefore one expects the total G ( t ) tobe stretched, the extent of which depends on the relativecontributions of second/outer-sphere to inner-sphere.The correlation time τ D can be predicted as such: τ D = d D W + D Gd = 94 τ T , (34)where the simulated diffusion coefficients of H O ( D W )and Gd ( D Gd ) in Gd –aqua are taken from Fig. 7,the results of which are listed in Table II. The hard-corediameter d is taken from the local maximum at ≃ . g ( r ) in Fig. 1 (which isattributed to second-sphere water). The resulting τ D ≃
91 ps is a factor ≃
10 larger than that of bulk water( τ D ≃ . d ≃ . d ≃ D W D Gd d τ D τ T (˚A /ps) (˚A /ps) (˚A) (ps) (ps)0.20 0.03 4.6 91 40TABLE II. Diffusion coefficients of H O ( D W ) and Gd ( D Gd ) in Gd –aqua, distance of closest approach ( d ) betweenGd and second-sphere H O according to g ( r ) (Fig. 1),and, resulting translational-diffusion correlation time τ D (=9 / τ T ) from Eq. 34. The value τ D ≃
91 ps can be compared with the distri-bution P ( τ ) in Fig. 6. More specifically, the translationalcorrelation time τ T ≃
40 ps is plotted in Fig. 6, wherethe relation τ D = 9 / τ T and the origin of the factor 9/4is explained in [49, 73]. The width of the main peak in P ( τ ) can therefore be accounted for by the sharp peak at τ in = 26 ps from inner-sphere water, plus a broad peakat τ T = 40 ps from second/outer-sphere water.We note that P ( τ ) in Fig. 6 has a small contribution atshort τ ≃ − ps, which is a result of the sharp drop in G ( t ) over the initial t ≃ . –aqua, which can perhaps be explained0using a two rotational-diffusion model such as found inbulk water [84].The r dispersion in Fig. 8 results in a mild increase inthe ratio T /T = r /r ≃ / f &
10 MHz, until f & T /T increases further. CombiningEqs. 18 and 19 with J ( ω e ) = 0 (i.e., slow-motion regime)and J ( ω ) = J (0) (i.e., fast-motion regime) accounts forthe factor T /T = 7/6 within the frequency range f ≃
10 MHz ↔ T /T = 7 / r in . C. Electron-spin relaxation
The difference between r meas and r below f . T e = T e (0) = T e (0) into account. As-suming that the correlation times P ( τ ) are uncorrelatedwith the electron-spin relaxation times, the following ex-pression results [3]: r e , (0) = 1[ W ] 209 ∆ ω h τ ′ i , (35)1 h τ ′ i = 1 h τ i + 1 T e . (36)This is equivalent to introducing an exponential decayterm exp( − t/T e ) inside the FFT integral of Eq. 17. Thefitted value of T e ≃
153 ps is determined by matching r e below f . . r meas . Theresulting T e ≃
153 ps is consistent with the publishedrange of T e ≃ ↔
160 ps [13, 15, 19, 20, 22, 29, 30].Investigations are underway to incorporate T e ( ω e ) and T e ( ω e ) dispersion [4, 22, 24, 29] over the full frequencyrange in r meas . IV. CONCLUSIONS
Atomistic MD simulations of H NMR relaxivity r forwater in Gd –aqua complex at 25 ◦ C show good agree-ment (within ≃ f & – H dipole-dipole relaxation for unlike spins. The resultsabove f & B & r in chelated Gd contrast-agents for clinical MRI, without relying on theSolomon-Bloembergen-Morgan (SBM) model and its 8(or more) free parameters.Water molecules more strongly bound to the inner-sphere are partitioned into “group A” based on timespent in the inner-sphere of the Gd complex during the 1.64 ns production run of the simulation. A resi-dence time of the order τ m ≃ .
28 ns is found for watermolecules expected to set the time scale for the rejuvena-tion of the inner-sphere population, which is consistentwith published values.The autocorrelation function G in ( t ) for inner-spherewater is consistent with the mono-exponential decayused in the SBM inner-sphere model. Furthermore, therotational-diffusion correlation time τ in ≃
26 ps for theGd –aqua complex, the average Gd – H separation r in ≃ .
09 ˚A, and the inner-sphere coordination q ≃ . g ( r ), areall consistent with previously published values whichused the SBM inner-sphere model.The total autocorrelation function G ( t ) shows a mildlystretched-exponential decay, with an average correlationtime h τ i ≃
36 ps. The mildly stretched-exponentialnature of G ( t ) is expected given that it includes bothinner-sphere and (stretched) second/outer-sphere contri-butions. A distance of closest approach (i.e. hard-corediameter) for Gd –H O of d ≃ . g ( r ) (attribute to second-sphere water), which together with the simulated diffu-sion coefficients of Gd and H O are used to estimatethe translational-diffusion correlation time τ D ≃
91 ps.The force-free Hwang-Freed (HF) model of second/outer-sphere relaxation, and the estimated value for τ D , helpexplain the distribution in correlation times P ( τ ) deter-mined from the inverse Laplace transform of G ( t ).The total relaxivity r at all f consists of ≃
62% con-tribution from inner-sphere water r in , plus a ≃
38% con-tribution from second/outer-sphere water. The contribu-tion from inner-sphere is lower than previously published,however note that we determine r in and r without anyadjustable parameters or models.Below f . r compared to measurements, which is reconciled by tak-ing the zero-field electron-spin relaxation time T e intoaccount. The resulting fitted value T e ≃
153 ps is con-sistent with published values.
ACKNOWLEDGMENTS
We thank Vinegar Technologies LLC, Chevron EnergyTechnology Company, and the Rice University Consor-tium on Processes in Porous Media for financial support.We gratefully acknowledge the National Energy ResearchScientific Computing Center, which is supported by theOffice of Science of the U.S. Department of Energy (No.DE-AC02-05CH11231) and the Texas Advanced Com-puting Center (TACC) at The University of Texas atAustin for high-performance computer time and support.We also thank Arjun Valiya-Parambathu, Xinglin Wang,and Thiago Pinheiro for helpful discussions.1 [1] Solomon, I. Relaxation processes in a system of two spins
Phys. Rev.
99 (2), 559–565[2] Bloembergen, N. Proton relaxation times in paramag-netic solutions
J. Chem. Phys.
27 (2), 572–573[3] Bloembergen, N.; Morgan, L. O. Proton relaxation timesin paramagnetic solutions. Effects of electron spin relax-ation
J. Chem. Phys.
34, 842–850[4] Kowalewski, J.; Nordenski¨old, L.; Benetis, N.; West-lund, P.-O. Theory of nuclear spin relaxation on para-magnetic systems in solution
Prog. Nucl. Magn. Reson.Spect.
17, 141–185[5] Helm, L. Relaxivity in paramagnetic systems: Theoryand mechanisms
Prog. Nucl. Magn. Reson. Spect.
49, 45–64[6] Morgan, L. O.; Nolle, A. W. Proton spin relaxationin aqueous solutions of paramagnetic ions, II. Cr +++ ,Mn ++ , Ni ++ , Cu ++ , and Gd +++ J. Chem. Phys.
31 (2), 365–368[7] Bernheim, R. A.; Brown, T. H.; Gutowsky, H. S.; Woess-ner, D. E. Temperature dependence of proton relaxationtimes in aqueous solutions of paramagnetic ions
J. Chem.Phys.
30 (4), 905–956[8] Fries, P. H.; Ferrante, G.; Belorizky, E.; Rast, S. Therotational motion and electronic relaxation of the Gd(III)aqua complex in water revisited through a full protonrelaxivity study of a probe solute
J. Chem. Phys.
Phys. Rev.
Dynamics of Solutions and Fluid Mix-tures by NMR
John Wiley & Son, Chichester [11] Kowalewski, J.; M¨aler, L.
Nuclear Spin Relaxation inLiquids: Theory, Experiments, and Applications
Taylor& Francis Group [12] Belorizky, E.; Fries, P. H.; Helm, L.; Kowalewski, J.;Kruk, D.; Sharp, R. R.; Westlund, P.-O. Comparison ofdifferent methods for calculating the paramagnetic relax-ation enhancement of nuclear spins as a function of themagnetic field
J. Chem. Phys. ions J. Chem. Phys.
63, 2279[14] Southwood-Jones, R. V.; Earl, W. L.; Newman, K. E.;Merbach, A. E. Oxygen-17 NMR and EPR studies ofwater exchange from the first coordination sphere ofgadolinium(III) aquoion and gadolinium(III) propylene-diaminetetraacetate
J. Chem. Phys.
73, 5909[15] Banci, L.; Bertini, I.; Luchinat, C. H NMR studies ofsolutions of paramagentic metal ions in etheleneglycol
Inorgan. Chimica Acta
Chem. Rev.
87, 901–927[17] Powell, D. H.; Merbach, A. E.; Gonz´alez, G.; Br¨ucher,E.; Micskei, K.; Ottaviani, M. F.; K¨ohler, K.; Zelewsky,A. V.; Grinberg, O. Y.; Lebedev, Y. S. Magnetic-Field-Dependent Electronic Relaxation of Gd in AqueousSolutions of the Complexes [Gd(H O) ] , [Gd(propane-1,3-diamine-N,N,N’,N’-tetraacetate)(H O) ] − , and [Gd(N,N’-bis[(N-methylcarbamoyl)methyl]-3-azapentane-1,5-diamine-3,N,N’-triaetate)(H O)] ofinterest in magnetic-resonance imaging
Helv. Chim.Acta
76, 2129–2146[18] Micskei, K.; Powell, D. H.; Helm, L.; Br¨ucher, E.;Merbach, A. E. Water Exchange on [Gd(H ] and[Gd(PDTA)(H O) ] − in Aqueous Solution: a Variable-Pressure, -Temperature and -Magnetic Field O NMRStudy
Magn. Res. Chem.
31, 1011–1020[19] Strandberg, E.; Westlund, P.-O. H NMRD profileand ESR lineshape calculation for an isotropic elec-tron spin system with S = 7/2 a generalized modifiedSolomon–Bloembergen–Morgan theory for nonextreme-narrowing conditions
J. Magn. Reson.
25, 261–266[20] Powell, D. H.; Dhubhghaill, O. M. N.; Pubanz, D.; Helm,L.; Lebedev, Y. S.; Schlaepfer, W.; Merbach, A. E. Struc-tural and dynamic parameters obtained from O NMR,EPR, and NMRD studies of monomeric and dimericGd complexes of interest in magnetic resonance imag-ing: An integrated and theoretically self-consistent ap-proach J. Amer. Chem. Soc.
Chem. Rev.
J. Chem.Phys.
Solution NMRof Paramagntic Molecules; Applications to Metallo-biomolecules and Models
Elsevier [24] Borel, A.; Yerly, F.; Helm, L.; Merbach, A. E. Multiex-ponential electronic spin relaxation and redfield’s limit inGd(III) complexes in solution: Consequences for O/ HNMR and EPR simultaneous analysis
J. Amer. Chem.Soc.
124 (9), 2042–2048[25] Zhou, X.; Caravan, P.; Clarkson, R. B.; Westlund, P.-O. On the philosophy of optimizing contrast agents. Ananalysis of H NMRD profiles and ESR lineshapes of theGd(III) complex MS-325 +HSA
J. Magn. Reson.
Spectrochim. Acta A
62, 76–82[27] Yazyev, O. V.; Helm, L. Gadolinium (III) ion in liquidwater: Structure, dynamics, and magnetic interactionsfrom first principles
J. Chem. Phys.
Eur. J. Inorg. Chem. O) +39 complex Phys. Chem.Chem. Phys.
11, 10368–10376[30] Luchinat, C.; Parigi, G.; Ravera, E. Can metal ion com-plexes be used as polarizing agents for solution DNP? Atheoretical discussion
J. Biomol. NMR
58, 239–249[31] S, A.; Botta, M.; Esteban-G´omez, D.; Platas-Iglesias, C.Characterisation of magnetic resonance imaging (MRI) contrast agents using NMR relaxometry Mol. Phys.
117 (7-8), 898–909[32] Fragai, M.; Ravera, E.; Tedoldi, F.; Luchinat, C.; Pa-rigi, G. Relaxivity of Gd-based MRI contrast agentsin crosslinked hyaluronic acid as a model for tissues
ChemPhysChem
20, 2204–2209[33] Li, H.; Meade, T. J. Molecular magnetic resonance imag-ing with Gd(III)-based contrast agents: Challenges andkey advances
J. Amer. Chem. Soc.
Chem. Rev.
J. Magn. Reson.
J. Magn. Reson.
123 (A), 95–104[37] Straley, C. A mechanism for the temperature dependenceof the surface relaxation rate in carbonates
Soc. CoreAnalysts
SCA2002-27[38] Zhang, G. Q.; Hirasaki, G. J.; House, W. V. Internalfield gradients in porous media
Petrophysics
44 (6),422–434[39] Mitchell, J.; Gladden, L. F.; Chandrasekera, T. C.; Ford-ham, E. J. Low-field permanent magnets for industrialprocess and quality control
Prog. Nucl. Magn. Reson.Spect.
76, 1–60[40] Faux, D. A.; Cachia, S.-H. P.; McDonald, P. J.; Bhatt,J. S.; Howlett, N. C.; Churakov, S. V. Model for theinterpretation of nuclear magnetic resonance relaxometryof hydrated porous silicate materials
Phys. Rev. E
91, 032311[41] Saidian, M.; Prasad, M. Effect of mineralogy on nuclearmagnetic resonance surface relaxivity: A case study ofMiddle Bakken and Three Forks formations
Fuel
Phys. Rev.
73 (7), 679–712[43] Torrey, H. C. Nuclear spin relaxation by translationaldiffusion
Phys. Rev.
92 (4), 962–969[44] Hwang, L.-P.; Freed, J. H. Dynamic effects of pair corre-lation functions on spin relaxation by translational diffu-sion in liquids
J. Chem. Phys.
63 (9), 4017–4025[45] Lipari, G.; Szabo, A. Model-free approach to the in-terpretation of nuclear magnetic resonance relaxation inmacromolecules. 1. Theory and range of validity
J. Amer.Chem. Soc.
J.Amer. Chem. Soc.
J. Magn.Reson.
SEG/AAPG/EAGE/SPE Researchand Development Petroleum Conference and Exhibition
J. Chem. Phys.
148 (16),164507[50] Asthagiri, D.; Chapman, W. G.; Hirasaki, G. J.; Singer,P. M. NMR H– H dipole relaxation in fluids: Relaxationof individual H– H pairs versus relaxation of molecularmodes
J. Phys. Chem. B
124 (47), 10802–10810[51] Singer, P. M.; Asthagiri, D.; Chapman, W. G.; Hi-rasaki, G. J. NMR spin-rotation relaxation and diffusionof methane
J. Chem. Phys.
148 (20), 204504[52] Singer, P. M.; Valiya Parambathu, A.; Wang, X.; Astha-giri, D.; Chapman, W. G.; Hirasaki, G. J.; Fleury, M.Elucidating the H NMR relaxation mechanism in poly-disperse polymers and bitumen using measurements, MDsimulations, and models
J. Phys. Chem. B n -heptane in a polymer matrix revealed by MD simulations J. Phys. Chem. B
124 (18), 3801–3810[54] Ren, P.; Ponder, J. W. Polarizable Atomic Multipole Wa-ter Model for Molecular Mechanics Simulation
J. Phys.Chem. B
J. Chem. Phys.
Journal of Chemical Theory and Computation
14, 5273–5289[57] Asthagiri, D.; Pratt, L. R.; Paulaitis, M. E.; Rempe,S. B. Hydration structure and free energy of biomolecu-larly specific aqueous dications, including Zn and firsttransition row metals J. Amer. Chem. Soc.
J. Chem. Phys.
97, 1990–2001[59] Lagardere, L.; Aviat, F.; Piquemal, J.-P. Pushing thelimits of multiple-time-step strategies for polarizablepoint dipole molecular dynamics
J. Phys. Chem. Lett.
10, 2593–2599[60] Nos´e, S. A molecular dynamics method for simulation inthe canonical ensemble
Mol. Phys.
52, 255–268[61] Hoover, W. G. Canonical dynamics: Equilibrium phase-space distribution
Phys. Rev. A
31, 1695–1697[62] Marcus, Y.
Ion solvation
John Wiley and Sons, Oxford,UK [63] Andersen, H. C. Molecular dynamics at constant pressureand/or temperature
J. Chem. Phys.
72, 2384–2393[64] Impey, R. W.; Madden, P. A.; McDonald, I. R. Hydrationand mobility of ions in solution
J. Phys. Chem.
Coord. Chem. Rev.
J. Phys. Chem. B [67] Garc´ıa, A. E.; Stiller, L. Computation of the mean resi-dence time of water in the hydration shell of biomolecules J. Comput. Chem.
14, 1396–1406[68] Merchant, S.; Asthagiri, D. Thermodynamically domi-nant hydration structures of aqueous ions
J. Chem. Phys.
Chem. Phys. Lett.
Coord. Chem. Rev.
Principles of Nuclear Magnetism
OxfordUniversity Press, International Series of Monographs onPhysics [72] McConnell, J.
The Theory of Nuclear Magnetic Relax-ation in Liquids
Cambridge University Press [73] Cowan, B.
Nuclear Magnetic Resonance and Relaxation
Cambridge University Press [74] Kimmich, R.
NMR Tomography, Diffusometry and Re-laxometry
Springer-Verlag [75] Peter, C.; Daura, X.; van Gunsteren, W. F. Calculationof NMR-relaxation parameters for flexible molecules frommolecular dynamics simulations
J. Biomol. NMR
20, 297–310[76] Case, D. A. Molecular dynamics and NMR spin relax-ation in proteins
Acc. Chem. Res.
35 (6), 325–331 [77] Venkataramanan, L.; Song, Y.-Q.; H¨urlimann, M. D.Solving fredholm integrals of the first kind with tensorproduct structure in 2 and 2.5 dimensions
IEEE Trans.Sig. Process.
50 (5), 1017–1026[78] Song, Y.-Q.; Venkataramanan, L.; H¨urlimann, M. D.;Flaum, M.; Frulla, P.; Straley, C. T - T correlation spec-tra obtained using fast two-dimensional laplace inversion J. Magn. Reson.
Diffusion Fundam.
29, 1–8[80] Singer, P. M.; Arsenault, A.; Imai, T.; Fujita, M.
La NMR investigation of the interplay between lat-tice, charge, and spin dynamics in the charge-orderedhigh- T c cuprate La . Ba . CuO Physical Review B
J. Phys.Chem. B
J. Chem. Phys. D - T relaxation-diffusion correlation of n -alkanes Appl. Magn. Reson.
47 (12), 1391–1408[84] Madhavi, W. A. M.; Weerasinghe, S.; Momot, K. I.Rotational-diffusion propagator of the intramolecularproton- proton vector in liquid water: A molecular dy-namics study
J. Phys. Chem. B2017