Computational Modeling of (De)-Solvation Effects and Protein Flexibility in Protein-Ligand Binding using Molecular Dynamics Simulations

Water is a crucial participant in virtually all cellular functions. Evidently, water molecules in the binding site contribute significantly to the strength of intermolecular interactions in the aqueous phase by mediating protein-ligand interactions, solvating and de-solvating both ligand and protein upon protein-ligand dissociation and association. Recently many published studies use water distributions in the binding site to retrospectively explain and rationalize unexpected trends in structure-activity relationships (SAR). However, traditional approaches cannot quantitatively predict the thermodynamic properties of water molecules in the binding sites and its associated contribution to the binding free energy of a ligand. We have developed and validated a computational method named WATsite to exploit high-resolution solvation maps and thermodynamic profiles to elucidate the water molecules’ potential contribution to protein-ligand and protein-protein binding. We have also demonstrated the utility of the computational method WATsite to help direct medicinal chemistry efforts by using explicit water de-solvation. In addition, protein conformational change is typically involved in the ligand-binding process which may completely change the position and thermodynamic properties of the water molecules in the binding site before or upon ligand binding. We have shown the interplay between protein flexibility and solvent reorganization, and we provide a quantitative estimation of the influence of protein flexibility on desolvation free energy and, therefore, protein-ligand binding. Different ligands binding to the same target protein can induce different conformational adaptations. In order to apply WATsite to an ensemble of different protein conformations, a more efficient implementation of WATsite based on GPU-acceleration and system truncation has been developed. Lastly, by extending the simulation protocol from pure water to mixed water-organic probes simulations, accurate modeling of halogen atom-protein interactions has been achieved

ATOMISTIC SIMULATIONS OF INTRINSICALLY DISORDERED PROTEIN FOLDING AND DYNAMICS

Intrinsically disordered proteins (IDPs) are crucial in biology and human diseases, necessitating a comprehensive understanding of their structure, dynamics, and interactions. Atomistic simulations have emerged as a key tool for unraveling the molecular intricacies and establishing mechanistic insights into how these proteins facilitate diverse biological functions. However, achieving accurate simulations requires both an appropriate protein force field capable of describing the energy landscape of functionally relevant IDP conformations and sufficient conformational sampling to capture the free energy landscape of IDP dynamics. These factors are fundamental in comprehending potential IDP structures, dynamics, and interactions. I first conducted explicit solvent simulations to assess the performance of two state-of-the-art protein force fields, namely CHARMM36m and a99SB-disp, in capturing the stability of small protein-protein interactions. To evaluate their accuracy, I selected a set of 46 amino acid backbone and side chain pairs with representative configurations and computed the free energy profiles of their interactions. The results demonstrated that CHARMM36m consistently predicted stronger protein-protein interactions compared to a99SB-disp. Notably, the most significant overestimation in CHARMM36m occurred in charged pairs involving Arg and Glu side chains, with an overestimation of up to 2.9 kcal/mol. Through free energy decomposition analysis, I determined that these overestimations were primarily driven by protein-water electrostatic interactions rather than van der Waals (vdW) interactions. Consequently, these findings suggest that careful rebalancing of electrostatic interactions should be considered in the further optimization of protein force fields. In order to enhance the conformational sampling of IDPs, I developed an integrated approach that combines an improved implicit solvent model called Generalized Born with molecular volume and solvent accessible surface area (GBMV2/SA) with a multiscale enhanced sampling (MSES) technique. To make this approach more efficient, I implemented it as a standalone OpenMM plugin on Graphics Processing Units (GPUs). The results demonstrated that the GPU-GBMV2/SA model achieved numerical equivalence to the original CPU-GBMV2/SA models, while providing a remarkable ~60x speedup on a single NVIDIA TITAN X (Pascal) graphics card for molecular dynamic simulations of both folded and unstructured proteins. This significant acceleration greatly facilitated the application of the approach in biomolecular simulations. In addition, I conducted an evaluation of the reliability of GBMV2/SA models in simulating both folded and unfolded proteins. The results revealed that the GBMV2/SA model accurately describes small proteins, but its applicability is limited when it comes to larger proteins such as KID and p53-TAD proteins. This limitation can be attributed to the absence of long-range solute-solvent dispersion interactions in the model. To address this issue, I introduced a comprehensive treatment of nonpolar solvation free energy called GBMV2/NP model. Unfortunately, the GBMV2/NP model exhibited a destabilizing effect on well-folded proteins, particularly larger ones, due to an inaccurate representation of the repulsive solvent accessible surface area (SASA) model caused by the utilization of unphysical van der Waals volume. This observation highlights the need for further improvements in accurately describing the nonpolar term in the model

