Protein-Ligand Binding from Distancefield Distances and Hamiltonian Replica Exchange Simulations

Methods for Classical-Mechanical Molecular Simulation in Chemistry: Achievements, Limitations, Perspectives

More than a half century ago it became feasible to simulate, using classical-mechanical equations of motion, the dynamics of molecular systems on a computer. Since then classical-physical molecular simulation has become an integral part of chemical research. It is widely applied in a variety of branches of chemistry and has significantly contributed to the development of chemical knowledge. It offers understanding and interpretation of experimental results, semiquantitative predictions for measurable and nonmeasurable properties of substances, and allows the calculation of properties of molecular systems under conditions that are experimentally inaccessible. Yet, molecular simulation is built on a number of assumptions, approximations, and simplifications which limit its range of applicability and its accuracy. These concern the potential-energy function used, adequate sampling of the vast statistical-mechanical configurational space of a molecular system and the methods used to compute particular properties of chemical systems from statistical-mechanical ensembles. During the past half century various methodological ideas to improve the efficiency and accuracy of classical-physical molecular simulation have been proposed, investigated, evaluated, implemented in general simulation software or were abandoned. The latter because of fundamental flaws or, while being physically sound, computational inefficiency. Some of these methodological ideas are briefly reviewed and the most effective methods are highlighted. Limitations of classical-physical simulation are discussed and perspectives are sketched.

A molecularly engineered, broad-spectrum anti-coronavirus lectin inhibits SARS-CoV-2 and MERS-CoV infection in vivo

"Pan-coronavirus" antivirals targeting conserved viral components can be designed. Here, we show that the rationally engineered H84T-banana lectin (H84T-BanLec), which specifically recognizes high mannose found on viral proteins but seldom on healthy human cells, potently inhibits Middle East respiratory syndrome coronavirus (MERS-CoV), severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) (including Omicron), and other human-pathogenic coronaviruses at nanomolar concentrations. H84T-BanLec protects against MERS-CoV and SARS-CoV-2 infection in vivo. Importantly, intranasally and intraperitoneally administered H84T-BanLec are comparably effective. Mechanistic assays show that H84T-BanLec targets virus entry. High-speed atomic force microscopy depicts real-time multimolecular associations of H84T-BanLec dimers with the SARS-CoV-2 spike trimer. Single-molecule force spectroscopy demonstrates binding of H84T-BanLec to multiple SARS-CoV-2 spike mannose sites with high affinity and that H84T-BanLec competes with SARS-CoV-2 spike for binding to cellular ACE2. Modeling experiments identify distinct high-mannose glycans in spike recognized by H84T-BanLec. The multiple H84T-BanLec binding sites on spike likely account for the drug compound's broad-spectrum antiviral activity and the lack of resistant mutants.

Free Energy Landscape of H2A-H2B Displacement From Nucleosome

Nucleosome reconstitution plays an important role in many cellular functions. As an initial step, H2A-H2B dimer displacement, which is accompanied by disruption of many of the interactions within the nucleosome, should occur. To understand how H2A-H2B dimer displacement occurs, an adaptively biased molecular dynamics (ABMD) simulation was carried out to generate a variety of displacements of the H2A-H2B dimer from the fully wrapped to partially unwrapped nucleosome structures. With regards to these structures, the free energy landscape of the dimer displacement was investigated using umbrella sampling simulations. We found that the main contributors to the free energy were the docking domain of H2A and the C-terminal of H4. There were various paths for the dimer displacement which were dependent on the extent of nucleosomal DNA wrapping, suggesting that modulation of the intra-nucleosomal interaction by external factors such as histone chaperons could control the path for the H2A-H2B dimer displacement. Key residues which contributed to the free energy have also been reported to be involved in the mutations and posttranslational modifications (PTMs) which are important for assembling and/or reassembling the nucleosome at the molecular level and are found in cancer cells at the phenotypic level. Our results give insight into how the H2A-H2B dimer displacement proceeds along various paths according to different interactions within the nucleosome.

