Multiscale methods for nanochemistry and biophysics in solution

Journal of Molecular Liquids 164 (2011) 101–112

Contents lists available at SciVerse ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Multiscale methods for nanochemistry and biophysics in solution Andriy Kovalenko a, b,⁎, Nikolay Blinov a, b a b

National Institute for Nanotechnology, 11421 Saskatchewan Dr., Edmonton, AB, Canada T6G 2M9 Department of Mechanical Engineering, University of Alberta, Edmonton, AB, Canada T6G 2G8

a r t i c l e

i n f o

Available online 12 October 2011 Keywords: 3D-RISM-KH molecular theory of solvation Molecular dynamics Amber package Ligand binding affinity Ion channel Solvation free energy Thiamine binding protein

a b s t r a c t Statistical–mechanical, 3D-RISM-KH molecular theory of solvation is an integral equation theory of molecular liquids capable of describing, from the first principles, solvation structure and thermodynamics in complex nanoscale systems. The 3D-RISM-KH theory accurately yields the solvation structure for biomolecular systems as large as chaperonins and ion channels. It predicts 3D maps of ligand binding affinity at once without using empirical scoring functions, and has significant potential for computer-aided drug design. Recently, the 3D-RISM-KH theory was coupled with MD simulations by contracting solvent degrees of freedom to treat solvated biomolecules, and the MD/3D-RISM-KH multiscale method was implemented in the Amber molecular dynamics package [T. Luchko et al., 2010, [29]]. The MD/3D-RISM-KH method allows one to study biomolecular processes on extremely long timescales, as the statistics of rare solvent and ligand events is accounted for analytically. We further replaced the MM/ GB(PB)SA post-processing which evaluates the free energy by empirical treatment of non-polar solvation effects with the MM/3D-RISM-KH method which employs 3D-RISM-KH to evaluate the solvation thermodynamics. In this paper, we briefly review the 3D-RISM-KH theory as well as the MD/3D-RISM-KH multiscale method implemented in the Amber package. We illustrate the MM/3D-RISM-KH method on binding of thiamine (vitamin of group B) to extracytoplasmic thiamine-binding lipoprotein MG289, and on the microscopic solvation properties of the bacterial Gloeobacter violaceus pentameric ligand-gated ion channel (GLIC) homologue. Crown Copyright © 2011 Published by Elsevier B.V. All rights reserved.

1. Introduction Molecular dynamics (MD) simulations with explicit solvent yield detailed modeling of solvated biomolecules for processes within accessible time scales [1]. A major computational burden comes from solvent molecules. Moreover, solvent enters pockets and inner cavities of proteins through a very slow process of their conformational changes, nearly as difficult to model as protein folding. The generalized Born model [2–7] represents the solvent polarization effects by a cavity in dielectric continuum (optionally, with Debye screening), whereas the non-electrostatic contributions are empirically parametrized. However, it bears the fundamental drawbacks of implicit solvation: hydrogen bonding and hydrophobic interactions are non-transferable to new complex systems (e.g. co-solvent or different buffer ions); the solvent accessible surface and volumetrics are not well defined, the entropic term is absent in continuum solvation. Also, the implicit solvation models do not account properly for the dispersion interactions and excluded volume effects [8,9] and require the phenomenological parameters (ionic radii, the surface tension coefficients) as input for modeling [5]. An alternative is statistical–mechanical, molecular theory of solvation, a.k.a. three-dimensional reference interaction site model [10–12], which starts from an explicit solvent model with all-atom force field ⁎ Corresponding author at: National Institute for Nanotechnology, 11421 Saskatchewan Dr., Edmonton, AB, Canada T6G 2M9. Tel.: +1780 641 1716; fax: +1780 641 1601. E-mail address: [email protected] (A. Kovalenko).

but operates with 3D correlation functions of species in a statistical ensemble rather than with individual molecular trajectories and predicts the solvation structure and thermodynamics of biomolecules from the first principles of statistical mechanics. Complemented with the Kovalenko–Hirata (KH) closure approximation [12], the 3D-RISM-KH theory properly accounts for chemical functionalities of biomolecules and solvent by representing both electrostatic and non-polar features of the solvation structure, such as hydrogen bonding, hydrophobicity, salt bridges, structural solvent, etc.; moreover, it analytically yields the solvation thermodynamics, including the solvation free energy, its energetic and entropic decomposition, and partial molar volume effects [12–16]. The approach was proven to be successful in description of the solvation structure and thermodynamics of proteins and protein complexes under various solvent conditions and composition [17–27]. It can be used for analysis of biomolecular systems as large and complex as the solvated chaperonin GroEL/ES protein complex [28]. Recently, we coupled the 3D-RISM-KH method contracting solvent degrees of freedom with molecular dynamics simulation of biomolecules in the Amber molecular dynamics package [29]. This included a number of accelerating schemes with several cutoffs for the interaction potentials and correlation functions, an iterative guess for the 3D-RISM solutions, and extrapolating solvent-induced forces and applying them in large multi-time steps (up to 20 fs) to enable simulation of large biomolecules. The multiscale MD/3D-RISM-KH method makes feasible modeling of biomolecular structures of practical interest and thus has tremendous potential for computer-aided drug design. It allows one

0167-7322/$ – see front matter. Crown Copyright © 2011 Published by Elsevier B.V. All rights reserved. doi:10.1016/j.molliq.2011.09.011

102

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

to study biomolecular systems involving solvent processes occurring on extremely long timescales, as the statistics of rare solvent events such as exchange of solvent molecules and ions, and ligand binding at pockets and inner spaces is accounted for statistically–mechanically. The 3D-RISM-KH theory combined with Amber molecular dynamics package [29], or in its standalone implementation, replaces the MM/GBSA or MM/PBSA post-processing suffering from the empirical treatment of non-polar contributions with the MM/3D-RISM-KH accurate evaluation of the solvation thermodynamics [29,27]. It also addresses the problem of molecular recognition [12,30,31] by yielding at once 3D maps of binding affinity without any phenomenological approximations [30–32], which can be used to improve the quality of 3D site binding maps used in fragment-based screening of drug compounds in most of the novel biomolecular software packages. It was shown recently that 3D-RISM-KH predicts binding maps of potential inhibitors of the pathological conversion prion proteins [33]. Below we briefly summarize the 3D molecular theory of solvation and the coupling scheme of the MD/3D-RISM-KH method with the Amber molecular dynamics package [29]. Then we employ the standalone version of the theory and the MM/3D-RISM-KH method to predict binding maps, binding affinity and microscopic solvation structure of thiamine at extracytoplasmic thiamine binding lipoprotein MG289 [34]. For illustration of the capabilities of the 3D-RISM-KH theory to deal with large-scale biomolecular systems, we analyze the solvation structure, permeation and selectivity of the bacterial Gloeobacter violaceus pentameric ligand-gated ion channel (GLIC) homologue [35]. 2. Theory 2.1. Three-dimensional molecular theory of solvation The solvation structure is represented by the probability density of finding site γ of solvent molecules at 3D space position r around the solute (supra) molecule, ργgγ(r), which is determined by the average number density ργ in the solution bulk times the 3D distribution function (normalized density distribution) gγ(r) of solvent site γ. The latter indicates site density enhancement when gγ(r) N 1 or depletion when gγ(r) b 1 relative to the average density at a distance from the solute in the solution bulk where gγ(r) → 1. The 3D distribution functions of solvent interaction sites are obtained from the 3D-RISM integral equation [12]

 0 ′ ′ ð1Þ hγ ðrÞ ¼ ∑ ∫ dr cα r−r Þχαγ r ; α

