<figure><img src="image_1.png" /></figure>

HAL
open science

# Simultaneous Identification of Multiple Binding Sites in Proteins: A Statistical Mechanics Approach

Patrice Koehl, Marc Delarue, Henri Orland

► To cite this version:

Patrice Koehl, Marc Delarue, Henri Orland. Simultaneous Identification of Multiple Binding Sites in Proteins: A Statistical Mechanics Approach. Journal of Physical Chemistry B, 2021, 125 (19), pp.5052-5067. 10.1021/acs.jpcb.1c02658. pasteur-03269759

HAL Id: pasteur-03269759

https://pasteur.hal.science/pasteur-03269759

Submitted on 3 Nov 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Simultaneous identification of multiple
binding sites in proteins:
A statistical mechanics approach

Patrice Koehl,*,† Marc Delarue,*,‡ and Henri Orland*,¶

†Department of Computer Science and Genome Center, University of California, Davis,
CA 95616, USA

‡Architecture et Dynamique des Macromolécules Biologiques, Department of Structural Biology and Chemistry, UMR 3528 du CNRS, Institut Pasteur, 75015 Paris, France

¶Institut de Physique Théorique, Université Paris-Saclay, CEA, 91191 Gif/Yvette Cedex,
France

E-mail: koehl@cs.ucdavis.edu; delarue@pasteur.fr; henri.orland@ipht.fr
We present an extension of the Poisson-Boltzmann model in which the solute of interest is immersed in an assembly of self-orienting Langevin water dipoles, anions, cations, and hydrophobic molecules, all of variable densities. Interactions between charges are controlled by electrostatics, while hydrophobic interactions are modeled with a Yukawa potential. We impose steric constraints by assuming that the system is represented on a cubic lattice. We also assume incompressibility, i.e. all sites of the lattice are occupied. This model, which we refer to as the Hydrophobic Dipolar Poisson Boltzmann Langevin (HDPBL) model, leads to a system of two equations whose solutions give the water dipole, salt, and hydrophobic molecule densities, all of them in the presence of the others in a self-consistent way. We use those to study the organization of the ions, co-solvent and solvent molecules around proteins. In particular, peaks of densities are expected to reveal, simultaneously, the presence of compatible binding sites of different kinds on a protein. We have tested and validated the ability of HDPBL to detect pockets in proteins that bind to hydrophobic ligands, polar ligands and charged small probes as well as to characterize the environment of membrane proteins.

# 1. Introduction

Proteins are unique among biomolecules in that their function is modulated by a wide variety of small molecules that use different types of interactions in their binding sites. These binding sites can bind natural substrates, as observed for example in enzyme active sites, as well as in allosteric regulatory sites. They are also the target of artificial or non natural ligands, developed by the pharmaceutical industry as drugs that can inhibit or activate protein function. The identification and characterization of those binding sites are therefore essential steps for understanding and controlling protein function by structure-based drug design, and as such have been central themes in structural biology. Many experimental methods have
been developed for this purpose, based on NMR,¹ X-ray crystallography,² as well as other techniques.³ Those techniques, however, are time-consuming and often expensive. This has led to parallel developments in computational structural biology with the same goals of identifying druggable binding sites in proteins.⁴⁻⁶ The corresponding methods are usually fast and easy to implement. Their success, however, is often mitigated by the fact that they rely on severe simplifications, such as ignoring solvent effects. A successful computational technique for identifying drug-binding sites requires adequate and efficient sampling of the conformation of the ligand and of the environment within the putative binding site, including the presence of ordered water molecules and salt. Molecular dynamics (MD) simulations provides a rigourous framework to implement such sampling.⁷⁻⁹ Of particular interest, the incorporation of co-solvents to mimic the ligand in such simulations improves the sampling and the identification of binding sites.¹⁰⁻¹⁵ Molecular dynamics, however, is time-consuming as proper sampling requires long simulation times. In this paper we propose a complementary technique based on statistical mechanics for identifying simultaneously one or several binding sites of different kinds in a given protein.

In an MD simulation interactions between atoms are usually described by semi-empirical "force fields", with different levels of approximations (for recent reviews, see¹⁶⁻¹⁹). Applications of these force fields imply that the positions of all atoms be known. While this seems to be a simple requirement, it is unfortunately difficult to meet when modeling solvated large biomolecular systems. This is mostly due to the inherent difficulties in accounting for the mobile solvent molecules and ions that surrounds the solute and of the size of the system that increases significantly when one includes solvent and ions in the simulation. To circumvent these difficulties, there has been continuous efforts to develop simplified models that remain physically accurate. Many of these models consider the solvent implicitly, reducing the solute-solvent interactions to their mean field characteristics. In such models the solvent is treated as a dielectric continuum; they are referred to as continuum dielectric models. The Poisson-Boltzmann theory is one such model, unquestionably the most popular, which pro-
vides a framework for calculating the electrostatics solvation free energy of a solute in such a dielectric continuum. The corresponding PB equation (PBE) is not however the remedy to all problems associated with characterizing the electrostatics interactions for a biomolecule: it remains a mean field approximation, with known limitations. Many improvements have been proposed to the PB equation. The novel technique proposed in this paper is one such extension. We first review the current theoretical developments around the PB model and their relationships with our new model.

PBE is only a mean-field approximation to the multibody problem of solvent-solute electrostatics interactions. It is based on several approximations that have proved to be limitations in some cases. For example, PBE does not include effects due to ion size or ion-ion correlations in its treatment (for review, see Grochowski and Trylska<sup>20</sup>). Solutions have been proposed to account for at least ion size using either a single size<sup>21</sup> or two different sizes,<sup>22</sup> yielding a size-modified Poisson-Boltzmann equation (SMPB), which is used to study ion channels<sup>23</sup> and electrolyte solutions.<sup>24</sup> Classical PBE handles solvent implicitly, with a fixed dielectric constant inside the biomolecule (usually set to 2-4) that abruptly jumps to 80 at the interface between the biomolecule and the solvent. This approximation leads to too much importance to the definition of this interface, usually set to the molecular surface or the vdW surface of the solute.<sup>25</sup> Nicholls and colleagues<sup>26</sup> proposed to solve this problem by using a Gaussian representation of the atoms of the solutes. Their solution, however, does not circumvent an even more important limitation of the PB equation: because of polarization effects in the vicinity of charges, a representation of the solvent as a homogeneous dielectric medium is bound to be erroneous close to the interface. We have developed an extension to the PB equation in which the solvent is described as an assembly of interacting dipoles on a lattice gas to account for the non uniform dielectric property of the solvent.<sup>27-31</sup> We referred to the corresponding equation as the Dipolar Poisson-Boltzmann-Langevin (DPBL) equations. Here we describe an extension to DPBL, which we refer to as HDPBL, in which we allow for the addition of an hydrophobic cosolvent in the environment of the solute of interest.
These cosolvent molecules interact with each other and with hydrophobic “charges” set on hydrophobic groups inside the solute (usually CH2 and CH3 groups). Those interactions are represented with a Yukawa potential. The electrostatic field and an additional attractive field associated with the Yukawa potential are found to be solution of a system of two PB-like partial differential equations (PDE). Those fields are used to compute the water, ions, and co-solvent densities that are then used to map the characteristics of pockets in the neighborhood of the protein of interest.

The HDPBL model proposed in this paper is an alternate approach to the cosolvent-enhanced MD simulations. Its novelty and possible advantages compared to these simulations are based on

i) *Improved sampling.* The conformational space and degrees of freedom of solvent, ions, and co-solvent molecules within and around the target biomolecule are sampled efficiently and their densities are computed self-consistently.

ii) *Multiprobe exploration.* HDPBL enables incorporation of multiple probe types in the analyses. The output of HDPBL are the densities of dipoles representing a polar solvent, anions, cations, and hydrophobic cosolvent molecules as a function of position, thereby allowing for the identification of polar, positively charged, negatively charged, and hydrophobic binding sites simultaneously.

iii) *Efficiency.* MD simulations solve Newton's equations at the atomic level numerically to construct a picture of the dynamics of the biomolecule of interest. Those simulations require significant computing time to provide adequate sampling. In contrast, HDPBL relies on solving once a system of coupled second order PDEs. We will show that this can be performed efficiently using standard techniques for solving elliptic PDEs.

This paper is organized as follows. The next section provides a complete overview of the HDPBL model. The following section describes our implementation of a solver of HDPBL in a program we refer to as AquaVit. AquaVit is a continuation of AquaSol;³¹ it is heavily based
on the package MG developed by Michael Holst.³² The next section provides examples of the usefulness of HDPBL for detecting and characterizing binding sites in proteins. We conclude the paper with a discussion on the limitations of HDPBL, of possible improvements that would circumvent those limitations, as well as on possible developments around HDPBL.

# 2. The HDPBL model

## 2.1 The system: solute, salt, solvent, and cosolvent

We consider a fixed solute molecule (protein) in a solution containing water, a $z : z$ salt (e.g. NaCl for which the valence $z$ is 1), and a cosolvent, i.e. neutral hydrophobic molecules. Those hydrophobic molecules will serve as probes to assess the hydrophobic environment within, and around the solute of interest. They may be considered as methane, xenon, krypton, or any other small hydrophobic molecule.

We assume that the system comprises $N_w$ water molecules, $N_s$ salt molecules (that is $N_s$ cations of charge $+ze_c$ and $N_s$ anions of charge $-ze_c$ where $e_c$ is the elementary electronic charge) and $N_h$ inert hydrophobic molecules. Those molecules are only present in the region outside of the solute, which is characterized by its molecular surface. The water molecules are modeled as dipoles with fixed dipole moment $p_0$. The charged ions and the water dipoles interact through the Coulomb interaction, and all particles are subject to steric repulsion. To model this steric repulsion, we assume that the system is represented on a cubic lattice gas of lattice spacing $a$ and that at each site of the lattice there is either no particle, or one ion or one water molecule or one hydrophobic molecule. Each mobile particle is thus modeled as a hard sphere with radius $a/2$. We denote respectively by $n_+(\mathbf{r})$, $n_-(\mathbf{r})$, $n_w(\mathbf{r})$, $n_h(\mathbf{r})$ the occupation numbers of site $\mathbf{r}$ by a cation, an anion, a water or a hydrophobic molecule, where each of these numbers are 1 or 0 depending whether the corresponding particle is present or absent. Note that $\mathbf{r}$ represents one of the cuboids of the cubic lattice. The steric constraint
imposes that there is at most one particle at any site $\mathbf{r}$,

$$n_+(\mathbf{r}) + n_-(\mathbf{r}) + n_w(\mathbf{r}) + n_h(\mathbf{r}) = 0 \text{ or } 1 \quad (1)$$

If the system is incompressible, the above sum is strictly equal to 1, meaning that there is exactly one particle at each site. The relationships between the occupation numbers of all species and their number of molecules are given by the additional constraints

