Gaussian fluctuations and linear response in an electron transfer protein

Potential distribution theory of alchemical transfer

We present an analytical description of the Alchemical Transfer Method (ATM) for molecular binding using the Potential Distribution Theory (PDT) formalism. ATM models the binding free energy by mapping the bound and unbound states of the complex by translating the ligand coordinates. PDT relates the free energy and the probability densities of the perturbation energy along the alchemical path to the probability density at the initial state, which is the unbound state of the complex in the case of a binding process. Hence, the ATM probability density of the transfer energy at the unbound state is first related by a convolution operation of the probability densities for coupling the ligand to the solvent and coupling it to the solvated receptor-for which analytical descriptions are available-with parameters obtained from maximum likelihood analysis of data from double-decoupling alchemical calculations. PDT is then used to extend this analytical description along the alchemical transfer pathway. We tested the theory on the alchemical binding of five guests to the tetramethyl octa-acid host from the SAMPL8 benchmark set. In each case, the probability densities of the perturbation energy for transfer along the alchemical transfer pathway obtained from numerical calculations match those predicted from the theory and double-decoupling simulations. The work provides a solid theoretical foundation for alchemical transfer, offers physical insights on the form of the probability densities observed in alchemical transfer calculations, and confirms the conceptual and numerical equivalence between the alchemical transfer and double-decoupling processes.

Explicit Projection of Stokes Shifts onto Solvent Motion in an Aqueous Liquid and Linear Response Theory

We investigate the molecular origin of the fluorescence Stokes shift in an aqueous liquid. By examining the speed of energy change, the solvation response function is explicitly projected onto the translational and rotational motions of water molecules for both nonequilibrium relaxation and equilibrium fluctuations. Molecular dynamics simulations of a tryptophan solution show that these two processes have highly consistent dynamics, not only for the total response function but also for the decomposed components in terms of specific molecular movements. We found that the rotational mode governs the relaxation of the Stokes shift, whereas the translational mode contributes non-negligibly with slower dynamics. This consistency implies the similarity of the underlying translational and rotational movements of water molecules as the system is far away from and at equilibrium, supporting the validity of the linear response theory at the molecular level. The decomposition methodology is also applicable to a rigid solvent.

Electronic Polarization Is Essential for the Stabilization and Dynamics of Buried Ion Pairs in Staphylococcal Nuclease Mutants

Buried charged residues play important roles in the modulation of protein stabilities and conformational dynamics and make crucial contributions to protein functions. Considering the generally nonpolar nature of protein interior, a key question concerns the contribution of electronic polarization to the stabilization and properties of buried charges. We answer this question by conducting free energy simulations using the latest polarizable CHARMM force field based on Drude oscillators for a series of Staphylococcal nuclease mutants that involve a buried Glu-Lys pair in different titration states and orientations. While a nonpolarizable model suggests that the ionized form of the buried Glu-Lys pair is more than 40 kcal/mol less stable than the charge-neutral form, the two titration states are comparable in stability when electronic polarization is included explicitly, a result better reconcilable with available experimental data. Analysis of free energy components suggests that additional stabilization of the ionized Glu-Lys pair has contributions from both the enhanced salt-bridge strength and stronger interaction between the ion-pair and surrounding protein residues and penetrated water. Despite the stronger direct interaction between Glu and Lys, the ion-pair exhibits considerably larger and faster structural fluctuations when polarization is included, due to compensation of interactions in the cavity. Collectively, observations from this work provide compelling evidence that electronic polarization is essential to the stability, hydration, dynamics, and therefore function of buried charges in proteins. Therefore, our study advocates for the explicit consideration of electronic polarization for mechanistic and engineering studies that implicate buried charged residues, such as enzymes and ion transporters.

The redox potential of a heme cofactor in Nitrosomonas europaea cytochrome c peroxidase: a polarizable QM/MM study

