Distant residues modulate conformational opening in SARS-CoV-2 spike protein.

Infection by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) involves the attachment of the receptor-binding domain (RBD) of its spike proteins to the ACE2 receptors on the peripheral membrane of host cells. Binding is initiated by a down-to-up conformational change in the spike protein, the change that presents the RBD to the receptor. To date, computational and experimental studies that search for therapeutics have concentrated, for good reason, on the RBD. However, the RBD region is highly prone to mutations, and is therefore a hotspot for drug resistance. In contrast, we here focus on the correlations between the RBD and residues distant to it in the spike protein. This allows for a deeper understanding of the underlying molecular recognition events and prediction of the highest-effect key mutations in distant, allosteric sites, with implications for therapeutics. Also, these sites can appear in emerging mutants with possibly higher transmissibility and virulence, and preidentifying them can give clues for designing pan-coronavirus vaccines against future outbreaks. Our model, based on time-lagged independent component analysis (tICA) and protein graph connectivity network, is able to identify multiple residues that exhibit long-distance coupling with the RBD opening. Residues involved in the most ubiquitous D614G mutation and the A570D mutation of the highly contagious UK SARS-CoV-2 variant are predicted ab initio from our model. Conversely, broad-spectrum therapeutics like drugs and monoclonal antibodies can target these key distant-but-conserved regions of the spike protein.

Predicting residue cooperativity during protein folding: A combined, molecular dynamics and unsupervised learning approach.

Allostery in proteins involves, broadly speaking, ligand-induced conformational transitions that modulate function at active sites distal to where the ligand binds. In contrast, the concept of cooperativity (in the sense used in phase transition theory) is often invoked to understand protein folding and, therefore, function. The modern view on allostery is one based on dynamics and hinges on the time-dependent interactions between key residues in a complex network, interactions that determine the free-energy profile for the reaction at the distal site. Here, we merge allostery and cooperativity, and we discuss a joint model with features of both. In our model, the active-site reaction is replaced by the reaction pathway that leads to protein folding, and the presence or absence of the effector is replaced by mutant-vs-wild type changes in key residues. To this end, we employ our recently introduced time-lagged independent component analysis (tICA) correlation approach [Ray et al. Proc. Natl. Acad. Sci. 118(43) (2021), e2100943118] to identify the allosteric role of distant residues in the folded-state dynamics of a large protein. In this work, we apply the technique to identify key residues that have a significant role in the folding of a small, fast folding-protein, chignolin. Using extensive enhanced sampling simulations, we critically evaluate the accuracy of the predictions by mutating each residue one at a time and studying how the mutations change the underlying free energy landscape of the folding process. We observe that mutations in those residues whose associated backbone torsion angles have a high correlation score can indeed lead to loss of stability of the folded configuration. We also provide a rationale based on interaction energies between individual residues with the rest of the protein to explain this effect. From these observations, we conclude that the tICA correlation score metric is a useful tool for predicting the role of individual residues in the correlated dynamics of proteins and can find application to the problem of identifying regions of protein that are either most vulnerable to mutations or-mutatis mutandis-to binding events that affect their functionality.

Force-Field-Dependent DNA Breathing Dynamics: A Case Study of Hoogsteen Base Pairing in A6-DNA.

The Hoogsteen (HG) base pairing conformation, commonly observed in damaged and mutated DNA helices, facilitates DNA repair and DNA recognition. The free energy difference between HG and Watson-Crick (WC) base pairs has been computed in previous studies. However, the mechanism of the conformational transition is not well understood. A detailed understanding of the process of WC to HG base pair transition can provide a deeper understanding of DNA repair and recognition. In an earlier study, we explored the free energy landscape for this process using extensive computer simulation with the CHARMM36 force field. In this work, we study the impact of force field models in describing the WC to HG base pairing transition using meta-eABF enhanced sampling, quasi-harmonic entropy calculation, and nonbonded energy analysis. The secondary structures of both base pairing forms and the topology of the free energy landscapes were consistent over different force field models, although the relative free energy, entropy, and the interaction energies tend to vary. The relative stability of the WC and HG conformations is dictated by a delicate balance between the enthalpic stabilization and the reduced entropy of the structurally rigid HG structure. These findings highlight the impact that subtleties in force field models can have on accurately modeling DNA base pair dynamics and should stimulate further computational investigations into other dynamically important motions in DNA.

