Structure and Spectromagnetic Properties of the Superoxide Radical

Aug 19, 2010 - (5, 6) This technique involves the addition of a short-lived free radical to a ... of DMPO in biological spin trapping are limited by t...
0 downloads 0 Views 3MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 11793–11803

11793

Structure and Spectromagnetic Properties of the Superoxide Radical Adduct of DMPO in Water: Elucidation by Theoretical Investigations Ce´line Houriez,*,† Nicolas Ferre´,† Didier Siri,† Paul Tordo,† and Michel Masella‡ UMR 6264 Laboratoire Chimie ProVence, Faculte´ des Sciences de Saint-Je´rome Case 521, AVenue Escadrille Normandie-Niemen, 13397 Marseille Cedex 20, France, and Laboratoire de Chimie du ViVant, SerVice d’inge´nierie mole´culaire des prote´ines, Institut de biologie et de technologies de Saclay, Commissariat a l’e´nergie atomique, Centre de Saclay, 91191 Gif-sur-YVette Cedex, France ReceiVed: April 14, 2010; ReVised Manuscript ReceiVed: July 26, 2010

In the field of spin trapping chemistry, the design of more efficient radical traps can be assisted by the development of theoretical methods able to give a quantitative evaluation of the electron paramagnetic resonance (EPR) spectrum features of the spin-adduct radical, even before initiating the experimental work. The superoxide radical adduct of the 5,5-dimethyl-1-pyrroline-N-oxide nitrone (DMPO-OOH) has been reported in a huge number of papers devoted to the study of the oxidative stress. Here, we present for the first time the theoretical study of DMPO-OOH in an explicit water solution, based on the combined QM/MM//MD protocol we recently proposed, featuring a full coupling between the solute and all the explicit water molecules. Our results show that the DMPO-OOH EPR spectrum, whose interpretation is still debated, can be explained in the light of two sites in chemical exchange, in agreement with the most recent experimental data. Moreover, we demonstrate that each site consists of an equilibrium between the two main 5-membered ring conformations of DMPO-OOH. We provide also an analysis of the solvent contribution to the hyperfine coupling constants (hcc’s) as well as an exhaustive study of the possible relationship between the hcc’s and the main structural characteristics of DMPO-OOH. Our QM/MM//MD protocol appears thus to be an accurate theoretical tool allowing the investigation of the magnetic properties of large nitroxide spin adducts in complex environments. I. Introduction The superoxide radical O2-• has been evidenced as one of the major reactive oxygen species playing a key role in numerous physiological processes and human diseases.1-3 Hence, the ability of detecting in real time the presence of O2-• in a biological milieu is crucial, and different techniques of identification have deserved much attention during the last years.4 Among them, the spin trapping technique coupled with electron paramagnetic resonance (EPR) spectroscopy is intensively used to detect this radical both in vitro and in vivo.5,6 This technique involves the addition of a short-lived free radical to a diamagnetic spin trap, leading to the formation of a spin-adduct radical that is persistent enough to be detected by EPR spectroscopy. Beyond the access to the specific signature of this new paramagnetic molecule, the spectrum analysis provides direct information concerning the spin adduct conformational behavior.7-9 Five-membered cyclic nitrones based on the pyrroline-Noxide moiety have been extensively used during the last 15 years.10-15 One of the most popular and widely used traps is the 5,5-dimethyl-1-pyrroline-N-oxide (DMPO), a 5-membered cyclic nitrone whose reaction with the superoxide radical leads to the formation of the DMPO-OOH nitroxide spin adduct (see Figure 1). The distinctive features of the DMPO-OOH EPR spectrum, especially as compared to the spectrum of DMPO-OH (resulting from trapping the free radical OH•), allow one to detect the presence of O2-•, and that has had a profound influence on the understanding of the role played by O2-• and OH•, respec* To whom correspondence should be addressed. Present address: CEA, DAM, DIF, 91297 Arpajon, France. † UMR 6264 Laboratoire Chimie Provence. ‡ Institut de biologie et de technologies de Saclay.

Figure 1. Spin trapping of the superoxide radical by the DMPO nitrone.

tively, in numerous biological processes involving an oxidative stress. However, applications of DMPO in biological spin trapping are limited by the short half-life time of the DMPO-OOH spin adduct (t1/2 = 60 s at pH 7.2) and by its rapid decay to the DMPO-OH spin adduct. Thus, a large number of pyrroline-N-oxides overcoming to some extent these limitations have been designed in the last 15 years.13 Nevertheless, their different behavior compared to DMPO has not been satisfactorily explained. The DMPO-OOH EPR spectrum, reported for the first time 35 years ago,16 is marked by a noticeable asymmetry whose origin has been discussed for years.17-21 To interpret this spectrum, two main alternative assumptions were proposed: (A) the superimposition of two 6-line EPR spectra characterized by different nitrogen and β-hydrogen hyperfine coupling constants (hcc’s) and no γ-hydrogen hcc18 and (B) the 2-site chemical exchange broadening of a 12-line spectrum involving hyperfine couplings with nitrogen and β- and γ-hydrogen nuclei.22 Different sets of parameters deduced from the best fit between the calculated and experimental spectra recorded in different experimental conditions have been proposed and are collected in Table 1. As already mentioned, model A is based on the superimposition hypothesis, while parameter sets B1,19 B2,22

10.1021/jp1033307  2010 American Chemical Society Published on Web 08/19/2010

11794

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Houriez et al.

TABLE 1: Four Sets of hcc Parameters aX (X ) N, Hβ, Hγ) Allowing One to Reproduce Satisfactorily the DMPO-OOH Experimental EPR Spectrum in Water Solution Based on Standard Spectrum Simulation Techniques and by Considering Two Sites c1 and c2 as Major Contributors to the Spectrum (the hcc’s of c2 Are Provided in Parentheses)a model

aN

aH β

aH γ

A B1 B2 B3

13.8 (13.8) 14.2 (14.1) 14.4 (14.2) 14.1 (14.1)

12.2 (9.8) 11.3 (11.8) 12.8 (10.6) 11.8 (10.9)

1.6 (0.0) 1.1 (1.4) 0.9 (1.4)

a Details concerning the models A,18 B1,19 B2,22 and B320 are provided in the text. All values are in Gauss.

and B320 have been obtained from the chemical exchange model. In 2005, the combination of quantum chemical calculations and experimental investigation of γ-deuterated DMPO-OOH species has evidenced a coupling with γ-hydrogen nuclei, validating the chemical exchange interpretation.19 However, up to now, the nature of the two sites (denoted by c1 and c2 hereafter) remains unknown. Their determination should be essential to fully understand the EPR features of DMPO-OOH and the origin of its short lifetime. At the exception of a few simple molecules for which the microwave technique is applicable,23 usually no experimental structural data is available for nitroxide spin adducts in solution, as in the case of DMPO-OOH. Hence, one may usefully consider molecular-based simulations for the elucidation of the conformational dependence of complex EPR spectra. Obviously, such simulation protocols have to be accurate enough to describe the solvation of the nitroxide, as well as to compute reliable hcc values, in particular on a statistical basis. Aside from the challenge of developing new spin traps more efficient than DMPO (which is characterized by a slow reactivity with O2-• and whose spin adduct exhibits an overall fast adduct decay rate), the second major challenge in the spin trapping field is to develop a theoretical method able to evaluate the EPR features of a new spin adduct, even before initiating the experimental work. With this aim, Villamena and co-workers investigated the ability of common quantum mechanical tools, mainly based on the density functional theory, to investigate the formation, the decay, and the spectromagnetic properties of some nitroxide spin adducts like DMPO-OOH. However, the molecular models they used ignored the solvent effects or accounted for them using the available continuum solvation methods, such as the polarizable continuum model (PCM).24 Moreover, the major part of their studies focused on a single nitroxide conformation, and, as in the case of DMPO-OOH, 25 the nitroxide optimized geometries considered for interpretation often exhibit intramolecular hydrogen bonds (HBs), whose existence in a protic medium such as water is questionable. Lastly, if they attempted recently to characterize the incidence of some conformational effects (involving hydroperoxide torsional angles) on hcc’s, they still ignored the possible influence of 5-membered ring nitroxide conformations on the hcc’s, which has been shown to be pivotal for understanding the EPR spectra of diphosporyl nitroxides, 26 for instance. To overcome the latter drawbacks, methods based on the investigation of nitroxide conformational behavior in explicit solvent have been proposed. They originally relied on Monte Carlo simulations and quantum computations.27,28 More recently, Pavone and co-workers29-32 proposed an integrated scheme combining a Car-Parrinello (CP) molecular dynamics (MD)

