
==== Front
PLoS One
PLoS One
plos
PLOS ONE
1932-6203
Public Library of Science San Francisco, CA USA

10.1371/journal.pone.0309996
PONE-D-24-21031
Research Article
Physical Sciences
Physics
Thermodynamics
Free Energy
Engineering and Technology
Mechanical Engineering
Rotors
Research and Analysis Methods
Simulation and Modeling
Physical Sciences
Chemistry
Computational Chemistry
Molecular Dynamics
Physical Sciences
Physics
Thermodynamics
Computer and Information Sciences
Artificial Intelligence
Machine Learning
Support Vector Machines
Medicine and Health Sciences
Pharmacology
Drug Research and Development
Drug Discovery
Physical Sciences
Materials Science
Material Properties
Density
Physical Sciences
Materials Science
Materials Physics
Density
Physical Sciences
Physics
Materials Physics
Density
Calculated hydration free energies become less accurate with increases in molecular weight
Calculated hydration free energies become less accurate with increases in molecular weight
https://orcid.org/0000-0002-6949-2328
Ivanov Stefan M. Conceptualization Data curation Formal analysis Visualization Writing – original draft Writing – review & editing *
Faculty of Pharmacy, Medical University of Sofia, Sofia, Bulgaria
Bhakat Soumendranath Editor
AlloTec Bio, UNITED STATES OF AMERICA
Competing Interests: The author has declared that no competing interests exist.

* E-mail: sivanov@ddg-pharmfac.net
19 9 2024
2024
19 9 e030999624 5 2024
22 8 2024
© 2024 Stefan M. Ivanov
2024
Stefan M. Ivanov
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

In order for computer-aided drug design to fulfil its long held promise of delivering new medicines faster and cheaper, extensive development and validation work must be done first. This pertains particularly to molecular dynamics force fields where one important aspect–the hydration free energy (HFE) of small molecules–is often insufficiently analyzed. While most benchmarking studies report excellent accuracies of calculated hydration free energies–usually within 2 kcal/mol of experimental values–we find that deeper analysis reveals significant shortcomings. Herein, we report a dependence of HFE prediction errors on ligand molecular weight–the higher the weight, the bigger the prediction error and the higher the probability the calculated result is erroneous by a large amount. We show that in the drug-like molecular weight region, HFE predictions can easily be off by 5 kcal/mol or more. This is likely to be highly problematic in a drug discovery and development setting. We make our HFE results and molecular descriptors freely and fully available in order to encourage deeper analysis of future molecular dynamics results and facilitate development of the next generation of force fields.

This study is financed by the European Union-NextGenerationEU through the National Recovery and Resilience Plan of the Republic of Bulgaria, project № BG-RRP-2.004-0004-C01. The in silico calculations were performed at the Centre of Excellence for Informatics and ICT, supported by the Science and Education for Smart Growth Operational Program and co-financed by the European Union through the European Structural and Investment Funds (Grant No. BG05M2OP001-1.001-0003). The funders played no role in the study design, data collection and analysis, decision to publish, or preparation of the manuscript. Data AvailabilityAll relevant data are within the manuscript and attached figures and S1 File.
Data Availability

All relevant data are within the manuscript and attached figures and S1 File.
==== Body
pmcIntroduction

For decades, computer-aided drug design (CADD) has held the promise of greatly accelerating drug discovery and bringing much needed pharmacotherapeutic solutions to patients and medical professionals. CADD practitioners like to point to success stories such as the development of raltegravir–an HIV integrase inhibitor [1]–and doramapimod–a highly potent kinase inhibitor–where modeling and simulation have been instrumental in bringing about new therapeutic agents [2]. In the case of raltegravir, all-atom molecular dynamics (MD) simulations [3] revealed a cryptic trench in HIV integrase and led to the discovery of a new class of antiretrovirals–the integrase inhibitors [2]. Similarly, explicit solvent MD simulations revealed a 10 Å shift of the F169 side chain of p38 MAP kinase, exposing a cryptic pocket in the presence of BIRB 796 –a novel ligand, later named doramapimod—which was evaluated by Boehringer Ingelheim in clinical trials for the treatment of inflammatory diseases [4]. As it samples and explores different conformational states, molecular dynamics is particularly well suited to identifying cryptic pockets which remain occluded in many crystal structures [5,6]. While raltegravir and doramapimod are excellent examples of molecular modeling and dynamics delivering new and actionable knowledge that leads to new drugs or drug-like molecules, there are hundreds of types of cancer [7] accounting for millions of deaths yearly [8]. Clearly, much remains to be desired. Indeed, while the rapid turnaround in silico hit discovery [9] provides and the unraveling of biomolecular association [10–16] mechanisms are certainly beneficial to drug design and discovery, there is no shortage of diseases and malignancies in need of new vaccines and treatments. The lack of vaccines for malaria [17] and hepatitis C [18], despite the decades of research, are two prominent examples. Moreover, even in certain cancers where there are therapeutic options available, such as pancreatic cancer, the survival rate is still very low [19], further exacerbating the need for new antineoplastic molecules [20].

One of the enabling–or limiting–factors in the utility of computer modeling and all-atom molecular dynamics simulations [5,21] for the purposes of drug design and discovery is the fidelity and accuracy of the potentials used in these campaigns. On the one end of the computational spectrum, we have molecular docking [22]—a popular and conceptually simple technique used to rapidly screen libraries of ligands against a target of interest. Docking typically employs unsophisticated, easily computable empirical or knowledge-based potentials [23] to evaluate different ligand conformations and rank all the ligands in a library. On the other extreme, we have the theoretically rigorous and computationally demanding absolute binding free energy calculations (ABFEs) [24]. Typically, docking tries to discern ligands that bind the target of interest from ones that do not by rapidly evaluating simple scoring functions while ABFEs try to accurately calculate, at least in principle, the absolute free energy of binding (ΔGbind or simply ΔG) between the target and every ligand being screened by covering the relevant conformational space and evaluating ΔGbind using potential energy functions referred to as force fields. In principle, ΔGbind is the difference between the free energy (G) of the bound state and the free energy of the unbound protein and unbound ligand in solution. In practice, ABFEs do not calculate absolute free energies (Gs) but only free energy changes or differences (ΔGs) by constructing a thermodynamic cycle and decoupling the ligand from the bound state and from solvent and taking the difference of these two terms [25]. The free energy change upon decoupling the ligand from solvent, i.e. desolvating the ligand, is the desolvation free energy or the solvation free energy multiplied by -1. In the case where the solvent is water, what is being calculated is the (de)hydration free energy (HFE). In contrast to ABFEs, where the ligand is decoupled from the system, in relative binding free energy (RBFE) calculations, the ligand is transformed into a different ligand, typically a close analog, in order to calculate the relative change in free energy (ΔΔG) between the two [26–28].

For a chemical compound to be a viable drug in clinical practice, its hydration free energy must be within a suitable range. If a molecule that is too hydrophilic is administered orally, e.g. as a tablet, it will remain mostly in the fluid inside the gastrointestinal tract and fail to distribute throughout the body [29]. Conversely, if a molecule that is too hydrophobic is administered orally, it will also fail to reach its intended target, as it fails to pass into the intestinal fluid [30]. The hydration free energy is a determinant of a drug molecule’s liberation from its pharmaceutical formulation and its absorption and distribution in the body, which are the subject of an entire scientific discipline–pharmacokinetics [31]. Moreover, highly lipophilic molecules are much more likely to be promiscuous binders and tend to behave as pan-assay interference compounds (PAINS) [32]. Further still, highly hydrophilic or hydrophobic compounds are often difficult to work with in terms of chromatography analysis, purification [33], and formulation [34]. Given the importance of HFEs in the development of pharmaceuticals, much effort has been devoted to developing methods for their rapid and accurate prediction [35,36], as evidenced in the SAMPL challenges [37].

(Δ)ΔG can be computed numerically via thermodynamic integration (TI [38]). In TI, the free energy difference between the two states of interest (also referred to as end states) is evaluated by connecting them through some functional form f(λ) of a nonphysical coordinate λ, calculating the derivative of the potential with respect to λ (∂U/∂λ), and estimating the area under the ∂U/∂λ versus λ curve [39–41]. TI makes use of the malleability of the potential energy function which can be calculated for the nonphysical, intermediate, mixed states, which is why it belongs to a class of free energy methods referred to as alchemical transformations. In this way, one ligand can be transmuted into another and the free energy difference between the two can be calculated. The transformation can be done by modifying van der Waals and electrostatic interactions simultaneously (the so-called one-step approach) or separately (the two-step approach) [42]. In the simplest case, f(λ) = λ, which is referred to as linear scaling [43]. To avoid singularities in the potential energy function when particles appear or disappear, i.e. are being inserted into or decoupled from the system, a modified or softcore potential is used [44].

Clearly, the quality of the force fields is one of the major factors in determining the outcomes in such campaigns. The calculated hydration free energy of a ligand is fully determined by two force fields: the water model and the small molecule force field, and their interplay. Over the past several decades, literally hundreds of water models have been proposed and developed, yet in biomolecular modeling, the rigid, fixed-charge, 3-point TIP3P [45] model remains among the most widely used, if not the most widely used [46], despite its many known shortcomings [47]. Similarly, the general Amber force field (GAFF [48]) and its subsequent refinements hold this position in simulating drug and drug-like molecules. Other notable families of force fields are CHARMM [49], the Open Force Field (OpenFF) [50], and OPLS [51]. More recently, a new 4-point rigid water model named OPC [52] was developed by optimizing the point charge distribution in water that reproduces its bulk properties far more accurately than TIP3P; a 3-point version named OPC3 [53] was developed as a compromise between speed and accuracy. It was hoped that an improvement in bulk property reproduction would also translate into an improvement in hydration free energy predictions. For small sets of organic molecules, the new water models did indeed outperform TIP3P [52,54]. However, those studies use only a small number of ligands to carry out the validation and do not report how they select the molecules they use out of the 642 available for benchmarking from the FreeSolv data set [55]–the largest and most comprehensive curated set of experimentally measured hydration free energies.

Here, we report a thorough analysis of the performance of hydration free energy calculations on the entire FreeSolv set for a combination of the GAFF2.11 force field with the TIP3P, OPC, and OPC3 water models, and compare and contrast all-atom molecular dynamics hydration free energy calculations to the more affordable two-dimensional quantitative structure–activity relationship (2D QSAR [56]) approach based on low-level molecular descriptors [57]. We discuss their merits and shortcomings from the standpoint of computer-aided drug discovery and development. Moreover, we also compare the performance of GAFF2.11 with the previous version of GAFF– 1.81 –demonstrating a robust means of validating edits to force fields that is far more relevant to drug discovery, design, and development than what physics focused research laboratories typically employ. We demonstrate that calculated HFEs become increasingly inaccurate with increases in ligand molecular weight and the number of rotatable bonds. Finally, we show that overall statistics like the coefficient of determination (R2) and root mean square error (RMSE) in benchmarking studies can often be misleading and overshadow significant shortcomings. We demonstrate that deeper analysis is needed to reveal and interpret these shortcomings and showcase such an analysis for HFE calculations on the FreeSolv set.