How the phage T4 injection machinery works including energetics, forces, and dynamic pathway.

The virus bacteriophage T4, from the family Myoviridae, employs an intriguing contractile injection machine to inject its genome into the bacterium Escherichia coli Although the atomic structure of phage T4 is largely understood, the dynamics of its injection machinery remains unknown. This study contributes a system-level model describing the nonlinear dynamics of the phage T4 injection machinery interacting with a host cell. The model employs a continuum representation of the contractile sheath using elastic constants inferred from atomistic molecular-dynamics (MD) simulations. Importantly, the sheath model is coupled to component models representing the remaining structures of the virus and the host cell. The resulting system-level model captures virus-cell interactions as well as competing energetic mechanisms that release and dissipate energy during the injection process. Simulations reveal the dynamical pathway of the injection process as a "contraction wave" that propagates along the sheath, the energy that powers the injection machinery, the forces responsible for piercing the host cell membrane, and the energy dissipation that controls the timescale of the injection process. These results from the model compare favorably with the available (but limited) experimental measurements.

Free Energy Landscape and Conformational Kinetics of Hoogsteen Base Pairing in DNA vs. RNA.

Genetic information is encoded in the DNA double helix, which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleo-base pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen (HG) base pairing pattern in naked DNA duplexes, and estimated its relative stability and lifetime. In contrast with DNA, such structures were not observed in RNA duplexes. Understanding HG base pairing is important because the underlying "breathing" motion between the two conformations can significantly modulate protein binding. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force-bias method) and Markov state modeling to obtain accurate free energy, kinetics, and the intermediates in the transition pathway between Watson-Crick and HG base pairs for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extrahelical intermediates in the seemingly simple process of base flipping in B-DNA. Extending our calculation to A-RNA, for which HG base pairing is not observed experimentally, resulted in relatively unstable, single-hydrogen-bonded, distorted Hoogsteen-like bases. Unlike B-DNA, the transition pathway primarily involved base paired and intrahelical intermediates with transition timescales much longer than that of B-DNA. The seemingly obvious flip-over reaction coordinate (i.e., the glycosidic torsion angle) is unable to resolve the intermediates. Instead, a multidimensional picture involving backbone dihedral angles and distance between hydrogen bond donor and acceptor atoms is required to gain insight into the molecular mechanism.

Weighted ensemble milestoning (WEM): A combined approach for rare event simulations.

To directly simulate rare events using atomistic molecular dynamics is a significant challenge in computational biophysics. Well-established enhanced-sampling techniques do exist to obtain the thermodynamic functions for such systems. However, developing methods for obtaining the kinetics of long timescale processes from simulation at atomic detail is comparatively less developed an area. Milestoning and the weighted ensemble (WE) method are two different stratification strategies; both have shown promise for computing long timescales of complex biomolecular processes. Nevertheless, both require a significant investment of computational resources. We have combined WE and milestoning to calculate observables in orders-of-magnitude less central processing unit and wall-clock time. Our weighted ensemble milestoning method (WEM) uses WE simulation to converge the transition probability and first passage times between milestones, followed by the utilization of the theoretical framework of milestoning to extract thermodynamic and kinetic properties of the entire process. We tested our method for a simple one-dimensional double-well potential, for an eleven-dimensional potential energy surface with energy barrier, and on the biomolecular model system alanine dipeptide. We were able to recover the free energy profiles, time correlation functions, and mean first passage times for barrier crossing events at a significantly small computational cost. WEM promises to extend the applicability of molecular dynamics simulation to slow dynamics of large systems that are well beyond the scope of present day brute-force computations.

Kinetics and free energy of ligand dissociation using weighted ensemble milestoning.

We consider the recently developed weighted ensemble milestoning (WEM) scheme [D. Ray and I. Andricioaei, J. Chem. Phys. 152, 234114 (2020)] and test its capability of simulating ligand-receptor dissociation dynamics. We performed WEM simulations on the following host-guest systems: Na+/Cl- ion pair and 4-hydroxy-2-butanone ligand with FK506 binding protein. As a proof of principle, we show that the WEM formalism reproduces the Na+/Cl- ion pair dissociation timescale and the free energy profile obtained from long conventional MD simulation. To increase the accuracy of WEM calculations applied to kinetics and thermodynamics in protein-ligand binding, we introduced a modified WEM scheme called weighted ensemble milestoning with restraint release (WEM-RR), which can increase the number of starting points per milestone without adding additional computational cost. WEM-RR calculations obtained a ligand residence time and binding free energy in agreement with experimental and previous computational results. Moreover, using the milestoning framework, the binding time and rate constants, dissociation constants, and committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant (kon) when binding is primarily diffusion driven. We show that the WEM method can efficiently calculate multiple experimental observables describing ligand-receptor binding/unbinding and is a promising candidate for computer-aided inhibitor design.

