Theoretical Determination of the Standard Reduction Potentials of

7210

J. Phys. Chem. B 2007, 111, 7210-7217

Theoretical Determination of the Standard Reduction Potentials of Pheophytin-a in N,N-Dimethyl Formamide and Membrane Nital Mehta and Sambhu N. Datta* Department of Chemistry, Indian Institute of TechnologysBombay, Powai, Mumbai 400076, India ReceiVed: NoVember 8, 2006; In Final Form: April 18, 2007

Quantum mechanical/molecular mechanics (QM/MM) calculations were performed on the neutral, anionic, and dianionic forms of Pheophytin-a (Pheo-a) in N,N-dimethyl formamide (DMF) in order to calculate the absolute free energy of reduction of Pheo-a in solution. The geometry of the solvated species was optimized by restricted open-shell density functional treatment (ROB3LYP) using the 6-31G(d) basis set for the molecular species while the primary solvent shell consisting of 45 DMF molecules was treated by the MM method using the universal force field (UFF). Electronic energies of the neutral, anionic, and dianionic species were obtained by carrying out single point density functional theory (DFT) calculations using the 6-311+G(2d,2p) basis set on the respective ONIOM optimized geometries. The CHARMM27 force field was used to account for the dynamical nature of the primary solvation shell of 45 DMF molecules. In the calculations using solvent shells, the required atomic charges for each solvent molecule were obtained from ROB3LYP/6-31G(d) calculation on a single isolated DMF molecule. Randomly sampled configurations generated by the Monte Carlo (MC) technique were used to determine the contribution of the primary shell to the free energy of solvation of the three species. The dynamical nature of the primary shell significantly corrects the free energy of solvation. Frequency calculations at the ROB3LYP/6-31G(d) level were carried out on the optimized geometries of truncated 47-atom models of the neutral and ionic species in vacuum so as to determine the differences in thermal energy and molecular entropy. The Born energy of ion-dielectric interaction, the Onsager energy of dipole-dielectric interaction, and the Debye-Hu¨ckel energy of ion-ionic cloud interaction for the pheophytin-solvent aggregate were added as perturbative corrections. The Born interaction also makes a large contribution to the absolute free energy of reduction. An implicit solvation model (DPCM) was also employed for the calculation of standard reduction potentials in DMF. Both the models were successful in reproducing the standard reduction potentials. An explicit solvent treatment(QM/MM/MC + Born + Onsager + Debye corrections) yielded the one electron reduction potential of Pheo-a as -0.92 ( 0.27 V and the two electron reduction potential as -1.34 ( 0.25 V at 298.15 K, while the implicit solvent treatment yielded the corresponding values as -1.03 ( 0.17 and -1.30 ( 0.17 V, respectively. The calculated values more or less agree with the experimental midpoint potentials of -0.90 and -1.25 V, respectively. Moreover, a numerical finite difference Poisson-Boltzmann solver (FDPB) along with the DPCM methodology was employed to obtain the reduction potential of pheophytin in the thylakoid membrane. The calculated reduction potential value of -0.58 V is in excellent agreement with the reported value -0.61 V.

I. Introduction Recent years have seen a surge of activity in the theoretical investigation of the electron transport route for the light phase of photosynthesis.1-4 The outcome of these investigations mainly confirms the redox properties of the involved species. The results also generate a physical picture of the environments, that is, the organization of the protein chains surrounding the species within the thylakoid membrane and of the associated solvent molecules in solution. There have been other advances as well. Berry and Rumberg5 calculated the rate of electron transport and trans-membrane proton translocation, the protonelectron stoichiometry, and the stationary trans-membrane pH difference by carrying out numerical simulation of a model of the photosynthetic electron transport chain under the steadystate condition. Newton et al. investigated electron-transfer reactions in condensed phases.6 Parson and Warshel prepared a density matrix model of photosynthetic electron transfer with * E-mail: [email protected].

microscopically based estimates of vibrational relaxation times.7 Our activity in this field has been based on two aspects. First, we calculated the standard redox potentials of some of the biomolecules involved in the Z-scheme of photosynthesis.8 Second, the computed total energies and molecular orbitals can be used to calculate the rate of electron transfer from one redox species to another in a semirigid condensed phase like the thylakoid membrane region.9 Pheophytin-a (Pheo-a) is one of the biomolecules involved in the electron-transfer process in the light phase of photosynthesis. In this work, we investigate the electronic structure of pheophytin in solvent N,N-dimethyl formamide (DMF) and determine its redox properties. The solvent DMF was chosen as it is one of the most common solvents in which polarographic and voltammetric experiments are conducted. Studies of Photosystem-II (PS II) reaction center, P680, by Klimov et al. have shown that pheophytin accepts an electron from the primary electron donor in PS II and donates an electron to a quinone.10 Experimental investigations using optical, Fourier transform

10.1021/jp067383t CCC: $37.00 © 2007 American Chemical Society Published on Web 05/31/2007

Standard Reduction Potentials of Pheophytin-a infrared (FTIR), and magnetic resonance methods11a-d have confirmed the presence of a charge-separated form. A number of recent studies suggest that the charge pair is Chl+-Pheo-awhere Chl is the true primary electron donor in P680.12 Therefore, the redox potentials of Pheo-a are of sufficient experimental and theoretical interest. Kazakova et al. studied the products of Pheo-a reduction and identified their reduction potentials from sampled dc-polarograms of these species under aprotic conditions in DMF.13 The experimental determination of redox potential is normally undertaken by using the calomel electrode as reference, and the observed data are converted to a scale where the normal hydrogen electrode (NHE) potential is zero. A study of the changes in the cyclic voltammetry curves at various pH values identified a number of products of pheophytin reduction (standard reduction potentials are given in brackets): Pheo‚(-0.90 V), Pheo2- (-1.25 V), PheoH- (-0.72 V), and PheoH2 (+0.05 V).13 The standard reduction potential of Pheo-a within the membrane is -0.61 V which is smaller in magnitude than -0.90 V, the standard potential of the species in DMF. This is so in spite of the fact that there is a hydrophobic polypeptide network around Pheo-a in PS II.14 The calculation of solvation energy becomes very important for the theoretical determination of the absolute redox potentials. If the total energy of the solvated species varies by 0.01 hartree, the calculated redox potential deviates by 0.27 V. We have employed both explicit and implicit solvation models in our calculations. An explicit treatment of a number of solvent molecules using the ONIOM (QM/MM) methodology followed by a dynamical consideration through Monte Carlo (MC) simulations of the solvent molecules arranged around the moiety is expected to provide a greater insight into the role played by the solvents in determining the redox potentials. It is also possible to treat the solvent effect by a polarized continuum model (PCM) approach.15 These aspects are made clear in the latter sections. As a part of our ongoing investigation on the redox processes occurring in the Z-scheme,5,6 we have theoretically determined the absolute reduction free energy of Pheo-a in DMF and in the membrane. The present work is on the theoretical determination of the one- and two-electron reduction potentials of Pheo-a in DMF by using (1) the QM/MM methodology followed by MC which is a “discrete” or “cluster” approach that maintains the molecular nature of the solvent and (2) by using the density functional theory-implicit solvation model (DFT-DPCM) method which is a “continuum” approach involving a macroscopic description of the solvent in terms of the dielectric constant. This work also theoretically determines the reduction potential of Pheo-a in the membrane by using the finite difference Poisson-Boltzmann solver (FDPB) along with the DFT-DPCM method. The electrostatic and nonelectrostatic contributions to the free energy change of the reduction reaction were obtained by the FDPB and DPCM methods, respectively. The molecular structure selected for our quantum chemical calculations is discussed in section II. A short description on the method of calculation is given in section III. Section IV presents the computed results and a summary of conclusions. II. Molecular Geometry The structure of Pheo-a is shown in Figure 1a. There are 141 atoms in this molecule. For computational purposes, we work on the truncated forms Pheo-47 and Pheo-86, species with 47 and 86 atoms respectively, as shown in Figure 1b and c. Pheo-47 has the vinyl group intact and all the other side chains

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7211

Figure 1. Structure of Pheophytin-a: (a) the whole molecule, (b) the truncated model Pheo-47 with all substituents except the vinyl group replaced by hydrogen (gas-phase optimized geometry), and (c) ethyl derivative Pheo-86 with all the side chains retained except for the phytyl chain that is replaced by an ethyl group (optimized geometry in the gas phase).

replaced by hydrogen atoms. Pheo-86 retains all the side chains intact except that the phytyl chain is replaced by an ethyl group. Quantum chemical calculations to determine the electronic energy values are carried out on Pheo-86, its anion and its dianion. Atomic coordinates are optimized for the isolated species. To obtain the solvated geometries, atomic coordinates are reoptimized by using a two-layer ONIOM technique16 with mechanical embedding. Pheo-86 constitutes the higher-level DFT layer while 45 DMF molecules are considered to form the lower-level layer as shown in Figure 2. The number of solvent molecules has been actually varied to study the corresponding effect, and the final calculations have been done with 45 molecules. The dynamical nature of the primary solvation shell is considered by adopting Monte Carlo (MC) simulation of the solvent layer consisting of the DMF molecules. A typical snapshot of the configuration of all molecules involved is exhibited in Figure 3. For the PCM calculations in DMF, we use the vacuumoptimized geometries. For the FDPB method, we use the crystal structure from the pdb file (2AXT entry in the Brookhaven Protein Data Bank) at 3.00 Å of resolution. The pdb geometry of Pheo-a with its environment at a radial distance of 4 Å from it is shown in Figure 4. The computation of thermal energies of the large number of atoms involved in the two ONIOM layers is a formidable task.

7212 J. Phys. Chem. B, Vol. 111, No. 25, 2007

Figure 2. System chosen for ONIOM geometry optimization. The DFT layer is in the wire form and is lightly shaded. The lower layer is in the tube form and is dark colored. The results of optimization are used for the calculation of the electronic energy of the neutral species and its anions.

Mehta and Datta

Figure 4. Section of the overall structure of the PSII model (2AXT entry in the Brookhaven Protein Data Bank). This figure shows all the proteins and cofactors at a radial distance of 4 Å from Pheo-a within the thylakoid membrane. The truncated model consists of the following residues: (a) 4 units of LEU, 4 ALA, 1 THR, 2 ILE, 3 PHE, 2 TYR, 1 GLN, 2 PRO, 1 MET, 1 GLY, 1 VAL, and 1 TRP (shown as narrow sticks), (b) 2 molecules of CLA (chlorophyll-a) and 1 MGE [(1S)-2(alpha-L-allopyranosyloxy)-1-[tridecanoyloxy)methyl] ethyl palmitate] (shown in the ball and stick model), and (c) 1 PHO (pheophytin-a) (shown as broad lines).

side chains do not change significantly on ionization. In fact, the thermal energy difference between Pheo-a and its anion and dianion is only a minor contributor to the free energy of reduction of the species. III. Method of Calculation

Figure 3. Snapshot of ethyl derivative of Pheo-a with the primary solvation shell of 45 DMF molecules. The primary solvation shell is treated by MC simulation and is shown in the ball and stick model. For the sake of clarity, hydrogen atoms are not shown here.

What we need is the difference in thermal energy of Pheo-a from its anion and dianion. Therefore, we separately optimize the geometry of bare Pheo-47, its anion and dianion, at the ROB3LYP/6-31G(d) level. Frequency calculations at the same level were performed on the newly optimized geometries. The thermal energy and the molecular entropy were obtained from a harmonic analysis of the vibrational frequencies, using an ideal gas approximation. This technique is expected to give rise to a reliable difference of thermal energies, because the additional charge in the anion and dianion is localized largely within the π frame, and consequently, the bond force constants for the

The calculation of the absolute redox free energy changes of biomolecules that are generally large species is extremely difficult. Different methods have been employed in the past to study the redox behavior of biomolecules. Noodleman et al. studied the redox potentials of Fe-S clusters in both protein and solvent media by employing DFT along with the broken symmetry (BS) technique to perform electronic structure calculations.2a,b Olsson et al.3 carried out investigations of plastocyanin and rusticyanin by treating the active site at the DFT level. The electron densities over the protein backbone were considered to be frozen relative to the nuclear coordinates, in order to incorporate the effect of the backbone atoms. They employed the semi-macroscopic PDLD/S-LRA method and classical all-atom simulations with and without a polarizable force field. Their calculations reproduced the difference between the reduction potentials of plastocyanin and rusticyanin. Corni et al.4a investigated the details of electron transfer in redox active azurin molecules. These authors performed plane wave pseudopotential density functional theory calculations of the protein active site in the two possible oxidation states of Cu(I) and Cu(II). Friesner et al.4b employed the self-consistent reaction field (SCRF) calculations at the Hartree-Fock level using the ab initio electronic structure code PSGVB [6-31G(d,p)] together

Standard Reduction Potentials of Pheophytin-a

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7213

with an accurate numerical Poisson-Boltzmann solver to obtain the relative reduction potential of bacteriopheophytin (BPh). The computed reduction potential was found to be equal to -2.655 V in DMF whereas the experimental value is around -4.00 V.4c,d Ryde et al.4e adopted the DFT-COSMOPCM methodology to obtain the reduction potential values of pheophytin as -2.58, -1.53 and -0.99 V at dielectric constants of 1, 4, and 80, respectively. In a previous work, we calculated the absolute redox free energy and the reduction potential of plastocyanin at pH 7.0, and estimated the entropy of reduction at pH 3.8.8g We also determined the absolute free energies and standard oxidation and reduction potentials of Chl-a in acetonitrile.8h An ONIOM methodology was adopted to optimize the solvated active site, with the active site being treated at the UB3LYP/6-31G level and the nearby solvents treated by the UFF force field. In the present work, we investigate the ethyl derivative of pheophytin-a, and the basic task is to compute the difference of electronic and thermal energies of the anionic and dianionic forms of Pheo-a from that of its neutral form in the presence of a solvent or any other polarizable medium. The calculated free energies of these species are directly used to evaluate the absolute reduction free energy change of Pheo-a in DMF and in the membrane. Geometry Optimization. The overall procedure for the optimization of molecular geometry requires a consideration of not only the solute but also the contribution of the solvent molecules around it, either by the PCM treatment or by the QM/ MM/MC method. The use of an ONIOM QM/MM methodology to optimize the pheophytin geometry also bears a great similarity to the method used later to determine energies of the associated solvent molecules by the MC method. In ONIOM, the total energy of the system is given by

EONIOM,2layer ) Ereal,MM + Emodel,QM - Emodel,MM

(1)

where the real system contains all the atoms and is calculated only at the MM level. The model system contains the part that is treated at the QM level. Both QM and MM calculations are carried out for the model system. The solvent molecules are modeled as the lower layer, while DFT is used on pheophytin that forms the higher layer. The lower layer is treated by the universal force field (UFF) methodology.17 The UFF uses a charge equilibration scheme (QEq), which enables charges to be generated dynamically in response to the environment and has been shown to reproduce atomic charges accurately for a variety of systems. For the geometry optimization, the computations on the higher layer are performed using the ROB3LYP/ 6-31G(d) methodology. The Gaussian 03 (G03) suite of programs18 is used in all the calculations. Both of the methods (QM/MM/MC and DFT-DPCM) have successfully reproduced the standard reduction potential values. Electronic and Thermal Energies. The optimized structures of Pheo-a are obtained by the ONIOM method at the ROB3LYP level employing the 6-31G(d) basis set while the solvent molecules are treated by the UFF method. Electronic energies of the neutral, anionic, and dianionic forms are finally obtained by carrying out single point calculations of the solvated species at the ROB3LYP/6-311+G(2d,2p) level. As mentined earlier, molecular characteristics such as thermal energy and entropy at the 6-31G(d) level are obtained for the optimized truncated 47-atom model of the neutral, anionic, and dianionic forms of Pheo-a. All thermal energy corrections include zero point vibrational corrections to energy.

Solvation Energy. The solvation free energy is especially difficult to estimate to the accuracy required for the determination of redox potentials. Each technique to represent the solvation presents its own peculiar weakness. For example, not all of the electronic aspects of hydrogen bonding between solvent and solute are correctly described by continuum models, whereas explicit descriptions, which are always size-limited, cannot fully take into account long-range (or bulk) effects.15 In the explicit solvation model, electronic energies of the neutral, anionic, and dianionic species under the influence of 45 DMF molecules are obtained by the QM/MM technique so that the computed energy includes all forms of electrostatic contributions of the surrounding solvent molecules. Furthermore, the primary solvation layer is treated to be dynamical, and its energy obtained from MC simulations accounts for not only corrections to solvation energy but also entropic effects. A constant temperature (298.15 K), constant volume (cubic unit cell of volume 25 × 25 × 25 Å3) (NVT) ensemble is adopted. Configurations of 45 DMF molecules around Pheo, Pheo-, and Pheo2- were initially optimized under the molecular mechanics (CHARMM27) force field by employing the steepest descent algorithm. The optimized configurations were allowed to relax at a regular max delta (maximum trial atomic displacement) of 0.05 Å and for a total of 1000 run steps to simulate an adequate description of the MC solvent dynamics. Each run step involved arbitrary translational and rotational movement of all the solvent molecules in separate microsteps. To mimic the infinite system and to reduce the finite boundary effects, periodic boundary conditions are employed. One and two Na+ ions are introduced as counterions near pheophytin at the center of the cubic box to neutralize the charges of the anion and dianion, respectively. The geometry of each Pheophytin moiety was kept fixed at the ONIOM optimized structure. MC simulations involved movement of only the DMF molecules in the presence of fixed Pheophytin and counterions. No attempt was made to ensure a reasonable pressure. The system was equilibrated for the first 300 run steps. Thereafter, data were collected after each run step. The remaining 700 data points were averaged by Boltzmann distribution. The standard deviation of the simulated energies was calculated for neutral, anion, and dianion species. Along with the electrostatic effects, other interactions such as dispersive interactions are considered through the CHARMM27 force field. Mulliken charges generated by ROB3LYP/6-31G(d) calculations for neutral, anionic, and dianionic forms of Pheo (not shown here) and for DMF (O -0.4682, C 0.3516, H 0.1120, N -0.3883, C -0.2972, C -0.3119, H 0.1673, H 0.1674, H 0.1674, H 0.1478, H 0.1759, H 0.1759) are used. We found simulations employing ESP-derived charges (CHelpG and Merz Kollman charges) to be unstable. Another set of simulations of 45 DMF molecules was performed in the presence of only the fixed counterions. Subtracting the average energy of the solvent layer in the presence of only the fixed counterions from that in the presence of both fixed Pheophytin and counterions gives the average energy of the solvent layer, which includes the free energy of interaction of Pheophytin with the primary solvation shell. Different positions of the counterion give nearly the same energy difference. Hyperchem Professional release 7 for Windows was used in these calculations.20 The reason for the use of a different force field in ONIOM and MC is that UFF is not available in Hyperchem and furthermore CHARMM27 has been designed to work well on biomolecules. DMF does not have a hydrogen-bonded network. Consequently, the entropic effects arising from a structure breakdown are limited to the region near the solute. Therefore, the remaining

7214 J. Phys. Chem. B, Vol. 111, No. 25, 2007

Mehta and Datta

solvent is treated as a structureless dielectric continuum. The Born free energy of the ion-dielectric interaction is given by

GBorn ) -

Q2 1 12a0 

(

)

(2)

where a0 is the radius and Q is the total charge of the soluteprimary layer complex. The quantity  is the dielectric constant of the bulk solvent. The Onsager free energy term is given by

GOnsager ) -