investigation of nitroxide in polar solvent, followed by the extraction of selected molecular clusters (the nitroxide and the two closest solvent molecules) whose instantaneous nitrogen hcc’s are computed at a DFT/PCM level of theory and then averaged. In the same line as Pavone and co-workers, Neugebauer and co-workers proposed a DFT//CP MD method in which the solvent molecules are treated as frozen densities able to polarize the nitroxide distribution of charges.33 However, as we recently highlighted in the case of some small nitroxides in water, the sampling of picosecond-scale trajectories usually leads to results which can be affected by large statistical biases. Accordingly, such results have to be considered with particular care.34-36 Moreover, because of the molecular size of systems such as DMPO-OOH, their study in solution based on CP MD techniques with long enough simulation times is far from being feasible. During the last years, we devised an alternative strategy based on the combination of classical MD, featuring an ad hoc parametrized polarizable force field, with large scale QM/MM calculations. We have shown34-36 that it is possible to provide statistically relevant nitrogen hcc’s with this method, reaching an agreement with experimental values as good as 0.3 G on average. In the current work, we report the application of such a scheme to the investigation of DMPO-OOH in aqueous solution. We particularly paid attention to the identification of the structural factors influencing hcc’s of nitrogen and β- and γ-hydrogen atoms of that nitroxide in water. Owing to both the strong sensitivity of the hcc values to the nitroxide structure (in particular to the 5-membered ring conformation) and the uncertainties inherent to quantum and force field computations, we defined several force field parameter sets and accordingly computed several sets of average hcc’s to check the reliability of our predictions. Moreover, we also computed the free energy profile of DMPO-OOH in explicit solvent as a function of the dihedral angle OOCN to establish the possible existence of different conformers of the DMPO-OOH in solution. Anticipating our results, this QM/MM// MD strategy allows us to confirm that the aqueous DMPO-OOH EPR spectrum can be rationalized from two sites in chemical exchange, mainly differing in the orientation of the hydroperoxyl moiety. Furthermore, we show that each site has to be seen as an average over two sets of 5-membered ring conformations, in contrast to previous theoretical results that completely ignored this degree of freedom. The paper is organized as follows. In the first part, the details of the force field we used are described, with special emphasis on key torsional degrees of freedom. In a second part, the nitroxide conformation influence on hcc’s is investigated in detail using high level quantum computations in the gas phase. Finally, the method is applied to investigate DMPO-OOH in explicit water boxes, and the corresponding structural and spectromagnetic data are analyzed and compared to available experimental data. II. Computational Details We describe here in brief our QM/MM//MD protocol devoted to compute nitroxide hcc values in aqueous solution. It involves three steps: (1) development of an accurate polarizable force field to model the intramolecular interactions within the nitroxide; (2) generation of 8-ns-long trajectories of the nitroxide embedded in a water box; (3) computation of the instantaneous hcc values of small nitroxide/water clusters extracted from the trajectories, using a suitable QM/MM level of theory.

Properties of the Superoxide Radical Adduct of DMPO

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11795

Averaging the latter set of values provides the theoretical estimate of the nitroxide hcc mean values in solution. A. Force Field Details. As in our previous studies, we consider the TCPEp polarizable force field37 based on the decomposition of the energy into five terms

U ) Urep + Uqq' + Urel + Upol + UHB

(1)

corresponding, respectively, to repulsive, classical charge-charge electrostatic, intramolecular relaxation, polarization, and specific hydrogen bond terms. The polarization term is based on an induced dipole moment approach (however, only the oxygen, carbon, and nitrogen atoms are considered as polarizable centers). The intramolecular relaxation term Urel includes classical harmonic stretching, bending, and torsional terms. With the exception of the degrees of freedom discussed below, all the parameters used in the present study are either provided as Supporting Information or available in previous studies.34,36 In the particular case of the COOH moiety, we assigned its force field intramolecular parameters to reproduce quantum B3LYP/6-311+G(d,p) results characterizing the methyl peroxide in its gas phase equilibrium geometry (characterized by a dihedral angle ∠HOOC of 133°). 1. Modeling the Hydrogen Bonds. The hydrogen bonding term UHB consists of a short-range anisotropic function

UHB )

∑ f(r)g(φ, ψ)

Figure 2. Definitions of the small nitroxide radicals MENO-OOH and IPENO-OOH. The arrow indicates the direction according to which the HβCNC (φ) dihedral angle is computed.

torsion energy term UST to the original TCPEp potential energy function ST 2 UST ) kCX[1 - cos(2ωST)](rCX - rCX )

(3)

(2)

Here, the sum runs over all the hydrogen-bonded atom pairs, and the definition of the r, φ, and ψ parameters is detailed in previous studies.34,36,37 The UHB parameters are defined to reproduce quantum results concerning only small hydrogenbonded dimers and trimers (see ref 34 in the particular case of NO · · · water HBs). In order to model the HBs taking place between the hydroperoxyl moiety and the water molecules, we altered the UHB analytical form specifically. In accordance with our own MP2 quantum computations, this was done by introducing a new angular dependence in the function g. This purpose-built function favors hydroperoxyl/water structures whose water hydrogen points toward only one hydroperoxyl oxygen, as compared to the hydroperoxyl/water bridged structure whose water hydrogen interacts with the two hydroperoxyl oxygen atoms. The details concerning such a term will be discussed in an upcoming work; however, the reliability of our improved potential may be assessed by comparing the geometrical and energetic results obtained by using this new potential with the quantum MP2/6311+G(2df,2p) ones, concerning 15 small CH3OOH/(H2O)n clusters (n ) 1-4) (see the Supporting Information). Indeed, on average, the agreement in energy is achieved within 0.5 kcal mol-1, and the root-mean-square deviation in the aggregate coordinates between both methods is about 0.2 Å. 2. Anomeric Effects Affecting the CHβ and CO Bonds. The DMPO-OOH nitroxide is affected by the so-called anomeric effects, that are usually interpreted as hyperconjugative interactions taking place mainly between the π/π* orbitals describing the 3-electron cloud of the NO bond on the one hand and the antibonding orbitals σ*(C-Hβ) and σ*(C-O) on the other hand. These effects are responsible for a noticeable increase of the C-X bond length (X ) Hβ, O) when the dihedral angle ∠XCNO is close to 90/270°. These geometrical effects are accounted for in our simulation protocol by adding a stretch-

where rCX and ωST are, respectively, the C-X bond length and the dihedral angle ∠XCNO and kCX and rST CX are two parameters. With this potential, the incidence of anomeric effect on chemical bonds is the strongest for ωST values close to 90/270°, as expected. The UST parameters were assigned together with those of Urep, UHB, and the torsion term Utor, to reproduce the potential energy surface (PES) of the DMNO-OOH radical defined by the dihedral angles ∠OCNO and ∠OOCN. The corresponding quantum and force field PESs as well as the CHβ and CO bond lengths are compared in the figures provided as Supporting Information. With our force field approach, it is possible to reproduce with a good accuracy the quantum results concerning the PES and, on average, the bond lengths. In particular, UST allows one to reproduce bond length differences of 10-2 (CHβ) and 4 × 10-2 (CO) Å for ωST values corresponding to weak and strong anomeric effects, whereas these bond lengths are predicted to be constant within less than 10-3 Å when omitting UST. 3. Key Dihedral Degrees of Freedom. Cremer and Pople demonstrated that the intrinsic flexibility of a 5-membered ring (ring puckering) can be characterized via the degrees of freedom q and ω, i.e., the amplitude and the phase of the ring deformation, respectively.38 The conformation sets of twist (T) and envelope (E) types can be found relatively to ω: they differ by the number of ring atoms lying out of the ring mean plane. In the case of the DMPO-OOH 5-membered ring, the following usual convention is adopted: a twist conformation is denoted by 3T4 if the ring carbon atoms 3 and 4, as shown in Figure 3, are, respectively, above and under the ring mean plane, and vice versa for 4T3. An envelope conformation implies that only one ring atom is out of the ring mean plane and is denoted by n E or En (with n being the atom number) according to the same convention. For the sake of readability, an equilibrium between a conformer characterized by an angle ω around 90° (i.e., E4, 3 T4, and 3E) and a conformer with ω around 270° (4E, 4T3, and E3) is referred to hereafter as a 3T4/4T3 equilibrium.