Chimeric Cellobiose Dehydrogenases Reveal the Function of Cytochrome Domain Mobility for the Electron Transfer to Lytic Polysaccharide Monooxygenase

The natural function of cellobiose dehydrogenase (CDH) to donate electrons from its catalytic flavodehydrogenase (DH) domain via its cytochrome (CYT) domain to lytic polysaccharide monooxygenase (LPMO) is an example of a highly efficient extracellular electron transfer chain. To investigate the function of the CYT domain movement in the two occurring electron transfer steps, two CDHs from the ascomycete Neurospora crassa (NcCDHIIA and NcCDHIIB) and five chimeric CDH enzymes created by domain swapping were studied in combination with the fungus’ own LPMOs (NcLPMO9C and NcLPMO9F). Kinetic and electrochemical methods and hydrogen/deuterium exchange mass spectrometry were used to study the domain movement, interaction, and electron transfer kinetics. Molecular docking provided insights into the protein–protein interface, the orientation of domains, and binding energies. We find that the first, interdomain electron transfer step from the catalytic site in the DH domain to the CYT domain depends on steric and electrostatic interface complementarity and the length of the protein linker between both domains but not on the redox potential difference between the FAD and heme b cofactors. After CYT reduction, a conformational change of CDH from its closed state to an open state allows the second, interprotein electron transfer (IPET) step from CYT to LPMO to occur by direct interaction of the b-type heme and the type-2 copper center. Chimeric CDH enzymes favor the open state and achieve higher IPET rates by exposing the heme b cofactor to LPMO. The IPET, which is influenced by interface complementarity and the heme b redox potential, is very efficient with bimolecular rates between 2.9 × 10⁵ and 1.1 × 10⁶ M^–1 s^–1.

Charge-Changing Perturbations and Path Sampling via Classical Molecular Dynamic Simulations of Simple Guest-Host Systems

Currently, two different methods dominate the field of biomolecular free-energy calculations for the prediction of binding affinities. Pathway methods are frequently used for large ligands that bind on the surface of a host, such as protein-protein complexes. Alchemical methods, on the other hand, are preferably applied for small ligands that bind to deeply buried binding sites. The latter methods are also widely known to be heavily artifacted by the representation of electrostatic energies in periodic simulation boxes, in particular, when net-charge changes are involved. Different methods have been described to deal with these artifacts, including postsimulation correction schemes and instantaneous correction schemes (e.g., co-alchemical perturbation of ions). Here, we use very simple test systems to show that instantaneous correction schemes with no change in the system net charge lower the artifacts but do not eliminate them. Furthermore, we show that free energies from pathway methods suffer from the same artifacts.

Advances in the calculation of binding free energies

In recent years, calculations of binding affinities from molecular simulations seem to have matured significantly. While the number of applications of such methods in drug design and biotechnology increases, the number of truly new methodological developments decreases. This review provides an overview of the current status of the field as reflected in recent publications. The focus is on the challenges that remain when using endstate, alchemical and pathway methods. For endstate methods this is the calculation of entropic contributions. For alchemical methods there are unsolved problems associated with the solvation of the active site, sampling slow degrees of freedom and when modifying the net charge. For pathway methods achieving sufficient sampling remains challenging. New trends are also highlighted, including the use of pathway methods for the quantification of protein-protein interactions.

Correcting electrostatic artifacts due to net-charge changes in the calculation of ligand binding free energies

Alchemically derived free energies are artifacted when the perturbed moiety has a nonzero net charge. The source of the artifacts lies in the effective treatment of the electrostatic interactions within and between the perturbed atoms and remaining (partial) charges in the simulated system. To treat the electrostatic interactions effectively, lattice-summation (LS) methods or cutoff schemes in combination with a reaction-field contribution are usually employed. Both methods render the charging component of the calculated free energies sensitive to essential parameters of the system like the cutoff radius or the box side lengths. Here, we discuss the results of three previously published studies of ligand binding. These studies presented estimates of binding free energies that were artifacted due to the charged nature of the ligands. We show that the size of the artifacts can be efficiently calculated and raw simulation data can be corrected. We compare the corrected results with experimental estimates and nonartifacted estimates from path-sampling methods. Although the employed correction scheme involves computationally demanding continuum-electrostatics calculations, we show that the correction estimate can be deduced from a small sample of configurations rather than from the entire ensemble. This observation makes the calculations of correction terms feasible for complex biological systems. To show the general applicability of the proposed procedure, we also present results where the correction scheme was used to correct independent free energies obtained from simulations employing a cutoff scheme or LS electrostatics. In this work, we give practical guidelines on how to apply the appropriate corrections easily.