Investigating substitutions in antibody–antigen complexes using molecular dynamics: a case study with broad-spectrum, Influenza A antibodies

In studying the binding of host antibodies to the surface antigens of pathogens, the structural and functional characterization of antibody–antigen complexes by X-ray crystallography and binding assay is important. However, the characterization requires experiments that are typically time consuming and expensive: thus, many antibody–antigen complexes are under-characterized. For vaccine development and disease surveillance, it is often vital to assess the impact of amino acid substitutions on antibody binding. For example, are there antibody substitutions capable of improving binding without a loss of breadth, or antigen substitutions that lead to antigenic escape? The questions cannot be answered reliably from sequence variation alone, exhaustive substitution assays are usually impractical, and alanine scans provide at best an incomplete identification of the critical residue–residue interactions. Here, we show that, given an initial structure of an antibody bound to an antigen, molecular dynamics simulations using the energy method molecular mechanics with Generalized Born surface area (MM/GBSA) can model the impact of single amino acid substitutions on antibody–antigen binding energy. We apply the technique to three broad-spectrum antibodies to influenza A hemagglutinin and examine both previously characterized and novel variant strains observed in the human population that may give rise to antigenic escape. We find that in some cases the impact of a substitution is local, while in others it causes a reorientation of the antibody with wide-ranging impact on residue–residue interactions: this explains, in part, why the change in chemical properties of a residue can be, on its own, a poor predictor of overall change in binding energy. Our estimates are in good agreement with experimental results—indeed, they approximate the degree of agreement between different experimental techniques. Simulations were performed on commodity computer hardware; hence, this approach has the potential to be widely adopted by those undertaking infectious disease research. Novel aspects of this research include the application of MM/GBSA to investigate binding between broadly binding antibodies and a viral glycoprotein; the development of an approach for visualizing substrate–ligand interactions; and the use of experimental assay data to rescale our predictions, allowing us to make inferences about absolute, as well as relative, changes in binding energy

Extending fragment-based free energy calculations with library Monte Carlo simulation: Annealing in interaction space

Pre-calculated libraries of molecular fragment configurations have previously been used as a basis for both equilibrium sampling (via "library-based Monte Carlo") and for obtaining absolute free energies using a polymer-growth formalism. Here, we combine the two approaches to extend the size of systems for which free energies can be calculated. We study a series of all-atom poly-alanine systems in a simple dielectric "solvent" and find that precise free energies can be obtained rapidly. For instance, for 12 residues, less than an hour of single-processor is required. The combined approach is formally equivalent to the "annealed importance sampling" algorithm; instead of annealing by decreasing temperature, however, interactions among fragments are gradually added as the molecule is "grown." We discuss implications for future binding affinity calculations in which a ligand is grown into a binding site

Exploring Mechanical Properties and Configurational Energetics of Toxbox Using Molecular Dynamics Simulation

All-atom Molecular Dynamics Simulations (MDS) have been performed to obtain 5 ns trajectory of the solvated, neutralized, and equilibrated toxbox system in NPT ensemble at 300 K. This trajectory data has been used to calculate the configurational entropy of toxbox by employing a quantum mechanical approach. The method is based on evaluating determinant of the covariance matrix, built from generalized coordinates of all atoms for each frame. The upper limit to the configurational entropy of toxbox has been calculated to be 30,030 J/mol-K. A preliminary investigation has been conducted to study the effects of sequence-dependent DNA conformation (DNA Crookedness) on the mechanical properties of toxbox by implementing constant force MDS. Results of this research may serve as the reference for studying ToxT – DNA interactions

Development and analysis of Tinker-OpenMM as a GPU-based free energy perturbation engine