11796

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Houriez et al.

Figure 3. Schematic definitions of the 3T4/4T3 equilibria corresponding to the DMPO-OOH sites c1 and c2 discussed in the present study. The DMPO-OOH/water structures correspond to snapshots extracted along the solvated trajectories T1(T1′) and T8(T8′). The arrows indicate the directions according to which the corresponding dihedral angles φ1 and φ2 are computed (see text for a definition).

As discussed in section III, some DMPO-OOH hcc’s are highly sensitive to the latter equilibrium ratio. Moreover, we estimate the quantum uncertainties (originating from both the basis set extension and the level of theory) affecting the difference in energy between 3T4 and 4T3 conformers to be included within 0.5 and 1 kcal mol-1 (see the Supporting Information). Hence, instead of considering a single set of torsional parameters to handle the 5-membered ring motion and therefore bias our conclusions, we chose to consider eight different force field parameter sets allowing us to generate trajectories corresponding to various DMPO-OOH 3T4/4T3 equilibrium ratios. Those parameter sets only differ by the torsional parameter values related to the dihedral angles ∠OCCC and ∠HβCCC (parameters are provided as Supporting Information). With this strategy, it is possible to thoroughly investigate the relation between the DMPO-OOH 5-membered ring interconversion and the hcc’s of this molecule. We identify other four dihedral degrees of freedom as key parameters to predict accurate hcc values of the HOOCHβNOC moiety: the ∠HβCNC, ∠OOCN, and ∠HOOC dihedral angles (φ is defined as 180°, ∠HβCNC and the other two angles are, respectively, labeled φ1 and φ2) and the ∠ONCC improper dihedral angle (θ) quantifying the nitrogen out-of-plane displacement. In line with our previous study,36 a harmonic improper dihedral potential is selected to describe ∠ONCC, i.e., Uimp ) kimpθ2. The force field torsional, improper, and repulsive parameters are assigned on the basis of the quantum PES as a function of (θ, φ) of the IPENO-OOH nitroxide (see Figure 2), computed at the B3LYP/6-311+G(d,p) level of theory. In Figure 4, that PES is compared to the one obtained using the most reliable set of force field parameters. The agreement between the two surfaces is particularly good in the region (150° < θ < 210°, 30° < φ < 120°), which turns out to be particularly explored along the solvated trajectories of DMPO-OOH. Concerning φ1, we refined the torsional parameters to reproduce the B3LYP/6-311+G(d,p) energy profiles as a function of that angle for DMPO-OOH in both 3T4 and 4T3

Figure 4. Potential energy surface of IPENO-OOH vs the improper dihedral angle θ and the dihedral angle φ, based on quantum B3LYP/ 6-311+G(d,p) (up) and force field (down) computations. The contour values are in kcal mol-1.

conformations (φ2 and θ values were set to 180°). The quantum profiles are shown in Figure 5. Whatever the conformer, the profiles are characterized by three minima corresponding to φ1 values of 90, 170, and 290°: the deepest minimum corresponds to φ1 ) 290° (with the OOH moiety above the 5-membered ring) and is more stable by 1-2 kcal mol-1, relative to the other two minima. Note that these latter two minima are almost isoenergetic for the 3T4 conformation, whereas, for 4T3, the minimum corresponding to φ1 ) 90° is more stable by 1 kcal mol-1, relative to the φ1 ) 170° one. Whatever the conformation, the first two minima are separated by an energy barrier of 2 kcal mol-1, whereas the barrier is at least twice larger between

Properties of the Superoxide Radical Adduct of DMPO

Figure 5. Potential energy and free energy profiles vs the dihedral angle OOCN for DMPO-OOH. (a and b) Energy profile in vacuo for DMPO-OOH in its 3T4 and 4T3 conformations, respectively (blue line, TCPEp results; black line, B3LYP/6-311+ G(d,p) results; dashed black line, MP2/6-311+G(2df,2p)//B3LYP-6-311+G(d,p) results). (c) Free energy profile in solution computed using an umbrella sampling approach. c1 and c2 refer to the two DMPO-OOH conformation sites in slow exchange discussed in the text.

the second and the last ones (4 and 6 kcal mol-1, respectively, for the 3T4 and 4T3 conformers). Compared to the MP2/6311+G(2df,2p) results, these barriers are weaker by 1-2 kcal mol-1; however, the energy differences between the minima are identical regardless of the level of theory. The close features observed for the profiles corresponding to 3T4 and 4T3 conformations clearly suggest that two forms of DMPO-OOH may exist in solution regardless of the 5-membered ring conformation: the conformer family c1 corresponding to a φ1 angle included within 60 and 210° and the c2 one corresponding to φ1 ≈ 290°. This hypothesis is discussed in detail below. Concerning the force field profiles, they agree with the quantum ones for both the 3T4 and 4T3 conformers. However, for the 3T4 conformer, the force field profile is shifted to larger φ1 values by 5°, relative to the B3LYP one. For the 4T3 conformer, our force field approach predicts the minima corresponding to φ1 ) 90° and φ1 ) 170° to be isoenergetic (similarly to the 3T4 conformer case), whereas quantum computations have shown the φ1 ) 90° minimum to be more stable by 1 kcal mol-1 than the second one. However, this small difference does not impact our conclusions about the study of DMPO-OOH in solution. Concerning φ2, we refined the torsional parameters similarly to reproduce the B3LYP/6-311+G(d,p) energy profiles as a function of that angle for DMPO-OOH in both 3T4 and 4T3 conformations (φ1 and θ values were set, respectively, to 80

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11797 and 180°; profiles are provided as Supporting Information). The 3 T4 and 4T3 quantum energy profiles are close: they present a deep minimum at φ2 ) 70-80° (corresponding to the intramolecular HB COOH · · · ON), which is more stable by 3 (3T4) and 5 (4T3) kcal mol-1, as compared to a flat energy region located at 180-270° (region principally explored along the solvated trajectories). Note that the 3T4 force field and quantum energy profiles agree within 0.1 kcal mol-1, assessing the quality of our force field. B. Molecular Dynamics Details, Free Energy Calculation, Extracted Clusters, and ESPF Method. The initial DMPO-OOH structures corresponding to various conformers for the solvated trajectories were randomly taken along gas phase trajectories of 10 ns. However, from the above discussion regarding the possible role played by the c1 and c2 conformation sites in solution, we take care to select eight structures corresponding to a dihedral angle of ∠OOCN e 210° and the other eight ones corresponding to ∠OOCN g 240°. These structures were embedded into cubic boxes containing about 1000 water molecules. Hence, two series (one for each conformer site) of eight MD trajectories of 8 ns duration were performed; they were generated using each of the eight force fields discussed above. MD details are almost identical to those of our recent studies concerning small solvated nitroxide34-36 (details are also available as Supporting Information); in particular, the equations of motion were solved using a multi-time-step propagator devoted to polarizable force fields.39 However, we relaxed the volume of the simulation boxes by performing a first series of NPT simulations of 2 ns duration (300 K, 1 atm). Once these initial runs were achieved, the cell volume converged: it increased by about 1% on average along the trajectories. The simulations were then restarted in the NVT ensemble by considering a cell volume averaged over the last 50 ps of NPT simulations. For a given force field parameter set, we have also computed the free energy profile of the DMPO-OOH in solution as a function of its dihedral angle ∠OOCN using an umbrella sampling protocol, made of 24 windows. Each window consists of a 4 ns NPT MD trajectory at ambient conditions, with the angle ∠OOCN restrained to a selected value (included within 0 and 345°), based on a harmonic potential (whose force constant is set to 10 kcal mol-1 rad-2, ensuring a sufficient overlap of the adjacent histograms40). The MD simulations were unbiased using a WHAM scheme41 integrated in our POLARIS(MD) code (the convergence criterion and the histogram bin size are, respectively, set to 10-6 kcal mol-1 and 2°). The solvated nitroxide structures taken into account in QM computations correspond to clusters containing the nitroxide and its first hydration shell, including the water molecules located at less than 3.6 Å from the three oxygen atoms of the NO and COOH moieties, i.e., the nitroxide atoms that can be involved in HBs with water. This distance corresponds to the average position of the first minimum located after the first peak of water oxygen radial distribution functions around the three oxygen atoms of the nitroxide. The influence of the bulk water molecules which do not belong to the clusters was accounted for in QM computations by considering the electrostatic potential Vi they generate on each cluster atom i, by means of the ESPF method.42 In brief, this method consists of altering the matrix elements of the Fock operator of a system to afford the influence of any kind of external electrostatic potential Vi. Contrary to the polarizable continuum model (PCM), it allows one to consider the instantaneous anisotropic electrostatic potential arising from the