String Method for Protein-Protein Binding Free-Energy Calculations

A powerful computational strategy to determine the equilibrium association constant of two macromolecules with explicit-solvent molecular dynamics (MD) simulations is the "geometric route", which considers the reversible physical separation of the bound complex in solution. Nonetheless, multiple challenges remain to render this type of methodology reliable and computationally efficient in practice. In particular, in one, formulation of the geometric route relies on the potential of mean force (PMF) for physically separating the two binding partners restrained along a straight axis, which must be selected prior to the calculation. However, practical applications indicate that the calculation of the separation PMF along the predefined rectilinear pathway may be suboptimal and slowly convergent. Recognizing that a rectilinear straight separation pathway is generally not representative of how the protein complex physically separates in solution, we put forth a novel theoretical framework for binding free-energy calculations, leaning on the optimal curvilinear minimum free-energy path (MFEP) determined from the string method. The proposed formalism is validated by comparing the results obtained using both rectilinear and curvilinear pathways for a prototypical host-guest complex formed by cucurbit[7]uril (CB[7]) binding benzene, and for the barnase-barstar protein complex. On the basis of multi-microsecond MD calculations, we find that the calculations following the traditional rectilinear pathway and the string-based curvilinear pathway agree quantitatively, but convergence is faster with the latter.

Binding Modes and Metabolism of Caffeine

A correct estimate of ligand binding modes and a ratio of their occupancies is crucial for calculations of binding free energies. The newly developed method BLUES combines molecular dynamics with nonequilibrium candidate Monte Carlo. Nonequilibrium candidate Monte Carlo generates a plethora of possible binding modes and molecular dynamics enables the system to relax. We used BLUES to investigate binding modes of caffeine in the active site of its metabolizing enzyme Cytochrome P450 1A2 with the aim of elucidating metabolite-formation profiles at different concentrations. Because the activation energies of all sites of metabolism do not show a clear preference for one metabolite over the others, the orientations in the active site must play a key role. In simulations with caffeine located in a spacious pocket above the I-helix, it points N3 and N1 to the heme iron, whereas in simulations where caffeine is in close proximity to the heme N7 and C8 are preferably oriented toward the heme iron. We propose a mechanism where at low caffeine concentrations caffeine binds to the upper part of the active site, leading to formation of the main metabolite paraxanthine. On the other hand, at high concentrations two molecules are located in the active site, forcing one molecule into close proximity to the heme and yielding metabolites theophylline and trimethyluretic acid. Our results offer an explanation of previously published experimental results.

Free Energy of Binding and Mechanism of Interaction for the MEEVD-TPR2A Peptide-Protein Complex