where hγ(r) is the 3D total correlation function of solvent site γ related to the 3D site distribution function by gγ(r) = hγ(r) + 1, and cγ(r) is the 3D direct correlation function which has the asymptotic of the solute–solvent site interaction potential, cγ(r) ∼−uγ(r)/(kBT); the site–site susceptibility of pure solvent χαγ(r) is an input to the 3D-RISM theory; and indexes α and γ enumerate all sites on all sorts of solvent species. Another relation between the 3D total and direct correlation functions, called a closure, is necessary to complement the 3D-RISM integral equation (1). The exact closure can be formally expressed as a series in terms of multiple integrals of the combinations of the total correlation functions [36]. However, it is extremely cumbersome, and in practice is replaced with amenable approximations. We use the 3D-KH closure approximation, proven to be appropriate to describe various association effects in complex liquids and electrolyte solutions [12], and in supramolecular synthetic organic [2–4] and biomolecular [5–11] systems in solution,

gγ ðrÞ ¼

f


 exp dγ ðrÞ

for dγ ðrÞ≤0;

1 þ dγ ðrÞfor dγ ðrÞN0; uγ ðrÞ þ hγ ðrÞ−cγ ðrÞ; dγ ðrÞ ¼ − kB T

ð2Þ

where uγ(r) is the 3D interaction potential between the whole solute and solvent site γ specified by the molecular force field, and kBT is the Boltzmann constant times the solution temperature. The closure (2) couples in a non-trivial way the so-called mean spherical approximation (MSA) and hypernetted chain (HNC) approximation [36], the former being applied to the spatial regions of solvent density enrichment (gγ(r) N 1) such as association peaks, and the latter to those of solvent density depletion (gγ(r) b 1) including the repulsive core [12]. The site–site susceptibility of solvent breaks up into the intra- and intermolecular terms, χαγ ðr Þ ¼ ωαγ ðr Þ þ ρα hαγ ðr Þ;

ð3Þ

where the intramolecular correlation function


 2 ωαγ ðrÞ ¼ δαγ δðrÞ þ 1−δαγ δ r−lαγ = 4πlαγ represents the geometry of solvent molecules with site–site separations lαγ specified by the molecular force field (z-matrix in quantum chemistry), and hαγ(r) is the radial total correlation function between sites α and γ. In advance to the 3D-RISM-KH calculation, hαγ(r) are obtained from the dielectrically consistent RISM theory [37] coupled with the KH closure (DRISM-KH) [12], applied to the bulk solvent with counter ions, co-solvent, and ligands at a given concentration. DRISM has an electrostatic bridge correction in the analytical form which was shown to be proper for consistently getting the dielectric constant of electrolyte solution through the dielectric and ionic routes [37]. The susceptibility (3) of the bulk solution is then input into the 3D-RISM integral equation (1). The solvation free energy of the solute supramolecule in multicomponent solvent following from the 3D-RISM-KH integral equations (1) and (2) is given by the closed analytical expression [12]  
 1 1 2 μsolv ¼ kB T ∑ ργ ∫dr hγ ðrÞΘ −hγ ðrÞ − hγ ðrÞcγ ðrÞ−cγ ðrÞ ; 2 2 γ

ð4Þ

where Θ(x) is the Heaviside step function. The potential of mean force Wγ(r) acting on solvent site γ near the biomolecule is defined in terms of the 3D site distribution function as Wγ(r) =−kBTlngγ(r). For a structural water molecule “w”, for example, used here as an illustration of a simple ligand molecule, the binding strength is determined as the difference between the potential of mean force of water oxygen O in the effective potential well at the biomolecule and that in the first peak of bulk water. The binding strength is thus expressed in terms of the water oxygen peaks at the biomolecule and in bulk water as binding

max at solute

max in bulk

Aw ¼W −WO
O  max at solute max in bulk ¼ −kB Tln gO ðrÞ=gO ðrÞ :

ð5Þ

The solvation free energy (4) and binding energy (5) obtained by 3D-RISM-KH can be decomposed into partial contributions of solute sites, providing a basis for spatial decomposition analysis (SDA) of association effects [38]. To properly treat electrostatic forces in electrolyte solution with polar molecular solvent and ionic species when evaluating the convolution in the 3D-RISM-KH integral equations (1)–(2), the radial correlations (3) of bulk solvent, and in the thermodynamics expressions (4) and (5), the electrostatic asymptotics of all the correlation functions (both the 3D and radial ones) are treated analytically [1,14]. The non-periodic electrostatic asymptotics are separated out in the direct and reciprocal space and the remaining short-range terms of the correlation functions are discretized on a rectangular grid in a non-periodic uniform box large enough to ensure decay of the short-range terms at the box boundaries [14]. The convolution of the short-range terms in Eq. (1) is calculated using 3D fast Fourier transform (3D-FFT) in double size box to prevent aliasing.

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

Accordingly, the electrostatic asymptotic terms in the thermodynamics integral (4) are handled analytically, reducing to one-dimensional integrals that are easy to compute [14]. This analytical treatment of the electrostatics eliminates the periodicity artifacts arising when using the Ewald summation and 3D-FFT supercell for the electrostatic terms, in particular, those in the values of gγ(k=0) coming from the cancelation or renormalization of the Coulomb singularities of cγ(k→0) [1,14]. Eqs. (1)–(2) are solved by using the modified direct inversion in the iterative subspace (MDIIS) accelerated numerical solver of integral equations of liquid state theory [12]. The convergence to a relative root-mean-square of the residual term of 10 − 4−10 − 5 provides a required level of accuracy for the solution. 2.2. MD coupled with 3D-RISM-KH In the coupled MD/3D-RISM-KH method implemented in Amber molecular dynamics package [29], MD is applied to the biomolecule while the solvation structure, free energy and associated forces are obtained by 3D-RISM-KH. The latter are derived by differentiating the Kirkwood thermodynamic integration formula analytically, also valid for the solvation free energy (4) in the 3D-KH approximation (2), FðRi Þ ¼ –

∂μsolv ¼ – ∑ ργ ∫drg ðrÞ ∂ui ðr−Ri Þ ; γ ∂Ri γ ∂R

ð6Þ

i

where uγ(r−Ri) is the pairwise interaction potential between biomolecule site i located at Ri and solvent site γ at r. To obtain meaningful sampling of solute conformations, we reduce the computational expense of 3D-RISM calculations by several optimization strategies [29]: (i) creating high-quality initial guesses to the direct correlation function cγ(r) from multiple previous solutions; (ii) accelerating the pre- and postprocessing of the solute–solvent potentials, long-range asymptotics, and forces using a cutoff scheme and minimal solvation box; and (iii) avoiding direct calculation of 3D-RISM-KH solvation forces fFgðkÞ at current time step k by interpolating them based of those from previous time steps, or force-coordinate extrapolation multi-time step procedure, ðkÞ

fFg

¼

N X

ðlÞ

akl fFg ; l∈3D−RISM steps;

ð7Þ

l¼1

where the weight coefficients akl are obtained as the best representation of the matrix of arrangement of solute atoms fRgðkÞ at current time step k in terms of its projections onto the “basis” of l = 1,N previous ones ðlÞ fRg obtained from the MD/3D-RISM-KH, by minimizing the norm of the difference  2 N   X  ðkÞ ðl Þ  akl fRg  : minimizefRg −   l¼1

