We construct a mean-field variational model to study how the dependence of dielectric coefficient (i.e., relative permittivity) on local ionic concentrations affects the electrostatic interaction in an ionic solution near a charged surface. The electrostatic free-energy functional of ionic concentrations, which is the key object in our model, consists mainly of the electrostatic potential energy and the ionic ideal-gas entropy. The electrostatic potential is determined by Poisson's equation in which the dielectric coefficient depends on the sum of concentrations of individual ionic species. This dependence is assumed to be qualitatively the same as that on the salt concentration for which experimental data are available and analytical forms can be obtained by the data fitting. We derive the first and second variations of the free-energy functional, obtain the generalized Boltzmann distributions, and show that the free-energy functional is in general nonconvex. To validate our mathematical analysis, we numerically minimize our electrostatic free-energy functional for a radially symmetric charged system. Our extensive computations reveal several features that are significantly different from a system modeled with a dielectric coefficient independent of ionic concentration. These include the non-monotonicity of ionic concentrations, the ionic depletion near a charged surface that has been previously predicted by a one-dimensional model, and the enhancement of such depletion due to the increase of surface charges or bulk ionic concentrations.

The solute-solvent interface that separates biological molecules from their surrounding aqueous solvent characterizes the conformation and dynamics of such molecules. In this work, we construct a solvent fluid dielectric boundary model for the solvation of charged molecules and apply it to study the stability of a model cylindrical solute-solvent interface. The motion of the solute-solvent interface is defined to be the same as that of solvent fluid at the interface. The solvent fluid is assumed to be incompressible and is described by the Stokes equation. The solute is modeled simply by the ideal-gas law. All the viscous force, hydrostatic pressure, solute-solvent van der Waals interaction, surface tension, and electrostatic force are balanced at the solute-solvent interface. We model the electrostatics by Poisson's equation in which the solute-solvent interface is treated as a dielectric boundary that separates the low-dielectric solute from the high-dielectric solvent. For a cylindrical geometry, we find multiple cylindrically shaped equilibrium interfaces that describe polymodal (e.g., dry and wet) states of hydration of an underlying molecular system. These steady-state solutions exhibit bifurcation behavior with respect to the charge density. For their linearized systems, we use the projection method to solve the fluid equation and find the dispersion relation. Our asymptotic analysis shows that, for large wavenumbers, the decay rate is proportional to wavenumber with the proportionality half of the ratio of surface tension to solvent viscosity, indicating that the solvent viscosity does affect the stability of a solute-solvent interface. Consequences of our analysis in the context of biomolecular interactions are discussed.

Near a charged surface, counterions of different valences and sizes cluster; and their concentration profiles stratify. At a distance from such a surface larger than the Debye length, the electric field is screened by counterions. Recent studies by a variational mean-field approach that includes ionic size effects and by Monte Carlo simulations both suggest that the counterion stratification is determined by the ionic valence-to-volume ratios. Central in the mean-field approach is a free-energy functional of ionic concentrations in which the ionic size effects are included through the entropic effect of solvent molecules. The corresponding equilibrium conditions define the generalized Boltzmann distributions relating the ionic concentrations to the electrostatic potential. This paper presents a detailed analysis and numerical calculations of such a free-energy functional to understand the dependence of the ionic charge density on the electrostatic potential through the generalized Boltzmann distributions, the role of ionic valence-to-volume ratios in the counterion stratification, and the modification of Debye length due to the effect of ionic sizes.

Protein–ligand binding is a key biological process at the molecular level. The identification and characterization of small-molecule binding sites on therapeutically relevant proteins have tremendous implications for target evaluation and rational drug design. In this work, we used the recently developed level-set variational implicit-solvent model (VISM) with the Coulomb field approximation (CFA) to locate and characterize potential protein–small-molecule binding sites. We applied our method to a data set of 515 protein–ligand complexes and found that 96.9% of the cocrystallized ligands bind to the VISM-CFA-identified pockets and that 71.8% of the identified pockets are occupied by cocrystallized ligands. For 228 tight-binding protein–ligand complexes (i.e, complexes with experimental pKd values larger than 6), 99.1% of the cocrystallized ligands are in the VISM-CFA-identified pockets. In addition, it was found that the ligand binding orientations are consistent with the hydrophilic and hydrophobic descriptions provided by VISM. Quantitative characterization of binding pockets with topological and physicochemical parameters was used to assess the “ligandability” of the pockets. The results illustrate the key interactions between ligands and receptors and can be very informative for rational drug design.