Advances in milestoning. II. Calculating time-correlation functions from milestoning using stochastic path integrals.

In the milestoning framework, and more generally in related transition interface sampling schemes, one significantly enhances the calculation of relaxation rates for complex equilibrium kinetics from molecular dynamics simulations between the milestones or interfaces. The goal of the present paper is to advance milestoning applications into the realm of non-equilibrium statistical mechanics, in particular, to calculate entire time correlation functions. In order to accomplish this, we introduce a novel methodology for obtaining the flux through a given milestone configuration as a function of both time and initial configuration and build upon it with a novel formalism describing autocorrelation for Langevin motion in a discrete configuration space. The method is then applied to three different test systems: a harmonic oscillator, which we solve analytically, a two-well potential, which is solved numerically, and an atomistic molecular dynamics simulation of alanine dipeptide.

Advances in milestoning. I. Enhanced sampling via wind-assisted reweighted milestoning (WARM).

The milestoning algorithm of Elber and co-workers creates a framework for computing the time scale of processes that are too long and too complex to be studied using simply brute force simulations. The fundamental objects involved in the milestoning algorithm are the first passage time distributions KAB (τ) between adjacent conformational milestones A and B. The method proposed herein aims to further enhance milestoning (or other interface based sampling methods) by employing an artificially applied force, akin to a wind that blows the trajectories from their initial to their final states, and by subsequently applying corrective weights to the trajectories to yield the true first passage time distributions KAB (τ) in a fraction of the computation time required for unassisted calculations. The re-weighting method is rooted in the formalism of stochastic path integrals. The theoretical basis for the technique and numerical examples are presented.

Electric-Field-Induced Protein Translocation via a Conformational Transition in SecDF: An MD Study.

SecDF is an important component of the Sec protein translocation machinery embedded in the bacterial membrane, which is associated with many functions, such as stabilizing other Sec translocon components within the membrane, maintaining the transmembrane (TM) potential, and facilitating the ATP-independent stage of the translocation mechanism. Related studies suggest that SecDF undergoes functionally important conformational changes that involve mainly its P1-head domain and that these changes are coupled with the proton motive force (Δp). However, there still is not a clear understanding of how SecDF functions, its exact role in the translocation machinery, and how its function is related to Δp. Here, using all-atom molecular dynamics simulations combined with umbrella sampling, we study the P1-head conformational change and how it is coupled to the proton motive force. We report potentials of mean force along a root-mean-square-distance-based reaction coordinate obtained in the presence and absence of the TM electrical potential. Our results show that the interaction of the P1 domain dipole moment with the TM electrical field considerably lowers the free-energy barrier in the direction of F-form to I-form transition.

Dynamic Model Exposes the Energetics and Dynamics of the Injection Machinery for Bacteriophage T4.

Bacteriophage T4 infects the bacterial host (Escherichia coli) using an efficient genomic delivery machine that is driven by elastic energy stored in a contractile tail sheath. Although the atomic structure of T4 is largely known, the dynamics of its fascinating injection machinery is not understood. This article contributes, to our knowledge, the first predictions of the energetics and dynamics of the T4 injection machinery using a novel dynamic model. The model employs an atomistic (molecular dynamics) representation of a fraction of the sheath structure to generate a continuum model of the entire sheath that also couples to a model of the viral capsid and tail tube. The resulting model of the entire injection machine reveals estimates for the energetics, timescale, and pathway of the T4 injection process as well as the force available for cell rupture. It also reveals the large and highly nonlinear conformational changes of the sheath whose elastic energy drives the injection process.

Transient Hoogsteen base pairs in canonical duplex DNA.