Combined with other accelerating techniques, such as propagation of solution of the 3D-RISM equations, speed-up of up to 10 times compared to the basic implementation of 3D-RISM-KH can be achieved [29]. Optimal choice of parameters for the force extrapolation depends on the system under consideration and type of calculations (accurate energy conserving dynamics vs efficient sampling of conformational space of a macromolecule in solution). As shown in Ref. [29], with this method 3D-RISM-KH calculations can be performed at up to 20 fs intervals, which corresponds to N = 10 with the base time step of 2 fs. 2.3. Force field parameters for 3D-RISM-KH calculations Details of setup of 3D-RISM-KH calculations depend on a particular application, but there are some common features summarized in

103

this section. In the current study, the 3D-RISM-KH Eqs. (1)–(2) for all the systems were solved in rectangular boxes on the grid with the resolution of 0.5 Å. Such a grid is fine enough to obtain the 3DRISM solvation free energies with an accuracy better than 0.5% for most of the systems. The box size for each system was chosen to be sufficiently large to avoid issues with boundary effects. The all-atom ff03 force field [39,40] from the Amber 9 molecular dynamics package [41,42] was used to parametrize the solute–solvent interaction potentials. For water, the SPC/E model [43] was used with the LennardJones parameters for the hydrogen site of water corrected to properly describe the thermodynamic and structural properties of bulk water, based on the one-dimensional dielectrically consistent RISM theory [37]. The partial charges for water oxygen and hydrogen are qO =−0.8476 e and qH =0.4238 e (e is the electron charge), the Lennard-Jones parameters for oxygen are σO =3.166 Å and εO =0.1554 kcal/mol, and those for hydrogen are σH =0.4 Å and εH =0.046 kcal/mol. All the thermodynamic properties were calculated for ambient water at temperature T=298.15 K, number density ρ=0.997 g/cm3, and dielectric constant ε=78.4. The temperature derivatives for the solvation entropy were calculated using the first-order central finite difference scheme with a temperature step of 2 K. 2.4. Molecular dynamics simulations and energy minimization To sample conformational space of the thiamine–MG289 complex in bound and unbound states, molecular dynamics trajectories were generated in all-atom explicit solvent molecular dynamics simulations by using Amber 9 molecular dynamics package [41,42] with the ff03 all-atom force field [39,40]. The following protocol was used. First, counter-ions were added at the points of extrema of the electrostatic potential to neutralize a system (two sodium ions were added in the case of the MG289 protein in the unbound state, one sodium ion was added for the thiamine–MG289 complex, and one chloride ion for the thiamine molecule in the unbound state). The initial structures were solvated with the SPC/E [43] water molecules in rectangular periodic boxes. A distance of at least 12 Å was allowed between the solute molecules and periodic boundaries. The above procedure was followed by two minimization runs to relax solvent and solute (the protein and thiamine) degrees of freedom. The production MD simulations were preceded by two thermalization/equilibration runs. First, the temperature was gradually raised from 0 to the target value of 298.15 K. Then, the constant pressure MD was performed to relax solvent degrees of freedom and to adjust the solvent density to approximately 1 g/cm 3. In the equilibration runs, the coordinates of the protein and thiamine were restrained to their minimum energy values. The MD step was set to 1 fs. In the production runs, all bonds involving hydrogen atoms were constrained with the SHAKE algorithm [44], which allowed us to set the MD step to 2 fs. Default value of 10 − 5 Å was used for relative geometrical tolerance for coordinate resetting in SHAKE algorithm. Langevin dynamics was used for temperature control with the collision frequency of 5 ps − 1. Isotropic position scaling was used to control pressure with pressure relaxation time set to 2.0 ps. The cutoff for nonbonded interactions was 16 Å. Other parameters for nonbonded interactions, including the parameters for the Particle Mesh Ewald (PME) method used to treat the electrostatics, were set to their default values [42]. 3. Results and discussion 3.1. Multiscale modeling based on 3D-RISM-KH: protein–ligand binding and free energy calculations The 3D-RISM-KH theory can be combined with different molecular modeling algorithms allowing one to treat variety of biophysical and biochemical problems on the same footing. In the case of protein–ligand binding the 3D-RISM-KH theory can be used for

104

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

molecular recognition and ligand mapping on protein surfaces, initial docking, refinements of binding poses of ligand in a binding pocket [30,19,31], studying the effects of structural solvation on binding, and calculations of the free energy of binding [45]. In this section, we employ the 3D-RISM-KH-based approach to study binding modes and free energy of binding of thiamine (the vitamin of group B) against extracytoplasmic thiamine binding lipoprotein MG289. Mycoplasma genitalium (MG) is one of the smallest free-living bacterium, and it is of interest in the context of studying smallest organisms capable of self-replication. The high resolution structure the thiamine–MG289 complex has been recently solved in X-ray experiments [34]. The binding site of thiamine is buried deep in the protein fold (Figs. 1 and 2) which makes the analysis based on the standard docking methods a challenging task. As discussed in Ref. [34], hydrophobic effects play an important role in binding, with the bulk of the interactions between the protein and thiamine being non-polar. Also, because of the non-zero charge of thiamine (+1 e), the desolvation effects and specificity of electrostatic interactions may play important role in binding. Thiamine forms a direct hydrogen bond with Asp41 and a water-mediate hydrogen bond with Asp269 residues. The pyrimidine ring of thiamine π-stacks with Trp275 [34], and possibly with Tyr307. The corresponding residues are shown in Panel (C) of Fig. 1. The initial structure of the complex was taken from the Protein Data Bank (PDB accession code 3MYU). Solvent (water and ions) as well as the ligand molecule (thiamine) was described with the 3D-

(A)

(B)

(C)

Fig. 1. Panels (A) and (C): thiamine in the binding pocket of the MG289 monomer [34]. Residues within a distance of 4 Å from thiamine are shown in the stick representation. OD2 of Asp41 forms a hydrogen bond with thiamine, Asp269 is involved in water-mediated hydrogen bonding (the water molecule is represented by the oxygen shown as a van der Waals sphere), Trp275 residue π-stacks with the pyrimidine ring of thiamine [34]. Panel (B): thiamine molecule and its fragments as used in the 3D-RISM-KH calculations. The figure was prepared with DiscoveryStudio 2.0 (Accelrys Inc.) software.