We introduce a software package for the analysis of biomolecular solvation. The package collects computer codes that implement numerical methods for a variational implicit-solvent model (VISM). The input of the package includes the atomic data of biomolecules under consideration and the macroscopic parameters such as solute-solvent surface tension, bulk solvent density and ionic concentrations, and the dielectric coefficients. The output includes estimated solvation free energies and optimal macroscopic solute-solvent interfaces that are obtained by minimizing the VISM solvation free-energy functional among all possible solute-solvent interfaces enclosing the solute atoms. We review the VISM with various descriptions of electrostatics. We also review our numerical methods that consist mainly of the level-set method for relaxing the VISM free-energy functional and a compact coupling interface method for the dielectric Poisson-Boltzmann equation. Such numerical methods and algorithms constitute the central modules of the software package. We detail the structure of the package, format of input and output files, workflow of the codes, and the postprocessing of output data. Our demo application to a host-guest system illustrates how to use the package to perform solvation analysis for biomolecules, including ligand-receptor binding systems. The package is simple and flexible with respect to minimum adjustable parameters and a wide range of applications. Future extensions of the package use can include the efficient identification of ligand binding pockets on protein surfaces.

We incorporate the Poisson-Boltzmann (PB) theory of electrostatics into our variational implicit-solvent model (VISM) for the solvation of charged molecules in an aqueous solvent. In order to numerically relax the VISM free-energy functional by our level-set method, we develop highly accurate methods for solving the dielectric PB equation and for computing the dielectric boundary force. We also apply our VISM-PB theory to analyze the solvent potentials of mean force and the effect of charges on the hydrophobic hydration for some selected molecular systems. These include some single ions, two charged particles, two charged plates, and the host-guest system Cucurbit[7]uril and Bicyclo[2.2.2]octane. Our computational results show that VISM with PB theory can capture well the sensitive response of capillary evaporation to the charge in hydrophobic confinement and the polymodal hydration behavior and can provide accurate estimates of binding affinity of the host-guest system. We finally discuss several issues for further improvement of VISM.

Diffusion is often the rate-determining step in many biological processes. Currently, the two main computational methods for studying diffusion are stochastic methods, such as Brownian dynamics, and continuum methods, such as the finite element method. A previous study introduced a new hybrid diffusion method that couples the strengths of each of these two methods, but was limited by the lack of interactions among the particles; the force on each particle had to be from an external field. This study further develops the method to allow charged particles. The method is derived for a general multidimensional system and is presented using a basic test case for a one-dimensional linear system with one charged species and a radially symmetric system with three charged species.

We design and implement numerical methods for the incompressible Stokes solvent flow and solute-solvent interface motion for nonpolar molecules in aqueous solvent. The balance of viscous force, surface tension, and van der Waals type dispersive force leads to a traction boundary condition on the solute-solvent interface. To allow the change of solute volume, we design special numerical boundary conditions on the boundary of a computational domain through a consistency condition. We use a finite difference ghost fluid scheme to discretize the Stokes equation with such boundary conditions. The method is tested to have a second-order accuracy. We combine this ghost fluid method with the level-set method to simulate the motion of the solute-solvent interface that is governed by the solvent fluid velocity. Numerical examples show that our method can predict accurately the blow up time for a test example of curvature flow and reproduce the polymodal (e.g., dry and wet) states of hydration of some simple model molecular systems.

Mean-field methods, such as the Poisson-Boltzmann equation (PBE), are often used to calculate the electrostatic properties of molecular systems. In the past two decades, an enhancement of the PBE, the size-modified Poisson-Boltzmann equation (SMPBE), has been reported. Here, the PBE and the SMPBE are reevaluated for realistic molecular systems, namely, lipid bilayers, under eight different sets of input parameters. The SMPBE appears to reproduce the molecular dynamics simulation results better than the PBE only under specific parameter sets, but in general, it performs no better than the Stern layer correction of the PBE. These results emphasize the need for careful discussions of the accuracy of mean-field calculations on realistic systems with respect to the choice of parameters and call for reconsideration of the cost-efficiency and the significance of the current SMPBE formulation.

Ligand-receptor binding and unbinding are fundamental biomolecular processes and particularly essential to drug efficacy. Environmental water fluctuations, however, impact the corresponding thermodynamics and kinetics and thereby challenge theoretical descriptions. Here, we devise a holistic, implicit-solvent, multimethod approach to predict the (un)binding kinetics for a generic ligand-pocket model. We use the variational implicit-solvent model (VISM) to calculate the solute-solvent interfacial structures and the corresponding free energies, and combine the VISM with the string method to obtain the minimum energy paths and transition states between the various metastable ("dry" and "wet") hydration states. The resulting dry-wet transition rates are then used in a spatially dependent multistate continuous-time Markov chain Brownian dynamics simulation and the related Fokker-Planck equation calculations of the ligand stochastic motion, providing the mean first-passage times for binding and unbinding. We find the hydration transitions to significantly slow down the binding process, in semiquantitative agreement with existing explicit-water simulations, but significantly accelerate the unbinding process. Moreover, our methods allow the characterization of nonequilibrium hydration states of pocket and ligand during the ligand movement, for which we find substantial memory and hysteresis effects for binding vs. unbinding. Our study thus provides a significant step forward toward efficient, physics-based interpretation and predictions of the complex kinetics in realistic ligand-receptor systems.