$$\begin{align*}
\sum_{\mathbf{r}} n_+(\mathbf{r}) &= N_s \\
\sum_{\mathbf{r}} n_-(\mathbf{r}) &= N_s \\
\sum_{\mathbf{r}} n_w(\mathbf{r}) &= N_w \\
\sum_{\mathbf{r}} n_h(\mathbf{r}) &= N_h
\end{align*}
\tag{2}$$

The corresponding particle densities are given by

$$\begin{align}
\rho_{\pm}(\mathbf{r}) &= \frac{n_{\pm}(\mathbf{r})}{a^3} \nonumber \\
\rho_w(\mathbf{r}) &= \frac{n_w(\mathbf{r})}{a^3} \nonumber \\
\rho_h(\mathbf{r}) &= \frac{n_h(\mathbf{r})}{a^3}
\end{align}
\tag{3}$$

We emphasize that all the occupation numbers of mobile particles are zero inside the solute, since we assume that mobile particles cannot penetrate the solute. Let $c_s$, $c_w$, and $c_h$ be the bulk concentrations of salt, water, and hydrophobic probes, respectively. We introduce the
volume fraction $\Phi$ for each type of particle

$$ 
\begin{aligned}
\Phi_s &= 2c_s a^3 \\
\Phi_w &= c_w a^3 \\
\Phi_h &= c_h a^3
\end{aligned}
\qquad (4) $$

In addition, we consider a volume fraction for vacancies:

$$ \Phi_v = 1 - (\Phi_s + \Phi_w + \Phi_h) \qquad (5) $$

$\Phi_v$ should be understood as follows. If the system is incompressible, $\Phi_v = 0$ and the concentrations of salt, water, and hydrophobic probes are necessarily dependent. Otherwise, $\Phi_v$ is positive, and vacancies are possible in the environment of the solute. In this case, $c_s$, $c_w$ and $c_h$ are independent, although they still need to satisfy $\Phi_s+\Phi_w+\Phi_h = 2c_s a^3+c_w a^3+c_h a^3 \le 1$.

We denote by $v(\mathbf{r}) = \frac{1}{4\pi\varepsilon_0 r}$ the Coulomb potential, where $\varepsilon_0$ is the dielectric permittivity of the vacuum and $r = |\mathbf{r}|$.

The hydrophobic interactions between the hydrophobic molecules are captured by an attractive Yukawa potential

$$ w(\mathbf{r}) = -\frac{w_0 e^{-\kappa r}}{4\pi r} \qquad (6) $$

where $\kappa = 1/l_h$ defines the range of the hydrophobic interaction, and $w_0 > 0$ its strength. The negative sign in Equation 6 denotes the attractive nature of the interaction.
## 2.2 An effective free energy for the system

The canonical partition function of the system on the lattice can be written as

$$ Z_c(N_s, N_w, N_h) = \sum_{\{n_{\pm}(\mathbf{r})=0,1\}} \sum_{\{n_w(\mathbf{r})=0,1\}} \sum_{\{n_h(\mathbf{r})=0,1\}} \delta\left(\sum_{\mathbf{r}} n_+(\mathbf{r}) - N_s\right) \delta\left(\sum_{\mathbf{r}} n_-(\mathbf{r}) - N_s\right) \\ \times \delta\left(\sum_{\mathbf{r}} n_w(\mathbf{r}) - N_w\right) \delta\left(\sum_{\mathbf{r}} n_h(\mathbf{r}) - N_h\right) \\ \times \int \prod_{\mathbf{r}} d\mathbf{p}_w(\mathbf{r}) \delta(\mathbf{p}_w^2(\mathbf{r}) - p_0^2) \\ \exp\left(-\frac{\beta}{2} \sum_{\mathbf{r},\mathbf{r}'} \rho_c(\mathbf{r})v(\mathbf{r}-\mathbf{r}')\rho_c(\mathbf{r}') - \frac{\beta}{2} \sum_{\mathbf{r},\mathbf{r}'} \rho_H(\mathbf{r})w(\mathbf{r}-\mathbf{r}')\rho_H(\mathbf{r}')\right) \quad (7) $$

where $\beta = 1/k_B T$ with T the temperature and $k_B$ the Boltzmann constant. The notation $\sum_{\{n_{\pm}(\mathbf{r})=0,1\}}$ denotes a sum over all possible combinations of values of the occupation numbers (0 or 1) at all lattice sites. The total charge density $\rho_c$ in (7) is the sum of the charge densities of the ions, of the point-like water dipoles and of the charges of the solute

$$ \rho_c(\mathbf{r}) = ze_c(\rho_+(\mathbf{r}) - \rho_-(\mathbf{r})) - \mathbf{p}_w(\mathbf{r}) \cdot \nabla \rho_w(\mathbf{r}) + \rho_f(\mathbf{r}) \quad (8) $$

where $\rho_f(\mathbf{r})$ is the charge density of the fixed charges of the solute at point $\mathbf{r}$ and $\mathbf{p}_w(\mathbf{r})$ is the dipole moment of the water molecule at the same point. These dipole moments have a fixed magnitude $p_0$ but can take all possible orientations. This is accounted for by the integral over all dipole orientations in (7).

Finally, $\rho_H(\mathbf{r})$ is the sum of the density $\rho_h(\mathbf{r})$ of mobile hydrophobic particles and of the density $\rho_p(\mathbf{r})$ of the hydrophobic sites of the fixed solute

$$ \rho_H(\mathbf{r}) = \rho_h(\mathbf{r}) + \rho_p(\mathbf{r}) \quad (9) $$

Going to the grand canonical ensemble and introducing the chemical potentials $\mu$ and
fugacities $\lambda$ of the various species

$$ \lambda_s = e^{\beta\mu_s} $$

$$ \lambda_w = e^{\beta\mu_w} $$

$$ \lambda_h = e^{\beta\mu_h} $$

the grand partition function can be written as

$$ \Xi = \sum_{N_s=0}^{+\infty} \sum_{N_w=0}^{+\infty} \sum_{N_h=0}^{+\infty} \frac{e^{2\beta\mu_s N_s + \beta\mu_w N_w + \beta\mu_h N_h}}{(N_s!)^2 N_w! N_h!} Z_c(N_s, N_w, N_h) \quad (10) $$

where the canonical partition function $Z_c$ is given in (7).

Following the formalism introduced in$^{27-30}$ we perform two Hubbard-Stratonovich transforms within Equation 10 and integrate over the dipole moments. After a few standard manipulations, the partition function can be written in an exact manner as an integral over two fields $\varphi(\mathbf{r})$ and $\psi(\mathbf{r})$ corresponding to the electrostatic and the Yukawa interactions, respectively:

$$ \begin{aligned} \Xi = & \int \mathcal{D}\varphi(\mathbf{r})\mathcal{D}\psi(\mathbf{r})e^{-\left(\frac{\beta}{2}\int d\mathbf{r}d\mathbf{r}'\phi(\mathbf{r})v^{-1}(\mathbf{r}-\mathbf{r}')\phi(\mathbf{r}')+\frac{\beta}{2}\int d\mathbf{r}d\mathbf{r}'\psi(\mathbf{r})w^{-1}(\mathbf{r}-\mathbf{r}')\psi(\mathbf{r}')\right)} \\ & \times \exp\left(-i\beta\int d\mathbf{r}\varphi(\mathbf{r})\rho_f(\mathbf{r}) - \beta\int d\mathbf{r}\psi(\mathbf{r})\rho_p(\mathbf{r})\right) \\ & \times \prod_{\mathbf{r}}\left(\lambda_v + 2\lambda_s \cosh(i\beta z e_c \varphi(\mathbf{r})) + \lambda_w \frac{\sinh(ip_0\beta|\nabla\varphi(\mathbf{r})|)}{ip_0\beta|\nabla\varphi(\mathbf{r})|} + \lambda_h e^{-\beta\psi(\mathbf{r})}\right)^{\gamma(\mathbf{r})} \end{aligned} \quad (11) $$

where we have introduce a pseudo fugacity for vacancies $\lambda_v$ such that $\lambda_v = 0$ if the system is incompressible, and $\lambda_v = 1$ otherwise, and $\gamma(\mathbf{r})$ is the indicator function for the points available to the mobile particles, namely $\gamma(\mathbf{r}) = 1$ outside the solute, and $\gamma(\mathbf{r}) = 0$ inside.

The operator $w^{-1}(\mathbf{r})$ is the inverse of the Yukawa interaction (6), and is given by

$$ w^{-1}(\mathbf{r} - \mathbf{r}') = -\frac{1}{w_0}(-\nabla^2 + \kappa^2)\delta(\mathbf{r} - \mathbf{r}') $$
To simplify the notations and the equations, we take the continuous limit $a \to 0$ in the lattice product above, and treat the fields $\varphi(\mathbf{r})$ and $\psi(\mathbf{r})$ as defined in the whole space. The sums are replaced by integrals according to

$$ \sum_{\mathbf{r}} = \frac{1}{a^3} \int d\mathbf{r} \qquad (12) $$

where the integral on the right side is over the whole 3D space, and we obtain

$$ \begin{aligned} \Xi = & \int \mathcal{D}\varphi(\mathbf{r}) \mathcal{D}\psi(\mathbf{r}) e^{\frac{\beta\varepsilon_0}{2} \int d\mathbf{r} (\nabla\varphi(\mathbf{r}))^2 - \frac{\beta}{2w_0} \int d\mathbf{r} ((\nabla\psi(\mathbf{r}))^2 + \kappa^2\psi^2(\mathbf{r}))} \\ & \times \exp \left( -i\beta \int d\mathbf{r} \varphi(\mathbf{r}) \rho_f(\mathbf{r}) - \beta \int d\mathbf{r} \psi(\mathbf{r}) \rho_p(\mathbf{r}) \right) \\ & \times \exp \left( \frac{1}{a^3} \int d\mathbf{r} \gamma(\mathbf{r}) \ln \left( \lambda_v + 2\lambda_s \cosh(i\beta z e_c \varphi(\mathbf{r})) + \lambda_w \frac{\sinh(ip_0\beta |\nabla\varphi(\mathbf{r})|)}{ip_0\beta |\nabla\varphi(\mathbf{r})|} + \lambda_h e^{-\beta\psi(\mathbf{r})} \right) \right) \end{aligned} \qquad (13) $$

The functional integral in Equation 13 is evaluated by the Saddle-Point Approximation. This method, which is also called Mean-Field Theory, consists of minimizing the exponent in the above equation with respect to the two fields $\varphi$ and $\psi$. The field $\varphi$ is pure imaginary at the saddle-point, while $\psi$ is real, and the exponent above can be written as an effective free energy as