RISM-KH method. Non-polar hydrogen atoms of thiamine were merged with adjacent carbon atoms. The force field parameters for the inter-protein and protein–solvent (which includes the thiamine molecule) interactions were assigned based on the Amber ff03 force field [39,40]. To generate conformations for the free energy calculations, all atom explicit solvent molecular dynamics simulations were performed with the Amber molecular dynamics package [41,42] as described in Section 2.4. 3.1.1. Protein–ligand binding: 3D-RISM-KH for mapping ligand binding sites Following the framework developed in Refs. [30,19,31], we treat ligand (the thiamine molecule) degrees of freedom as a part of the complex solvent which also includes water molecules and chloride counter-ions. The 3D density distributions for all solvent sites (the oxygen and hydrogen sites of water, chloride ions and 21 atomic sites of thiamine, which excludes non-polar thiamine hydrogens) have been calculated with the 3D-RISM-KH method for the experimental X-ray structure of the protein. The densities for selected atomic sites of the thiamine molecule (pyrimidine C carbon and N2 nitrogen, thiazole sulfur and hydroxyl oxygen sites) are represented with the isosurfaces corresponding to isosurface level of 2 (excess of number density relative to its bulk value) in Figs. 2 and 3. In the right panel of Fig. 2, the front parts of the density distributions have been clipped to reveal a protrusion of these densities into the binding pocket. The high density pockets of the thiamine atomic sites can be related to the most probable locations of these sites in the bound state. As discussed in the previous section, the density distributions of ligand atomic sites can be used to define a potential of mean force (PMF) which fully accounts for all the physical protein–ligand interactions as well as for the solvation effects, including the favorable hydrophobic interactions and the desolvation penalty. Such PMFs can be used in scoring functions in the docking algorithms to find/select binding conformations with the lowest energy. Compared to the conventional grid-based docking algorithms implemented, for example, in the AutoDock program [46], the new docking approach is characterized by the high level of transferability of its solvation part, which can account for the competition between different ligand species, ligand concentration and solvent composition effects as well as for the thermodynamic state of a system under consideration. The 3D-RISM-KH based docking approach also provides direct information on solvation enthalpic and entropic contributions to the binding energy. The latter is especially important because the solvation entropic effects are related to hydrophobicity. Additionally, the molecular theory of solvation accounts for the specificity of ligand site–protein interactions through the geometric restraints imposed on ligand degrees of freedom during calculations of 3D density distributions for atomic sites of a ligand molecule. For further illustration, we show the 3D density distributions for the sulfur atomic site (located at the thiazole ring of thiamine) and the OH group of thiamine in the binding pocket of MG289 in Fig. 4. The isosurfaces shown here are for isolevel of 1.5 (excess of number density over the bulk value) which tentatively corresponds to the first “solvation” shells of thiamine atomic sites in the proximity of the protein. While there are some correlations between these distributions and locations of the corresponding sites of thiamine in binding conformation, the higher density peaks are more specific and better describe the possible locations of ligand sites in the binding pocket. We present such densities by isosurfaces with isovalue of 2 for C carbon and N2 nitrogen sites of the pyrimidine ring of thiamine on the left and right panels of Fig. 5, respectively. The carbon and nitrogen densities are clearly confined in the plane of the pyrimidine ring of thiamine in the experimental conformation, with a good agreement between the locations of C and N2 sites in the binding pocket from one hand, and the density distributions for the same sites, from the other. This planar density pocket is clearly seen in

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

(A)

105

(B)

Fig. 2. Panel (A): Thiamine in binding pocket of MG289. Panel (B): Isosurface representation of the 3D density distributions of the pyrimidine carbon (in cyan) and thiazole sulfur (in yellow) atomic sites. Isosurfaces correspond to isovalue of 2. Front parts of the isosurfaces have been clipped to show the binding pocket. In both panels, the solvent accessible surface of the protein is shown in gray. The image was produced by using the VMD program [63].

Fig. 5 just below Trp275 residue. While such planar structure is a characteristic of 3D density distributions for all ligand sites, energetically most favorable orientation of the thiamine molecule will be when one of its aromatic rings commensurates with this planar structure, thus maximizing inter-molecular van der Waals and, possibly, hydrophobic interactions. There will be steric clashes between the protein side chains and thiamine in the case of the thiazole ring of thiamine located below Trp275. Also, the thiazole sulfur has preferable locations outside that density pocket. Thus, the planar structure from Fig. 5 gives a probable location of the pyrimidine ring of thiamine in the binding pocket. The density distributions for C carbon and N2 nitrogen sites are also consistent with the orientation of the thiamine molecule in the experimental binding pocket. There are two major reasons why there is no “one-to-one” agreement between locations of the actual atomic sites of the thiamine molecule and the densities profiles. First, the 3D-RISM-KH method relies on 3D density distributions for atomic sites of a solvent molecule (which includes the thiamine molecule in this particular case), and these densities are averaged over all orientations of a solvent molecule, which smooths out actual six dimensional distributions of solvent molecules around a protein. Secondly, in our calculations the thiamine molecule has been represented by three independent fragments (Fig. 1) to account for its flexibility. Thus, the additional geometrical restraints should be imposed in 3D-RISM based docking algorithms to select real binding conformations. This can be done, for example, by combining the PMF generated with the 3D-RISM-KH approach with the intra-ligand potential terms (bonded and nonbonded). 3.1.2. Free energy of binding: MM/3D-RISM-KH approach To calculate free energy of binding for the thiamine–MG289 complex, we use the approach which combines the 3D-RISM-KH method to calculate the solvation free energy with molecular mechanics for the gas-phase (intra- and inter-molecular) energy [27]. The protein and ligand conformational ensembles are sampled in the explicit solvent molecular dynamics simulations. Previously, a similar approach

was used for calculation of the binding energy in Ref. [45]. The approach can be regarded as an alternative to the MM/PBSA(GBSA) method [41,47,48] where the Poisson–Boltzmann (PB) equation, or the generalized Born (GB) approximation combined with the surface area (SA) term, have been replaced with the 3D molecular theory of solvation. This allows one to account rigorously for the non-polar (including hydrophobic) solvation effects, and it also makes the theory highly transferable and applicable to different solvents at different thermodynamic conditions, being free from any adjustable phenomenological parameters. The free energy of binding ΔΔGbind is defined as
 ΔΔGbind ¼ ΔGcomplex − ΔGprotein þ ΔGligand ;

ð8Þ

where ΔGcomplex is the free energy of a protein–ligand complex, ΔGprotein and ΔGligand are the free energies of a protein and a ligand in unbound states. Each term in the above equation can be represented as

〈

〉 〈

〉 〈 〉

ΔG ¼ ΔHgas þ ΔGsolv −T ΔS ;

ð9Þ

where 〈ΔHgas〉 is the gas-phase energy, 〈ΔGsolv〉 is the solvation free energy and − T〈ΔS〉 is a contribution of the conformational entropy of a solute molecule to the free energy of a system. Angle brackets represent a canonical average over instantaneous (calculated for every molecular conformation) quantities. As discussed in Section 2.4, we employ explicit solvent all-atom molecular dynamics simulations with the Amber 9 molecular dynamics package [41,42] to generate canonical ensembles of protein and ligand conformations in bound and unbound states for calculations of average quantities in Eq. (9). The above approach requires sampling of conformational space of protein and ligand in both bound and unbound states. At the same time, a simplification is commonly used which replaces true ΔGprotein and ΔGligand with the quantities calculated using conformations sampled for a protein–ligand complex in the bound state only, thus eliminating the necessity of additional MD runs for a protein and a ligand

106

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

(A)

(B)

(C)

(D)

Fig. 3. Isosurface representation of 3D density distributions for pyrimidine C carbon and N2 nitrogen atomic sites (Panels (A) and (B)), thiazole sulfur and hydroxyl oxygen (Panels (C) and (D)) atomic sites. Isosurfaces correspond to isolevel of 2. The solvent accessible surface of the protein is shown in gray. The image was produced by using the VMD program [63].

in the unbound states. Obviously, this is a reasonable approximation only in the case when protein–ligand binding does not lead to major conformational changes of participating molecules. To assess the validity of such approximation in the particular case of the thiamine–MG289 complex, we compared the root mean square displacement (RMSD) and the mean square fluctuations (B-factors) of backbone α-carbon atoms of the protein in bound and unbound states. The results are represented in Figs. 6 and 7, respectively. The

complex and protein in the unbound state are stable with the evolution of RMSDs in the similar ranges at the time scale of simulations. Also, the mean square fluctuations, which characterize flexibility of a system, show similar patterns with only minor variations. It means that calculations of the free energy of binding using protein and ligand conformations sampled from an equilibrium ensemble of a complex only may be a quite reasonable approximation for the thiamine–MG289 system. Also, in this particular case the entropic

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

(A)

107

(B)

