Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems

Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems View ORCID Profile Shivani Gupta , Qinfang Sun , View ORCID Profile Ronald M. Levy doi: https://doi.org/10.1101/2025.10.19.683321 Shivani Gupta 1 Center for Biophysics and Computational Biology, Temple University , Philadelphia, PA, USA 2 Department of Chemistry, Temple University , Philadelphia, PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Shivani Gupta Qinfang Sun 1 Center for Biophysics and Computational Biology, Temple University , Philadelphia, PA, USA 2 Department of Chemistry, Temple University , Philadelphia, PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site Ronald M. Levy 1 Center for Biophysics and Computational Biology, Temple University , Philadelphia, PA, USA 2 Department of Chemistry, Temple University , Philadelphia, PA, USA Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Ronald M. Levy For correspondence: ronlevy{at}temple.edu Abstract Full Text Info/History Metrics Preview PDF ABSTRACT Mutation analysis of single and double mutations of proteins including for the two widely studied proteins: Staphylococcus nuclease (S. nuclease) and T4 lysozyme, provide critical insights into their stability and fitness. Furthermore, such analysis facilitates a deeper understanding of the complex interplay between a protein’s sequence, structural conformation, and functional output. Free Energy Perturbation (FEP) is an alchemical, all-atom molecular dynamics based computational approach that determines the free energy change ( ΔΔG ) from wild type to mutant states. Two widely adopted software platforms used for this purpose are Schrödinger and GROMACS. We have compared the results of FEP simulations for mutations using these platforms, employing the OPLS4 force field implemented in Schrödinger, and the Amber99SB-ILDN force field implemented in GROMACS for this work. For the 38 single mutants of S. nuclease, the Pearson r between the experimental and the calculated free energy change ( ΔΔG ) was 0.86 using Schrödinger and 0.87 using GROMACS. The reliability of the FEP method using the two software platforms with the specified force fields was further demonstrated by its performance on 24 single mutants of T4 lysozyme, yielding strong correlation between predicted and experimental ΔΔG values, with Pearson r of 0.80 for Schrödinger and 0.85 for GROMACS. Additionally, the computed folding free energy changes for 45 double mutations in S. nuclease ( ) using Schrödinger correlated well with both experimental measurements (Pearson r = 0.74) and previously reported GROMACS values (Pearson r = 0.71). Correspondingly, the nonadditivities( ) of the double mutations derived using Schrödinger for these 45 double mutants of S. nuclease were also found to be in good agreement with the experimental values (Pearson r = 0.79), as well as with the previously reported GROMACS results (Pearson r = 0.61). A good correlation was also observed between computed values from Schrödinger and GROMACS, with a Pearson r of 0.71 for double mutants and 0.61 for their nonadditivities. Collectively, these findings establish the efficiency, accuracy, and reliability of Schrödinger FEP+ and GROMACS in predicting mutation induced free energy changes in S. nuclease and T4 lysozyme. The integration of FEP based computational methodologies with experimental validation provides a framework for quantifying mutation induced changes in protein thermostability, which can be used as a tool for protein engineering and design. INTRODUCTION Mutational analysis serves as a powerful tool for elucidating the complex interdependencies among protein sequence, structure, and function. 1 By systematically assessing the impact of specific amino acid substitutions on protein thermostability 2 and consequently on functional activity, researchers can uncover fundamental principles governing protein folding, molecular interactions, and evolutionary fitness. 3 , 4 This approach provides critical insights into the intrinsic trade-offs between stability and function, thereby enabling rational protein engineering and aiding the interpretation of pathogenic or disease associated variants. The impact of mutations on protein thermostability has been extensively studied across diverse protein systems. 5 – 9 A protein’s ability to remain properly folded and functionally active at elevated temperatures is governed by its thermostability, making it a critical parameter for understanding degradation, aggregation, and denaturation over time. Thermostable proteins are of significant industrial relevance due to their broad-spectrum applications. Their pivotal roles as biocatalysts in chemical synthesis, components in biofuel production, agents in food processing, and additives in detergents are inherently dependent on their structural and functional stability under thermal stress. 10 In the pharmaceutical industry, thermostable therapeutics such as monoclonal antibodies, enzymes, and cytokines exhibit enhanced stability during manufacturing, storage, and transportation, making them more suitable for clinical and commercial applications. 11 Accurate prediction of protein thermostability through computational approaches provide a potentially efficient complement to experimental techniques. 12 In this context, various physics based approaches such as Molecular Mechanics Generalized Born Surface Area (MM/GBSA) 13 , 14 and Free Energy Perturbation 15 – 17 methods have gained attention. Unlike data driven machine learning models, which may be constrained by the quality, bias, or representativeness of their training datasets, these methods are founded on fundamental physical principles and provide predictive insights that are independent of empirical training data. Free energy perturbation (FEP) 18 – 21 is an alchemical free energy simulation method, firmly rooted in statistical thermodynamics that can be used to estimate free energy differences, provided sufficiently long simulation times, together with enhanced sampling techniques and an accurate potential energy function are employed. It is an explicit solvent, molecular dynamics (MD) based method, offering a framework for evaluating the thermodynamic impact of amino acid substitutions on protein stability. The term ‘alchemical’ refers to the generation of a series of nonphysical intermediate states connecting two physical end states typically the wild-type (initial) and mutant (final) forms, facilitating the estimation of the difference between the free energy to fold the wild-type and mutant proteins. The transformation from wild-type to mutant residues, for single mutants, has been performed in both the folded and unfolded states of the protein that follows a thermodynamic cycle, as illustrated in Figure 1 of the Methodology section. Similarly, the thermodynamic cycle representing double mutations and their corresponding nonadditivity has been depicted in Figure 2 of the Methodology section. Download figure Open in new tab Figure 1. Thermodynamic Cycle to calculate Stability Change of Single Mutant T22C ( ΔΔ G). The green color highlights wild-type residue and blue represents the mutated residue, with rest colored as red, being same in both states. The folding free energy change upon single amino acid mutation is ΔΔ G = Δ G Mut - Δ G WT - Δ G 1 - Δ G 2 Download figure Open in new tab Figure 2. Double mutant cycle with three different pathways (blue, orange, and green) to calculate the double1nutation free energy change and the corresponding nonadditivity ( δ ). Here in , the Subscript WT = Reference/initial state, Superscript A = Mutated/Final state and same for others. In principle free energy change fro1n all three routes should be the same. If then mutation A and B are nonadditive, highly corelated and exhibit a therrnodyna1nic coupling. Staphylococcal nuclease (S. nuclease) is a small extracellular protein of 149 amino acid residue length with no cysteines, secreted by Staphylococcus aureus . Its ability to degrade both DNA and RNA in the presence of free Ca 2+ ions, finds various applications in nucleic acid chemistry since 1956. 22 , 23 It’s simple solution behavior being soluble in both native and denatured states has made it a benchmark system used for studies on protein folding 24 , 25 and NMR spectroscopy. 26 The experimental values for both single and double mutants in S. nuclease were sourced from the extensive work of Shortle et al in S. nuclease. 27 – 29 A second benchmark protein we focus on in the present study is T4 lysozyme (164 residues), 30 Schellman and co-workers 31 have measured the entropy, enthalpy and free energy change for the wildtype protein and a large number of mutants. 32 A diverse set of experimentally characterized single mutants for T4 lysozyme, along with their corresponding ΔΔG values indicating changes in protein stability, has been compiled from the literature. 33 – 41 and considered for the present study. The present study investigates a total of 62 single mutants (SMs) across two proteins: S. nuclease (38 mutants) and T4 lysozyme (24 mutants), and 45 double mutations of S. nuclease using the Schrödinger (FEP+) and GROMACS packages. These codes are both widely used for studying free energy changes associated with protein-ligand binding in explicit solvent, and more recently for studying the effects of mutations on protein stability. There have been few efforts to directly compare the performance of these codes on the same benchmark systems. This study provides the first systematic comparison of free energy change predictions obtained from the Schrödinger and GROMACS simulation platforms across a comprehensive set of single and double mutations in S. Nuclease and a set of single mutations in T4 lysozyme, for which extensive experimental stability data are available. Despite methodological distinctions and differences in the underlying potential energy functions between the two platforms, both yield results that are broadly consistent with experimental measurements and in close agreement with each other. METHODS Dataset Selection and Preparation Single mutants for two proteins under investigation: S. nuclease and T4 lysozyme were curated from the studies by Duan et al (2020) 42 and de Groot et al. (2021). 43 The initial dataset comprised of 27 single mutants (SMs) of S. nuclease (PDB ID: 1EY0) and 26 single mutants of T4 lysozyme (PDB IDs: 2LZM and 1L63), with all experimental measurements performed at pH 7 ± 1. An additional 18 SMs from the dataset associated with PDB: 1STN 43 were incorporated, augmenting the initial 27 mutations and resulting in a total of 45 single mutations analyzed for S. nuclease. After excluding three overlapping mutations, the final dataset comprised of 42 unique SMs. Among these, only 38 mutations were successfully simulated using the Schrödinger and GROMACS platforms, while the remaining four were omitted due to technical issues involving proline residues. Furthermore, a set of 45 double mutants (DMs) 43 for the S. nuclease was also analyzed. These DMs were simulated using only Schrödinger to estimate the free energy changes ( ΔΔG ) and assess their corresponding nonadditive effects. For T4 lysozyme, although the initial dataset contained 26 SMs, two proline-related mutations were excluded due to similar technical constraints, resulting in a final set of 24 single mutations. No double mutant simulations were performed for T4 lysozyme owing to the lack of reliable experimental data for benchmarking. The experimental values were sourced from comprehensive studies by Shortle (for S. nuclease) 27 – 29 and Matthews (for T4 lysozyme), 34,35,37–41 as mentioned above, serving as benchmark for evaluating the calculated free energy changes from Schrödinger and GROMACS simulations across SM and DM datasets for these two proteins. Before discussing the methodology for FEP calculations using Schrödinger and GROMACS, it is important to first understand the underlying thermodynamic cycle and the associated coupling parameter that are common to both platforms. Description of the Alchemical Thermodynamic Cycles The alchemical thermodynamic cycles provide a conceptual framework for the transformations from the wild-type to single and double mutants, as illustrated in Figure 1 and Figure 2 , respectively. The corresponding nonadditive effects associated with the double mutants are also depicted in Figure 2 . The transformations include both physical and alchemical (non-physical) states and since free energy is a state function, its value depends only on the initial and final states and hence independent of the pathway taken. The objective is to compute the folding free energy change resulting from a single mutation, transitioning from the wild-type (WT) to the mutant (Mut) state. Rather than directly simulating the folding process for each state, which is computationally intensive, an alchemical approach has been employed. As illustrated in Figure 1 , the physical pathways represent the folding of the WT and Mut state from their respective unfolded structure, with the associated folding free energy denoted as ΔG WT (path WT) and ΔG Mut (path Mut) respectively. The alchemical pathways consist of the free energy associated with mutating WT to Mut in the folded state denoted as ΔG 1 (path 1) and in the unfolded state denoted as ΔG 2 (path 2). Since the free energy is a state function, the sum of the ΔG around the thermodynamic cycle is zero. Consequently, the effect of the mutation on protein stability ( ΔΔG ) can be estimated using Equation 1 , which relates the differences in folding free energies of the WT and mutant states via the alchemical transformations. The unfolded state of the protein, for both the wild-type and mutant forms, is represented by a capped tripeptide, which preserves the original sequence environment by placing the residue of interest in the central position. The use of a capped tripeptide offers a more realistic and chemically accurate representation for alchemical transformations compared to an uncapped peptide, by mimicking the local backbone environment and minimizing end effects. 44 For double mutants (DMs) an additional thermodynamic cycle has been constructed, as shown in Figure 2 . In this cycle, the transformation from the wild-type (WT) to the double mutant (Mut AB) final state can proceed via three distinct routes. In route 1, the transformation proceeds from WT to Mut A, followed by a second mutation to reach state Mut AB. In Route 2, the path goes from WT to Mut B, and then to Mut AB. Route 3 represents a direct transformation from WT to Mut AB, where both mutations A and B are introduced simultaneously. Since free energy is a state function, the total free energy change associated with the double mutation, denoted as ( in Equation 2 ), must be independent of the transformation path. Therefore, all three routes should theoretically yield identical values. However, in practice, Route 3 (the direct transformation) is employed because experimentally determined PDB structures for the mutated states A or B are not available In the case of double mutations, the combined effect may deviate from the sum of the individual single mutation effects, a phenomenon referred to as nonadditivity ( in Equation 3 ). When , the two mutations are energetically coupled and therefore nonadditive. Conversely, if , the mutations are additive and independent. Within the framework of alchemical FEP, the contribution from the unfolded state is assumed to be additive. Consequently, it cancels out during the calculation of nonadditivity, simplifying the analysis. As a result, the nonadditivity is evaluated solely based on the free energy differences between the single and double mutations in the folded state as formally represented in Equation 3 . Coupling parameter (λ) In alchemical free energy simulations, the coupling parameter (λ) is a dimensionless variable used to denote a transition between two thermodynamic states commonly, state A (wild-type) and state B (mutant). 18 As λ progresses from 0 to 1, the system’s Hamiltonian is gradually modified in a series of steps from H A to H b . Mathematically, the interpolated Hamiltonian is represented in Equation 4 . This parameter governs the pattern of the transformation and serves as a foundational component for computing free energy differences using FEP methods. With the theoretical framework established, the following sections provide a detailed description of the methodologies employed on each platform. First, we outline the workflow implemented in Schrödinger, followed by the methodology employed in GROMACS for FEP calculations. FEP calculations using Schrödinger (FEP+) The FEP+ module from the Schrödinger Release 2024-3 (FEP+, Schrödinger, LLC, New York, NY, 2024) was employed to perform FEP calculations. Initial protein structures for different proteins were retrieved from the Protein Data Bank (PDB) and prepared using Maestro’s Protein Preparation Wizard. The crystal water molecules were retained, hydrogen atoms were added, and hydrogen-bonding networks were optimized. PropKa was used to predict the ionization states of residues at pH 7.0. Subsequently, the structures underwent restrained minimization of heavy atoms using the OPLS4 force field, 45 and the process was terminated once the root mean square deviation (RMSD) of heavy atoms reached 0.3 Å. Each system was solvated using a cubic SPC water box (default water model), extending 5 Å from the protein in all directions. To simulate physiological conditions, counterions Na + or Cl − were added to both folded and unfolded systems, followed by the addition of excess Na + and Cl − ions to achieve a final salt concentration of 0.15 M NaCl. All ions were randomly distributed throughout the simulation box. For unfolded models, capped tripeptides extracted from the crystal structures of the corresponding native proteins, containing the mutation site at the central position in the peptide were used. 42 The capping groups consisted of an acetyl group at the N-terminus and an N-methyl group at the C-terminus. The protein FEP+ calculations were performed using a dual topology algorithm, in which multiple MD simulations were performed to transform the system from wild-type (WT) to the mutated state. During this alchemical transformation, the electrostatic and van der Waals interactions of the WT side chain gradually turned off, while those of the mutant side chain were simultaneously turned on. The mutated residue was included in the Replica Exchange Solute Tempering (REST2) region, 46 , 47 which effectively raises the temperature for enhanced sampling of that residue. The FEP+/REST2 framework combines solute tempering with replica exchange to improve sampling of sidechain conformational flexibility near the mutation site, typically within 5 Å. The solvated systems were then relaxed and equilibrated using the default protocol of the Desmond Molecular Dynamics System (D. E. Shaw Research, New York, NY) 48 as implemented in Maestro. This protocol involves a series of energy minimizations and short MD simulations with restraints. Each perturbation was simulated over 16 to 24 λ windows (where λ denotes the alchemical intermediate states or coupling factor as explained above), depending on the mutation type. The charge conserving mutations were simulated using 16 λ windows while charge altering mutations used 24 λ windows. Electrostatic and van der Waals transformations were integrated across the same λ schedule. Each λ-window was simulated for 10 ns which was extended up to 100 ns to ensure convergence for challenging cases. For charge changing mutations, a co-alchemical water approach was employed to mitigate long range electrostatic artifacts. The electrostatic interactions were handled using the Particle Mesh Ewald (PME) approach 49 , 50 with a real-space cutoff of 1.0 nm. Van der Waals interactions were also gradually shifted to zero at this same cutoff distance. Final production simulations were conducted with the default ensemble using the Desmond engine and free energy differences between alchemical states were computed using Multistate Bennett Acceptance Ratio (MBAR) 51 in Schrödinger. FEP calculations using GROMACS The alchemical FEP calculations for both protein systems were also performed using GROMACS 2022.6 52 , 53 which included 38 SMs from S. nuclease and 23 SMs from T4 lysozyme. The dual topology and the hybrid coordinate files for the mutated protein were generated using the pmx server. 54 Similarly, the tripeptide structures, extracted from the full length protein using Maestro, were processed by pmx server to obtain their corresponding dual topology representations and hybrid tripeptides. The AMBER99sb * ILDN force field 55 – 57 has been applied for these transformation in both the cases as this force field has been widely adopted for FEP calculations, as supported by several studies. 15 , 58 The crystal waters from the protein structure were retained during topology processing of the hybrid protein. Each hybrid system (protein and tripeptide) was solvated in a TIP3P 59 water box with a 10 Å buffer. The counterions Na + or Cl − were added to both the systems to neutralize the charge and achieve the final salt concentration of 0.15 M, consistent with standard FEP simulation practices. The retention of crystallographic water and the use of 0.15 M salt concentration align with established FEP protocols in both Schrödinger and GROMACS frameworks. Then the hybrid protein and tripeptide systems underwent energy minimization using the steepest descent algorithm across 27 λ windows (from λ= 0 to λ = 26. The 27 λ windows were partitioned into electrostatic decoupling (10 λ windows) and van der Waals decoupling (17 λ windows), with soft-core potentials applied to enable smooth interpolation of non-bonded interactions. This was followed by an additional minimization using the l-bfgs integrator for further relaxation. Subsequently, the systems were equilibrated under the NVT ensemble for 10 ps using the velocity-rescaling thermostat. 60 This step was followed by NPT equilibration ensemble, performed sequentially using the Berendsen and Parrinello-Rahman barostats 61 each for 30 ps. The electrostatic interactions were treated using the Particle Mesh Ewald (PME) method 49 , 50 and a real space cutoff of 1.0 nm was applied. The van der Waals interactions were gradually shifted to zero at the same cutoff distance (1.0 nm). All bonds involving hydrogen atoms were constrained using the LINCS 50 algorithm (order = 6). Following equilibration, replicas were simulated for 50 ns with NPT ensemble for both the hybrid protein and the tripeptide systems. The free energy differences between the wild-type and the mutant forms were computed using the Bennett Acceptance Ratio (BAR) 62 , 63 method, implemented via the gmx_bar module in GROMACS. In summary, to assess the impact of mutations on protein free energy landscapes, we performed a comparative analysis of single and double mutations FEP protocols implemented in GROMACS and Schrödinger’s FEP+ framework. Although both platforms utilize alchemical transformation methods, they differ markedly in system setup, topology handling, and sampling strategies. This study aims to evaluate the consistency and reliability of predicted ΔΔG values across these two distinct, widely used, computational platforms. RESULTS and DISCUSSION The experimental determination of protein thermostability changes due to mutations typically using chemical or thermal denaturation techniques remains the gold standard for validation. However, these methods are often time-consuming, resource-intensive, and low-throughput, making them more difficult for large-scale mutational screening. Moreover, while such methods yield valuable thermodynamic data, they offer limited insight into the atomic level interactions underlying protein stability changes. In contrast, computational methods such as FEP, which leverage all-atom MD simulations, potentially provide a faster, scalable, and cost-effective alternative provided they have sufficient accuracy. These methods not only enable high-throughput in-silico mutagenesis but also deliver atomistic insights into stability mechanisms, allowing simulations of various protein conformational states to better understand folding landscapes and thermostability. Recent advances in high performance computing including GPU acceleration as well as improvements in force fields and simulation algorithms, have enhanced the predictive accuracy of these computational approaches. In this work, we employed two widely used MD platforms: Schrödinger and GROMACS, for FEP simulations, each offering distinct advantages and limitations. GROMACS is an open-source, freely available MD engine noted for its flexibility, extensive customizability, and support for a broad range of force fields including AMBER, CHARMM, OPLS and GROMOS. However, its implementation for FEP workflows requires manual setup, command-in-line proficiency, and greater time investment, making it less accessible to non-expert users. In contrast, Schrödinger FEP+ is a commercial, closed-source platform that supports only the modified OPLS force field. While it offers less flexibility, its user-friendly graphical interface, built-in FEP workflow as FEP+ module, and enhanced sampling methods such as REST (Replica Exchange with Solute Tempering) streamline the simulation process and significantly reduce setup complexity and computational overhead. Despite the widespread use of both platforms for FEP-based protein stability predictions, a systematic comparison of their performance using experimental data as a benchmark has been lacking. In this work, we address this gap by directly evaluating and comparing free energy changes for single and double mutations using Schrödinger and GROMACS against experimental thermostability data, thereby providing insights into their relative accuracy, usability, and computational efficiency. The free energy changes ( ΔΔG ) for 38 single mutations of S. nuclease were predicted using the two MD platforms: Schrödinger and GROMACS. These values are provided in the Supporting Information (Table S1). Correlation plots comparing experimentally determined free energy changes with predictions calculated using Schrödinger and GROMACS are shown in Figure 3A and 3B respectively. Both platforms demonstrated state-of-the-art accuracy in reproducing experimental measurements, with Pearson correlation coefficients of 0.86 for Schrödinger ( Figure 3A ) and 0.87 for GROMACS ( Figure 3B ) between calculated and experimental ΔΔG values, indicating strong linear agreement. Additionally, the mean absolute difference also referred to as the average unsigned error (AUE), was further calculated to evaluate predictive accuracy. The AUE values were found to be 0.99 kcal/mol for Schrödinger and 0.61 kcal/mol for GROMACS, both within the generally accepted threshold of ∼1 kcal/mol for reliable free energy predictions. Notably, these AUE values are comparable to those reported for FEP+ protocol in predicting small molecule protein relative binding affinities, 64 , 65 supporting the robustness of our results. The strong Pearson r of 0.87 between the previously reported 42 and our calculated folding free energy changes ( ΔΔG ) highlights the reliability and robustness of these calculations, as shown in Figure S1 and Table S1. The high correlation and low AUE observed in this study are consistent with previously reported FEP studies, 42 , 66 where correlations between computed and experimental ΔΔG values ranged from 0.64 to 0.82, further validating the reliability of both platforms. Overall, these findings establish both Schrödinger and GROMACS as accurate and reliable tools for predicting protein thermostability changes induced by single mutations using FEP-based approaches. Download figure Open in new tab Figure 3. Correlation plots of the free energy changes for 38 single mutations for S. nuclease between the experi111ental values and calculated ones fro111 (A) Schrodinger and (B) GROMACS. Pearson co1Telation coefficient (r) and coefficient of deten11ination (R 2 ) are indicated. The data plotted here is provided in Table S1. The analysis was extended to evaluate folding free energy change for 45 double mutants (DMs) of S. nuclease represented as . The values calculated using Schrödinger and previously reported values using GROMACS 43 alongside experimental measurements, which continue to serve as the benchmark, were summarized in Table S2. The calculated using Schrödinger yielded a Pearson correlation coefficient of 0.74 when compared with the experimental data ( Figure 4A ). A similarly strong correlation (r = 0.71) was observed between the previously reported GROMACS values and experimental results, as shown in Figure 4B . The AUE was found to be 0.79 kcal/mol for Schrödinger and 0.76 kcal/mol for GROMACS, both within the generally accepted threshold of 1 kcal/mol. Moreover, we compared previously reported GROMACS values with those calculated using Schrödinger, yielding a Pearson correlation of 0.71 (Figure S2). Comparable accuracy and low AUE values demonstrate that both Schrödinger (OPLS4 force field) and GROMACS (AMBER99sb * ILDN force field) gave consistent and reliable results for double mutations in the S. nuclease protein. Download figure Open in new tab Figure 4. Correlation plots of the free energy changes for S. nuclease’s 45 double mutants between the experi1nental values and calculated ones with (A) Schrodinger and (B) GROMACS. Pearson correlation coefficient (r) and coefficient of detennination (R 2 ) are indicated. The data plotted here is provided in Table S2. In the case of double mutants (DMs), the free energy changes can deviate significantly from the sum of the individual single mutations, a phenomenon known as nonadditivity ( ) . The corresponding nonadditivities for the 45 DMs have been calculated employing the FEP+ protocol in Schrödinger. For comparison, previously reported free energy changes calculated using GROMACS 43 were included. Experimental values continue to serve as the benchmark for this study, consistent with the analyses of single and double mutations. Schrödinger predictions of nonadditivity yielded a Pearson r of 0.79 with experimental data ( Figure 5A ), whereas the previously reported GROMACS values showed a correlation of 0.63 ( Figure 5B ). The AUE was found to be 0.31 kcal/mol for Schrödinger and 0.35 kcal/mol for GROMACS relative to experimental values. Additionally, a Pearson r of 0.61 was observed for the nonadditivity values ( ) of 45 double mutants of S. nuclease between previously reported GROMACS predictions and our calculated Schrödinger values, as illustrated in Figure S3. All these nonadditivity values obtained from experimental measurements, GROMACS and Schrödinger calculations were tabulated in Table S3. Download figure Open in new tab Figure 5. Correlation plots for the corresponding nonadditivity for S. nuclease’s 45 double mutants between experi,nental values and calculated ones with (A) Schrodinger and (B) GROMACS. Pearson correlation coefficient (r) and coefficient of detennination (R 2 ) are indicated. The data plotted here is provided in Table S3. We have extended our analysis further to another widely studied protein system, T4 lysozyme in order to gain more confidence in thermostability predictions using computational approaches. The dataset consists of 24 single mutations for T4 lysozyme protein, for which free energy changes ( ΔΔG ) were calculated using both the FEP+ module of Schrödinger and the GROMACS. However, double mutants were not included for the T4 lysozyme system, as experimental values were not available. To assess the agreement between calculated and experimental values, correlation analyses were performed. Figure 6A shows the correlation between experimental ΔΔG values and those calculated using Schrödinger (Pearson r = 0.85), while Figure 6B presents the correlation between experimental values and GROMACS predictions (Pearson r = 0.80). The corresponding AUEs were 0.63 kcal/mol for Schrödinger and 0.55 kcal/mol for GROMACS. A direct comparison of computed ΔΔG values with Schrödinger and GROMACS yielded a Pearson r of 0.81 ( Figure 6C ). Furthermore, a Pearson r of 0.96 was observed between previously reported 42 Schrödinger results and the ΔΔG values calculated using Schrödinger in this study (Figure S4). All free energy changes: experimental, calculated via GROMACS and Schrödinger, and previously reported Schrödinger values were summarized in Table S4. Collectively, these results demonstrate the effectiveness in modeling free energy changes arising from single mutations in T4 lysozyme 42 by GROMACS and Schrödinger. Download figure Open in new tab Figure 6. Correlation plots of the free energy changes for T4 lysozy1ne 24 single mutations between the experi1nental and calculated values with (A) Schrodinger and (B) GROMACS and (C) Schrodinger vs GROMACS. Pearson correlation coefficient (r) and coefficient of detennination (R 2 ) are indicated. The data plotted here is provided in Table S4. CONCLUSION FEP methods are increasingly employed to estimate changes in protein thermostability caused by mutations. In this study, we compared two widely used MD platforms: GROMACS and Schrödinger’s FEP+, to perform FEP-based calculations for single mutants of S. nuclease and T4 lysozyme proteins. Additionally, FEP-based predictions were extended to DMs of S. nuclease, including the assessment of their corresponding nonadditivity effects. Our results showed strong correlation between experimental and computed free energy changes ( ΔΔG ) for 62 single mutants in two different proteins using GROMACS and Schrödinger. Schrödinger accurately captured the free energy changes and their corresponding nonadditivities ( ) for 45 double mutants of S. nuclease protein, showing good correlation with experimental values. A strong Pearson correlation for free energy changes and nonadditivity for 45 double mutants of S. nuclease, has also been observed between the previously reported GROMACS values and those calculated in this study using Schrödinger. The findings from the present study demonstrate the robustness of both Schrödinger and GROMACS in modeling single and double mutations, accurately predicting their free energy changes and nonadditivity. The strong correlation between experimental and computed free energy changes further highlights the consistency and accuracy of these platforms in capturing the thermodynamic effects of mutations. The choice between two platforms often depends on factors such as the complexity of the system under study, the cost of the software involved, preference towards a specific force field or enhanced sampling technique, and of course user familiarity. While FEP is not a substitute for experimental validation, it offers a rapid, scalable, and mechanistically informative approach to mutation analysis in proteins. Importantly, FEP enables atomic level insights into the effects of mutations, which are not directly accessible through experimental denaturation assays alone. Therefore, integrating FEP-based computational predictions (either by GROMACS or Schrödinger) with experimental measurements presents a powerful and complementary strategy for protein engineering and thermostability optimization. Supporting Information The supporting information includes detailed folding free energy changes ( ΔΔG ) for 38 single mutants (SMs) of the S. nuclease protein (Table S1), comparing experimental values, previously reported data using Schrödinger or GROMACS, and values calculated in this study using both platforms. Figure S1 illustrates correlation plot of free energy changes ( ΔΔG ) for 38 SMs in the S. nuclease protein, comparing previously reported values using Schrödinger or GROMACS with values calculated in this study using Schrödinger. Pearson correlation coefficient (r) and coefficient of determination (R 2 ) are indicated. Table S2 gives information on folding free energy changes ( ) in kcal/mol for 45 double mutants (DMs) of the S. nuclease protein. The experimental values and previously reported values using GROMACS alongside values calculated in this study using Schrödinger are shown for comparison. Figure S2 illustrates correlation plot of free energy changes ( ) for 45 DMs in the S. nuclease protein, comparing previously reported values from GROMACS with values calculated in this study using Schrödinger. Pearson correlation coefficient (r) and coefficient of determination (R 2 ) are shown. Table S3 gives nonadditivity values ( ) in kcal/mol for 45 double mutants (DMs) of the S. nuclease protein. The experimental and previously reported values obtained using GROMACS are shown alongside values calculated in this study using Schrödinger for comparison. Figure S3 illustrates correlation plot of nonadditivity ( ) for 45 DMs in the S. nuclease protein, comparing previously reported values from GROMACS with values calculated in this study using Schrödinger. Pearson correlation coefficient (r) and coefficient of determination (R 2 ) are shown. Table S4 highlights Folding free energy changes ( ΔΔG ) in kcal/mol for 24 single mutants (SMs) of the T4 lysozyme protein. The experimental values and previously reported values using Schrödinger alongside values calculated in this study using GROMACS and Schrödinger are provided for comparison. Figure S4 illustrated correlation plot of free energy changes ( ΔΔG ) for 24 SMs in the T4 lysozyme protein, comparing previously reported values with calculated values from Schrödinger. Pearson correlation coefficient (r) and coefficient of determination (R 2 ) are shown. Acknowledgements This work was supported in part by NIH Grant R35 GM132090. Footnotes Competing Interest Statement: The authors declare no competing interests. References 1. ↵ Soskine , M. & Tawfik , D. S. Mutational effects and the evolution of new protein functions . Nat. Rev. Genet . 11 , 572 – 582 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 2. ↵ Nisthal , A. , Wang , C. Y. , Ary , M. L. & Mayo , S. L. Protein stability engineering insights revealed by domain-wide comprehensive mutagenesis . Proc. Natl. Acad. Sci . 116 , 16367 – 16377 ( 2019 ). OpenUrl Abstract / FREE Full Text 3. ↵ Mannige , R. Dynamic New World: Refining Our View of Protein Structure, Function and Evolution . Proteomes 2 , 128 – 153 ( 2014 ). OpenUrl PubMed 4. ↵ Dill , K. A. & MacCallum , J. L. The protein-folding problem, 50 years on . Science vol. 338 1042 – 1046 at doi: 10.1126/science.1219021 ( 2012 ). OpenUrl Abstract / FREE Full Text 5. ↵ Pack , S. P. & Yoo , Y. J. Protein thermostability: structure-based difference of amino acid between thermophilic and mesophilic proteins . J. Biotechnol . 111 , 269 – 277 ( 2004 ). OpenUrl CrossRef PubMed 6. Duan , J. , Nilsson , L. & Lambert , B. Structural and functional analysis of mutations at the human hypoxanthine phosphoribosyl transferase (HPRT1) locus . Hum. Mutat . 23 , 599 – 611 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 7. Stefl , S. , Nishi , H. , Petukh , M. , Panchenko , A. R. & Alexov , E. Molecular Mechanisms of Disease-Causing Missense Mutations . J. Mol. Biol . 425 , 3919 – 3936 ( 2013 ). OpenUrl CrossRef PubMed 8. Petukh , M. , Kucukkal , T. G. & Alexov , E. On Human Disease-Causing Amino Acid Variants: Statistical Study of Sequence and Structural Patterns . Hum. Mutat . 36 , 524 – 534 ( 2015 ). OpenUrl CrossRef PubMed 9. ↵ Ponnuswamy , P. K. , Muthusamy , R. & Manavalan , P. Amino acid composition and thermal stability of proteins . Int. J. Biol. Macromol . 4 , 186 – 190 ( 1982 ). OpenUrl CrossRef Web of Science 10. ↵ Rigoldi , F. , Donini , S. , Redaelli , A. , Parisini , E. & Gautieri , A. Review: Engineering of thermostable enzymes for industrial applications . APL Bioeng . 2 , ( 2018 ). 11. ↵ McConnell , A. D. et al. A general approach to antibody thermostabilization . MAbs 6 , 1274 – 1282 ( 2014 ). OpenUrl CrossRef PubMed 12. ↵ Scarabelli , G. , Oloo , E. O. , Maier , J. K. X. & Rodriguez-Granillo , A. Accurate Prediction of Protein Thermodynamic Stability Changes upon Residue Mutation using Free Energy Perturbation . J. Mol. Biol . 434 , 167375 ( 2022 ). OpenUrl CrossRef PubMed 13. ↵ Wang , C. et al. Investigation of endosome and lysosome biology by ultra pH-sensitive nanoprobes . Adv. Drug Deliv. Rev . 113 , 87 – 96 ( 2017 ). OpenUrl CrossRef PubMed 14. ↵ Srinivasan , J. , Miller , J. , Kollman , P. A. & Case , D. A. Continuum Solvent Studies of the Stability of RNA Hairpin Loops and Helices . J. Biomol. Struct. Dyn . 16 , 671 – 682 ( 1998 ). OpenUrl CrossRef PubMed Web of Science 15. ↵ Gapsys , V. , Michielssens , S. , Seeliger , D. & de Groot , B. L. Accurate and Rigorous Prediction of the Changes in Protein Free Energies in a Large-Scale Mutation Scan . Angew. Chemie Int . Ed. 55 , 7364 – 7368 ( 2016 ). OpenUrl 16. Dang , L. X. , Merz , K. M. & Kollman , P. A. Free energy calculations on protein stability: Thr-157 .fwdarw. Val-157 mutation of T4 lysozyme . J. Am. Chem. Soc . 111 , 8505 – 8508 ( 1989 ). OpenUrl 17. ↵ Jespers , W. et al. QresFEP: An Automated Protocol for Free Energy Calculations of Protein Mutations in Q . J. Chem. Theory Comput . 15 , 5461 – 5473 ( 2019 ). OpenUrl PubMed 18. ↵ Jorgensen , W. L. & Ravimohan , C. Monte Carlo simulation of differences in free energies of hydration . J. Chem. Phys . 83 , 3050 – 3054 ( 1985 ). OpenUrl CrossRef Web of Science 19. Beveridge , D. L. & DiCapua , F. M. Free Energy Via Molecular Simulation: Applications to Chemical and Biomolecular Systems . Annu. Rev. Biophys. Biophys. Chem . 18 , 431 – 492 ( 1989 ). OpenUrl CrossRef PubMed Web of Science 20. Pohorille , A. , Jarzynski , C. & Chipot , C. Good Practices in Free-Energy Calculations . J. Phys. Chem. B 114 , 10235 – 10253 ( 2010 ). OpenUrl CrossRef PubMed 21. ↵ Zwanzig , R.W. High-Temperature Equation of State by a Perturbation Method. I. Nonpolar Gases . J. Chem. Phys . 22 , 1420 – 1426 ( 1954 ). OpenUrl CrossRef PubMed Web of Science 22. ↵ Flanagan , J. M. , Kataoka , M. , Shortle , D. & Engelman , D. M. Truncated staphylococcal nuclease is compact but disordered . Proc. Natl. Acad. Sci . 89 , 748 – 752 ( 1992 ). OpenUrl Abstract / FREE Full Text 23. ↵ Shortle , D. Staphylococcal Nuclease: A Showcase of m-Value Effects . in Advances in Protein Chemistry vol. 46 217 – 247 ( 1995 ). OpenUrl 24. ↵ Tucker , P. W. , Hazen , E. E. & Cotton , F. A. Staphylococcal nuclease reviewed: A prototypic study in contemporary enzymology . Mol. Cell. Biochem . 23 , 131 – 141 ( 1979 ). OpenUrl PubMed 25. ↵ Anfinsen , C. B. Principles that Govern the Folding of Protein Chains . Science (80- . ). 181 , 223 – 230 ( 1973 ). OpenUrl FREE Full Text 26. ↵ NUCLEAR MAGNETIC RESONANCE AS A STRUCTURAL METHOD by Oleg Jarde fzky and Norma G . Wade-Jarde fzky . 27. ↵ Shortle , D. , Stites , W. E. & Meeker , A. K. Contributions of the Large Hydrophobic Amino Acids to the Stability of Staphylococcal Nuclease . Biochemistry 29 , 8033 – 8041 ( 1990 ). OpenUrl CrossRef PubMed Web of Science 28. Green , S. M. , Meeker , A. K. & Shortle , D. Contributions of the Polar, Uncharged Amino Acids to the Stability of Staphylococcal Nuclease: Evidence for Mutational Effects on the Free Energy of the Denatured State . Biochemistry 31 , 5717 – 5728 ( 1992 ). OpenUrl CrossRef PubMed Web of Science 29. ↵ Green , S. M. & Shortle , D. Patterns of Nonadditivity between Pairs of Stability Mutations in Staphylococcal Nuclease-*-Only a Limited Number of Studies of Paired Mutant Effects on Protein Stability Have Been Reported, with the General Conclusion Being That Effects Often Appear to B . Shortle & Meeker vol. 32 https://pubs.acs.org/sharingguidelines ( 1993 ). 30. ↵ Alber , T. & Matthews , B. W. Structure and thermal stability of phage T4 lysozyme . In Methods in Enzymology vol. 154 511 – 533 ( 1987 ). OpenUrl 31. ↵ Schellman , J. A. , Lindorfer , M. , Hawkes , R. & Grutter , M. Mutations and protein stability . Biopolymers 20 , 1989 – 1999 ( 1981 ). OpenUrl CrossRef PubMed Web of Science 32. ↵ Becktel , W. J. & Baase , W. A. Thermal denaturation of bacteriophage T4 lysozyme at neutral pH . Biopolymers 26 , 619 – 623 ( 1987 ). OpenUrl CrossRef PubMed Web of Science 33. ↵ Cornish , V. W. , Kaplan , M. I. , Veenstra , D. L. , Kollman , P. A. & Schultz , P. G. Stabilizing and Destabilizing Effects of Placing β-Branched Amino Amino Acids in Protein α-Helixes . Biochemistry 33 , 12022 – 12031 ( 1994 ). OpenUrl PubMed 34. Nicholson , H. , Tronrud , D. E. , Becktel , W. J. & Matthews , B. W. Analysis of the effectiveness of proline substitutions and glycine replacements in increasing the stability of phage T4 lysozyme . Biopolymers 32 , 1431 – 1441 ( 1992 ). OpenUrl CrossRef PubMed Web of Science 35. Matsumura , M. , Becktel , W. J. & Matthews , B. W. Hydrophobic stabilization in T4 lysozyme determined directly by multiple substitutions of Ile 3 . Nature 334 , 406 – 410 ( 1988 ). OpenUrl CrossRef PubMed 36. Heinz , D. W. et al. Accommodation of Amino Acid Insertions in an α-Helix of T4 Lysozyme . J. Mol. Biol . 236 , 869 – 886 ( 1994 ). OpenUrl CrossRef PubMed Web of Science 37. Daopin , S. , Sauer , U. , Nicholson , H. & Matthews , B. W. Contributions of Engineered Surface Salt Bridges to the Stability of T4 Lysozyme Determined by Directed Mutagenesis . Biochemistry 30 , 7142 – 7153 ( 1991 ). OpenUrl CrossRef PubMed 38. Bell , J. A. , Becktel , W. J. , Sauer , U. , Baase , W. A. & Matthews , B. W. Dissection of Helix Capping in T4 Lysozyme by Structural and Thermodynamic Analysis of Six Amino Acid Substitutions at Thr 59 . Biochemistry 31 , 3590 – 3596 ( 1992 ). OpenUrl CrossRef PubMed Web of Science 39. Matthews , B. W. , Nicholson , H. & Becktel , W. J. Enhanced protein thermostability from site-directed mutations that decrease the entropy of unfolding . Proc. Natl. Acad. Sci. U. S. A . 84 , 6663 – 6667 ( 1987 ). OpenUrl Abstract / FREE Full Text 40. Heinz , D. W. , Baase , W. A. & Matthews , B. W. Folding and function of a T4 lysozyme containing 10 consecutive alanines illustrate the redundancy of information in an amino acid sequence . Proc. Natl. Acad. Sci. U. S. A . 89 , 3751 – 3755 ( 1992 ). OpenUrl Abstract / FREE Full Text 41. ↵ Nicholson , H. , Anderson , D. E. , Dao-pin , S. & Matthews , B. W. Analysis of the interaction between charged side chains and the alpha-helix dipole using designed thermostable mutants of phage T4 lysozyme . Biochemistry 30 , 9816 – 28 ( 1991 ). OpenUrl CrossRef PubMed Web of Science 42. ↵ Duan , J. , Lupyan , D. & Wang , L. Improving the Accuracy of Protein Thermostability Predictions for Single Point Mutations . Biophys. J . 119 , 115 – 127 ( 2020 ). OpenUrl CrossRef PubMed 43. ↵ Werner , M. , Gapsys , V. & de Groot , B. L. One Plus One Makes Three: Triangular Coupling of Correlated Amino Acid Mutations . J. Phys. Chem. Lett . 12 , 3195 – 3201 ( 2021 ). OpenUrl CrossRef PubMed 44. ↵ Markthaler , D. , Kraus , H. & Hansen , N. Overcoming Convergence Issues in Free-Energy Calculations of Amide-to-Ester Mutations in the Pin1-WW Domain . J. Chem. Inf. Model . 58 , 2305 – 2318 ( 2018 ). OpenUrl PubMed 45. ↵ Lu , C. et al. OPLS4: Improving Force Field Accuracy on Challenging Regimes of Chemical Space . J. Chem. Theory Comput . 17 , 4291 – 4300 ( 2021 ). OpenUrl CrossRef PubMed 46. ↵ Wang , L. et al. Modeling Local Structural Rearrangements Using FEP/REST: Application to Relative Binding Affinity Predictions of CDK2 Inhibitors . J. Chem. Theory Comput . 9 , 1282 – 1293 ( 2013 ). OpenUrl PubMed 47. ↵ Wang , L. , Berne , B. J. & Friesner , R. A. On achieving high accuracy and reliability in the calculation of relative protein–ligand binding affinities . Proc. Natl. Acad. Sci . 109 , 1937 – 1942 ( 2012 ). OpenUrl Abstract / FREE Full Text 48. ↵ Bowers , K. J. et al. Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters . in ACM/IEEE SC 2006 Conference (SC’06) 43 – 43 ( IEEE , 2006 ). doi: 10.1109/SC.2006.54 . OpenUrl CrossRef 49. ↵ Darden , T. , York , D. & Pedersen , L. Particle mesh Ewald: An N · log( N ) method for Ewald sums in large systems . J. Chem. Phys . 98 , 10089 – 10092 ( 1993 ). OpenUrl CrossRef PubMed Web of Science 50. ↵ Essmann , U. et al. A smooth particle mesh Ewald method . J. Chem. Phys . 103 , 8577 – 8593 ( 1995 ). OpenUrl CrossRef Web of Science 51. ↵ Shirts , M. R. & Chodera , J. D. Statistically optimal analysis of samples from multiple equilibrium states . J. Chem. Phys . 129 , ( 2008 ). 52. ↵ Van Der Spoel , D. et al. GROMACS: Fast, flexible, and free . J. Comput. Chem . 26 , 1701 – 1718 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 53. ↵ Kohnke , B. , Kutzner , C. & Grubmüller , H. A GPU-Accelerated Fast Multipole Method for GROMACS: Performance and Accuracy . J. Chem. Theory Comput . 16 , 6938 – 6949 ( 2020 ). OpenUrl PubMed 54. ↵ Gapsys , V. & de Groot , B. L. pmx Webserver: A User Friendly Interface for Alchemistry . J. Chem. Inf. Model . 57 , 109 – 114 ( 2017 ). OpenUrl CrossRef PubMed 55. ↵ Lindorff-Larsen , K. et al. Improved side-chain torsion potentials for the Amber ff99SB protein force field . Proteins Struct. Funct. Bioinforma . 78 , 1950 – 1958 ( 2010 ). OpenUrl 56. Best , R. B. & Hummer , G. Optimized Molecular Dynamics Force Fields Applied to the Helix-Coil Transition of Polypeptides . J. Phys. Chem. B 113 , 9004 – 9015 ( 2009 ). OpenUrl CrossRef PubMed 57. ↵ Hornak , V. et al. Comparison of multiple Amber force fields and development of improved protein backbone parameters . Proteins Struct. Funct. Bioinforma . 65 , 712 – 725 ( 2006 ). OpenUrl CrossRef 58. ↵ Aldeghi , M. , Gapsys , V. & de Groot , B. L. Accurate Estimation of Ligand Binding Affinity Changes upon Protein Mutation . ACS Cent. Sci . 4 , 1708 – 1718 ( 2018 ). OpenUrl PubMed 59. ↵ Jorgensen , W. L. , Chandrasekhar , J. , Madura , J. D. , Impey , R. W. & Klein , M. L. Comparison of simple potential functions for simulating liquid water . J. Chem. Phys . 79 , 926 – 935 ( 1983 ). OpenUrl CrossRef PubMed Web of Science 60. ↵ Bussi , G. , Donadio , D. & Parrinello , M. Canonical sampling through velocity rescaling . J. Chem. Phys . 126 , ( 2007 ). 61. ↵ Parrinello , M. & Rahman , A. Polymorphic transitions in single crystals: A new molecular dynamics method . J. Appl. Phys . 52 , 7182 – 7190 ( 1981 ). OpenUrl CrossRef PubMed Web of Science 62. ↵ Bennett , C. H. Efficient estimation of free energy differences from Monte Carlo data . J. Comput. Phys . 22 , 245 – 268 ( 1976 ). OpenUrl CrossRef 63. ↵ Shirts , M. R. & Pande , V. S. Comparison of efficiency and bias of free energies computed by exponential averaging, the Bennett acceptance ratio, and thermodynamic integration . J. Chem. Phys . 122 , ( 2005 ). 64. ↵ Abel , R. , Wang , L. , Harder , E. D. , Berne , B. J. & Friesner , R. A. Advancing Drug Discovery through Enhanced Free Energy Calculations . Acc. Chem. Res . 50 , 1625 – 1632 ( 2017 ). OpenUrl CrossRef PubMed 65. ↵ Abel , R. , Wang , L. , Mobley , D. L. & Friesner , R. A. A Critical Review of Validation, Blind Testing, and Real-World Use of Alchemical Protein-Ligand Binding Free Energy Calculations . Curr. Top. Med. Chem . 17 , 2577 – 2585 ( 2017 ). OpenUrl CrossRef PubMed 66. ↵ Steinbrecher , T. et al. Predicting the Effect of Amino Acid Single-Point Mutations on Protein Stability—Large-Scale Validation of MD-Based Relative Free Energy Calculations . J. Mol. Biol . 429 , 948 – 963 ( 2017 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted October 20, 2025. Download PDF Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems Shivani Gupta , Qinfang Sun , Ronald M. Levy bioRxiv 2025.10.19.683321; doi: https://doi.org/10.1101/2025.10.19.683321 Share This Article: Copy Citation Tools Benchmarking Free Energy Calculations: Analysis of Single and Double Mutations Across Two Simulation Software Platforms for Two Protein Systems Shivani Gupta , Qinfang Sun , Ronald M. Levy bioRxiv 2025.10.19.683321; doi: https://doi.org/10.1101/2025.10.19.683321 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Biophysics Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41911) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13371) Ecology (19887) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40364) Molecular Biology (17163) Neuroscience (88537) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15129) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9817) Zoology (2269)