$$ \begin{aligned} \mathcal{F} = & -\frac{\varepsilon_0}{2} \int d\mathbf{r} (\nabla\varphi(\mathbf{r}))^2 + \frac{1}{2w_0} \int d\mathbf{r} ((\nabla\psi(\mathbf{r}))^2 + \kappa^2\psi^2(\mathbf{r})) \\ & + \int d\mathbf{r} \varphi(\mathbf{r}) \rho_f(\mathbf{r}) + \int d\mathbf{r} \psi(\mathbf{r}) \rho_p(\mathbf{r}) \\ & - \frac{k_B T}{a^3} \int d\mathbf{r} \gamma(\mathbf{r}) \ln \left( \lambda_v + 2\lambda_s \cosh(\beta z e_c \varphi(\mathbf{r})) + \lambda_w \frac{\sinh(ip_0\beta |\nabla\varphi(\mathbf{r})|)}{ip_0\beta |\nabla\varphi(\mathbf{r})|} + \lambda_h e^{-\beta\psi(\mathbf{r})} \right) \end{aligned} \qquad (14) $$
## 2.3 Optimizing the free energy

The mean field equations are the Euler-Lagrange equations obtained by minimizing 14 with respect to $\varphi$ and $\psi$

$$ 
\begin{aligned}
-\varepsilon_0 \nabla^2 \varphi(\mathbf{r}) &= \rho_f(\mathbf{r}) - \frac{2\lambda_s}{a^3} z e_c \gamma(\mathbf{r}) \frac{\sinh(\beta z e_c \varphi(\mathbf{r}))}{\mathcal{D}(\mathbf{r})} + \frac{p_0}{a^3} \lambda_w \gamma(\mathbf{r}) \nabla \cdot \left( \frac{\nabla \varphi(\mathbf{r})}{|\nabla \varphi(\mathbf{r})|} \frac{g(p_0 \beta |\nabla \varphi(\mathbf{r})|)}{\mathcal{D}(\mathbf{r})} \right) \\
\frac{1}{w_0} (-\nabla^2 + \kappa^2) \psi(\mathbf{r}) &= -\rho_p(\mathbf{r}) - \frac{\lambda_h}{a^3} \gamma(\mathbf{r}) \frac{e^{-\beta \psi(\mathbf{r})}}{\mathcal{D}(\mathbf{r})} \\
\mathcal{D}(\mathbf{r}) &= \lambda_v + 2\lambda_s \cosh(\beta z e_c \varphi(\mathbf{r})) + \lambda_w \frac{\sinh(p_0 \beta |\nabla \varphi(\mathbf{r})|)}{p_0 \beta |\nabla \varphi(\mathbf{r})|} + \lambda_h e^{-\beta \psi(\mathbf{r})}
\end{aligned}
\quad (15) $$

where

$$ g(x) = \frac{\cosh x}{x} - \frac{\sinh x}{x^2} $$

Note that $\varphi(\mathbf{r}) \to 0$ and $\psi(\mathbf{r}) \to \psi_0$ as $\mathbf{r} \to +\infty$, i.e. in the bulk, far from the solute. All coefficients in those equations are computed either from physical constants, or from input information describing the system, with the exception of the fugacities and $\psi_0$, which we derive now.

The fugacities are determined by the equations

$$ 
\begin{aligned}
-\lambda_s \frac{\partial(\beta F)}{\partial \lambda_s} &= N_s \\
-\lambda_w \frac{\partial(\beta F)}{\partial \lambda_w} &= N_w \\
-\lambda_h \frac{\partial(\beta F)}{\partial \lambda_h} &= N_h
\end{aligned}
\quad (16) $$
These equations translate into

$$ 
\begin{aligned}
\int d\mathbf{r} \frac{2\lambda_s}{\mathcal{D}(\mathbf{r})} \cosh (\beta z e_c \varphi(\mathbf{r})) &= a^3 N_s \\
\int d\mathbf{r} \frac{\lambda_w}{\mathcal{D}(\mathbf{r})} \frac{\sinh (p_0 \beta |\nabla \varphi(\mathbf{r})|)}{p_0 \beta |\nabla \varphi(\mathbf{r})|} &= a^3 N_w \\
\int d\mathbf{r} \frac{\lambda_h}{\mathcal{D}(\mathbf{r})} e^{-\beta \psi(\mathbf{r})} &= a^3 N_h
\end{aligned}
\qquad (17) $$

Assuming that the volume of the solution is large compared to that of the solute,

$$ 
\begin{aligned}
\frac{\lambda_s}{\lambda_v + 2\lambda_s + \lambda_w + \lambda_h e^{-\beta\psi_0}} &= a^3 c_s = \frac{\Phi_s}{2} \\
\frac{\lambda_w}{\lambda_v + 2\lambda_s + \lambda_w + \lambda_h e^{-\beta\psi_0}} &= a^3 c_w = \Phi_w \\
\frac{\lambda_h e^{-\beta\psi_0}}{\lambda_v + 2\lambda_s + \lambda_w + \lambda_h e^{-\beta\psi_0}} &= a^3 c_H = \Phi_h
\end{aligned}
\qquad (18) $$

where $\Phi_s$, $\Phi_w$ and $\Phi_h$ the volume fractions defined in Equation 4. Recall that $\Phi_v = 1.0 - \Phi_s - \Phi_w - \Phi_h$ is the molar fraction of vacancies. We consider two cases:

i) *Compressible system*

In a compressible system, we allow for vacancies and $\Phi_v \neq 0$ and $\lambda_v = 1$. Solving the system in Equation 18, we get:

$$ 
\begin{aligned}
2\lambda_s &= \frac{\Phi_s}{\Phi_v} \\
\lambda_w &= \frac{\Phi_w}{\Phi_v} \\
\lambda_h e^{-\beta\psi_0} &= \frac{\Phi_h}{\Phi_v}
\end{aligned}
\qquad (19) $$

ii) *Incompressible system*

In the case of an incompressible system, there are no vacancies, and $\Phi_s + \Phi_w + \Phi_h = 1$. The fugacities are not independent, and it is possible to choose for example $\lambda_w = 1$.
The fugacities are then given by

$$
\begin{align}
\lambda_w &= 1 \nonumber \\
2\lambda_s &= \frac{\Phi_s}{\Phi_w} \nonumber \\
\lambda_h e^{-\beta\psi_0} &= \frac{\Phi_h}{\Phi_w} \nonumber \\
\Phi_w &= 1 - \Phi_s - \Phi_h \tag{20}
\end{align}
$$

The bulk value $\psi_0$ is given by

$$
\psi_0 = - \frac{w_0 \Phi_h}{\kappa^2 a^3} \qquad (21)
$$

Using the values derived above for the fugacities, In all cases, the meanfield equations 15
can be rewritten as

$$
-\epsilon_0 \nabla^2 \varphi(\mathbf{r}) = \rho_f(\mathbf{r}) - \frac{\Phi_s}{a^3} z e_c \gamma(\mathbf{r}) \frac{\sinh(\beta z e_c \varphi(\mathbf{r}))}{\mathcal{D}_a(\mathbf{r})} + \frac{p_0}{a^3} \Phi_w \gamma(\mathbf{r}) \nabla \cdot \left( \frac{\nabla \varphi(\mathbf{r})}{|\nabla \varphi(\mathbf{r})|} \frac{g(p_0 \beta |\nabla \varphi(\mathbf{r})|)}{\mathcal{D}_a(\mathbf{r})} \right) \quad (22a)
$$

$$
\frac{1}{w_0} (-\nabla^2 + \kappa^2) \psi(\mathbf{r}) = -\rho_p(\mathbf{r}) - \frac{\Phi_h}{a^3} \gamma(\mathbf{r}) \frac{e^{-\beta(\psi(\mathbf{r}) - \psi_0)}}{\mathcal{D}_a(\mathbf{r})} \quad (22b)
$$

where

$$
\mathcal{D}_a(\mathbf{r}) = \Phi_v + \Phi_s \cosh(\beta z e_c \varphi(\mathbf{r})) + \Phi_w \frac{\sinh(p_0 \beta |\nabla \varphi(\mathbf{r})|)}{p_0 \beta |\nabla \varphi(\mathbf{r})|} + \Psi_H e^{-\beta \psi(\mathbf{r})} \quad (23)
$$

Note that $\mathcal{D}_a(\mathbf{r}) \to 1$ when $\mathbf{r} \to +\infty$.

2.4 The HDPBL system of equations

The meanfield equations 22 given above fully describe the system under study. Equation 22a
is a Dipolar Poisson-Boltzmann Langevin (DPBL) equation,²⁷,²⁹,³¹ while equation 22b is a
Poisson-Boltzmann-like equation that relates to the hydrophobic interactions involving the
hydrophobic probes in the solvent and the hydrophobic charges on the solute. As a consequence, we refer to this system of equations as the Hydrophobic Dipolar Poisson-Boltzmann Langevin equations, or HDPBL equations in short. In the following, we provide modified equations involving dimensionless potentials that are more amenable to a numerical solution, and we derive all constants necessary for those equations.

### 2.4.1 Revisiting the DPBL-like equation 22a

The electrostatic potential $\varphi(\mathbf{r})$ and the field $\psi(\mathbf{r})$ are expressed in volts and Joules in the SI unit system, respectively. It is common to consider instead the dimensionless potentials $u(\mathbf{r})$ and $v(\mathbf{r})$ defined as:

$$u(\mathbf{r}) = \frac{e_c \varphi(\mathbf{r})}{k_B T} = \beta e_c \varphi(\mathbf{r}) \qquad (24a)$$

$$v(\mathbf{r}) = \frac{\psi(\mathbf{r}) - \psi_0}{k_B T} = \beta(\psi(\mathbf{r}) - \psi_0) \qquad (24b)$$

(24c)

Equation 22a can then be rewritten as

$$-\varepsilon_0 \nabla^2 u(\mathbf{r}) = \beta e_c^2 \rho_{fd}(\mathbf{r}) - \frac{\Phi_s}{a^3} \beta z e_c^2 \gamma(\mathbf{r}) \frac{\sinh(zu(\mathbf{r}))}{\mathcal{D}_1(\mathbf{r})} + \frac{\beta e_c^2 p_e}{a^3} \Phi_w \gamma(\mathbf{r}) \nabla \cdot \left( \frac{\nabla u(\mathbf{r})}{|\nabla u(\mathbf{r})|} \frac{g(p_e |\nabla u(\mathbf{r})|)}{\mathcal{D}_1(\mathbf{r})} \right) \quad (25)$$

where we have defined $p_e = p_0/e_c$, and

$$\mathcal{D}_1(\mathbf{r}) = \Phi_v + \Phi_s \cosh(zu(\mathbf{r})) + \Phi_w \frac{\sinh(p_e |\nabla u(\mathbf{r})|)}{p_e |\nabla u(\mathbf{r})|} + \Phi_h e^{-v(\mathbf{r})} \qquad (26)$$

In equation 25, $\rho_{fd}$ is the density of fixed charges expressed as fraction of electrons, hence the $e_c^2$ as a coefficient, where the first $e_c$ comes from the change to a dimensionless potential, and the second $e_c$ is factored from the solute charges.
Let us introduce the function $F(x) = \frac{g(x)}{x}$. Note that

$$F(x) = \frac{\cosh(x)}{x^2} - \frac{\sinh(x)}{x^3} = \frac{\sinh(x)}{x^2} \mathcal{L}(x) \quad (27)$$