The association between the MEEVD C-terminal peptide from the heat shock protein 90 (Hsp90) and tetratricopeptide repeat A (TPR2A) domain of the heat shock organizing protein (Hop) is a useful prototype to study the fundamental molecular details about the Hop-Hsp90 interaction. We study here the mechanism of binding/unbinding and compute the standard binding free energy and potential of mean force for the association of the MEEVD peptide to the TPR2A domain using the Adaptive Biasing Force (ABF) methodology. We observe conformational changes of the peptide and the protein receptor induced by binding. We measure the binding free energy of -8.4 kcal/mol, which is consistent with experimental estimates. The simulations achieve multiple unbinding and rebinding events along a consistent pathway connecting the binding site to solvent. The MEEVD peptide slowly dissociates disrupting the hydrogen bonds first, then tilting on the side while preserving the interaction with the side chain of residue Asp 5 of the peptide. After this initial displacement, the peptide completely dissociates and moves into the solvent. Rebinding of the MEEVD peptide from the solvent to the receptor binding site occurs slowly through the portal of entry. Unbinding and rebinding go through intermediate states characterized by the peptide interacting with a lateral helix, helix A1, of the receptor with mainly Asp 5, Val 4, and Glu 3 of the peptide. This newly discovered intermediate structure is characterized by numerous contacts with the receptor which lead to complete formation of the bound complex. The structure of the bound complex obtained after rebinding is structurally very similar to the crystal structure of the complex (0.48 Å root-mean square deviation). The residues Asp 5, Val 4, and Glu 3 adopt conformations and intermolecular contacts with excellent structural similarity with the native ones. Finally, the dissociation and reassociation of MEEVD induce hydration/dehydration transitions, which provide insights on the role of desolvation and solvation processes in protein-peptide binding.

Exploring the binding pathways of the 14-3-3ζ protein: Structural and free-energy profiles revealed by Hamiltonian replica exchange molecular dynamics with distancefield distance restraints

The 14-3-3 protein family performs regulatory functions in eukaryotic organisms by binding to a large number of phosphorylated protein partners. Whilst the binding mode of the phosphopeptides within the primary 14-3-3 binding site is well established based on the crystal structures of their complexes, little is known about the binding process itself. We present a computational study of the process by which phosphopeptides bind to the 14-3-3ζ protein. Applying a novel scheme combining Hamiltonian replica exchange molecular dynamics and distancefield restraints allowed us to map and compare the most likely phosphopeptide-binding pathways to the 14-3-3ζ protein. The most important structural changes to the protein and peptides involved in the binding process were identified. In order to bind phosphopeptides to the primary interaction site, the 14-3-3ζ adopted a newly found wide-opened conformation. Based on our findings we additionally propose a secondary interaction site on the inner surface of the 14-3-3ζ dimer, and a direct interference on the binding process by the flexible C-terminal tail. A minimalistic model was designed to allow for the efficient calculation of absolute binding affinities. Binding affinities calculated from the potential of mean force along the binding pathway are in line with the available experimental estimates for two of the studied systems.

Predicting Binding Free Energies: Frontiers and Benchmarks

Binding free energy calculations based on molecular simulations provide predicted affinities for biomolecular complexes. These calculations begin with a detailed description of a system, including its chemical composition and the interactions among its components. Simulations of the system are then used to compute thermodynamic information, such as binding affinities. Because of their promise for guiding molecular design, these calculations have recently begun to see widespread applications in early-stage drug discovery. However, many hurdles remain in making them a robust and reliable tool. In this review, we highlight key challenges of these calculations, discuss some examples of these challenges, and call for the designation of standard community benchmark test systems that will help the research community generate and evaluate progress. In our view, progress will require careful assessment and evaluation of new methods, force fields, and modeling innovations on well-characterized benchmark systems, and we lay out our vision for how this can be achieved.

Free-energy calculations of residue mutations in a tripeptide using various methods to overcome inefficient sampling

Previous free-energy calculations have shown that the seemingly simple transformation of the tripeptide KXK to KGK in water holds some unobvious challenges concerning the convergence of the forward and backward thermodynamic integration processes (i.e., hysteresis). In the current study, the central residue X was either alanine, serine, glutamic acid, lysine, phenylalanine, or tyrosine. Interestingly, the transformation from alanine to glycine yielded the highest hysteresis in relation to the extent of the chemical change of the side chain. The reason for that could be attributed to poor sampling of φ2 /ψ2 dihedral angles along the transformation. Altering the nature of alanine's Cβ atom drastically improved the sampling and at the same time led to the identification of high energy barriers as cause for it. Consequently, simple strategies to overcome these barriers are to increase simulation time (computationally expensive) or to use enhanced sampling techniques such as Hamiltonian replica exchange molecular dynamics and one-step perturbation. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Multistate λ-local-elevation umbrella-sampling (MS-λ-LEUS): method and application to the complexation of cations by crown ethers

