Bridging solvent molecules mediate RNase A - Ligand binding.

Due to its high catalytic activity and readily available supply, ribonuclease A (RNase A) has become a pivotal enzyme in the history of protein science. Moreover, this great interest has carried over to computational chemistry and molecular dynamics, where RNase A has become a model system for various types of studies, all the while being an important drug design target in its own right. Here, we present a detailed molecular dynamics study of RNase-ligand binding involving 22 compounds, spanning nearly five orders of magnitude in affinity, and totaling 8.8 μs of sampling with the standard Amber parameters and an additional 8.8 μs of sampling with a modified potential. We show that short-lived, solvent-mediated bridging interactions are crucial to RNase-ligand binding. We characterize the behavior of bridging solvent molecules, uncovering a power-law dependence between the lifetime of a solvent bridge and the probability of its occurrence. We also demonstrate that from an energetic perspective, bridging solvent in RNase A-ligand binding behaves like part of the enzyme, rather than the ligands. Moreover, we describe an automated pipeline for the detection and processing of bridging interactions, and offer an independent assessment of the performance of the state-of-the-art fixed-charge force fields. Thus, our work has broad implications for drug design and computational chemistry in general.

Introduction

Ribonuclease A (RNase A) is a member of the endonuclease family of enzymes that degrade RNA. It and its homologs are involved in a multitude of functions in both the healthy and diseased state [1,2]. Some of the more prominent of these homologs are the eosinophil-derived neurotoxin (RNase 2), and eosinophil cationic protein (RNase 3), which have been implicated in hypereosinophilic and allergic conditions, and angiogenin (RNase 5), which induces neovascularization during normal organ growth, cancer, and in vascular and rheumatoid disorders [3]. Furthermore, RNase A has been shown to have the highest catalytic activity of all known family members. In part due to its abundance and high catalytic activity, bovine pancreatic ribonuclease became the first enzyme to have its catalytic mechanism described [4]. It is, therefore, unsurprising that bovine pancreatic ribonuclease A has been the subject of many pioneering studies in protein chemistry and enzymology, and is commonly used in laboratories. Moreover, apart from experimental studies, it has also served as a model system for molecular dynamics investigations focusing on enthalpy-entropy compensation, ligand binding entropy, protein engineering for thermal stability, and the role of protein dynamics in enzymatic catalysis [5]. Finally, RNase A and its homologs have become important targets for drug design [6,7]. Previous crystallographic studies [8] have suggested bridging water molecules [9,10] are involved in RNase–ligand binding, possibly forming a network of water-mediated interactions [11]. The extent of solvent involvement in RNase—ligand binding, however, remains uncertain.

Here, we present a molecular dynamics investigation of the binding between full-length bovine pancreatic RNase A and 22 nucleotide, nucleoside, and pyrophosphate-containing ligands, spanning nearly five orders of magnitude in affinity (Fig 1; note that in the figure phosphate groups are drawn protonated, whereas they have been simulated in a fully deprotonated state). We show that bridging water molecules and ions are an essential component of RNase A–ligand binding, and characterize this interaction at resolution and in detail that is unavailable to structural studies. We demonstrate that a far greater number of water molecules mediate the binding process than what is visible in published RNase–ligand structures [6-15], and characterize their behavior and contribution to the binding process. Our work involves extensive sampling–nearly 10 μs with the standard Amber parameters and nearly 10 μs with a modified potential. This serves as an independent test of the latest version of the general Amber force field–GAFF 2.11 [16]. Moreover, our work has broad implications for drug design campaigns targeting RNases, as well as for the broader field of computer-aided drug design (CADD) [17] in general.

Fig 1

Ligand molecules examined in the present study and their maximum common substructure.

Below each ligand, we give its name and its pK value for RNase A. All pK values are pKis, except for PAX, which is a pKd. Therefore, throughout the present report, we refer to these values simply as pKs. Note that phosphate groups are drawn protonated, whereas they have been simulated in a fully deprotonated state.

Methods

System preparation

Initial coordinates for the complexes between RNase A and 22 different ligands (Fig 2) were obtained from the 2018 refined set of the PDBbind database [18,19]. Water molecules up to 4 Å away from ligand atoms were retained to ensure that all waters mediating protein–ligand interactions from the published structures are included in the simulations. As all 22 structures included the full-length protein with no chain breaks, caps were not needed. Protein chains were protonated and solvated in a truncated octahedral box with TIP3P water [20] with the tleap program from the Amber18 package with a minimal wall distance of 15 Å, for a total system size of around 30 000 atoms. 0.150 M of NaCl with Joung and Cheatham parameters [21] were added to approximate a physiological salt concentration while ensuring charge neutrality. Ligand protonation states and net charges were set to those distributed by the PDBbind curators; only for ligands with an odd number of bonding electrons were charges manually adjusted to their correct values, corresponding to a pH of 7. Ligand parameters were obtained using the general Amber force field (GAFF 2.11) [16] with AM1-BCC charges [22] using antechamber.

Fig 2

The 22 starting protein–ligand crystal and NMR structures aligned by Cα.

The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and ligands are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. The view is centered on the ligands in the binding groove of the enzyme. Also labeled are several lysine and arginine residues, which are further discussed in the text.

Molecular dynamics simulations

The solvated systems were subjected to a series of 1000 steps of energy minimization with the steepest descent method and harmonic restraints of 3 kcal*mol−1 *Å-2 applied to all heavy atoms of the solutes (ligand and protein), followed by 1000 steps of conjugate gradient minimization with identical restraints. The systems were heated from 0 to 300 K over a period of 1 ns at constant volume with 3 kcal*mol−1 *Å-2 harmonic restraints on solute heavy atoms, followed by 1 ns of constant pressure density equilibration with restraints. The systems were then equilibrated for 1 ns without any restraints and simulated for 100 ns of production dynamics under constant pressure (1 bar) and temperature (300 K) conditions, maintained with the Berendsen barostat [23] and the Langevin thermostat, respectively [24]. Collision frequencies for temperature coupling were set to 2 ps−1; the pressure relaxation time was set to 2 ps, using isotropic position scaling. All systems were simulated in four independent replicas with the ff14SB force field under periodic boundary conditions [25]. A 12.0 Å cutoff was used for both van der Waals and electrostatic interactions; long-range electrostatics beyond the real-space cutoff were computed with the particle-mesh Ewald (PME) scheme [26]. During heating, density equilibration, preproduction, and production dynamics, bonds to hydrogen were constrained using the SHAKE algorithm [27], allowing for a 2 fs time step; only during energy minimization bonds to hydrogen were not constrained. During production dynamics, frames were saved every 100 ps (0.1 ns) for a total of 1000 per trajectory, to be used in subsequent analysis.

Trajectory processing, identification of bridging solvent, and residence time analysis

Rototranslational alignment, as well as ligand heavy atom and protein carbon alfa (Cα) root-mean-square deviation (RMSD), and root-mean-square fluctuation (RMSF) calculations (with respect to the starting coordinates) were performed with cpptraj V4.14.0 [28]. For each production trajectory, bridging solvent molecules and ions were identified with the hbond command for cpptraj using the nointramol keyword, which excludes intramolecular bridging waters and ions from consideration. The distance and angle cutoffs for hbond were 3 Å and 135 degrees, respectively. After all intermolecular bridging agents were identified, they, along with the protein and ligand, were stripped into a separate trajectory file. Corresponding topologies for these trajectories were generated with the ante-MMPBSA.py utility of AmberTools19.

For each production dynamics run, the number of bridging solvent molecules and ions was calculated, as were the total number of bridging interactions formed during the 100 ns simulations, and the residence times of the solvent molecules. Herein, we define a bridging water molecule as one that is simultaneously hydrogen bonded to the protein and ligand in at least one of the 1000 frames per production dynamics run. When counting the number of bridges per trajectory, we defined each bridge as a single, continuous interaction. For example, if a given water molecule is involved in bridging interactions in frames 300, 390, 391, 392, 911, and 912, we record three separate bridges with lifetimes of 0.1, 0.3, and 0.2 ns, respectively. We stress that a given water molecule or ion can be involved in multiple bridging interactions. Hence, the number of bridges in a molecular dynamics simulation is typically much larger than the number of bridging molecules. For example, an ion may become involved in a single, long-residence time bridging interaction where it remains located between the ligand and protein for several nanoseconds without interruption. Then, it can leave this location and return to it multiple times for short periods of time. Hence, we record one bridging ion and several bridging interactions–one long-lived and multiple short-lived bridges. Due to the ambiguity in defining “residence time” throughout the literature, we point out that in the present report, we use the terms “lifetime of a bridging interaction” and “residence time of a water molecule/ion” interchangeably.

