How to Sample Dozens of Substitutions per Site with λ Dynamics.

Alchemical free energy methods are useful in computer-aided drug design and computational protein design because they provide rigorous statistical mechanics-based estimates of free energy differences from molecular dynamics simulations. λ dynamics is a free energy method with the ability to characterize combinatorial chemical spaces spanning thousands of related systems within a single simulation, which gives it a distinct advantage over other alchemical free energy methods that are mostly limited to pairwise comparisons. Recently developed methods have improved the scalability of λ dynamics to perturbations at many sites; however, the size of chemical space that can be explored at each individual site has previously been limited to fewer than ten substituents. As the number of substituents increases, the volume of alchemical space corresponding to nonphysical alchemical intermediates grows exponentially relative to the size corresponding to the physical states of interest. Beyond nine substituents, λ dynamics simulations become lost in an alchemical morass of intermediate states. In this work, we introduce new biasing potentials that circumvent excessive sampling of intermediate states by favoring sampling of physical end points relative to alchemical intermediates. Additionally, we present a more scalable adaptive landscape flattening algorithm for these larger alchemical spaces. Finally, we show that this potential enables more efficient sampling in both protein and drug design test systems with up to 24 substituents per site, enabling, for the first time, simultaneous simulation of all 20 amino acids.

A λ-Dynamics Investigation of Insulin Wakayama and Other A3 Variant Binding Affinities to the Insulin Receptor.

Insulin Wakayama is a clinical insulin variant where a conserved valine at the third residue on insulin's A chain (ValA3) is replaced with a leucine (LeuA3), weakening insulin receptor (IR) binding by 140-500-fold. This severe impact on binding from a subtle modification has posed an intriguing problem for decades. Although experimental investigations of natural and unnatural A3 mutations have highlighted the sensitivity of insulin-IR binding at this site, atomistic explanations of these binding trends have remained elusive. We investigate this problem computationally using λ-dynamics free energy calculations to model structural changes in response to perturbations of the ValA3 side chain and to calculate associated relative changes in binding free energy (ΔΔGbind). The Wakayama LeuA3 mutation and seven other A3 substitutions were studied in this work. The calculated ΔΔGbind results showed high agreement compared to experimental binding potencies with a Pearson correlation of 0.88 and a mean unsigned error of 0.68 kcal/mol. Extensive structural analyses of λ-dynamics trajectories revealed that critical interactions were disrupted between insulin and the insulin receptor as a result of the A3 mutations. This investigation also quantifies the effect that adding an A3 Cδ atom or losing an A3 Cγ atom has on insulin's binding affinity to the IR. Thus, λ-dynamics was able to successfully model the effects of mutations to insulin's A3 side chain on its protein-protein interactions with the IR and shed new light on a decades-old mystery: the exquisite sensitivity of hormone-receptor binding to a subtle modification of an invariant insulin residue.

A λ-dynamics investigation of insulin Wakayama and other A3 variant binding affinities to the insulin receptor.

Insulin Wakayama is a clinical insulin variant where a conserved valine at the third residue on insulin's A chain (ValA3) is replaced with a leucine (LeuA3), impairing insulin receptor (IR) binding by 140-500 fold. This severe impact on binding from such a subtle modification has posed an intriguing problem for decades. Although experimental investigations of natural and unnatural A3 mutations have highlighted the sensitivity of insulin-IR binding to minor changes at this site, an atomistic explanation of these binding trends has remained elusive. We investigate this problem computationally using λ-dynamics free energy calculations to model structural changes in response to perturbations of the ValA3 side chain and to calculate associated relative changes in binding free energy (ΔΔGbind). The Wakayama LeuA3 mutation and seven other A3 substitutions were studied in this work. The calculated ΔΔGbind results showed high agreement compared to experimental binding potencies with a Pearson correlation of 0.88 and a mean unsigned error of 0.68 kcal/mol. Extensive structural analyses of λ-dynamics trajectories revealed that critical interactions were disrupted between insulin and the insulin receptor as a result of the A3 mutations. This investigation also quantifies the effect that adding an A3 Cδ atom or losing an A3 Cγ atom has on insulin's binding affinity to the IR. Thus, λ-dynamics was able to successfully model the effects of subtle modifications to insulin's A3 side chain on its protein-protein interactions with the IR and shed new light on a decades-old mystery: the exquisite sensitivity of hormone-receptor binding to a subtle modification of an invariant insulin residue.