Methods

System setup for thermodynamic integration simulations

To generate starting structures for our simulations, we used the sdf-format ligand structures as provided by FreeSolv. All 642 neutral FreeSolv ligands were processed with antechamber from Amber18 to generate Amber-compatible mol2 files with the appropriate atom types and AM1_BCC partial charges [58]. All ligands were parameterized with GAFF versions 2.11 and 1.81. GAFF2.11-parameterized ligands were solvated with tleap in cubic boxes of TIP3P, OPC, and OPC3 water with a wall distance of 24 Å; GAFF 1.81-parameterized ligands were solvated with TIP3P water only.

Simulation protocol for thermodynamic integration simulations

The solvated systems were subjected to 2000 steps of energy minimization. The systems were then heated from 100 to 300 K over a period of 100 ps at constant volume, followed by 100 ps of density equilibration, 100 ps of constant pressure and temperature (NPT) equilibration, and were finally subjected to NPT production runs of 250 ps under 1 bar, 300 K with periodic boundary conditions. Constant temperature and pressure were maintained with the the Langevin thermostat [59] and Berendsen barostat [60], respectively. Collision frequencies for temperature coupling were 2 ps−1; the pressure relaxation time was set to 2 ps. A cutoff of 12.0 Å was used for nonbonded interactions. Long-range electrostatics beyond the real space cutoff were computed with the particle-mesh Ewald (PME) scheme [61]. The time step was set to 1 fs to keep the simulations stable during alchemical transformations (decoupling); bonds to hydrogen were not constrained.

Thermodynamic integration hydration free energy calculations

Simulations were carried out with one-step thermodynamic integration with softcore potentials and linear scaling [42,44] using the pmemd.cuda MD engine from Amber18 [62]. Ligands were decoupled from the solvation box in 21 evenly spaced λ-windows ranging from 0.0 to 1.0 in intervals of 0.05. The scalpha and scbeta parameters, which control the softness of the potential, were set to 0.5 and 12 Å2, respectively, as in previous work [12]. To avoid any Hamiltonian lag [63], energy minimization, heating, density equilibration, NPT preproduction equilibration, and production dynamics were all carried out with potential energy functions, corresponding to the λ-value of every λ-window. Hydration free energies, error estimates, and convergence metrics [64,65] were computed from the decorrelated dV/dl values (the derivative of the potential with respect to λ) from the Amber output files with the alchemlyb [66] TI estimator.

2D QSAR hydration free energy calculations

Molecular descriptors were calculated with RDKit as described previously [12]. Descriptors were also normalized by molecular weight; parameters with zero variance were removed from consideration, leaving around 280 descriptors in total, which were scaled from 0 to 1. We then trained a Gaussian kernel support vector regressor (SVR [67,68]) with scikit-learn [69] on the experimental hydration free energies. The 642 ligands were split into a training and test set in a 70:30 ratio; the scaled descriptors (the independent variables) and hydration free energies (the dependent variables) for the training compounds were used to build an SVR model to predict hydration free energies; the model was then tested against the HFEs from the test set. To choose the optimal model hyperparameters, during training we performed a grid search with 5-fold crossvalidation for the regularization parameter C varying it from 10−9 to 109 in 10-fold increments and the kernel coefficient gamma using ‘scale’, ’auto’, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 5, 10, 15, and 25 as possible options. The best performing model was then used to predict the HFEs of the test set ligands and the prediction errors (ΔGexperimental− ΔGcalculated, in kcal/mol) and relative errors (ΔGcalculated*100/ΔGexperimental, in per cent (%)) were recorded. This procedure was repeated 10,000 times for random train/test splits, all in a 70:30 ratio.

Results

Thermodynamic integration hydration free energy calculations

The hydration free energy calculations using the TIP3P water model and GAFF2.11 force field resulted in a pleasing R2 of 0.84 and an RMSE of 1.78 ± 0.05 kcal/mol when compared to experiment (Fig 1).

10.1371/journal.pone.0309996.g001 Fig 1 Calculated versus experimental free energies.

a Scatter plot of calculated versus experimental free energies for GAFF2.11 with TIP3P water. The coefficient of determination (R2) and root mean square error (RMSE) are given in the figure legend. Error estimates for experimental and calculated HFEs are given as horizontal and vertical error bars, respectively (note that most of the calculated error estimates are too small to see at this scale). b Scatter plot of calculated versus experimental free energies for GAFF2.11 with OPC water. The R2 and RMSE are given in the figure legend. Error estimates for experimental and calculated HFEs are given as horizontal and vertical error bars, respectively (note that most of the calculated error estimates are too small to see at this scale). c Scatter plot of calculated versus experimental free energies for GAFF2.11 with OPC3 water. The R2 and RMSE are given in the figure legend. Error estimates for experimental and calculated HFEs are given as horizontal and vertical error bars, respectively (note that most of the calculated error estimates are too small to see at this scale). d Error (ΔGexperimental− ΔGcalculated) distributions for the three water models The x = 0 kcal/mol axis indicates the position of perfect predictions; the regions between the vertical solid and dashed lines indicate the location of predictions within 1 and 2 kcal/mol of experimental values, respectively. e Relative error (ΔGcalculated*100/ΔGexperimental) distributions for the three water models. The x = 100% axis indicates the position of perfect predictions; the x = 0% axis separates ligands whose calculated HFE is of the correct sign (right of the axis) from ligands whose HFE sign has been mispredicted (left of the axis).

Somewhat surprisingly, the OPC model produced a lower R2 and a higher RMSE (R2 = 0.67, RMSE = 2.54 ± 0.05 kcal/mol), whereas OPC3 performed similarly to TIP3P (R2 = 0.81, RMSE = 2.18 ± 0.07 kcal/mol). Notably, there is near perfect agreement between the 3-point water models TIP3P and OPC3 (R2 = 0.97, RMSE = 0.97 ± 0.07 kcal/mol, S1 Fig), far more than between OPC and OPC3 (R2 = 0.79, RMSE = 2.02 ± 0.06 kcal/mol). It is also notable that the majority of points lie above the x = y identity line for all three water models. Consequently, the density distribution curves for the prediction errors or residuals (ΔGexperimental− ΔGcalculated) are left-shifted from the x = 0 kcal/mol line (Fig 1D). In the case of a near-perfect prediction, the density distribution curve is tall and narrow, tightly centered around x = 0 kcal/mol; the poorer the predictions, the lower the peak and the wider the distribution. Among the three water models, TIP3P has the tallest and narrowest error distribution. Again, the distribution for the 3-point OPC3 is very similar in shape, only slightly shorter and wider, whereas the OPC distribution is appreciably different with larger errors (residuals) being much more prevalent than with TIP3P and OPC3 where the majority of predictions fall within 2 kcal/mol of the experimental values (the region between the dashed black lines in Fig 1D).

Moreover, we analyze prediction accuracy in terms of relative errors given in percentages as ΔGcalculated*100/ΔGexperimental. This is important as it allows clearly and easily identifying cases where the workflow mispredicts the sign of the hydration free energy (i.e. HFE has been estimated to be positive when it is on fact negative or vice versa), which is a more substantial misprediction than merely getting the value wrong by a certain amount. Such points lie to the left of the x = 0% axis (identified with a solid red line in Fig 1E). In the case of a perfect predictor, all data points would stack at x = 100% (solid black line in Fig 1E), i.e. the prediction perfectly matches experiment in sign and in magnitude. For x > 0% and x ≠ 100%, the predictor has correctly estimated the sign of the HFE, but not its magnitude, whereas for x < 0%, the HFE sign has been mispredicted. In terms of relative errors, the three water models have much more similar distributions with the three modes lying close to each other. Again, TIP3P has the tallest and narrowest distribution, with the second tallest this time being OPC rather than OPC3, the latter two being very similar in shape and position. The TIP3P curve also has the smallest area lying to the left of x = 0%, i.e. TIP3P has the highest number of correctly predicted HFE signs (574/642 or 89%), followed by OPC (547/642 or 85%) and OPC3 (530/642 or 83%).

Further, we present the error distributions as a function of molecular weight (Fig 2A). We see that the molecular mechanics-based predictions exhibit a rightward fanning out of the error distributions, i.e. errors become larger as molecular weight increases. Moreover, we stress that not only does the distribution become wider to the right, the proportion of errors within the -2 to 2 kcal/mol region between the dashed black lines becomes smaller and smaller as one moves to the right of the x axis. Near the beginning of the x axis, the majority of predictions lie within 2 kcal/mol of experimental values. Moving to the right, such predictions become the minority. This means that not only do errors become larger with increasing molecular weight, large errors become more likely with increasing molecular weight. Interestingly, an opposite trend appears when examining relative errors as a function of molecular weight (Fig 2B)–the distribution is wider on the left of the x axis, up to around 200 g/mol, and narrows after that. However, relative errors behave similarly to absolute errors in that smaller residuals are more prevalent in the low-molecular-weight region of the plot and become less likely in the high-weight region.

10.1371/journal.pone.0309996.g002 Fig 2 Hydration free energy prediction errors versus molecular weight.

a GAFF2.11 HFE errors (ΔGexperimental− ΔGcalculated) for the three water models plotted against molecular weight. The regions between the horizontal solid and dashed lines indicate the location of predictions within 1 and 2 kcal/mol of experimental values, respectively. b GAFF2.11 HFE relative errors (ΔGcalculated*100/ΔGexperimental) for the three water models plotted against molecular weight. The y = 100% axis indicates the position of perfect predictions.

For the TIP3P + GAFF2.11 combination, we performed energy minimization, heating, density equilibration, preproduction equilibration, and production dynamics in four independent replicas. Overall results for the hydration free energy estimates, computed from the production stages of the four replicas, are nearly indistinguishable–the 4 replicas agree with each other up to an R2 = 0.99 (S2 Fig). Finally, we compare the latest version of GAFF– 2.11 –to version 1.81 using the TIP3P water model. We find a minor difference in performance–GAFF2.11 produced an R2 of 0.84 and an RMSE of 1.78 ± 0.05 kcal/mol; GAFF1.81 produced and R2 of 0.88 and an RMSE of 1.70 ± 0.05 kcal/mol with a slightly narrower error distribution (S3 Fig).

2D QSAR hydration free energy calculations

We now examine the performance of the support vector machine (SVM) machine learning (ML) models on HFE predictions. We generated 10,000 random train/test splits in a 70:30 ratio which means that each split has 193 ligands in the test set and we have 1,930,000 predictions in total. We generated multiple models from different train/test splits to gauge the performance of such models much more reliably than if we were to use a single model. Overlaying the SVM and GAFF2.11 + TIP3P error distributions plotted against molecular weight highlights the difference in performance (Fig 3). The machine learning models have a narrow, symmetric error distribution tightly centered around y = 0 kcal/mol and y = 100% for the absolute and relative errors, respectively. Notably, the error distributions from the ML predictions do not get wider to the right of the molecular weight axis. This contrasts starkly with the TIP3P ΔG error distribution which becomes wider to the right; conversely, TIP3P relative errors tend to be largest around a molecular weight of 100 g/mol and become narrower to the extremes of the weight range. In contrast to the ML ΔG errors which have a narrow, symmetric, elongated distribution centered on the y = 0 kcal/mol (perfect prediction) axis, TIP3P ΔG errors are downshifted and centered around y = -1 kcal/mol. Similarly, the TIP3P relative error distribution is much wider than its ML counterpart, albeit not as downshifted from the y = 100% (perfect prediction) axis as the ΔG error distribution. Crucially, the ML models mispredicted the signs of the HFEs in around 2% of cases compared to around 11% for the TI-based workflow with TIP3P water, and around 15 and 17% with OPC and OPC3 water, respectively.