Sequence-directed variations in the canonical DNA double helix structure that retain Watson-Crick base-pairing have important roles in DNA recognition, topology and nucleosome positioning. By using nuclear magnetic resonance relaxation dispersion spectroscopy in concert with steered molecular dynamics simulations, we have observed transient sequence-specific excursions away from Watson-Crick base-pairing at CA and TA steps inside canonical duplex DNA towards low-populated and short-lived A•T and G•C Hoogsteen base pairs. The observation of Hoogsteen base pairs in DNA duplexes specifically bound to transcription factors and in damaged DNA sites implies that the DNA double helix intrinsically codes for excited state Hoogsteen base pairs as a means of expanding its structural complexity beyond that which can be achieved based on Watson-Crick base-pairing. The methods presented here provide a new route for characterizing transient low-populated nucleic acid structures, which we predict will be abundant in the genome and constitute a second transient layer of the genetic code.

Automated placement of interfaces in conformational kinetics calculations using machine learning.

Several recent implementations of algorithms for sampling reaction pathways employ a strategy for placing interfaces or milestones across the reaction coordinate manifold. Interfaces can be introduced such that the full feature space describing the dynamics of a macromolecule is divided into Voronoi (or other) cells, and the global kinetics of the molecular motions can be calculated from the set of fluxes through the interfaces between the cells. Although some methods of this type are exact for an arbitrary set of cells, in practice, the calculations will converge fastest when the interfaces are placed in regions where they can best capture transitions between configurations corresponding to local minima. The aim of this paper is to introduce a fully automated machine-learning algorithm for defining a set of cells for use in kinetic sampling methodologies based on subdividing the dynamical feature space; the algorithm requires no intuition about the system or input from the user and scales to high-dimensional systems.

Dendrimers in Nanoscale Confinement: The Interplay between Conformational Change and Nanopore Entrance.

Hyperbranched dendrimers are nanocarriers for drugs, imaging agents, and catalysts. Their nanoscale confinement is of fundamental interest and occurs when dendrimers with bioactive payload block or pass biological nanochannels or when catalysts are entrapped in inorganic nanoporous support scaffolds. The molecular process of confinement and its effect on dendrimer conformations are, however, poorly understood. Here, we use single-molecule nanopore measurements and molecular dynamics simulations to establish an atomically detailed model of pore dendrimer interactions. We discover and explain that electrophoretic migration of polycationic PAMAM dendrimers into confined space is not dictated by the diameter of the branched molecules but by their size and generation-dependent compressibility. Differences in structural flexibility also rationalize the apparent anomaly that the experimental nanopore current read-out depends in nonlinear fashion on dendrimer size. Nanoscale confinement is inferred to reduce the protonation of the polycationic structures. Our model can likely be expanded to other dendrimers and be applied to improve the analysis of biophysical experiments, rationally design functional materials such as nanoporous filtration devices or nanoscale drug carriers that effectively pass biological pores.

Free-energy landscape and characteristic forces for the initiation of DNA unzipping.

DNA unzipping, the separation of its double helix into single strands, is crucial in modulating a host of genetic processes. Although the large-scale separation of double-stranded DNA has been studied with a variety of theoretical and experimental techniques, the minute details of the very first steps of unzipping are still unclear. Here, we use atomistic molecular-dynamics simulations, coarse-grained simulations, and a statistical-mechanical model to study the initiation of DNA unzipping by an external force. Calculation of the potential of mean force profiles for the initial separation of the first few terminal basepairs in a DNA oligomer revealed that forces ranging between 130 and 230 pN are needed to disrupt the first basepair, and these values are an order of magnitude larger than those needed to disrupt basepairs in partially unzipped DNA. The force peak has an echo of ∼50 pN at the distance that unzips the second basepair. We show that the high peak needed to initiate unzipping derives from a free-energy basin that is distinct from the basins of subsequent basepairs because of entropic contributions, and we highlight the microscopic origin of the peak. To our knowledge, our results suggest a new window of exploration for single-molecule experiments.

Slowdown of Interhelical Motions Induces a Glass Transition in RNA.