In silico λ-dynamics predicts protein binding specificities to modified RNAs.

RNA modifications shape gene expression through a smorgasbord of chemical changes to canonical RNA bases. Although numbering in the hundreds, only a few RNA modifications are well characterized, in part due to the absence of methods to identify modification sites. Antibodies remain a common tool to identify modified RNA and infer modification sites through straightforward applications. However, specificity issues can result in off-target binding and confound conclusions. This work utilizes in silico λ-dynamics to efficiently estimate binding free energy differences of modification-targeting antibodies between a variety of naturally occurring RNA modifications. Crystal structures of inosine and N6-methyladenosine (m6A) targeting antibodies bound to their modified ribonucleosides were determined and served as structural starting points. λ-Dynamics was utilized to predict RNA modifications that permit or inhibit binding to these antibodies. In vitro RNA-antibody binding assays supported the accuracy of these in silico results. High agreement between experimental and computed binding propensities demonstrated that λ-dynamics can serve as a predictive screen for antibody specificity against libraries of RNA modifications. More importantly, this strategy is an innovative way to elucidate how hundreds of known RNA modifications interact with biological molecules without the limitations imposed by in vitro or in vivo methodologies.

Fast free energy estimates from λ-dynamics with bias-updated Gibbs sampling.

Relative binding free energy calculations have become an integral computational tool for lead optimization in structure-based drug design. Classical alchemical methods, including free energy perturbation or thermodynamic integration, compute relative free energy differences by transforming one molecule into another. However, these methods have high operational costs due to the need to perform many pairwise perturbations independently. To reduce costs and accelerate molecular design workflows, we present a method called λ-dynamics with bias-updated Gibbs sampling. This method uses dynamic biases to continuously sample between multiple ligand analogues collectively within a single simulation. We show that many relative binding free energies can be determined quickly with this approach without compromising accuracy. For five benchmark systems, agreement to experiment is high, with root mean square errors near or below 1.0 kcal mol-1. Free energy results are consistent with other computational approaches and within statistical noise of both methods (0.4 kcal mol-1 or less). Notably, large efficiency gains over thermodynamic integration of 18-66-fold for small perturbations and 100-200-fold for whole aromatic ring substitutions are observed. The rapid determination of relative binding free energies will enable larger chemical spaces to be more readily explored and structure-based drug design to be accelerated.

Chemically induced partial unfolding of the multifunctional Apurinic/apyrimidinic endonuclease 1.

Targeting of the multifunctional enzyme apurinic/apyrimidinic endonuclease I/redox factor 1 (APE1) has produced small molecule inhibitors of both its endonuclease and redox activities. While one of the small molecules, the redox inhibitor APX3330, completed a Phase I clinical trial for solid tumors and a Phase II clinical trial for Diabetic Retinopathy/Diabetic Macular Edema, the mechanism of action for this drug has yet to be fully understood. Here, we demonstrate through HSQC NMR studies that APX3330 induces chemical shift perturbations (CSPs) of both surface and internal residues in a concentration-dependent manner, with a cluster of surface residues defining a small pocket on the opposite face from the endonuclease active site of APE1. Furthermore, APX3330 induces partial unfolding of APE1 as evidenced by a time-dependent loss of chemical shifts for approximately 35% of the residues within APE1 in the HSQC NMR spectrum. Notably, regions that are partially unfolded include adjacent strands within one of two beta sheets that comprise the core of APE1. One of the strands comprises residues near the N-terminal region and a second strand is contributed by the C-terminal region of APE1, which serves as a mitochondrial targeting sequence. These terminal regions converge within the pocket defined by the CSPs. In the presence of a duplex DNA substrate mimic, removal of excess APX3330 resulted in refolding of APE1. Our results are consistent with a reversible mechanism of partial unfolding of APE1 induced by the small molecule inhibitor, APX3330, defining a novel mechanism of inhibition.

Identification of nonhistone substrates of the lysine methyltransferase PRDM9.