10.1371/journal.pone.0309996.g003 Fig 3 Hydration free energy prediction errors versus molecular weight.

a Bivariate kernel density estimate plot of HFE errors (ΔGexperimental− ΔGcalculated) for GAFF2.11 with the TIP3P water model versus support vector machine (SVM) predictors against molecular weight. The regions between the horizontal solid and dashed lines indicate the location of predictions within 1 and 2 kcal/mol of experimental values, respectively. b Bivariate kernel density estimate plot of HFE relative errors (ΔGcalculated*100/ΔGexperimental) for GAFF2.11 with the TIP3P water model and SVM predictors against molecular weight. The y = 100% axis indicates the position of perfect predictions.

Notably, some of the descriptors show significant correlation or anticorrelation (given here as the correlation coefficient R) with the error from HFE calculations (all RDKit descriptors and HFEs are given in S1 File; descriptor definitions can be found in the RDKit documentation and the references therein). When looking at signed errors from the TIP3P calculations, the descriptor with the strongest anticorrelation to error is the VSA_Estate3 molecular surface descriptor [70] normalized by molecular weight (VSA_Estate3/MolWt, R = -0.38). Conversely, the molecular surface descriptors PEOE_VSA14 and SlogP_VSA3 [71] have the strongest positive correlation to HFE error (R = 0.41 and R = 0.46, respectively). The strongest correlations between descriptors and HFE errors, however, appear when taking the absolute values of the errors. In this case, the descriptor with the largest correlation is the number of heteroatoms (R = 0.53) followed by the VSA_Estate3 molecular surface descriptor (R = 0.46).

Discussion

We report a thorough, rigorous benchmark of the TIP3P, OPC, and OPC3 water models with the popular general Amber force field for small molecules. Due to a limitation on the computational resources available to us, we present results based on fairly limited sampling–for every ligand, the HFE estimate is based on 5.25 ns of production dynamics (21 windows x 0.25 ns of sampling per window). Admittedly, this is not sufficient in duration for our automated pipeline to achieve convergence for every ligand. However, as we compare different water models under identical conditions, we can still make fair and valid comparisons and draw meaningful, useful conclusions. In S4 Fig, we plot the distributions of the A_c convergence metric [64,65] from the simulations using the TIP3P, OPC, and OPC3 water models. An A_c value of 1 indicates perfect convergence across the 21 λ windows; a value of 0 indicates that the simulations have not even begun to converge. Interestingly, the TIP3P distribution is by far the right-most, followed by OPC3, demonstrating that the 3-point models have achieved the most convergence given a fixed amount of simulation time. It could thus be expected that 3-point water models reach convergence faster and require less simulation time than 4- and especially 5-point models. Moreover, with all three water models examined here, there is a discernible dependence between molecular weight and convergence–the higher the weight, the lower the A_c metric (S5 Fig), indicating that larger ligands require longer simulations to converge. Further still, there is significant agreement between the convergence estimates among the different water models (S5 Fig).

Another important property that can be expected to strongly affect the convergence in molecular simulations is the number of rotatable bonds (henceforth also referred to as rotors for brevity) a molecule has. In Fig 4, we show a breakdown of the A_c convergence metric and molecular weight by the number of rotatable bonds.

10.1371/journal.pone.0309996.g004 Fig 4 A_c convergence metric and molecular weight distributions broken down by number of rotatable bonds.

a Bivariate kernel density estimate plot of the A_c convergence metric for GAFF2.11 with the TIP3P water model plotted against molecular weight. Ligands have been grouped in 5 separate panels based on the number of rotatable bonds– 0 or 1, 2 or 3, 4 or 5, 6 or 7, and 8 or more. For ligands with 0 or 1 or 6 or 7 rotatable bonds, the molecular weight regions around 400 g/mol have been highlighted with blue and red rectangles, respectively. b Bivariate kernel density estimate plot of the A_c convergence metric for GAFF2.11 with the OPC water model plotted against molecular weight. Ligands have been grouped in 5 separate panels based on the number of rotatable bonds– 0 or 1, 2 or 3, 4 or 5, 6 or 7, and 8 or more. For ligands with 0 or 1 or 6 or 7 rotatable bonds, the molecular weight regions around 400 g/mol have been highlighted with blue and red rectangles, respectively. c Bivariate kernel density estimate plot of the A_c convergence metric for GAFF2.11 with the OPC3 water model plotted against molecular weight. Ligands have been grouped in 5 separate panels based on the number of rotatable bonds– 0 or 1, 2 or 3, 4 or 5, 6 or 7, and 8 or more. For ligands with 0 or 1 or 6 or 7 rotatable bonds, the molecular weight regions around 400 g/mol have been highlighted with blue and red rectangles, respectively.

We see that as the number of rotors increases, the A_c distributions shift downward to lower values, whereas the weight distributions move to the right, indicating that molecules with more rotors tend to have larger molecular weights and tend to exhibit lower convergence scores. With all three water models examined here, nearly all molecules with 8 or more rotatable bonds have an A_c value below 0.4, whereas nearly all ligands with 0 or 1 rotor have an A_c value above 0.4.

We uncover and report a connection between molecular size and calculation error. We focus on molecular weight because it is a direct measure of how much there is to describe in a molecule. However, molecular weight is not the only important factor that should be examined, the number of rotors perhaps being the next logical choice. However, we note that the number of rotatable bonds, being an integer, is a far more coarse descriptor than molecular weight. Moreover, significant chemical complexity can be introduced with few or no rotatable bonds, e.g. with large aromatic or (spiro)cyclic systems. For example, Fig 4 shows that there are compounds with 0 or 1 rotatable bonds but with large molecular weight (around the 400 g/mol region, blue rectangles) and that most of these tend to have lower A_c values than compounds with the same number of rotors but lower weight (the regions above and to the left of the blue rectangles in Fig 4). Moreover, there are ligands with considerably more rotatable bonds and lower weight that have higher convergence metrics (these are the ligands in the 2 or 3, 4 or 5, and 6 or 7 rotors panels lying above the level of the blue rectangles). Finally, Fig 4 shows that there are ligands with 6 or 7 rotatable bonds that also have molecular weights in the 400 g/mol region (red rectangles in the figure). These have similar convergence values to the ligands with 0 or 1 rotor in the 400 g/mol region (blue rectangles) indicating that molecular weight is a better predictor of convergence than the number of rotatable bonds. Indeed, molecular weight is one of the few descriptors that begin to approach significance in correlations to the magnitude of the calculated error (R = 0.38, R = 0.36, and R = 0.43 for the TIP3P, OPC, and OPC3 models, respectively), i.e. the higher the weight, the larger the error in the calculated result, whereas the number of rotors has no correlation to error (R < 0.1 for all three water models). Interestingly, the descriptor that exhibited the highest correlation to the magnitude of the error is the number of heteroatoms–R = 0.53 with TIP3P, R = 0.49 with OPC3; the correlation is much lower for OPC (R = 0.30).

No single descriptor is likely to account for a majority in the variance in the results. Moreover, we find that more sophisticated surface area descriptors such as PEOE_VSA14 and SlogP_VSA3 [71] are among the most highly correlated to error in HFE calculations and that they correlate to error more strongly than simpler descriptors such as total polar area or total molecular surface. PEOE_VSA and SlogP_VSA are the van der Waals atomic surfaces in a molecule that fall within a specific range of partial charge [72] and logP (the octanol/water partition coefficient) values, respectively; they are specific subdivisions of the molecular surface. Our feature engineering efforts, where we normalize properties by molecular weight, also proved productive in that some of the resulting features also showed significant (anti)correlation to HFE errors. One such example is VSA_Estate3 –the electrotopological state [70] within a specific range of van der Waals surface area values–normalized by molecular weight (VSA_Estate3/MolWt, all descriptors are available in S1 File). One could envisage further feature engineering schemes, e.g. normalizing by the number of rotatable bonds or even molecular weight divided by the number of rotatable bonds, to name but two. With the present report, we hope to stimulate further investigations into this matter and, more broadly, to stimulate deeper analysis into results coming out of computational chemistry pipelines. To this end, we make our data fully and freely available (S1 File).

Performing 3 more independent replicas with the TIP3P water model yielded nearly identical results–the 4 replicas agree with each other up to an R2 of 0.99 and three of the four replicas have the same R2 with the experimental HFEs (0.84). Only replica 3 has an R2 of 0.83, indicating that few new ligand conformations are sampled by the additional replicas or that if new conformations are sampled, these are isoenergetic with previously sampled conformations. While it is generally recognized that multiple short simulations are more efficient in covering conformational space than one long simulation [73], our results demonstrate that this general rule has to be applied judiciously–there exists a threshold below which additional sampling becomes inefficient, failing to reach new energy basins from the conformational landscape, all the while adding computational costs.

We emphasize that our goal here is not to draw definitive conclusions about the different force fields and water models and recommend one combination over another but to demonstrate how they should be rigorously and thoroughly interrogated and validated. Ideally, new water models, force fields or versions of force fields should demonstrate appreciably improved performance in terms of error distributions, both absolute and relative, centered around x = 0 kcal/mol and x = 100% (perfect predictions), respectively, being taller and narrower than their predecessor or the previous gold standard, rather than merely improving parameters for one group or another or marginally improving R2 in benchmarking. This does not appear to be the case for GAFF2.11 and 1.81, at least when using TIP3P water (S3 Fig). However, we reserve judgement until sufficient sampling can be attained, and merely suggest how this matter should be interrogated.

While our simulations are fairly short, our results are nearly identical to those of Matos et al. [74], who used a very similar protocol to our own with 20 λ windows and 5 ns of sampling per window– 20-fold greater than what we have used. Running nearly 20 times more computing with TIP3P and GAFF1.7 achieves a modest improvement in accuracy–R2 = 0.87 (from 100 ns of total sampling per ligand) vs 0.84 from our results (5.25 ns of total sampling); RMSE = 1.53 ± 0.07 kcal/mol (from 100 ns of sampling) vs 1.78 ± 0.05 kcal/mol from our work (5.25 ns of sampling). Indeed, our results from 5.25 ns of sampling per ligand are virtually identical to their results from 100 ns of sampling (R2 = 0.97). Further details become evident when comparing their results to the ones we report here with a plot of errors as a function of molecular weight (Fig 5).

10.1371/journal.pone.0309996.g005 Fig 5 Hydration free energy prediction errors versus molecular weight.

HFE errors (ΔGexperimental− ΔGcalculated) from this work and from Matos et al. [74] plotted against molecular weight. The regions between the horizontal solid and dashed lines indicate the location of predictions within 1 and 2 kcal/mol of experimental values, respectively.