To gain further insight into the global behavior of the bridging waters and ions over the course of the 400 ns of production dynamics for each complex, we generated their radial distribution functions (RDFs) around the ligands with the Gromacs rdf tool. We have performed RDF analysis for bridging water and bridging ions separately.

MM-PB(GB)SA calculations and analysis

For each complex trajectory, the enthalpy of interaction between the protein and the bound ligand (ΔH) was computed with the MMPBSA.py script [29], part of the Amber18 package. Briefly, molecular mechanics—Poisson-Boltzmann surface area (MM-PBSA) calculations are an implicit solvent, end-state free energy method [30]. A thermodynamic cycle is constructed and the energy of interaction is evaluated as the difference between the complex free energy and the sum of the individual components' free energies: where the angle brackets denote an ensemble average of snapshots taken from a (usually explicit solvent) trajectory [31]. For each component, ΔG is evaluated as where Egas is a gas phase energy, calculated from the trajectory in accordance with the force field parameters, ΔGsolvation is the solvation free energy, computed with an implicit solvent model, and TS is the entropic contribution to binding, estimated via normal mode analysis or with the quasiharmonic approximation [29,32]. The solvation energy, in turn, comprises an electrostatic and nonpolar component: The former is computed by modeling the solute as a set of spheres with appropriate charges and radii, embedded in a structureless continuum (the solvent and ions dissolved in it) and numerically solving the Poisson-Boltzmann equation [33]. Poisson-Boltzmann calculations were performed using the internal PBSA solver in sander. The nonpolar part was computed as a sum of two terms–a repulsive term stemming from the formation of a cavity in the solvent, which accommodates the solute, as well as from repulsive solute–solvent interactions, and an attractive term, arising from favorable solute–solvent interactions [34].

It has been shown that the attractive term can be approximated by the van der Waals interactions between solute and solvent [35,36] which are modeled with the Lennard-Johnes 6–12 potential. Finally, the repulsive component of the nonpolar free energy of solvation is obtained via a simple linear relationship with the solvent accessible surface area (SASA) of the solute: where γ is a surface tension coefficient and c is a cavity offset corresponding to ΔGrepulsive for a solute of zero volume [34].

In our subsequent analysis, we also included the similar in spirit but computationally cheaper molecular mechanics—generalized Born surface area calculations (MM-GBSA). Like MM-PBSA, MM-GBSA is also an end state, post-processing method for evaluating binding free energies. Here, polar solvation energies are computed from pairwise summations over charge–charge interactions, scaled in accord with effective atomic burial or the “Born radius” of the atoms [37]. Born radii were obtained from the work of Onufriev et al. [38] by setting the igb parameter to 5. SASA is calculated with the linear combinations of pairwise overlaps (LCPO) method from atomic radii [39].

MM-PBSA and MM-GBSA calculations were performed using mbondi2 radii and default settings for the nonpolar decomposition scheme, surface tension, cavity offset, and external and internal dielectric constants. We only adjusted the default setting for the ionic strength (0.0 M) to the one used during production dynamics– 0.150 M.

MM-PBSA and MM-GBSA calculations were performed on the complex trajectories (the so-called “one trajectory approach”) on the complexes between the protein and the ligands, as well as complexes including the protein, ligand, and all bridging waters and ions.Per-residue energy decompositions were also performed for every individual frame, adding 1–4 energy terms to internal energy terms; the energies for every residue were averaged over the 1000 frames of production dynamics. This allowed us to assess each individual residue’s contribution to binding over the course of an entire trajectory–favorable, unfavorable, or indifferent.

As we were primarily interested in relative binding energies, rather than their absolute values, and bearing in mind the considerable approximations and inaccuracies involved in computing entropy [30,40,41], as well the considerable time and computational effort required, we chose to omit entropy calculations from our analysis. Given that entropies were not explicitly calculated, i.e. the entropy term has not been included in Eq (2), our MM-PB(GB)SA calculations produced the enthalpy of binding (ΔH), which can be calculated rather reliably, as opposed to the free energy (ΔG) and entropy (ΔS), which cannot.

Results

Complex stability

The protein components of the complexes remained stable throughout the 100 ns dynamics runs in all replicas. The mean Cα RMSDs typically varied below 2.5 Å, with only individual trajectories exceeding 3 Å (see S1 Fig and S1 File). A more detailed analysis of the proteins was performed by examining the individual Cα fluctuations (S2 Fig and S1 File). This analysis revealed that the proteins were highly stable, with only unstructured regions, primarily the N-terminal loop, exhibiting high fluctuations (see Figs 2 and S2 and S1 File). The ligands exhibited somewhat varied behavior–most were stably bound, as assessed by their heavy atom RMSDs, although some exhibited greater mobility. This pertains to some of the smaller compounds, which tended to explore more of the relatively broad and shallow binding groove of the RNase (see, for example, Fig 3A and 3B), as well as the large, pyrophosphate (also known as diphosphate)-containing compounds, which reorient themselves to make extensive phosphate–lysine/arginine and pyrophosphate–lysine/arginine interactions with the protein (Fig 3C and 3D). As lysine and arginine side chains bear terminal amino groups, these interactions are of the amine–phosphate type. Most (pyro)phosphate–amine interactions involved K7, R10, R39, and K41 (Fig 3C and 3D), which form the binding groove of the enzyme. During dynamics, the amine–phosphate contacts often form clusters of salt bridges and salt-linked triads, where two salt bridges share a residue in between [42,43] (see the salt-linked triad featuring K7 in Fig 3D). Such a rearrangement is also observed with the small, phosphate-containing compounds where the crystal structures have little or no direct amine–phosphate contacts (see Fig 4A and 4C). After 100 ns of explicit solvent dynamics, the compounds reorient themselves to make extensive amine–phosphate contacts. In the C5P simulations, for example, the ligand becomes involved in salt-linked triads with R39 and K41 (compare Fig 4A and 4B). An even greater clustering occurs in the A3P simulations, as this compound has two phosphate groups, and forms a cluster of amine-phosphate contacts with K7, K10, K47, and K41 (compare Fig 4C and 4D).

Fig 3

Comparison between crystal and simulated structures.

(A) The RNase A–C3P complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray, the ligand is given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. (B) The RNase A–C3P complex after 100 ns of production dynamics. (C) The RNase A–PAX complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled; polar contacts are shown as dotted lines. (D) The RNase A–PAX complex after 100 ns of production dynamics.

Fig 4

Comparison between crystal and simulated structures.

(A) The RNase A–C5P complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. (B) The RNase A–C5P complex after 100 ns of production dynamics; polar contacts are shown as dotted lines. (C) The RNase A–A3P complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled; polar contacts are shown as dotted lines. (D) The RNase A–A3P complex after 100 ns of production dynamics.

Solvent-mediated bridging interactions

Trajectory processing revealed that in each protein–ligand system, hundreds to thousands of different solvent molecules participated in at least one bridging interaction (Fig 5A). In the majority of trajectories, this corresponded to bridging interactions numbering in the thousands. Phosphate-containing compounds tended to attract more bridging waters than compounds that did not have phosphate groups. Furthermore, the larger, pyrophosphate-containing compounds tended to attract more bridging waters than the small compounds that have only one phosphate group (see Figs 1 and 5A). Moreover, examining the composition of the bridging agents showed that no chloride ions are involved in bridging interactions, i.e. only water and sodium ions mediate RNase–ligand binding. We then examined the duration of these bridging interactions. It was found that the vast majority of bridges (98 to over 99%) have rather short lifetimes–below 1 ns. However, among the different RNase–ligand systems, we also observed tens of bridges with lifetimes between 1 and 10 ns, and a small number of bridges, single-digit numbers for most compounds, that lasted 10 or more nanoseconds. The three longest-lasting water bridges in our data set were observed in the PUA–RNase, C3P –RNase, and PAX–RNase systems, and had lifetimes of 55.5, 35.6 and 23.3 ns, respectively.

Fig 5

Average numbers of bridging water molecules and ions, and bridging water–ligand and bridging Na+—ligand RDFs.

(A) For every compound, we give the average number of bridging molecules (blue columns) and bridges (red columns) for the four production runs; error bars represent standard deviations. (B) Bridging water–ligand radial distribution functions computed from the four 100 ns production dynamics simulations for every compound. (C) Bridging Na+–ligand radial distribution functions computed from the four 100 ns production dynamics simulations for every compound.

