"Martinizing" the Variational Implicit Solvent Method (VISM): Solvation Free Energy for Coarse-Grained Proteins.

Solvation is a fundamental driving force in many biological processes including biomolecular recognition and self-assembly, not to mention protein folding, dynamics, and function. The variational implicit solvent method (VISM) is a theoretical tool currently developed and optimized to estimate solvation free energies for systems of very complex topology, such as biomolecules. VISM's theoretical framework makes it unique because it couples hydrophobic, van der Waals, and electrostatic interactions as a functional of the solvation interface. By minimizing this functional, VISM produces the solvation interface as an output of the theory. In this work, we push VISM to larger scale applications by combining it with coarse-grained solute Hamiltonians adapted from the MARTINI framework, a well-established mesoscale force field for modeling large-scale biomolecule assemblies. We show how MARTINI-VISM (MVISM) compares with atomistic VISM (AVISM) for a small set of proteins differing in size, shape, and charge distribution. We also demonstrate MVISM's suitability to study the solvation properties of an interesting encounter complex, barnase-barstar. The promising results suggest that coarse-graining the protein with the MARTINI force field is indeed a valuable step to broaden VISM's and MARTINI's applications in the near future.

Introduction

Solvation, or more specifically hydration, is a fundamental driving force in essentially all biological processes. Water and ions, and their interactions, play a huge role in protein folding, dynamics, and function,[1−3] self-assembly of membranes and proteins,[4] and molecular recognition and binding,[2,5] to name only a few. Many computational methods attempt to provide an accurate description of the solvation effects underlying such processes. The most efficient among these are the so-called dielectric boundary implicit solvation methods. By representing the solvent environment as a continuum medium separated from the solute by a well-defined solvation boundary, they avoid the recurrent issue of statistical sampling of solvent configurations. Instead, the average forces and thermal fluctuations giving rise to the polarity of the solvent are implicitly represented by a uniform dielectric constant.[6]

The majority of implicit solvation methods are based on electrostatics—described by the Poisson–Boltzmann[7−9] or generalized Born theory[10−12]—with hydrophobic, van der Waals, and first-shell solvation effects taken into account by means of empirical relations between solvation free energy and accessible surface area, also known as SASA (solvent accessible surface area) methods.[6,13] While successful in many cases, their performance is ultimately limited by the fact that these methods rely on fixed and pre-established solute–solvent interfaces—normally guessed as van der Waals-based surfaces. Therefore, they fail to capture subtle though important behavior deriving from heterogeneous solvation patterns, such as polymodal hydration and dewetting,[14] which are observed in explicit solvent simulations.[4,15,16] Moreover, Poisson–Boltzmann calculations are extremely sensitive to the chosen dielectric boundary,[17] so that a poorly guessed interface can lead to very significant errors.

The variational implicit solvent method (VISM) is a solvation free energy method that avoids guessing an a priori dielectric solvation boundary. Instead, VISM expresses the solvation contributions from hydrophobic effects, van der Waals interactions, and electrostatics as a functional of all possible solute–solvent interfaces (see eq ).[18,19] VISM’s variational formulation makes it possible to identify the dielectric interface—or solvation state—that minimizes the solvation free energy, producing the optimal interface as an output of the theory. Moreover, the surface energy, van der Waals interactions, and electrostatics are not independent but coupled together through the solvation interface. In this way, sophisticated hydrophilic–hydrophobic compensation effects can be captured.[20,21]

VISM is currently developed to predict solvation free energies for systems of complex topology by means of a robust level-set method,[21−23] combined with a user-defined choice of the electrostatic formulation that varies from the Coulomb field approximation (CFA)[24,25] to the nonlinear Poisson–Boltzmann (PB) theory including ionic effects.[20] Besides predicting accurate solvation free energies of small solutes and small globular proteins, several studies have highlighted VISM’s unique ability in capturing nontrivial solvation effects such as capillary evaporation,[20,21] dry–wet transitions in ligand binding events,[25−29] dewetting of druggable binding cavities,[30] and, since it is based on a true Hamiltonian approach, being extended to solvent equilibrium fluctuations.[31]