( - 1) µ2 (2 + 1) a03

(3)

where µ is the ground state dipole moment of the solute-layer complex. We have used a value of 38.7 for  of DMF in our work. Gaussian 03 was used to obtain the values of a0 and µ of the neutral, anionic, and dianionic forms of solvated Pheo-86 at the ROB3LYP/6-31G(d) level (shown in the Supporting Information, Table S5). The computed a0 values are a little less than half of the side length of the box considered for periodic boundary conditions in the MC simulations. The additional stabilization provided by the medium through the ion-ionic atmosphere interaction is estimated by using Debye-Hu¨ckel theory. We write

EDH ) -

Q2κ 2

(4)

The Debye-Hu¨ckel reciprocal length is given by

κ)

[ ∑( ) ] 4πe02 kBT

1/2

Nci

i

1000

zi

2

(5)

where e0 is the electronic charge, kB is the Boltzmann constant, T is the temperature of the system in Kelvin, N is the Avogadro number, ci is the concentration of the ith species in molar units, and zie0 is the charge of the ith ion.21 We take 1 M concentrations for each type of ion (ionic strength: 1 for anionic solution and 4 for dianionic solution). The standard Gibbs free energy G° of a species is given by

G° ) Eel(ONIOM) + Ethermal + (GMC + GBorn + GOnsager + EDH) + PV - TS (6) where Eel(ONIOM) ) (E)model,QM. Thus, for a specific process, ∆G° per molecule is written as