where $\mathcal{L}(x) = \frac{1}{\tanh(x)} - \frac{1}{x}$ is the Langevin function and $F(x) \to 1/3$ when $x \to 0$.

After a few standard manipulations, taking into account that $\nabla\gamma(\mathbf{r}) = \mathbf{0}$, equation 25 can be rewritten as:

$$-\nabla \cdot \left( \left( 1 + \gamma(\mathbf{r}) C_w \Phi_w \frac{F(p_e |\nabla u(\mathbf{r})|)}{\mathcal{D}_1(\mathbf{r})} \right) \nabla u(\mathbf{r}) \right) + \gamma(\mathbf{r}) C_s \Phi_s \frac{\sinh(zu(\mathbf{r}))}{\mathcal{D}_1(\mathbf{r})} = C_f \rho_{fd}(\mathbf{r}) \quad (28)$$

where we have introduced the three constants:

$$C_w = \frac{4\pi l_B p_e^2}{a^3} \quad (29a)$$

$$C_s = \frac{z4\pi l_B}{a^3} \quad (29b)$$

$$C_f = 4\pi l_B \quad (29c)$$

where $l_B$ is the Bjerrum length in vacuum, namely

$$l_B = \frac{\beta e_c^2}{4\pi\varepsilon_0} \quad (30)$$

As written, equation 28 is a PB-like second order differential equation, with a relative permittivity $\varepsilon(\mathbf{r}, u, v)$ defined as:

$$\varepsilon(\mathbf{r}, u, v) = 1 + \gamma(\mathbf{r}) C_w \Phi_w \frac{F(p_e |\nabla u(\mathbf{r})|)}{\mathcal{D}_1(\mathbf{r})} \quad (31)$$

$\varepsilon(\mathbf{r}, u, v)$ is 1 inside the solute, and depends on both the position $\mathbf{r}$ and on the values of the potentials $u(\mathbf{r})$ and $v(\mathbf{r})$ at that position, for $\mathbf{r}$ outside the molecule. Equation 31 gives a self-consistent representation of the dielectric permittivity of the system.
### 2.4.2 Revisiting the PB-like equation 22b

Using the dimensionless potentials $u(\mathbf{r})$ and $v(\mathbf{r})$ defined in equations 24, equation 22b becomes

$$ \frac{1}{w_0} (-\nabla^2 + \kappa^2) (v(\mathbf{r}) + v_0) = -\beta \rho_p(\mathbf{r}) - \frac{\beta \Phi_h}{a^3} \gamma(\mathbf{r}) \frac{e^{-v(\mathbf{r})}}{D_1(\mathbf{r})} \quad (32) $$

where

$$ v_0 = \beta \psi_0 = - \frac{\beta w_0 \Phi_h}{\kappa^2 a^3} \quad (33) $$

Note that $\beta w_0$ is a length, which we write as $l_Y$, and

$$ v_0 = \frac{l_Y \Phi_h}{\kappa^2 a^3} \quad (34) $$

We rewrite Equation 32 as:

$$ -\nabla^2 v(\mathbf{r}) + \kappa^2 (v(\mathbf{r}) + v_0) + \gamma(\mathbf{r}) \frac{l_Y}{a^3} \frac{e^{-v(\mathbf{r})}}{D_1(\mathbf{r})} = -l_Y \rho_p(\mathbf{r}) \quad (35) $$

## 2.5 The densities or water dipoles, ions, and hydrophobic probes

Once the dimensionless fields $u^{MF}(\mathbf{r})$ and $v^{MF}(\mathbf{r})$ have been derived as solutions of the HDPBL system of equations, the densities of the various molecules are given by

i) Anions and cations:

$$ \rho_{\pm}(\mathbf{r}) = \frac{1}{a^3} \frac{\Phi_s e^{\mp z u^{MF}(\mathbf{r})}}{2 D_1^{MF}(\mathbf{r})} \quad (36) $$

ii) Salt:

$$ \rho_s(\mathbf{r}) = \frac{1}{a^3} \frac{\Phi_s \cosh(z u^{MF}(\mathbf{r}))}{D_1^{MF}(\mathbf{r})} \quad (37) $$
iii) Water dipoles

$$ \rho_w(\mathbf{r}) = \frac{1}{a^3} \frac{\Phi_w}{\mathcal{D}_1^{MF}(\mathbf{r})} \frac{\sinh(p_e |\nabla u^{MF}(\mathbf{r})|)}{p_e |\nabla u(\mathbf{r})|} \quad (38) $$

iv) Hydrophobic probes

$$ \rho_h(\mathbf{r}) = \frac{1}{a^3} \frac{\Phi_h}{\mathcal{D}_1^{MF}(\mathbf{r})} e^{-v(\mathbf{r})} \quad (39) $$

where $\mathcal{D}_1^{MF}(\mathbf{r}) = \Phi_v + \Phi_s \cosh(zu^{MF}(\mathbf{r})) + \Phi_w \frac{\sinh(p_e |\nabla u^{MF}(\mathbf{r})|)}{p_e |\nabla u^{MF}(\mathbf{r})|} + \Phi_h e^{-v^{MF}(\mathbf{r})}$.

Equation 38 defines the local densities of water around the solute. In parallel, we can also compute the polarisation density

$$ \mathbf{P}(\mathbf{r}) = \frac{p_0}{a^3} \Phi_w \frac{g(p_0 \beta |\mathbf{E}^{MF}(\mathbf{r})|)}{\mathcal{D}_1^{MF}(\mathbf{r})} \hat{\mathbf{E}}^{MF}(\mathbf{r}) \quad (40) $$

where $\hat{\mathbf{E}}$ is the unit vector of the electric field $\mathbf{E}$. Using the expression for the water density, we have

$$ \mathbf{P}(\mathbf{r}) = \rho_w(\mathbf{r}) \mathcal{L}(p_0 \beta |\mathbf{E}^{MF}(\mathbf{r})|) \hat{\mathbf{E}}^{MF}(\mathbf{r}) \quad (41) $$

where $\mathcal{L}$ is the usual Langevin function (see above). For small electric field $\mathbf{E}$, eq. 41 becomes the standard linear relation

$$ \mathbf{P}(\mathbf{r}) = \frac{p_0^2 \beta c_w}{3} \mathbf{E}^{MF}(\mathbf{r}) = \alpha \mathbf{E}^{MF}(\mathbf{r}) \quad (42) $$

Note that if we use the expression for the polarisation density $\mathbf{P}(\mathbf{r})$ we can rewrite the DPBL-like equation 28 as

$$ \nabla \cdot (\varepsilon_0 \mathbf{E}(\mathbf{r}) + \mathbf{P}(\mathbf{r})) = \rho_f(\mathbf{r}) + \rho_{\text{ions}}(\mathbf{r}) \quad (43) $$
### 3. Numerical solutions to the HDPBL system of equations

Let us first recall the system of equations HDPBL:

$$ -\nabla \cdot (\varepsilon(\mathbf{r}, u, v)\nabla u(\mathbf{r})) + \gamma(\mathbf{r})C_s \Phi_s \frac{\sinh(zu(\mathbf{r}))}{D_1(\mathbf{r}, u, v)} = C_f \rho_{fd}(\mathbf{r}) \quad (44a) $$

$$ -\nabla^2 v(\mathbf{r}) + \kappa^2(v(\mathbf{r}) + v_0) + \gamma(\mathbf{r}) \frac{l_Y}{a^3} \Phi_h \frac{e^{-v(\mathbf{r})}}{D_1(\mathbf{r}, u, v)} = -l_Y \rho_p(\mathbf{r}) \quad (44b) $$

where $D_1(\mathbf{r}, u, v)$ and $\varepsilon(\mathbf{r}, u, v)$ are functions of the position, $\mathbf{r}$, and of the fields $u$ and $v$:

$$ \varepsilon(\mathbf{r}, u, v) = 1 + \gamma(\mathbf{r})C_w \Phi_w \frac{F(p_e |\nabla u(\mathbf{r})|)}{D_1(\mathbf{r}, u, v)} $$

$$ D_1(\mathbf{r}, u, v) = \Phi_v + \Phi_s \cosh(zu(\mathbf{r})) + \Phi_w \frac{\sinh(p_e |\nabla u(\mathbf{r})|)}{p_e |\nabla u(\mathbf{r})|} + \Phi_h e^{-v(\mathbf{r})} \quad (45) $$

The two equations in this system are dependent as they both involve the two fields $u$ and $v$. They are PB-like, but cannot be solved directly using a PB solver. Equation 44 for example is a second order elliptic non linear PDE, like PB; however, its response coefficients (the coefficients in the divergence term) are not constant, as they are non linear functions of the two fields $u$ and $v$. For equation 44, it is the Helmholtz term that is a non linear function of the two fields $u$ and $v$. Instead of considering a new specific solver for each of those equations, we propose to use a standard inexact Newton method developed for the PB equation by Hold and Saied$^{33}$ within a self-consistent algorithm for solving the HDPBL system, as described in algorithm 1.

This algorithm alternatively solves for $u$ (part 1), given $v$, and then solves for $v$ (part 2), given $u$, until convergence, i.e. until those fields do not change anymore. To solve for $u$, step (1b) sets the diffusion coefficient $\varepsilon_n$ independent of the fields $u$ and $v$. Similarly, step (1c) defines a Helmholtz-like term $H_m$ whose value at position $\mathbf{r}$ only depends on the value of the potential at that position. The PDE in step (1d) is then a PB equation that can be solved directly with a PB algorithm without modification. The updates in step (1e) and (2d) follow a typical trick for self-consistent methods that removes oscillations in the
**Algorithm 1** AquaVit: Self Consistent Newton method for solving the DPBL equation

Initialize $u_0(\mathbf{r}) = 0$ and $v_0(\mathbf{r}) = 0, \forall \mathbf{r}$

for $n = 0, \dots$ until convergence do

(1) Solve equation 44a for $u(\mathbf{r})$ with $v(\mathbf{r})$ fixed

(1a) Initialize a field $\psi(\mathbf{r}) = 0, \forall \mathbf{r}$;

Define $F(\mathbf{r}, \psi) = -\nabla \cdot (\varepsilon(\mathbf{r}, \psi, v_n)\nabla\psi(\mathbf{r})) + \gamma(\mathbf{r})C_s\Phi_s \frac{\sinh(z\psi(\mathbf{r}))}{D_1(\mathbf{r}, \psi, v_n)} - C_f\rho_{fd}(\mathbf{r})$

for $m = 0, \dots$ until convergence do

(1b) Set $\varepsilon_m(\mathbf{r}) = \varepsilon(\mathbf{r}, \psi_m, v_n)$

(1c) Set $D_m(\mathbf{r}) = D_1(\mathbf{r}, \psi_m, v_n)$

Define $H_m(\mathbf{r}, \psi) = \gamma(\mathbf{r})C_s\Phi_s \frac{\sinh(z\psi(\mathbf{r}))}{D_m(\mathbf{r})}$