Given the proven potential and performance of VISM in describing the complex solvation phenomena in relatively small systems, it is urgently needed now to broaden its application range toward larger scale processes, such as multiprotein association and self-assembly. For modeling consistency and a not too large scale-separation, VISM can thus be connected with efficient coarse-grained (CG) models, desirably those that are efficiently coarse but keep chemical specificity and have been developed already for the use in a coarse-grained or implicit solvent environment. A popular state-of-the art CG model of that kind is provided by the MARTINI family of models.

The MARTINI force field (FF) is a well-established coarse-grained model originally developed to simulate lipids and surfactants in lipid bilayers[32,33] and later extended to proteins and other biomolecules.[34−38] Its parametrization philosophy strongly relies on reproducing water/oil or water/membrane partition coefficients. Although the majority of MARTINI applications to proteins involve membrane-embedded or membrane-anchored proteins,[39] there are a few studies focused on water-soluble proteins.[40,41] Among these, Stark et al. found that MARTINI tends to overestimate protein–protein interactions in aqueous environment, but the discrepancy to experiments could be counteracted by a simple downscaling of the Lennard-Jones energy parameter, ε.[40] Originally, MARTINI was developed for the use with a CG water model. The issues and performance problems of the latter have been discussed and a “dry” Martini version with an efficient, fully implicit water has been introduced[42] by adapting the nonbonded parameter matrix. Still, this implicit-solvent model in dry-Martini has the same limitations as the traditional (pre-VISM) implicit solvents discussed above.

In this work, we combined the current best of two worlds and adapted VISM to work with a for our purposes amended MARTINI FF to construct an efficient and quantitative model to attack new problems with VISM, such as biomolecular assembly on large scales. As a first demonstration, we predict the solvation free energies of some selected proteins. We tested MARTINI-VISM (MVISM) against atomistic VISM (AVISM) for a set of six proteins differing in shape, size, and charge distribution. With a few adjustments and simple corrections in the Lennard-Jones parameters, MVISM is capable of reproducing atomistic solvation free energies to within 10% for 5 out of 6 proteins, besides producing reasonable electrostatic potentials at the interacting surfaces of the proteins. Finally, both MVISM and AVISM seem to correctly predict the nature of dry–wet transitions in the encounter of barnase and barstar. These results pave the way for using MVISM to address more challenging problems such as solvation dynamics and many-protein interactions, in the near future.

Theory and Methods

The Free Energy Functional

To estimate the solvation free energy of a protein (or other biomolecule) in the aqueous environment (water + ions), we start by dividing the system into three regions, as displayed in Figure . The solute region, Ωm, contains all of the N atoms belonging to the solute molecule, which are located at x1, ..., x inside Ωm and carry point charges Q1, ..., Q, respectively. The solvent region, Ωw, is treated implicitly as a continuum. The third part consists of the solute–solvent interface, Γ, which geometrically separates the regions Ωm and Ωw.

Figure 1

Schematic view of a solvation system with an implicit solvent. A solute–solvent interface, Γ, separates the solvent region, ΩW, from the solute region, Ωm. The solute atoms are located at x1, ..., x and carry point charges Q1, ..., Q, respectively.

In VISM, the solvation free energy is expressed as a functional of the solvation interface, Γ (eq ). By minimizing this functional against all possible solvation interfaces, one can find the local minima that correspond to the stable hydration states of the system.The first term in eq is purely geometrical and accounts for the hydrophobic effect by integrating the surface tension along the solvation interface. Because for systems of nanometer scale the surface tension strongly depends on the curvature, we define the local surface tension aswhere γ0 is the constant macroscopic surface tension for a planar liquid–vapor interface, H is the mean curvature defined as the average of the two principal curvatures, and τ is a curvature correction coefficient, which essentially accounts for the relative size of the solvent molecules with respect to the solute local curvature.[18]

