Field-dependent paramagnetic relaxation enhancement in solutions of Ni(II): What happens above the NMR proton frequency of 1 GHz?

An extended set of paramagnetic relaxation enhancement (PRE) data, up to the field of 32.9 Tesla, is reported for protons in an acidified aqueous solution of a Ni(II) salt in the presence and in the absence of added glycerol. For the 55% w/w glycerol sample, a distinct maximum in the PRE vs magnetic field curve is observed for the first time. The data are analysed using the Swedish slow-motion theory, including both the intramolecular (inner-sphere) and intermolecular (outer-sphere) contributions. The results indicate that estimating the outer-sphere part in the presence of the more efficient inner-sphere term is a difficult task.


Introduction
Studies of nuclear spin relaxation in paramagnetic systems in solution are an important field of research, described in recent monographs [1,2] and reviews [3,4]. The relaxation measurements are often carried out as a function of the magnetic field, in experiments known as nuclear magnetic relaxation dispersion (NMRD). About a decade ago, we reported an investigation of proton NMRD for aqueous solutions of acidified nickel(II) trifluoromethanesulfonate (an anion with low tendency to form ion pairs [5]), with or without added glycerol [6]. That study will be denoted as I in the following. The present work extends that study in two ways. First, we expand the range of magnetic fields up to about 33 Tesla, using the resistive magnet at Laboratoire National des Champs Magnétiques Intenses (LNCMI) in Grenoble [7]. In this way the very broad range of magnetic fields from the limit of zero-field splitting (ZFS) of Ni(II) dominating the proton Zeeman interaction to the opposite limit when ZFS is dominated by the Zeeman coupling has been covered.
Second, we interpret the data making use of a protocol capable to fit a model including both the inner sphere [8,9] and the outer sphere [10] paramagnetic relaxation enhancement (PRE). The inner-sphere part of the code was some time ago compared to a mathematically very different (but physically equivalent) approach proposed by the group in Grenoble. It was found that the two computational approaches agreed very well with each other [11].
Spin relaxation processes in low pH aqueous solutions of Ni(II) salts have been studied many times, starting with the classical work reported by Morgan and Nolle [12]. They found no field dependence for the proton spin lattice relaxation time in the field range 0.05-1.4 Tesla. Bloembergen and Morgan [13] explained this observation introducing a theory of rapid electron spin relaxation, caused by fluctuations of the zero-field splitting, originating from distortions of the hydration sphere of the ion by collisions. In late eighties, we reported a study of low-pH aqueous solution of Ni(II) at a series of magnetic field up to 11.7 Tesla [14] and interpreted the experimental data using a theoretical model of Westlund et al. [15], again assuming modulation of the ZFS by distortions of the coordination sphere. The same data sets were re-interpreted by Svoboda et al. [16], using a slightly more sophisticated variety of the theory. The experiments show that the PRE increases with increasing magnetic field in the high field range. This is in agreement with the theory, which predicts that the NMRD curve should reach a maximum and then turn down. In I, we tried to reach the maximum by extending the measurements to even higher magnetic field, 21.1 Tesla (900 MHz 1 H resonance frequency). In addition, the NMRD profiles were measured for a Ni(II) salt dissolved in a mixture of water and glycerol, deuterated in the chain. However, the expected decrease of the relaxation rate at high field was still not reached at that time. Here, this PRE reduction at high field was indeed observed experimentally.