An extension of the λ-local-elevation umbrella-sampling (λ-LEUS) scheme [ Bieler et al. J. Chem. Theory Comput. 2014 , 10 , 3006 ] is proposed to handle the multistate (MS) situation, i.e. the calculation of the relative free energies of multiple physical states based on a single simulation. The key element of the MS-λ-LEUS approach is to use a single coupling variable Λ controlling successive pairwise mutations between the states of interest in a cyclic fashion. The Λ variable is propagated dynamically as an extended-system variable, using a coordinate transformation with plateaus and a memory-based biasing potential as in λ-LEUS. Compared to other available MS schemes (one-step perturbation, enveloping distribution sampling and conventional λ-dynamics) the proposed method presents a number of important advantages, namely: (i) the physical states are visited explicitly and over finite time periods; (ii) the extent of unphysical space required to ensure transitions is kept minimal and, in particular, one-dimensional; (iii) the setup protocol solely requires the topologies of the physical states; and (iv) the method only requires limited modifications in a simulation code capable of handling two-state mutations. As an initial application, the absolute binding free energies of five alkali cations to three crown ethers in three different solvents are calculated. The results are found to reproduce qualitatively the main experimental trends and, in particular, the experimental selectivity of 18C6 for K(+) in water and methanol, which is interpreted in terms of opposing trends along the cation series between the solvation free energy of the cation and the direct electrostatic interactions within the complex.

Exploring biomolecular dynamics and interactions using advanced sampling methods

Molecular dynamics (MD) and Monte Carlo (MC) simulations have emerged as a valuable tool to investigate statistical mechanics and kinetics of biomolecules and synthetic soft matter materials. However, major limitations for routine applications are due to the accuracy of the molecular mechanics force field and due to the maximum simulation time that can be achieved in current simulations studies. For improving the sampling a number of advanced sampling approaches have been designed in recent years. In particular, variants of the parallel tempering replica-exchange methodology are widely used in many simulation studies. Recent methodological advancements and a discussion of specific aims and advantages are given. This includes improved free energy simulation approaches and conformational search applications.

Identification of Possible Binding Sites for Morphine and Nicardipine on the Multidrug Transporter P-Glycoprotein Using Umbrella Sampling Techniques

The multidrug transporter P-glycoprotein (P-gp) is central to the development of multidrug resistance in cancer. While residues essential for transport and binding have been identified, the location, composition, and specificity of potential drug binding sites are uncertain. Here molecular dynamics simulations are used to calculate the free energy profile for the binding of morphine and nicardipine to P-gp. We show that morphine and nicardipine primarily interact with key residues implicated in binding and transport from mutational studies, binding at different but overlapping sites within the transmembrane pore. Their permeation pathways were distinct but involved overlapping sets of residues. The results indicate that the binding location and permeation pathways of morphine and nicardipine are not well separated and cannot be considered as unique. This has important implications for our understanding of substrate uptake and transport by P-gp. Our results are independent of the choice of starting structure and consistent with a range of experimental studies.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Net charge changes in the calculation of relative ligand-binding free energies via classical atomistic molecular dynamics simulation

The calculation of binding free energies of charged species to a target molecule is a frequently encountered problem in molecular dynamics studies of (bio-)chemical thermodynamics. Many important endogenous receptor-binding molecules, enzyme substrates, or drug molecules have a nonzero net charge. Absolute binding free energies, as well as binding free energies relative to another molecule with a different net charge will be affected by artifacts due to the used effective electrostatic interaction function and associated parameters (e.g., size of the computational box). In the present study, charging contributions to binding free energies of small oligoatomic ions to a series of model host cavities functionalized with different chemical groups are calculated with classical atomistic molecular dynamics simulation. Electrostatic interactions are treated using a lattice-summation scheme or a cutoff-truncation scheme with Barker-Watts reaction-field correction, and the simulations are conducted in boxes of different edge lengths. It is illustrated that the charging free energies of the guest molecules in water and in the host strongly depend on the applied methodology and that neglect of correction terms for the artifacts introduced by the finite size of the simulated system and the use of an effective electrostatic interaction function considerably impairs the thermodynamic interpretation of guest-host interactions. Application of correction terms for the various artifacts yields consistent results for the charging contribution to binding free energies and is thus a prerequisite for the valid interpretation or prediction of experimental data via molecular dynamics simulation. Analysis and correction of electrostatic artifacts according to the scheme proposed in the present study should therefore be considered an integral part of careful free-energy calculation studies if changes in the net charge are involved.