Furthermore, we examined the radial distribution functions of the bridging waters and bridging Na+ around the ligands throughout the production dynamics runs. The radial distribution function provides information how likely it is for a molecule or chemical group to be found at a given distance from another given group or molecule, as opposed to finding it in bulk solvent, where RDF = 1. In other words, RDF values below 1 indicate that a given molecule is less likely to be found around a given molecule at a given distance compared to being found in bulk solvent; RDF values above 1 indicate the opposite. Examining the RDFs for bridging waters and bridging Na+ reveals an interesting pattern. It is evident that a bridging water molecule is more likely to be found in bulk solvent (Fig 5B and S1 File), rather than within 3 Å of the ligands, mediating their interactions with the protein. Conversely, for most ligands, sodium ions are more likely to be found in the vicinity of the compounds, rather than bulk solvent (Fig 5C and S1 File). Notably, the Na+- ligand radial distribution functions for the pyrophosphate (or equivalently, diphosphate) compounds have high peaks, implying that bridging sodium ions are tens of times more likely to be found near the ligand, rather than in bulk solution (see Figs 1 and 5 and S1 File). Conversely, most of the compounds whose Na+ RDFs do not exhibit peaks lack phosphate groups (see Figs 1 and 5 and S1 File).

We monitored the MM-PBSA and MM-GBSA energies over the course of the production dynamics runs to assess their convergence. When bridging solvent is excluded from the computations, most enthalpies converge (i.e., reach a plateau) well before the end of the simulation, around the 50 ns mark (S4 Fig and S1 File). Including bridging solvent in the calculations makes the systems larger and more complex. Correspondingly, most of the energies take longer to converge (S5 Fig and S1 File); this is particularly the case with the larger, pyrophosphate-containing compounds, e.g. PAP, PUA, and PAX. We note that a pyrophosphate group, known also as a diphosphate group, is simply two phosphate groups condensed together. These compounds tended to attract the greatest number of bridging solvent molecules, numbering in the thousands (see also Figs 1 and 5A). Nevertheless, there is no significant difference between results performed on the entire production runs and results derived only from the latter half of each trajectory (see the data in S1 File). Therefore, for completeness, from this point forward, we present and discuss only the former.

When excluding bridging particles from the MM-PBSA and MM-GBSA calculations and analysis, i.e. performing calculations on only the protein and the ligand, no correlation is observed between calculated ΔH values and experimental affinities, expressed as pK (Fig 6 and S1 File). Moreover, the lines of best fit for the ΔH/pK correlations from the four replicas cross each other at large angles. Conversely, when including bridging molecules in the calculations, the correlations become considerable, with the slopes of the lines of best fit from the four replicas becoming much more similar. A more detailed inspection of the ligands in Fig 1, the plots in Fig 6C and 6D, and the molecular dynamics trajectories reveals that the greatest outliers are phosphate or pyrophosphate containing molecules whose interactions with the protein are dominated by (pyro)phosphate–lysine or (pyro)phosphate–arginine interactions. These compounds have greatly overestimated affinities for the enzyme, noticeably eroding R2.

Fig 6

ΔH/pK correlations.

(A) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-PBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black); bridging waters and ions are not included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color. (B) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-GBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black); bridging waters and ions are not included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color. (C) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-PBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black); bridging waters and ions are included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color. (D) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-GBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black); bridging waters and ions are included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color.

As PAX follows the same general ΔH/pK trend as the other ligands, we conclude that its pKd value is directly comparable to the pKi values of the remaining compounds.

Discussion

We chose to utilize the PDBbind database, as it contains a large, manually-curated compilation of high-quality protein–ligand structures. Moreover, the refined set contains only pKi and pKd values, which are directly comparable to each other, unlike pIC50 values, which depend critically on experimental conditions and can only be compared across identical assays. The PDBbind database contains only molecules consisting of C, N, O, P, S, F, Cl, Br, I, and H atoms and excludes ligands with unusual chemistries, e.g. compounds containing Be, B, Si or metals. Thus, it is particularly well suited to the needs of drug design.

As in previous work [44,45], we decomposed ΔGnonpolar into a dispersive (attractive) and cavitation (repulsive) term [34] in the MM-PBSA calculations, as this scheme was shown to provide much better agreement between computed [44] and isothermal calorimetry results [46] for ΔH. The alternative scheme, where ΔGnonpolar is linearly dependent on SASA, has been shown to grossly overestimate ΔH [44].

The current implementation of MM-PB(GB)SA calculations in MMPBSA.py [29] requires that the receptor and ligand be explicitly defined. When including bridging waters and ions in the calculations, this creates an ambiguity–one can define these particles either as part of the receptor or the ligand. Previous work on a diverse range of systems, including topoisomerase–camptothecin derivatives, α-thrombin–benzothiophene and benzopiperidine derivatives, penicillopepsin–peptide, penicillopepsin–naphthalene derivative complexes, and avidin–biotin analogues [47], has shown that computed ΔH correlates with affinity to a significant degree only when the bridging solvent is treated as part of the receptor. Our results on the complexes between bovine pancreatic RNase and 22 nucleotides, nucleosides, and their analogs corroborate this finding. Indeed, the R2 values presented in Fig 6C and 6D were obtained by including bridging water molecules and ions in the receptor; including them in the ligand produces no significant correlation (see the data in S1 File). It is, therefore, evident that this behavior is consistent across different systems, rather than being system-dependent. Moreover, this behavior is perhaps unexpected and certainly warrants rationalization. High-resolution structural studies on the oligopeptide binding protein (OppA) complexed to different tripeptides of the KXK type, where only the middle residue is allowed to vary, show that different side chains allow for a different water content in the binding cavity of the protein, but the waters that remain in common occupy the same positions from one peptide to another [48]. This has prompted some authors to suggest that in a sense, water acts as part of the protein, moving around to change the shape of the binding cavity by selecting from a predefined set of conserved positions [49]. Intriguingly, we arrive at a similar conclusion after approaching the problem from a different angle. Our independent study, therefore, complements previous structural work, adding an energetic perspective, and corroborates this hypothesis.

Our results on bridging waters are in semiquantitative agreement with published structural studies on ligand binding to bovine pancreatic RNase, where a small number (usually < 10) of bridging water molecules is observed, thereby lending credence to the TIP3P water model, the ff14SB, and general Amber force field 2.11, and their compatibility. It is also instructive to compare published structural data with the behavior of the complexes in our molecular dynamics simulations. For example, crystal structures of the PUA–RNase [7] and PAX–RNase [50] complexes show that the two compounds form bridging interactions with the protein on both ends, through their uracyl and adenosyl, and thymine and adenosyl groups, respectively (Fig 7A and 7D).We observed that during dynamics, these water molecules rapidly exchange with others, while generally maintaining similar networks of interactions to those from the crystal structures (compare the starting structure in Fig 7A, 7B and 7C). Moreover, in certain trajectories, the ligands tended to reorient themselves to make more amine–phosphate interactions with the protein, which was accompanied with reorganization of the bridging interactions (compare 7D to 7E and F, also see Fig 3C and 3D). Interestingly, sodium ions can also become involved in the bridging network formed by the adenosyl group (Fig 7C). Additionally, in the course of dynamics, sodium ions tended to approach the phosphates and form long-lived bridging interactions among the phosphates and negatively charged or polar side chains from the protein, along with backbone oxygen atoms (see, for example, Fig 7B, 7C, 7E, 7F, 7H and 7I). One such example is the RNase–C3P complex [8], where the ligand drifted from its starting position in the center of the groove to its edge, forming a long-lived interaction with a sodium ion and E86, along with nearby backbone and side chain oxygen atoms (see Figs 7G, 7H, 7I and 3A and 3B). Trajectory analysis, together with the data presented in Fig 5, paints a picture of numerous (hundreds to thousands) bridging water molecules, which rapidly exchange with each other, and much less numerous (usually below 10) sodium ions, which tend to form longer-lived bridging interactions. Indeed, the lack of any bridging Cl- in our entire data set likely points to the importance of the phosphate groups for binding to RNase A. Moreover, the complete lack of sodium ions, bridging or otherwise, in the 22 published structures we have used to perform molecular dynamics, taken together with the presence of citrate or sulfate ions in some of the structures, and the difficulty of distinguishing Na+ from water in crystallography [51,52], suggests that the number of sodium ions in the structural databases may be underestimated.

Fig 7

Comparison between crystal and simulated structures.