The second term in the functional accounts for dispersion attraction and Pauli repulsion, or van der Waals interactions between solute and solvent. For each atom type of the solute, we define a pairwise interaction energy with the solvent, U. As traditionally used in molecular dynamics (MD) simulations, U is modeled by a Lennard-Jones (LJ) function (eq ):For the ith atom of the solute, located at x, the interaction energy is calculated as a volume integral over the solvent region, which is coarse-grained represented by grid points located at x, with a density prefactor, ρw.

Finally, the third term of the functional, Gelec(Γ), is the electrostatic part of the solvation free energy. In the current VISM implementation, this term can be calculated with CFA or PB theory (in its linearized or nonlinear forms). Detailed descriptions of this term are given in ref (24) (for CFA) or ref (21) (for PB).

VISM Calculations

Atomistic protein structures were obtained from the Protein Data Bank,[43] and hydrogens were added. VISM input files were generated by combining the crystallographic atomic coordinates with LJ parameters and point charges borrowed from the CHARMM36 FF[44] (Figure A). LJ parameters for water were obtained from the TIP3P model and combined with the solute parameters by using the standard Lorentz–Berthelot rules:Throughout our AVISM calculations, we fixed T = 298 K, γ0 = 0.1315 kBT/Å2, τ = 1.0 Å, ρw = 0.0333 Å–3, ϵw = 78.0, and ϵm = 2.0.[21]

Figure 2

Atomistic (A) and MARTINI (B) representation of ubiquitin, with corresponding VISM input files. For the sake of clarity, only parts of the input files are shown, with σ and ε in Å and kBT, respectively.

The electrostatic component of the free energy functional was calculated with CFA theory during the largest part of the minimization procedure and then refined in the last steps with linearized PB theory. We included +1 and −1 ionic species at a concentration of 0.1 M. We also pulled the final VISM interface closer to the solute molecules by δ = 1 Å (or δ = 2.5 Å, in MVISM calculations) before calculating the final electrostatic energies. This approach was proved to significantly improve the accuracy of electrostatic VISM energies in previous studies.[20]

VISM Calculations with MARTINI

For MVISM calculations, we submitted the atomistic PDB files to the martinize.py script with MARTINI 2.1,[34] which produced the corresponding coarse-grained PDB files. We then generated VISM input files by combining the bead coordinates from the coarse-grained PDB files with LJ parameters and point charges from MARTINI 2.1 (Figure B). Because MARTINI provides the Lennard-Jones parameters in the precombined form (U), we had to adapt the VISM algorithm to skip the LJ combination rules (see eqs and 5). We also used MARTINI’s LJ potential in its shifted form from rshift = 9 Å to rcut = 12 Å, smoothly vanishing to zero beyond rcut. Since one MARTINI water bead occupies the volume of approximately four TIP3P water molecules (data obtained from test simulations with MARTINI water and TIP3P water, at 298 K), we fixed ρw = 0.0333/4 = 0.008325 Å–3 in our MVISM calculations. Unless specified otherwise, all remaining parameters in MVISM calculations were the same used in the AVISM calculations.

In both AVISM and MVISM calculations, the grid resolution ranged from approximately 0.3 to 0.5 Å, depending on the size of the protein (∼0.3 for ubiquitin, E3, and barstar; ∼0.4 Å for barnase and BLIP; and ∼0.5 Å for β-lactamase). Each calculation was performed in a single CPU, and calculation times ranged from less than 1 h to ∼32 h depending on the system size/geometry and the choice of the initial surfaces. Calculations starting from tight initial surfaces tend to converge significantly faster than calculations that start from loose initial surfaces (check ref (21) for more details).[21]

Results and Discussion

To test and adjust MVISM performance, we chose a test set of six proteins differing in size, shape, and charge distribution (Table ). These proteins also form pairwise complexes and do not undergo large conformational changes upon binding. We started by performing MVISM calculations and comparing the resulting solvation energies with those from AVISM calculations.

Table 1

Test Set of Proteins