Lysine methylation is a dynamic, posttranslational mark that regulates the function of histone and nonhistone proteins. Many of the enzymes that mediate lysine methylation, known as lysine methyltransferases (KMTs), were originally identified to modify histone proteins but have also been discovered to methylate nonhistone proteins. In this work, we investigate the substrate selectivity of the KMT PRDM9 to identify both potential histone and nonhistone substrates. Though normally expressed in germ cells, PRDM9 is significantly upregulated across many cancer types. The methyltransferase activity of PRDM9 is essential for double-strand break formation during meiotic recombination. PRDM9 has been reported to methylate histone H3 at lysine residues 4 and 36; however, PRDM9 KMT activity had not previously been evaluated on nonhistone proteins. Using lysine-oriented peptide libraries to screen potential substrates of PRDM9, we determined that PRDM9 preferentially methylates peptide sequences not found in any histone protein. We confirmed PRDM9 selectivity through in vitro KMT reactions using peptides with substitutions at critical positions. A multisite λ-dynamics computational analysis provided a structural rationale for the observed PRDM9 selectivity. The substrate selectivity profile was then used to identify putative nonhistone substrates, which were tested by peptide spot array, and a subset was further validated at the protein level by in vitro KMT assays on recombinant proteins. Finally, one of the nonhistone substrates, CTNNBL1, was found to be methylated by PRDM9 in cells.

Addressing Intersite Coupling Unlocks Large Combinatorial Chemical Spaces for Alchemical Free Energy Methods.

Alchemical free energy methods are playing a growing role in molecular design, both for computer-aided drug design of small molecules and for computational protein design. Multisite λ dynamics (MSλD) is a uniquely scalable alchemical free energy method that enables more efficient exploration of combinatorial alchemical spaces encountered in molecular design, but simulations have typically been limited to a few hundred ligands or sequences. Here, we focus on coupling between sites to enable scaling to larger alchemical spaces. We first discuss updates to the biasing potentials that facilitate MSλD sampling to include coupling terms and show that this can provide more thorough sampling of alchemical states. We then harness coupling between sites by developing a new free energy estimator based on the Potts models underlying direct coupling analysis, a method for predicting contacts from sequence coevolution, and find it yields more accurate free energies than previous estimators. The sampling requirements of the Potts model estimator scale with the square of the number of sites, a substantial improvement over the exponential scaling of the standard estimator. This opens up exploration of much larger alchemical spaces with MSλD for molecular design.

Optimizing Multisite λ-Dynamics Throughput with Charge Renormalization.

With the ability to sample combinations of alchemical perturbations at multiple sites off a small molecule core, multisite λ-dynamics (MSλD) has become an attractive alternative to conventional alchemical free energy methods for exploring large combinatorial chemical spaces. However, current software implementations dictate that combinatorial sampling with MSλD must be performed with a multiple topology model (MTM), which is nontrivial to create by hand, especially for a series of ligand analogues which may have diverse functional groups attached. This work introduces an automated workflow, referred to as msld_py_prep, to assist in the creation of a MTM for use with MSλD. One approach for partitioning partial atomic charges between ligands to create a MTM, called charge renormalization, is also presented and rigorously evaluated. We find that msld_py_prep greatly accelerates the preparation of MSλD ready-to-use files and that charge renormalization can provide a successful approach for MTM generation, as long as bookending calculations are applied to correct small differences introduced by charge renormalization. Charge renormalization also facilitates the use of many different force field parameters with MSλD, broadening the applicability of MSλD for computer-aided drug design.

Improved Cell-Potent and Selective Peptidomimetic Inhibitors of Protein N-Terminal Methyltransferase 1.

Protein N-terminal methyltransferase 1 (NTMT1) recognizes a unique N-terminal X-P-K/R motif (X represents any amino acid other than D/E) and transfers 1-3 methyl groups to the N-terminal region of its substrates. Guided by the co-crystal structures of NTMT1 in complex with the previously reported peptidomimetic inhibitor DC113, we designed and synthesized a series of new peptidomimetic inhibitors. Through a focused optimization of DC113, we discovered a new cell-potent peptidomimetic inhibitor GD562 (IC50 = 0.93 ± 0.04 µM). GD562 exhibited improved inhibition of the cellular N-terminal methylation levels of both the regulator of chromosome condensation 1 and the oncoprotein SET with an IC50 value of ~50 µM in human colorectal cancer HCT116 cells. Notably, the inhibitory activity of GD562 for the SET protein increased over 6-fold compared with the previously reported cell-potent inhibitor DC541. Furthermore, GD562 also exhibited over 100-fold selectivity for NTMT1 against several other methyltransferases. Thus, this study provides a valuable probe to investigate the biological functions of NTMT1.

Exploring Unconventional SAM Analogues To Build Cell-Potent Bisubstrate Inhibitors for Nicotinamide N-Methyltransferase.