Fig. 4. Isosurface representation of thiazole sulfur (right panel), and hydroxyl oxygen and hydrogen (left panel) atomic sites of thiamine. The isosurfaces correspond to isolevel of 1.5 (number density excess over the bulk value) which can be tentatively referred as first “solvation” shells of ligand atomic sites around the protein. The image was produced by using the VMD program [63].

conformational penalty might be less important compared to the other components of the free energy of binding (the gas-phase energy and the solvation free energy). Both thiamine and MG289 have non-zero charges at normal conditions and neutral pH. Most likely, under the physiological conditions the molecules will be associated with ions. The effects of counter-ions (structural ions) on binding may be important to account for to make a connection between theoretical calculations for binding energy and experimental results. The locations of ions in

(A)

the explicit solvation schemes should be sampled with MD or Monte Carlo methods along with a solute and water degrees of freedom. Accounting for explicit ions in the free energy calculations might be also useful from the point of view of the Ewald summation because it guarantees the neutrality of the periodic cell. In low ionic concentration environments, it can be computationally challenging to obtain statistically significant results for ion distributions. As an alternative, ions can be placed at the locations of extrema of the electrostatic potential, or they can be treated implicitly with the PB or

(B)

Fig. 5. Isosurface representation of thiamine C carbon (left panel) and N2 nitrogen (right panel) sites from the pyrimidine ring of thiamine. The isosurfaces are for isolevel of 2. The image was produced by using the VMD program [63].

108

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

Fig. 6. Root mean square displacements of α-carbon atoms from the initial structures in the production MD runs for thiamine–MG289 complex and MG289 in unbound state.

GB methods. In the framework of the 3D-RISM-KH-based approaches, ions are considered as separate solvent species, with account for all the correlations relevant to ions, including ion–ion correlations. In addition to the Debye screening accounted for in the implicit solvation methods, 3D-RISM-KH also describes microscopic solvation structural effects, including specific interactions between ions and biomolecules. For charged solutes, the simulations box should be large enough to accommodate asymptotics of the distribution functions (the electrostatic potential in the case of PB calculations). In the case of the 3DRISM-KH method, the analytical treatment of the long distance behavior of the total and direct correlation functions has been proposed [49]. It allows one to use 3D-RISM-KH for calculations of the solvation free energy of charged solutes for different solvent compositions under different thermodynamic conditions. For relatively simple solvents (water and counter ions), realistic estimates for the electrostatic part of the solvation energy can be obtained with PB or GB methods, and the 3D-RISM-KH method can complement these methods by providing accurate results for the non-polar part of the solvation free energy and the entropic contribution to the solvation free energy. The gas-phase and the solvation free energy contributions to the binding energy of the thiamine–MG289 complex along with their physical components are summarized in Table 1. The binding energy

ASP41 TYR167 HIS42 TYR176 ASP269 TYR275

TYR307 ASP308 ILE343

was calculated based on Eqs. (8) and (9). The averages of the instantaneous quantities were calculated using 100 snapshots taken at equal separations from the 0.5 ns equilibrated parts of all-atom explicit solvent trajectories of 10 ns for the thiamine–MG289 complex, MG289 monomer and thiamine in an unbound state, simulated with Amber 9 [41,42] molecular dynamics package and ff03 all-atom force field [39,40]. The solvation free energy and its physical components for the complex, protein and ligand in sodium chloride water solution at 0.208 mol/kg concentration were calculated by using the 3D-RISM-KH method. We also decomposed the solvation part of the binding energy into electrostatic and non-electrostatic parts. They are 43 (11) and 11 (3) kcal/mol, respectively (the numbers in the parentheses giving the statistical errors). The convergence of the MM/3D-RISM-KH free energy as a function of number of snapshots was discussed elsewhere [27]. Both the electrostatic and non-electrostatic solvation effects are the factors unfavorable for thiamine binding against MG289. This is a common situation for association of oppositely charged solutes. The favorable electrostatic gas-phase (inter-molecular) part of the binding energy overcompensates the solvation electrostatic contribution to the binding energy. Electrostatic desolvation penalty is a consequence of removing solvent particles from polar solvent exposed surfaces of a protein and ligand upon formation of the protein–ligand interface. This results in decrease of favorable direct electrostatic interactions between a charged solute and atomic sites of solvent (polarized in the proximity of a protein). Additionally, binding reduces the amount of hydrogen bonds formed between a charged solute (protein and ligand molecules) and water molecules. The non-electrostatic contribution to the desolvation penalty is mostly due to excluding solvent particles from (favorable) van der Waals interactions between solvent and solute molecules. The unfavorable solvation effects are partly compensated by formation of specific hydrophobic contacts, water-mediated hydrogen bonds and salt bridges. For example, it follows from Table 1 that the solvation entropic effects, closely related to hydrophobicity, play a favorable role in thiamine–MG289 binding, which supports the suggestions from Ref. [34]. The exact balance between different factors and, as a consequence, a sign and extent of changes in the solvation free energy upon protein–ligand binding depends on both details of interactions between all species involved and also solvent conditions (such as, for example, temperature and level of pH).

3.1.3. Desolvation effects in protein–ligand binding Displacement and change in orientation of solvent particles upon binding can be represented by the incremental density distribution defined as a difference between solvent densities (for each solvent site) of a protein–ligand complex and a protein in the unbound state. Because solvent particles in binding pocket (including possible structural ions and water molecules) can directly affect pathways of binding and binding energetics (binding affinity), the incremental densities can provide valuable information on binding. In Figs. 8 and 9 we display incremental densities for the oxygen and hydrogen sites of water, and sodium and chloride ions in the binding pocket of thiamine at MG289 calculated for the experimental X-ray structure of the complex. The isosurfaces shown on the left panel of Fig. 8 with isovalue of −1 (negative values correspond to the depletion of solvent densities as a result of binding) tentatively represent the volume

Table 1 Physical components of binding energy for MG289–thiamine complex. Statistical errors are given in parentheses. All values are in kcal/mol. Fig. 7. Mean square fluctuations (B-factors) of α-carbon atoms resolved on per-residue basis for the thiamine–MG289 complex and MG289 in the unbound state. Vertical lines show locations of the residues in the binding pocket. B-factors are normalized as follows: Bi = 8π2/3〈Δri2〉, where angle brackets denote an average over MD conformations.

〈ΔΔHgas〉 − 108 (12) 〈ΔΔHelect〉 〈ΔΔHvdW〉 − 72 (11) − 37 (3.5)

〈ΔΔHbonded〉 0.9 (0.7)

〈ΔΔGsolv〉 54 (11) 〈ΔΔHsolv〉 142 (15)

− T〈ΔΔSsolv〉 − 88 (11)

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

(A)

109

(B)

Fig. 8. Incremental densities for water oxygen (in red) and hydrogen (in blue) atomic sites in the thiamine binding pocket of MG289. Panel (A): density depletion (desolvation) upon thiamine binding is illustrated by the isosurfaces corresponding to isolevel of − 1. Panel (B): isosurface representation of the excess water density of the thiamine–MG289 complex relative to MG289 in the unbound state. The isosurfaces are for isolevel of 1. Thiamine in the binding pocket is shown with van der Waal spheres. In both panels, binding site residues (Asp41, His42, Tyr167, Tyr176, Asp269, Tyr275, Tyr307, Asp308 and Ile343) are shown in the stick representation. The image was produced by using the VMD program [63].