Redox reactions are crucial to biological processes that protect organisms against oxidative stress. Metalloenzymes, such as peroxidases which reduce excess reactive oxygen species into water, play a key role in detoxification mechanisms. Here we present the results of a polarizable QM/MM study of the reduction potential of the electron transfer heme in the cytochrome c peroxidase of Nitrosomonas europaea. We have found that environment polarization does not substantially affect the computed value of the redox potential. Particular attention has been given to analyzing the role of electrostatic interactions within the protein environment and the solvent on tuning the redox potential of the heme co-factor. We have found that the electrostatic interactions predominantly explain the fluctuations of the vertical ionization/attachment energies of the heme for the sampled configurations, and that the long range electrostatic interactions (up to 40 Å) contribute substantially to the absolute values of the vertical energy gaps.

Reverse Protonation of Buried Ion-Pairs in Staphylococcal Nuclease Mutants

Although buried titratable residues in protein cavities are often of major functional importance, it is generally challenging to understand their properties such as the ionization state and factors of stabilization based on experimental studies alone. A specific set of examples involve buried Glu-Lys pairs in a series of variants of Staphylococcal nuclease, for which recent structural and thermodynamic studies appeared to suggest that both the stability and the ionization state of the buried Glu-Lys pair are sensitive to its orientation (i.e., Glu23-Lys36 vs Lys23-Glu36). To further clarify the situation, especially ionization states of the buried Glu-Lys pairs, we have conducted extensive molecular dynamics simulations and free energy computations. Microsecond molecular dynamics simulations show that the hydration level of the cavity depends on the orientation of the buried ion-pair therein as well as its ionization state; free energy simulations recapitulate the relative stability of Glu23-Lys36 (EK) vs Lys23-Glu36 (KE) mutants measured experimentally, although the difference is similar in magnitude regardless of the ionization state of the Glu-Lys pair. A complementary set of free energy simulations strongly suggests that, in contrast to the original suggestion in the experimental analysis, the Glu and Lys residues prefer to adopt their charge-neutral rather than the ionized states. This result is consistent with the low dielectric constant computed for water in the cavity, which makes it difficult for the protein cavity to stabilize a pair of charged Glu-Lys residues, even with water penetration. The current study highlights the role of free energy simulations in understanding the ionization state of buried titratable residues and the relevant energetic contributions, forming the basis for the rational design of buried charge networks in proteins.

Linear-Response and Nonlinear-Response Formulations of the Instantaneous Marcus Theory for Nonequilibrium Photoinduced Charge Transfer

Instantaneous Marcus theory (IMT) offers a way for capturing the time-dependent charge transfer (CT) rate coefficient in nonequilibrium photoinduced CT processes, where the system was photoexcited from its equilibrated ground state vertically to the excitonic state, followed by an electronic transition to the CT state. As derived from the linearized semiclassical nonequilibrium Fermi's golden rule (LSC NE-FGR), the original IMT requires expensive all-atom nonequilibrium molecular dynamics (NEMD) simulations. In this work, we propose computationally efficient linear-response and nonlinear-response formulations for IMT rate calculations, which only require equilibrium molecular dynamics simulations. The linear- and nonlinear-response IMT methods were tested to predict the transient behavior in the photoinduced CT dynamics of the carotenoid-porphyrin-C₆₀ molecular triad solvated in explicit tetrahydrofuran. Our result demonstrated that the nonlinear-response IMT is in excellent agreement with the benchmark NEMD for all cases investigated here, whereas the linear-response IMT predicts the correct trend for all cases but overestimates the transient CT rate in one case involving a significant nonequilibrium relaxation. This mild breakdown of linear-response IMT is due to neglecting the higher-order terms in the exact nonlinear-response IMT. Taking advantage of time translational symmetry, the linear- and nonlinear-response approaches were demonstrated to be able to reduce the computational cost by 80% and 60% compared with NEMD simulations, respectively. Thus, we highly recommend the readily applicable and accurate nonlinear-response IMT approach for simulating nonequilibrium CT processes in complex molecular systems in the condensed phase.

Theoretical Insights on Nonlinear Response Theory of Fluorescence Spectroscopy in Liquids