Nicotinamide N-methyltransferase (NNMT) methylates nicotinamide and has been associated with various diseases. Herein, we report the first cell-potent NNMT bisubstrate inhibitor II399, demonstrating a Ki of 5.9 nM in a biochemical assay and a cellular IC50 value of 1.9 µM. The inhibition mechanism and cocrystal structure confirmed II399 engages both the substrate and cofactor binding pockets. Computational modeling and binding data reveal a balancing act between enthalpic and entropic components that lead to II399's low nM binding affinity. Notably, II399 is 1 000-fold more selective for NNMT than closely related methyltransferases. We expect that II399 would serve as a valuable probe to elucidate NNMT biology. Furthermore, this strategy provides the first case of introducing unconventional SAM mimics, which can be adopted to develop cell-potent inhibitors for other SAM-dependent methyltransferases.

Generalizing the Discrete Gibbs Sampler-Based λ-Dynamics Approach for Multisite Sampling of Many Ligands.

In this work, the discrete λ variant of the Gibbs sampler-based λ-dynamics (d-GSλD) method is developed to enable multiple functional group perturbations to be investigated at one or more sites of substitution off a common ligand core. The theoretical framework and special considerations for constructing discrete λ states for multisite d-GSλD are presented. The precision and accuracy of the d-GSλD method is evaluated with three test cases of increasing complexity. Specifically, methyl → methyl symmetric perturbations in water, 1,4-benzene hydration free energies and protein-ligand binding affinities for an example HIV-1 reverse transcriptase inhibitor series are computed with d-GSλD. Complementary MSλD calculations were also performed to compare with d-GSλD's performance. Excellent agreement between d-GSλD and MSλD is observed, with mean unsigned errors of 0.12 and 0.22 kcal/mol for computed hydration and binding free energy test cases, respectively. Good agreement with experiment is also observed, with errors of 0.5-0.7 kcal/mol. These findings support the applicability of the d-GSλD free energy method for a variety of molecular design problems, including structure-based drug design. Finally, a discussion of d-GSλD versus MSλD approaches is presented to compare and contrast features of both methods.

Gas-Phase Ion/Ion Chemistry for Structurally Sensitive Probes of Gaseous Protein Ion Structure: Electrostatic and Electrostatic to Covalent Cross-Linking.

Intramolecular interactions within a protein are key in maintaining protein tertiary structure and understanding how proteins function. Ion mobility-mass spectrometry (IM-MS) has become a widely used approach in structural biology since it provides rapid measurements of collision cross sections (CCS), which inform on the gas-phase conformation of the biomolecule under study. Gas-phase ion/ion reactions target amino acid residues with specific chemical properties and the modified sites can be identified by MS. In this study, electrostatically reactive, gas-phase ion/ion chemistry and IM-MS are combined to characterize the structural changes between ubiquitin electrosprayed from aqueous and denaturing conditions. The electrostatic attachment of sulfo-NHS acetate to ubiquitin via ion/ion reactions and fragmentation by electron-capture dissociation (ECD) provide the identification of the most accessible protonated sites within ubiquitin as the sulfonate group forms an electrostatic complex with accessible protonated side chains. The protonated sites identified by ECD from the different solution conditions are distinct and, in some cases, reflect the disruption of interactions such as salt bridges that maintain the native protein structure. This agrees with previously published literature demonstrating that a high methanol concentration at low pH causes the structure of ubiquitin to change from a native (N) state to a more elongated A state. Results using gas-phase, electrostatic cross-linking reagents also point to similar structural changes and further confirm the role of methanol and acid in favoring a more unfolded conformation. Since cross-linking reagents have a distance constraint for the two reactive sites, the data is valuable in guiding computational structures generated by molecular dynamics. The research presented here describes a promising strategy that can detect subtle changes in the local environment of targeted amino acid residues to inform on changes in the overall protein structure.

Mutant thermal proteome profiling for characterization of missense protein variants and their associated phenotypes within the proteome.