(A) The RNase A–PUA complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres, polar contacts are shown as dashed lines. The locations of the urydyl and adenosyl moieties are indicated with a red and a blue arrow, respectively. (B) The RNase A–PUA complex after 42 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, Q10, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues and Q10 are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines. The locations of the urydyl and adenosyl moieties are indicated with a red and a blue arrow, respectively. (C) The RNase A–PUA complex after 100 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, Q10, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues and Q10 are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines. The locations of the urydyl and adenosyl moieties are indicated with a red and a blue arrow, respectively. (D) The RNase A–PAX complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres, polar contacts are shown as dashed lines. The locations of the thymine and adenosyl moieties are indicated with a red and a blue arrow, respectively. (E) The RNase A–PAX complex after 60 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines. The locations of the thymine and adenosyl moieties are indicated with a red and a blue arrow, respectively. (F) The RNase A–PAX complex after 100 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines. The locations of the thymine and adenosyl moieties are indicated with a red and a blue arrow, respectively. (G) The RNase A–C3P complex from a crystal structure. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres, polar contacts are shown as dashed lines. (H) The RNase A–C3P complex after 60 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines. (I) The RNase A–C3P complex after 100 ns of production dynamics. The protein backbone is in cartoon representation, colored in gray. Lysines, arginines, and the ligand are given in stick representation with carbon atoms in white, nitrogens in blue, oxygens in red, and phosphorus in orange. Certain lysine and arginine residues are labeled. Water oxygens are represented as red spheres; sodium ions are represented as magenta spheres; polar contacts are shown as dashed lines.

Our work demonstrates that apart from the long-lived, high residence time bridging interactions, there exist 2–3 orders of magnitude more abundant short-lived (< 1 ns) bridging interactions. We show that residence times follow a power law distribution–there exist a single-digit number of bridging interactions with residence times on the order of tens of nanoseconds, tens of bridging interactions with lifetimes between 1 and 10 ns, and hundreds to thousands of such interactions with residence times below 1 ns. In this respect, bridging solvent behaves similarly to hydrating solvent [53,54].

Crucially, a simple, semiquantitative consideration of bridging solvent shows that for many RNase ligands, the energetic contribution to binding of the short-lived bridges is comparable or greater than that of the long-lived bridges. A convenient example is provided by the RNase A–C3P system where in replica 2, all bridging interactions have lifetimes below 1 ns, whereas replica 4 has multiple long-lived bridging interactions– 3 above 10 ns and 18 between 1 and 10 ns. However, the MM-GBSA and MM-PBSA results for replica 2 are lower than for replica 4 (-54.7 and -40.7 kcal/mol versus -19.7 and -31.6 kcal/mol, respectively.) While this is a crude, semiquantitative treatment of the matter (as other differences in between replicas may also influence the computed results), it is useful in highlighting that long-lived bridges may not be the dominant energetic contributor. Indeed, long residence times do not indicate particularly strong protein–water or water–ligand interactions, but rather a topography that prevents the water molecule from exchanging by a cooperative mechanism [55], such as exchange-mediated orientational randomization (EMOR) at high confinement [56]. Therefore, preserving bridging interactions from crystallographic structures in a drug design campaign targeting a particular protein or protein–protein interaction [57] might not be sufficient. One also needs to account for short-lived solvent bridges. It is likely that molecular dynamics can be a powerful tool in this regard. This line of reasoning is further corroborated by the observation that a strong, positive correlation (R2 = 0.55, see S1 File) exists between ligand potency and the average number of bridging molecules, observed in our simulations.

The ΔH/pK plots excluding solvent show no significant correlation. Moreover, the lines of best fit exhibit greatly differing slopes, i.e. they cross each other at large angles. Conversely, the corresponding plots with bridging solvent tend to exhibit significant correlations and smaller angles in between the lines of best fit. These findings can be interpreted as suggesting that the calculations without solvent omit the key factors, determining RNase–ligand binding, producing random, largely orthogonal patterns. Conversely, the calculations with bridging solvent include in themselves most of the determinants governing binding, and consistently produce reliable, reproducible patterns.

The variations in ΔH and R2 across replicas indicate that the individual 100 ns simulations do not visit all conformations, relevant to binding. This is likely to be particularly pronounced for the pyrophosphate (or, equivalently, diphosphate) series of ligands, which also contact the protein’s unstructured regions, including the highly flexible N-terminal loop (see Figs 1, 2, and S2). Indeed, out of all lysine and arginine residues, K1 exhibits the greatest rotamer diversity from the 22 published structures (Fig 2). For loops, a great number of conformational states is possible, especially for terminal loops, which are bound only on one end and have vastly more degrees of freedom than loops which are bound on both ends. During the simulations, we observed conformations where K1 contacts the ligands and others where it is far removed from them. Moreover, we observed conformations where K7, R39, and K41 contact the ligands, as well as conformations where they do not (Fig 7). Finally, although more distant residues, such as K66 and R85 (see Fig 2) did not participate substantially in ligand binding in our simulations, we cannot rule out that they are involved in vivo, where sampling of conformational space is not limited to 400 ns. Indeed, this is a vast conformational space we cannot hope to sample exhaustively. It has previously been demonstrated that in such a scenario, performing several shorter replicas is far more efficient in terms of sampling the relevant conformational space than performing a lesser number of longer replicas [58]. This particular system offers another, highly illustrative example of why this is the case. Some of the smaller ligands, being more mobile, tended to leave the binding cleft near the end of the simulations in certain replicas. In order to observe multiple binding and unbinding events for these complexes, we would need to increase the length of our simulations by at least an order of magnitude, likely two or more, which is presently not feasible. Therefore, a more efficient strategy to obtain different conformations, relevant to protein–ligand binding, is to perform several shorter replicas, which is the approach we have adopted. Indeed, although the ΔH and R2 values vary in between replicas, R2 consistently points to a significant correlation between affinity and bridging interactions, indicating that this is the case throughout all of conformational space. In other words, the strong correlations and the reproducible patterns we observe in our simulations across a wide spectrum of affinity and in four independent replicas indicate that our results and conclusions are robust.

It is noteworthy that the greatest outliers in the plots in Fig 6C and 6D are (pyro)phosphate-containing ligands which during dynamics make extensive amine–phosphate contacts with the protein through its arginine and lysine residues. Indeed, RNase A is abundant in such residues, mostly within contact distance of the ligands around the binding groove (Fig 2). For example, in Fig 6C, C5P and A3P consistently lie well below the lines of best fit, meaning that their affinity for RNase A is greatly overestimated, eroding overall R2. As previously demonstrated, these ligands make little direct amine–phosphate contacts with the protein in the crystal structures, but reorient themselves to make extensive contacts of this nature during molecular dynamics (Fig 4). Moreover, examining the per-residue energies reveals that lysine and arginine residues are indeed the greatest contributors to binding (Fig 8).

Fig 8

Sources of affinity assessed via energetics analysis based on protein–ligand complex trajectories.

(A) The RNase A–C5P complex after 100 ns of molecular dynamics colored by per-residue energies of interaction calculated from the entire trajectory. Blue indicates a favorable contribution to binding, red indicates an unfavorable contribution. R39, K41, K66, R85, and D121 are shown in stick representation and explicitly labeled; the ligand is also shown in stick representation. Note the frame in the image is the same as in Fig 4B. (B) The RNase A–A3P complex after 100 ns of molecular dynamics colored by per-residue energies of interaction calculated from the entire trajectory. Blue indicates a favorable contribution to binding, red indicates an unfavorable contribution. K1, E2, K7, R10, K37, and K41 are shown in stick representation and explicitly labeled; the ligand is also shown in stick representation. Note the frame in the image is the same as in Fig 4D.

The only negative contribution to binding affinity is made by acidic residues in the vicinity of the phosphates, whereas the greatest positive contributions come from arginines and lysines, clustered together with the phosphates (compare Figs 4 and 8), highlighting once more the key role of the phosphates. We also note the importance of salt-linked triads to protein–ligand binding in this system, as this is a cooperative interaction that has previously been shown to be involved in protein folding [42,43] and protein–protein recognition and binding [44].

