Benoı̂t Roux

IMPLICIT SOLVENT MODELS

Implicit solvent models for biomolecular simulations are reviewed and their underlying statistical mechanical basis is discussed. The fundamental quantity that implicit models seek to approximate is the solute potential of mean force, which determines the statistical weight of solute conformations, and which is obtained by averaging over the solvent degrees of freedom. It is possible to express the total free energy as the reversible work performed in two successive steps. First, the solute is inserted in the solvent with zero atomic partial charges; second, the atomic partial charges of the solute are switched from zero to their full values. Consequently, the total solvation free energy corresponds to a sum of non-polar and electrostatic contributions. These two contributions are often approximated by simple geometrical models (such as solvent exposed area models) and by macroscopic continuum electrostatics, respectively. One powerful route is to approximate the average solvent density distribution around the solute, i.e. the solute-solvent density correlation functions, as in statistical mechanical integral equations. Recent progress with semi-analytical approximations makes continuum electrostatics treatments very efficient. Still more efficient are fully empirical, knowledge-based models, whose relation to explicit solvent treatments is not fully resolved, however. Continuum models that treat both solute and solvent as dielectric continua are also discussed, and the relation between the solute fluctuations and its macroscopic dielectric constant(s) clarified.

A simple polarizable model of water based on classical Drude oscillators

A simple polarizable water model is developed and optimized for molecular dynamics simulations of the liquid phase under ambient conditions. The permanent charge distribution of the water molecule is represented by three point charges: two hydrogen sites and one additional M site positioned along the HOH bisector. Electronic induction is represented by introducing a classical charged Drude particle attached to the oxygen by a harmonic spring. The oxygen site carries an equal and opposite charge, and is the center of an intermolecular Lennard-Jones interaction. The HOH gas-phase experimental geometry is maintained rigidly and the dipole of the isolated molecule is 1.85 D, in accord with experiment. The model is simulated by considering the dynamics of an extended Lagrangian in which a small mass is attributed to the Drude particles. It is parametrized to reproduce the salient properties of liquid water under ambient conditions. The optimal model, refered to as SWM4-DP for “simple water model with four site...

Extension to the weighted histogram analysis method: combining umbrella sampling with free energy calculations

The Weighted Histogram Analysis Method (WHAM) of Kumar et al. (J. Comput. Chem. 13 (1992) 1011), is used to combine free energy perturbations with umbrella sampling calculations. The formulation is general and allows optimal calculation of the free energies from a set of molecular dynamics simulations generated in the presence of arbitrary biasing umbrella sampling window potentials. The method yields the free energy associated with a given simulation as well as the probability distribution of the molecular system configurations by extracting the information contained in all the biased simulations (the windows) in an optimal way. The method presents some advantages compared to the standard free energy perturbation (FEP) and thermodynamic integration (TI) methods, because the window potential can be used for restricting the conformational space to specific regions during free energy calculations.

Computations of standard binding free energies with molecular dynamics simulations.

An increasing number of studies have reported computations of the standard (absolute) binding free energy of small ligands to proteins using molecular dynamics (MD) simulations and explicit solvent molecules that are in good agreement with experiments. This encouraging progress suggests that physics-based approaches hold the promise of making important contributions to the process of drug discovery and optimization in the near future. Two types of approaches are principally used to compute binding free energies with MD simulations. The most widely known is the alchemical double decoupling method, in which the interaction of the ligand with its surroundings are progressively switched off. It is also possible to use a potential of mean force (PMF) method, in which the ligand is physically separated from the protein receptor. For both of these computational approaches, restraining potentials may be activated and released during the simulation for sampling efficiently the changes in translational, rotational, and conformational freedom of the ligand and protein upon binding. Because such restraining potentials add bias to the simulations, it is important that their effects be rigorously removed to yield a binding free energy that is properly unbiased with respect to the standard state. A review of recent results is presented, and differences in computational methods are discussed. Examples of computations with T4-lysozyme mutants, FKBP12, SH2 domain, and cytochrome P450 are discussed and compared. Remaining difficulties and challenges are highlighted.

Modeling induced polarization with classical Drude oscillators: Theory and molecular dynamics simulation algorithm

A simple treatment for incorporating induced polarization in computer simulations is formulated on the basis of the classical Drude oscillator model. In this model, electronic induction is represented by the displacement of a charge-carrying massless particle attached to a polarizable atom under the influence of the local electric field. The traditional self-consistent field (SCF) regime of induced polarization is reproduced if these auxiliary particles are allowed to relax instantaneously to their local energy minima for any given fixed configuration of the atoms in the system. In practice, such treatment is computationally prohibitive for generating molecular dynamics trajectories because the electric field must be recalculated several times iteratively to satisfy the SCF condition, and it is important to seek a more efficient way to simulate the classical Drude oscillator model. It is demonstrated that a close approximation to the SCF regime can be simulated efficiently by considering the dynamics of a...