Temperature-sensitive (TS) missense mutants have been foundational for characterization of essential gene function. However, an unbiased approach for analysis of biochemical and biophysical changes in TS missense mutants within the context of their functional proteomes is lacking. We applied MS-based thermal proteome profiling (TPP) to investigate the proteome-wide effects of missense mutations in an application that we refer to as mutant thermal proteome profiling (mTPP). This study characterized global impacts of temperature sensitivity-inducing missense mutations in two different subunits of the 26S proteasome. The majority of alterations identified by RNA-Seq and global proteomics were similar between the mutants, which could suggest that a similar functional disruption is occurring in both missense variants. Results from mTPP, however, provide unique insights into the mechanisms that contribute to the TS phenotype in each mutant, revealing distinct changes that were not obtained using only steady-state transcriptome and proteome analyses. Computationally, multisite λ-dynamics simulations add clear support for mTPP experimental findings. This work shows that mTPP is a precise approach to measure changes in missense mutant-containing proteomes without the requirement for large amounts of starting material, specific antibodies against proteins of interest, and/or genetic manipulation of the biological system. Although experiments were performed under permissive conditions, mTPP provided insights into the underlying protein stability changes that cause dramatic cellular phenotypes observed at nonpermissive temperatures. Overall, mTPP provides unique mechanistic insights into missense mutation dysfunction and connection of genotype to phenotype in a rapid, nonbiased fashion.

Exploring pH Dependent Host/Guest Binding Affinities.

When the electrostatic environment surrounding binding partners changes between unbound and bound states, the net uptake or release of a proton is possible by either binding partner. This process is pH-dependent in that the free energy required to uptake or release the proton varies with pH. This pH-dependence is typically not considered in conventional free energy methods where the use of fixed protonation states is the norm. In the present paper, we apply a simple two-step approach to calculate the pH-dependent binding free energy of a model cucubit[7]uril host/guest system. By use of λ-dynamics with an enhanced sampling protocol, adaptive landscape flattening, pKa shifts and reference binding free energies upon complexation were determined. This information enables the construction of pH-dependent binding profiles that accurately capture the pKa shifts and reproduce binding free energies at the different pH conditions that were observed experimentally. Our calculations illustrate a general framework for computing pH-dependent binding free energies but also point to some issues in modeling the molecular charge distributions within this series of molecules with CGenFF. However, by introducing some minor charge modifications to the CGenFF force field, we saw significant improvement in accuracy of the calculated pKa shifts.

Accelerated CDOCKER with GPUs, Parallel Simulated Annealing, and Fast Fourier Transforms.

Fast Fourier transform (FFT)-based protein ligand docking together with parallel simulated annealing for both rigid and flexible receptor docking are implemented on graphical processing unit (GPU) accelerated platforms to significantly enhance the throughput of the CDOCKER and flexible CDOCKER - the docking algorithms in the CHARMM program for biomolecule modeling. The FFT-based approach for docking, first applied in protein-protein docking to efficiently search for the binding position and orientation of proteins, is adapted here to search ligand translational and rotational spaces given a ligand conformation in protein-ligand docking. Running on GPUs, our implementation of FFT docking in CDOCKER achieves a 15 000 fold speedup in the ligand translational and rotational space search in protein-ligand docking problems. With this significant speedup it becomes practical to exhaustively search ligand translational and rotational space when docking a rigid ligand into a protein receptor. We demonstrate in this paper that this provides an efficient way to calculate an upper bound for docking accuracy in the assessment of scoring functions for protein-ligand docking, which can be useful for improving scoring functions. The parallel molecular dynamics (MD) simulated annealing, also running on GPUs, aims to accelerate the search algorithm in CDOCKER by running MD simulated annealing in parallel on GPUs. When utilized as part of the general CDOCKER docking protocol, acceleration in excess of 20 times is achieved. With this acceleration, we demonstrate that the performance of CDOCKER for redocking is significantly improved compared with three other popular protein-ligand docking programs on two widely used protein ligand complex data sets: the Astex diverse set and the SB2012 test set. The flexible CDOCKER is similarly improved by the parallel MD simulated annealing on GPUs. Based on the results presented here, we suggest that the accelerated CDOCKER platform provides a highly competitive docking engine for both rigid-receptor and flexible-receptor docking studies and will further facilitate continued improvement in the physics-based scoring function employed in CDOCKER docking studies.

Absolute Free Energy of Binding Calculations for Macrophage Migration Inhibitory Factor in Complex with a Druglike Inhibitor.