11798

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Houriez et al.

bulk water molecules along the trajectories. In the context of hcc computations, we have shown that this method coupled with a polarizable force field yields solvent shift estimates in agreement with experiment.36 C. Quantum Computations of Average hcc’s. The isotropic hcc’s aX of a monoradical (corresponding to the Fermi contact terms) are computed by considering them proportional to the electronic spin density FsX evaluated at the nucleus X

aX )

8π g µ gXµXFX 3 e e n n s

(4)

where gnX and µnX are, respectively, the Lande´ g factor and the Bohr magneton of the nucleus and ge and µe are the g factor and the magneton of the electron.43 Our previous studies34-36 have shown that experimental nitroxide nitrogen hcc’s can be reproduced with a good accuracy by considering the PBE0/6-31+G(d) level of theory. Compared to the much more resource-demanding QCISD level (used in conjunction with purposely tailored basis sets, such as the Chipman basis set44 or the EPR-II one45), this effective level of theory is able to reproduce within 0.5 G on average the nitrogen hcc’s of small solvated nitroxides. In the present study, the PBE0/6-31+G(d) level of theory has been systematically considered to compute the hcc’s of the nitroxide structures extracted along the solvated trajectories. However, for comparison purposes, the CCSD/EPR-II level of theory was also considered to compute the hcc’s of the MENO-OOH radical (see Figure 2), in particular to investigate the conformational dependence of the HOOCHβNO moiety hcc’s. All quantum computations were performed using the Gaussian 03 package of programs,46 in which we locally incorporated the ESPF method (via an additional link). The average hcc values were computed on the basis of statistical ensembles including 2000 structures, extracted along the solvated trajectories. The sampling interval was 1 ps, as an interval of 500 fs provided a set of temporally uncorrelated structures for nitrogen hcc computation in the case of the dimethyl nitroxide.34 The sampling starts once an equilibration period of 4 ns is achieved and goes on along a trajectory segment of 2 ns. III. Results and Discussion A. Geometrical Parameters Inpacting the HOOCHβNO Moiety hcc’s: A Gas Phase Study. The identification of the key intramolecular degrees of freedom influencing the DMPO-OOH hcc’s was initiated by studying the dependence of the gas phase hcc’s of the MENO-OOH model radical (see Figure 2) on the dihedral angles ∠HβCNC and ∠OOCN, as well as on the improper dihedral angle ∠ONCC. As above, φ is defined as 180° - ∠HβCNC and the other two angles are denoted by φ1 and θ, respectively. The hcc’s have been computed at the CCSD/EPR-II level of theory, considering three sets of radical structures corresponding to φ1 values of 80, 180, and 300°, respectively. In the three sets, the considered structures correspond to regularly spaced θ and φ values (within 130 and 230° and within 0 and 180°, respectively). The dihedral angle ∠HOOC was always set to 180°. In accordance with our previous studies, all the corresponding geometries have been optimized at the B3LYP/6-311+G(d,p) level of theory. Selected hcc contour plots (corresponding to nitrogen, hydrogen Hβ, and one of the γ-hydrogen atoms) are shown in Figure 6. Note that the contour plots corresponding to φ1 ) 80 and 300° are very close.

Figure 6. Nitrogen (up), β-hydrogen (middle), and γ-hydrogen (down) hcc values vs the improper dihedral angle θ and the dihedral angle φ for MENO-OOH, computed at the CCSD/EPR-II//B3LYP/6-311+G(d,p) level of theory. On the left, the values correspond to φ1 ) 180°, and on the right, to φ1 ) 80°. For φ1 ) 300°, the contour plots are very close to those corresponding to φ1 ) 80° (see the Supporting Information). The contour values are in Gauss.

The nitrogen hcc depends mostly on the sole θ value. Whatever φ, the nitrogen hcc profile as a function of θ is similar to the ones already reported for small nitroxides, such as the dimethyl nitroxide:47 the hcc values are minimum (about 11 G) for θ values close to 180° (planar nitroxide), while they increase up to 21 G when θ evolves toward 130 or 230° (bent nitroxide). The hydrogen Hβ hcc depends on both θ and φ angles. The largest hcc values (>16 G) are found in a region located around θ ) 165° and φ ) 105°, whereas the smallest ones (as low as 0 G) are obtained for φ values smaller than 30°. However, an overall small dependence of the Hβ hcc on φ1 exists in the (140° < θ < 190°, 90° < φ < 120°) domain, corresponding to the largest Hβ hcc values. Interestingly, the region in which this hcc is larger than 16 G is smaller for φ1 ) 80° than for φ1 ) 180°. Actually, this area corresponds to Hβ in a pseudoaxial position, with respect to the ONCC mean plane, resulting in the strongest hyperconjugative interactions between the orbitals π/π* of the NO moiety and the antibonding orbital σ*(C-Hβ) (anomeric effect). As a consequence, the electronic spin density at the Hβ nucleus is maximum in this (θ, φ) region, leading to a large Hβ hcc value. On the other hand, the above-mentioned orbital interactions are weak in the region φ < 30° (i.e., Hβ in a pseudoequatorial position), explaining the negligible Hβ hcc value observed there.

Properties of the Superoxide Radical Adduct of DMPO Among the three γ-hydrogen atoms of MENO-OOH, two are almost silent, i.e., their hcc’s are very weak, in particular in the domain (150° < θ < 210°, 30° < φ < 120°), which turns out to be particularly explored along the solvated trajectories of DMPO-OOH. In absolute value, they are all smaller than 0.5 G, and they are usually of 0.2 G. These two γ-hydrogen atoms correspond to dihedral angles ∠HγCCN of (60°. On the contrary, the hcc of the third γ-hydrogen atom (whose dihedral angle HγCCN is 180°) can reach significant values (up to 2 G). Similar to Hβ, the hcc of the latter Hγ depends on both θ and φ angles and the largest hcc values are found in a region centered at (θ ) 150°, φ ) 30°). Such high hcc values for a hydrogen atom in the γ position can be interpreted as conjugated singly occupied molecular orbital delocalization and spin polarization effects (see the interesting discussion in ref 48). The influence of the dihedral angle ∠HOOC (φ2) on the latter hcc values has been investigated by computing the hcc profiles at the CCSD/EPR-II level of theory as a function of φ2. To this end, we consider a series of MENO-OOH conformations corresponding to Hβ in a pseudoaxial position (i.e., ∠HβCNC ) 90°). Again, the structures (corresponding to φ1 ) 80° and to regularly spaced φ2 values within 0 and 350°) were optimized at the B3LYP/6-311+G(d,p) level of theory (profiles provided as Supporting Information). In the 60° < φ2 < 120° domain, which is characterized by an intramolecular HB OOH · · · ON, the hcc’s of the nitrogen and the third hydrogen Hγ depend weakly on φ2, whereas, in the 120° < φ2 < 270° region, the dependence is much less pronounced and may be considered as negligible. As the latter region is mainly explored along the solvated trajectories of DMPO-OOH, it appears thus that the dihedral angle ∠HOOC is not expected to influence the hcc values of the solvated DMPO-OOH. To check the ability of the cheaper PBE0/6-31+G(d) level of theory to give hcc’s of similar quality on average to the ones computed at the CCSD/EPR-II high level of theory, the abovementioned hcc surfaces and profiles have been recalculated (they are provided as Supporting Information). Overall, the main features of the CCSD surfaces are well reproduced. Moreover, this DFT level of theory is able to provide accurate enough hcc’s: the maximum deviation between both the latter levels of theory is 2 G (e.g., for Hβ in the axial position), while, most of the time, the discrepancy is included within 0.2-0.6 G (in absolute values). Together, all of these results demonstrate the reliability of the PBE0/6-31+G(d) level to compute effective hcc’s, not only for the nitrogen nucleus but also for the hydrogen ones. B. DMPO-OOH Conformations in Aqueous Solution. All the statistical averages and distribution functions discussed in the present section were computed along 4 ns segments of the 16 NVT MD trajectories of the solvated DMPO-OOH, denoted by T1 to T8 for c1-type starting structures and by T1′ to T8′ for c2-type ones. Note that deriving averages from 2 ns segments or 4 ns ones leads to indistinguishable results. However, along two trajectories corresponding to a c2-type conformer starting structure, a transition of the DMPO-OOH structure toward a c1-type form was observed (the transition occurs at different simulation times). Such a c1-c2 transition is not observed along the other 14 trajectories. Hence, a new series of eight trajectories corresponding to c2-type starting structures was generated, by adding to the potential energy function a harmonic potential kφ1(φ1 - φ10)2 to prevent any c1-c2 transition. kφ1 was set to 10 kcal mol-1, which allows one to reproduce the φ1 spread observed along the unrestrained trajectories. φ10 was set to the

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11799