Correlations between the nonequilibrium solvation dynamics upon the photon excitation of the chromophore and a system's equilibrium fluctuations are deeply studied. As the linear response of the solvent has been linked with Gaussian statistics of the energy fluctuations in the literature, we specifically explore the cases beyond the regime of the linear response theory due to deviation from Gaussian fluctuations. As a continuation of our previous work, an analytical formalism is presented to project the energy shift with various order moments, where the non-Gaussian statistics arise from the overlap of the energy basins on the perturbed potential energy surface. It is shown that the nonequilibrium dynamics still correlate with the spontaneous regressions at equilibrium and are controlled by the decay rates of those higher order components with the prevailing contributions to the energy shift. Molecular dynamics simulations were performed in the protein Staphylococcus nuclease, in which even the dynamics of the high order moments are available. The results further verify the above relationship. Our scheme is used to evaluate Stokes shift using the information on non-Gaussian statistics at equilibrium, thus presenting a broad picture on the correlation between the nonequilibrium process and equilibrium properties in liquids.

A molecular density functional theory approach to electron transfer reactions

Beyond the dielectric continuum description initiated by Marcus theory, the standard theoretical approach to study electron transfer (ET) reactions in solution or at interfaces is to use classical force field or ab initio molecular dynamics simulations. We present here an alternative method based on liquid-state theory, namely molecular density functional theory, which is numerically much more efficient than simulations while still retaining the molecular nature of the solvent. We begin by reformulating molecular ET theory in a density functional language and show how to compute the various observables characterizing ET reactions from an ensemble of density functional minimizations. In particular, we define within that formulation the relevant order parameter of the reaction, the so-called vertical energy gap, and determine the Marcus free energy curves of both reactant and product states along that coordinate. Important thermodynamic quantities such as the reaction free energy and the reorganization free energies follow. We assess the validity of the method by studying the model Cl⁰ → Cl⁺ and Cl⁰ → Cl^- ET reactions in bulk water for which molecular dynamics results are available. The anionic case is found to violate the standard Marcus theory. Finally, we take advantage of the computational efficiency of the method to study the influence of a solid-solvent interface on the ET, by investigating the evolution of the reorganization free energy of the Cl⁰ → Cl⁺ reaction when the atom approaches an atomistically resolved wall.

Termination of Biological Function at Low Temperatures: Glass or Structural Transition?

Energy of life is produced by electron transfer in energy chains of respiration or photosynthesis. A small input of free energy available to biology puts significant restrictions on how much free energy can be lost in each electron-transfer reaction. We advocate the view that breaking ergodicity, leading to violation of the fluctuation-dissipation theorem (FDT), is how proteins achieve high reaction rates without sacrificing the reaction free energy. Here we show that a significant level of nonergodicity, represented by a large extent of the configurational temperature over the kinetic temperature, is maintained in the entire physiological range for the cytochrome c electron transfer protein. The protein returns to the state consistent with the FDT below the crossover temperature close to the temperature of the protein glass transition. This crossover leads to a sharp increase in the activation barrier of electron transfer and is displayed by a kink in the Arrhenius plot for the reaction rate constant.

BAR-based optimum adaptive sampling regime for variance minimization in alchemical transformation: the nonequilibrium stratification

Following the previously proposed equilibrate-state sampling based adaptive sampling regime Optimum Bennett Acceptance Ratio (OBAR), we introduce its nonequilibrium extension, the Optimum Crooks' Equation (OCE) in the current work. The efficiency of the NonEquilibrium Work (NEW) stratification is improved by adaptively manipulating the significance of each nonequilibrium realization followed by importance sampling. As is exhibited in the equilibrium case, the nonequilibrium extension outperforms the simple equal time rule used in nonequilibrium stratification in the sense of minimizing the total variance of the free energy estimate. The speedup of this non-equal time rule is more than 1-fold. The Time Derivative of total Variance (TDV) proposed for the OBAR criterion is extended to determine the importance of each nonequilibrium transformation, which is linearly dependent on the variance. The TDV in the nonequilibrium case gives a totally different importance rank from the standard errors of the free energy differences and OBAR TDV due to the duration of nonequilibrium pulling being added into the OCE equation. The performance of the OCE workflow is demonstrated in the solvation of several small molecules with a series of lambda increments and relaxation times between successive perturbations. To the best of our knowledge, such a nonequilibrium adaptive sampling regime in alchemical transformation is unprecedented.