Calculation of the absolute free energy of binding (ΔGbind) for a complex in solution is challenging owing to the need for adequate configurational sampling and an accurate energetic description, typically with a force field (FF). In this study, Monte Carlo (MC) simulations with improved side-chain and backbone sampling are used to assess ΔGbind for the complex of a druglike inhibitor (MIF180) with the protein macrophage migration inhibitory factor (MIF) using free energy perturbation (FEP) calculations. For comparison, molecular dynamics (MD) simulations were employed as an alternative sampling method for the same system. With the OPLS-AA/M FF and CM5 atomic charges for the inhibitor, the ΔGbind results from the MC/FEP and MD/FEP simulations, -8.80 ± 0.74 and -8.46 ± 0.85 kcal/mol, agree well with each other and with the experimental value of -8.98 ± 0.28 kcal/mol. The convergence of the results and analysis of the trajectories indicate that sufficient sampling was achieved for both approaches. Repeating the MD/FEP calculations using current versions of the CHARMM and AMBER FFs led to a 6 kcal/mol range of computed ΔGbind. These results show that calculation of accurate ΔGbind for large ligands is both feasible and numerically equivalent, within error limits, using either methodology.

Overcoming Challenging Substituent Perturbations with Multisite λ-Dynamics: A Case Study Targeting ß-Secretase 1.

Alchemical free energy calculations have made a dramatic impact upon the field of structure-based drug design by allowing functional group modifications to be explored computationally prior to experimental synthesis and assay evaluation, thereby informing and directing synthetic strategies. In furthering the advancement of this area, a series of 21 ß-secretase 1 (BACE1) inhibitors developed by Janssen Pharmaceuticals were examined to evaluate the ability to explore large substituent perturbations, some of which contain scaffold modifications, with multisite λ-dynamics (MSλD), an innovative alchemical free energy framework. Our findings indicate that MSλD is able to efficiently explore all structurally diverse ligand end-states simultaneously within a single MD simulation with a high degree of precision and with reduced computational costs compared to the widely used approach TI/MBAR. Furthermore, computational predictions were shown to be accurate to within 0.5-0.8 kcal/mol when CM1A partial atomic charges were combined with CHARMM or OPLS-AA-based force fields, demonstrating that MSλD is force field independent and a viable alternative to FEP or TI approaches for drug design.

Fast Solver for Large Scale Multistate Bennett Acceptance Ratio Equations.

The multistate Bennett acceptance ratio method (MBAR) and unbinned weighted histogram analysis method (UWHAM) are widely employed approaches to calculate relative free energies of multiple thermodynamic states that gain statistical precision by employing free energy contributions from configurations sampled at each of the simulated λ states. With the increasing availability of high throughput computing resources, a large number of configurations can be sampled from hundreds or even thousands of states. Combining sampled configurations from all states to calculate relative free energies requires the iterative solution of large scale MBAR/UWHAM equations. In the current work, we describe the development of a fast solver to iteratively solve these large scale MBAR/UWHAM equations utilizing our previous findings that the MBAR/UWHAM equations can be derived as a Rao-Blackwell estimator. The solver is implemented and distributed as a Python module called FastMBAR. Our benchmark results show that FastMBAR is more than 2 times faster than the widely used solver pymbar, when it runs on a central processing unit (CPU) and more than 100 times faster than pymbar when it runs on a graphical processing unit (GPU). The significant speedup achieved by FastMBAR running on a GPU is useful not only for solving large scale MBAR/UWHAM equations but also for estimating uncertainty of calculated free energies using bootstrapping where the MBAR/UWHAM equations need to be solved multiple times.

Approaching protein design with multisite λ dynamics: Accurate and scalable mutational folding free energies in T4 lysozyme.

The estimation of changes in free energy upon mutation is central to the problem of protein design. Modern protein design methods have had remarkable success over a wide range of design targets, but are reaching their limits in ligand binding and enzyme design due to insufficient accuracy in mutational free energies. Alchemical free energy calculations have the potential to supplement modern design methods through more accurate molecular dynamics based prediction of free energy changes, but suffer from high computational cost. Multisite λ dynamics (MSλD) is a particularly efficient and scalable free energy method with potential to explore combinatorially large sequence spaces inaccessible with other free energy methods. This work aims to quantify the accuracy of MSλD and demonstrate its scalability. We apply MSλD to the classic problem of calculating folding free energies in T4 lysozyme, a system with a wealth of experimental measurements. Single site mutants considering 32 mutations show remarkable agreement with experiment with a Pearson correlation of 0.914 and mean unsigned error of 1.19 kcal/mol. Multisite mutants in systems with up to five concurrent mutations spanning 240 different sequences show comparable agreement with experiment. These results demonstrate the promise of MSλD in exploring large sequence spaces for protein design.