(1d) Solve the PB-like PDE:

$$\nabla \cdot (\varepsilon_n(\mathbf{r})\nabla\psi(\mathbf{r})) + H_n(\mathbf{r}, \psi(\mathbf{r})) = C_f\rho_{fd}(\mathbf{r})$$

for $\psi_{sol}$, using the inexact Newton method of Holst and Saied³³

(1e) Update $\psi$:

$$\psi_{m+1}(\mathbf{r}) = \lambda\psi_{sol}(\mathbf{r}) + (1 - \lambda)\psi_m(\mathbf{r}), \quad \forall \mathbf{r}$$

(1f) Check for convergence: if $\frac{\sum_{\mathbf{r}} \|\mathbf{F}(\mathbf{r}, \psi_{m+1})\|}{\sum_{\mathbf{r}} \|\mathbf{F}(\mathbf{r}, \psi_0)\|} < \text{TOL}$, stop

end for

(1g) Update $u_{n+1}(\mathbf{r}) = \psi_{m+1}(\mathbf{r}), \quad \forall \mathbf{r}$

(2) Solve equation 44b for $v(\mathbf{r})$ with $u(\mathbf{r})$ fixed

(2a) Initialize a field $\phi(\mathbf{r}) = 0, \forall \mathbf{r}$;

(2b) Define $G(\mathbf{r}, \phi) = \kappa^2(\phi(\mathbf{r}) + v_0) + \gamma(\mathbf{r}) \frac{l_Y}{a^3} \Phi_h \frac{e^{-\phi(\mathbf{r})}}{D_1(\mathbf{r}, u_{n+1}, \phi)}$

(2c) Solve the PB-like PDE:

$$-\nabla^2 v(\mathbf{r}) + G(\mathbf{r}, \phi(\mathbf{r})) = -l_Y \rho_p(\mathbf{r})$$

for $\phi_{sol}$, using the inexact Newton method of Holst and Saied³³

(2d) Update $v_{n+1}(\mathbf{r}) = \lambda\phi_{sol}(\mathbf{r}) + (1 - \lambda)v_n(\mathbf{r}), \quad \forall \mathbf{r}$

(3) Check for convergence: if $\sum_{\mathbf{r}} |v_{n+1}(\mathbf{r}) - v_n(\mathbf{r})| < \text{TOL}$, stop

end for
convergence behavior. Note that here is no need to solve the PDE in step (1d) exactly. As its solution $\psi(\mathbf{r})$ is used as a correction for the solution of the DPBL equation (step 1e), it is appropriate to use an approximation: the number of total iterations may then be higher, but this is compensated by the fact that the amount of work per iteration is smaller.

The algorithm described above was implemented in a software package, AquaVit. AquaVit is written in Fortran and is designed specifically to solve the HDPBL system of equations. AquaVit is mostly inspired from AquaSol,³¹ and uses many routines from the Fortran package MG developed by Michael Holst.³² More information on the implementation is available in the paper describing Aquasol in length.³¹

# 4. Methods

## 4.1 System set up

The coordinates of the atoms of the solute(s) as well as their partial charges are read from a single file under the PQR format used by APBS. For large biomolecules, PQR files can be readily generated from the correspondent PDB³⁴ files using the service PDB2PQR.³⁵ For all examples described below, we used the PARSE parameter dataset³⁶ to assign charges. The PQR file may contain several molecules. The PQR file was then modified to add hydrophobic charges to all atoms. While those are not charge per se, they enable interactions between the hydrophobic probe in the solvent and the solute. Using the nomenclature of CHARMM, all atoms identified as CT (aliphatic carbon), CH1E (extended carbon, with 1 hydrogen), CH2E (extended carbon, with 2 hydrogens), CH3E (extended carbon, with 3 hydrogens), S (sulfur) and SH1E (extended atom S with one hydrogen) were assigned a non-zero hydrophobic charge of +2, while all other atoms had no non polar charges.

AquaVit starts by building a regular mesh around the solutes. The mesh is positioned such that its center matches with the center of the solute. The user provides the number of points and the mesh spacing in each directions. AquaVit checks that there is at least a
distance of $2l_B^w$, ($l_B^w$ being the Bjerrum length in water at 300K, i.e. approximately 7 Å) from
any point on the surface of the solute to the closest edge of the mesh; if this condition is not
met, the mesh size is adjusted accordingly.

The interface between the interior and exterior of the solutes is modeled based on their
molecular surface. The molecular surface is the lower envelope obtained by rolling a water
probe of radius $R_{probe}$ on the vdW surface of the molecule. A full description of how this
surface is computed can be found in Ref. ³¹

Classical treatment of electrostatics assigns a point charge to each atom, usually located
at the center of the sphere representing this atom. The mesh considered in AquaVit is
Cartesian; as such, the centers of the atoms of the solute(s) will most likely not coincide
with its vertices. We therefore need to project the atomic charges on the vertices of the
mesh; we have used trilinear interpolation.

## 4.2 Parameterizing the HDPBL system

The HDPBL system of equations include 10 parameters: the number of vertices in the mesh
in each dimension, the lattice size $a$, the temperature $T$, the concentrations of water $c_w$, salt,
$c_s$, and hydrophobic probes, $c_h$, the valence $z$ of the anions and cations from the salt, the
strength of the water dipole, $p_0$, and the parameters of the Yukawa potential, $l_Y(=\beta w_0)$ and
$l_h(=1/\kappa)$.

The mesh was set with 193 vertices in each direction, $x, y$, and $z$. Those vertices are
equally spaced, and the distance between two vertex is computed automatically based on
the size of the solute and the fact that the borders of the mesh are set to be at least $2l_B^w$
away from the solute.

Assuming incompressibility, in the presence of pure water, we expect $\Phi_w = 1$, i.e. that
$c_w a^3 = 1$, where $c_w$ is the concentration of bulk water, namely 55 M. This leads to $a = 3.11$.
In all the simulations described in the result section, we have considered monovalent (i.e.
$z = 1$) salt at 0.2 M and hydrophobic probes at 1 M. This amount of hydrophobic probes
is equivalent to the amount used in ligand mapping molecular dynamics programs such as SILCS³⁷ and SWISH.³⁸ Again, assuming incompressibility, the concentration of water is then fixed as we have the relation (see above):

$$2c_s a^3 + c_w a^3 + c_h a^3 = 1$$

Using the prescribed concentrations of salt and hydrophobic probes, and the lattice size $a = 3.11$, this gives us an apparent concentration of water $c_w = 53.6$ M.

The parameter $l_h$ defines the range of the Yukawa potential. We have set it equal to the size of the lattice, i.e. $l_h = 3.1$ Å. $l_Y$ is a characteristic length for the Yukawa potential that directly relate to its strength. We have set it to $l_Y = 4$ Å.

The temperature is set to 300K.

The experimental dipole moment of water is 1.85D in the gas phase. We have observed previously that with this value for $p_0$, we could not obtain the correct dielectric permittivity of water in the DPBL model.²⁹ As HDPBL is based on DPBL, we follow the recommendation of increasing $p_0$ which we set at 2.8D.

In Table 1, we give the values of the three constants $C_w$, $C_s$, $C_f$, and $p_e$ for a typical run of AquaVit with $T = 300$K and the values of the different parameters given above.

<table><caption>**Table 1: Typical values for the constants in equation 28**</caption><thead><tr><th>Name</th><th>Expression</th><th>value<sup>(a)</sup></th><th>unit</th></tr></thead><tbody><tr><td>$l_B$</td><td>$\frac{\beta e_c^2}{4\pi\varepsilon_0}$</td><td>556.99995</td><td>Å</td></tr><tr><td>$C_w$</td><td>$\frac{4\pi l_B p_e^2}{a^3}$</td><td>78.51664</td><td>Dimensionless</td></tr><tr><td>$C_s$</td><td>$\frac{z4\pi l_B}{a^3}$</td><td>231.83464</td><td>Å<sup>-2</sup></td></tr><tr><td>$C_f$</td><td>$4\pi l_B$</td><td>6999.46779</td><td>Å</td></tr><tr><td>$p_e$</td><td>$\frac{p_0}{e_c}$</td><td>0.58195</td><td>Å</td></tr></tbody></table>

<footnote>(a) $T = 300$K, $a = 3.11$ Å (size of the lattice), $p_D = 2.8$D (dipole moment of water), and $z = 1$ (valence of salt ions).</footnote>
## 4.2 Output format

The output of AquaVit are the maps corresponding to the densities of the different species in the solvent, namely the water dipoles, the anions, cations, and the hydrophobic probes. While those maps are initially expressed as relative densities with respect to the bulk concentrations of the respective species, we have re-expressed those maps as Z-maps, using:

$$Z(\mathbf{r}) = \frac{\rho(\mathbf{r}) - \mu_{\rho}}{\sigma_{\rho}}$$

where $\mu_{\rho}$ and $\sigma_{\rho}$ are the mean and standard deviation of the densities computed over all positions $\mathbf{r}$ outside of the solute, respectively.

# 5. Results and Discussion

HDPBL differs from a standard PB model, as we have included hydrophobic probes around the solute of interest. Our primary goal is to use those probes to detect hydrophobic pockets in this solute. We have therefore tested the ability of HDPBL to detect pockets that bind to hydrophobic ligands as well as to characterize the environment of membrane proteins.

HDPBL also provides the densities of anions, cations, and water dipole around the solute. We have included a test on detecting binding sites characterized with electrostatics interactions to illustrate their usefulness. Finally, we will describe limitations of HDPBL when it comes to detecting cryptic binding sites.

## 5.1 Detecting hydrophobic pockets in proteins

The set of proteins used for validation includes four proteins from different families and functions: a dihydrofolate reductase (DHFR) from Staphylococcus Aureus (PDB code 2W9H$^{39}$), a maize lipid-transfer protein (PDB code 1FK2$^{40}$), a human retinol binding protein 1 (PDB code 5HBS$^{41}$), and an insect Takeout 1 protein (PDB code 3E8W$^{42}$). The PDB structures
<figure><img src="image_3.png" /><figcaption>Lipid transport protein (PDB code 1F2K)</figcaption></figure>

Figure 1: Hydrophobic occupancy maps, black mesh, superimposed onto the PDB structures for DHFR (PDB code 2W9H) and a lipid binding protein (LBP, PDB code 1FK2). The maps are derived from the densities of hydrophobic probes computed by AquaVit, and represented at +20 σ. Those maps are derived from the apo structure of the protein, i.e. in the absence of all crystallographic ligands and water molecules. In panels (a) and (c) we show the hydrophobic maps for DHFR and LBP superposed to the apo PDB structure, while in panels (b) and (d) we visualize the hydrophobic ligand (trimethoprim for DHFR and myristic acid for LBD) in magenta. All images were generated using Pymol.⁴³

correspond to those proteins bound to a hydrophobic ligand. All ligands, ions, and crystallographic water molecules, however, were removed prior to running AquaVit. In Figures 1 and 2, we show the resulting hydrophobic probe occupancy maps superimposed on the PDB structures for all four proteins, with and without the ligand.