Extensive Assessment of Various Computational Methods for Aspartate's pKa Shift

A series of computational methods for pK_a shift prediction are extensively tested on a set of benchmark protein systems, aiming at identifying pitfalls and evaluating their performance on high variants. Including 19 ASP residues in 10 protein systems, the benchmark set consists of both residues with highly shifted pK_a values as well as those varying little from the reference value, with an experimental RMS free energy differences of 2.49 kcal/mol with respect to blocked amino acid, namely the RMS pK_a shift being 1.82 pK_a units. The constant pH molecular dynamics (MD), alchemical methods, PROPKA3.1, and multiconformation continuum electrostatics give RMSDs of 1.52, 2.58, 1.37, and 3.52 pK_a units, respectively, on the benchmark set. The empirical scoring method is the most accurate one with extremely low computational cost, and the pH-dependent model is also able to provide accurate results, while the accuracy of MD sampling incorporating alchemical free energy simulation is prohibited by convergence achievement and the performance of conformational search incorporating multiconformation continuum electrostatics is bad. Former research works did not define statistical uncertainty with care and yielded the questionable conclusion that alchemical methods perform well in most benchmarks. In this work the traditional alchemical methods are thoroughly tested for high variants. We also performed the first application of nonequilibrium alchemical methods to the pK_a cases.

Efficient Criterion To Evaluate Linear Response Theory in Optical Transitions

The role of the Gaussian statistics on the solvation dynamics upon the photon excitation of the chromophore is deeply explored. The linear response theory for the fluorescence Stokes shift is investigated. An analytical formulism is presented to recast Stokes shift into the contributions of the equilibrium time correlation functions of the solute-solvent interactions on the excited-state surface, and the latter is further reformed and depicted by the time relaxation of the moment. As the first application of the formulism in the molecular dynamics simulations, it is verified that the efficiency of the linear response theory relies on the Gaussian characteristics of the dominant moments in terms of the Stokes shift, which is identified by the same relaxation dynamics between those moments and the linear order one. The comparisons between the above observations on the linearity of Stokes shift and the explanations in the literature are discussed. The key finding is the development of explicit criterion to measure the appropriateness of applying linear response theory.

BAR-based optimum adaptive sampling regime for variance minimization in alchemical transformation

The efficiency of alchemical free energy simulations with a staging strategy is improved by adaptively manipulating the significance of each ensemble followed by importance sampling.

All-Atom Continuous Constant pH Molecular Dynamics With Particle Mesh Ewald and Titratable Water

Development of a pH stat to properly control solution pH in biomolecular simulations has been a long-standing goal in the community. Toward this goal recent years have witnessed the emergence of the so-called constant pH molecular dynamics methods. However, the accuracy and generality of these methods have been hampered by the use of implicit-solvent models or truncation-based electrostatic schemes. Here we report the implementation of the particle mesh Ewald (PME) scheme into the all-atom continuous constant pH molecular dynamics (CpHMD) method, enabling CpHMD to be performed with a standard MD engine at a fractional added computational cost. We demonstrate the performance using pH replica-exchange CpHMD simulations with titratable water for a stringent test set of proteins, HP36, BBL, HEWL, and SNase. With the sampling time of 10 ns per replica, most pK_a's are converged, yielding the average absolute and root-mean-square deviations of 0.61 and 0.77, respectively, from experiment. Linear regression of the calculated vs experimental pK_a shifts gives a correlation coefficient of 0.79, a slope of 1, and an intercept near 0. Analysis reveals inadequate sampling of structure relaxation accompanying a protonation-state switch as a major source of the remaining errors, which are reduced as simulation prolongs. These data suggest PME-based CpHMD can be used as a general tool for pH-controlled simulations of macromolecular systems in various environments, enabling atomic insights into pH-dependent phenomena involving not only soluble proteins but also transmembrane proteins, nucleic acids, surfactants, and polysaccharides.