We see that in the low-weight region of the plot, the results from their ~100 ns simulations lie closer to the x = 0 kcal/mol error axis (perfect prediction) than our results from ~5 ns of sampling. However, in the high-weight region, the additional sampling has not always produced more accurate results, as is evident from Fig 5. For some of the ligands, the results from Matos et al. are less accurate then our own, despite the 20-fold increase in sampling time. Regrettably, Matos et al. do not provide an estimate for the convergence in their simulations; it would be very interesting to see how convergence scales with simulation time.

Given the fact that 100 ns of sampling are often not enough to achieve accurate results for larger ligands, a conservative estimate would posit a requirement of one microsecond of sampling for ligands with drug- or drug-like molecular weights, which could be loosely defined as 400 g/mol or more. For single-digit numbers of drug-like ligands, this amount of sampling is still trivial in terms of cost. For large sets of ligands, however, it becomes advantageous to turn to ML models, which also appear to have more favorable error profiles–SVM models exhibit a much narrower error distribution across the entire molecular mass range than TI-based estimates. Crucially, the key advantage of ML HFE models is that they exhibit no real dependence between accuracy and ligand molecular weight–the ΔG error distribution for the ML predictors does not broaden moving to the right, unlike the distribution for TI-based estimates.

Conclusions

While herein we focus primarily on molecular weight and to a lesser extent on rotatable bonds, we merely view this work as laying out a blueprint for interrogating results from HFE calculations and in computational chemistry more broadly. That is to say that not only should novel descriptors [75] be explored but also feature engineering with existing descriptors. We examine molecular weight and rotatable bonds only as the first properties that should be looked at while interrogating results from HFE, ABFE or RBFE calculations; we encourage the reader to think carefully about what other properties might be highly relevant to their particular data set(s) and offer potential feature engineering pathways for the reader to explore–normalization by weight, rotors, or both; other avenues could certainly prove fruitful. Further, we hope to inspire more rigorous validation of HFE and MD results in order to avoid situations where overall statistics look favorable, but that is only because the data sets being used are dominated by molecules with properties that are easy to (retro)predict. Fig 5 offers a striking example of this–the overwhelming majority of data points lie in the low weight region below 200 g/mol and have their HFEs predicted to a high degree of accuracy. This makes overall RMSE and R2 values appear favorable, but overshadows the fact that HFE calculations seem to struggle beyond 350 or 400 g/mol–the very region they are most needed. The lack of negative data in published computational chemistry work, particularly ABFE and RBFE studies, is an even more egregious flaw we plan to address in a future publication. If computer-aided drug design is to truly usher in the next generation of pharmaceuticals, it has to go beyond the bay and sail the sea.

Supporting information

S1 Fig Experimental and computed hydration free energies calculated with GAFF2.11 and the TIP3P, OPC, and OPC3 water models.

The density distributions and scatter plots between the four properties along with the corresponding R2 values and equations of best fit are given.

(TIFF)

S2 Fig Experimental and computed hydration free energies calculated with GAFF2.11 and the TIP3P water model in four independent replicas.

The density distributions and scatter plots between the five properties along with the corresponding R2 values and equations of best fit are given.

(TIFF)

S3 Fig Hydration free energy prediction errors.

a HFE error (ΔGexperimental− ΔGcalculated) distributions for GAFF versions 2.11 and 1.81 calculated with the TIP3P water model. The x = 0 kcal/mol axis indicates the position of perfect predictions; the regions between the vertical solid and dashed lines indicate the location of predictions within 1 and 2 kcal/mol of experimental values, respectively. b HFE relative error (ΔGcalculated*100/ΔGexperimental) distributions for GAFF versions 2.11 and 1.81 calculated with the TIP3P water model. The x = 100% axis indicates the position of perfect predictions; the x = 0% axis separates ligands whose calculated HFE is of the correct sign (right of the axis) from ligands whose HFE sign has been mispredicted (left of the axis).

(TIFF)

S4 Fig Convergence metric (A_c) distributions from the TIP3P, OPC, and OPC3 simulations with GAFF2.11.

(TIFF)

S5 Fig Molecular weights for the FreeSolv ligands and convergence metrics (A_c) from their corresponding simulations with the TIP3P, OPC, and OPC3 water models using GAFF2.11.

The density distributions and scatter plots between the four properties along with the corresponding R2 values and equations of best fit are given.

(TIFF)

S1 File A spreadsheet of all molecules used in this study.

For ease of comparison, the parm column gives the molecular ID from the corresponding Matos et al. paper [74], followed by the SMILES, name, the experimental HFE, the experimental uncertainty, the HFE values calculated by Matos et al. and their uncertainties, followed by our results. The ti_A_c_tip3p_gaff2.11 column contains the A_c convergence metric values with GAFF2.11 and TIP3P, the ti_tip3p_gaff2.11 column gives the calculated TI values, and the ti_uncertainty_tip3p_gaff2.11 column contains our error estimates from the calculations. What follows are analogous columns for the calculations with GAFF2.11 and OPC3 and GAFF2.11 with OPC; the GAFF1.81 with TIP3P calculations, followed by the RDKit canonical SMILES, and the RDKit molecular descriptors (their definitions can be found in the RDKit documentation and the references therein), followed by the same descriptors normalized by molecular weight. Finally, the tip3p_errors column lists the errors from the HFE calculations using GAFF2.11 and TIP3P, analogously for opc3_errors and opc_errors; mobley_errors contains the errors from the Matos et al. calculations; the last six columns give the A_c and HFE values from the 3 additional replicas with GAFF2.11 and TIP3P from this work. The author declares that no competing interests exist.

(CSV)

10.1371/journal.pone.0309996.r001
Decision Letter 0
Bhakat Soumendranath Academic Editor
© 2024 Soumendranath Bhakat
2024
Soumendranath Bhakat
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version0
18 Jun 2024

PONE-D-24-21031Calculated hydration free energies become less accurate with increases in molecular weightPLOS 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. Your manuscript has been reviewed by three independent reviewers and they raised valuable questions which require further clarifications and should be reflected in the revised manuscript. I will be delighted to consider the revised version of the manuscript. 

Please submit your revised manuscript by Aug 02 2024 11:59PM. If you will need more time than this to complete your revisions, please reply to this message or contact the journal office at plosone@plos.org. When you're 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.

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). You should upload this letter as a separate file labeled 'Response to Reviewers'.

A marked-up copy of your manuscript that highlights changes made to the original version. You should upload this as a separate file labeled 'Revised Manuscript with Track Changes'.

An unmarked version of your revised paper without tracked changes. You should upload this as a separate file labeled 'Manuscript'.

If you would like to make changes to your financial disclosure, please include your updated statement in your cover letter. Guidelines for resubmitting your figure files are available below the reviewer comments at the end of this letter.

If applicable, we recommend that you deposit your laboratory protocols in protocols.io to enhance the reproducibility of your results. Protocols.io assigns your protocol its own identifier (DOI) so that it can be cited independently in the future. For instructions see: https://journals.plos.org/plosone/s/submission-guidelines#loc-laboratory-protocols. Additionally, PLOS ONE offers an option for publishing peer-reviewed Lab Protocol articles, which describe protocols hosted on protocols.io. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols.

We look forward to receiving your revised manuscript.

Kind regards,

Soumendranath Bhakat

Academic Editor

PLOS ONE

Journal Requirements:

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

1. 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 

https://journals.plos.org/plosone/s/file?id=wjVg/PLOSOne_formatting_sample_main_body.pdf and 

https://journals.plos.org/plosone/s/file?id=ba62/PLOSOne_formatting_sample_title_authors_affiliations.pdf

2. Please note that PLOS ONE has specific guidelines on code sharing for submissions in which author-generated code underpins the findings in the manuscript. In these cases, we expect all author-generated code to be made available without restrictions upon publication of the work. 

Please review our guidelines at https://journals.plos.org/plosone/s/materials-and-software-sharing#loc-sharing-code and ensure that your code is shared in a way that follows best practice and facilitates reproducibility and reuse.

3. We note that the grant information you provided in the ‘Funding Information’ and ‘Financial Disclosure’ sections do not match. 

When you resubmit, please ensure that you provide the correct grant numbers for the awards you received for your study in the ‘Funding Information’ section.

4. Thank you for stating the following in the Acknowledgments Section of your manuscript: 

"Acknowledgments. This study is financed by the European Union-NextGenerationEU through the National Recovery and Resilience Plan of the Republic of Bulgaria, project № BG-RRP-2.004-0004-C01. The in silico calculations were performed at the Centre of Excellence for Informatics and ICT, supported by the Science and Education for Smart Growth Operational Program and co-financed by the European Union through the European Structural and Investment Funds (Grant No. BG05M2OP001-1.001-0003)."

Please note that funding information should not appear in the Acknowledgments section or other areas of your manuscript. We will only publish funding information present in the Funding Statement section of the online submission form. Please remove any funding-related text from the manuscript. 