The utilization of computational technologies for the lead optimization process is one of the biggest challenges in the computational chemistry field. In this dissertation, I describe the addition of GPU-based absolute and relative free energy calculation methods using polarizable force field AMOEBA to Tinker-OpenMM. I then proceed to test the capabilities of this platform by studying the binding free energy and binding structures of derivatives of the MELK inhibitor IN17. Also, I present the implementation of virial-based pressure control to the Tinker-OpenMM platform that is needed for performing isobaric simulations.Cellular and Molecular Biolog

Investigation of Structural Dynamics of Enzymes and Protonation States of Substrates Using Computational Tools.

This review discusses the use of molecular modeling tools, together with existing experimental findings, to provide a complete atomic-level description of enzyme dynamics and function. We focus on functionally relevant conformational dynamics of enzymes and the protonation states of substrates. The conformational fluctuations of enzymes usually play a crucial role in substrate recognition and catalysis. Protein dynamics can be altered by a tiny change in a molecular system such as different protonation states of various intermediates or by a significant perturbation such as a ligand association. Here we review recent advances in applying atomistic molecular dynamics (MD) simulations to investigate allosteric and network regulation of tryptophan synthase (TRPS) and protonation states of its intermediates and catalysis. In addition, we review studies using quantum mechanics/molecular mechanics (QM/MM) methods to investigate the protonation states of catalytic residues of β-Ketoacyl ACP synthase I (KasA). We also discuss modeling of large-scale protein motions for HIV-1 protease with coarse-grained Brownian dynamics (BD) simulations

AutoDock-SS: AutoDock for Multiconformational Ligand-based Virtual Screening

Ligand-based virtual screening (LBVS) can be pivotal for identifying potential drug leads, especially when the target protein’s structure is unknown. However, current LBVS methods are limited in their ability to consider the ligand conformational flexibility. This study presents AutoDock-SS (Similarity Searching), which adapts protein–ligand docking for use in LBVS. AutoDock-SS integrates novel ligand-based grid maps and AutoDock-GPU into a novel three-dimensional LBVS workflow. Unlike other approaches based on pregenerated conformer libraries, AutoDock-SS’s built-in conformational search optimizes conformations dynamically based on the reference ligand, thus providing a more accurate representation of relevant ligand conformations. AutoDock-SS supports two modes: single and multiple ligand queries, allowing for the seamless consideration of multiple reference ligands. When tested on the Directory of Useful Decoys─Enhanced (DUD-E) data set, AutoDock-SS surpassed alternative 3D LBVS methods, achieving a mean AUROC of 0.775 and an EF1% of 25.72 in single-reference mode. The multireference mode, evaluated on the augmented DUD-E+ data set, demonstrated superior accuracy with a mean AUROC of 0.843 and an EF1% of 34.59. This enhanced performance underscores AutoDock-SS’s ability to treat compounds as conformationally flexible while considering the ligand’s shape, pharmacophore, and electrostatic potential, expanding the potential of LBVS methods

The kth nearest neighbor method for estimation of entropy changes from molecular ensembles

All processes involving molecular systems entail a balance between associated enthalpic and entropic changes. Molecular dynamics simulations of the end-points of a process provide in a straightforward way the enthalpy as an ensemble average. Obtaining absolute entropies is still an open problem and most commonly pathway methods are used to obtain free energy changes and thereafter entropy changes. The kth nearest neighbor (kNN) method has been first proposed as a general method for entropy estimation in the mathematical community 20 years ago. Later, it has been applied to compute conformational, positional–orientational, and hydration entropies of molecules. Programs to compute entropies from molecular ensembles, for example, from molecular dynamics (MD) trajectories, based on the kNN method, are currently available. The kNN method has distinct advantages over traditional methods, namely that it is possible to address high-dimensional spaces, impossible to treat without loss of resolution or drastic approximations with, for example, histogram-based methods. Application of the method requires understanding the features of: the kth nearest neighbor method for entropy estimation; the variables relevant to biomolecular and in general molecular processes; the metrics associated with such variables; the practical implementation of the method, including requirements and limitations intrinsic to the method; and the applications for conformational, position/orientation and solvation entropy. Coupling the method with general approximations for the multivariable entropy based on mutual information, it is possible to address high dimensional problems like those involving the conformation of proteins, nucleic acids, binding of molecules and hydration. This article is categorized under: Molecular and Statistical Mechanics > Free Energy Methods Theoretical and Physical Chemistry > Statistical Mechanics Structure and Mechanism > Computational Biochemistry and Biophysics