Molecular dynamics of the KcsA K(+) channel in a bilayer membrane.

Molecular dynamics (MD) simulations of an atomic model of the KcsA K(+) channel embedded in an explicit dipalmitoylphosphatidylcholine (DPPC) phospholipid bilayer solvated by a 150 mM KCl aqueous salt solution are performed and analyzed. The model includes the KcsA K(+) channel, based on the recent crystallographic structure of, Science. 280:69-77), 112 DPPC, K(+) and Cl(-) ions, as well as over 6500 water molecules for a total of more than 40,000 atoms. Three K(+) ions are explicitly included in the pore. Two are positioned in the selectivity filter on the extracellular side and one in the large water-filled cavity. Different starting configurations of the ions and water molecules in the selectivity filter are considered, and MD trajectories are generated for more than 4 ns. The conformation of KcsA is very stable in all of the trajectories, with a global backbone root mean square (RMS) deviation of less than 1.9 A with respect to the crystallographic structure. The RMS atomic fluctuations of the residues surrounding the selectivity filter on the extracellular side of the channel are significantly lower than those on the intracellular side. The motion of the residues with aromatic side chains surrounding the selectivity filter (Trp(67), Trp(68), Tyr(78), and Tyr(82)) is anisotropic with the smallest RMS fluctuations in the direction parallel to the membrane plane. A concerted dynamic transition of the three K(+) ions in the pore is observed, during which the K(+) ion located initially in the cavity moves into the narrow part of the selectivity filter, while the other two K(+) ions move toward the extracellular side. A single water molecule is stabilized between each pair of ions during the transition, suggesting that each K(+) cation translocating through the narrow pore is accompanied by exactly one water molecule, in accord with streaming potential measurements (, Biophys. J. 55:367-371). The displacement of the ions is coupled with the structural fluctuations of Val(76) and Gly(77), in the selectivity filter, as well as the side chains of Glu(71), Asp(80), and Arg(89), near the extracellular side. Thus the mechanical response of the channel structure at distances as large as 10-20 A from the ions in the selectivity filter appears to play an important role in the concerted transition.

Ion Permeation and Selectivity of OmpF Porin: A Theoretical Study Based on Molecular Dynamics, Brownian Dynamics, and Continuum Electrodiffusion Theory

Three different theoretical approaches are used and compared to refine our understanding of ion permeation through the channel formed by OmpF porin from Escherichia coli. Those approaches are all-atom molecular dynamics (MD) in which ions, solvent, and lipids are represented explicitly, Brownian dynamics (BD) in which ions are represented explicitly, while solvent and lipids are represented as featureless dielectrics, and Poisson-Nernst-Planck (PNP) electrodiffusion theory in which both solvent and local ion concentrations are represented as a continuum. First, the ability of the different theoretical approaches in reproducing the equilibrium average ion density distribution in OmpF porin bathed by a 1M KCl symmetric salt solution is examined. Under those conditions the PNP theory is equivalent to the non-linear Poisson-Boltzmann (PB) theory. Analysis shows that all the three approaches are able to capture the important electrostatic interactions between ions and the charge distribution of the channel that govern ion permeation and selectivity in OmpF. The K(+) and Cl(-) density distributions obtained from the three approaches are very consistent with one another, which suggests that a treatment on the basis of a rigid protein and continuum dielectric solvent is valid in the case of OmpF. Interestingly, both BD and continuum electrostatics reproduce the distinct left-handed twisted ion pathways for K(+) and Cl(-) extending over the length of the pore which were observed previously in MD. Equilibrium BD simulations in the grand canonical ensemble indicate that the channel is very attractive for cations, particularly at low salt concentration. On an average there is 1.55 K(+) inside the pore in 10mM KCl. Remarkably, there is still 0.17 K(+) on average inside the pore even at a concentration as low as 1microM KCl. Secondly, non-equilibrium ion flow through OmpF is calculated using BD and PNP and compared with experimental data. The channel conductance in 0.2M and 1M KCl calculated using BD is in excellent accord with the experimental data. The calculations reproduce the experimentally well-known conductance-concentration relation and also reveal an asymmetry in the channel conductance (a larger conductance is observed under a positive transmembrane potential). Calculations of the channel conductance for three mutants (R168A, R132A, and K16A) in 1M KCl suggest that the asymmetry in the channel conductance arises mostly from the permanent charge distribution of the channel rather than the shape of the pore itself. Lastly, the calculated reversal potential in a tenfold salt gradient (0.1:1M KCl) is 27.4(+/-1.3)mV (BD) and 22.1(+/-0.6)mV (PNP), in excellent accord with the experimental value of 24.3mV. Although most of the results from PNP are qualitatively reasonable, the calculated channel conductance is about 50% higher than that calculated from BD probably because of a lack of some dynamical ion-ion correlations.