RNA function depends crucially on the details of its dynamics. The simplest RNA dynamical unit is a two-way interhelical junction. Here, for such a unit--the transactivation response RNA element--we present evidence from molecular dynamics simulations, supported by nuclear magnetic resonance relaxation experiments, for a dynamical transition near 230 K. This glass transition arises from the freezing out of collective interhelical motional modes. The motions, resolved with site-specificity, are dynamically heterogeneous and exhibit non-Arrhenius relaxation. The microscopic origin of the glass transition is a low-dimensional, slow manifold consisting largely of the Euler angles describing interhelical reorientation. Principal component analysis over a range of temperatures covering the glass transition shows that the abrupt slowdown of motion finds its explanation in a localization transition that traps probability density into several disconnected conformational pools over the low-dimensional energy landscape. Upon temperature increase, the probability density pools then flood a larger basin, akin to a lakes-to-sea transition. Simulations on transactivation response RNA are also used to backcalculate inelastic neutron scattering data that match previous inelastic neutron scattering measurements on larger and more complex RNA structures and which, upon normalization, give temperature-dependent fluctuation profiles that overlap onto a glass transition curve that is quasi-universal over a range of systems and techniques.

Membrane Protein Structure, Function, and Dynamics: a Perspective from Experiments and Theory.

Membrane proteins mediate processes that are fundamental for the flourishing of biological cells. Membrane-embedded transporters move ions and larger solutes across membranes; receptors mediate communication between the cell and its environment and membrane-embedded enzymes catalyze chemical reactions. Understanding these mechanisms of action requires knowledge of how the proteins couple to their fluid, hydrated lipid membrane environment. We present here current studies in computational and experimental membrane protein biophysics, and show how they address outstanding challenges in understanding the complex environmental effects on the structure, function, and dynamics of membrane proteins.

Rate turnover in mechano-catalytic coupling: A model and its microscopic origin.

A novel aspect in the area of mechano-chemistry concerns the effect of external forces on enzyme activity, i.e., the existence of mechano-catalytic coupling. Recent experiments on enzyme-catalyzed disulphide bond reduction in proteins under the effect of a force applied on the termini of the protein substrate reveal an unexpected biphasic force dependence for the bond cleavage rate. Here, using atomistic molecular dynamics simulations combined with Smoluchowski theory, we propose a model for this behavior. For a broad range of forces and systems, the model reproduces the experimentally observed rates by solving a reaction-diffusion equation for a "protein coordinate" diffusing in a force-dependent effective potential. The atomistic simulations are used to compute, from first principles, the parameters of the model via a quasiharmonic analysis. Additionally, the simulations are also used to provide details about the microscopic degrees of freedom that are important for the underlying mechano-catalysis.

Structural ensemble and dynamics of toroidal-like DNA shapes in bacteriophage ϕ29 exit cavity.

In the bacteriophage ϕ29, DNA is packed into a preassembled capsid from which it ejects under high pressure. A recent cryo-EM reconstruction of ϕ29 revealed a compact toroidal DNA structure (30-40 basepairs) lodged within the exit cavity formed by the connector-lower collar protein complex. Using multiscale models, we compute a detailed structural ensemble of intriguing DNA toroids of various lengths, all highly compatible with experimental observations. In particular, coarse-grained (elastic rod) and atomistic (molecular dynamics) models predict the formation of DNA toroids under significant compression, a largely unexplored state of DNA. Model predictions confirm that a biologically attainable compressive force of 25 pN sustains the toroid and yields DNA electron density maps highly consistent with the experimental reconstruction. The subsequent simulation of dynamic toroid ejection reveals large reactions on the connector that may signal genome release.

A general method for constructing atomic-resolution RNA ensembles using NMR residual dipolar couplings: the basis for interhelical motions revealed.

The ability to modulate alignment and measure multiple independent sets of NMR residual dipolar couplings (RDCs) has made it possible to characterize internal motions in proteins at atomic resolution and with time scale sensitivity ranging from picoseconds up to milliseconds. The application of such methods to the study of RNA dynamics, however, remains fundamentally limited by the inability to modulate alignment and by strong couplings between internal and overall motions that complicate the quantitative interpretation of RDCs. Here, we address this problem by showing that RNA alignment can be generally modulated, in a controlled manner, by variable elongation of A-form helices and that the information contained within the measured RDCs can be extracted even in the presence of strong couplings between motions and overall alignment via structure-based prediction of alignment. Using this approach, four RDC data sets, and a broad conformational pool obtained from a 8.2 µs molecular dynamics simulation, we successfully construct and validate an atomic resolution ensemble of human immunodeficiency virus type I transactivation response element RNA. This ensemble reveals local motions in and around the bulge involving changes in stacking and hydrogen-bonding interactions, which are undetectable by traditional spin relaxation and drive global changes in interhelical orientation. This new approach broadens the scope of using RDCs in characterizing the dynamics of nucleic acids.