Polarizability of the active site of cytochrome c reduces the activation barrier for electron transfer

Enzymes in biology’s energy chains operate with low energy input distributed through multiple electron transfer steps between protein active sites. The general challenge of biological design is how to lower the activation barrier without sacrificing a large negative reaction free energy. We show that this goal is achieved through a large polarizability of the active site. It is polarized by allowing a large number of excited states, which are populated quantum mechanically by electrostatic fluctuations of the protein and hydration water shells. This perspective is achieved by extensive mixed quantum mechanical/molecular dynamics simulations of the half reaction of reduction of cytochrome c. The barrier for electron transfer is consistently lowered by increasing the number of excited states included in the Hamiltonian of the active site diagonalized along the classical trajectory. We suggest that molecular polarizability, in addition to much studied electrostatics of permanent charges, is a key parameter to consider in order to understand how enzymes work.

Equilibrium fluctuation relations for voltage coupling in membrane proteins

A general theoretical framework is developed to account for the effects of an external potential on the energetics of membrane proteins. The framework is based on the free energy relation between two (forward/backward) probability densities, which was recently generalized to non-equilibrium processes, culminating in the work-fluctuation theorem. Starting from the probability densities of the conformational states along the "voltage coupling" reaction coordinate, we investigate several interconnected free energy relations between these two conformational states, considering voltage activation of ion channels. The free energy difference between the two conformational states at zero (depolarization) membrane potential (i.e., known as the chemical component of free energy change in ion channels) is shown to be equivalent to the free energy difference between the two "equilibrium" (resting and activated) conformational states along the one-dimensional voltage couplin reaction coordinate. Furthermore, the requirement that the application of linear response approximation to the free energy functionals of voltage coupling should satisfy the general free energy relations, yields a novel closed-form expression for the gating charge in terms of other basic properties of ion channels. This connection is familiar in statistical mechanics, known as the equilibrium fluctuation-response relation. The theory is illustrated by considering the coupling of a unit charge to the external voltage in the two sites near the surface of membrane, representing the activated and resting states. This is done using a coarse-graining (CG) model of membrane proteins, which includes the membrane, the electrolytes and the electrodes. The CG model yields Marcus-type voltage dependent free energy parabolas for the response of the electrostatic environment (electrolytes etc.) to the transition from the initial to the final configuratinal states, leading to equilibrium free energy difference and free energy barrier that follow the trend of the equilibrium fluctuation relation and the Marcus theory of electron transfer. These energetics also allow for a direct estimation of the voltage dependence of channel activation (Q-V curve), offering a quantitative rationale for a correlation between the voltage dependence parabolas and the Q-V curve, upon site-directed mutagenesis or drug binding. Taken together, by introducing the voltage coupling as the energy gap reaction coordinate, our framework brings new perspectives to the thermodynamic models of voltage activation in voltage-sensitive membrane proteins, offering an a framework for a better understating of the structure-function correlations of voltage gating in ion channels as well as electrogenic phenomena in ion pumps and transporters. Significantly, this formulation also provides a powerful bridge between the CG model of voltage coupling and the conventional macroscopic treatments.

Interplay of Electrostatics and Hydrophobic Effects in the Metamorphic Protein Human Lymphotactin

The human lymphotactin (hLtn) is a protein that features two native states both of which are physiologically relevant: it is a monomer (hLtn10) at 10 °C with 200 mM salt and a dimer (hLtn40) at 40 °C and without salt. Here we focus on the networks of electrostatic and hydrophobic interactions that display substantial changes upon the conversion from hLtn10 to hLtn40 since they are expected to modulate the relative stability of the two folds. In addition to the Arg 23-Arg 43 interaction discussed in previous work, we find several other like-charge pairs that are likely important to the stability of hLtn10. Free energy perturbation calculations are carried out to explicitly evaluate the contribution of the Arg 23-Arg 43 interaction to the hLtn10 stability. hLtn40 features a larger number of salt bridges, and a set of hydrophobic residues undergo major changes in the solvent accessible surface area between hLtn10 and hLtn40, pointing to their importance to the relative stability of the two folds. We also discuss the use of explicit and implicit solvent simulations for characterizing the conformational ensembles under different solution conditions.