occupied by the water molecules displaced from the binding pocket. In contrast, the right panel of the same figure shows the water brought to the binding pocket by the thiamine molecule. As expected, this mostly includes water molecules close to the (positively charged) thiazole ring of thiamine. It is clearly seen from this figure that the hydrogen and oxygen sites of water are affected differently by the

(A)

binding. The implication is that simple desolvation potentials commonly used in the standard docking software packages might not account properly for solvent polarization effects. Also, a pattern in the distributions of water oxygen and hydrogen sites from the right panel of Fig. 8 suggests a possibility of formation of hydrogen bonds between water molecule and thiamine. The difference between

(B)

Fig. 9. Incremental densities for sodium (in blue) and chloride (in orange) ions in the thiamine binding pocket of MG289. Panel (A): density depletion (desolvation) upon thiamine binding is illustrated by the isosurfaces with isovalue of − 0.5. Panel (B): isosurface representation of the excess ions densities of the thiamine–MG289 complex relative to MG289 in the unbound state. The isosurfaces are for isolevel of 0.5. In both panels, the binding site residues (Asp41, His42, Tyr167, Tyr176, Asp269, Tyr275, Tyr307, Asp308 and Ile343) are shown in the stick representation. The image was produced with the VMD program [63].

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

incremental densities for sodium and chloride ions shown in Fig. 9 can be explained by the presence of acidic residues in the binding pocket (Asp41, Asp269 and Asp308 are all located within 4 Å distance from the thiamine molecule in the binding pocket) deprotonated at neutral pH. Only one residue, His42, which may be positively charged at normal conditions is located close to the binding pocket. As a result, in the unbound state there is an excess density of sodium ions to be expelled from the binding pocket by the thiamine molecule. This is illustrated on the left panel of Fig. 9 where the incremental densities of sodium and chloride ions are represented by the isosurfaces with isovalue of −0.5, which corresponds to the depletion of the densities. On the contrary, in the bound state, there is an excess density of chloride ions, mostly close to the solvent exposed positively charged thiazole ring of thiamine (the right panel of Fig. 9). Thus, in the particular case of thiamine–MG289 binding, the desolvation effects may include displacement of cations (possibly, structural) from the binding pocket, which may have important implications for the pathway of binding and the free energy of binding under the physiological conditions.

1.1 Experiment DRISM

1.05 1

Activity coefficient

110

0.95

NaCl in water

0.9 0.85 0.8 0.75 0.7 0.65 0.6 0

0.2

0.4

0.6

0.8

1

Concentration, mol/kg Fig. 10. Mean activity coefficient for sodium chloride water solution at normal conditions plotted as a function of NaCl concentration. Comparison between the theoretical (1D-DRISM-KH) and experimental [60] results.

3.2. Microscopic solvation and selectivity of GLIC ion channels Ionic channels belong to an important class of biomolecular systems where solvent (ions, water, protons and possible ligands) performs its biological functions in the restricted geometry. Solvents confined in restricted geometries can have properties quite different from those in the bulk state, which may have interesting implications for biologically relevant phenomena [50]. Because passing of ions through a channel is a rare event (on the time scale of relaxation of the bulk solvent), the explicit solvent simulations [51–53] require extensive computational resources to achieve statistically reliable results. Additionally, studying of selectivity mechanisms and permeation properties of channels under the physiological conditions may require modeling of complex solvent conditions, including low ionic concentration environments, which is currently beyond the scope of explicit solvent molecular dynamics simulations. To overcome these difficulties, the coarse-grained models of a channel were proposed and successfully used to study the physical properties of ionic channels [54–56]. While such methods provide an important physical insight into functioning of channels (in particular, the importance of protein polarization and side chains mobility on the selectivity of sodium channels have been recently elucidated [54]), they do not account for the atomic structure of a real ionic channel. Also, they do not allow one to account in a transferable manner for more complex solvent compositions, where solvent includes different cofactors, ligands and other small molecules. This can be done based on the 3D molecular theory of solvation. Thus, the 3D-RISM-KH method provides a natural bridge between the coarse-grained description and all-atom explicit solvent molecular dynamics simulations. It was proven that the method is capable of describing solvent in restricted geometries, including protein pores and cavities [28,32,57,30,58,18,26]. 3.2.1. Activity coefficient The 3D-RISM-KH method requires bulk solvent 1D distribution functions as input. These functions can be obtained based on the 1D dielectrically consistent RISM approach (1D-DRISM). Here we show that 1D-DRISM gives a quantitatively accurate description of the activity coefficient for sodium chloride water solution in a wide range of salt concentrations [58,30]. The mean activity coefficient is plotted as a function of concentration for sodium chloride water solution at normal condition in Fig. 10. The coefficient is defined as pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ðγCl− γNaþ Þ, where γX = exp(ΔμX/kBT−ΔμX0/kBT) is the activity coefficient for ion X. Here, ΔμX and ΔμX0 are the solvation free energies of ion X in salt solution and in water at infinite dilution, respectively. The activity describes a deviation of the chemical potentials of the

mixture of chemical substances from the ideal behavior. The theoretical results represented in Fig. 10, were obtained with the modified OPLS (Optimized Potentials for Liquid Simulations) force field [59] parameters (the van der Waals radii and energies for chloride atoms are σCl − = 3.771 Å, εCl − = 0.2650 kcal/mol and those for sodium ions are σNa + = 3.328 Å, εNa + = 0.0028 kcal/mol). As in the previous sections, the SPC/E model [43] was used for water molecules. The comparison with the experimental data [60] shows the quantitative agreement between the experimental and the calculated activity coefficients in the wide range of salt concentrations. 3.2.2. Solvation structure and selectivity of GLIC ion channel In this section, we employ the 3D-RISM-KH molecular theory of solvation to study microscopic solvation structure and selectivity of the pentameric ligand-gated ion channel embedded in the cellular membrane. The structure of the bacterial G. violaceus pentameric ligand-gated ion channel homologue (GLIC) in an open conformation was solved in the X-ray experiments [35]. It was also demonstrated that this channel is sensitive to general anesthetics [61]. Recently, the crystal structures of GLIC in the complex with propofol and desflurane ligands were identified in X-ray experiments [62]. Here, we present the results of modeling the microscopic solvation structure of the channel without ligands to demonstrate the capability of the 3D-RISM-KH method to predict selectivity and permeation properties of ion channels. The initial coordinates of the channel protein were taken from the Protein Data Bank (PDB accession code 3EAM [35]). Similar to the previous theoretical modeling of the system [53,62], the protein was inserted into 1-palmitoyl-2-oleoyl-sn-glycero-3phosphocholine (POPC) lipid bilayer with the help of VMD software [63] (Fig. 11). Parameters of the solute (the channel protein and the lipid bilayer)–solvent interaction potential needed as an input for 3D-RISM calculations were taken from the Amber ff03 force field [39,40]. Equations of the 3D molecular theory of solvation (1) and (2) were solved for GLIC–membrane complex in sodium chloride water solution at salt concentration of 0.208 mol/kg for the energy minimized structure of the complex. Minimization was performed with the Amber 9 molecular dynamics package [41,42], the all-atom ff03 force field [39,40] and the implicit generalized Born solvation model. The 3D density distributions for the oxygen site of water as well as for chloride and sodium ions around the channel protein embedded into cellular membrane are represented by isosurfaces in the left and right panels of Fig. 12. It is clearly seen from the figure that the negative chloride ions are mostly expelled from the pore of this

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

111