DHFR is the enzyme responsible for the NADPH-dependent reduction of dihydrofolate
to tetrahydrofolate, an essential cofactor in the synthesis of purines, methionine, and other key metabolites. Because of its importance in a wide range of cellular functions, DHFR has been the subject of much research targeting the enzyme for anticancer, antibacterial, and antimicrobial agents. Clinically used compounds targeting DHFR include methotrexate for the treatment of cancer and trimethoprim (TMP) for the treatment of bacterial infections. The active site of DHFR is comprised of a large hydrophobic pocket which serves as the folate-binding site. This pocket was successfully detected by AquaVit, as illustrated in Figure 1(a). This high density region of hydrophobic probes overlaps well with the TMP ligand that was co-crystallized with DHFR and bound within the hydrophobic active site³⁹ (see Figure 1(b)).

Lipid binding proteins (LBP) facilitate the transfer of lipids between membrane. We consider a non specific LBP from maize, whose hydrophobic cavity can accommodate various lipids from C10 to C18. AquaVit was successful in identifying this cavity, as illustrated in Figure 1(c). Interestingly, as we superimpose the ligand found in the PDB structure 1FK2, namely myristic acid, we find that the high density region of hydrophobic probes identified by AquaVit extends beyond this ligand (see Figure 1(d)). Myristic acid is C14, while the pocket found by AquaVit indicates that the hydrophobic cavity of the maize LBP can accommodate larger ligand.⁴⁰

The cellular retinol-binding protein 1(CRBP1) is another example of a protein with a large hydrophobic active site. CRBP1 is important in regulating the uptake, storage and metabolism of retinoids (vitamin A and its derivatives). Its vitamin A binding site is lined with hydrophobic residues that are well conserved among retinol binding proteins. Those residues provide the non-polar interactions that stabilize the retinoid ligand. This pocket was successfully detected by AquaVit, as illustrated in Figure 2(a). The high density region of hydrophobic probes overlaps well with the all trans retinol ligand that was co-crystallized with CRBP1,⁴¹ as observed in Figure 2(b).

Takeout (To) proteins are found exclusively in insects in which they have been proposed to
<figure><img src="image_4.png" /><figcaption>Takeout 1 (PDB code 3E8W)</figcaption></figure>

Figure 2: Hydrophobic occupancy maps, black mesh, superimposed onto the PDB structures for a retinol binding protein (RBP, PDB code 5HBS) and a TakeOut 1 protein (To1, PDB code 3E8W). The maps are derived from the densities of hydrophobic probes computed by AquaVit, and represented at +20 σ. Those maps are derived from the apo structure of the protein, i.e. in the absence of all crystallographic ligands and water molecules. In panels (a) and (c) we show the hydrophobic maps for RBP and To1 superposed to the apo PDB structure, while in panels (b) and (d) we visualize the hydrophobic ligand (trans retinol for RBP and 8-ubiquinone for To1) in magenta. All images were generated using Pymol.⁴³

play important roles in their physiology and behavior.⁴⁴ Of particular interest to us, they have been suggested to bind to hydrophobic ligands. We considered the To1 protein from *Epiphyas postvittana*, a light brown apple moth, whose structure was solved by crystallography at 1.3 Å⁴² in the presence of ubiquinone-8. The crystal structure revealed a 45 Ålong hydrophobic
internal tunnel that extends to the full length of the protein. This pocket was successfully detected by AquaVit, as illustrated in Figure 2(c). Note that the same pocket was originally characterized based on geometry only.⁴⁴ AquaVit provides the additional information that this pocket is amenable to interaction with an hydrophobic ligand. This is confirmed as the high density region of hydrophobic probes overlaps well with the hydrophobic ubiquinone ligand, see Figure 2(d).

## 5.2 Environments of membrane proteins

The subsection above illustrates that AquaVit is able to locate hydrophobic pockets in proteins. However, the solution of the HDPBL system of equations is more comprehensive and also provides information on dipolar density in the presence of a given salt concentration. Here we assess its ability to characterize the polar and non polar environments of membrane proteins. We consider two types of such protein, a member of the G protein coupled receptor (GPCR) family, with a transmembrane domain that consists of a helical bundle (PDB code 2RH1⁴⁵), and a porin that consist of a β-pleated sheet (PDB code 2POR⁴⁶). All ligands, ions, and crystallographic water molecules were removed prior to running AquaVit on those structures. In Figure 3 we show the resulting hydrophobic probe occupancy maps and water dipole occupancy maps, superimposed on the PDB structures, including the crystallographic water molecules.

β₂-adrenergic receptors (β₂AR) are members of the GPCR family that reside predominantly in smooth muscle. Their antagonists are used in particular in the treatment of asthma.⁴⁷ To study the structure of this membrane protein, Cherezov et al designed studied a chimera consisting of β₂AR and the T4 Lysozyme (T4L).⁴⁵ The crystal structure of this chimera reveals a fold for β₂AR composed of a transmembrane domain with a 7 helix bundle, and a standard all-helix fold for the T4L. Interactions between β₂AR and T4L are minimal. AquaVit provides a consistent image of the environment of this chimera, with a mostly hydrophobic environment for the transmembrane domain, and a mostly hydrophilic
<figure><img src="image_5.png" /><figcaption>(b) Porin (PDB code 2POR)</figcaption></figure>

<figure><img src="image_6.png" /></figure>

Figure 3: Hydrophobic occupancy maps, black mesh, and water dipole occupancy maps, blue mesh superimposed onto the PDB structures for (a) the β₂ adrenergic GPCR (PDB code 5HBS) and (b) a porin (PDB code 2POR). The crystallographic waters are shown as cyan spheres. The maps are derived from the densities of hydrophobic probes and water dipoles computed by AquaVit, respectively, and represented at +10 σ for the hydrophobic probes and at +0.3 σ for the water dipoles. Those maps are derived from the apo structure of the proteins, i.e. in the absence of all crystallographic ligands and water molecules.

environment for T4L (see Figure 3 (a)). Interestingly, the water dipole occupancy map superimposes well with the crystallographic water molecules detected in the structure.

Porins are integral membrane proteins that are found in the outer membrane of Gram-negative bacteria, mitochondria and chloroplasts. The crystal structure of the porin from *Rhodobacter capsulatus* reveals a 16-stranded β-barrel, with all β-strands antiparallel and connected to their neighbors.⁴⁶ AquaVit reveals a hydrophobic environment on the outside of the β-pleated sheet, with the loops between the strands in a more hydrophilic environment (Figure 3 (b)). This is consistent with the positions of the crystallographic water molecules that map with the water dipole density map.
### 5.3 Multiprobe analysis of a complex active site

All examples presented above focused on the identification of hydrophobic pockets and characterization of the hydrophobic environment of proteins, based on the analyses of the hydrophobic occupancy maps generated by AquaVit. However, much akin to MixMD⁴⁸ and mLMMMD,⁴⁹ AquaVit can be seen as a multiprobe analysis of the environment of a protein. It generates the densities of water dipole, anion, cation, and hydrophobic molecules surrounding the protein that serve as probes to identify and characterize pockets in the protein. We use the ferredoxin oxidoreductase system to illustrate this functionality of AquaVit.

Ferredoxin-dependent bilin reductases (FBDR) are enzymes that are involved in the reduction of biliverdin (BV) to form phycobilins used for light-perception or light harvesting in plants and cyanobacteria.⁵¹ Several members of this family have been identified in multiple species.⁵² Among those, PebA reduces BV at its C15-C16 double bond to produce 15-16 dihydrobiliverdin (DHBV). The structure of PebA from the cyanobacterium *Synechococcus* sp. WH8020 with its substrate BV was determined at 1.55 Å (PDB code 1X9O⁵⁰). This structure consists of a seven-stranded antiparallel β-sheet surrounded by six α helices (Figure 4(a)). We used AquaVit to solve the HDPBL system of equations that characterize the environment of PeBA. The calculation was performed on the protein structure alone, i.e. in the absence of the BV ligand and crystallographic water molecules. Figure 4(b) illustrates the densities of water dipoles, anions, cations, and hydrophobic probes around the structure of PebA in the form of occupancy maps. The whole protein is surrounded with water dipoles. This is in agreement with the fact that many crystallographic water molecules have been identified. Those molecules superimpose well with the water dipole occupancy map Figure 4(c). Of significant interest are the anion, cation, hydrophobic probe occupancy maps in the active site of PeBA; those maps are illustrated in Figure 4(d). The superposition of the ligand structure and of the two residues (Asp84 and Asn67) that define the central polar centering “pin” of the binding pocket⁵⁰ on those maps shows that the anion densities (in red) map well with the two carboxyl groups on BV (red arrows). We also observe an anion pocket
<figure><img src="image_7.png" /><figcaption>(b)</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>(c)</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Figure 4: (a) Overall structure of PcbA, a ferredoxin-dependent bilin reductase (cartoon mode) in presence of its substrate, biliverdin (BV): PDB code 1X9O.⁵⁰ (b) Superimposition of the occupancy maps of water (+0.3 σ, blue), hydrophobic probes (+10 σ, grey), anion (+20 σ, red) and cation (+20 σ, cyan) on the structure of PcbA, in the absence of BV. (c) The PcbA structure, the water occupancy map, and the crystallographic water molecules (shown as blue spheres). (d) Superposition of the occupancy maps of hydrophobic probes (grey), anion (red), and cation (cyan) with the structure of BV (sticks). We also show the position of residues Asp84 and Asn67 that are known to play an important role in placing the ligand in its binding pocket.⁵⁰ Note the presence of anions density overlapping with the carboxyl groups of BV (red arrows) and forming a pocket in front of the N of the terminal group of Asn67. Note also the presence of cations around the carboxyl group of Asp84 (blue arrow).</figcaption></figure>
in front of the Nitrogen of the terminal group of Asn67 (black arrow). Note that the densities
are expected to match with the chemical structure of the ligand, and be complementary to
the chemical properties of the solute. The hydrophobic probe density map superimposes
well with the hydrophobic parts of BV. Finally, we observe a cation occupancy in the region
surrounding the 4 nitrogen of the rings of BV, as well as a pocket in from of the carboxyl
moiety of Asp84 (blue arrow). These observations highlight the advantage of AquaVit to
account for the multiple species that form the environment of a protein.

**5.4 Computing time**