Revealing the functional states in the active site of BLUF photoreceptors from electrochromic shift calculations

Photoexcitation with blue light of the flavin chromophore in BLUF photoreceptors induces a switch into a metastable signaling state that is characterized by a red-shifted absorption maximum. The red shift is due to a rearrangement in the hydrogen bond pattern around Gln63 located in the immediate proximity of the isoalloxazine ring system of the chromophore. There is a long-lasting controversy between two structural models, named Q63A and Q63J in the literature, on the local conformation of the residues Gln63 and Tyr21 in the dark state of the photoreceptor. As regards the mechanistic details of the light-activation mechanism, rotation of Gln63 is opposed by tautomerism in the Q63A and Q63J models, respectively. We provide a structure-based simulation of electrochromic shifts of the flavin chromophore in the wild type and in various site-directed mutants. The excellent overall agreement between experimental and computed data allows us to evaluate the two structural models. Compelling evidence is obtained that the Q63A model is incorrect, whereas the Q63J is fully consistent with the present computations. Finally, we confirm independently that a keto-enol tautomerization of the glutamine at position 63, which was proposed as molecular mechanism for the transition between the dark and the light-adapted state, explains the measured 10 to 15 nm red shift in flavin absorption between these two states of the protein. We believe that the accurateness of our results provides evidence that the BLUF photoreceptors absorption is fine-tuned through electrostatic interactions between the chromophore and the protein matrix, and finally that the simplicity of our theoretical model is advantageous as regards easy reproducibility and further extensions.

An Overview of Electrostatic Free Energy Computations for Solutions and Proteins

Free energy simulations for electrostatic and charging processes in complex molecular systems encounter specific difficulties owing to the long-range, 1/r Coulomb interaction. To calculate the solvation free energy of a simple ion, it is essential to take into account the polarization of nearby solvent but also the electrostatic potential drop across the liquid-gas boundary, however distant. The latter does not exist in a simulation model based on periodic boundary conditions because there is no physical boundary to the system. An important consequence is that the reference value of the electrostatic potential is not an ion in a vacuum. Also, in an infinite system, the electrostatic potential felt by a perturbing charge is conditionally convergent and dependent on the choice of computational conventions. Furthermore, with Ewald lattice summation and tinfoil conducting boundary conditions, the charges experience a spurious shift in the potential that depends on the details of the simulation system such as the volume fraction occupied by the solvent. All these issues can be handled with established computational protocols, as reviewed here and illustrated for several small ions and three solvated proteins.

Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects

The calculation of a protein-ligand binding free energy based on molecular dynamics (MD) simulations generally relies on a thermodynamic cycle in which the ligand is alchemically inserted into the system, both in the solvated protein and free in solution. The corresponding ligand-insertion free energies are typically calculated in nanoscale computational boxes simulated under periodic boundary conditions and considering electrostatic interactions defined by a periodic lattice-sum. This is distinct from the ideal bulk situation of a system of macroscopic size simulated under non-periodic boundary conditions with Coulombic electrostatic interactions. This discrepancy results in finite-size effects, which affect primarily the charging component of the insertion free energy, are dependent on the box size, and can be large when the ligand bears a net charge, especially if the protein is charged as well. This article investigates finite-size effects on calculated charging free energies using as a test case the binding of the ligand 2-amino-5-methylthiazole (net charge +1 e) to a mutant form of yeast cytochrome c peroxidase in water. Considering different charge isoforms of the protein (net charges -5, 0, +3, or +9 e), either in the absence or the presence of neutralizing counter-ions, and sizes of the cubic computational box (edges ranging from 7.42 to 11.02 nm), the potentially large magnitude of finite-size effects on the raw charging free energies (up to 17.1 kJ mol(-1)) is demonstrated. Two correction schemes are then proposed to eliminate these effects, a numerical and an analytical one. Both schemes are based on a continuum-electrostatics analysis and require performing Poisson-Boltzmann (PB) calculations on the protein-ligand system. While the numerical scheme requires PB calculations under both non-periodic and periodic boundary conditions, the latter at the box size considered in the MD simulations, the analytical scheme only requires three non-periodic PB calculations for a given system, its dependence on the box size being analytical. The latter scheme also provides insight into the physical origin of the finite-size effects. These two schemes also encompass a correction for discrete solvent effects that persists even in the limit of infinite box sizes. Application of either scheme essentially eliminates the size dependence of the corrected charging free energies (maximal deviation of 1.5 kJ mol(-1)). Because it is simple to apply, the analytical correction scheme offers a general solution to the problem of finite-size effects in free-energy calculations involving charged solutes, as encountered in calculations concerning, e.g., protein-ligand binding, biomolecular association, residue mutation, pKa and redox potential estimation, substrate transformation, solvation, and solvent-solvent partitioning.