Figure 7. Cremer-Pople ring phase ω distribution functions along selected c1 trajectories corresponding to five different force field parameter sets.

average value of φ1 observed along the c2 trajectories, before any c1-c2 transition occurs. In Figure 7, the 1D distribution functions gω corresponding to the phase ω degree of freedom quantifying the 5-membered ring conformation have been plotted for five selected c1 trajectories. Note that very close gω distributions have been obtained from c2 trajectories. They are all characterized by two well-defined peaks at 70-90 and 250-270°, and by a flat and almost null minimum at 180° (the distribution spread in the vicinity of the maxima goes up to 60°). Moreover, a clear relationship between the maxima appears: if the first peak for some trajectories is located around 70°, the second one is rather around 270°, whereas for some other trajectories the maxima are around 90 and 260°. Hence, a 3T4/E3 or a 3E/4T3 ring conformer equilibrium emerges from this analysis: the first one is observed along trajectories T1(T1′) to T4(T4′) and the second one for T5(T5′) to T8(T8′). As mentioned in section II.A.3, these two equilibria are referred to a 3T4/4T3 equilibrium in the following. The populations (p3T4, p4T3) corresponding to this latter equilibria were evaluated by integrating the gω profiles along the trajectories: they vary from about (5%, 95%) [T1(T1′)] to about (95%, 5%) [T8 (T8′)]. The p3T4 values are summarized in Table 2. Moreover, on the basis of the temporal fluctuations of the phase ω along the trajectories T4 and T5′, we estimate for DMPO-OOH the average time spent in a 3T4-like form to be about 10 ps before interconversion toward a 4T3-like form (and vice versa). The 3T4/4T3 equilibrium echoes in the shape of the 2D distribution functions g2D(θ, φ) (corresponding to angles ∠ONCC and ∠HβCNC), computed along the trajectories. For instance, the distributions corresponding to three sets of trajectories, namely, (T1,T1′), (T4,T5′), and (T8,T8′), are plotted in Figure 8. The g2D(θ, φ) is centered at about (180°, 45°) for the first trajectory set (p3T4 ≈ 95%), whereas it is centered at about (180°, 75°) for the third set (p3T4 ) 5-8%; the distribution spread is larger for θ in this latter case). Hence, Hβ can be considered as being in an almost axial position along the trajectories (T1,T1′), while it is equatorial in (T8,T8′). Concerning the improper dihedral angle ONCC, its 1D distribution functions correspond to Gaussian profiles centered at 180°, with almost negligible contributions for θ < 150° and θ > 210°, similar to the case of the PROXYL nitroxide we recently investigated.36 The larger distribution spread for p3T4 e 8%, relative to p3T4 ≈ 95%, can be readily explained by considering the IPENO-OOH PES plotted as a function of (θ, φ) in Figure 4: for φ values greater than 60°, the PES is flatter than for smaller angle values. Finally, in the case of the intermediate ring conformer distribution (p3T4 ≈ 50% along T4 and T5′),

11800

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Houriez et al.

TABLE 2: Average hcc Values along the Eight Solvated Trajectories Corresponding to the Sites c1 and c2 (Respectively Labeled T1-T8 and T1′-T8′)a trajectory T1 T2 T3 T4 T5 T6 T7 T8

T1′ T2′ T3′ T4′ T5′ T6′ T7′ T8′

p3T4

ajHβ

ajHγ

ajN

0 4.9 17.2 26.3 45.7 59.3 78.3 85.5 96.4 100 0 6.2 18.3 25.1 41.4 54.1 63.0 80.9 92.0 100

15.9 [14.8] 2.2 15.7 [14.6] 2.2 13.7 [12.7] 1.8 12.7 [11.8] 1.7 10.8 [10.0] 1.6 9.5 [8.8] 1.4 7.4 [6.8] 1.1 6.5 [6.0] 0.9 5.0 [4.6] 0.7 4.8 [4.5] 0.7 15.6 [14.4] 3.7 15.3 [14.3] 3.6 14.0 [13.0] 3.2 13.5 [12.6] 3.2 12.4 [11.6] 3.0 11.4 [10.6] 2.9 10.4 [9.6] 2.5 9.0 [8.3] 2.2 8.1 [7.4] 2.0 7.3 [6.7] 1.9

0.2 [0.0] -0.1 0.3 [0.1] -0.1 0.5 [0.3] 0.0 0.8 [0.6] 0.0 1.1 [0.9] 0.1 1.5 [1.3] 0.3 1.8 [1.6] 0.3 2.1 [1.9] 0.4 2.2 [2.0] 0.4 2.3 [2.1] 0.4 -0.1 [-0.3] -0.1 0.1 [-0.1] -0.1 0.4 [0.2] 0.0 0.5 [0.3] 0.0 0.8 [0.6] 0.1 1.1 [0.9] 0.3 1.3 [1.1] 0.3 1.7 [1.5] 0.4 2.0 [1.8] 0.4 2.1 [1.9] 0.4

14.2 [14.9] 2.3 14.2 [14.9] 2.4 14.0 [14.7] 2.4 13.9 [14.6] 2.4 13.6 [14.3] 2.5 13.5 [14.1] 2.5 13.3 [14.0] 2.6 13.2 [13.8] 2.6 13.2 [13.8] 2.5 13.1 [13.7] 2.6 15.2 [15.9] 3.2 15.0 [15.7] 3.3 14.9 [15.6] 3.1 14.7 [15.3] 3.2 14.0 [14.7] 3.2 14.4 [15.1] 3.2 14.1 [14.8] 3.1 14.0 [14.7] 3.1 14.0 [14.6] 3.2 13.8 [14.4] 3.2

a p3T4: 3T4-like conformer population values in percent along each trajectory. All hcc ajX values are in Gauss, as estimated from PBE0/ 6-31+G(d) computations (in brackets, CCSD/EPR-II extrapolated values; in italics, the solvent contribution δS(ajX) to the hcc’s based on PBE0/6-31+G(d) results). For p3T4 ) 0/100%: extrapolated hcc values based on linear regressions of the trajectory averages.

Figure 8. 2D distribution functions g2D(θ, φ) computed along c1 and c2 solvated trajectories corresponding to populations p3T4, respectively, of about 5% (up), 50% (middle), and 95% (down). Orange and blue contours correspond, respectively, to weak and high g2D(θ, φ) values.

g2D(θ, φ) can be considered as the average of the other two distributions. To further characterize the DMPO-OOH structures in solution, some other 1D distribution functions were computed