Our structural and energetics analysis suggest that the overestimated ΔH values stem from the overly attractive amine–phosphate potential–a previously recognized problem in the CHARMM and Amber series of force fields [59]. Indeed, given the abundance of amine-containing side chains on the protein side and the abundance of phosphates on the ligand side, RNase–ligand simulations are likely to be heavily burdened by the inaccuracies in the parameterization of the amine–phosphate interaction. We, therefore, applied NBFIX (nonbonded fix) corrections by increasing the σ parameter in the Lennard-Jones potential for the N–O = P (amine nitrogen–phosphate oxygen) pair [59,60]. NBFIX parameters are pair-specific corrections, i.e. they affect only a specific atom pair, e.g. amine nitrogen–carboxylate oxygen or amine nitrogen–phosphate oxygen. The NBFIX corrections represent an increase in the σ parameter for a given atom pair, making it less attractive. We also applied NBFIX corrections to the CT–CT atom pair (the interaction between two alyphatic carbons) [61], as these have also been shown to be overly stabilizing, artificially reducing the radius of gyration of simulated proteins in comparison to experiment [62], as well as NBFIX corrections to the Na+–Cl- (sodium ion–chloride ion) and Na+–O = P (sodium ion–phosphate oxygen) pair [63]. After performing analogous analysis on the systems simulated with NBFIX corrections, we found that compounds that lack phosphate groups (TXS, U1S, U2S, 0G0, 0FT, and 0EY) exhibit very similar ΔH values to their respective non-NBFIX behavior (compare S6A and S6C and S6B and S6D Fig, note that the Y axis scales are different; also see S1 File), confirming that the amine–phosphate interactions are the dominant difference in between the NBFIX and standard parameter simulations. For the phosphate-containing compounds, the modified potential resulted in a reduced propensity for direct amine–phopshate contacts during production dynamics. For the pyrophosphate derivatives, this also resulted in a decrease in the number of salt-linked triads and clusters formed by phosphate oxygens and positively charged side chains. This corresponded to an upshift in ΔH values for the phosphate-containing compounds, as compared to the non-NBFIX simulations. The NBFIX corrections, however, reduced R2, rather than increase it. This implies that the previously published parameters [59,60] are unsuitable for the current protein–ligand data set, highlighting the problem of parameter transferability in molecular dynamics. Parameterizing the vast chemical space that drug and drug-like molecules constitute is certainly a daunting task. A recent analysis on the charging and dispersive-repulsive contributions to solvation free energies in an analysis of GAFF and OPLSA-AA with the TIP3P, SPCE, and OPC3 water models has demonstrated that while GAFF generally performs well, further improvements are more likely to be obtained through “adjusting and tuning the available atomic charge calculation protocols, namely AM1/BCC for GAFF2 and 1.14*CM1A or 1.14*CM1A-LBCC for OPLS-AA” [64]. Our extensive, fully independent study complements this work by examining a different property in a different system, leading us to concur with Vasetti et al. [64].

Finally, we discuss two methodologically important points. First, we have opted to use a relatively large solvation shell (15 Å instead of the usual 10–11–12) because we noticed an imaging artifact in certain trajectories when using smaller solvation shells where bridging interactions could not be detected by cpptraj, leading to an underestimation of the magnitude of ΔH. Second, we point out that the computationally much cheaper MM-GBSA calculations offer similar performance to the more theoretically rigorous MM-PBSA approach in terms of R2 using the standard Amber parameters (Fig 6), despite being a much cruder approximation, at a lower level of theory. Thus, in scenarios with limited computational resources, particularly limited CPU resources, MM-GBSA offers a viable alternative to the more demanding (typically, 1–2 orders of magnitude in terms of compute time) MM-PBSA analysis.

Conclusions

We have performed and presented extensive free energy calculations on a protein–ligand data set spanning 22 compounds and nearly five orders of magnitude in affinity. This work presents an independent test of the state-of-the-art fixed-charge force fields. Crucially, we describe an automated workflow for the detection of bridging interactions in protein–ligand binding–an often important but neglected factor in intermolecular interactions. Moreover, our workflow is extensible and amenable to modification to accommodate other types of biologically important (macro)molecules such as nucleic acids, saccharides, lipids, and polyamines, as well as multicomponent systems of these, where complex, multifactor bridging interactions are likely to be crucial for an accurate, atomic-level description of the biology of interest [65,66]. The workflow we describe is also likely to be applicable to four-dimensional, time-dependent quantitative-structure activity relationship (4D-QSAR) studies [45] in drug design and optimization.

Cα RMSDs.

(A) Cα RMSDs for the protein over the course of production dynamics in replica 1. (B) Cα RMSDs for the protein over the course of production dynamics in replica 2. (C) Cα RMSDs for the protein over the course of production dynamics in replica 3. (D) Cα RMSDs for the protein over the course of production dynamics in replica 4.

(EPS)

Click here for additional data file.

Cα RMSFs.

(A) Cα RMSFs for the protein over the course of production dynamics in replica 1. (B) Cα RMSFs for the protein over the course of production dynamics in replica 2. (C) CαRMSFs for the protein over the course of production dynamics in replica 3. (D) Cα RMSFs for the protein over the course of production dynamics in replica 4.

(EPS)

Click here for additional data file.

Ligand heavy atom RMSDs.

(A) Ligand heavy atom RMSDs over the course of production dynamics in replica 1. (B) Ligand heavy atom RMSDs over the course of production dynamics in replica 2. (C) Ligand heavy atom RMSDs over the course of production dynamics in replica 3. (D)Ligand heavy atom RMSDs over the course of production dynamics in replica 4.

(EPS)

Click here for additional data file.

Enthalpy convergence without bridging solvent.

(A) Enthalpy (ΔH) convergence from the MM-PBSA calculations from replica 1 (top), replica 2, replica 3, and replica 4 (bottom); bridging waters and ions are not included in the calculations. (B) Enthalpy (ΔH) convergence from the MM-GBSA calculations from replica 1 (top), replica 2, replica 3, and replica 4 (bottom); bridging waters and ions are not included in the calculations.

(PDF)

Click here for additional data file.

Enthalpy convergence with bridging solvent.

(A) Enthalpy (ΔH) convergence from the MM-PBSA calculations from replica 1 (top), replica 2, replica 3, and replica 4 (bottom); bridging waters and ions are included in the calculations. (B) Enthalpy (ΔH) convergence from the MM-GBSA calculations from replica 1 (top), replica 2, replica 3, and replica 4 (bottom); bridging waters and ions are included in the calculations.

(EPS)

Click here for additional data file.

ΔH/pK correlations with NBFIX parameters.

(A) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-PBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black) with NBFIX corrections; bridging waters and ions are included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color. (B) Correlations between computed enthalpy (ΔH) and experimental affinity from the MM-GBSA calculations from replica 1 (blue), replica 2 (red), replica 3 (green), and replica 4 (black) with NBFIX corrections; bridging waters and ions are included in the calculations. For each replica, the line of best fit and the corresponding equation are also given in the respective color.

(EPS)

Click here for additional data file.

Supplementary data.

Sheet 1 contains the experimental affinity data between RNase A and the 22 ligands, expressed as pKi values, and the computed enthalpies from the four replicas. Sheet 2 contains the average Cα RMSD values; sheet 3 –the Cα RMSFs; sheet 4 –the ligand heavy atom RMSDs; sheet 5 –the convergence of the MM-PBSA and MM-GBSA values without bridging solvent; sheet 6 –the convergence of the MM-PBSA and MM-GBSA values with bridging solvent; sheet 7 –the average number of bridging waters and ions across the four replicas for every compound and the corresponding standard deviation, followed by the average number of bridges and the standard deviations; sheet 8 –the radial distribution functions for bridging water and bridging sodium ions around the ligands calculated from the four replicas.

(XLSX)

Click here for additional data file.

11 Sep 2019

[EXSCINDED]

PONE-D-19-23148

Bridging solvent molecules mediate RNase A – ligand binding

PLOS ONE

Dear Dr Ivanov,

Thank you for submitting your manuscript to PLOS ONE. After careful consideration, we feel that it has merit but does not fully meet PLOS ONE’s publication criteria as it currently stands. Therefore, we invite you to submit a revised version of the manuscript that addresses the points raised during the review process.

==============================

The two reviewers combine for 16 important points. While some of them are minor, most of them are significant enough to make me unable to determine the validity of this work.

I encourage you to address all of the points so that we can fully assess this manuscript.

==============================

We would appreciate receiving your revised manuscript by Oct 26 2019 11:59PM. When you are ready to submit your revision, log on to https://www.editorialmanager.com/pone/ and select the 'Submissions Needing Revision' folder to locate your manuscript file.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter.

To enhance the reproducibility of your results, we recommend that if applicable you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. For instructions see: http://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols

Please include the following items when submitting your revised manuscript:

A rebuttal letter that responds to each point raised by the academic editor and reviewer(s). This letter should be uploaded as separate file and labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. This file should be uploaded as separate file and labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. This file should be uploaded as separate file and labeled 'Manuscript'.

Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

We look forward to receiving your revised manuscript.

Kind regards,

Freddie Salsbury , Jr, PhD

Academic Editor

PLOS ONE

Journal Requirements:

1.  When submitting your revision, we need you to address these additional requirements.

Please ensure that your manuscript meets PLOS ONE's style requirements, including those for file naming. The PLOS ONE style templates can be found at

http://www.journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and http://www.journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

[Note: HTML markup is below. Please do not edit.]

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Partly

Reviewer #2: Partly

**********

2. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: No

Reviewer #2: Yes

**********

3. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: No

Reviewer #2: Yes

**********

4. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors presented a MD study for RNase A with different ligand

compounds using the AMBER force field. Their key result includes

bridging water analysis, heavy atoms RMSD, and energetic analysis for

the ligands with and without water and ions. While the outcome that

inclusion of bridging waters lead to higher correlations with