In Figure 5 we report the computing times for solving the HDPBL system of equations
using AquaVit for three systems, DHFR (PDB code 2W9H), GPCR (PDB code 2RH1), and
PEBA (PDB code 2X9O) under the standard conditions defined in the Methods section. In
addition to the inputs specific to the system under consideration (structure of the solute,
concentrations of water dipoles, salt, and hydrophobic probes, dipole moment of the water
dipole, and parameters of the Yukawa potential for hydrophobic interactions), AquaVit relies
on the parameters of algorithm 1 to solve the HDPBL system. These parameters include
on the size of the Cartesian grid used to solve the PDEs, the tolerance TOL that serves
as stopping criteria when solving the system self consistently, and the parameter λ used
to update the fields u and v (see Algorithm 1 presented in section 3). In all calculations
presented above, TOL is set to 10⁻⁴ and λ is set to 0.7 (this number is somewhat arbitrary
and could in fact be optimized for each system considered). We have also used consistently
a grid of size 193³, but in figure 5 we report the computing time of AquaVit as a function of
this grid size. We observe a near linear dependence of the computing times of AquaVit with
respect to the total number of points in the grid, as expected.³¹ These computing time do not
differ significantly between the three proteins, despite the fact that they are very different
in size (157 residues for 2W9H, 237 residues for 2X9O, and 465 residues for 2RH1). The
average computing time for a grid of 193³ points is 1240 s, i.e. approximately 21 minutes.
<figure><img src="image_11.png" /><figcaption>Figure 5: The computing time required to solve the system of equations HDPBL using AquaVit with a residual error lower than $TOL = 10^{-4}$ (see text for details) is plotted versus the total number of points in the mesh for the three proteins DHFR (PDB code 2W9H), GPCR (PDB code 2RH1), and PEBA (PDB code 2X9O). All computations were performed on an Intel Core i7 processor with 8 cores running at 4.00GHz, and 64GB of memory, although AquaVit was compiled without parallelization options.</figcaption></figure>

We note that the computing times reported are both CPU and clock time, i.e. AquaVit ran on a single core and does not benefit from parallelization. We acknowledge that there is much room for improvement in the implementation of algorithm 1 in AquaVit.

## 5.5 Limitations of the HDBPL model: absence of flexibility

We have shown above that the HDPBL model enables the discovery and characterization of binding sites on proteins. In addition, HDPBL is fast (see Figure 5) or at least competitive with respect to computing time compared to the ligand mapping molecular dynamics simulations. It has one major limitation, however, as it considers the structure of the protein to be static. While it was not an issue in the examples described in the previous sections, we illustrate here a case in which dynamics matter.

Polo-like kinases (PLKs) are a family of serine/threonine kinases related to the polo gene product of *Drosophila melanogaster*. Most PLKs have multiple functions which map
<figure><img src="image_12.png" /><figcaption>Figure 6: (a) Superimposition of the occupancy maps of hydrophobic probes (+20 σ, blue) on the structure of the polo binding domain (PBD) of human polo like kinase 1 (PLK1) (PDB code 1Q4K). The occupancy map is computed using AquaVit, on the PBD structure in the absence of the ligand and crystallographic water molecules. (b) Same as in (a), but with the ligand (the phosphopeptide MetGlnSerpThrProLeu) shown in magenta. Note that the ligand overlaps well with one of the hydrophobic pockets identified by AquaVit. (c) Superimposition of the occupancy maps of hydrophobic probes (+20 σ, blue) on the structure of the polo binding domain (PBD) of human polo like kinase 1 (PLK1) (PDB code 3P37K). The occupancy map is computed using AquaVit, on the PBD structure in the absence of the ligand and crystallographic water molecules. Note the longer hydrophobic pocket compared to 1Q4K, as a result of the change of conformation of Tyr417 and Tyr481. (d) Same as in (a), but with the ligand (the phosphopeptide PheAspProProLeuHisSerpThrAla) shown in magenta. The ligand overlaps well with the long hydrophobic pocket identified by AquaVit.</figcaption></figure>
with their organization in domains. Their C-terminal regions, for example, contain the polo
box domain (PBD), which helps in their sub-cellular localization by binding to serine- or
threonine- phosphorylated sequences on target proteins. Cheng et al have determined the
structure of the PBD domain of the human PLK1 in the presence of a phosphopeptide with
sequence MetGlnSerpThrProLeu, where pThr indicates that the threonin is phosphorylated
(PDB code 1Q4K⁵³). The phosphopeptide was found to bind on the surface of the protein,
in an hydrophobic pocket. We ran AquaVit on a single copy of the PBD domain in the
absence of the peptide ligand and of all crystallographic water molecules. We can identify
the ligand binding hydrophobic pocket of PBD from the resulting hydrophobic probe density
map, as illustrated in Figure 6(a) and (b).

Interestingly, the binding of MetGlnSerpThrProLeu onto PBD was found to induce very
little conformational changes.⁵³ Two subsequent independent studies,⁵⁴,⁵⁵ however, identi-
fied a cryptic hydrophobic binding site that is close to the original phosphate binding site
identified by Cheng et al.⁵³ The opening of this cryptic binding site is the result of a change
in the conformation of the side chains of Tyr417 and Tyr481. This opening enables PBD
to bind longer phosphopeptides. Sledz et al⁵⁵ for example presented the structure of
the complex of PBD with the phosphopeptide PheAspProProLeuHisSerpThrAla (PDB code
3P37). Interestingly, when we ran AquaVit on a single copy of the PBD domain in the
absence of the peptide ligand and of all crystallographic waters in its configuration from
3P37, we were able to identify the longer ligand binding hydrophobic pocket of PBD which
includes the cryptic pockets, which can bind this longer phosphopeptide. This is illustrated
in Figure 6(c) and (d). AquaVit, however, was not able to detect this longer pocket from
the structure in 1Q4K, as residues Tyr417 and Tyr481 were closing the cryptic pocket.
6. Conclusion

We have developed the HDBPL model and implemented it in the program AquaVit as a tool
for identifying and characterizing binding sites in proteins. A protein of interest is immersed
in a lattice gas containing water dipoles, anions, cations, and hydrophobic molecules. The
charged molecules interact between themselves and with the solute charges through electro-
static interactions, while the hydrophobic molecules (including the non polar groups on the
solute) interact based on a Yukawa potential. We impose steric constraints on the lattice,
as well as incompressibility, i.e. all sites of the lattice are occupied. The system (protein of
interest and its environment) is then characterized with an effective free energy that depends
on two fields φ(r) and ψ(r) corresponding to the electrostatics and hydrophobic interactions,
respectively. The Euler-Lagrange equations obtained by minimizing the free energy with
respect to those two fields form a system of two Poisson-Boltzmann like PDEs, HDPBL,
which we solve using a self consistent approach, implemented in the program AquaVit. The
outputs of HDPBL are the densities of the different species, and peaks of densities are ex-
pected to reveal the presence of compatible binding sites. We have tested and validated the
ability of HDPBL to detect pockets that bind to hydrophobic ligands (the DHFR, a lipid
binding protein, HIV protease, and a retinol binding protein), polar ligands (biliverdin), as
well as to characterize the environment of membrane proteins such as GPCR and a porin.

The HDPBL equations form a system of two coupled second order elliptic nonlinear PDE.
While those equations are akin to the PB equation, they cannot be solved directly with a
PB solver, mostly because the two equations are strongly dependent. We have proposed,
however, an algorithm that makes use of a standard PB solver by solving those equations
self consistently. The same approach was used previously to solve the DBPL equation,³¹ as
well as the YULP system of equation.³⁰ This algorithm is relatively simple to implement and
can be adapted to any PB solver, including those based on finite elements methods, which
we did not consider here. We are currently developing a more versatile solver based on those
methods.
In the current implementation of Aquavit, the lattice is populated with point-like electric charges, dipoles and hydrophobic molecules. One could also include finite-size dipoles made of an electric and an hydrophobic charge, as well as an electric dipole attached to a hydrophobic charge. In this case, these entities will react to an hydrophobic field $∇Ψ(\mathbf{r})$, using the same formalism. This could be well adapted to study the solvation of proteins in the presence of large and more complex co-solvents such as acetonitrile or DMSO.

AquaVit, our program for solving the HDPBL system of equations, is fast and as such compares favorably with the ligand mapping molecular dynamics simulations that have been designed with the same goal of detecting and characterizing binding sites in proteins. The latter, however, have the significant advantage that they account for the dynamics of the solute. As such, they have been shown to detect cryptic bind sites in proteins,<sup>11,14,15,49</sup> namely sites that are not accessible unless a structural change occur in the protein. In its current formulation, the HDPBL system of equations assumes that the protein is static; such conformational changes are then inaccessible. As an extension to HDBPL, one could model the dynamics of the solute protein with its low frequency normal modes, using for example a simple elastic model to compute those modes.<sup>56-58</sup> Ultimately, we want to develop AquaVit as a tool for structure-based as well as dynamics-based drug design.

## Acknowledgement

The work discussed here originated from a visit by P.K. at the Institut de Physique Théorique, CEA Saclay, France, during the fall of 2019. He thanks them for their hospitality and financial support. PK acknowledges support from the University of California Multicampus Research Programs and Initiatives (grant No. M21PR3267).
References

(1) Shuker, S.; Hadjuk, P.; Meadows, R.; Fesik, S. Discovering high-affinity ligands for proteins: SAR by NMR. *Science* **1996**, *274*, 1531–1534.

(2) Mattos, C.; Ringe, D. Locating and characterizing binding sites on proteins. *Nat. Biotechnol.* **1996**, *14*, 595–599.

(3) Hardy, J.; Wells, J. Searching for new allosteric sites in enzymes. *Curr. Opin. Struct. Biol.* **2004**, *14*, 706–715.

(4) Henrich, S.; Salo-Ahen, O.; Huang, B.; Rippmann, F.; Cruciani, G.; Wade, R. Computational approaches to identifying and characterizing protein binding sites for ligand design. *J. Molec. Recognit.* **2010**, *23*, 209–219.

(5) Konc, J.; Janežič, D. Binding site comparison for function prediction and pharmaceutical discovery. *Curr. Opin. Struct. Biol.* **2014**, *25*, 34–39.

(6) Zhao, J.; Cao, Y.; Zhang, L. Exploring the computational methods for protein-ligand binding site prediction. *Comput. Struct. Biotechnol. J.* **2020**, *18*, 417–426.

(7) Ivetac, A.; McCammon, J. A molecular dynamics ensemble-based approach for the mapping of druggable binding sites. *Methods Mol. Biol.* **2018**, *819*, 3–12.

(8) Feng, T.; Barakat, K. Molecular dynamics simulation and prediction of druggable binding sites. *Methods Mol. Biol.* **2018**, *1762*, 87–103.

(9) Śledź, P.; Caflisch, A. Protein structure-based drug design: from docking to molecular dynamics. *Curr. Opin. Struct. Biol.* **2018**, *48*, 93–102.

(10) Basse, N.; Kaar, J.; Settanni, G.; Joerger, A.; Rutherford, T.; Fersht, A. Toward the rational design of p53-stabilizing drugs: Probing the surface of the oncogenic Y220C mutant. *Chem. Biol.* **2010**, *17*, 46–56.
(11) Tan, Y.; Sledz, P.; Lang, S.; Stubbs, C.; Spring, D.; Abell, C.; Best, R. Using ligand-mapping simulations to design a ligand selectively targeting a cryptic surface pocket of Polo-Like kinase 1. *Angew. Chem., Int. Ed.* **2012**, *51*, 10078-10081.