along each solvated trajectory, corresponding to the dihedral angle OOCN (φ1), the dihedral angle HOOC (φ2), and the distance RHO between the hydroperoxyl hydrogen atom and the nitroxide oxygen one (plots of these functions are provided as Supporting Information). Concerning the φ1 distributions computed along the c1 trajectories, they are marked by a single sharp peak at about 80°, with a spread of 60° showing that the geometry of the OOCN moiety does not depend on the ring conformation. In the gas phase, and regardless of the ring conformation, the energy difference between the structures with φ1 values of 80 and 170° is at most 1 kcal mol-1. In solution, the structure corresponding to φ1 ) 80° is largely stabilized due to the existence of a stable COOH · · · water · · · ON HB network observed along all the c1 trajectories (see Figure 3). Note that the above-mentioned disagreement between our force field approach and the quantum computations concerning these two structures in the gas phase (our approach predicts them to be isoenergetic whatever the ring conformation) is not expected to affect the reliability of our conclusions concerning DMPO-OOH in water. Along the c2 trajectories, as expected by considering a harmonic potential to restrain the φ1 angle (see above), all the φ1 distributions present a single Gaussian-like peak centered at about 300°, with a spread of 50°. To evaluate the stability of the c1/c2 conformers in solution (especially for the c2-type conformers), we have computed the free energy profile of DMPO-OOH as a function of φ1 (in the domain 0° < φ1 < 345°) according to the umbrella sampling protocol presented in section 2.B. We considered here the force field parameter set leading to a 1:1 3T4/4T3 equilibrium ratio for the two conformer sites c1 and c2 (i.e., the parameter set used to generate the trajectories T4 and T4′). The free energy profile is plotted in Figure 5. It shows that the conformer sites c1 and c2 are separated by two energy barriers (corresponding to φ1 ≈ 10 and 230°) estimated in solution to be of 7 kcal mol-1. Compared to vacuum energy profiles for conformers 3T4 and 4 T3, the free energy barrier at 230° is larger in solution by 1-3 kcal mol-1. Moreover, in water, the conformations of c1 corresponding to φ1 ≈ 80° are slightly more stable than the c2 ones, by 0.5 kcal mol-1. The second c1 conformation, corresponding to φ1 ≈ 180°, is strongly destabilized in solution (see the discussion above). Hence, these results demonstrate that two conformation sites of the DMPO-OOH have to be mainly observed in solution. Because of the height of the energy barrier connecting them in solution, the c1 and c2 sites are in slow chemical exchange. Concerning the HOOC (φ2) distribution functions, their profiles are very close whatever the trajectory and the conformer family: they present two peaks located around 140 and 220°, separated by a mimimum at about 180-200°. In c1 trajectories, the peak heights are close, whereas, in c2 trajectories, the height of the first peak is 4 times greater than the second one. As already discussed, the nitrogen and hydrogen hcc dependence on φ2 in that latter domain is negligible. As mentioned in section I, Villamena and co-workers based their studies about nitroxide-type spin adducts in water on the assumption that some intramolecular HBs may favor the corresponding conformations. This assumption being rather counterintuitive, the percentage of DMPO-OOH structures in which an intramolecular HB occurs in aqueous solution has been estimated from the RHO distributions (see the Supporting Information). Actually, a HB is assumed to exist when a simple geometric criterion is met: RHO < 2.5 Å. By integrating the latter distribution functions, we note that such an intramolecular bond

Properties of the Superoxide Radical Adduct of DMPO

Figure 9. c1 (empty symbols) and c2 (plain symbols) hcc values vs the population p3T4: (a) nitrogen (circles) and hydrogen Hβ (diamonds) hcc’s; (b) H1γ (circles) and H2γ (diamonds). Plain and dashed lines: linear regression fits.

is observed at most in 4% of the solvated DMPO-OOH structures along the c1 trajectories. However, such HBs are observed only along the trajectories corresponding to a p3T4 value smaller than 50%, whereas it corresponds to a rare event for the remaining trajectories. Hence, in aqueous solution, the less constrained COOH · · · water · · · ON HB network is in any case favored against any DMPO-OOH intramolecular HB, invalidating Villamena’s assumption. In conclusion, together, all of these results demonstrate that the DMPO-OOH structure in an aqueous medium can be readily interpreted as a chemical exchange between the two conformation sites c1 and c2 corresponding to a dihedral angle ∠OOCN, respectively, of 80 and 300°. Moreover, each family involves a fast interconversion process between the two 5-membered ring major conformations, i.e., the 3T4 and 4T3 conformers shown in Figure 3. These conformers do not present any intramolecular COOH · · · ON HB in solution, even if the dihedral angle ∠OOCN value is 80° on average in the c1 family. This structural analysis is completely new, and it will lead to providing reliable arguments to interpret the EPR spectrum of DMPO-OOH in water, as highlighted in the next sections. C. DMPO-OOH hcc’s in Solution. 1. hcc Values Ws DMPO-OOH Conformer Equilibrium Ratios. The average hcc values ajN, ajHβ, and ajHγ1, computed according to our QM/MM protocol at the PBE0/6-31+G(d) level of theory, are summarized in Table 2 for both types of trajectories. They are also plotted as a function of p3T4 in Figure 9, together with the hcc values ajHγ2 of the second γ-hydrogen (see Figure 3 for the γ-hydrogen labeling). The root-mean-square deviations corresponding to these average values are included within 4 and 6 G (cf. also the distributions of nitrogen, hydrogen β and γ1 hcc values along the trajectories corresponding to sites c1 and c2 provided as Supporting Information). Assuming that the sets of extracted solvated nitroxide structures are temporally uncorrelated, the uncertainty affecting these average values obtained from statistical ensembles including 2000 structures is estimated to be at most 0.2 G. Note that we also computed average hcc values

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11801 along a single trajectory by considering a twice larger statistical ensemble (and extracted along a 4 ns segment): the new values agree with those derived from 2000 structures within less than 0.2 G. All the average hcc values depend linearly on p3T4: the corresponding correlation coefficients are all larger than 0.98. Hence, these results demonstrate that the magnitude of the DMPO-OOH spin densities in aqueous solution is intimately related to the properties of the 3T4/4T3 equilibrium described in section III.B, regardless of the c1 and c2 sites. The strongest dependence on this ratio concerns ajHβ: it varies from 15.9 to 4.8 G (trajectories c1) and from 15.6 to 7.3 G (trajectories c2), respectively, when p3T4 increases from 0 to 100%. A weaker dependence is observed for aj N, aj Hγ1 , and aj Hγ2 : they vary by 1-2 G. A second noticeable feature concerns the nitrogen and β-hydrogen hcc values: the c1 hcc’s are systematically smaller compared to the c2 ones if p3T4 values are larger than 5%. Concerning the nitrogen nucleus, the hcc difference between both sites is almost constant (about 0.9 G), whereas it linearly increases for the β-hydrogen nucleus from about 0.0 to 2.5 G as p3T4 varies from 0 to 100%. Concerning the γ-hydrogen nuclei, the opposite trend is observed: the c1 hcc’s are larger than the c2 ones. The difference in Hγ1 hcc values is almost constant (about 0.3 G), whereas it linearly increases from 0.1 to 0.7 G for Hγ2 . Note that only Hγ1 gives rise to a significant hcc value (>1 G), in agreement with the gas phase analysis presented in section III.A. Its position is characterized by a dihedral angle ∠HγCCN close to 180° (W arrangement), only in 3T4-like DMPO-OOH conformers (regardless of whether they are of c1 or c2 type). The above differences in β-hydrogen and nitrogen hcc’s originate from the properties of the hcc surfaces shown in Figure 6 and from those of the 2D distribution function g2D(θ, φ) shown in Figure 8. As compared to c1 trajectories, the distributions g2D(θ, φ) spread over a wider θ value domain and they are shifted to larger φ values (by about 5-10°) in the case of large p3T4 values along c2 trajectories. This can be interpreted as resulting (1) from weaker COOH-NO interactions in site c2 than in c1, reducing the intramolecular constraints affecting the motion of the oxygen atom of the NO moiety, and (2) from strong interactions between the COOH moiety and the 5-membered ring occurring in site c2, which in turn affect the motion of the β-hydrogen atom. In our previous works, the uncertainty related to the quantum level of theory selected to compute the hcc’s was estimated. Here, the same procedure is applied, by calculating the quantum uncertainties δaXquantum as the expected difference in hcc values between the effective PBE0/6-31+G(d) and the CCSD/EPR-II levels of theory using an interpolation scheme. As discussed above, the dihedral angle OOCN value is, respectively, 80 and 300° on average along the c1 and c2 trajectories (and it does not depend on the DMPO-OOH ring conformation). Therefore, the CCSD hcc’s of the model nitroxide MENO-OOH in the gas phase have been interpolated using a bicubic interpolation scheme based on the instantaneous values of the improper ∠ONCC and dihedral ∠HβCNC angles observed along the solvated DMPO-OOH trajectories and on the hcc surfaces at both CCSD/EPR-II and PBE0/6-31+G(d) levels of theory corresponding to ∠OOCN ) 80 and 300° (the surfaces are plotted in Figure 6 or provided as Supporting Information). The difference between both sets of MENO-OOH interpolated hcc values provides the δaXquantum values, then allowing one to compute the estimated corrected CCSD hcc’s of DMPO-OOH