∆G° ) ∆(E)model,QM + ∆Ethermal + (∆GMC + ∆GBorn + ∆GOnsager + ∆EDH) + kBT∆ng - T∆S (7) where ∆ng is the change in the number of moles in the gas phase. The PV term is negligibly small for the solvated species. Hence, for a specific process in a solvent medium,

∆G° ) ∆(E)model,QM + ∆Ethermal + (∆GMC + ∆GBorn + ∆GOnsager + ∆EDH) - T∆S (8) We have also employed the DPCM methodology in order to obtain the free energy of reduction of Pheo-a in solution and within the membrane. The FDPB method along with DPCM solvation model at the ROB3LYP/6-31G(d) level were employed to obtain the total free energies of Pheo-a and Pheo-awithin the thylakoid membrane. The dielectric constant of the thylakoid membrane (in) is taken as 4.0, the value that was

used with success earlier.22 The dielectric constant of the external medium (out) is taken as 80.0.23 The electrostatic and nonelectrostatic contributions to the free energy change during the reduction process are obtained by employing the FDPB and DPCM methods, respectively. We then obtain

∆G°red ) ∆Gel + ∆Gnon-el

(9)

Standard Reduction Potential. The standard reduction potential is written as

E°red (in V) ) -

1 [∆G°red (in eV) - ∆G°H+/(1/2)H2(g) (in eV)] ne (10)