system	surface areaa (Å2)	heavy atoms/beads	largest cross-section (Å)	net charge	charged + polar residuesb (%)	PDB ID	reference
E3 ligase	318	350/93	∼25	–2	33.3 + 0 = 33.3	2OOB	Peschard et al.[45]
ubiquitin	413	574/156	∼25	–1	16.6 + 25.0 = 41.6	2OOB
barstar	487	699/195	∼30	–6	33.3 + 16.6 = 50.0	1BRS	Buckle et al.[46]
barnase	592	839/248	∼35	+2	50.0 + 33.3 = 83.3	1BRS
BLIP	795	1235/351	∼40	–2	23.3 + 34.6 = 60.0	1JTG	Lim et al.[47]
β-lactamase	1092	2022/550	∼50	–6	29.0 + 29.0 = 60.0	1JTG

Obtained from AVISM calculations starting from tight initials.

At the binding interfaces.

Figure reports the solvation free energy and its decomposition in surface (hydrophobic), Lennard-Jones (van der Waals), and electrostatic components, for the six individual proteins, distributed according to their size (in ascendant order). As expected, both the surface and the LJ components correlate well with protein size, as they are both dependent on the protein surface area. We also found that MVISM calculations correctly capture the overall trend in the solvation free energy and in its energy components but tend to systematically overestimate the magnitude of the LJ energies by ∼40% with respect to the atomistic calculations.

Figure 3

Solvation free energy and correspondent surface, Lennard-Jones, and electrostatic components computed with AVISM (in black) or MVISM (in blue) for six different proteins, distributed according to protein size (ascendant order).

Overestimation of LJ interactions is a known trait of the of MARTINI FF[40] and likely a deliberate way of compensating (enthalpically) for the loss of degrees of freedom—and consequent loss of entropically driven hydrophobic effects—in MD simulations with the coarse-grained model.[39] In VISM, since hydrophobic effects are taken into account implicitly in the geometric term of the free energy functional, we feel justified to downscale the MARTINI-based LJ interactions in MVISM.

Downscaling LJ Interactions by Tuning the ε Parameter

We chose to tune the LJ parameter, ε, by interpolating the well-depth of the LJ potential between (i) the value set in the original MARTINI FF and (ii) the value employed by MARTINI to describe nominally “repulsive” interactions, as originally proposed by Stark et al.[40] The degree of scaling is controlled by a scaling parameter, α, with α = 1 corresponding to the original MARTINI FF and α = 0 corresponding to the use of the weakest possible bead–bead interaction type for all bead–bead interactions (for which MARTINI uses ε = 2.0 kJ mol–1). For a given value of α, the new εα is determined usingCalculations performed with different α values reveal that a good agreement between atomistic and coarse-grained LJ energies is obtained with α = 0.5 (Figure A and Table S1). Interestingly, the significant improvement in the Lennard-Jones energies does not systematically improve the solvation free energy for all proteins, as shown in Figure B. This is a direct effect of VISM’s theoretical formulation, which allows for different energy components to couple through the solvation interface. Indeed, a more detailed analysis of each energy term reveals a compensation between the LJ and electrostatic components, such that downscaling the coarse-grained LJ interactions concomitantly allows for the coarse-grained electrostatic energies to become more negative and closer to atomistic values (Figure C). Overall, downscaling the LJ interactions in MVISM not only improves the description of the LJ energies, but it also produces a more accurate distribution of the solvation free energy among each of its energy components. We thus decided to fix α = 0.5 in all subsequent MVISM calculations.

Figure 4

Lennard-Jones (A) and solvation free energies (B) for MVISM calculations (colored lines) with different levels of ε scaling in comparison with AVISM (in black). (C) Detailed analysis of the effects of downscaling ε on each energy term. The energies depicted as dashed lines correspond to the coarse-grained energies with different levels of ε scaling. For comparison, the corresponding atomistic energies are depicted as solid lines.

