Implicit solvation (sometimes termed continuum solvation) is a method to represent solvent as a continuous medium instead of individual “explicit” solvent molecules, most often used in molecular dynamics simulations and in other applications of molecular mechanics. The method is often applied to estimate free energy of solute-solvent interactions in structural and chemical processes, such as folding or conformational transitions of proteins, DNA, RNA, and polysaccharides, association of biological macromolecules with ligands, or transport of drugs across biological membranes.
The implicit solvation model is justified in liquids, where the potential of mean force can be applied to approximate the averaged behavior of many highly dynamic solvent molecules. However, the interfaces and the interiors of biological membranes or proteins can also be considered as media with specific solvation or dielectric properties. These media are not necessarily uniform, since their properties can be described by different analytical functions, such as “polarity profiles” of lipid bilayers.
There are two basic types of implicit solvent methods: models based on accessible surface areas (ASA) that were historically the first, and more recent continuum electrostatics models, although various modifications and combinations of the different methods are possible. The accessible surface area (ASA) method is based on experimental linear relations between Gibbs free energy of transfer and the surface area of a solute molecule. This method operates directly with free energy of solvation, unlike molecular mechanics or electrostatic methods that include only the enthalpic component of free energy. The continuum representation of solvent also significantly improves the computational speed and reduces errors in statistical averaging that arise from incomplete sampling of solvent conformations, so that the energy landscapes obtained with implicit and explicit solvent are different. Although the implicit solvent model is useful for simulations of biomolecules, this is an approximate method with certain limitations and problems related to parameterization and treatment of ionization effects.
Accessible surface area-based method
The free energy of solvation of a solute molecule in the simplest ASA-based method is given by:
where is the accessible surface area of atom i, and is solvation parameter of atom i, i.e., a contribution to the free energy of solvation of the particular atom i per surface unit area. The needed solvation parameters for different types of atoms (carbon (C), nitrogen (N), oxygen (O), sulfur (S), etc.) are usually determined by a least squares fit of the calculated and experimental transfer free energies for a series of organic compounds. The experimental energies are determined from partition coefficients of these compounds between different solutions or media using standard mole concentrations of the solutes.
Notably, solvation energy is the free energy needed to transfer a solute molecule from a solvent to vacuum (gas phase). This energy can supplement the intramolecular energy in vacuum calculated in molecular mechanics. Thus, the needed atomic solvation parameters were initially derived from water-gas partition data. However, the dielectric properties of proteins and lipid bilayers are much more similar to those of nonpolar solvents than to vacuum. Newer parameters have thus been derived from octanol-water partition coefficients or other similar data. Such parameters actually describe transfer energy between two condensed media or the difference of two solvation energies.
or (in mks):
where represents the position-dependent dielectric, represents the electrostatic potential, represents the charge density of the solute, represents the concentration of the ion i at a distance of infinity from the solute, is the valence of the ion, q is the charge of a proton, k is the Boltzmann constant, T is the temperature, and is a factor for the position-dependent accessibility of position r to the ions in solution (often set to uniformly 1). If the potential is not large, the equation can be linearized to be solved more efficiently.
Although this equation has solid theoretical justification, it is computationally expensive to calculate without approximations. A number of numerical Poisson-Boltzmann equation solvers of varying generality and efficiency have been developed, including one application with a specialized computer hardware platform. However, performance from PB solvers does not yet equal that from the more commonly used generalized Born approximation.
Generalized Born model
The Generalized Born (GB) model is an approximation to the exact (linearized) Poisson-Boltzmann equation. It is based on modeling the solute as a set of spheres whose internal dielectric constant differs from the external solvent. The model has the following functional form:
where is the permittivity of free space, is the dielectric constant of the solvent being modeled, is the electrostatic charge on particle i, is the distance between particles i and j, and is a quantity (with the dimension of length) termed the effective Born radius. The effective Born radius of an atom characterizes its degree of burial inside the solute; qualitatively it can be thought of as the distance from the atom to the molecular surface. Accurate estimation of the effective Born radii is critical for the GB model.
With accessible surface area
The Generalized Born (GB) model augmented with the hydrophobic solvent accessible surface area (SA) term is GBSA. It is among the most commonly used implicit solvent model combinations. The use of this model in the context of molecular mechanics is termed MM/GBSA. Although this formulation has been shown to successfully identify the native states of short peptides with well-defined tertiary structure, the conformational ensembles produced by GBSA models in other studies differ significantly from those produced by explicit solvent and do not identify the protein's native state. In particular, salt bridges are overstabilized, possibly due to insufficient electrostatic screening, and a higher-than-native alpha helix population was observed. Variants of the GB model have also been developed to approximate the electrostatic environment of membranes, which have had some success in folding the transmembrane helixes of integral membrane proteins.
Ad hoc fast solvation models
Another possibility is to use ad hoc quick strategies to estimate solvation free energy. A first generation of fast implicit solvents is based on the calculation of a per-atom solvent accessible surface area. For each of group of atom types, a different parameter scales its contribution to solvation ("ASA-based model" described above).
The reference solvation free energy of i corresponds to a suitably chosen small molecule in which group i is essentially fully solvent-exposed. The integral is over the volume Vj of group j and the summation is over all groups j around i. EEF1 additionally uses a distance-dependent (non-constant) dielectric, and ionic side-chains of proteins are simply neutralized. It is only 50% slower than a vacuum simulation. This model was later augmented with the hydrophobic effect and called Charmm19/SASA.
Hybrid implicit-explicit solvation models
It is possible to include a layer or sphere of water molecules around the solute, and model the bulk with an implicit solvent. Such an approach is proposed by M. J. Frisch and coworkers and by other authors. For instance in Ref. the bulk solvent is modeled with a Generalized Born approach and the multi-grid method used for Coulombic pairwise particle interactions. It is reported to be faster than a full explicit solvent simulation with the particle mesh Ewald summation (PME) method of electrostatic calculation. There are a range of hybrid methods available capable of accessing and acquiring information on solvation.
Effects unaccounted for
The hydrophobic effect
Models like PB and GB allow estimation of the mean electrostatic free energy but do not account for the (mostly) entropic effects arising from solute-imposed constraints on the organization of the water or solvent molecules. This is termed the hydrophobic effect and is a major factor in the folding process of globular proteins with hydrophobic cores. Implicit solvation models may be augmented with a term that accounts for the hydrophobic effect. The most popular way to do this is by taking the solvent accessible surface area (SASA) as a proxy of the extent of the hydrophobic effect. Most authors place the extent of this effect between 5 and 45 cal/(Å2 mol). Note that this surface area pertains to the solute, while the hydrophobic effect is mostly entropic in nature at physiological temperatures and occurs on the side of the solvent.
Implicit solvent models such as PB, GB, and SASA lack the viscosity that water molecules impart by randomly colliding and impeding the motion of solutes through their van der Waals repulsion. In many cases, this is desirable because it makes sampling of configurations and phase space much faster. This acceleration means that more configurations are visited per simulated time unit, on top of whatever CPU acceleration is achieved in comparison to explicit solvent. It can, however, lead to misleading results when kinetics are of interest.
Viscosity may be added back by using Langevin dynamics instead of Hamiltonian mechanics and choosing an appropriate damping constant for the particular solvent. In practical bimolecular simulations one can often speed-up conformational search significantly (up to 100 times in some cases) by using much lower collision frequency . Recent work has also been done developing thermostats based on fluctuating hydrodynamics to account for momentum transfer through the solvent and related thermal fluctuations. One should keep in mind, though, that the folding rate of proteins does not depend linearly on viscosity for all regimes.
Hydrogen bonds with solvent
Solute-solvent hydrogen bonds in the first solvation shell are important for solubility of organic molecules and especially ions. Their average energetic contribution can be reproduced with an implicit solvent model.
Problems and limitations
All implicit solvation models rest on the simple idea that nonpolar atoms of a solute tend to cluster together or occupy nonpolar media, whereas polar and charged groups of the solute tend to remain in water. However, it is important to properly balance the opposite energy contributions from different types of atoms. Several important points have been discussed and investigated over the years.
Choice of model solvent
It has been noted that wet 1-octanol solution is a poor approximation of proteins or biological membranes because it contains ~2M of water, and that cyclohexane would be a much better approximation. Investigation of passive permeability barriers for different compounds across lipid bilayers led to conclusion that 1,9-decadiene can serve as a good approximations of the bilayer interior, whereas 1-octanol was a very poor approximation. A set of solvation parameters derived for protein interior from protein engineering data was also different from octanol scale: it was close to cyclohexane scale for nonpolar atoms but intermediate between cyclohexane and octanol scales for polar atoms. Thus, different atomic solvation parameters should be applied for modeling of protein folding and protein-membrane binding. This issue remains controversial. The original idea of the method was to derive all solvation parameters directly from experimental partition coefficients of organic molecules, which allows calculation of solvation free energy. However, some of the recently developed electrostatic models use ad hoc values of 20 or 40 cal/(Å2 mol) for all types of atoms. The non-existent “hydrophobic” interactions of polar atoms are overridden by large electrostatic energy penalties in such models.
Strictly speaking, ASA-based models should only be applied to describe solvation, i.e., energetics of transfer between liquid or uniform media. It is possible to express van der Waals interaction energies in the solid state in the surface energy units. This was sometimes done for interpreting protein engineering and ligand binding energetics, which leads to “solvation” parameter for aliphatic carbon of ~40 cal/(Å2 mol), which is 2 times bigger than ~20 cal/(Å2 mol) obtained for transfer from water to liquid hydrocarbons, because the parameters derived by such fitting represent sum of the hydrophobic energy (i.e., 20 cal/Å2 mol) and energy of van der Waals attractions of aliphatic groups in the solid state, which corresponds to fusion enthalpy of alkanes. Unfortunately, the simplified ASA-based model cannot capture the "specific" distance-dependent interactions between different types of atoms in the solid state which are responsible for clustering of atoms with similar polarities in protein structures and molecular crystals. Parameters of such interatomic interactions, together with atomic solvation parameters for the protein interior, have been approximately derived from protein engineering data. The implicit solvation model breaks down when solvent molecules associate strongly with binding cavities in a protein, so that the protein and the solvent molecules form a continuous solid body. On the other hand, this model can be successfully applied for describing transfer from water to the fluid lipid bilayer.
Importance of extensive testing
More testing is needed to evaluate the performance of different implicit solvation models and parameter sets. They are often tested only for a small set of molecules with very simple structure, such as hydrophobic and amphiphilic alpha helixes (α). This method was rarely tested for hundreds of protein structures.
Treatment of ionization effects
Ionization of charged groups has been neglected in continuum electrostatic models of implicit solvation, as well as in standard molecular mechanics and molecular dynamics. The transfer of an ion from water to a nonpolar medium with dielectric constant of ~3 (lipid bilayer) or 4 to 10 (interior of proteins) costs significant energy, as follows from the Born equation and from experiments. However, since the charged protein residues are ionizable, they simply lose their charges in the nonpolar environment, which costs relatively little at the neutral pH: ~4 to 7 kcal/mol for Asp, Glu, Lys, and Arg amino acid residues, according to the Henderson-Hasselbalch equation, ΔG = 2.3RT (pH - pK). The low energetic costs of such ionization effects have indeed been observed for protein mutants with buried ionizable residues. and hydrophobic α-helical peptides in membranes with a single ionizable residue in the middle. However, all electrostatic methods, such as PB, GB, or GBSA assume that ionizable groups remain charged in the nonpolar environments, which leads to grossly overestimated electrostatic energy. In the simplest accessible surface area-based models, this problem was treated using different solvation parameters for charged atoms or Henderson-Hasselbalch equation with some modifications. However even the latter approach does not solve the problem. Charged residues can remain charged even in the nonpolar environment if they are involved in intramolecular ion pairs and H-bonds. Thus, the energetic penalties can be overestimated even using the Henderson-Hasselbalch equation. More rigorous theoretical methods describing such ionization effects have been developed, and there are ongoing efforts to incorporate such methods into the implicit solvation models.
- Marsh D (Jul 2001). "Polarity and permeation profiles in lipid membranes". Proceedings of the National Academy of Sciences of the United States of America. 98 (14): 7777–82. Bibcode:2001PNAS...98.7777M. doi:10.1073/pnas.131023798. PMC 35418. PMID 11438731.
- Richards FM (1977). "Areas, volumes, packing and protein structure". Annual Review of Biophysics and Bioengineering. 6: 151–76. doi:10.1146/annurev.bb.06.060177.001055. PMID 326146.
- Roux B, Simonson T (Apr 1999). "Implicit solvent models". Biophysical Chemistry. 78 (1–2): 1–20. doi:10.1016/S0301-4622(98)00226-9. PMID 17030302.
- Zhou R (Nov 2003). "Free energy landscape of protein folding in water: explicit vs. implicit solvent". Proteins. 53 (2): 148–61. doi:10.1002/prot.10483. PMID 14517967. S2CID 22977210.
- Ben-Naim AY (1980). Hydrophobic interactions. New York: Plenum Press. ISBN 978-0-306-40222-7.
- Holtzer A (Jun 1995). "The "cratic correction" and related fallacies" (Free full text). Biopolymers. 35 (6): 595–602. doi:10.1002/bip.360350605. PMID 7766825. S2CID 45118648.
- Ooi T, Oobatake M, Némethy G, Scheraga HA (May 1987). "Accessible surface areas as a measure of the thermodynamic parameters of hydration of peptides" (Free full text). Proceedings of the National Academy of Sciences of the United States of America. 84 (10): 3086–90. Bibcode:1987PNAS...84.3086O. doi:10.1073/pnas.84.10.3086. PMC 304812. PMID 3472198.
- Eisenberg D, McLachlan AD (Jan 1986). "Solvation energy in protein folding and binding". Nature. 319 (6050): 199–203. Bibcode:1986Natur.319..199E. doi:10.1038/319199a0. PMID 3945310. S2CID 21867582.
- Fogolari F, Brigo A, Molinari H (Nov 2002). "The Poisson-Boltzmann equation for biomolecular electrostatics: a tool for structural biology". Journal of Molecular Recognition. 15 (6): 377–92. doi:10.1002/jmr.577. PMID 12501158. S2CID 17184352.
- Shestakov AI, Milovich JL, Noy A (Mar 2002). "Solution of the nonlinear Poisson-Boltzmann equation using pseudo-transient continuation and the finite element method". Journal of Colloid and Interface Science. 247 (1): 62–79. Bibcode:2002JCIS..247...62S. doi:10.1006/jcis.2001.8033. PMID 16290441.
- Lu B, Zhang D, McCammon JA (Jun 2005). "Computation of electrostatic forces between solvated molecules determined by the Poisson-Boltzmann equation using a boundary element method" (PDF). The Journal of Chemical Physics. 122 (21): 214102. Bibcode:2005JChPh.122u4102L. doi:10.1063/1.1924448. PMID 15974723. S2CID 23731263.
- Baker NA, Sept D, Joseph S, Holst MJ, McCammon JA (Aug 2001). "Electrostatics of nanosystems: application to microtubules and the ribosome". Proceedings of the National Academy of Sciences of the United States of America. 98 (18): 10037–41. Bibcode:2001PNAS...9810037B. doi:10.1073/pnas.181342398. PMC 56910. PMID 11517324.
- Höfinger S (Aug 2005). "Solving the Poisson-Boltzmann equation with the specialized computer chip MD-GRAPE-2". Journal of Computational Chemistry. 26 (11): 1148–54. doi:10.1002/jcc.20250. PMID 15942918. S2CID 19378083.
- Koehl P (Apr 2006). "Electrostatics calculations: latest methodological advances". Current Opinion in Structural Biology. 16 (2): 142–51. doi:10.1016/j.sbi.2006.03.001. PMID 16540310.
- Still WC, Tempczyk A, Hawley RC, Hendrickson T (1990). "Semianalytical treatment of solvation for molecular mechanics and dynamics". J Am Chem Soc. 112 (16): 6127–6129. doi:10.1021/ja00172a038.
- Onufriev A, Case DA, Bashford D (Nov 2002). "Effective Born radii in the generalized Born approximation: the importance of being perfect". Journal of Computational Chemistry. 23 (14): 1297–304. CiteSeerX 10.1.1.107.962. doi:10.1002/jcc.10126. PMID 12214312. S2CID 12244353.
- Ho BK, Dill KA (Apr 2006). "Folding very short peptides using molecular dynamics". PLOS Computational Biology. 2 (4): e27. Bibcode:2006PLSCB...2...27H. doi:10.1371/journal.pcbi.0020027. PMC 1435986. PMID 16617376.
- Im W, Feig M, Brooks CL (Nov 2003). "An implicit membrane generalized born theory for the study of structure, stability, and interactions of membrane proteins". Biophysical Journal. 85 (5): 2900–18. Bibcode:2003BpJ....85.2900I. doi:10.1016/S0006-3495(03)74712-2. PMC 1303570. PMID 14581194.
- Wesson L, Eisenberg D (Feb 1992). "Atomic solvation parameters applied to molecular dynamics of proteins in solution" (Free full text). Protein Science. 1 (2): 227–35. doi:10.1002/pro.5560010204. PMC 2142195. PMID 1304905.
- Lazaridis T, Karplus M (May 1999). "Effective energy function for proteins in solution". Proteins. 35 (2): 133–52. doi:10.1002/(SICI)1097-0134(19990501)35:2<133::AID-PROT1>3.0.CO;2-N. PMID 10223287.
- Ferrara P, Apostolakis J, Caflisch A (Jan 2002). "Evaluation of a fast implicit solvent model for molecular dynamics simulations". Proteins. 46 (1): 24–33. CiteSeerX 10.1.1.25.1195. doi:10.1002/prot.10001. PMID 11746700. S2CID 17671699.
- TA Keith, MJ Frisch (1994). "Chapter 3: Inclusion of Explicit Solvent Molecules in a Self-Consistent-Reaction Field Model of Solvation". In Smith D (ed.). Modeling the hydrogen bond. Columbus, OH: American Chemical Society. ISBN 978-0-8412-2981-5.
- Lee MS, Salsbury FR, Olson MA (Dec 2004). "An efficient hybrid explicit/implicit solvent method for biomolecular simulations". Journal of Computational Chemistry. 25 (16): 1967–78. doi:10.1002/jcc.20119. PMID 15470756. S2CID 16996683.
- Marini A, Muñoz-Losa A, Biancardi A, Mennucci B (Dec 2010). "What is solvatochromism?". The Journal of Physical Chemistry B. 114 (51): 17128–35. doi:10.1021/jp1097487. PMID 21128657.
- Skyner RE, McDonagh JL, Groom CR, van Mourik T, Mitchell JB (Mar 2015). "A review of methods for the calculation of solution free energies and the modelling of systems in solution". Physical Chemistry Chemical Physics. 17 (9): 6174–91. Bibcode:2015PCCP...17.6174S. doi:10.1039/C5CP00288E. PMID 25660403.
- Sharp KA, Nicholls A, Fine RF, Honig B (Apr 1991). "Reconciling the magnitude of the microscopic and macroscopic hydrophobic effects". Science. 252 (5002): 106–9. Bibcode:1991Sci...252..106S. doi:10.1126/science.2011744. PMID 2011744.
- Schlick T (2002). Molecular Modeling and Simulation: An Interdisciplinary Guide Interdisciplinary Applied Mathematics: Mathematical Biology. New York: Springer. ISBN 978-0-387-95404-2.
- Anandakrishnan R, Drozdetski A, Walker RC, Onufriev AV (March 2015). "Speed of conformational change: comparing explicit and implicit solvent molecular dynamics simulations". Biophysical Journal. 108 (5): 1153–64. Bibcode:2015BpJ...108.1153A. doi:10.1016/j.bpj.2014.12.047. PMC 4375717. PMID 25762327.
- Wang Y, Sigurdsson JK, Brandt E, Atzberger PJ (August 2013). "Dynamic implicit-solvent coarse-grained models of lipid bilayer membranes: fluctuating hydrodynamics thermostat". Physical Review E. 88 (2): 023301. arXiv:1212.0449. Bibcode:2013PhRvE..88b3301W. doi:10.1103/PhysRevE.88.023301. PMID 24032960. S2CID 118409341.
- Zagrovic B, Pande V (September 2003). "Solvent viscosity dependence of the folding rate of a small protein: distributed computing study". Journal of Computational Chemistry. 24 (12): 1432–6. doi:10.1002/jcc.10297. PMID 12868108. S2CID 16666252.
- Lomize AL, Pogozheva ID, Mosberg HI (Apr 2011). "Anisotropic solvent model of the lipid bilayer. 1. Parameterization of long-range electrostatics and first solvation shell effects". Journal of Chemical Information and Modeling. 51 (4): 918–29. doi:10.1021/ci2000192. PMC 3089899. PMID 21438609.
- Lomize AL, Pogozheva ID, Mosberg HI (Apr 2011). "Anisotropic solvent model of the lipid bilayer. 2. Energetics of insertion of small molecules, peptides, and proteins in membranes". Journal of Chemical Information and Modeling. 51 (4): 930–46. doi:10.1021/ci200020k. PMC 3091260. PMID 21438606.
- Radzicka A, Wolfenden R (1988). "Comparing the polarities of the amino acids: side-chain distribution coefficients between the vapor phase, cyclohexane, 1-octanol, and neutral aqueous solution". Biochemistry. 27 (5): 1664–1670. doi:10.1021/bi00405a042.
- Mayer PT, Anderson BD (Mar 2002). "Transport across 1,9-decadiene precisely mimics the chemical selectivity of the barrier domain in egg lecithin bilayers". Journal of Pharmaceutical Sciences. 91 (3): 640–6. doi:10.1002/jps.10067. PMID 11920749.
- Walter A, Gutknecht J (1986). "Permeability of small nonelectrolytes through lipid bilayer membranes". The Journal of Membrane Biology. 90 (3): 207–17. doi:10.1007/BF01870127. PMID 3735402. S2CID 23972379.
- Lomize AL, Reibarkh MY, Pogozheva ID (Aug 2002). "Interatomic potentials and solvation parameters from protein engineering data for buried residues" (Free full text). Protein Science. 11 (8): 1984–2000. doi:10.1110/ps.0307002. PMC 2373680. PMID 12142453.
- Eriksson AE, Baase WA, Zhang XJ, Heinz DW, Blaber M, Baldwin EP, Matthews BW (Jan 1992). "Response of a protein structure to cavity-creating mutations and its relation to the hydrophobic effect". Science. 255 (5041): 178–83. Bibcode:1992Sci...255..178E. CiteSeerX 10.1.1.461.7843. doi:10.1126/science.1553543. PMID 1553543. S2CID 12278880.
- Funahashi J, Takano K, Yutani K (Feb 2001). "Are the parameters of various stabilization factors estimated from mutant human lysozymes compatible with other proteins?" (Free full text). Protein Engineering. 14 (2): 127–34. doi:10.1093/protein/14.2.127. PMID 11297670.
- Lomize AL, Pogozheva ID, Mosberg HI (Oct 2004). "Quantification of helix-helix binding affinities in micelles and lipid bilayers" (Free full text). Protein Science. 13 (10): 2600–12. doi:10.1110/ps.04850804. PMC 2286553. PMID 15340167.
- Lomize AL, Pogozheva ID, Lomize MA, Mosberg HI (Jun 2006). "Positioning of proteins in membranes: a computational approach" (Free full text). Protein Science. 15 (6): 1318–33. doi:10.1110/ps.062126106. PMC 2242528. PMID 16731967.
- Dao-pin S, Anderson DE, Baase WA, Dahlquist FW, Matthews BW (Dec 1991). "Structural and thermodynamic consequences of burying a charged residue within the hydrophobic core of T4 lysozyme". Biochemistry. 30 (49): 11521–9. doi:10.1021/bi00113a006. PMID 1747370.
- Caputo GA, London E (Mar 2003). "Cumulative effects of amino acid substitutions and hydrophobic mismatch upon the transmembrane stability and conformation of hydrophobic alpha-helices". Biochemistry. 42 (11): 3275–85. doi:10.1021/bi026697d. PMID 12641459.
- Schaefer M, van Vlijmen HW, Karplus M (1998). "Electrostatic contributions to molecular free energies in solution". Linkage Thermodynamics of Macromolecular Interactions. Advances in Protein Chemistry. Vol. 51. pp. 1–57. doi:10.1016/S0065-3233(08)60650-6. ISBN 978-0-12-034251-8. PMID 9615168.
- García-Moreno EB, Fitch CA (2004). "Structural interpretation of pH and salt-dependent processes in proteins with computational methods". Energetics of Biological Macromolecules, Part E. Methods in Enzymology. Vol. 380. pp. 20–51. doi:10.1016/S0076-6879(04)80002-8. ISBN 978-0-12-182784-7. PMID 15051331.