experimental data is useful there are signficant problems with the

manuscript. Comments follow.

1) Importanly, the authors need to discuss published experimental

structural studies in relation to their work validation beyond two

lines on page 15. Based on text on page 2 (references 8 to 11) there

is experimental information available.

2) In Figure 1 the phosphates of the ligands are shown to be fully

protonated (ie. neutral). If this is the protonation state used in

the simualtions, then the results are rather meaningless as this state

will not exist in the vicinity of pH 7.0.

3) What is the need to use the modified parameters for additional 8.8 or

10 microseconds? What are the modified parameters and to they impact

the results at all. This is not discussed in the methods, but appears

to be associated with the use of NBFIX parameters that is presented

later in the manuscript. This issue needs significant clarification.

4) The MM-PBSA/GBSA as presented as enthalpies of binding; however, PB

and GB are giving the free energies of solvation, such that using free

energy of binding may be more appropraite. This needs to be

clarified.

5) The authors need to perform some structural analysis (maybe

spatial/angular/radial distribution) where the bridging functional

groups between protein and the ligand are involved. Or role of the

phosphate group in the binding mechanism of salt bridges

(phosphate-lysine etc). The detailed structural analysis will show the

importance of those functional groups that get involved in the most

likely bridging through the entire trajectory.

6) On page 4 the authors state "Average Cα RMSD is 0.07 Å." This is

clearly incorrect based on the RMSD figures

7) The RMSD/RMSF figures should be moved to the Supplementary Material

8) Figures 7 and 8 may also be moved to the Supplementary Material.

It would be nice to provide some interpretation in relation to Figure

6. Why a certain fragment has less or more number of bridges in

comparison to another ligand? Also, some more details/discussion in

the calculation to the lifetimes of the ligand or all the all bound

for the duration of the simulations. Please clarify.

9) Clarify the scientific information derived from Figures 9? We know the

importance of water and ions in the formation of the bridges. How

about having a somewhat larger variation in the DH for a different

replica in particularly for PAX, PUA, PAP (Figure 9C and D)?

10) Indeed, the Figures 9 and 10 should be moved to the Supplementary

Material and the results included in the main text as a table. This

will make it much easier for the reader to access those results. This

is important as the main result is the improved correlation with

experiment when water is included in the enthalpy of binding

calculation.

11) The following needs some structural analysis to support the statement:

“the plots in Fig 9C and D, and the molecular dynamics trajectories

reveals that the greatest outliers are phosphate or pyrophosphate

containing molecules whose interactions with the protein are dominated

by phosphate – lysine or phosphate – arginine interactions.”

Reviewer #2: The manuscript PONE-D-19-23148 by Stefan Ivanov et al. studies the water bridging interactions in RNase-ligand binding using explicit molecular dynamics simulations. Each of the 22 compounds bound to RNase A is simulated in four independent replicas with a length of 100 ns for each replica. Some analyses are performed on bridging water molecules to show the importance of their role in RNase-ligand binding. MM-PBSA, MM-GBSA are used to compute the binding energies of the compounds. In order to explain the overestimated dH values, the same systems are simulated with nonbonded fix (NBFIX) corrections to the force fields. And similar analyses are performed on the second round of simulations. However, there are a number of issues that I hope the authors could address.

Major points are as follows:

1. One advantage of explicit solvent MD simulation is to provide structural information at atomic resolution. However, I couldn't see any detailed structural information of the bridging water molecules in the binding site. For example, where are the water molecules in the binding sites? Is there any hydration site in the binding pocket that is critical for the binding? What are the configurations of the residues, ligands and water molecules in the binding site? Overall, I think the authors should include more detailed discussions about bridging interactions with atomic structures.

2. In the paragraph starting from line 186, "As we were primarily interested in relative binding energies ...", entropy calculations are omitted, however, I don’t think the assumption that the entropy changes are similar across different systems in this work is valid without any proof. Some of the 22 compounds simulated in this work (Fig 1) are quite different from each other in terms of rotational bonds, molecular weights. The water molecule involved in the bridging interactions may also have a considerable contribution to the entropy. Could you show any validation that the entropy is negligible in computing the relative free energies?

3. Although the enthalpy changes computed in the manuscript are related to the relative binding energies, they are not the same. Could you explain more about the choice of computing dHs instead of computing the ddGs?

4. I am puzzled by the plots in Fig. 9 and Fig. 10, which is trying to correlate between the computed relative binding energies and the absolute experimental binding affinities. Shouldn’t the correlation be between the computed ddG vs experimental ddG?

5. I also have concerns about the statistical stability of computed binding free energies. There are large variations between the dHs computed from four different replicas and the slops of fitted linear functions in Fig. 9 and Fig. 10 also show large variations, which may indicate that the conformational samplings from the simulations aren’t enough for the binding energy calculations.

Some minor points:

1. Line 179, “MM-PBA” should be “MM-PBSA”

2. In Fig. 1, the pKi value of U2S is in the wrong place.

3. Reference 56 is in the wrong line.

**********

6. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

[NOTE: If reviewer comments were submitted as an attachment file, they will be attached to this email and accessible via the submission site. Please log into your account, locate the manuscript record, and check for the action link "View Attachments". If this link does not appear, there are no attachment files to be viewed.]

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, https://pacev2.apexcovantage.com/. PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Registration is free. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at figures@plos.org. Please note that Supporting Information files do not need this step.

24 Sep 2019

We thank the reviewers for their positive and constructive comments. Below, we provide a point-by-point response to each of their concerns.

Reviewer #1: The authors presented a MD study for RNase A with different ligand compounds using the AMBER force field. Their key result includes bridging water analysis, heavy atoms RMSD, and energetic analysis for the ligands with and without water and ions. While the outcome that inclusion of bridging waters lead to higher correlations with experimental data is useful there are signficant problems with the manuscript. Comments follow.

1) Importanly, the authors need to discuss published experimental structural studies in relation to their work validation beyond two lines on page 15. Based on text on page 2 (references 8 to 11) there is experimental information available.

We thank the reviewer for his comment, which we feel is very constructive. Moreover, it aligns nicely with several other comments by both reviewers. Accordingly, as requested by the reviewer, we have included structural analysis on the behavior of several complexes throughout production dynamics. We believe the added details, along with further structural and analysis details we have included throughout the manuscript, will be highly instructive for the PLOS ONE readership. More details are provided below.

2) In Figure 1 the phosphates of the ligands are shown to be fully protonated (ie. neutral). If this is the protonation state used in the simualtions, then the results are rather meaningless as this state will not exist in the vicinity of pH 7.0.

The reviewer is indeed correct to point out that the phosphate groups will not be protonated at pH 7. Accordingly, in our original submission, the last line of the Figure 1 caption explicitly states that the phosphate groups have been drawn protonated, but have been simulated in a fully deprotonated state - “Note that phosphate groups are drawn protonated, whereas they have been simulated in a fully deprotonated state. ” We have deliberately placed this sentence at the very end of the caption to make it stand out. In our revised manuscript, the Figure 1 caption remains unchanged. We thank the reviewer for pointing out that this information can be hard to notice. Therefore, to help guide the reader, we have added this notice in the main text as well.

3) What is the need to use the modified parameters for additional 8.8 or 10 microseconds? What are the modified parameters and to they impact the results at all. This is not discussed in the methods, but appears to be associated with the use of NBFIX parameters that is presented later in the manuscript. This issue needs significant clarification.

The NBFIX parameters are pair-specific corrections to the potential energy function. It has previously been noted by other authors that certain interactions in the CHARMM and AMBER force fields are significantly overestimated. For example, if glycine is simulated in water, it will aggregate in simulations, whereas in reality it is highly soluble. This is because the amine – acetate interaction between the N-terminal amine and the C-terminal COO- group is too attractive in the current force fields. This is corrected by increasing the σ (sigma) parameter in the potential energy function for the amine nitrogen (N) – carboxylate oxygen (O) atom pair. More details about NBFIX parameters and parametrization can be found in the papers we reference. We have included simulations with NBFIX corrections because testing how general and useful they are is of great interest to the drug design and computational chemistry community. We have 22 complexes, for each of which we have performed 400 ns (4 replicas x 100 ns) of production dynamics = 8.8 microseconds in total. For a fair comparison, we have performed the exact same amount of sampling with NBFIX corrections, otherwise any comparison might be misleading. Whether or not the NBFIX parameters are applicable to our data set, however, is of secondary importance in this work; the main focus is on bridging interactions. Therefore, we briefly discuss NBFIX in the Discussion section. We have chosen not to include the NBFIX simulations in the Methods section, but to simply state that we have performed the same type of analysis on simulations performed with NBFIX corrections, in order to maintain the focus of this work on bridging interactions. We thank the reviewer for pointing out this can be seen as a little unclear. We have added more details and explanations on the NBFIX work to make this section clearer and easier to digest.