The corresponding solvation surfaces produced by VISM are illustrated in Figure , along with the original molecular complexes. The coarse-grained interfaces strongly resemble the atomistic ones, although they are slightly larger and less curved due to the size and shape of MARTINI beads. Because in VISM one minimizes the free energy functional with respect to the solvation interface, the final results depend on the initial guess of the solvation surface. To account for different solvation states, it is common to start VISM calculation from tight and loose surface initials.[20,21,24,26] A loose initial corresponds to a large surface loosely encompassing all of the solute atoms, while a tight initial corresponds to a van der Waals surface of the solute atoms. In this way, VISM calculations starting from loose initials are more likely to produce “dry” states—with water expulsion near hydrophobic patches—whereas tight initials tend to produce “wet”, fully solvated states. For the six proteins tested, the solvation energies obtained from tight or loose initials are extremely similar, in both AVISM and MVISM calculations (Figure S1 and Table S2), as are the corresponding solvation interfaces (Figures S2, S3, and S4). Therefore, these proteins do not appear to display multiple hydration states in the apo-state.

Figure 5

Atomistic (magenta and violet) or coarse-grained (yellow and cyan) solvation surfaces obtained with AVISM and MVISM, respectively. These calculations started from loose initial surfaces.

The Curvature Correction Coefficient, τ

The microscopic curvature correction coefficient (τ) is a parametric coefficient used in the geometric, or hydrophobic, part of the VISM functional (see eqs and 2). It implicitly accounts for how sensitive the organization of the solvent molecules is with respect to the local curvature of the solvation surface, which is related to the solvent size. In atomistic VISM calculations, τ is typically set to 1 Å.[21] Since combining groups of atoms into beads affects the protein local curvatures, we next investigate how adjustments in τ can affect the MVISM performance with respect to atomistic calculations.

Several AVISM and MVISM calculations were performed with τ values ranging from 0.5 to 2.0 Å, using the ubiquitin and barnase proteins as test cases. As expected, we found that the τ parameter mainly affects the surface term, with small concomitant changes in the Lennard-Jones or electrostatic terms (<20 kBT, not shown). For both proteins, the surface energies decrease with increasing τ (Figure A). This is consistent with the fact that globular proteins display predominantly convex surfaces, where a larger τ decreases the local surface tension and, therefore, also the surface energy. In concave regions—such as binding pockets—a larger τ has the opposite effect of increasing the local surface tension and therefore enhancing the hydrophobic effect.[27] As a consequence, larger τ values tend to produce solvation surfaces with less pronounced pockets. This effect is illustrated in Figure , which shows the barstar binding cavity located in the barnase protein, in atomistic (B) and coarse-grained (C) representations. Comparison of the surfaces obtained with small (0.5 Å) or large (2.0 Å) τ shows that the more pronouncedly concave regions of the binding cavity become shallower with larger τ values (see, for instance, the cyan patch in Figure C). A physical interpretation for this is that a larger τ makes the binding cavities more hydrophobic, with a higher tendency to expel water (dewetting effect).

Figure 6

Effect of the curvature correction coefficient on the surface energies and solvation interfaces. (A) Surface energies with different τ values for ubiquitin (left) and barnase (right). (B) Atomistic solvation interfaces obtained with AVISM, using τ = 0.5 Å (magenta) or τ = 2.0 Å (yellow) for barnase. (C) Coarse-grained solvation interfaces obtained with MVISM, using τ = 0.5 Å (purple) or τ = 2.0 Å (cyan) for barnase.

Quantitative analysis of the surface energies shows that, in order to match the MVISM surface energies to the atomistic (AVISM) energies obtained with τ = 1 Å, one should use slightly larger τ values. In the case of ubiquitin, a good agreement is obtained with τ = 1.5 Å, whereas, for barnase, a good agreement is obtained with τ = 1.2 Å (see dotted arrows in Figure A). Since barnase is larger than ubiquitin, our results suggest that the need to adjust τ might decrease with protein size. For the remaining MVISM calculations in this work, we used τ = 1 Å, for simplicity.

Electrostatic Energies and the Shifting Parameter, δ

