Modeling and analysis of ensemble average solvation energy and solute–solvent interfacial fluctuations

2.1 Background

We consider a solute–solvent system with a fixed biomolecular structure contained in a bounded Lipschitz domain Ω ⊂ R 3 by using a (statistical) grand canonical ensemble. The temperature, chemical potentials, and volume of the system are kept constant. Suppose that there are N c ion species in Ω and the system contains N c + 2 types of particles, i.e., the solute and solvent molecules and N c ion species. We assume that the solute atoms are centered at x 1 , … , x N m ∈ Ω .

Define a probability space ( S , ℱ , P ) , where the sample space S = { S α : α ∈ A } denotes the set of all possible states of the system with A being the index set of the states. ℱ is a σ -algebra of S and P is the probability measure. Because a biomolecular structure is fixed, each state S α is uniquely determined by a function:

where C α , j ( x ) is the (number) concentration of the jth ion species at x for j = 1 , … , N c , or the (number) concentration of solvent molecules for j = 0 . In state S α , a point x ∈ Ω is in the solute phase if C α ( x ) = ( 0 , … , 0 ) , otherwise in the solvent phase. A point that is not occupied by any particle is typically considered to be in the vacuum. However, due to the similar dielectric constants of the solute and the vacuum, it is reasonable to assume that such points can be treated as being in the solute phase.

It is important to note that different states can result in the same solute and solvent phases. For instance, consider the scenario where a solvent molecule and an ion are interchanged, which are originally located at different positions, within a given state. Allowing for a slight abuse of notation, we refer to the subset of states within S that share the same solute and solvent phases as a microstate of the solute–solvent system. Assume that the system undergoes K microstates: ℳ 1 , … , ℳ K ∈ ℱ ; and each ℳ k occurs with a probability p k = P ( ℳ k ) . For notational brevity, we put K = { 1 , 2 , … , K } . This leads to a decomposition S = ⨆ k ∈ K ℳ k ; and thus ∑ k ∈ K p k = 1 .

We define several random variables that will be extensively used in this article.

• X : S → K , where X = ∑ k ∈ K k χ ℳ k . Thus, X ( S α ) indicates the microstate that S α belongs to.

• C j : S → L 1 ( Ω ) , where C j ( S α ) = C α , j with j ∈ { 0 , 1 , … , N c } . Thus, when j ∈ { 1 , … , N c } , C j ( S α ) ( x ) is the (number) concentration of the jth ion species at x ∈ Ω in S α ; C 0 ( S α ) ( x ) is the (number) concentration of the solvent molecule at x ∈ Ω in S α .

Similarly, if one considers ℳ k as a (statistic) ensemble, the ensemble average of Y in ℳ k is defined by

Throughout the remainder of this article, unless explicitly stated otherwise, the ensemble average of a quantity Y always refers to the object defined by (2.1).

2.2 A sharp interface VISM

In this section, we will introduce the solvation energy defined in a microstate ℳ k with a fixed sharp solute–solvent interface. Assume that in microstate ℳ k , the solute phase is represented by a Caccioppoli subset D k ⋐ Ω , or equivalently, U k = Ω \ D ¯ k is the solvent phase. Then the solute–solvent interface is given by ∂ D k . See Figure 1(a) for an illustration.

Figure 1

(a) Illustration of a microstate ℳ k : D k is the solute region; U k is the solvent region; and ∂ D k is the solute–solvent interface. (b) Domain decomposition for a grand canonical ensemble: Ω i is the region occupied by the solute in all microstates; Ω e denotes the region occupied by the solvent in all microstates; and Ω t denotes the transition region where 0 ≤ u ≤ 1 .

In this article, solute atoms are treated as hard spheres. More precisely, we consider the atom centered at x i , i = 1 , … , N m , as a sphere of radius σ i > 0 , i.e., B ( x i , σ i ) , for which solvent molecules and ions cannot enter. On the other hand, the solvent and ion concentrations are the same as their bulk concentrations in regions sufficiently far away from the solute. As a consequence, such regions cannot be contained in the solute phase in any microstate. On the basis of these observations, we may assume that there are two open subsets, Ω i and Ω e , of Ω with ⋃ i = 1 N m B ¯ ( x i , σ i ) ⊂ Ω i ⋐ Ω \ Ω ¯ e and ∂ Ω ⊂ ∂ Ω e such that