characterized by the enhanced densities of the anions. The sodium distribution inside the channel pore correlates with the channel radius and the pattern in the distributions of the hydrophobic and hydrophilic residues in the channel pore. Our results show how the charge distribution along the pore and its geometry contribute to the selectivity and permeation properties of the ion channel. To the best of our knowledge, the above constitutes the first application of the statistical–mechanical integral equation theory of molecular liquids to ion channels, although the method was previously used in the similar context to study internal solvation and permeation properties of the aquaporin channels [57,30].

4. Conclusions and summary

Fig. 11. Schematic representation of the GLIC channel (PDB accession code 3EAM [35]) in the POPC lipid bilayer. The extracellular and intracellular spaces are above and below the lipid membrane, respectively. The figure was prepared with DiscoveryStudio 2.0 (Accelrys Inc.) software.

cationic channel. In contrast, there are non-zero densities of sodium ions throughout the channel pore. This illustrates the capability of the 3D-RISM-KH method to predict selectivity of ionic channels in agreement with the experimental data. The solvent distributions in the proximity of the protein are characterized by some interesting features. For example, the chloride ions' preferable solvation sites are located along the external (facing the extracellular space above the membrane) ring of the channel “funnel”. Inside the funnel, right before the entrance to the channel pore, there is complementarity between the distributions of water and sodium ions. There is a pocket of chloride ions which separates the upper and lower parts of the funnel (the area above the level of the lipid bilayer on the left panel of Fig. 12). The lower part, being adjacent to the channel pore, is characterized by the preferable solvation of sodium ions as much as the rest of the channel pore (the right panel of Fig. 12). The exit from the channel pore from the intracellular side (below the membrane) is surrounded by the solvation sites of chloride ions. In general, both the extra and intracellular surfaces of the lipid bilayer are

(A)

The three-dimensional molecular theory of solvation (also known as 3D-RISM-KH [12]) is based on the rigorous statistical–mechanical formulation of integral theory of molecular liquids [36]. The theory provides a reliable basis for studying nanoscale phenomena in complex environments, and can be used as an integral part of multiscale methods for nanochemistry and biophysics in solution. 3D-RISM-KH describes solvent based on the statistical–mechanical distribution functions, thus avoiding conformational sampling of the solvent degrees of freedom. The theory is fully transferable and can be used for molecular liquids of different compositions and complexities. It yields an analytical representation for the solvation free energy in terms of the solute–solvent correlation functions which allows one to skip thermodynamic integration in calculations of the solvation free energy. In this paper, we give a brief overview of the 3D molecular theory of solvation and the multiscale method that couples the 3D-RISM-KH theory treating the solvent with molecular dynamics applied to the biomolecule and is implemented in the Amber molecular dynamics package. We postpone application of the MD/3D-RISM-KH multiscale method for biomolecular and other nanoscale systems to our future publications. For illustration of the capability of the 3D-RISM-KH theory to deal with problems in complex biophysical systems, the MM/3D-RISM-KH post-processing method was used to study the binding modes and the free energy of binding of thiamine (vitamin of group B) against extracytoplasmic thiamine binding lipoprotein MG289 [34]. We show how 3D-RISM-KH predicts the binding maps,

(B)

Fig. 12. Ions and water distributions in the proximity of the GLIC channel. The channel protein is shown in schematic representation with β-sheet (in yellow) and α-helices (in magenta). The lipid monomers shown in the ball-and-stick representation. The orange, blue and gray colors are used to represent density isosurfaces corresponding to chloride ions, sodium ions, and water oxygen sites, respectively. Panel (A): chloride ions and water distributions in the proximity of GLIC embedded in the cellular membrane. Panel (B): distributions of sodium ions and water around the channel protein and membrane.

112

A. Kovalenko, N. Blinov / Journal of Molecular Liquids 164 (2011) 101–112

binding affinity and microscopic solvation structure of the thiamine– MG289 complex. Efficient contraction of solvent degrees of freedom allows one to use 3D-RISM-KH for a statistically meaningful description of slow process and rare effects such as, for example, permeation of ions and small molecules through channel proteins [57]. As a case study, we have analyzed the solvation structure, permeation and selectivity of the bacterial G. violaceus pentameric ligand-gated ion channel (GLIC) homologue [35]. We demonstrate that the 3D-RISM-KH method is capable of predicting selectivity of ionic channels. The results we obtained are in agreement with experimental data. Acknowledgments This work was supported by the Alberta Prion Research Institute (APRI) and the National Research Council (NRC) of Canada. References [1] M. Karplus, Accounts of Chemical Research 35 (6) (2002) 321–323. [2] W.C. Still, A. Tempczyk, R.C. Hawley, T. Hendrickson, Journal of the American Chemical Society 112 (16) (1990) 6127–6129. [3] V. Tsui, D. Case, Journal of the American Chemical Society 122 (11) (2000) 2489–2498. [4] J. Srinivasan, M.W. Trevathan, P. Beroza, D.A. Case, Theor. Chim. Acta 101 (1999) 426–434. [5] H. Gohlke, D.A. Case, Journal of Computational Chemistry 25 (2) (2004) 238–250. [6] A. Onufriev, D. Bashford, D.A. Case, Proteins: Structure, Function, and Bioinformatics 55 (2) (2004) 383–394. [7] M. Feig, A. Onufriev, M.S. Lee, W. Im, D.A. Case, C.L. Brooks III, Journal of Computational Chemistry 25 (2) (2004) 265–284. [8] E. Gallicchio, Ronald M. Levy, Linda Y. Zhang, A.K. Felts, Journal of the American Chemical Society 125 (2003) 9523–9530. [9] J.A. Wagoner, N.A. Baker, Proceedings of the National Academy of Sciences of the United States of America 103 (22) (2006) 8331–8336. [10] D. Beglov, B. Roux, The Journal of Physical Chemistry B 101 (39) (1997) 7821–7826. [11] A. Kovalenko, F. Hirata, Journal of Chemical Physics 110 (20) (1999) 10095–10112. [12] A. Kovalenko, in: F. Hirata (Ed.), Molecular Theory of Solvation, Understanding Chemical Reactivity, Vol. 24, Kluwer Academic Publishers, Dordrecht, 2003, pp. 169–275. [13] J.G. Moralez, J. Raez, T. Yamazaki, R.K. Motkuri, A. Kovalenko, H. Fenniri, Journal of the American Chemical Society 127 (23) (2005) 8307–8309. [14] R.S. Johnson, T. Yamazaki, A. Kovalenko, H. Fenniri, Journal of the American Chemical Society 129 (17) (2007) 5735–5743. [15] G. Tikhomirov, T. Yamazaki, A. Kovalenko, H. Fenniri, Langmuir 24 (9) (2008) 4447–4450. [16] T. Yamazaki, N. Blinov, D. Wishart, A. Kovalenko, Biophysical Journal 95 (10) (2008) 4540–4548. [17] T. Imai, S. Ohyama, A. Kovalenko, F. Hirata, Protein Science 16 (9) (2007) 1927–1933. [18] T. Imai, R. Hiraoka, A. Kovalenko, F. Hirata, Proteins: Structure, Function, and Bioinformatics 66 (4) (2007) 804–813. [19] T. Imai, R. Hiraoka, T. Seto, A. Kovalenko, F. Hirata, The Journal of Physical Chemistry B 111 (39) (2007) 11585–11591. [20] T. Imai, Y. Harano, M. Kinoshita, A. Kovalenko, F. Hirata, Journal of Chemical Physics 126 (22) (2007) 225102. [21] T. Imai, Y. Harano, M. Kinoshita, A. Kovalenko, F. Hirata, Journal of Chemical Physics 125 (2) (2006) 024911. [22] Y. Harano, T. Imai, A. Kovalenko, M. Kinoshita, F. Hirata, Journal of Chemical Physics 114 (21) (2001) 9506–9511. [23] P. Drabik, S. Gusarov, A. Kovalenko, Biophysical Journal 92 (2) (2007) 394–403. [24] T. Imai, H. Isogai, T. Seto, A. Kovalenko, F. Hirata, The Journal of Physical Chemistry B 110 (24) (2006) 12149–12154. [25] T. Imai, A. Kovalenko, F. Hirata, Molecular Simulation 32 (10&11) (2006) 817–824. [26] T. Imai, R. Hiraoka, A. Kovalenko, F. Hirata, Journal of the American Chemical Society 127 (44) (2005) 15334–15335. [27] N. Blinov, L. Dorosh, D. Wishart, A. Kovalenko, Biophysical Journal 98 (2) (2010) 282–296. [28] M.C. Stumpe, N. Blinov, D. Wishart, A. Kovalenko, V.S. Pande, The Journal of Physical Chemistry B 115 (2) (2011) 319–328.