Materials and methods
The nickel(II) trifluoromethanesulfonate solutions without glycerol and with 55% w/w glycerol were the same as used in I [6]. Similarly to I, the measurements up to 1.066 Tesla were performed on 10 mm samples using a Stellar relaxometer. The range 2.1-18.8 Tesla was covered by commercial high resolution spectrometers. Precautions were taken to avoid radiation damping problems. At some fields, the saturation-recovery and inversion-recovery obtained at different laboratories were compared and found to give good agreement with each other. The high field (above 20 Tesla) measurements were performed at LNCMI using inversionrecovery method on 1 mm capillary in a 24 MW resistive magnet. The high-field measurements followed the protocol for T 1 experiments on a different system published few years ago [17]. All measurements were performed at 295 K.
The paramagnetic relaxation enhancements were obtained by subtracting the field-independent water proton relaxation rates from the measured relaxation rates for the paramagnetic solutions. With the relatively high nickel(II) concentrations, the diamagnetic corrections were always less than 1% of the paramagnetic rates. The PREs were converted into relaxivities (PRE at unit concentration of paramagnetic ions) using the sample composition data in Table 1 of I. Using relaxivities in the fitting procedures is advantageous when both the inner-sphere and the outer-sphere contributions are included in the calculations.
The relaxivity, r, at a given field is given as a sum of the innerand outer-sphere contributions [10]: c Ni is the Ni(II) ion concentration (in mMol units), R 1;I is the incomplex spin-lattice relaxation rate for nuclear spin I (here 1 H), p is the ratio of molar concentrations of paramagnetic species and the free ligands, q is the number of ligand molecules in the coordination shell. R 1;outer is the outer-sphere relaxation rate. R 1;I and R 1;outer are the quantities obtained by theoretical calculations, here making use of the so-called Swedish slow motion theory [3,[8][9][10]. Very briefly, the method is based on a mixed quantum mechanical and classical description of the lattice, which is handled in the spirit of the stochastic Liouville equation. The approach retains its validity even beyond the perturbation limit (the Redfield limit). The relaxation rate R 1;I is expressed in terms of a composite spectral density function taken at the nuclear spin (here proton) Larmor frequency. Calculation of this spectral density involves setting up and inversion of a very large matrix, representing the Liouville superoperator in a basis set containing direct products of Wigner rotation matrices and appropriate spin-space operators. The relaxation rate depends on the model parameters: interspin distance in the complex (r IS ), rotational and distortional correlation times (s R ; s v , respectively), static and transient ZFS parameters ðD S ; D T Þ.
The static ZFS refers to the part of the ZFS which is independent of time in the molecular frame, while the transient ZFS characterizes the fluctuations caused by distortions of the hydration sphere. Computation of the outer-sphere relaxation rate is based on taking numerical integral of a product of a function related to the innersphere spectral density and a translational diffusion correlation function. R 1;outer depends on the same correlation times and ZFS parameters as above and, in addition, on the distance of closest approach d between the spins in the intermolecular case and the mutual translational diffusion coefficient D tr .
Fitting the model parameters to the experimental data was performed iteratively using the slow-motion Fortran program interfaced to the Minuit program package [18]. The fitting strategy for the glycerol-containing solution involved least-squares optimization of five inner-sphere parameters ðr IS ; s R ; s v ; D S ; D T Þ for given pairs ðd; D tr Þ. The minimum in the seven-parameter space was localized comparing the sum of square deviations for different ðd; D tr Þ pairs. For the solution in pure water, we assumed that the static ZFS vanished, in analogy with I. The outer-sphere calculations are computationally demanding and we have developed a version of the program running on 32 cores of a parallel computer.
Computations have been carried out on a local computational cluster, consisting of 4 Dell PowerEdge 815 nodes. A single fiveparameters fitting run took few days of the computer time, while less than a day was required for the D S ¼ 0 case.

Results and discussion
The measured field-dependent relaxivities are shown in Fig. 1. For convenience of the reader, the relaxivities are also listed in Table A1 of the appendix. We notice that the case with 55% glycerol shows quite clearly the relaxivity maximum at high field, while for the aqueous sample there is only a weak indication of the maximum.
In I, we have essentially used the fitting to an inner-sphere model, complemented by estimates of the outer-sphere contributions (which were subsequently found to be underestimated due to erroneous scaling). Here, we use a fitting approach including both mechanisms. We initiate the analysis of the data for both samples by an inner-only fitting, which provides a sum-of-square deviations to be used as reference point in the following inner + outer modelling.

Sample without glycerol
The inner-only fit leads to the parameter set shown in the first line of Table 1. It may be of interest to compare the parameters with those reported in I as global fit at 323 K. The interspin distance in the present work is somewhat shorter than 240 pm given in I. The rotational correlation time is longer than the earlier value (s R ¼ 6:8 ps ), which is expected because of the temperature difference. On the other hand, the distortional correlation time is close to the value from I (s v ¼ 2:3 ps), which is in line with suggestion presented there that, because of the crudeness of the model, this parameter does not need to be assumed temperature dependent. Finally, the transient ZFS in Table 1 is somewhat larger than the value D T ¼ 5:1 cm À1 in the earlier study.
The results of the present work based on the inner + outer sphere model (the second line in Table 1 Table 1 (inner + outer A) correspond to the lowest point in the v 2 -grid. The corresponding v 2 value goes down by about 10% with respect to the inner-only reference. The minimum found in the surface is very shallow and poorly defined. We can notice that the inner-sphere parameters are affected to a small extent by inclusion of the outer-sphere contribution. This goes well together with conclusion of the previous study concerning the dominant role of the innersphere mechanism.
In addition to this procedure, we also report the result of fiveparameter fit, where both the ZFS parameters are kept at fixed values, while D tr and d are optimized along with r IS and the two correlation times (inner + outer B). This leads to a very slight improvement of the v 2 -value and a rather large change in the parameters characterizing the outer sphere contribution, in particular D tr . The latter quantity comes out too small in both fits. We have therefore also performed fitting with a fixed value of D tr = 2.3Á10 À9 m 2 s À1 , the diffusion coefficient for pure water [19]. This is shown in Table 1 as ''inner + outer C". The change in parameters is dramatic, the r IS becomes much longer (and the innersphere relaxation rate is proportional to r IS À6 ) while d goes down strongly, which makes the outer sphere contribution dominant over the inner sphere. In our opinion, this is not physically reasonable and we rather believe in the lines (inner + outer A) or (inner + outer B) and the interpretation of D tr as some kind of ''effective" parameter.
The calculated curve based on the ''inner + outer, A" scheme is shown in Fig. 2, together with the experimental points.

Sample with 55% glycerol
The results of the fitting for the glycerol-containing sample are summarized in Table 2. Here, we assume that the glycerol molecule (deuterated in the chain) can enter the first coordination shell, distorting the average symmetry of the complex. We assume that the number of protons in the first shell is 11. As a consequence, the static ZFS is no longer zero. The inner-only fit parameters are, again, in reasonable agreement with I. The most significant differences are in the r IS and s v . Including the outer-sphere mechanism results in about 7% reduction of the target function (v 2 ) and in large changes of the best-fit parameters. The variation of the v 2 as a function of the outer-sphere parameters is shown in Fig. A2 in the appendix. We notice that the r IS increases very substantially, which implies a reduced importance of the inner-shell component. At the same time, the distance of closest approach in the outer-sphere contribution comes out very short, increasing the importance to the outer-sphere part, in analogy with the last line of Table 1. Mathematically, this opposite behaviour makes sense but the physical picture is unclear. The third line in Table 2 (inner + outer B) corresponds to the five-parameter fitting based on fixed values of the ZFS parameters and simultaneous fit of r IS , the two correlation times, d and D tr . The results are in fair agreement with the gridtype calculations. The change of the diffusion coefficient with respect to the aqueous solution (fits ''inner + outer A" and ''inner + outer B" in Table 1) is contrary to expectations, but the value as such is in reasonable agreement with about 4Á10 À10 m 2 s À1 estimated through interpolation from experimental data [20,21]. Thus, we did not see any reason to perform fits with fixed value of D tr . A possible conclusion to draw is that the fitting process may contain too many parameters and that the force-free diffusion approach gives an oversimplified description of the real physics.
As a possible tool to find a physically more satisfying description of the PRE profile, we have decided to try introducing a constraint on some parameters that can be estimated by independent means. First, we performed molecular dynamics (MD) simulations for a system of similar composition and used the trajectory to compute the time-correlation function (tcf) for the rank-2 spherical harmonics describing the orientation of the nickel-to-aqueous proton (or the proton of coordinated OH-group of glycerol) axis with respect to the z-axis of the laboratory frame. The details of the MD simulation are provided in the appendix and the tcf's for different protons in the first coordination sphere are shown in Fig. 3. These functions are shown in red colour. The black line is the averaged autocorrelation function approximated by a single exponential (green). The timescale is in picoseconds. The  Theoretical and experimental relaxivities r (in s À1 mMol À1 ) for the sample without glycerol as a function of the magnetic field for the best-fit parameters (inner + outer A) in Table 1. corresponding correlation time is about 485 ps. With the assumption of this value of s R , we have performed a six-parameter fit shown in Table 2 as ''Inner + outer C". Again, the fitting leads to a small inner-sphere contribution, not substantially different from the ''inner + outer B" case.
Another parameter which can be estimated independently is the static ZFS. This was achieved by means of ab initio calculations (n-electron valence state perturbation theory with a minimal active space of orbitals NEVPT2(8,5)) similar to the earlier work from one of our laboratories [22] at a series the geometric configurations of a Ni(II)(H 2 O) 5 (glycerol) obtained from a MD trajectory mentioned above. The details of the calculations are given in the appendix. Here, we present some important features of the results. First, we show in Fig. 4 the ZFS energy levels calculated along the trajectory. We chose to present the results in this form, rather than in terms of the usual D and E parameters, to avoid the complication of frequent sign changes in D. Clearly, the level structure is such that the most negative and the most positive level take on the role of the D zz = D parameter and the middle level is relatively close to zero. We can use the splittings along the trajectory to estimate the magnitude of the static D S j j ¼ 5:06 cm À1 .
The results of the fitting using the fixed values of s R and D S are shown in Table 2 as ''inner + outer D". Again, similarly to the ''Inner + outer C", the essential features of the fitted parameters (large outer sphere contribution and small inner sphere part) remain the same. The ZFS data along the trajectory (Fig. 4) suggest a possible modification of the model used for fitting: rather than the cylindrically-symmetric ZFS tensor, we might assume the three levels to be equidistant, so that E = D/3 (the fully rhombic case). This way of thinking was applied long time ago for the case of Ni The calculated curve based on the ''inner + outer, A" scheme for the 55% glycerol sample is shown in Fig. 5 together with the experimental points. a Assumed and kept fixed. b Upper bound of the allowed range.   Theoretical and experimental relaxivities r (in s À1 mMol À1 ) for the 55% glycerol as a function of the magnetic field for the best-fit parameters (inner + outer A) in Table 2.
(II) in aqueous solution [16]. An important property of this model is that it keeps the number of fitting parameters unchanged. The results are displayed in Table 2 as ''inner + outer E". Once more, the target function goes up with respect to the ''inner + outer A" or ''B", while the problem of the dominant outer sphere contribution remains. As a last attempt to resolve this problem, we make a rather dramatic assumption D S = 0. The result of these fits, with and without locking the rotational correlation time, are shown in Table 2 as ''inner + outer F" and ''inner + outer G". These fits result in a large value of the translational diffusion coefficient, long intermolecular distance of closest approach and a high value of the target function.

Concluding remarks
In this study, we were finally able to observe a high-field maximum in the NMRD profile for a solution of Ni(II) salt, characterized by fast electron spin relaxation (the maxima are well-known for paramagnetic ions with more symmetric electronic structure and with slower electron relaxation, such as Gd(III) or Mn(II) [1]). The maximum is well-defined for the solution containing 55% w/w of glycerol, while it is only weakly indicated for the aqueous sample. The data were interpreted using a model allowing for both the intramolecular (inner-shell) and intermolecular (outer shell) contributions, considered within the Swedish slow-motion approach [8][9][10]. The static and transient ZFS tensors are (in most cases) assumed cylindrically symmetric with the principal axis coinciding with the principal axis of the intramolecular dipole-dipole interaction. For the 55% glycerol sample, we also tried the possibility of locking certain parameters and of the fully rhombic static ZFS tensor, i.e. E S = D S /3. The fitted parameters for the 55% glycerol sample indicate, for cases A-E, a significantly reduced intramolecular component, compensated by a probably overestimated outer-sphere contribution. For the last two cases in Table 2, where D S ¼ 0 is assumed, the intramolecular part is more up to expectations, but the fit quality (measured by the target function) is rather poor.
For the sample without glycerol, the minimum in the v 2 surface is very shallow and the outer-sphere component appears to be a small correction, except for the case of fixed value of the diffusion coefficient (inner + outer C). The results indicate that estimating the outer-sphere part in the presence of the more efficient innersphere term is a difficult task. One should also be aware that the concept of outer-sphere relaxation has been introduced to take into account the relaxation contribution originating from solvent molecules not belonging to the coordination sphere, but the assumption that these molecules perform free translation dynamics is an oversimplification.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Relaxation data and details of molecular dynamics simulations
Relaxation data Relaxivities as a function of magnetic field are collected in Table A1.
Molecular dynamics simulation details MD simulations were performed using the MDynaMix [23] package. A dilute solution of Ni(II) ion in 72 wt% glycerol-water mixture (1 ion, 170 glycerol and 340 water molecules in the simulation cell) was studied. Counter-ions were not included because they are expected to disturb the nickel complex structure and dynamics. All simulations were carried out in an isothermal-  isobaric (NPT) ensemble at 298 K and atmospheric pressure. The temperature was kept constant by the Nose-Hoover method [24,25], and the pressure was regulated by the Hoover barostat [25]. The equations of motion were solved using the double time step algorithm by Tuckerman et al. [26,27] with a long time step of 10 short time steps of 0.2 fs. The long-range Coulomb forces were treated by the Ewald summation method. The systems considered are not electro-neutral but the neutrality requirement in the Ewald algorithm was satisfied by adding a uniform background charge density that cancels the ion charge [28]. The corresponding corrections are implemented in the MDynaMix package.
Interaction potentials Any classical simulation methodology is always as good as the used potential models describing the interactions between the molecules. In the present study water was simulated with the flexible SPC-F model by Toukan and Rahman [29]. To describe glycerol interactions the slightly modified version of Blieck et al. [30] potential was used. The full description of model parameters employed is given in Ref. [31]. To describe the nickel ion-water interaction an effective 2-body potential proposed by Chillemi et al. [32] was used. To our knowledge, there is no potential especially adopted to describe the Ni(II)-glycerol interactions available in the literature. In the present simulations the standard sum of Coulomb and Lennard-Jones 6-12 potentials was used. Cation Lennard-Jones parameters r Ni = 2.525 Å, e Ni = 0.0628 kJ/mol proposed by Rappe et al. [33] were employed.
Nickel cation solvation The model system was equilibrated during the 3 ns simulation. The following production run was 7 ns. It is necessary to mention here, that nickel complexes both with water and glycerol are characterized by slow ligand exchange (corresponding lifetimes are 4 Ä 5Á10 À5 s for water [34] and about 2.5Á10 À4 s for glycerol [35]). In this case the standard MD simulations even on a few microseconds time scale could not achieve complete Boltzmann sampling. It was found that in 72 wt% glycerol-water mixture the first coordination sphere of the nickel ion consisted of five waters and one glycerol and remained stable over the whole simulation time of 7 ns. The glycerol molecule acts as a monodentate ligand and is coordinated with respect to the Ni(II) ion by central OHgroup.
Computational details of the ZFS calculations From the MD simulations of aqueous Ni(II) salt glycerol solutions (72 wt%) we sampled 889 configurations each separated by 200 fs. The cubic simulation cell with a cell side of 31.40 Å contained 1 nickel, 340 H 2 O and 170 glycerols, from which we extracted clusters of the local nickel coordination, which was Ni (II) (H 2 O) 5 (glycerol) during the entire trajectory. Furthermore, the glycerol molecule was coordinated with respect to the Ni(II) by central OH-group.
The zero-field splitting of the ground state triplet was obtained using multi-configuration quantum chemical calculations in ORCA 4.1.2 [36] within the Douglas-Kroll-Hess (DKH) formalism [37,38] in which SOC effects are added to scalar relativistic calculations of states with pure electron spin multiplicities. Multiconfigurational scalar relativistic calculations were performed using the DKH-def2-TZVPP basis set from ORCA with exponents from def2-TZVPP [39], re-contracted for 2nd order DKH. Using a minimal active space containing the nickel 3d orbitals, quasidegenerate perturbation theory based on NEVPT2(8,5) calculations of 10 triplets and 15 singlets with spin-orbit coupling was used for calculation of the zero-field splitting of the triplet ground state in Ni(II), as in our previous studies of the aqueous Ni(II) complex [22]. For efficiency, we employed the resolution of identity and chain of spheres (RIJCOSX) approximation using resolution of identity for the coulomb integrals and COSX numerical integration of exchange integrals. The AutoAux generation procedure [40] was used to generate the auxiliary basis sets. Evaluation of the RIJCOSX approximation and of the basis set dependence by increase to the QZVPP basis set or the ANO-RCC-VTZP basis set showed stability within 0.07 cm À1 for the zero-field splitting.