For instance, one can take Ω i to be the region enclosed by a smoothed solvent excluded surface and Ω e to be region outside a smoothed solvent accessible surface. Let Ω t ≔ Ω \ ( Ω ¯ i ∪ Ω ¯ e ) be the transition region. By Assumption (2.3), in all microstates,

• ions are located outside Ω i , and

• the solute–solvent interfaces are located inside Ω ¯ t .

The hydrophobicity of the amino acids in the biomolecule varies from position to position. To account for this phenomenon, we introduce a positive variable surface tension function θ ∈ C 1 ( Ω ¯ t ) . Without loss of generality, we may extend θ to be a function in C 1 ( Ω ¯ ) , still denoted by θ , such that

for some constants 0 < θ 0 < θ 1 . By [10, Lemma 1],

See Definition C.1 for the definition of B V -functions. It is clear that, for any f ∈ B V ( Ω ) with f ≡ 0 in Ω e and f ≡ 1 in Ω i , the value of ∫ Ω θ ( x ) d ∣ D f ( x ) ∣ is independent of how θ is extended outside Ω ¯ t .

As proposed in [24,25], the solvation free energy in microstate ℳ k predicted by a sharp-interface VISM is the value of the following energy functional evaluated at χ D k :

where ψ : Ω → R is the electrostatic potential. The two components I np and I p are termed the nonpolar and polar portion of the solvation energy, respectively. For every fixed D k , ψ = ψ χ D k is chosen to maximize I p ( χ D k , ψ ) among all ψ taking a predetermined Dirichlet boundary value ψ D on ∂ Ω .

The nonpolar solvation energy consists of three parts:

When θ ≡ γ for some constant γ > 0 , the first term reduces to γ Per ( D k ; Ω ) with Per ( D k ; Ω ) being the perimeter (in Ω ) of the biomolecule region D k . This term measures the disruption of intermolecular and/or intramolecular bonds that occurs when a surface is created. In addition, Vol ( D k ) represents the volume of D k ; P h is the (constant) hydrodynamic pressure. Therefore, P h Vol ( D k ) is the mechanical work of creating the biomolecular size vacuum in the solvent. In the last integral, ρ s is the (constant) solvent bulk density; and U vdW represents the Lennard-Jones potential [57]; as such U vdW ∈ C ∞ ( Ω ¯ \ { x 1 , … , x N m } ) .

Currently, one of the most widely used polar solvation models is based on the PB theory [1,28,29,35,40]. In the framework of the classical PB theory, the polar energy is expressed as follows:

where ρ is an L ∞ -approximation of the solute partial charges supported in ⋃ i = 1 N m B ( x i , σ i ) . ε p and ε s are the dielectric constants of the solute and the solvent, respectively. Usually, ε p ≈ 1 for the protein and ε s ≈ 80 for the water. q j is the charge of ion species j , j = 1 , 2 , … , N c , and β = 1 ⁄ k B T , where k B is the Boltzmann constant, T is the (constant) absolute temperature. c j ∞ is the (constant) bulk concentration of the jth ionic species. For notational brevity, we define

In addition, we assume the charge neutrality condition

It is important to observe that B ( 0 ) = 0 and, by (2.7), B ′ ( 0 ) = 0 and B ′ ( ± ∞ ) = ± ∞ . Further, B ″ ( s ) > 0 . We thus conclude that B ( 0 ) = min s ∈ R B ( s ) and B is strictly convex.

To sum up, in microstate ℳ k , the solvation energy is given by

where ψ = ψ χ D k solves the PB equation:

in the admissible set

The boundary data ψ D is fixed throughout this article.

3.1 Classic PB theory

In this subsection, we will first consider the classic PB theory, in which ions and solvent molecules are assumed to be point-like and have no correlation between their concentrations.