The Active Site Loop Modulates the Reorganization Energy of Blue Copper Proteins by Controlling the Dynamic Interplay with Solvent

Understanding the factors governing the rate of electron transfer processes in proteins is crucial not only to a deeper understanding of redox processes in living organisms but also for the design of efficient devices featuring biological molecules. Here, molecular dynamics simulations performed on native azurin and four chimeric cupredoxins allow for the calculation of the reorganization energy and of structure-related quantities that were used to clarify the molecular determinants to the dynamics/function relationship in blue copper proteins. We find that the dynamics of the small, metal-binding loop region controls the outer-sphere reorganization energy not only by determining the exposure of the active site to solvent but also through the modulation of the redox-dependent rearrangement of the whole protein scaffold and of the surrounding water molecules.

Charging free energy calculations using the Generalized Solvent Boundary Potential (GSBP) and periodic boundary condition: a comparative analysis using ion solvation and oxidation free energy in proteins

Free energy simulations using a finite sphere boundary condition rather than a periodic boundary condition (PBC) are attractive in the study of very large biomolecular systems. To understand the quantitative impact of various approximations in such simulations, we compare charging free energies in both solution and protein systems calculated in a linear response framework with the Generalized Solvent Boundary Potential (GSBP) and PBC simulations. For simple ions in solution, we find good agreements between GSBP and PBC charging free energies, once the relevant correction terms are taken into consideration. For PBC simulations with the particle-mesh-Ewald for long-range electrostatics, the contribution (ΔG(P-M)) due to the use of a particle rather than molecule based summation scheme in real space is found to be significant, as pointed out by Hünenberger and co-workers. For GSBP, when the inner region is close to be charge neutral, the key correction is the overpolarization of water molecules at the inner/outer dielectric boundary; the magnitude of the correction (ΔG(s-pol)), however, is relatively small. For charging (oxidation) free energy in proteins, the situation is more complex, although good agreement between GSBP and PBC can still be obtained when care is exercised. The smooth dielectric boundary approximation inherent to GSBP tends to make significant errors when the inner region is featured with a high net charge. However, the error can be corrected with Poisson-Boltzmann calculations using snapshots from GSBP simulations in a straightforward and robust manner. Because of the more complex charge and solvent distributions, the magnitudes of ΔG(P-M) and ΔG(s-pol) in protein simulations appear to be different from those derived for solution simulations, leading to uncertainty in directly comparing absolute charging free energies from PBC and GSBP simulations for protein systems. The relative charging/oxidation free energies, however, are robust. With the linear response approximation, for the specific protein system (CueR) studied, the effect of freezing the protein structure in the outer region is found to be small, unless a very small (8 Å) inner region is used; even in the latter case, the result is substantially improved when the nearby metal binding loop is allowed to respond to metal oxidation. The implications of these results to the applicability of GSBP to complex biomolecules and in ab initio QM/MM simulations are discussed.

Reorganization energy for internal electron transfer in multicopper oxidases