Gating charge displacement in voltage-gated ion channels involves limited transmembrane movement

Voltage-gated ion channels are responsible for generating electrical impulses in nerves and other excitable cells. The fourth transmembrane helix (S4) in voltage-gated channels is the primary voltage-sensing unit that mediates the response to a changing membrane electric field. The molecular mechanism of voltage sensing, particularly with respect to the magnitude of the transmembrane movement of S4, remains controversial. To determine the extent of this transmembrane movement, we use fluorescent resonance energy transfer between the S4 domain and a reference point in the lipid bilayer. The lipophilic ion dipicrylamine distributes on either side of the lipid bilayer depending on the membrane potential, and is used here as a resonance-energy-transfer acceptor from donor molecules attached to several positions in the Shaker K+ channel. A voltage-driven transmembrane movement of the donor should produce a transient fluorescence change because the acceptor also translocates as a function of voltage. In Shaker K+ channels no such transient fluorescence is observed, indicating that the S4 segment does not translocate across the lipid bilayer. Based on these observations, we propose a molecular model of voltage gating that can account for the observed 13e gating charge with limited transmembrane S4 movement.

Ions and counterions in a biological channel: a molecular dynamics simulation of OmpF porin from Escherichia coli in an explicit membrane with 1 M KCl aqueous salt solution.

A 5 ns all-atom molecular dynamics trajectory of Escherichia coli OmpF porin embedded in an explicit dimyristoyl-phosphatidylcholine (DMPC) bilayer bathed by a 1 M [KCl] aqueous salt solution is generated to explore the microscopic details of the mechanism of ion permeation. The atomic model includes the OmpF trimer, 124 DMPC, 13470 water molecules as well as 231 K+ and 201 Cl-, for a total of 70,693 atoms. The structural and dynamical results are in excellent agreement with the X-ray data. The global root-mean-square deviation of the backbone atoms relative to the X-ray structure is 1.4 A. A cluster of three fully charged arginine (Arg42, Arg82, and Arg132) facing two acidic residues (Asp113 and Glu117) on L3 in the narrowest part of the aqueous pore is observed to be very stable in the crystallographic conformation. In this region of the pore, the water molecules are markedly oriented perpendicular to the channel axis due to the strong transversal electrostatic field arising from those residues. On average the size of the pore is smaller during the simulation than in the X-ray structure, undergoing small fluctuations. No large movements of loop L3 leading to a gating of the pore are observed. Remarkably, it is observed that K+ and Cl- follow two well-separated average pathways spanning over nearly 40 A along the axis of the pore. In the center of the monomer, the two screw-like pathways have a left-handed twist, undergoing a counter-clockwise rotation of 180 degrees from the extracellular vestibule to the pore periplasmic side. In the pore, the dynamical diffusion constants of the ions are reduced by about 50% relative to their value in bulk solvent. Analysis of ion solvation across the channel reveals that the contributions from the water and the protein are complementary, keeping the total solvation number of both ions nearly constant. Unsurprisingly, K+ have a higher propensity to occupy the aqueous pore than Cl-, consistent with the cation selectivity of the channel. However, further analysis suggests that ion-ion pairs play an important role. In particular, it is observed that the passage of Cl- occurs only in the presence of K+ counterions, and isolated K+ can move through the channel and permeate on their own. The presence of K+ in the pore screens the negative electrostatic potential arising from OmpF to help the translocation of Cl- by formation of ion pairs.

Generalized solvent boundary potential for computer simulations

A general approach has been developed to allow accurate simulations of a small region part of a large macromolecular system while incorporating the influence of the remaining distant atoms with an effective boundary potential. The method is called the Generalized Solvent Boundary Potential (GSBP). By representing the surrounding solvent as a continuum dielectric, both the solvent-shielded static field from the distant atoms of the macromolecule and the reaction field from the dielectric solvent acting on the atoms in the region of interest are included. The static field is calculated once, using the finite-difference Poisson–Boltzmann (PB) equation, and the result is stored on a discrete grid for efficient simulations. The solvent reaction field is developed using a basis-set expansion whose coefficients correspond to generalized electrostatic multipoles. A matrix representing the reaction field Green’s function between those generalized multipoles is calculated only once using the PB equation and stored ...