4) The MM-PBSA/GBSA as presented as enthalpies of binding; however, PB and GB are giving the free energies of solvation, such that using free energy of binding may be more appropraite. This needs to be clarified.

Polar solvation free energies are obtained by solving the linearized Poisson-Boltzmann equation or the Generalized Born equation. This is where the “PB” and “GB” in MM-PBSA and MM-GBSA come from. However, Poisson-Boltzmann and Generalized Born calculations are not the only calculations in MM-PBSA and MM-GBSA – they also include surface area (or molecular volume) calculations for nonpolar solvation, gas phase energies, and van der Waals terms, calculated in accord with the force field parameters. Combined, these calculations are known as molecular mechanics – Poisson-Boltzmann surface area (MM-PBSA) and molecular mechanics – Generalized Born surface area (MM-GBSA) calculations. MM-PBSA and MM-GBSA, by themselves, produce the enthalpy of interaction (DH). In order to obtain the free energy of interaction (DG), as opposed to the enthalpy, one also needs to calculate the entropy of interaction (DS). This, however, is not possible at present, perhaps except for the very simplest of model systems, which is certainly not the case with our data set. Although calculations such as normal mode analysis (NMA) are sometimes performed, they are too expensive and too unreliable to obtain accurate DGs. Therefore, we focus our attention on the only part of the binding free energy that can be calculated with a certain amount of reliability – enthalpy (DH). Entropy calculations await further methodological and computational developments. More details on MM-PB(GB)SA can be found in the papers we reference.

5) The authors need to perform some structural analysis (maybe spatial/angular/radial distribution) where the bridging functional groups between protein and the ligand are involved. Or role of the phosphate group in the binding mechanism of salt bridges (phosphate-lysine etc). The detailed structural analysis will show the importance of those functional groups that get involved in the most likely bridging through the entire trajectory.

We thank the reviewer for this comment, as we feel it is particularly constructive and has led to interesting insights. We have performed additional analysis on the bridging waters and ions, examining their radial distribution around the ligands. We find that short-lived bridging is dominated by water, whereas long-lived bridging is dominated by sodium ions, which interact mostly with the ligands’ phosphate groups and oxygen atoms from protein side chains, as well as backbone oxygens. Moreover, the complete lack of bridging chloride ions (Cl-) in our entire data set clearly points to the preeminence of the phosphates in RNase – ligand binding. Moreover, even compounds that lack phosphate groups do not employ bridging Cl-, only Na+ and water, although, typically, to a lesser extent than phosphate-containing compounds. We have included more structural detail and discussion on this to the Results and Discussion sections.

6) On page 4 the authors state "Average Cα RMSD is 0.07 Å." This is clearly incorrect based on the RMSD figures

We thank the reviewer for pointing out that this figure and the caption information can be somewhat confusing. The figure presents the starting crystal and NMR structures we have used to perform MD. The 0.07 Å number for the Cα RMSD is given for the RMSD between the starting crystal and NMR structures, not the average RMSD during dynamics. The 0.07 number can easily be verified by loading chain A from PDB structures 4g90, 1rnm, 4g8y, 3d6p, 1w4p, 1jvu, 1qhc, 1jn4, 1w4o, 1afl, 1o0f, 1afk, 3d6o, 1u1b, 1o0h, 1o0n, 4g8v, 1o0m, 3d8z, 1z6s, 1rpf, 1w4q in Pymol and aligning everything to 4g90 – the output is 0.073 Å for average Cα RMSD. In order to avoid confusion, we have removed the number from the caption and have modified the figure title.

7) The RMSD/RMSF figures should be moved to the Supplementary Material

As suggested by the reviewer, these and several other figures have been moved to SI. We thank the reviewer for this suggestion, as this has opened up more space in the manuscript for more figures and discussion on structural analysis, which has improved the quality of the manuscript.

8) Figures 7 and 8 may also be moved to the Supplementary Material. It would be nice to provide some interpretation in relation to Figure 6. Why a certain fragment has less or more number of bridges in comparison to another ligand? Also, some more details/discussion in the calculation to the lifetimes of the ligand or all the all bound for the duration of the simulations. Please clarify.

As suggested by the reviewer, we have moved these figures to SI. In order to clarify what we mean by “bridging molecule” and “bridging interaction,” we have added more illustrative examples in the relevant Methods section. In the Results section, we note that phosphate-containing ligands tend to attract more bridging waters and ions than compounds that do not contain phosphates, and that the large, pyrophosphate-contining compounds attract more bridging waters than the small compounds that contain only one phosphate group. We note that a pyrophosphate (or diphosphate) group is simply two phosphate groups condensed together. It becomes clear, then, that phosphates tend to attract more bridging waters and ions. The more such waters and ions become involved, the larger the complexes become, and the longer it takes for the relevant MM-PBSA and MM-GBSA calculations to converge. This is clarified in the Results section.

9) Clarify the scientific information derived from Figures 9? We know the importance of water and ions in the formation of the bridges. How about having a somewhat larger variation in the DH for a different replica in particularly for PAX, PUA, PAP (Figure 9C and D)?

This comment aligns nicely with comment #5 from Reviewer 2. The variations in DH and R2 indicate that the individual simulations are not long enough to sample all conformations relevant to binding. Indeed, as we point out in the manuscript, some of the ligands contact the protein’s unstructured regions, which also bear lysine and arginine residues, especially the N-terminal loop, which bears lysine 1. Because this loop is free on the N terminus, it has a vast number of degrees of freedom, far more than the other loops, which are bound on both ends. This is reflected in the RMSF values of the loops. In fact, K1 has the greatest root-mean-square fluctuations of all residues, as evident from the relevant figure. In some conformations, K1 contacts the ligands, in others it is far removed from them. The same applies to many of the other lysines and arginines in the structure, e.g. K7, R39, K41, etc. We can not sample all possible conformations. In such a scenario, one should run multiple replicas to sample as much of conformation space as possible. Although there are variations in between replicas, the correlations are highly significant in every replica, which indicates that bridging waters are indeed a key determinant. Extended discussion of this has been included in the Discussion section.

10) Indeed, the Figures 9 and 10 should be moved to the Supplementary Material and the results included in the main text as a table. This will make it much easier for the reader to access those results. This is important as the main result is the improved correlation with experiment when water is included in the enthalpy of binding calculation.

As suggested by the reviewer, Figure 9 has been moved to SI. However, we feel that our data set is far too large to be included in the body of the manuscript in table form. We fully agree with the reviewer that the reader should have access to this data to examine it in detail. Therefore, we have included this data, along with more of our results, in a supplementary file the reader can download and examine in great detail.

11) The following needs some structural analysis to support the statement:

“the plots in Fig 9C and D, and the molecular dynamics trajectories reveals that the greatest outliers are phosphate or pyrophosphate containing molecules whose interactions with the protein are dominated by phosphate – lysine or phosphate – arginine interactions.”

As requested by the reviewer, we have performed ample structural and energetics analysis supporting our conclusions and have placed additional information on the relevant methods and discussion in the appropriate sections of the manuscript. Briefly, in our MM-PBSA and MM-GBSA calculations, similarly to our previously published work, we have enabled per-residue energy decompositions. This allows one to assess the contribution of each individual residue to binding during the molecular dynamics simulations. In other words, this analysis reveals whether a given residue makes a favorable or unfavorable contribution to binding or is largely indifferent. Together with our structural analysis, this confirmed that the overestimation in DH likely comes mostly from the inaccuracies in the parametrization of the amine – phosphate interaction. The relevant discussion has been placed in the Discussion section of the manuscript.

Reviewer #2: The manuscript PONE-D-19-23148 by Stefan Ivanov et al. studies the water bridging interactions in RNase-ligand binding using explicit molecular dynamics simulations. Each of the 22 compounds bound to RNase A is simulated in four independent replicas with a length of 100 ns for each replica. Some analyses are performed on bridging water molecules to show the importance of their role in RNase-ligand binding. MM-PBSA, MM-GBSA are used to compute the binding energies of the compounds. In order to explain the overestimated dH values, the same systems are simulated with nonbonded fix (NBFIX) corrections to the force fields. And similar analyses are performed on the second round of simulations. However, there are a number of issues that I hope the authors could address.

Major points are as follows:

1. One advantage of explicit solvent MD simulation is to provide structural information at atomic resolution. However, I couldn't see any detailed structural information of the bridging water molecules in the binding site. For example, where are the water molecules in the binding sites? Is there any hydration site in the binding pocket that is critical for the binding? What are the configurations of the residues, ligands and water molecules in the binding site? Overall, I think the authors should include more detailed discussions about bridging interactions with atomic structures.