Besides affecting the Lennard-Jones and surface energies, coarse-graining the protein with MARTINI 2.1 also affects the electrostatic energies, since it rearranges atomic partial charges into unit integer negative or positive bead charges (+1 or −1). We next investigated to which extent coarse-graining the charges impacts the electrostatic potential at binding interfaces. Figure shows the electrostatic potential at the binding interfaces of barnase and barstar obtained with AVISM (A) or MVISM (B). Clearly, the electrostatic potential produced by MVISM is a coarse-grained one as compared to the atomistic, as would be expected, since the nonpolarizable version of MARTINI does not use partial charges. However, the overall electrostatic complementarity is maintained, with a large negatively charged patch in barstar matching a largely positive binding cavity in barnase. Similar results were obtained with ubiquitin–E3 and β-lactamase–BLIP complexes (Figures S5 and S6, respectively), which is reassuring for future applications of MVISM in protein binding and assembly.

Figure 7

Electrostatic potential at the binding interface of the barnase and barstar proteins, calculated with AVISM (A) or MVISM (B). Positive and negative regions are displayed in blue and red, respectively, with the interacting interfaces indicated by dotted circles.

Martinizing the protein also has the effect of burying the point charges deeper in the solute region than they would be in atomistic representations of the proteins, due to the larger size of the beads and the location of the charge sites in MARTINI 2.1 beads. This causes the electrostatic attraction between protein charges and the (polar) solvent medium to decrease, leading to overestimated (less negative) electrostatic energies. One way to counteract this artifact is by simply adjusting the final shift (δ) that is applied to the solvation boundary, prior to the calculation of the final PB electrostatic energies. In atomistic VISM calculations, a 1 Å shift of the final interface toward the protein has been proved to favorably account for the fact that the SASA-like VISM interface is slightly different from the dielectric boundary interface.[20,21] In MVISM, this parameter can also be tuned to counteract the deeper burial of point charges in MARTINI 2.1. We thus report the electrostatic energies obtained with MVISM for the six proteins, using different shifting values (Figure ). As we increase the shifting of the VISM interface toward the solute atoms, the electrostatic energies become more negative, reflecting the favorable interaction between charged solute groups and the aqueous environment. As shown in Figure , a good agreement between coarse-grained and atomistic electrostatic energies is systematically obtained by performing MVISM calculations with a slightly larger shift of 2.5 Å.

Figure 8

Electrostatic energies calculated with MVISM using different shifting (δ) values. As a reference, atomistic electrostatic energies obtained with δ = 1 Å are indicated by pink dotted lines.

It is worth noting that there is still room to further improve MVISM electrostatic energies (and corresponding electrostatic potentials) by employing more sophisticated versions of the MARTINI force field. The polarizable MARTINI 2.2P, for instance, can improve the electrostatic description of polar proteins by including auxiliary partially charged sites for polar residues.[48] Whether MVISM calculations of highly polar proteins such as barnase could benefit from the polarizable MARTINI 2.2P remains to be tested, but it is certainly an interesting direction should one need more accurate electrostatics in MVISM applications.

The Parametrized MVISM Method

After parametrization, we suggest applying the MVISM method with (i) σ and (2-fold downscaled) ε Lennard-Jones parameters from the MARTINI 2.1 FF for the van der Waals part of the free energy functional; (ii) τ = 1 Å for the surface term in the free energy functional (though one could use a larger τ to get slightly more accurate agreement of the surface energies); and (iii) applying a slightly larger shift of the solvation boundary (δ = 2.5 Å) prior to the PB electrostatic calculations.

An Application Example

One of the most interesting aspects of the VISM method is its ability to capture different solvation states at regions where an accurate solvation description is critical, as in binding pockets or interfaces. As an application example, we decided to test whether MVISM can capture dry–wet transitions during a molecular encounter of the barnase and barstar proteins. Barstar and barnase form a tight complex which has achieved its extremely fast kinetics of binding by means of optimized electrostatic interactions.[49] As such, both barnase and barstar binding interfaces are rich in charged residues and polar residues, as shown in Table .

To create a reasonable encounter pathway for this complex, we separated the two proteins along the axis formed by their geometrical centers. Configurations were saved every 2 Å, with separation distances ranging from 0 (native bound complex) to 15 Å (unbound proteins). For each configuration, we performed VISM calculations starting from both tight and loose initials. Besides the solvation energies, we kept track of the resulting surfaces, from which we could see whether the calculations produced “wet” or “dry” encounter complexes.