where ne is the number of electrons transferred in the half-cell reaction and the ∆G° values are quantities per molecule. We used ∆G°[H+(aq)f(1/2)H2(g)] ) -4.36 eV.19 IV. Computed Results and Discussion We used Gaussian 03 software and an Intel Xeon 3.2 GHz server with 4 GB of virtual memory space for each calculation. The energy of each isolated species with ROB3LYP/6-31G(d) optimized geometry is given in the Supporting information (Table S1). To decide the number of solvent molecules to be used for the ONIOM calculations, we first employed the 3-21G basis set and performed the ONIOM calculation on the neutral species. The molecular geometry of neutral Pheophytin-a in 25, 50, and 100 DMF molecules was optimized by employing the ONIOM (ROB3LYP/3-21G:UFF) methodology with mechanical embedding. For each total number of solvent molecules, the position of the solvent molecules in the ONIOM-optimized aggregate was slightly changed to obtain a new solvent conformation. The new conformation was then taken as the starting point for another ONIOM optimization, and the process was repeated several times. The energy of the active site [E2 ) (E)model,QM] was noted (Supporting Information, Table S2). Three conformers were considered for 25 DMF molecules along with Pheo-a, six conformers were investigated for 50 DMF molecules, and one conformer was considered for 100 DMF molecules. The energy difference (E2)ONIOM - Eiso oscillated widely when 25 DMF molecules are considered. In the case of 50 DMF molecules, the energy difference varies slightly around 0.0028 au. The 100 DMF calculation leads to more or less the same energy difference. Therefore, we settled for 45 DMF molecules in all subsequent calculations. We also checked the effect of the variation of basis sets while using 45 DMF molecules. The calculated energy differences in atomic units are as follows: 0.0016, 0.0023 (STO-3G); 0.0026, 0.0029 (3-21G); 0.0078, 0.0075 (6-31G(d)) (Table S3). The difference increases with the basis size. Therefore, the final calculations were done with the fairly large 6-311+G(2d,2p) basis set. DMF causes a slight puckering of the protoporphyrin ring, thereby forcing it to adopt a nonplanar geometry as shown in Figure 5. The final ONIOM 6-31G(d) optimizations took 96, 101, and 102 CPU h for neutral, anionic, and dianionic species. The single point calculations with 6-311+G(2d,2p) bases needed 78, 80, and 88 CPU hours. The corresponding energies are given in the Supporting Information (Table S1). The ROB3LYP/6-31G(d) thermal energy calculations took 25, 42, and 59 CPU h for neutral, anionic, and dianionic species with 47 atoms. The results are given in the Supporting Information (Table S4). The volume (V) and the dipole moment (µ) of the solvated species, that is, the ONIOM optimized Pheophytin molecule with 45 DMF