Calculating the binding free energies of charged species based on explicit-solvent simulations employing lattice-sum methods: an accurate correction scheme for electrostatic finite-size effects

The calculation of a protein-ligand binding free energy based on molecular dynamics (MD) simulations generally relies on a thermodynamic cycle in which the ligand is alchemically inserted into the system, both in the solvated protein and free in solution. The corresponding ligand-insertion free energies are typically calculated in nanoscale computational boxes simulated under periodic boundary conditions and considering electrostatic interactions defined by a periodic lattice-sum. This is distinct from the ideal bulk situation of a system of macroscopic size simulated under non-periodic boundary conditions with Coulombic electrostatic interactions. This discrepancy results in finite-size effects, which affect primarily the charging component of the insertion free energy, are dependent on the box size, and can be large when the ligand bears a net charge, especially if the protein is charged as well. This article investigates finite-size effects on calculated charging free energies using as a test case the binding of the ligand 2-amino-5-methylthiazole (net charge +1 e) to a mutant form of yeast cytochrome c peroxidase in water. Considering different charge isoforms of the protein (net charges -5, 0, +3, or +9 e), either in the absence or the presence of neutralizing counter-ions, and sizes of the cubic computational box (edges ranging from 7.42 to 11.02 nm), the potentially large magnitude of finite-size effects on the raw charging free energies (up to 17.1 kJ mol(-1)) is demonstrated. Two correction schemes are then proposed to eliminate these effects, a numerical and an analytical one. Both schemes are based on a continuum-electrostatics analysis and require performing Poisson-Boltzmann (PB) calculations on the protein-ligand system. While the numerical scheme requires PB calculations under both non-periodic and periodic boundary conditions, the latter at the box size considered in the MD simulations, the analytical scheme only requires three non-periodic PB calculations for a given system, its dependence on the box size being analytical. The latter scheme also provides insight into the physical origin of the finite-size effects. These two schemes also encompass a correction for discrete solvent effects that persists even in the limit of infinite box sizes. Application of either scheme essentially eliminates the size dependence of the corrected charging free energies (maximal deviation of 1.5 kJ mol(-1)). Because it is simple to apply, the analytical correction scheme offers a general solution to the problem of finite-size effects in free-energy calculations involving charged solutes, as encountered in calculations concerning, e.g., protein-ligand binding, biomolecular association, residue mutation, pKa and redox potential estimation, substrate transformation, solvation, and solvent-solvent partitioning.

Navigating traditional chinese medicine network pharmacology and computational tools

The concept of "network target" has ushered in a new era in the field of traditional Chinese medicine (TCM). As a new research approach, network pharmacology is based on the analysis of network models and systems biology. Taking advantage of advancements in systems biology, a high degree of integration data analysis strategy and interpretable visualization provides deeper insights into the underlying mechanisms of TCM theories, including the principles of herb combination, biological foundations of herb or herbal formulae action, and molecular basis of TCM syndromes. In this study, we review several recent developments in TCM network pharmacology research and discuss their potential for bridging the gap between traditional and modern medicine. We briefly summarize the two main functional applications of TCM network models: understanding/uncovering and predicting/discovering. In particular, we focus on how TCM network pharmacology research is conducted and highlight different computational tools, such as network-based and machine learning algorithms, and sources that have been proposed and applied to the different steps involved in the research process. To make network pharmacology research commonplace, some basic network definitions and analysis methods are presented.