11802

J. Phys. Chem. B, Vol. 114, No. 36, 2010

Houriez et al.

in solution (they are summarized in Table 2). Compared to PBE0 values, the CCSD ones are larger (N) and weaker (Hβ and Hγ1 ) by about 0.7, 0.9, and 0.2 G, respectively. Note that the linear dependence of the hcc values on p3T4 is preserved, demonstrating that the conclusions derived from PBE0 results are not affected by the quantum level of theory used to compute the hcc’s. Again, the aqueous DMPO-OOH hcc’s characterizing the two c1 and c2 sites depend linearly on the 3T4/4T3 equilibrium ratio. 2. SolWent Contribution to the hcc’s. The influence of the solvent on the DMPO-OOH hcc values can be evaluated by computing the hcc mean values ajXn along each trajectory while considering the same statistical ensembles of 2000 structures from which the water molecules are discarded. Accordingly, the solvent contributions δS(ajX) to the average hcc’s ajX in solution are calculated as follows

δS(ajX) ) ajX - ajXn

(5)

The computed values are reported in Table 2. Except for Hγ, all these values are positive; the solvent is thus responsible for an increase of the hcc’s. However, the δS(ajX) magnitude depends on the nucleus nature, on the site type c1/c2, and on the 3T4/4T3 equilibrium ratio, with the exception of N, for which the solvent contribution appears to be constant whatever the trajectory: about 2.5 and 3.2 G in sites c1 and c2, respectively. These values are very close to those we reported recently for six small nitroxides.35,36 The contribution δS(ajHβ) decreases from 2.2 to 0.7 G (site c1) and from 3.7 to 1.9 G (site c2) as p3T4 increases from 0 to 100%, whereas the opposite trend is observed for δS(ajHγ1): it increases from -0.1 up to 0.4 G for both sites c1 and c2. Note that, for Hγ2 , the solvent contribution is constant and very weak: about -0.1 G. The local solvent structure in the vicinity of the nitroxide nuclei has been analyzed in terms of interatomic radial distribution functions (some examples are provided as Supporting Information). No particular relationship has been evidenced between the solvent structure and the magnitude of the solvent contribution to the hcc’s. For instance, the mean number of water oxygen atoms around the NO oxygen within less than 3.6 Å is 3.1 ( 0.1 whatever the trajectory and the conformer site. However, some differences in the solvent angular distribution functions at the vicinity of the nitrogen and β-hydrogen nuclei exist between the two DMPO-OOH sites, originating from the specific orientation of the nitroxide COOH moiety in both confomer types. A detailed analysis of the incidence of solvent structure on the differences of electrostatic potential affecting the DMPO-OOH nuclei (which are the main origin of the solvent contribution to the nitroxide hcc’s) is currently under investigation, and the conclusions will be presented elsewhere. D. Comparison to Experiment. As shown above, whatever the considered site, the DMPO-OOH nitrogen, β-hydrogen, and γ-hydrogen hcc’s linearly depend on the 3T4/4T3 equilibrium ratio. Given a particular force field parameter set to handle the intramolecular nitroxide degrees of freedom, the c1 and c2 3T4/ 4 T3 equilibrium ratios are usually found to be close within 5%. Accordingly, with these c1 and c2 sites being energetically separated by a large enough barrier (about 7 kcal mol-1), one may rely on their distinct existences in solution as the major contributors to the EPR spectrum of DMPO-OOH. In other words, our simulations support the chemical exchange interpretation of the DMPO-OOH EPR spectrum asymmetry. More interestingly, the asymmetry is certainly caused by different hyperfine coupling constants at all the considered nuclei, even

at N, a case that was not considered in the numerical models used to fit the experimental spectrum. On the basis of this hypothesis and keeping in mind the uncertainties affecting our hcc estimates, the inspection of the plots shown in Figure 9 shows that a 3T4/4T3 equilibrium ratio close to 1:1 for each site c1 and c2 leads to nitrogen, β-hydrogen, and γ-hydrogen hcc values in good agreement with the experimental models B1, B2, and B3 (see Table 1). However, it is worth mentioning that there is no reason to suppose that both sites should be characterized by strictly equal 3T4/4T3 ratios. More precisely, a 1:1 equilibrium ratio can be described by c1 and c2 nitrogen hcc’s of about 14 G (but different in each site), by β-hydrogen hcc’s of about 11 and 10 G, respectively, and by detectable γ-hydrogen hcc’s around 1.0 G (concerning these hcc’s, the PBE0 values and the extrapolated CCSD ones agree within 0.5 G on average). It is noteworthy that our computations predict the c2 nitrogen and β-hydrogen hcc’s to be larger than the c1 ones and the contrary for the γ-hydrogen hcc’s, agreeing with the above-mentioned experimental models. A very small free energy difference (e.g., ∆F = 0.25 kcal mol-1) can alter a 50/50% equilibrium ratio in a 60/40% one. This ∆F value is twice smaller than the quantum uncertainty affecting the differences in energy ∆E between 3T4 and 4T3 conformers (corresponding to various dihedral angle ∠OOCN values): about 0.5 kcal mol-1. This estimate has been computed by comparing ∆E values based on B3LYP and MP2 levels of theory,usedinconjunctionwith6-311+G(d,p)and6-311+G(2df,2p) basis sets, respectively. As our force field parameters are assigned from quantum results, this explains why we are not able to propose a more precise statistical weight for each conformation family. However, we have compared the latter quantum ∆E values to those computed using the eight force field parameter sets used to generate the c1 and c2 solvated trajectories. Interestingly, the best agreement between quantum and molecular modeling values is obtained for force field parameter sets leading to a 3T4/4T3 equilibrium in solution whose ratio is included within 50/50% and 80/20% (see the details provided as Supporting Information). Note that the best agreement for the ∆E between force field and quantum MP2 computations is observed for the parameter set leading to a ratio of 60/40%, whereas the best agreement with B3LYP computations corresponds to a force field parameter set leading to a ratio of 80/20%. This insures the reliability of our QM/MM protocol, which is thus able to provide a realistic and useful guess of parameters to fit the EPR spectra of large nitroxides in complex environments. IV. Conclusion In this work, we reported the first complete theoretical investigation of the structures and hyperfine coupling constants of the DMPO-OOH nitroxide spin adduct in explicit water. Our multiscale approach combines the statistical investigations of the solvated nitroxide potential energy surface, by using an anisotropic polarizable force field, and the effective computation of isotropic hcc’s evaluated at different nuclei, by means of QM/ MM calculations based on the ESPF coupling method. In agreement with the main experimental interpretation but contrary to previous theoretical works, our results clearly show that there is no unique (i.e., long-lived) DMPO-OOH structure in aqueous solution on which the spectromagnetic properties can rely. Rather, a chemical exchange between two conformation sites has to be considered, each site involving mainly an equilibrium between 3T4-like and 4T3-like conformations of the DMPO-OOH 5-membered ring. Accordingly, the apparent

Properties of the Superoxide Radical Adduct of DMPO DMPO-OOH hcc’s obtained from the EPR spectrum can be interpreted as the averages of the pure 3T4 and 4T3 hcc’s at each site. Owing to a probable 1:1 distribution of the sites, hcc values are in very good agreement with the experimentally fitted ones. Because some room still exists to improve the agreement between the theoretical and experimental interpretations of the DMPO-OOH EPR spectrum, the current work can be used as a new guess for hcc parameters that can be tested against the available spectra. This work establishes the reliability of the simulation protocol we set up to get quantitative estimates of nitroxide hcc’s in an aqueous medium. It will thus be applied to new challenging nitroxide spin adducts, such as the EMPO-OOH49 and DEPMPO-OOH22 ones, whose EPR spectra usually feature several contributions originating from different diasteroisomers. More generally, this work opens the road to the accurate determination of spectromagnetic properties of large radicals in a complex environment. Acknowledgment. All of the computations were performed on the large scale facilities of the “Centre de Calcul et de Recherche Technologique” (CCRT) of the French Nuclear Agency (CEA) and of the “Centre Informatique National de l’Enseignement Supe´rieur” (CINES) and at the “Centre Re´gional de Compe´tences en Mode´lisation Mole´culaire” (CRCMM) in Marseille. This work was granted access to the HPC resources of CCRT and CINES under the allocation 2009082226 made by GENCI (Grand E´quipement National de Calcul Intensif). Supporting Information Available: Force field details (concerning hydrogen bonding, anomeric effect, and intramolecular key degrees of freedom), energetical properties of DMPO-OOH conformers in vacuo, and structural properties of DMPO-OOH in solution. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Halliwell, B.; Guttenridge, J. Free Radicals in Biology and Medicine, 3rd ed.; Clarendon Oxford University Press: Oxford, U.K., 1999. (2) Sies, H. OxidatiVe stress, oxidants and antioxidants; Acad. Press: London, 1991. (3) Favier, A.; Cadet, J.; Kalyanaraman, B.; Fontecave, M.; Pierre, J. L. Analysis of Free Radicals in Biological Systems; Birkhauser Verlag: Basel, Boston, Berlin, 1995. (4) Wardman, P. Free Radical Biol. Med. 2007, 43, 995–1022. (5) Tordo, P. In Electron Paramagnetic Resonance; Gilbert, B. C., Atherton, N. M., Davies, M. J., Eds.; RSC: Cambridge, U.K., 1998; Vol. 16, Chapter “Spin trapping: recent developments and applications”, pp 116144. (6) Davies, M. J. In Electron Paramagnetic Resonance; Gilbert, B. C., Davies, M. J., Murphy, D. M., Eds.; RSC: Cambridge, U.K., 2002; Vol. 16, Chapter “Recent developments in spin trapping”, pp 47-73. (7) Berliner, L. Spin Labelling: Theory and Applications; Academic Press: New York: 1976. (8) Berliner, L. Spin Labelling: Theory and Applications II; Academic Press: New York, 1979. (9) Hubbel, W.; Cafiso, D.; Altenbach, C. Nat. Struct. Biol. 2000, 7, 735. (10) Kocherginsky, N.; Swartz, H.; Sentjure, M. Nitroxide Spin Labels: Reactions in Biology and Chemistry; CRC Press, Inc.: Boca Raton, FL, 1995.

J. Phys. Chem. B, Vol. 114, No. 36, 2010 11803 (11) Berliner, L. J. Biological Magnetic Resonnace: Spin Labelling in the next Millenium; Plenum Press: New York, 1998. (12) M. Engstrom, J. V.; Schimmelpfennig, B.; Ågren, H. J. Phys. Chem. B 2002, 106, 12354–12360. (13) Clement, J.-L.; Tordo, P. Electron Paramagnetic Resonance; RSC: Cambridge, U.K., 2007; Vol. 20, Chapter 2, pp 29-49. (14) Likhtenshtein, G.; Yamauchi, J.; Nakatsuji, S.; Smirnov, A.; Tamura, R. Nitroxides: Applications in Chemistry, Biomedecine and Materials Science; Wiley-VCH GmbH and Co.: Weinheim, 2008. (15) Brustolon, M., Giamello, E., Eds.; Electron Paramagnetic Resonance; J. Wiley and Sons, Inc.: New York, 2009. (16) Harbour, J.; Chow, V.; Bolton, J. Can. J. Chem. 1974, 52, 3549. (17) Buettner, G. Free Radical Res. Commun. 1990, 10, 11. (18) Rosen, G.; Beselman, A.; Tsai, P.; Pou, S.; Mailer, C.; Ichikawa, K.; Robinson, B.; Nielen, R.; Halpern, H.; MacKerell, A. J. Org. Chem. 2004, 69, 1321–1330. (19) Cle´ment, J.-L.; Ferre´, N.; Siri, D.; Karoui, H.; Rockenbauer, A.; Tordo, P. J. Org. Chem. 2005, 70, 1198. (20) Dikalov, S.; Jiang, J.; Mason, R. Free Radical Res. 2005, 39, 825. (21) Villamena, F. A.; Merle, J. K.; Hadad, C. M.; Zweier, J. L. J. Phys. Chem. A 2005, 112, 12607–12615. (22) Fre´javille, C.; Karoui, H.; Tuccio, B.; Moigne, F. L.; Culcasi, M.; Pietri, S.; Lauricella, R.; Tordo, P. J. Med. Chem. 1995, 38, 258–265. (23) Malloy, T.; Bauman, L.; Carreira, L. Top. Stereochem. 1979, 11, 97–185. (24) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027–2094. (25) Villamena, F. A.; Liu, Y.; Zweier, J. L. J. Phys. Chem. A 2008, 109, 6089. (26) Rockenbauer, A.; Gaudel-Siri, A.; Siri, D.; Berchadsky, Y.; Moigne, F. L.; Olives, G.; Tordo, P. J. Phys. Chem. A 2003, 107, 3851–3857. (27) Takase, H.; Kikuchi, O. Chem. Phys. 1994, 181, 57. (28) Yagi, T.; Kikuchi, O. J. Phys. Chem. A 1999, 103, 9132. (29) Pavone, M.; Benzi, C.; De Angelis, F.; Barone, V. Chem. Phys. Lett. 2004, 395, 120. (30) Pavone, M.; Cimino, P.; De Angelis, F.; Barone, V. J. Am. Chem. Soc. 2006, 128, 4338. (31) Pavone, M.; Cimino, P.; Crescenzi, O.; Sillanpa¨a¨, A.; Barone, V. J. Phys. Chem. B 2007, 111, 8928. (32) Pavone, M.; Sillanpaa, A.; Cimino, P.; Crescenzi, O.; Barone, V. J. Phys. Chem. B 2007, 110, 16189. (33) Neugebauer, J.; Louwerse, M. J.; Belanzoni, P.; Wesolowski, T. A.; Baerends, E. J. J. Chem. Phys. 2005, 123, 114101. (34) Houriez, C.; Ferre´, N.; Masella, M.; Siri, D. J. Chem. Phys. 2008, 128, 244504. (35) Houriez, C.; Ferre, N.; Masella, M.; Siri, D. THEOCHEM 2009, 898, 49. (36) Houriez, C.; Ferre´, N.; Siri, D.; Masella, M. J. Phys. Chem. B 2009, 113, 15047. (37) Masella, M.; Cuniasse, P. J. Chem. Phys. 2003, 119, 1866. (38) Cremer, D.; Pople, J. J. Am. Chem. Soc. 1975, 97, 1358. (39) Masella, M. Mol. Phys. 2006, 104, 415. (40) Souaille, M.; Roux, B. Comput. Phys. Commun. 2001, 135, 40– 57. (41) Kumar, S.; Bouzida, D.; Swendsen, R.; Kollman, P.; Rosenberg, J. J. Comput. Chem. 1992, 13, 1011. ´ ngya´n, J. G. Chem. Phys. Lett. 2002, 356, 331. (42) Ferre´, N.; A (43) Kaupp, M., Buhl, M., Malkin, V. G., Eds.; Calculation of NMR and EPR parameters; Wiley-VCH: Weinheim, 2004. (44) Chipman, D. M. Theor. Chim. Acta 1992, 82, 93. (45) Adamo, C.; Barone, V.; Fortunelli, A. J. Chem. Phys. 1995, 102, 384. (46) Frisch, M. J.; et al. Gaussian 03, revision.02; Gaussian, Inc.: Wallingford, CT, 2004. (47) Improta, R.; Barone, V. Chem. ReV. 2004, 104, 1231. (48) Barone, V.; Adamo, C.; Grand, A.; Subra, R. Chem. Phys. Lett. 1995, 246, 53–58. (49) Olive, G.; Mercier, A.; Moigne, F. L.; Rockenbauer, A.; Tordo, P. Free Radical Biol. Med. 2000, 28, 403–408.

JP1033307