A Fresh Look at the Normal Mode Analysis of Proteins: Introducing Allosteric Co-Vibrational Modes.

We propose a new way of utilizing normal modes to study protein conformational transitions. Instead of considering individual modes independently, we show that a weighted mixture of low-frequency vibrational modes can reveal dynamic information about the conformational mechanism in more detail than any single mode can. The weights in the mixed mode, termed the allosteric covibrational mode, are determined using a simple model where the conformational transition is viewed as a perturbation of the coupled harmonic oscillator associated with either of the two conformations. We demonstrate our theory in a biologically relevant example of high pharmaceutical interest involving the V617F mutation of Janus 2 tyrosine kinase (JAK2).

Comprehensive Approach to Simulating Large Scale Conformational Changes in Biological Systems Utilizing a Path Collective Variable and New Barrier Restraint.

Conformational sampling of complex biomolecules is an emerging frontier in drug discovery. Advances in lab-based structural biology and related computational approaches like AlphaFold have made great strides in obtaining static protein structures for biologically relevant targets. However, biology is in constant motion, and many important biological processes rely on conformationally driven events. Conventional molecular dynamics (MD) simulations run on standard hardware are impractical for many drug design projects, where conformationally driven biological events can take microseconds to milliseconds or longer. An alternative approach is to focus the search on a limited region of conformational space defined by a putative reaction coordinate (i.e., path collective variable). The search space is typically limited by applying restraints, which can be guided by insights about the underlying biological process of interest. The challenge is striking a balance between the degree to which the system is constrained and still allowing for natural motions along the path. A plethora of restraints exist to limit the size of conformational search space, although each has drawbacks when simulating complex biological motions. In this work, we present a three-stage procedure to construct realistic path collective variables (PCVs) and introduce a new kind of barrier restraint that is particularly well suited for complex conformationally driven biological events, such as allosteric modulations and conformational signaling. The PCV presented here is all-atom (as opposed to C-alpha or backbone only) and is derived from all-atom MD trajectory frames. The new restraint relies on a barrier function (specifically, the scaled reciprocal function), which we show is particularly beneficial in the context of molecular dynamics, where near-hard-wall restraints are needed with zero tolerance to restraint violation. We have implemented our PCV and barrier restraint within a hybrid sampling framework that combines well-tempered metadynamics and extended-Lagrangian adaptive biasing force (meta-eABF). We use three particular examples of high pharmaceutical interest to demonstrate the value of this approach: (1) sampling the distance from ubiquitin to a protein of interest within the supramolecular cullin-RING ligase complex, (2) stabilizing the wild-type conformation of the oncogenic mutant JAK2-V617F pseudokinase domain, and (3) inducing an activated state of the stimulator of interferon genes (STING) protein observed upon ligand binding. For examples 2 and 3, we present statistical analysis of meta-eABF free energy estimates and, for each case, code for reproducing this work.

Predicting the structural basis of targeted protein degradation by integrating molecular dynamics simulations with structural mass spectrometry.

Targeted protein degradation (TPD) is a promising approach in drug discovery for degrading proteins implicated in diseases. A key step in this process is the formation of a ternary complex where a heterobifunctional molecule induces proximity of an E3 ligase to a protein of interest (POI), thus facilitating ubiquitin transfer to the POI. In this work, we characterize 3 steps in the TPD process. (1) We simulate the ternary complex formation of SMARCA2 bromodomain and VHL E3 ligase by combining hydrogen-deuterium exchange mass spectrometry with weighted ensemble molecular dynamics (MD). (2) We characterize the conformational heterogeneity of the ternary complex using Hamiltonian replica exchange simulations and small-angle X-ray scattering. (3) We assess the ubiquitination of the POI in the context of the full Cullin-RING Ligase, confirming experimental ubiquitinomics results. Differences in degradation efficiency can be explained by the proximity of lysine residues on the POI relative to ubiquitin.

Structure-Based Analysis of Cryptic-Site Opening.

Many proteins in their unbound structures have cryptic sites that are not appropriately sized for drug binding. We consider here 32 proteins from the recently published CryptoSite set with validated cryptic sites, and study whether the sites remain cryptic in all available X-ray structures of the proteins solved without any ligand bound near the sites. It was shown that only few of these proteins have binding pockets that never form without ligand binding. Sites that are cryptic in some structures but spontaneously form in others are also rare. In most proteins the forming of pockets is affected by mutations or ligand binding at locations far from the cryptic site. To further explore these mechanisms, we applied adiabatic biased molecular dynamics simulations to guide the proteins from their ligand-free structures to ligand-bound conformations, and studied the distribution of druggability scores of the pockets located at the cryptic sites.

Global optimization of additive potential energy functions: predicting binary Lennard-Jones clusters.

We present a method for minimizing additive potential-energy functions. Our hidden-force algorithm can be described as an intricate multiplayer tug-of-war game in which teams try to break an impasse by randomly assigning some players to drop their ropes while the others are still tugging until a partial impasse is reached, then, instructing the dropouts to resume tugging, for all teams to come to a new overall impasse. Utilizing our algorithm in a non-Markovian parallel Monte Carlo search, we found 17 new putative global minima for binary Lennard-Jones clusters in the size range of 90-100 particles. The method is efficient enough that an unbiased search was possible; no potential-energy surface symmetries were exploited. All new minima are comprised of three nested polyicosahedral or polytetrahedral shells when viewed as a nested set of Connolly surfaces (though the shell structure has previously gone unscrutinized, known minima are often qualitatively similar). Unlike known minima, in which the outer and inner shells are comprised of the larger and smaller atoms, respectively, in 13 of the new minima, the atoms are not as clearly separated by size. Furthermore, while some known minima have inner shells stabilized by larger atoms, four of the new minima have outer shells stabilized by smaller atoms.