We have calculated the reorganization energy for the intramolecular electron transfer between the reduced type 1 copper site and the peroxy intermediate of the trinuclear cluster in the multicopper oxidase CueO. The calculations are performed at the combined quantum mechanics and molecular mechanics (QM/MM) level, based on molecular dynamics simulations with tailored potentials for the two copper sites. We obtain a reorganization energy of 91-133 kJ/mol, depending on the theoretical treatment. The two Cu sites contribute by 12 and 22 kJ/mol to this energy, whereas the solvent contribution is 34 kJ/mol. The rest comes from the protein, involving small contributions from many residues. We have also estimated the energy difference between the two electron-transfer states and show that the reduction of the peroxy intermediate is exergonic by 43-87 kJ/mol, depending on the theoretical method. Both the solvent and the protein contribute to this energy difference, especially charged residues close to the two Cu sites. We compare these estimates with energies obtained from QM/MM optimizations and QM calculations in a vacuum and discuss differences between the results obtained at various levels of theory.

Simulation of solution phase electron transfer in a compact donor-acceptor dyad

Charge separation (CS) and charge recombination (CR) rates in photosynthetic architectures are difficult to control, yet their ratio can make or break photon-to-current conversion efficiencies. A rational design approach to the enhancement of CS over CR requires a mechanistic understanding of the underlying electron-transfer (ET) process, including the role of the environment. Toward this goal, we introduce a QM/MM protocol for ET simulations and use it to characterize CR in the formanilide-anthraquinone dyad (FAAQ). Our simulations predict fast recombination of the charge-transfer excited state, in agreement with recent experiments. The computed electronic couplings show an electronic state dependence and are weaker in solution than in the gas phase. We explore the role of cis-trans isomerization on the CR kinetics, and we find strong correlation between the vertical energy gaps of the full simulations and a collective solvent polarization coordinate. Our approach relies on constrained density functional theory to obtain accurate diabatic electronic states on the fly for molecular dynamics simulations, while orientational and electronic polarization of the solvent is captured by a polarizable force field based on a Drude oscillator model. The method offers a unified approach to the characterization of driving forces, reorganization energies, electronic couplings, and nonlinear solvent effects in light-harvesting systems.

Microscopic pKa analysis of Glu286 in cytochrome c oxidase (Rhodobacter sphaeroides): toward a calibrated molecular model

As stringent tests for the molecular model and computational protocol, microscopic pK(a) calculations are performed for the key residue, Glu286, in cytochrome c oxidase (CcO) using a combined quantum mechanical/molecular mechanical (QM/MM) potential and a thermodynamic integration protocol. The impact of the number of water molecules in the hydrophobic cavity and protonation state of several key residues (e.g., His334, Cu(B)-bound water, and PRD(a3)) on the computed microscopic pK(a) values of Glu286 has been systematically examined. To help evaluate the systematic errors in the QM/MM-based protocol, microscopic pK(a) calculations have also been carried out for sites in a soluble protein (Asp70 in T4 lysozyme) and a better-characterized membrane protein (Asp85 in bacteriorhodopsin). Overall, the results show a significant degree of internal consistency and reproducibility that support the effectiveness of the computational framework. Although the number of water molecules in the hydrophobic cavity does not greatly influence the computed pK(a) of Glu286, the protonation states of several residues, some of which are rather far away, have more significant impacts. Adopting the standard protonation state for all titratable residues leaves a large net charge on the system and a significantly elevated pK(a) for Glu286, highlighting that any attempt to address the energetics of proton transfers in CcO at a microscopic level should carefully select the protonation state of residues, even those not in the immediate neighborhood of the active site. The calculations indirectly argue against the deprotonation of His334 for the proton pumping process, although further studies that explicitly compute its pK(a) are required for a more conclusive statement. Finally, the deprotonated Glu286 is found to be in a stable water-mediated connection with PRD(a3) for at least several nanoseconds when this presumed pumping site is protonated. This does not support the proposed role of Glu286 as a robust gating valve that prevents proton leakage, although a conclusive statement awaits a more elaborate characterization of the Glu286-PRD(a3) connectivity with free energy simulations and a protonated PRD(a3). The large sets of microscopic simulations performed here have provided useful guidance to the establishment of a meaningful molecular model and effective computational protocol for explicitly analyzing the proton transfer kinetics in CcO, which is required for answering key questions regarding the pumping function of this fascinating and complex system.