Standard Reduction Potentials of Pheophytin-a

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7215

Figure 5. Optimized geometries of (a) Pheo-a in vacuum, (b) Pheo-a, (c) Pheo-a-, and (d) Pheo-a2- in DMF. The ONIOM optimized geometries b, c, and d display nonplanarity and slight puckering of the protoporphyrin ring of Pheo in DMF.

molecules (Figure 2), are obtained by employing the ROB3LYP/ 6-31G(d) methodology. Ideally, one should calculate the volume and dipole moment of the solvated species at different steps of MC and, then, use the average values. However, this exercise would take an enormous amount of CPU time. We took the ONIOM optimized geometry (the initial geometry of MC) and the final geometry of MC. These two geometries gave quite similar volumes at the ROB3LYP/6-31G(d) level (3526.945 and 3532.665 Å3 for the neutral species) and similar dipole moments (8.924 and 9.043 D). The Born and Onsager corrections for the initial geometry of MC and the geometry after the final step of MC would differ by an amount on the order of 10-5 and 10-6 au only. Therefore, we used the volume and dipole moment of the ONIOM optimized geometry. These are given in the Supporting Information (Table S5). The effective radius (a0) is then obtained by employing the formula V ) (4/3)πa03. The MC calculations required 48-72 CPU h each. Calculated energy values are given in the Supporting Information (Table S6). The PCM calculations took 89, 102, and 132 CPU h for neutral Pheo86, its anion, and dianion, respectively. The EDPCM values are given for the dielectric constant 38.7 in Table S7. The medium interaction energy plays a critical role in deciding the redox potentials. Electronic energy differences are

found to be of the same order as the solvation energy differences. Thermal energy differences are much smaller in comparison, and the use of the in-vacuo thermal energies is more or less valid. We take two different radius values, a0 and (a0 + a0DMF)/ 2, for the calculation of uncertainties associated with both the Born and Onsager energy terms as the normal practice is to add half of the radius of the solvent. The average of the energies corresponding to radius a0 and (a0 + a0DMF)/2 is taken in each case. The deviation of the energy with radius a0 from the average energy is taken here as the error spread. The calculation of the total solvation energy is shown stepwise in the Supporting Information (Table S6). The dynamical nature of the primary layer significantly corrects the solvation free energy. Thus, the influence of the MC dynamics turns out to be a decisive factor. The Born free energy term also contributes significantly. The reduction reactions considered are as follows:

Pheo(DMF) + e-(g) f Pheo-(DMF)

(11)

Pheo(DMF) + 2e-(g) f Pheo2-(DMF)

(12)

The calculation of ∆H° is given in detail in Table 1. The ∆G°red calculated from eq 8 and the reduction potential calculated from

7216 J. Phys. Chem. B, Vol. 111, No. 25, 2007

Mehta and Datta

TABLE 1: Calculation of ∆H° by the QM/MM/MC Method reactiona Pheo(DMF) + (g) f Pheo‚ (DMF) Pheo(DMF) + 2e-(g) f Pheo2-(DMF) e-

-

∆(E)model,QMb (eV)

∆Ethermalc (eV)

∆Esolvd (eV)

∆H° e (eV)

-1.918 -0.890

-0.079 -0.158

-1.415 ( 0.266 -4.977 ( 0.499

-3.412 ( 0.266 -6.025 ( 0.499

a The energy of the electron is nearly zero and neglected. b (E)model,QM is the energy of the pheophytin moiety at the QM level under the influence of 45 DMF molecules. (E)model,QM is the ONIOM single point energy at the ROB3LYP/6-311+G(2d,2p) level on the ONIOM (ROB3LYP/6-31G(d): UFF) optimized geometry. c For Pheo-47 species at the ROB3LYP/6-31G(d) level. d For Pheo-86 species. The error is calculated as the square root of the sum of squares of the individual standard deviations. e ∆H° ) ∆(E)model,QM + ∆Ethermal + ∆Esolv.

TABLE 2: Calculation of the Absolute Free Energy of One- and Two-Electron Reduction (∆G°red) and E°red of Pheophytin-a in DMF by the QM/MM/MC Methoda midpoint E°red (V) reaction Pheo(DMF) + (g) f Pheo‚ (DMF) Pheo(DMF) + 2e-(g) f Pheo2-(DMF) e-

a

-

∆H° (eV)

-T∆Smol (eV)

∆G°red (eV)

calcd

obsd

-3.412 ( 0.266 -6.025 ( 0.499

-0.027 -0.019

-3.440 ( 0.266 -6.044 ( 0.499

-0.920 ( 0.266 -1.338 ( 0.250

-0.90 -1.25

b

We have used 1 au ) 27.21(16) eV and 1 eV ) 23.06(05) kcal mol-1. b ∆G°red ) ∆H° - T∆Smol.

TABLE 3: Calculation of the Absolute Free Energy of Reduction (∆G°red) and E°red of Pheophytin-aa midpoint E°red (V) reaction Pheo(DMF) + f Pheo(DMF) + 2e-(g) f Pheo2-(DMF) e-(g)

c

Pheo‚-(DMF)

∆EDPCM (eV)

∆Ethermal (eV)

-T∆Smol (eV)

∆G°redc (eV)

calcd

obsd

-3.222 ( -5.935 ( 0.173

-0.079 -0.158

-0.027 -0.019

-3.328 ( 0.173 -6.112 ( 0.173

-1.032 ( 0.173 -1.304 ( 0.173

-0.9 -1.25

0.173b

a We have used 1 au ) 27.2116 eV. b The inherent uncertainty associated in the application of the QM method on large molecules ∼ 4 kcal/mol. ∆G°red ) ∆EDPCM + ∆Ethermal - T∆Smol.

eq 10 are given in Table 2. The calculated value of the oneelectron reduction potential for Pheo (E°red ) -0.920 ( 0.266 V) is in good agreement with the experimentally observed midpoint value -0.90 V. The calculated value for two-electron reduction (E°red ) -1.338 ( 0.250 V) is also more or less in agreement with the observed value -1.25 V. The slight difference between the calculated and observed potential can be ascribed to the uncertainties in the calculational method and the reorganization entropy of the solvent continuum. The DFTDPCM calculation on the vacuum-optimized geometry (Supporting Information, Table S7) gave the one- and two-electron reduction potentials of Pheo as -1.032 ( 0.173 and -1.304 ( 0.173 V at 298.15 K, in agreement with the experimental result (Table 3). Similar results were obtained by employing the DFTCPCM methodology, but the results are not shown here. All the nonelectrostatic terms were included while calculating the total free energy of the Pheophytin moiety in solution. In this paper, we have shown that both explicit and implicit treatment of solvent molecules through QM/MM/MC + Born + Onsager + Debye-Hu¨ckel corrections and DFT-DPCM methodology can lead to the calculated reduction potentials for Pheo-a in DMF in agreement with the experimentally determined values. The difference between the electronic energy of the charged and neutral species is realistically generated only with a large basis set like 6-311+G(2d,2p) and only when some of the post Hartree-Fock effects are incorporated through the ROB3LYP method. The success is also due to a fairly accurate determination of the solvation free energy of the involved species in DMF. Calculation of E°red within the Membrane. The electrostatic contribution to the free energy change of the reduction reaction is obtained with the Delphi program package24 which solves the Poisson-Boltzmann equation25,26 when a distribution of charges and dielectric constants is given. Coordinates. All atomic coordinates were taken from the crystal structure of Pheophytin-a (PDB 2AXT). All cofactors and proteins at a radial distance of 4 Å from Pheo were considered (Figure 4). In the truncated model, the cleaved

TABLE 4: Electrostatic Contribution to the Free Energy of the System (i.e., Pheophytin and its Neighbors which are at a Radial Distance of 4 Å from it) Computed Using the FDPB Solver free energy of electrostatic interaction (Gel)a in eV

a

Pheo-a

Pheo-a-

∆Gel (eV)

-12.080

-15.861

-3.781

Computed using Delphi.

peptide linkages were converted to free carboxylic (COO-) or free amine (NH3+) groups as required. Both polar and nonpolar hydrogen atoms were inserted into the proteins while only polar hydrogen atoms were inserted into the cofactors [(CLA (chlorophyll a), 1 MGE [(1S)-2-(alpha-L-allopyranosyloxy)-1-[tridecanoyloxy)methyl] ethyl palmitate] and PHO (pheophytin-a)] by the WHAT IF software.27 Hydrogen atom positions were energetically optimized by using the MM+ force field. A radius was assigned to all the protein atoms from the Delphi database (CHARMM22 parameter set). The partial charges of the protein atoms were taken from the CHARMM22 parameter set. DFTESP (CHelpG) charges at the ROB3LYP/6-31G(d) level were considered for the cofactors. The grid spacing was set to 0.29 Å (i.e., 3.5 grid points per angstrom), and the truncated model filled 70% of the grid box. Changing the grid size did not affect the result. A probe radius of 1.4 Å, which defines the water-accessible surface, was chosen. Periodic boundary conditions were calculated with the option “coulombic”, and the convergence criterion was chosen to be 0.001 kT/e. Both linear and nonlinear iterations were kept to 1000. The dielectric constant of the truncated model (in) was set to 4.0,22 and the dielectric constant of the solvent (out) was set to 80.0.23 The ∆Gel value is given in Table 4. The nonelectrostatic contribution to the free energy change of the reduction reaction was obtained by employing the DPCM methodology at ROB3LYP/6-31G(d) on the structure of Pheophytin-a extracted from the pdb file (PDB 2AXT), hydrogen atoms inserted, and MM+ optimized. The above calculation gave the same value of Gnon-el (i.e., cavitation + dispersion + repulsion energy) for both Pheo-a and Pheo-a-. The free energy

Standard Reduction Potentials of Pheophytin-a

J. Phys. Chem. B, Vol. 111, No. 25, 2007 7217

TABLE 5: Calculation of the Absolute Free Energy of Reduction (∆G°red) and E°red of Pheophytin-a within the Membranea midpoint E°red (V) ∆Gelb (eV)

∆Gnon-el (eV)

∆G°redc (eV)

calcd

obsd

-3.781

0.0

-3.781

-0.58

-0.61

a

∆Gnon-el was computed employing the DPCM solvation method at the ROB3LYP/6-31G(d) level on the pdb geometry. b Computed taking in ) 4.0 and out ) 80.0. c ∆G°red ) ∆Gel + ∆Gnon-el.

change during the reaction can be ascribed solely due to the change in the electrostatic contribution to the free energy term (Table 5). The reduction potential calculated (E°red) has a value of -0.58 V which is in agreement with the reported value -0.61 V. Discussion. There are mainly three sources of errors in our calculations. The corresponding error spread has been estimated as follows. (1) The Born and Onsager energy corrections have been considered for radii a0 and (a0 + a0DMF)/2. The average energy is taken, and the deviation of the energy corresponding to radius a0 from the mean, as shown in Table S6, is taken as the error spread associated with the Born and Onsager energy correction terms. (2) The standard deviation of all the simulated energy terms corresponds to the uncertainty associated with the EMC term. (3) The inherent uncertainty of applying QM on large molecules about 4 kcal/mol (included in Table 3). Our next objective is to theoretically determine the rate of electron transfer of P680 f Pheo-a in membrane. Meanwhile, work on the determination of the reduction potentials of PheoHand PheoH2 in DMF is in progress, and the results will be communicated elsewhere. Acknowledgment. N.M. gratefully acknowledges financial support from Council of Scientific and Industrial Research. Supporting Information Available: Results of single point 6-311+G(2d,2p) calculations, frequency, and volume at ROB3LYP/6-31G(d) and total solvation energy of the species. Increase in the energy of Pheo-a and its ions in DMF with respect to that in vacuum is provided. The effects of the change in number and orientation of solvent molecules along with a change in the basis set are also supplied. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) (a) Svensson, B.; Etchebest, C.; Tuffery, P.; van Kan, P.; Smith, J.; Styring, S. Biochemistry 1996, 36, 14486. (b) Olsson, M. H. M.; Ryde, U. J. Biol. Inorg. Chem. 1999, 4, 654. (c) O’Malley, P. J.; Collins, S. J. J. Am. Chem. Soc. 2001, 123, 11042. Ryde, U.; Olsson, M. H. M.; Pierloot, K. Theor. Comp. Chem. 2001, 9, 1. (d) Sigfridsson, E.; Olsson, M. H. M.; Ryde, U. J. Phys. Chem. B 2001, 105, 5546. (e) Torres, R. A.; Lovell, T.; Noodleman, L.; Case, D. A. J. Am. Chem. Soc. 2003, 125, 1923. (f) Olsson, M. H. M.; Hong, G.; Warshel, A. J. Am. Chem. Soc. 2003, 125, 5025. (2) (a) Torres, R. A.; Lovell, T.; Noodleman, L.; Case, D. A. J. Am. Chem. Soc. 2003, 125, 1923. (b) Noodleman, L. J. Chem. Phys. 1981, 7, 5737. (3) Olsson, M. H. M.; Hong, G.; Warshel, A. J. Am. Chem. Soc. 2003, 125, 5025.

(4) (a) Corni, S.; Rienzo, F. D.; Felice, D. E.; Molinari, E. Theor. Biochem. Biophys. 2004, 102, 328. (b) Zhang, L. Y.; Friesner, R. A. J. Phys. Chem. 1995, 99, 16479. (c) Fajer, J.; Davis, M. S.; Brune, D. C.; Spaulding, L.; Borg, D.; Forman, A. BrookhaVen Symp. Biol. 1976, 28, 74. (d) Davis, M. S.; Forman, A.; Hanson, L. K.; Fajer, J. J. Phys. Chem. 1979, 83, 3325. (e) Heimdal, J.; Jensen, K P.; Devarajan, A.; Ryde, U. J. Biol. Inorg. Chem. 2007, 12, 49. (5) Berry, S.; Rumberg, B. Bioelectrochemistry 2000, 53, 35. (6) (a) Kornyshev, A.; Newton, M.; Ulstrup, J.; Sanderson, B. Chem. Phys. 2005, 319, 1. (b) Newton, M. D. Annu. ReV. Phys. Chem. 1984, 35, 437. (7) Parson, W. W.; Warshel, A. Chem. Phys. 2004, 296, 201; J. Phys. Chem. B 2004, 108, 10474. (8) (a) Datta, S. N.; Priyadarshy, S. Chem. Phys. Lett. 1990, 173, 360. (b) Mallik, B.; Datta, S. N. Int. J. Quantum Chem. 1994, 52, 629. (c) Datta, S. N.; Deshpande, R. Ind. J. Pure Appl. Phys. 1997, 35, 483. (d) Datta, S. N.; Mallik, B. J. Phys. Chem. B 1997, 101, 4171. Datta, S. N.; Nehra, V.; Jha, A. J. Chem. B 1999, 103, 8768. (e) Datta, S. N.; Prabhakar, B. G. S.; Nehra, V. J. Phys. Chem. B 1999, 103, 2291. (f) Datta, S. N.; Parendekar, P. V.; Lochan, R. C. J. Phys. Chem. B 2001, 105, 1442. (g) Datta, S. N.; Sudhamsu, J.; Pandey, A. J. Phys. Chem. B 2004, 108, 8007. (h) Pandey, A.; Datta, S. N. J. Phys. Chem. B 2005, 109, 9066. (9) Datta, S. N.; Mallik, B. Int. J. Quantum Chem. 1995, 53, 37; 1997, 61, 865. (10) Klimov, V. V.; Krasnovskii, A. A. Biofizika 1982, 27, 179. (11) (a) Forman, A.; Davis., M. S.; Fujita, I.; Hanson, L. K.; Smith, K. M.; Fajer, J. Isr. J. Chem. 1981, 21, 265. (b) Mantele, W. G.; Wollenweber, A. M.; Nabedryk, E.; Breton, J. Proc. Natl. Acad. Sci. U.S.A. 1988, 85, 8468. (c) Lubitz, W.; Issacson, R. A.; Okamura, M. Y.; Abresch, E. C.; Plato, M. Biochim. Biophys. Acta 1989, 977, 227. (d) Naberryk, E.; Andrianambnintsoa, S.; Berger, G.; Leonhard, M.; Mantele, W.; Breton, J. Biochim. Biophys. Acta 1990, 1016, 49. (12) Ernst-Walter K; Hiroshi, I.; Jacek, B.; Bernhard, L.; Wolfram, S. Angew. Chem. Int. Ed. 2006, 45, 1964. (13) Kazakova, A. A.; Kisselev, B. A.; Kozlov, Y. N. Bioelectrochem. Bioenerg. 1989, 21, 367. (14) Scheer, H. Chlorophylls; CRC: London, 1990. (15) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999. (16) Vreven, T.; Byun, K. S.; Koma´romi, I.; Dapprich, S.; Montgomery, J. A.; Morokuma, K.; Frisch, M. J. J. Chem. Theory Comput. 2006, 2, 815. (17) Rappe´, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A., III; Skiff, W. M. J. Am. Chem. Soc. 1992, 114, 10024. Rappe´, A. K.; Goddard, W. A., III. J. Phys. Chem. 1991, 95, 3358. (18) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A.; Gaussian 03, revision C. 02; Gaussian, Inc.: Wallingford, CT, 2004. (19) Lewis, A.; Bumpus, J. A.; Truhlar, D. G.; Cramer, C. J. J. Chem. Ed. 2004, 81, 596. (20) HyperChem Professional, release 7.01 for Windows; Hypercube Inc.: Gainesville, FL, 2002. (21) Bockris, J. O’M.; Reddy, A. K. N. Modern Electrochemistry; Plenum Press: New York and London, 1998. (22) Ishikita, H.; Knapp, E. FEBS Lett. 2005, 579, 3190. (23) Ishikita, H.; Knapp, E. J. Am. Chem. Soc. 2005, 127, 1963. (24) Delphi V.4, release 1.0; Accelrys Inc.: Pisa, Italy, 2001. (25) Rocchia, W.; Alexov, E.; Honig, B. J. Phys. Chem. B 2001, 105, 6507. (26) Rocchia, W.; Sridharan, S.; Nicholls, A, Alexov, E.; Chiabrera, A; Honig, B. J. Comp. Chem. 2002, 23, 128. (27) Vriend, G. J. Mol. Graph. 1990, 8, 52.