Mechanism of Na+/H+ antiporting.

Na+/H+ antiporters are central to cellular salt and pH homeostasis. The structure of Escherichia coli NhaA was recently determined, but its mechanisms of transport and pH regulation remain elusive. We performed molecular dynamics simulations of NhaA that, with existing experimental data, enabled us to propose an atomically detailed model of antiporter function. Three conserved aspartates are key to our proposed mechanism: Asp164 (D164) is the Na+-binding site, D163 controls the alternating accessibility of this binding site to the cytoplasm or periplasm, and D133 is crucial for pH regulation. Consistent with experimental stoichiometry, two protons are required to transport a single Na+ ion: D163 protonates to reveal the Na+-binding site to the periplasm, and subsequent protonation of D164 releases Na+. Additional mutagenesis experiments further validated the model.

Analyzing the performance of conformational search programs on compound databases.

We have studied the sampling performance of conformational search programs using geometric and energetic criteria. Ideally, a conformational search algorithm should identify the largest possible number of low-energy structures (energy criterion) covering the widest possible range of molecular shapes (geometric criterion). Geometric analysis consisted in comparing the distribution of conformations within the generated ensembles by multidimensional scaling and by analysing the eigenvalue structure of the pairwise coordinate covariance matrices. The energetic comparison was carried out by assessing the energy distribution of conformers after minimizing them all using the same semi-empirical quantum mechanics optimization protocol. The present investigation focused on five conformational search programs: DGEOM, QXP, ROTATE, LMOD and OMEGA. We have applied these methodologies to a maximally diverse 604-compound subset of the LeadQuest library. The program LMOD performs best according to the energetic criterion, whereas a wider range of geometrically diverse conformations is sampled by the other programs, at the cost of higher median conformer energies. In terms of speed, OMEGA is fastest. We recommend the use of LMOD or OMEGA for high-quality conformational search applications.

Application of the frozen atom approximation to the GB/SA continuum model for solvation free energy.

The generalized Born/surface area (GB/SA) continuum model for solvation free energy is a fast and accurate alternative to using discrete water molecules in molecular simulations of solvated systems. However, computational studies of large solvated molecular systems such as enzyme-ligand complexes can still be computationally expensive even with continuum solvation methods simply because of the large number of atoms in the solute molecules. Because in such systems often only a relatively small portion of the system such as the ligand binding site is under study, it becomes less attractive to calculate energies and derivatives for all atoms in the system. To curtail computation while still maintaining high energetic accuracy, atoms distant from the site of interest are often frozen; that is, their coordinates are made invariant. Such frozen atoms do not require energetic and derivative updates during the course of a simulation. Herein we describe methodology and results for applying the frozen atom approach to both the generalized Born (GB) and the solvent accessible surface area (SASA) parts of the GB/SA continuum model for solvation free energy. For strictly pairwise energetic terms, such as the Coulombic and van-der-Waals energies, contributions from pairs of frozen atoms can be ignored. This leaves energetic differences unaffected for conformations that vary only in the positions of nonfrozen atoms. Due to the nonlocal nature of the GB analytical form, however, excluding such pairs from a GB calculation leads to unacceptable inaccuracies. To apply a frozen-atom scheme to GB calculations, a buffer region within the frozen-atom zone is generated based on a user-definable cutoff distance from the nonfrozen atoms. Certain pairwise interactions between frozen atoms in the buffer region are retained in the GB computation. This allows high accuracy in conformational GB comparisons to be maintained while achieving significant savings in computational time compared to the full (nonfrozen) calculation. A similar approach for using a buffer region of frozen atoms is taken for the SASA calculation. The SASA calculation is local in nature, and thus exact SASA energies are maintained. With a buffer region of 8 A for the frozen-atom cases, excellent agreement in differences in energies for three different conformations of cytochrome P450 with a bound camphor ligand are obtained with respect to the nonfrozen cases. For various minimization protocols, simulations run 2 to 10.5 times faster and memory usage is reduced by a factor of 1.5 to 5. Application of the frozen atom method for GB/SA calculations thus can render computationally tractable biologically and medically important simulations such as those used to study ligand-receptor binding conformations and energies in a solvated environment.

Superspreading driven by Marangoni flow.

The spontaneous spreading (called superspreading) of aqueous trisiloxane ethoxylate surfactant solutions on hydrophobic solid surfaces is a fascinating phenomenon with several practical applications. For example, the ability of trisiloxane ethoxylate surfactants to enhance the spreading of spray solutions on waxy weed leaf surfaces, such as velvetleaf (Abutilion theophrasti), makes them excellent wetting agents for herbicide applications. The superspreading ability of silicone surfactants has been known for decades, but its mechanism is still not well understood. In this paper, we suggest that the spreading of trisiloxane ethoxylates is controlled by a surface tension gradient, which forms when a drop of surfactant solution is placed on a solid surface. The proposed model suggests that, as the spreading front stretches, the surface tension increases (the surfactant concentration becomes lower) at the front relative to the top of the droplet, thereby establishing a dynamic surface tension gradient. The driving force for spreading is due to the Marangoni effect, and our experiments showed that the higher the gradient, the faster the spreading. A simple model describing the phenomenon of superspreading is presented. We also suggest that the superspreading behavior of trisiloxane ethoxylates is a consequence of the molecular configuration at the air/water surface (i.e. small and compact hydrophobic part), as shown by molecular dynamics modeling. We also found that the aggregates and vesicles formed in trisiloxane solutions do not initiate the spreading process and therefore these structures are not a requirement for the superspreading process.