[29] T. Luchko, S. Gusarov, D.R. Roe, C. Simmerling, D.A. Case, J. Tuszynski, A. Kovalenko, Journal of Chemical Theory and Computation 6 (3) (2010) 607–624. [30] N. Yoshida, T. Imai, S. Phongphanphanee, A. Kovalenko, F. Hirata, The Journal of Physical Chemistry B 113 (4) (2009) 873–886. [31] T. Imai, K. Oda, A. Kovalenko, F. Hirata, A. Kidera, Journal of the American Chemical Society 131 (34) (2009) 12430–12440. [32] Y. Kiyota, N. Yoshida, F. Hirata, Journal of Molecular Liquids 159 (1) (2011) 93–98. [33] N. Blinov, L. Dorosh, D. Wishart, A. Kovalenko, Molec. Simul. 37 (2011) 718–728. [34] K.H. Sippel, B. Venkatakrishnan, S.K. Boehlein, B. Sankaran, J.G. Quirit, L. Govindasamy, M. Agbandje-McKenna, S. Goodison, C.J. Rosser, R. McKenna, Proteins: Structure, Function, and Bioinformatics 79 (2) (2011) 528–536. [35] N. Bocquet, H. Nury, M. Baaden, C. Le Poupon, J.-P. Changeux, M. Delarue, P.-J. Corringer, Nature 457 (2009) 111–114. [36] J.P. Hansen, I.R. McDonald, Theory of Simple Liquids, Academic Press, London, 1986. [37] J.S. Perkyns, B.M. Pettitt, Journal of Chemical Physics 97 (1992) 7656–7666. [38] T. Yamazaki, A. Kovalenko, Journal of Chemical Theory and Computation 5 (7) (2009) 1723–1730. [39] Y. Duan, C. Wu, S. Chowdhury, M.C. Lee, G.X.W. Zhang, R. Yang, P. Cieplak, R. Luo, T. Lee, J. Caldwell, J. Wang, P. Kollman, Journal of Computational Chemistry 24 (16) (2003) 1999–2012. [40] M.C. Lee, Y. Duan, Proteins: Structure, Function, and Bioinformatics 55 (3) (2004) 620–634. [41] D.A. Case, T.E. Cheatham III, T. Darden, H. Gohlke, R. Luo, K.M. Merz Jr., A. Onufriev, C. Simmerling, B. Wang, R.J. Woods, Journal of Computational Chemistry 26 (16) (2005) 1668–1688. [42] D.A. Case, T. Darden, T.E. Cheatham III, C. Simmerling, J. Wang, R.E. Duke, R. Luo, K.M. Merz, D.A. Pearlman, M. Crowley, R.C. Walker, W. Zhang, B. Wang, S. Hayik, A. Roitberg, G. Seabra, K.F. Wong, F. Paesani, X. Wu, S. Brozell, V. Tsui, H. Gohlke, L. Yang, C. Tan, J. Mongan, V. Hornak, G. Cui, P. Beroza, D.H. Mathews, C. Schafmeister, W.S. Ross, P.A. Kollman, AMBER 9, University of California, San Francisco, 2006. [43] H.J.C. Berendsen, J.R. Grigera, T.P. Straatsma, Journal of Physical Chemistry 91 (24) (1987) 6269–6271. [44] M.P. Allen, D.J. Tildesley, Computer Simulation of Liquids, Clarendon Press, New York, NY, USA, 1989. [45] S. Genheden, T. Luchko, S. Gusarov, A. Kovalenko, U. Ryde, The Journal of Physical Chemistry B 114 (25) (2010) 8505–8516. [46] G.M. Morris, R. Huey, W. Lindstrom, M.F. Sanner, R.K. Belew, D.S. Goodsell, A.J. Olson, Journal of Computational Chemistry 30 (16) (2009) 2785–2791. [47] J. Srinivasan, T. Cheatham, P. Cieplak, P. Kollman, D. Case, Journal of the American Chemical Society 120 (37) (1998) 9401–9409. [48] M.R. Lee, Y. Duan, P.A. Kollman, Proteins: Structure, Function, and Genetics 39 (2000) 309–316. [49] A. Kovalenko, F. Hirata, Journal of Chemical Physics 112 (23) (2000) 10391–10402. [50] D. Lucent, V. Vishal, V.S. Pande, Proceedings of the National Academy of Sciences of the United States of America 104 (25) (2007) 10430–10434. [51] M.S.P. Sansom, I.H. Shrivastava, K.M. Ranatunga, G.R. Smith, Trends in Biochemical Sciences 25 (8) (2000) 368–374. [52] S.Y. Noskov, S. Berneche, B. Roux, Nature 431 (2004) 830–834. [53] H. Nury, F. Poitevin, C. Van Renterghem, J.-P. Changeux, P.-J. Corringer, Delarue, M. Baaden, Proc. Natl. Acad. Sci. 107 (2010) 6275–6280. [54] D. Boda, W. Nonner, M. Valisko, D. Henderson, B. Eisenberg, D. Gillespie, Biophysical Journal 93 (6) (2007) 1960–1980. [55] D. Boda, M. Valiskó, B. Eisenberg, W. Nonner, D. Henderson, D. Gillespie, The Journal of Chemical Physics 125 (3) (2006) 034901. [56] D. Boda, M. Valiskó, B. Eisenberg, W. Nonner, D. Henderson, D. Gillespie, Physical Review Letters 98 (16) (2007) 168102. [57] S. Phongphanphanee, N. Yoshida, F. Hirata, The Journal of Physical Chemistry B 114 (23) (2010) 7967–7973. [58] N. Yoshida, S. Phongphanphanee, F. Hirata, The Journal of Physical Chemistry B 111 (17) (2007) 4588–4595. [59] W.L. Jorgensen, J. Tirado-Rives, Journal of the American Chemical Society 110 (6) (1988) 1657–1666. [60] W.J. Moore, Physical Chemistry, 5th edition Longman Publishing Group, March 1998. [61] Y. Weng, L. Yang, P.-J. Corringer, J.M. Sonner, Anesthesia and Analgesia 110 (2010) 59–63. [62] H. Nury, C. Van Renterghem, Y. Weng, A. Tran, M. Baaden, V. Dufresne, J.-P. Changeux, J.M. Sonner, M. Delarue, P.-J. Corringer, Nature 469 (2011) 428–431. [63] W. Humphrey, A. Dalke, K. Schulten, Journal of Molecular Graphics 14 (1) (1996) 33–38.