5. In this instance it seems there may be acceptable restrictions in place that prevent the public sharing of your minimal data. However, in line with our goal of ensuring long-term data availability to all interested researchers, PLOS’ Data Policy states that authors cannot be the sole named individuals responsible for ensuring data access (http://journals.plos.org/plosone/s/data-availability#loc-acceptable-data-sharing-methods).

Data requests to a non-author institutional point of contact, such as a data access or ethics committee, helps guarantee long term stability and availability of data. Providing interested researchers with a durable point of contact ensures data will be accessible even if an author changes email addresses, institutions, or becomes unavailable to answer requests.

Before we proceed with your manuscript, please also provide non-author contact information (phone/email/hyperlink) for a data access committee, ethics committee, or other institutional body to which data requests may be sent. If no institutional body is available to respond to requests for your minimal data, please consider if there any institutional representatives who did not collaborate in the study, and are not listed as authors on the manuscript, who would be able to hold the data and respond to external requests for data access? If so, please provide their contact information (i.e., email address). Please also provide details on how you will ensure persistent or long-term data storage and availability.

6. Please include captions for your Supporting Information files at the end of your manuscript, and update any in-text citations to match accordingly. Please see our Supporting Information guidelines for more information: http://journals.plos.org/plosone/s/supporting-information. 

Additional Editor Comments:

Please provide Github link to scripts and data used in this study

[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: Yes

Reviewer #2: Partly

Reviewer #3: Partly

**********

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

Reviewer #1: Yes

Reviewer #2: No

Reviewer #3: No

**********

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: Yes

Reviewer #2: No

Reviewer #3: No

**********

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: No

Reviewer #3: 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: This manuscript focuses on the impact of molecular properties & choice of water model on hydration free energy – a key component in determination ligand binding free energy. The findings and discussion presented in this study add value to the area of computational drug discovery efforts. I recommend this manuscript for publication after the following comments have been addressed by the author.

1. The author has measured molecular descriptors using RDKit. However, the discussion on the measured descriptors is missing in the discussion. I would recommend adding the list of the descriptors. Did the author find any correlation between one or more of those descriptors and the RMSE in HFE? Did any of those descriptors show correlation with the errors in the sign of the HFE?

2. The results presented in this manuscript show that the 3-point water models outperform the 4-point water model. Could this finding be extended to charged molecules? Should the choice of water model be determined by the physicochemical property of the molecules? A comment on this aspect could be useful if the author has any insights obtained from their analysis.

3. While the author has presented an interesting observation regarding the impact of molecular weight on the accuracy of hydration free energy prediction, it will be useful to add a section on why the errors get larger with increase in molecular weight. Specifically, does it have any correlation with the accessible surface area and total volume of the molecule? Do linear molecules show a different trend than cyclic or branched molecules?

Reviewer #2: See Attached Document for formatted review. Plaintext version pasted below:

Review for :

“Calculated hydration free energies become less accurate with increases in molecular weight”

Summary: This article seeks to address and benchmark long-standing problem of how accurately computation is able to predict the hydration free energy (HFE) of drug-like molecules. By using an established dataset of experimentally measured HFE values, the FreeSolv database, the author uses AMBER MD engine to carry out free energy perturbation methods to compute HFE estimations. The author then generates an machine-learned (ML) model trained on this data and assess the performance of said model relative to HFE predictions. The analysis of accuracy focuses on some select molecular properties, such as molecular weight, but does not do a complete benchmarking of other potentially relevant properties. As such, the article could benefit from a greater in-depth discussion on the results, why certain features were chosen, and more discussion on HFE’s relevance to drug discovery. As discussed below, other parts require further contextualization, and the discussion section requires a major rewrite as it is a single large paragraph.

Major (required) revisions:

1. The introduction would benefit from a greater direct contextualization of HFE calculations and their relevance in drug discovery. While Free Energy Calculations (FECs) are known to be incredibly relevant to drug discovery and inhibitor improvement, the introduction does not make it clear that HFEs also have their own utility. It would be beneficial if the introduction provided additional citations and additional direct contextualization about where HFEs are valuable in the drug discovery process and what point they’re relevant. Obviously, solubility is an important component of the ligand design process, but there is no discussion about additional methods for measuring solubility like logS or how they are utilized.

2. The discussion requires a major rewrite. It is currently one singular large paragraph that spans multiple pages, making it incredibly difficult to read and parse. It is also not clear how the discussion points are connected to the data or are further speculation without the data. It will be important to separate the discussion to clearly demarcate both of those types of paragraphs.

3. Additional discussion on why comparing error to molecular weight and no other properties would benefit the manuscript. In general, I think it is known that HFEs will scale with respect to molecular weight, and with increased ligand size there will be increased parameters to consider and increased time to convergence with simulation-based methods. Thus, it seems consistent with expectation that HFE variance will increase per ligand which would be remedied by increased replicates.

4. Consistent with the above discussion that convergence would require additional sampling with larger molecular weight due to a variety of chemical factors, it would be useful to compare how different numbers of replicates impact this error while still achieving convergence. It would be incredibly useful to identify if there is a scaling for molecular weight, number of atoms/bonds, and the amount of sampling needed.

5. Lastly, it would be useful to provide additional contextualization, testing, and out-of-data comparisons for the Machine Learning model that was constructed. The model tested using these train-test splits never tested against another orthogonal dataset that might be out of distribution. Testing only on data within the dataset allows the model to learn more similar chemistries between the across the train-test split without any guarantee that the model would indeed be able to extrapolate to new chemical topologies.

6. Additionally, a more thorough of the train-test characterization would be useful. It is currently not clear whether the train-test splits contain similar chemical identities in both the training and test sides of the split, which would make it difficult to test how well the dataset is able to extrapolate to new chemistries. Given the importance of the FreeSolv dataset to the ML model, but its small size, it would be useful to do more thorough chemically driven splitting.

7. Given that this was an effort done using established open-source libraries on openly available databases, it would good for there to be an associated github for sharing the results of this data and their weights.

8. Given the interesting nature of the data and the results, a conclusion section would benefit the manuscript and greatly improve readability.

Minor revisions:

1. Citations are needed at the following points:

a. Page 2, line 31, “an HIV integrase inhibitor”

b. Page 2, line 32-33, “modeling and simulation have been instrumental in bringing aobut new therapeutic agents”

c. Page 4, line 86 and 87, “one-step approach” and “two-step approach”

d. Page 5, line 97, “despite its many known shortcomings”

e. Page 5, line 96, “if not the most widely used”

2. The following phrases are not clear and could benefit from rewording:

a. Page 3, line 43-44

3. Start a new paragraph at Page 5 line 108 for clarity.

4. The end of the introduction is riddled with sharp transitions between sentences – some transition words and rephrasing can improve the flow for the reader here.

5. Given the historically outdated nature of the Berendsen Barostat (see papers such as: https://doi.org/10.1016/j.molliq.2022.120116 and https://doi.org/10.1016/j.bbamem.2016.02.004), it would be useful to provide some context for why the Berendsen was used in this simulation over other barostat methods. Alternatively a characterization across different barostats would be useful to see.

6. It would be useful to provide increased contextualization in the text for building the 2D QSAR models

7. Page 8, line 166-167: It would be useful to provide a description what these parameters of zero-variance were that were removed from consideration, and what the other 284 descriptors were in the SI.

8. Page 8, line 169: Please clarify why a 70:30 ratio was chosen for train-test splitting (or provide a citation)

9. Figure would benefit from larger fonts and heading text to improve readability

10. Page 11, line 232: Perhaps a more quantitative description of what it means that the three modes are lying close to each other?

11. Page 14, line 301: Provide a rationale/citation for why 5.25 ns of production dynamics was used, or a more robust sampling was done.

12. Page 15, line 314 appears to have a typo – it should be referring to S5 Fig. if I’m reading correctly?

Reviewer #3: This work compares the experimental Hydration Free Energy (HFE) in the Free Solve dataset to a) values calculated using alchemical free energy methods with thermodynamic integration (TI) estimator and b) to an ML model trained on what I assume is the experimental data. The main conclusion is that the error in the prediction from the alchemical estimation increases with increasing molecular weight.

Overall

There has been considerable effort in producing a useful set of calculations on an important area of computational chemistry, I thank the authors for their efforts. The main points of criticism are that:

1. The implications and reasons for the lack of convergence of the HFE estimates are not adequately explored. This affects the validity of the conclusions drawn.

2. The inclusion of the ML analysis isn’t fully justified.

3. The discussion areas should be a lot more focused on the topic of the paper.

4. There should be more references to relevant work.

There is much here that is interesting and a refocusing of the analysis would be welcome to explore the convergence properties of these calculations.

Introduction

- The introduction is well written and gives an overview of computer aided drug discovery.

- I believe it is too long and not focused enough on the specific area covered by the work and fails to make the case for why this work is necessary.

- The have been numerous works looking at computational predictions of HFE e.g. (non-exhaustive list), The SAMPL challenges or https://doi.org/10.1021/acs.jcim.0c00600, https://pubs.acs.org/doi/abs/10.1021/acs.jcim.0c00285 which should be mentioned.

- There have also been many papers on ML methods for QSAR, see https://paperswithcode.com/sota/molecular-property-prediction-on-freesolv for a ‘leaderboard’ of methods on the Free Solv database.

- Other forcefields e.g., CHARMM small molecule forcefield and the recent Open Force Field, Sage 2.0, and Machine Learned forcefields were also not mentioned.

Methods

- The methods were mostly clearly explained with the following exceptions: The use of the TI estimator was not fully justified given the success of other estimators and potential drawbacks of TI, e.g., MBAR. See https://pubs.acs.org/doi/abs/10.1021/acs.jcim.0c00285, https://doi.org/10.1063/1.5041835, and https://doi.org/10.1021/acs.jced.7b00104.

- The QSAR variables would be better listed in the SI, rather than given as a reference. While implied, the training data was not explicitly identified as the experimental values. Given there is value and precedence in fitting ML models to ABFE data this should be explicitly mentioned.

- The normalization by the molecular weight, while not wrong, was not justified as inclusion of the molecular weight itself should account for the effect of molecular weight on the predictions and regression coefficients.

- The fitting procedure was rigorous but the hyperparameter turning curve should be given in the SI.

- The metric denoted as ‘relative error’ is not consistent with general use of that term (which would include subtraction by 1) please either subtract 1 from the values or use a different term.

Results

The results are generally well reported and clear.

- The vertical lines in figure 1 where not explained in the caption.

- The discussion of the distribution of prediction errors would benefit from being quantitative (using terms like, bias, standard deviation, kurtosis, etc.) rather qualitative (e.g. ‘TIP3P has the tallest and narrowest error distribution’).

- In line 236 please convert these values into percentages.

- Please clarify what you mean by (line 242) ‘not only do the errors become larger with increasing MW but…’ as the error distribution looks to have an approximate mean of 0 for larger weight. Your note of increasing range and kurtosis looks accurate though.

- It is hard to draw meaningful conclusions from figure 2 due to its format (scatter plot with different colours). Plotting an estimate of the mean and range / standard deviation etc. of the errors vs MW would be more informative. One could use a Gaussian process, LOWESS smoother or even categorise the molecular weight into ranges and plot box plots. LOWESS smoothers are available in the Seaborn in the regplot function.

- Figure 3 is quite confusing. It would be ideal if you could keep the format of the comparison the same as Figure 2.

Discussion

- I do not believe that the conclusions you draw here are adequately supported by the data. This is because the convergence of the free energy estimates using the FE method drops significantly with increasing molecular weight.

- It’s not clear why MW has been singled out as the factor influencing accuracy given that number of rotatable bonds must also be very influential.

- I would like to see an analysis of convergence wrt to MW stratified by number of rotatable bonds.

- I would also like some investigation into the reasons for the lack of convergence. It’s not clear what you mean by ‘independent’ replicas in line 324. If they are not different configurations, perhaps perform replicas on some of the least converged molecules with different starting configurations.

- Comparisons between forcefields and water models is valid but only with converged estimates.

- Line 329: sampling multiple short trajectories are only useful if the starting configurations are drawn from the equilibrium configurational distribution (see comment earlier as well).

- In line 362 you say it is not the purpose of the paper to draw definitive conclusions yet the title of the paper is very definitive.

- The discussion is very wide ranging and could do with being shortened and restricted to the main points of the paper.

- The inclusion of the SVM models in the study was not justified.

**********

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

Reviewer #3: 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.]

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 PLOS at figures@plos.org. Please note that Supporting Information files do not need this step.

Attachment Submitted filename: 2024-june-hfe-plosOne.pdf

10.1371/journal.pone.0309996.r002
Author response to Decision Letter 0
Submission Version1
24 Jul 2024

I thank the reviewers for their positive and constructive comments. Below, I address the concerns they have raised.

Reviewer #1: This manuscript focuses on the impact of molecular properties & choice of water model on hydration free energy – a key component in determination ligand binding free energy. The findings and discussion presented in this study add value to the area of computational drug discovery efforts. I recommend this manuscript for publication after the following comments have been addressed by the author.

1. The author has measured molecular descriptors using RDKit. However, the discussion on the measured descriptors is missing in the discussion. I would recommend adding the list of the descriptors. Did the author find any correlation between one or more of those descriptors and the RMSE in HFE? Did any of those descriptors show correlation with the errors in the sign of the HFE?

I have made the complete set of descriptors and calculated HFEs fully and freely available with the supplementary information to this paper. As suggested by the reviewer, I have explored individual descriptors in the Results and Discussion sections.

2. The results presented in this manuscript show that the 3-point water models outperform the 4-point water model. Could this finding be extended to charged molecules? Should the choice of water model be determined by the physicochemical property of the molecules? A comment on this aspect could be useful if the author has any insights obtained from their analysis.

Currently, alchemical HFE calculations such as the ones we describe here are not capable of handling charged ligands because during decoupling the system’s net charge would change. This introduces large errors of a complex nature; the relevant theory is described in more detail in https://doi.org/10.1063/1.4826261. This is why FreeSolv is composed mostly of molecules that have a charge of 0 near neutral pH.

3. While the author has presented an interesting observation regarding the impact of molecular weight on the accuracy of hydration free energy prediction, it will be useful to add a section on why the errors get larger with increase in molecular weight. Specifically, does it have any correlation with the accessible surface area and total volume of the molecule? Do linear molecules show a different trend than cyclic or branched molecules?

I thank the reviewer for this comment, as I feel it is particularly productive. It also aligns with some of the comments made by the other reviewers. Molecular surface descriptors have been explored in more depth in the paper and they do indeed show correlation to the errors in the HFE calculations. Interestingly, the more sophisticated RDKit surface descriptors such as PEOE_VSA and SlogP_VSA have higher correlations than simple descriptors such as total polar surface or total surface area.

Reviewer #2: See Attached Document for formatted review. Plaintext version pasted below:

Review for :

“Calculated hydration free energies become less accurate with increases in molecular weight”

Summary: This article seeks to address and benchmark long-standing problem of how accurately computation is able to predict the hydration free energy (HFE) of drug-like molecules. By using an established dataset of experimentally measured HFE values, the FreeSolv database, the author uses AMBER MD engine to carry out free energy perturbation methods to compute HFE estimations. The author then generates an machine-learned (ML) model trained on this data and assess the performance of said model relative to HFE predictions. The analysis of accuracy focuses on some select molecular properties, such as molecular weight, but does not do a complete benchmarking of other potentially relevant properties. As such, the article could benefit from a greater in-depth discussion on the results, why certain features were chosen, and more discussion on HFE’s relevance to drug discovery. As discussed below, other parts require further contextualization, and the discussion section requires a major rewrite as it is a single large paragraph.

Major (required) revisions:

1. The introduction would benefit from a greater direct contextualization of HFE calculations and their relevance in drug discovery. While Free Energy Calculations (FECs) are known to be incredibly relevant to drug discovery and inhibitor improvement, the introduction does not make it clear that HFEs also have their own utility. It would be beneficial if the introduction provided additional citations and additional direct contextualization about where HFEs are valuable in the drug discovery process and what point they’re relevant. Obviously, solubility is an important component of the ligand design process, but there is no discussion about additional methods for measuring solubility like logS or how they are utilized.

As suggested by the reviewer, I have added an entire paragraph to the Introduction contextualizing HFEs and their relevance in drug design.

2. The discussion requires a major rewrite. It is currently one singular large paragraph that spans multiple pages, making it incredibly difficult to read and parse. It is also not clear how the discussion points are connected to the data or are further speculation without the data. It will be important to separate the discussion to clearly demarcate both of those types of paragraphs.

As suggested by the reviewer, the Discussion has undergone a major rewrite making it shorter and more focused on the work presented here. This comment also nicely aligns with comments by the other reviewers; I have taken this opportunity to add discussion on some of the descriptors to the Discussion section.

3. Additional discussion on why comparing error to molecular weight and no other properties would benefit the manuscript. In general, I think it is known that HFEs will scale with respect to molecular weight, and with increased ligand size there will be increased parameters to consider and increased time to convergence with simulation-based methods. Thus, it seems consistent with expectation that HFE variance will increase per ligand which would be remedied by increased replicates.

I thank the reviewer for this comment as it is particularly constructive. Indeed, my intention for this manuscript was to notify the community of the issue I report and to inspire more investigations into the connection between HFE results (and computational chemistry results in general) and molecular properties. I am well aware that molecular weight is far from the whole story. It is thanks to this comment that I realized that that is never explicitly stated in the text; it is not even suggested. Therefore, I have explicitly stated in the manuscript that it only aims to inspire further analysis; the paper does not claim to be comprehensive. Indeed, this subject can likely fill many more papers-worth of material, likely many books-worth. Here, I have included analysis and discussion on some of the more salient features – weight, number of rotatable bonds, and a few surface area descriptors.

4. Consistent with the above discussion that convergence would require additional sampling with larger molecular weight due to a variety of chemical factors, it would be useful to compare how different numbers of replicates impact this error while still achieving convergence. It would be incredibly useful to identify if there is a scaling for molecular weight, number of atoms/bonds, and the amount of sampling needed.

The reviewer is quite correct to point out that it would be very interesting to see how convergence scales with different factors. Sadly, Matos et al. have not included an estimate for the convergence metric from their 100 ns simulations in the supplementary information of their paper. Therefore, we cannot estimate how convergence scales with simulation time. As for scaling with molecular weight, supplementary figure 5 contains such scatter plots where the R2 between weight and the scaling parameter is around 0.30. However, at this point I am reluctant to make any major claims other than to say there certainly is a relationship. Quantifying it more precisely is the subject of future work, both by myself and other authors. That and the other questions the reviewer poses here would have to be the subject of future investigations as more compute time is presently not available.

5. Lastly, it would be useful to provide additional contextualization, testing, and out-of-data comparisons for the Machine Learning model that was constructed. The model tested using these train-test splits never tested against another orthogonal dataset that might be out of distribution. Testing only on data within the dataset allows the model to learn more similar chemistries between the across the train-test split without any guarantee that the model would indeed be able to extrapolate to new chemical topologies.

I address this comment together with the following one, see below.

6. Additionally, a more thorough of the train-test characterization would be useful. It is currently not clear whether the train-test splits contain similar chemical identities in both the training and test sides of the split, which would make it difficult to test how well the dataset is able to extrapolate to new chemistries. Given the importance of the FreeSolv dataset to the ML model, but its small size, it would be useful to do more thorough chemically driven splitting.

The machine learning section is used only as a story-telling device to introduce the chemical descriptors that will be used to analyze the outcomes from the HFE calculations. Therefore, ML is not the focus of the work.

7. Given that this was an effort done using established open-source libraries on openly available databases, it would good for there to be an associated github for sharing the results of this data and their weights.

All descriptors and HFE results are now included in the supplementary information of this paper for the computational chemistry community to use and analyze. Again, I hope to inspire more papers of this sort, potentially deriving better descriptors, models, and hopefully force fields.

8. Given the interesting nature of the data and the results, a conclusion section would benefit the manuscript and greatly improve readability.

As suggested by the reviewer, a Conclusions section has been added where I explicitly state the goals of this work and its implications.

Minor revisions:

1. Citations are needed at the following points:

a. Page 2, line 31, “an HIV integrase inhibitor”

b. Page 2, line 32-33, “modeling and simulation have been instrumental in bringing aobut new therapeutic agents”

c. Page 4, line 86 and 87, “one-step approach” and “two-step approach”

d. Page 5, line 97, “despite its many known shortcomings”

e. Page 5, line 96, “if not the most widely used”

As requested by the reviewer, all of the above references have been added.

2. The following phrases are not clear and could benefit from rewording:

a. Page 3, line 43-44

This has been simplified.

3. Start a new paragraph at Page 5 line 108 for clarity.

As requested by the reviewer, this is now in a separate paragraph.

4. The end of the introduction is riddled with sharp transitions between sentences – some transition words and rephrasing can improve the flow for the reader here.

The introduction has been rephrased slightly to improve readability.

5. Given the historically outdated nature of the Berendsen Barostat (see papers such as: https://doi.org/10.1016/j.molliq.2022.120116 and https://doi.org/10.1016/j.bbamem.2016.02.004), it would be useful to provide some context for why the Berendsen was used in this simulation over other barostat methods. Alternatively a characterization across different barostats would be useful to see.

Here, I use the Berendsen barostat simply because it is a popular choice and has been a part of our pipeline for a long time. I cannot make any claims about what the results would look like with other barostats, other than to say that I do not expect them to change much. Force fields should be the major determinant of the quality of the calculated results, not the barostat. After all, force fields are the field that is undergoing extensive development, barostats and thermostats – not so much. Presumably, barostats and thermostats are mature enough.

6. It would be useful to provide increased contextualization in the text for building the 2D QSAR models

Again, I do not want to place too much emphasis on the ML section; my primary focus is on the analysis of the HFE results.

7. Page 8, line 166-167: It would be useful to provide a description what these parameters of zero-variance were that were removed from consideration, and what the other 284 descriptors were in the SI.

All descriptors are now available in the SI of the paper.

8. Page 8, line 169: Please clarify why a 70:30 ratio was chosen for train-test splitting (or provide a citation)

70:30 is just a personal preference. Most often researches use 75:25. Indeed, this is the default in many packages, including scikit-learn, and often people use it without even realizing it. I tend to dislike it somewhat because I feel it artificially inflates results – if there is too much data in the training set and too little data in the test set, results may come out looking good simply because there is too little test data to make an error on. However, there is no split that is set in stone, it is entirely up to the user to choose the ratio.

9. Figure would benefit from larger fonts and heading text to improve readability

10. Page 11, line 232: Perhaps a more quantitative description of what it means that the three modes are lying close to each other?

While it would be straightforward to add exact numbers, I believe the text is already fairly heavy on numbers and data. Rather, the goal here is to guide the reader and help them interpret the data rather than to give exact details about the distributions.

11. Page 14, line 301: Provide a rationale/citation for why 5.25 ns of production dynamics was used, or a more robust sampling was done.

The rationale here is quite simple – this is the most I could afford. What is notable is that results differ marginally from 100 ns of sampling, as in the Matos paper (see Fig 5). I stress that the point of the paper is NOT that I can reproduce their results but that my results and their results clearly show a sharp deterioration of performance for larger, more complex ligands. This is something the community will need to acknowledge and address going forward.

12. Page 15, line 314 appears to have a typo – it should be referring to S5 Fig. if I’m reading correctly?

Indeed, the reviewer is correct, this should be S5 Fig. I thank the reviewer for pointing this out.

Reviewer #3: This work compares the experimental Hydration Free Energy (HFE) in the Free Solve dataset to a) values calculated using alchemical free energy methods with thermodynamic integration (TI) estimator and b) to an ML model trained on what I assume is the experimental data. The main conclusion is that the error in the prediction from the alchemical estimation increases with increasing molecular weight.

Overall

There has been considerable effort in producing a useful set of calculations on an important area of computational chemistry, I thank the authors for their efforts. The main points of criticism are that:

1. The implications and reasons for the lack of convergence of the HFE estimates are not adequately explored. This affects the validity of the conclusions drawn.

As requested by reviewers 2 and 3, more discussion on convergence has been added.

2. The inclusion of the ML analysis isn’t fully justified.

This comment aligns with one of the comments from reviewer #2. As I have noted in my response to reviewer #2, the ML section is used only as a story-telling device to introduce the chemical descriptors that will be used to analyze the outcomes from the HFE calculations. Again, ML is not the focus of the work.

3. The discussion areas should be a lot more focused on the topic of the paper.

This comment also nicely aligns with comments from the other reviewers

Attachment Submitted filename: Response to Reviewers.doc

10.1371/journal.pone.0309996.r003
Decision Letter 1
Bhakat Soumendranath Academic Editor
© 2024 Soumendranath Bhakat
2024
Soumendranath Bhakat
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Submission Version1
23 Aug 2024

Calculated hydration free energies become less accurate with increases in molecular weight

PONE-D-24-21031R1

Dear Dr. Ivanov,

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

Within one week, you’ll receive an e-mail detailing the required amendments. When these have been addressed, you’ll receive a formal acceptance letter and your manuscript will be scheduled for publication.

An invoice will be generated when your article is formally accepted. Please note, if your institution has a publishing partnership with PLOS and your article meets the relevant criteria, all or part of your publication costs will be covered. Please make sure your user information is up-to-date by logging into Editorial Manager at Editorial Manager® and clicking the ‘Update My Information' link at the top of the page. If you have any questions relating to publication charges, 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 help maximize its impact. If they’ll be preparing press materials, please inform our press team as soon as possible -- 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.

Kind regards,

Soumendranath Bhakat

Academic Editor

PLOS ONE

Additional Editor Comments (optional):

Dear Dr. Ivanov,

After carefully reviewing the reviewers comments, I am glad to accept your paper PONE-D-24-21031R1 for publication in PLOS ONE. Thanks for your patience during this process.

Best regards,

Soumendranath Bhakat

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: (No Response)

**********

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 author has satisfactorily addressed all the questions. Therefore, I recommend the manuscript for publication.

Reviewer #2: The author has addressed all comments in my previous review and revised the manuscript appropriately.

**********

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

**********

10.1371/journal.pone.0309996.r004
Acceptance letter
Bhakat Soumendranath Academic Editor
© 2024 Soumendranath Bhakat
2024
Soumendranath Bhakat
https://creativecommons.org/licenses/by/4.0/ This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
10 Sep 2024

PONE-D-24-21031R1

PLOS ONE

Dear Dr. Ivanov,

I'm pleased to inform you that your manuscript has been deemed suitable for publication in PLOS ONE. Congratulations! Your manuscript is now being handed over to our production team.

At this stage, our production department will prepare your paper for publication. This includes ensuring the following:

* All references, tables, and figures are properly cited

* All relevant supporting information is included in the manuscript submission,

* There are no issues that prevent the paper from being properly typeset

If revisions are needed, the production department will contact you directly to resolve them. If no revisions are needed, you will receive an email when the publication date has been set. At this time, we do not offer pre-publication proofs to authors during production of the accepted work. Please keep in mind that we are working through a large volume of accepted articles, so please give us a few weeks to review your paper and let you know the next and final steps.

Lastly, if your institution or institutions have a press office, please let them know about your upcoming paper now to help maximize its impact. If they'll be preparing press materials, 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.

If we can help with anything else, please email us at customercare@plos.org.

Thank you for submitting your work to PLOS ONE and supporting open access.

Kind regards,

PLOS ONE Editorial Office Staff

on behalf of

Dr. Soumendranath Bhakat

Academic Editor

PLOS ONE
==== Refs
References

1 Scarsi KK , Havens JP , Podany AT , Avedissian SN , Fletcher C V . HIV-1 Integrase Inhibitors: A Comparative Review of Efficacy and Safety. Drugs. 2020;80 (16 ):1649–76. doi: 10.1007/s40265-020-01379-9 32860583
2 Schames JR , Henchman RH , Siegel JS , Sotriffer CA , Ni H , McCammon JA . Discovery of a Novel Binding Trench in HIV Integrase. J Med Chem. 2004;47 (8 ):1879–81. doi: 10.1021/jm0341913 15055986
3 Hospital A , Goñi JR , Orozco M , Gelpi J . Molecular dynamics simulations: Advances and applications. Adv Appl Bioinforma Chem. 2015;8 (10 ):37–47. doi: 10.2147/AABC.S70333 26604800
4 Frembgen-Kesner T , Elcock AH . Computational Sampling of a Cryptic Drug Binding Site in a Protein Receptor: Explicit Solvent Molecular Dynamics and Inhibitor Docking to p38 MAP Kinase. J Mol Biol. 2006;359 (1 ):202–14. doi: 10.1016/j.jmb.2006.03.021 16616932
5 Durrant J , McCammon JA . Molecular dynamics simulations and drug discovery. BMC Biol. 2011;9 (71 ):1–9. doi: 10.1186/1741-7007-9-71 21214944
6 Oleinikovas V , Saladino G , Cossins BP , Gervasio FL . Understanding Cryptic Pocket Formation in Protein Targets by Enhanced Sampling Simulations. J Am Chem Soc. 2016;138 (43 ):14257–63. doi: 10.1021/jacs.6b05425 27726386
7 Kundra R , Zhang H , Sheridan R , Sirintrapun SJ , Wang A , Ochoa A , et al . OncoTree: A Cancer Classification System for Precision Oncology. JCO Clin Cancer Informatics. 2021;(5 ):221–30. doi: 10.1200/CCI.20.00108 33625877
8 Sung H , Ferlay J , Siegel RL , Laversanne M , Soerjomataram I , Jemal A , et al . Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71 (3 ):209–49. doi: 10.3322/caac.21660 33538338
9 Atanasova M , Dimitrov I , Ivanov S , Georgiev B , Berkov S , Zheleva-Dimitrova D , et al . Virtual Screening and Hit Selection of Natural Compounds as Acetylcholinesterase Inhibitors. Molecules. 2022;27 (10 ):1–19. doi: 10.3390/molecules27103139 35630613
10 Ivanov SM , Huber RG , Warwicker J , Bond PJ . Energetics and Dynamics Across the Bcl-2-Regulated Apoptotic Pathway Reveal Distinct Evolutionary Determinants of Specificity and Affinity. Structure. 2016;24 (11 ):2024–33. doi: 10.1016/j.str.2016.09.006 27773689
11 Ivanov SM , Cawley A , Huber RG , Bond PJ , Warwicker J . Protein-protein interactions in paralogues: Electrostatics modulates specificity on a conserved steric scaffold. PLoS One. 2017;12 (10 ):1–16. doi: 10.1371/journal.pone.0185928 29016650
12 Ivanov SM , Huber RG , Alibay I , Warwicker J , Bond PJ . Energetic Fingerprinting of Ligand Binding to Paralogous Proteins: The Case of the Apoptotic Pathway. J Chem Inf Model. 2019;59 (1 ):245–61. doi: 10.1021/acs.jcim.8b00765 30582811
13 Ivanov SM , Dimitrov I , Doytchinova IA . Bridging solvent molecules mediate RNase A–Ligand binding. PLoS One. 2019;14 (10 ): e0224271. doi: 10.1371/journal.pone.0224271 31644593
14 Ivanov SM , Atanasova M , Dimitrov I , Doytchinova IA . Cellular polyamines condense hyperphosphorylated Tau, triggering Alzheimer’s disease. Sci Rep. 2020;10 (1 ):1–13.31913322
15 Spassov DS , Atanasova M . Inhibitor Trapping in Kinases. Int J Mol Sci. 2024;25 (6 ):3249. doi: 10.3390/ijms25063249 38542228
16 Spassov DS . Binding Affinity Determination in Drug Design: Insights from Lock and Key, Induced Fit, Conformational Selection, and Inhibitor Trapping Models. Int J Mol Sci. 2024;25 (13 ):7124. doi: 10.3390/ijms25137124 39000229
17 Draper SJ , Sack BK , King CR , Nielsen CM , Rayner JC , Higgins MK , et al . Malaria Vaccines: Recent Advances and New Horizons. Cell Host Microbe. 2018;24 (1 ):43–56. doi: 10.1016/j.chom.2018.06.008 30001524
18 The Lancet Gastroenterology & Hepatology. The hunt for a vaccine for hepatitis C virus continues. Lancet Gastroenterol Hepatol. 2021;6 (4 ):253. doi: 10.1016/S2468-1253(21)00073-X 33714362
19 Cabasag CJ , Arnold M , Rutherford M , Bardot A , Ferlay J , Morgan E , et al . Pancreatic cancer survival by stage and age in seven high-income countries (ICBP SURVMARK-2): a population-based study. Br J Cancer. 2022;126 (12 ):1774–82. doi: 10.1038/s41416-022-01752-3 35236937
20 Wu Q , Qian W , Sun X , Jiang S . Small-molecule inhibitors, immune checkpoint inhibitors, and more: FDA-approved novel therapeutic drugs for solid tumors from 1991 to 2021. Vol. 15 , Journal of Hematology and Oncology. BioMed Central; 2022. 1–63 p. doi: 10.1186/s13045-022-01362-9 36209184
21 Bowers KJ , Sacerdoti FD , Salmon JK , Shan Y , Shaw DE , Chow E , et al . Molecular dynamics—Scalable algorithms for molecular dynamics simulations on commodity clusters. Proc 2006 ACM/IEEE Conf Supercomput—SC ‘06 2006:84 .
22 Gaba M , Punam G , Sarbjot S , Gupta G . An overview on Molecular Docking International Journal of Drug Development & Research. Int J Drug Dev Res. 2015;2 (February 2010):2.
23 Bello M , Martínez-Archundia M , Correa-Basurto J . Automated docking for novel drug discovery. Expert Opin Drug Discov. 2013;8 (7 ):821–34.23642085
24 Aldeghi M , Bluck JP , Biggin PC . Absolute alchemical free energy calculations for ligand binding: A beginner’s guide. Vol. 1762 , Methods in Molecular Biology. 2018. 199–232 p.
25 Williams-Noonan BJ , Yuriev E , Chalmers DK . Free Energy Methods in Drug Design: Prospects of “alchemical Perturbation” in Medicinal Chemistry. J Med Chem. 2018;61 (3 ):638–49. doi: 10.1021/acs.jmedchem.7b00681 28745501
26 Zhang S , Hahn DF , Shirts MR , Voelz VA . Expanded Ensemble Methods Can be Used to Accurately Predict Protein-Ligand Relative Binding Free Energies. J Chem Theory Comput. 2021. doi: 10.1021/acs.jctc.1c00513 34516130
27 Cournia Z , Allen B , Sherman W . Relative Binding Free Energy Calculations in Drug Discovery: Recent Advances and Practical Considerations. J Chem Inf Model. 2017;57 (12 ):2911–2937.29243483
28 Hahn DF , Gapsys V , de Groot BL , Mobley DL , Tresadern G . Current State of Open Source Force Fields in Protein-Ligand Binding Affinity Predictions. J Chem Inf Model. 2024. doi: 10.1021/acs.jcim.4c00417 38895959
29 Bowe CL , Mokhtarzadeh L , Venkatesan P , Babu S , Axelrod HR , Sofia MJ , et al . Design of compounds that increase the absorption of polar molecules. Proc Natl Acad Sci U S A. 1997;94 (22 ):12218–23. doi: 10.1073/pnas.94.22.12218 9342389
30 Gaucher G , Fournier E , Le Garrec D , Khalid MN , Hoarau D , Sant V , et al . Delivery of hydrophobic drugs through self-assembling nanostructures. 2004 Int Conf MEMS, NANO Smart Syst ICMENS 2004. 2004;56–7.
31 Fan J , De Lannoy IAM . Pharmacokinetics. Biochem Pharmacol. 2014;87 (1 ):93–120. doi: 10.1016/j.bcp.2013.09.007 24055064
32 Baell JB , Holloway GA . New substructure filters for removal of pan assay interference compounds (PAINS) from screening libraries and for their exclusion in bioassays. J Med Chem. 2010;53 (7 ):2719–40. doi: 10.1021/jm901137j 20131845
33 Condron MM , Monien BH , Bitan G . Synthesis and Purification of Highly Hydrophobic Peptides Derived from the C-Terminus of Amyloid β-Protein. Open Biotechnol J. 2008;2 (1 ):87–93.19898686
34 Vanommeslaeghe K. , Hatcher E. , Acharya C. , Kundu S. , Zhong S. , Shim J. , Darian E. , Guvench O. , Lopes P. , Vorobyov ADMJ I. Techniques to enhance solubility of hydrophobic drugs: An overview. Asian J Pharm. 2016;10 (2 ): S67–75.
35 Scheen J , Wu W , Mey ASJS , Tosco P , Mackey M , Michel J . Hybrid alchemical free Energy/Machine-Learning methodology for the computation of hydration free energies. J Chem Inf Model. 2020;60 (11 ):5331–9. doi: 10.1021/acs.jcim.0c00600 32639733
36 Tsai HC , Tao Y , Lee TS , Merz KM , York DM . Validation of free energy methods in AMBER. J Chem Inf Model. 2020;60 (11 ):5296–300.32551593
37 Bergazin TD , Tielker N , Zhang Y , Mao J , Gunner MR , Francisco K , et al . Evaluation of logP, pKa, and logD predictions from the SAMPL7 blind challenge. Vol. 35, Journal of Computer-Aided Molecular Design. Springer International Publishing; 2021. 771–802 p.
38 Jorge M , Garrido NM , Queimada AJ , Economou IG , MacEdo EA . Effect of the integration method on the accuracy and computational efficiency of free energy calculations using thermodynamic integration. J Chem Theory Comput. 2010;6 (4 ):1018–27.20467461
39 Giese TJ , York DM . A GPU-Accelerated Parameter Interpolation Thermodynamic Integration Free Energy Method. J Chem Theory Comput. 2018;14 (3 ):1564–82. doi: 10.1021/acs.jctc.7b01175 29357243
40 Klimovich P V. , Shirts MR, Mobley DL. Guidelines for the analysis of free energy calculations. J Comput Aided Mol Des. 2015;29 (5 ):397–411.25808134
41 Lee T-S , Allen BK , Giese TJ , Guo Z , Li P , Lin C , et al . Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery. J Chem Inf Model. 2020. doi: 10.1021/acs.jcim.0c00613 32936637
42 Steinbrecher T , Joung I , Case DA . Soft-core potentials in thermodynamic integration: Comparing one-and two-step transformations. J Comput Chem. 2011;32 (15 ):3253–63. doi: 10.1002/jcc.21909 21953558
43 Zou J , Tian C , Simmerling C . Blinded prediction of protein-ligand binding affinity using Amber thermodynamic integration for the 2018 D3R grand challenge 4. J Comput Aided Mol Des. 2018;(Md):1–24.29204945
44 Hornak V , Simmerling C . Development of softcore potential functions for overcoming steric barriers in molecular dynamics simulations. J Mol Graph Model. 2004;22 (5 ):405–13. doi: 10.1016/j.jmgm.2003.12.007 15099836
45 Jorgensen WL , Chandrasekhar J , Madura JD , Impey RW , Klein ML . Comparison of simple potential functions for simulating liquid water. J Chem Phys. 1983;79 (1983 ):926.
46 Kadaoluwa Pathirannahalage SP , Meftahi N , Elbourne A , Weiss ACG , McConville CF , Padua A , et al . Systematic Comparison of the Structural and Dynamic Properties of Commonly Used Water Models for Molecular Dynamics Simulations. J Chem Inf Model. 2021;61 (9 ):4521–36. doi: 10.1021/acs.jcim.1c00794 34406000
47 Kiss PT , Baranyai A . Sources of the deficiencies in the popular SPCE and TIP3P models of water. J Chem Phys. 2011;134(5). doi: 10.1063/1.3548869 21303091
48 Wang JM , Wolf RM , Caldwell JW , Kollman PA , Case DA . Development and Testing of a General Amber Force Field. J Comput Chem. 2004;25 (9 ):1157–74. doi: 10.1002/jcc.20035 15116359
49 Vanommeslaeghe K , Hatcher E , Acharya C , Kundu S , Zhong S , Shim J , et al . CHARMM General Force Field (CGenFF): A force field for drug-like molecules compatible with the CHARMM all-atom additive biological force fields. J Comput Chem. 2009;31 (4 ):671–90.
50 Boothroyd S , Behara PK , Madin OC , Hahn DF , Jang H , Gapsys V , et al . Development and Benchmarking of Open Force Field 2.0.0: The Sage Small Molecule Force Field. J Chem Theory Comput. 2023;19 (11 ):3251–75. doi: 10.1021/acs.jctc.3c00039 37167319
51 Lu C , Wu C , Ghoreishi D , Chen W , Wang L , Damm W , et al . OPLS4: Improving force field accuracy on challenging regimes of chemical space. J Chem Theory Comput. 2021;17 (7 ):4291–300. doi: 10.1021/acs.jctc.1c00302 34096718
52 Izadi S , Anandakrishnan R , Onufriev A V . Building water models: A different approach. J Phys Chem Lett. 2014;5 (21 ):3863–71. doi: 10.1021/jz501780a 25400877
53 Izadi S , Onufriev A V . Accuracy limit of rigid 3-point water models. J Chem Phys. 2016;145(7). doi: 10.1063/1.4960175 27544113
54 Vassetti D , Pagliai M , Procacci P . Assessment of GAFF2 and OPLS-AA General Force Fields in Combination with the Water Models TIP3P, SPCE, and OPC3 for the Solvation Free Energy of Druglike Organic Molecules. J Chem Theory Comput. 2019;15 (3 ):1983–95. doi: 10.1021/acs.jctc.8b01039 30694667
55 Mobley DL , Guthrie JP . FreeSolv: A database of experimental and calculated hydration free energies, with input files. J Comput Aided Mol Des. 2014;28 (7 ):711–20. doi: 10.1007/s10822-014-9747-x 24928188
56 Lewis RA , Wood D . Modern 2D QSAR for drug discovery. Wiley Interdiscip Rev Comput Mol Sci. 2014;4 (6 ):505–22.
57 Xue L , Bajorath J . Molecular Descriptors in Chemoinformatics, Computational Combinatorial Chemistry, and Virtual Screening. Comb Chem High Throughput Screen. 2012;3 (5 ):363–72.
58 Jakalian A , Bush BL , Jack DB , Bayly CI . Fast, efficient generation of high-quality atomic charges. AM1-BCC model: I. Method. J Comput Chem. 2000;21 (2 ):132–46.
59 Adelman SA , Doll JD . Generalized Langevin equation approach for atom/solid‐surface scattering: Collinear atom/harmonic chain model. J Chem Phys. 1974;61 (10 ):4242–5.
60 Berendsen HJ . C, Postma JPM, van Gunsteren WF, DiNola A, Haak JR. Molecular dynamics with coupling to an external bath. J Chem Phys. 1984;81 (1984 ):3684–90.
61 Darden T , York D , Pedersen L . Particle mesh Ewald: An N⋅log(N) method for Ewald sums in large systems. J Chem Phys. 1993;98 (1993 ):10089.
62 Kaus JW , Pierce LT , Walker RC , McCammon JA . Improving the efficiency of free energy calculations in the Amber molecular dynamics package. J Chem Theory Comput. 2013;9 (9 ):1–8. doi: 10.1021/ct400340s 26589006
63 Kollman P. Free-Energy Calculations—Applications to Chemical and Biochemical Phenomena. Chem Rev. 1993;93 (7 ):2395–417.
64 Fan S , Iorga BI , Beckstein O . Prediction of octanol-water partition coefficients for the SAMPL6- log P molecules using molecular dynamics simulations with OPLS-AA, AMBER and CHARMM force fields. J Comput Aided Mol Des. 2020;34 (5 ):543–60.31960254
65 Fan S , Nedev H , Vijayan R , Iorga BI , Beckstein O . Precise force-field-based calculations of octanol-water partition coefficients for the SAMPL7 molecules. J Comput Aided Mol Des. 2021;35 (7 ):853–70. doi: 10.1007/s10822-021-00407-4 34232435
66 Chodera JD . A Simple Method for Automated Equilibration Detection in Molecular Simulations. J Chem Theory Comput. 2016;12 (4 ):1799–805. doi: 10.1021/acs.jctc.5b00784 26771390
67 Platt J , others. Probabilistic outputs for support vector machines and comparisons to regularized likelihood methods. Adv large margin Classif. 1999;10 (3 ):61–74.
68 Chang CC , Lin CJ . LIBSVM: A Library for support vector machines. ACM Trans Intell Syst Technol. 2011;2 (3 ):1–40.
69 Pedregosa F , Varoquaux G , Gramfort A , Michel V , Thirion B , Grisel O , et al . Scikit-learn: Machine learning in Python. J Mach Learn Res. 2011; 12 :2825–30.
70 Hall LH , Mohney B , Kier LB . The electrotopological state: structure information at the atomic level for molecular graphs. J Chem Inf Model. 1991;31 (1 ):76–82.
71 Labute P. A widely applicable set of descriptors. J Mol Graph Model. 2000;18 (4–5 ):464–77. doi: 10.1016/s1093-3263(00)00068-1 11143563
72 Gasteiger J , Marsili M . Iterative partial equalization of orbital electronegativity-a rapid access to atomic charges. Tetrahedron. 1980;36 (22 ):3219–28.
73 Knapp B , Ospina L , Deane CM . Avoiding False Positive Conclusions in Molecular Simulation: The Importance of Replicas. J Chem Theory Comput. 2018;14 (12 ):6127–38. doi: 10.1021/acs.jctc.8b00391 30354113
74 Duarte Ramos Matos G , Kyu DY , Loeffler HH , Chodera JD , Shirts MR , Mobley DL . Approaches for Calculating Solvation Free Energies and Enthalpies Demonstrated with an Update of the FreeSolv Database. J Chem Eng Data. 2017;62 (5 ):1559–69. doi: 10.1021/acs.jced.7b00104 29056756
75 Dutt R , Madan AK . Development and application of novel molecular descriptors for predicting biological activity. Med Chem Res. 2017;26 (9 ):1988–2006.