This comment nicely aligns with several of the comments by Reviewer 1. We have included ample and illustrative structural information and analysis throughout the manuscript see above), which we believe would be highly instructive and very useful to the PLOS ONE readership.

2. In the paragraph starting from line 186, "As we were primarily interested in relative binding energies ...", entropy calculations are omitted, however, I don’t think the assumption that the entropy changes are similar across different systems in this work is valid without any proof. Some of the 22 compounds simulated in this work (Fig 1) are quite different from each other in terms of rotational bonds, molecular weights. The water molecule involved in the bridging interactions may also have a considerable contribution to the entropy. Could you show any validation that the entropy is negligible in computing the relative free energies?

First, we note that we are not performing relative free energy calculations (DDG) but enthalpy calculations (DH). The reviewer is indeed correct to point out that the ligands are significantly dissimilar and entropy cancelation should not be expected. Indeed, we chose to omit entropy calculations not because we believe that these terms will cancel out, but because at the given time it is impossible to calculate entropy reliably. It has been shown in previous work by us and others that enthalpies can be calculated reliably, but entropies can not. Given that entropy calculations are highly expensive but likely to introduce more noise than signal, we chose to not perform them. Moreover, given that we can demonstrate an important role for bridging waters – they key focus of this work – without entropy calculations, we feel that attempting to calculate entropy would be redundant. The phrase “As we were primarily interested in relative binding energies” pertains to the R2 values between the computed enthalpies, not the exact values of DH. Our 2016 paper shows that computed DH values are of similar magnitude to experimental DH results from ITC. This is explicitly referenced and discussed, as we feel that this is highly instructive for the reader.

3. Although the enthalpy changes computed in the manuscript are related to the relative binding energies, they are not the same. Could you explain more about the choice of computing dHs instead of computing the ddGs?

This relates to the reviewer’s previous comment. Given that with MM-PB(GB)SA we can not reliably compute entropies, we are left with the enthalpies of interaction, i.e. it is impossible to reliably calculate DG with this technique. Although other techniques, such as thermodynamic integration (TI), circumvent the problem of explicit entropy calculation, they suffer from other methodological weaknesses, such as dependence of the computed DG on box size, ligand net charge (calculated DGs for charged ligands are highly inaccurate), etc. This is a significant limitation in computational chemistry that is widely discussed in the literature, partly in our 2019 paper (referenced in this manuscript) as well, and the references therein. Stated briefly – it is currently not possible to reliably calculate DG for large and complex ligands, especially charged ones – this is a well known shortcoming of MD and computational chemistry in general. Given that it is impossible to calculate DG reliably but that it is possible to calculate DH with reasonable accuracy, we focus our analysis on enthalpies (DHs) rather than binding free energies (DGs).

4. I am puzzled by the plots in Fig. 9 and Fig. 10, which is trying to correlate between the computed relative binding energies and the absolute experimental binding affinities. Shouldn’t the correlation be between the computed ddG vs experimental ddG?

Figures 9 and 10 correlate the enthalpy of binding (DH) between RNase A and the 22 ligands with experimental affinities. We calculate only DH because it is practically not possible to calculate DG for these compounds. We then relate computed DH values to experimental affinity data from the PDBbind database, which comes in the form of pKi values. While experimental DG can be obtained from the pKi values, there is no need to do so because we are not attempting to calculate DG. We stress that we are not attempting to correlate calculated with experimental DG because we have not calculated DG.

5. I also have concerns about the statistical stability of computed binding free energies. There are large variations between the dHs computed from four different replicas and the slops of fitted linear functions in Fig. 9 and Fig. 10 also show large variations, which may indicate that the conformational samplings from the simulations aren’t enough for the binding energy calculations.

The reviewer is quite correct to point out that the variations in DH suggest that our simulations are not long enough to sample all relevant conformations – this is stated in our original submission as well. While performing longer simulations is certainly an option, it has been convincingly shown by others that more replicas offer much more efficient sampling than an equivalent amount of longer simulations. In our particular case, some of the smaller and more mobile ligands leave the binding cleft near the end of the simulations in some of the replicas. Doubling or even tripling the duration of the simulations will likely not be enough to observe multiple binding and unbinding events, which would require at least an order of magnitude more simulation time (likely several microseconds, at least). Therefore, a much more efficient strategy in this case is to perform several shorter simulations, rather that one long one, which is precisely the approach we have adopted. We have included more discussion on this in the manuscript, as well as a highly instructive reference the interested reader can examine.

Some minor points:

1. Line 179, “MM-PBA” should be “MM-PBSA”

2. In Fig. 1, the pKi value of U2S is in the wrong place.

3. Reference 56 is in the wrong line.

We thank the reviewer for his careful examination of our work and for noticing these errors. These have now been corrected.

Submitted filename: point_by_point_response.docx

Click here for additional data file.

10 Oct 2019

Bridging solvent molecules mediate RNase A – ligand binding

PONE-D-19-23148R1

Dear Dr. Ivanov,

We are pleased to inform you that your manuscript has been judged scientifically suitable for publication and will be formally accepted for publication once it complies with all outstanding technical requirements.

Within one week, you will receive an e-mail containing information on the amendments required prior to publication. When all required modifications have been addressed, you will receive a formal acceptance letter and your manuscript will proceed to our production department and be scheduled for publication.

Shortly after the formal acceptance letter is sent, an invoice for payment will follow. To ensure an efficient production and billing process, please log into Editorial Manager at https://www.editorialmanager.com/pone/, click the "Update My Information" link at the top of the page, and update your user information. If you have any billing related questions, please contact our Author Billing department directly at authorbilling@plos.org.

If your institution or institutions have a press office, please notify them about your upcoming paper to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, you must inform our press team as soon as possible and no later than 48 hours after receiving the formal acceptance. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information, please contact onepress@plos.org.

With kind regards,

Freddie Salsbury , Jr, PhD

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Reviewers' comments:

Reviewer's Responses to Questions

Comments to the Author

1. If the authors have adequately addressed your comments raised in a previous round of review and you feel that this manuscript is now acceptable for publication, you may indicate that here to bypass the “Comments to the Author” section, enter your conflict of interest statement in the “Confidential to Editor” section, and submit your "Accept" recommendation.

Reviewer #1: All comments have been addressed

Reviewer #2: All comments have been addressed

**********

2. Is the manuscript technically sound, and do the data support the conclusions?

The manuscript must describe a technically sound piece of scientific research with data that supports the conclusions. Experiments must have been conducted rigorously, with appropriate controls, replication, and sample sizes. The conclusions must be drawn appropriately based on the data presented.

Reviewer #1: Yes

Reviewer #2: Yes

**********

3. Has the statistical analysis been performed appropriately and rigorously?

Reviewer #1: Yes

Reviewer #2: Yes

**********

4. Have the authors made all data underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data—e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

**********

5. Is the manuscript presented in an intelligible fashion and written in standard English?

PLOS ONE does not copyedit accepted manuscripts, so the language in submitted articles must be clear, correct, and unambiguous. Any typographical or grammatical errors should be corrected at revision, so please note any specific errors here.

Reviewer #1: Yes

Reviewer #2: Yes

**********

6. Review Comments to the Author

Please use the space provided to explain your answers to the questions above. You may also include additional comments for the author, including concerns about dual publication, research ethics, or publication ethics. (Please upload your review as an attachment if it exceeds 20,000 characters)

Reviewer #1: The authors have addressed my concerns. xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx

Reviewer #2: I appreciate the efforts from authors to address the critiques from the previous review. The revised manuscript is improved and contains additional information that helps to understand the overall approach. I have no more major objections to its publication PLOS ONE. There are two minor points I want to point out here.

1. I think “relative binding energies” in Line 203 (revised manuscript) is very confusing. Considering change it to some more appropriate term like “binding enthalpies”

2. The subplots in Figure 6 look very crowded, especially the legends are overlapping with each other. I think a better way of showing the information is fitting a single line between the averaged dH to the experimental values. And move Fig 6 to SI.

**********

7. PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

14 Oct 2019

PONE-D-19-23148R1

Bridging solvent molecules mediate RNase A – ligand binding

Dear Dr. Ivanov:

I am pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now with our production department.

If your institution or institutions have a press office, please notify them about your upcoming paper at this point, to enable them to help maximize its impact. If they will be preparing press materials for this manuscript, please inform our press team within the next 48 hours. Your manuscript will remain under strict press embargo until 2 pm Eastern Time on the date of publication. For more information please contact onepress@plos.org.

For any other questions or concerns, please email plosone@plos.org.

Thank you for submitting your work to PLOS ONE.

With kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Freddie Salsbury , Jr

Academic Editor

PLOS ONE