(12) Tan, Y.; Spring, D.; Verma, C. The use of chlorobenzene as a probe molecule in molecular dynamics simulations. *J. Chem. Inf. Model.* **2014**, *54*, 1821-1827.

(13) Kalenkiewicz, A.; Grant, B.; Yang, C.-Y. Enrichment of druggable conformations from apo protein structures using cosolvent-accelerated molecular dynamics. *Biology* **2015**, *4*, 344-366.

(14) Kimura, S.; Hu, H.; Ruvinsky, A.; Sherman, W.; Favia, A. Deciphering cryptic binding sites on proteins by mixed-solvent molecular dynamics. *J. Chem. Inf. Model.* **2017**, *57*, 1388-1401.

(15) Schmidt, D.; Boehm, M.; McClendon, C.; Torella, R.; Gohlke, H. Cosolvent-enhanced sampling and unbiased identification of cryptic pockets suitable for structure-based drug design. *J. Chem. Theory Comput.* **2019**, *15*, 3331-3343.

(16) Nerenberg, P.; Head-Gordon, T. New developments in force fields for biomolecular simulations. *Curr. Opin. Struct. Biol.* **2018**, *49*, 129-138.

(17) Huang, J.; Jr., A. D. M. Force field development and simulations of intrinsically disordered proteins. *Curr. Opin. Struct. Biol.* **2018**, *48*, 40-48.

(18) Inakollu, V.; Geerke, D.; Rowley, C.; H, H. Y. Polarisable force fields: what do they add in biomolecular simulations? *Curr. Opin. Struct. Biol.* **2020**, *61*, 182-190.

(19) van der Spoel, D. Systematic design of biomolecular force fields. *Curr. Opin. Struct. Biol.* **2020**, *67*, 18-24.
(20) Grochowski, P.; Trylska, J. Continuum molecular electrostatics, salt effects, and counterion binding - A review of the Poisson-Boltzmann theory and its modifications. *Biopolymers* **2007**, *89*, 93-113.

(21) Borukhov, I.; Andelman, D.; Orland, H. Steric effects in electrolytes: A modified Poisson-Boltzmann equation. *Phys. Rev. Lett.* **1997**, *79*, 435-438.

(22) Chu, V.; Bai, Y.; Lipfert, J.; Herschlag, D.; Doniach, S. Evaluation of ion binding to DNA duplexes using a size-modified Poisson-Boltzmann theory. *Biophys. J.* **2007**, *93*, 3202-3209.

(23) Xie, D.; Audi, S.; Dash, R. A size modified Poisson-Boltzmann ion channel model in a solvent of multiple ionic species: application to voltage-dependent anion channel. *J. Comp. Chem.* **2020**, *41*, 218-230.

(24) Stein, C.; Herbert, J.; Head-Gordon, M. The Poisson-Boltzmann model for implicit solvation of electrolyte solutions: Quantum chemical implementation and assessment via Sechenov coefficients. *J. Chem. Phys.* **2019**, *151*, 224111.

(25) Pang, X.; Zhou, H.-X. Poisson-Boltzmann calculations: van der Waals or molecular surface? *Commun. Comput. Phys.* **2013**, *13*, 1-12.

(26) Grant, J.; Pickup, B.; Nicholls, A. A smooth permittivity function for Poisson-Boltzmann solvation methods. *J. Comp. Chem.* **2001**, *22*, 608-640.

(27) Azuara, C.; Lindahl, E.; Koehl, P.; Orland, H.; Delarue, M. Incorporating dipolar solvents with variable density in the Poisson-Boltzmann treatment of macromolecule electrostatics. *Nucleic Acids Res.* **2006**, *34*, W34-W42.

(28) Abrashkin, A.; Andelman, D.; Orland, H. Dipolar Poisson-Boltzmann equation: ions and dipoles close to charge interfaces. *Phys. Rev. Lett.* **2007**, *99*, 77801.
(29) Azuara, C.; Orland, H.; Bon, M.; Koehl, P.; Delarue, M. Incorporating dipolar solvents with variable density in Poisson-Boltzmann electrostatics. *Biophys. J.* **2008**, *95*, 5587–5605.

(30) Koehl, P.; Orland, H.; Delarue, M. Beyond Poisson-Boltzmann: modeling biomolecule-water and water-water interactions. *Phys. Rev. Lett.* **2009**, *102*, 087801.

(31) Koehl, P.; Delarue, M. Aquasol: an efficient solver for the dipolar Poisson-Boltzmann-Langevin equation. *J. Chem. Phys.* **2010**, *132*, 064101.

(32) Holst, M. Multilevel Methods for the Poisson-Boltzmann Equation. Ph.D. thesis, University of Illinois at Urbana-Champaign, USA, 1993.

(33) Holst, M.; Saied, F. Numerical solution of the nonlinear Poisson-Boltzmann equation: developing more robust and efficient methods. *J. Comp. Chem.* **1995**, *16*, 337–364.

(34) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. *Nucl. Acids. Res.* **2000**, *28*, 235–242.

(35) Dolinsky, T.; Nielsen, J.; Cammon, J. M.; Baker, N. PDB2PQR: an automated pipeline for the setup of Poisson-Boltzmann electrostatics calculations. *Nucl. Acids. Res.* **2004**, *32*, W665–W667.

(36) Sitkoff, D.; Sharp, K.; Honig, B. Accurate calculation of hydration free energies using macroscopic solvent models. *J. Phys. Chem.* **1994**, *98*, 1978–1988.

(37) Guvench, O.; MacKerell Jr, A. Computational fragment-based binding site identification by ligand competitive saturation. *PLoS Comput. Biol.* **2009**, *5*, e1000435.

(38) Oleinikovas, V.; Saladino, G.; Cossins, B.; Gervasio, F. Understanding cryptic pocket formation in protein targets by enhanced sampling simulations. *J. Am. Chem. Soc.* **2016**, *138*, 14257–14263.
(39) Heaslet, H.; Harris, M.; Fahnoe, K.; Sarver, R.; Putz, H.; Chang, J.; Subramanian, C.; Barreiro, G.; Miller, J. Structural comparison of chromosomal and exogenous dihydrofolate reductase from *Staphylococcus aureus* in complex with the potent inhibitor trimethoprim. *Proteins: Struct. Func. Bioinfo.* **2009**, *76*, 706–717.

(40) Han, G.; Lee, J.; Song, H.; Chang, C.; Min, K.; Moon, J.; Shin, D.; Kopka, M.; Sawaya, M.; Yuan, H. et al. Structural basis of non-specific lipid binding in maize lipid transfer protein complexes revealed by high-resolution X-ray crystallography. *J. Mol. Biol.* **2001**, *308*, 263–278.

(41) Silvaroli, J.; Arne, J.; Chelstowska, S.; Kiser, P.; Banerjee, S.; Golczak, M. Ligand binding induces conformational changes in human cellular retinol-binding protein 1 (CRBP1) revealed by atomic resolution crystal structures. *J. Biol. Chem.* **2016**, *291*, 8528–8540.

(42) Hamiaux, C.; Stanley, D.; Greenwood, D.; Baker, E.; Newcomb, R. D. Crystal structure of *Epiphyas postvittana* takeout 1 with bound ubiquinone supports a role as ligand carriers for takeout proteins in insects. *J. Biol. Chem.* **2009**, *284*, 3496–3503.

(43) Schrödinger, LLC, The PyMOL Molecular Graphics System, Version 1.8. 2015.

(44) Saurabh, S.; Vanaphan, N.; Wen, W.; Dauwalder, B. High functional conservation of takeout family members in a courtship model system. *PLoS One* **2018**, *13*, e0204615.

(45) Cherezov, V.; Rosenbaum, D.; Hanson, M.; Rasmussen, S.; Thian, F.; Kobilka, T.; Choi, H.; Kuhn, P.; Weis, W.; Kobilka, B. et al. High-resolution crystal structure of an engineered human beta2-adrenergic G protein-coupled receptor. *Science* **2007**, *318*, 1258–1265.

(46) Weiss, M.; Schulz, G. Structure of porin refined at 1.8 Å resolution. *J. Mol. Biol.* **1992**, *493–509*.
(47) Bai, T. Beta 2 adrenergic receptors in asthma: a current perspective. *Lung* **1992**, *170*, 125–141.

(48) Graham, S.; Leja, N.; Carlson, H. MixMD Probeview: robust binding site prediction from cosolvent simulations. *J. Chem. Inf. Model.* **2018**, *58*, 1426–1433.

(49) Tan, Y.; Verma, C. Straightforward incorporation of multiple ligand types into molecular dynamics simulations for efficient binding site detection and characterization. *J. Chem. Theory Comput.* **2020**, *16*, 6633–6644.

(50) Busch, A.; Reijerse, E.; Lubitz, W.; Frankerberg-Dinkel, N.; Hofmann, E. Structural and mechanistic insight into the ferredoxin-mediated two-electron reduction of bilins. *Biochem. J.* **2011**, *439*, 257–264.

(51) Dammeyer, T.; Frankenberg-Dinkel, N. Function and distribution of bilin biosynthesis enzymes in photosynthetic organisms. *Photochem. Photobiol. Sci.* **2008**, *7*, 1121–1130.

(52) Frankenberg, N.; Mukougawa, K.; Kohchi, T.; Lagarias, J. Functional genomics analysis of the HY2 family of ferredoxin-dependent bilin reductases from oxyhenic photosynthetic organisms. *Plant Cell* **2001**, *13*, 965–978.

(53) Cheng, K.-Y.; Lowe, E.; Sinclair, J.; Nigg, E.; Johnson, L. The crystal structure of the human polo-like kinase-1 polo box domain and its phospho-peptide complex. *EMBO J.* **2003**, *22*, 5757–5768.

(54) Liu, F.; Park, J.; Qian, W.-J.; Lim, D.; Graeber, M.; Berg, T.; Yaffe, M.; Lee, K.; Burke Jr, T. Serendipitous alkylation of a Plk1 ligand uncovers a new binding channel. *Nat. Chem. Biol.* **2011**, *7*, 595–601.

(55) Śledź, P.; Stubbs, C.; Lang, S.; Yang, Y.-Q.; MacKenzie, G.; Venkitaraman, A.; Hyvönen, M.; Abell, C. From crystal packing to molecular recognition: prediction and dis-
covery of a binding site on the surface of polo-like kinase 1. *Angew. Chem. Int. Ed.* **2011**, *50*, 4003–4006.

(56) Tirion, M. Large amplitude elastic motions in proteins from a single parameter, atomic analysis. *Phys. Rev. Lett.* **1996**, *77*, 1905–1908.

(57) Sanejouand, Y. Elastic network models: theoretical and empirical foundations. *Methods Mol. Biol.* **2013**, *914*, 601–616.

(58) Sinitskiy, A.; Voth, G. Coarse-graining of proteins based on elastic network models. *Chem. Phys.* **2013**, *422*.
Graphical TOC Entry

<figure><img src="image_13.png" /></figure>