Figure A displays the energetic profiles obtained for the association between barnase and barstar, starting from tight or loose initials. We found that AVISM calculations converge at separation distances <3 Å and >10 Å regardless of the initial surfaces. At less than 3 Å, AVISM calculations necessarily produce “dry” states because water cannot fit in the space between the interaction interfaces, being squeezed out. At separation distances larger than 10 Å, AVISM always produces “wet” states because the hydrophobic cost of producing large concave surfaces in the space between the proteins becomes too high. In MVISM calculations, convergence between tight and loose calculations is achieved at slightly larger separation distances. This is likely due to the larger size of water beads—which are squeezed out at separations smaller than 5 Å—and due to less favorable electrostatic interactions with water—which allow for dewetting to occur at separations as large as 11.4 Å.

Figure 9

Dry–wet transitions in the encounter complex of barstar and barnase. (A) AVISM (blue) and MVISM (purple) solvation energies obtained from tight (triangles) or loose (circles) initial surfaces. (B) Different solvation states captured by AVISM at 7.6 Å. (C) Different solvation states captured by MVISM at 11.4 Å.

At intermediate separation distances, however, both AVISM and MVISM capture distinctive “dry” or “wet” solvation states depending on the initial surface shape. This is illustrated by the solvation surfaces displayed in Figure B (atomistic) and C (coarse-grained), and also evident from the hysteresis in the energy profiles in Figure A. Interestingly, partial desolvation has been suggested to occur during the binding encounter of barstar and barnase, at separation distances up to 7 Å.[50] It has been argued that, without the solvent to screen the electrostatic interactions between the two proteins, strong long-range electrostatic interactions pull the two domains together. This effect would contribute favorably to the fast association of the proteins, even if destabilizing the final thermodynamics of binding.[51] Still, both AVISM and MVISM calculations reveal that the “wet” encounter pathway is significantly favored by approximately 150–180 kBT with respect to the “dry” pathway, which makes sense as this is a highly hydrophilic complex. Not surprisingly, the difference in free energy is dominated by the electrostatic component, since the interacting surfaces are highly charged and thus display very favorable electrostatic interactions with water, favoring the “wet” states.

Overall, MVISM seems to correctly predict the nature of dry–wet transitions in the encounter between barnase and barstar. Combining the VISM model with coarse-grained molecular dynamics simulations is one of our ultimate goals to study protein–protein interactions. This application example thus serves to highlight the potential of using MVISM to account for solvation effects in future coarse-grained simulations with implicit solvent.

Conclusions

In this work, we adapted the VISM method to work with the MARTINI 2.1 coarse-grained force field. Using atomistic VISM energies as a reference, the main adjustment of the MARTINI Hamiltonian to produce VISM solvation energies consisted of tuning the LJ parameter ε to correct for overestimated van der Waals interactions. MARTINI’s overestimation of LJ interactions arises in part as an enthalpic compensation to the loss of hydrophobic interactions due to reduced degrees of freedom of coarse-grained water. In VISM, such compensation is not required because hydrophobic interactions are accounted for implicitly in the surface term of the free energy functional.

After downscaling of LJ interactions, MVISM displayed good agreement with atomistic (AVISM) calculations, not only capturing the trends in solvation free energy of six different proteins but also correctly partitioning the free energy into hydrophobic, van der Waals, and electrostatic components. Despite the use of simplified coarse-grained charges, MVISM is capable of producing reasonable electrostatic complementarity between binding interfaces, which is a necessary feature for subsequent protein–protein binding studies. In terms of electrostatics, MVISM could be further improved with the use of more sophisticated versions of MARTINI, such as the polarizable MARTINI 2.2P.[48]

In the present formulation, MVISM can qualitatively capture and reproduce the dry–wet transitions in the binding of barnase and barstar, as observed in atomistic VISM calculations and supported by previous theoretical studies.[50] As such, the results reported herein are the first step toward the development of a hybrid approach combining the coarse-grained simulations of solute atoms and our VISM description of completely implicit solvent. We expect this work to broaden VISM applicability to study complex solvation properties of important protein assemblies in the near future.