Search

Monoclonal antibodies are emerging as a viable treatment for the coronavirus disease 19 (COVID-19). However, newly evolved variants of the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) can reduce the efficacy of currently available antibodies and can diminish vaccine-induced immunity. Here, we demonstrate that the microscopic dynamics of neutralizing monoclonal antibodies can be profoundly modified by the mutations present in the spike proteins of the SARS-COV-2 variants currently circulating in the world population. The dynamical perturbations within the antibody structure, which alter the thermodynamics of antigen recognition, are diverse and can depend both on the nature of the antibody and on the spatial location of the spike mutation. The correlation between the motion of the antibody and that of the spike receptor binding domain (RBD) can also be changed, modulating binding affinity. Using protein-graph-connectivity networks, we delineated the mutant-induced modifications in the information-flow along allosteric pathway throughout the antibody. Changes in the collective dynamics were spatially distributed both locally and across long-range distances within the antibody. On the receptor side, we identified an anchor-like structural element that prevents the detachment of the antibodies; individual mutations there can significantly affect the antibody binding propensity. Our study provides insight into how virus neutralization by monoclonal antibodies can be impacted by local mutations in the epitope via a change in dynamics. This realization adds a new layer of sophistication to the efforts for rational design of monoclonal antibodies against new variants of SARS-CoV2, taking the allostery in the antibody into consideration.

We present six advanced tutorials instructing users in the best practices of using key new features and plugins/extensions of the WESTPA 2.0 software package, which consists of major upgrades for enabling applications of the weighted ensemble (WE) path sampling strategy to even larger systems and/or slower processes. The tutorials demonstrate the use of the following key features: (i) a generalized resampler module for the creation of “binless” schemes, (ii) a minimal adaptive binning scheme for more efficient surmounting of free energy barriers, (iii) streamlined handling of large simulation datasets using an HDF5 framework, (iv) two different schemes for more efficient rate-constant estimation, (v) a Python API for simplified analysis of WE simulations, and (vi) plugins/extensions for Markovian Weighted Ensemble Milestoning and WE rule-based modeling at the system biology level. Applications of the tutorials range from atomistic to residue-level to non-spatial models, and include complex processes such as protein folding and the membrane permeability of a drug-like molecule. Users are expected to already have significant experience with running conventional molecular dynamics simulations and completed the previous suite of WESTPA tutorials.

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.

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.

We introduce a rare-event sampling scheme, named Markovian Weighted Ensemble Milestoning (M-WEM), which inlays a weighted ensemble framework within a Markovian milestoning theory to efficiently calculate thermodynamic and kinetic properties of long-timescale biomolecular processes from short atomistic molecular dynamics simulations. M-WEM is tested on the Müller-Brown potential model, the conformational switching in alanine dipeptide, and the millisecond timescale protein-ligand unbinding in a trypsin-benzamidine complex. Not only can M-WEM predict the kinetics of these processes with quantitative accuracy, but it also allows for a scheme to reconstruct a multidimensional free energy landscape along additional degrees of freedom which are not part of the milestoning progress coordinate. For the ligand-receptor system, the experimental residence time, association and dissociation kinetics, and binding free energy could be reproduced using M-WEM within a simulation time of a few hundreds of nanoseconds, which is a fraction of the computational cost of other currently available methods, and close to four orders of magnitude less than the experimental residence time. Due to the high accuracy and low computational cost, the M-WEM approach can find potential application in kinetics and free-energy based computational drug design.

We consider the recently developed weighted ensemble milestoning (WEM) scheme [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 (BUT) ligand with FK506 binding protein (FKBP). As proof or 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 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 constant and the committor probabilities could also be calculated at a low computational cost. We also present an analytical approach for estimating the association rate constant ($k_{\text{on}}$) 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.

Genetic information is encoded in the DNA double helix which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleobase pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen base pairing pattern in naked DNA duplexes and estimated its relative stability and lifetime. In contrast, HG transitions in RNA were not observed. Understanding Hoogsteen (HG) base pairing is important because the underlying "breathing" can modulate significantly DNA/RNA recognition by proteins. 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 WC and HG base pair for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extra-helical intermediates in this seemingly simple process of base pair conformation switching in B-DNA. Extending our calculation to A-RNA, for which HG base pair is not observed experimentally, resulted in relatively unstable single hydrogen bonded distorted Hoogsteen like base pair. Unlike B-DNA the transition pathway primarily involved base paired and intra-helical intermediates with transition timescales much higher than that of B-DNA. The seemingly obvious flip-over reaction coordinate, i.e., the glycosidic torsion angle is unable to resolve the intermediates; so a multidimensional picture, involving backbone dihedral angles and distance between atoms participating in hydrogen bonds, is required to gain insight into the molecular mechanism.SIGNIFICANCEFormation of unconventional Hoogsteen (HG) base pairing is an important problem in DNA biophysics owing to its key role in facilitating the binding of DNA repairing enzymes, proteins and drugs to damaged DNA. X-ray crystallography and NMR relaxation experiments revealed the presence of HG base pair in naked DNA duplex and protein-DNA complex but no HG base pair was observed in RNA. Molecular dynamics simulations could reproduce the experimental free energy cost of HG base pairing in DNA although a detailed mechanistic insight is still missing. We performed enhanced sampling simulation and Markov state modeling to obtain accurate free energy, kinetics and the intermediates in the transition pathway between WC and HG base pair for both B-DNA and A-RNA.

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.

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.

Contractile tails are key components of the biological nanomachinery involved in cell membrane puncturing, where they provide a means to deliver molecules and ions inside cells. Two intriguing examples of contractile tails are those from bacteriophage T4 and R2-pyocin. Although the two systems are different in terms of biological activity, they share a fascinatingly similar injection mechanism, during which the tail sheaths of both systems contract from a so-called extended state to around half of their length (the contracted state), accompanied by release of elastic energy originally stored in the sheath. Despite the great prevalence and biomedical importance of contractile delivery systems, many fundamental details of their injection machinery and dynamics are still unknown. In this work, we calculate the bending and torsional stiffness constants of a helical tail sheath strand of bacteriophage T4 and R2-pyocin, in both extended and contracted states, using molecular dynamics simulations of about one-sixth of the entire sheath. Differences in stiffness constants between the two systems are rationalized by comparing their all-atom monomer structures, changes in sheath architecture on contraction, and differences in interstrand interactions. The calculated coefficients indicate that the T4 strand is stiffer for both bending and torsion than the corresponding R2-pyocin strands in both extended and contracted conformations. The sheath strands also have greater stiffness in the contracted state for both systems. As the main application of this study, we describe how the stiffness constants can be incorporated in a model to simulate the dynamics of contractile nanoinjection machineries.

The B-DNA double helix can dynamically accommodate G-C and A-T base pairs in either Watson-Crick or Hoogsteen configurations. Here, we show that G-C(+) (in which + indicates protonation) and A-U Hoogsteen base pairs are strongly disfavored in A-RNA. As a result,N(1)-methyladenosine and N(1)-methylguanosine, which occur in DNA as a form of alkylation damage and in RNA as post-transcriptional modifications, have dramatically different consequences. Whereas they create G-C(+) and A-T Hoogsteen base pairs in duplex DNA, thereby maintaining the structural integrity of the double helix, they block base-pairing and induce local duplex melting in RNA. These observations provide a mechanism for disrupting RNA structure through post-transcriptional modifications. The different propensities to form Hoogsteen base pairs in B-DNA and A-RNA may help cells meet the opposing requirements of maintaining genome stability, on the one hand, and of dynamically modulating the structure of the epitranscriptome, on the other.

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.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.

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.

We study the conformational equilibrium between B-to-A forms of ds-DNA adsorbed onto a single-walled carbon nanotube (SWNT) using free energy profile calculations based on all-atom molecular dynamics simulations. The potential of mean force (PMF) of the B-to-A transition of ds-DNA in the presence of an uncharged (10,0) carbon nanotube for two dodecamers with poly-AT or poly-GC sequences is calculated as a function of a root-mean-square-distance (ΔRMSD) difference metric for the B-to-A transition. The calculations reveal that in the presence of a SWNT DNA favors B-form DNA significantly in both poly-GC and poly-AT sequences. Furthermore, the poly-AT DNA:SWNT complex shows a higher energy penalty for adopting an A-like conformation than poly-GC DNA:SWNT by several kcal/mol. The presence of a SWNT on either poly-AT or poly-GC DNA affects the PMF of the transition such that the B form is favored by as much as 10 kcal/mol. In agreement with published data, we find a potential energy minimum between A and B-form DNA at ΔRMSD ≈ −1.5 Å and that the presence of the SWNT moves this minimum by as much as ΔRMSD = 3 Å.

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.

The genetic material in living cells is organized into complex structures in which DNA is subjected to substantial contortions. Here we investigate the difference in structure, dynamics, and flexibility between two topological states of a short (107 base pair) DNA sequence in a linear form and a covalently closed, tightly curved circular DNA form. By employing a combination of all-atom molecular dynamics (MD) simulations and elastic rod modeling of DNA, which allows capturing microscopic details while monitoring the global dynamics, we demonstrate that in the highly curved regime the microscopic flexibility of the DNA drastically increases due to the local mobility of the duplex. By analyzing vibrational entropy and Lipari–Szabo NMR order parameters from the simulation data, we propose a novel model for the thermodynamic stability of high-curvature DNA states based on vibrational untightening of the duplex. This novel view of DNA bending provides a fundamental explanation that bridges the gap between classical models of DNA and experimental studies on DNA cyclization, which so far have been in substantial disagreement.

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.

The development of methods for predicting NMR chemical shifts with high accuracy and speed is increasingly allowing use of these abundant, readily accessible measurements in determining the structure and dynamics of proteins. For nucleic acids, however, despite the availability of semiempirical methods for predicting (1)H chemical shifts, their use in determining the structure and dynamics has not yet been examined. Here, we show that (1)H chemical shifts offer powerful restraints for RNA structure determination, allowing discrimination of native structure from non-native states to within 2-4 Å, and3 Å when highly flexible residues are ignored. Theoretical simulations shows that although (1)H chemical shifts can provide valuable information for constructing RNA dynamic ensembles, large uncertainties in the chemical shift predictions and inherent degeneracies lead to higher uncertainties as compared to residual dipolar couplings.

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.

Bacteriophage T4 is one of the most common and complex of the tailed viruses that infect host bacteria using an intriguing contractile tail assembly. Despite extensive progress in resolving the structure of T4, the dynamics of the injection machinery remains largely unknown. This paper contributes a first model of the injection machinery that is driven by elastic energy stored in a structure known as the sheath. The sheath is composed of helical strands of protein that suddenly collapse from an energetic, extended conformation prior to infection to a relaxed, contracted conformation during infection. We employ Kirchhoff rod theory to simulate the nonlinear dynamics of a single protein strand coupled to a model for the remainder of the virus, including the coupled translation and rotation of the head (capsid), neck and tail tube. Doing so provides an important building block towards the future goal of modeling the entire sheath structure which is composed of six interacting helical protein strands. The resulting numerical model exposes fundamental features of the injection machinery including the time scale and energetics of the infection process, the nonlinear conformational change experienced by the sheath, and the contribution of hydrodynamic drag on the head (capsid).

Enhanced sampling techniques are used to increase the frequency of “rare events” during computer simulations of complex molecules. Although methods exist that allow accurate thermodynamics to be recovered from enhanced simulations, recovering kinetics proves to be more challenging. Here we present an extrapolation approach that allows reliable kinetics to be recovered from potential-scaled MD simulations. The approach, based on Kramers’ rate theory, is simple and computationally efficient, and allows kinetics to be recovered without defining reaction coordinates. To test our approach, we use it to determine the kinetics of barrier crossing between two metastable states on the 2D-Muller potential and the C7eq to αR transition in alanine dipeptide. The mean first passage time estimates obtained are in excellent agreement with reference values obtained from direct simulations on the unscaled potentials performed over times that are orders of magnitude longer.

Bacteriophage T4 is one of the most common and complex of the tailed viruses that infect host bacteria using an intriguing contractile tail assembly. Despite extensive progress in resolving the structure of T4, the dynamics of the injection machinery remains largely unknown. This paper contributes a first model of the injection machinery that is driven by elastic energy stored in a structure known as the sheath. The sheath is composed of helical strands of protein that suddenly collapse from an energetic, extended conformation prior to infection to a relaxed, contracted conformation during infection. We employ Kirchhoff rod theory to simulate the nonlinear dynamics of a single protein strand coupled to a model for the remainder of the virus, including the coupled translation and rotation of the head (capsid), neck, and tail tube. Doing so provides an important building block toward the future goal of modeling the entire sheath structure which is composed of six interacting helical protein strands. The resulting numerical model exposes fundamental features of the injection machinery including the time scale and energetics of the infection process, the nonlinear conformational change experienced by the sheath, and the contribution of hydrodynamic drag on the head (capsid).

Sequence-specific DNA flexibility plays a key role in a variety of cellular interactions that are critical for gene packaging, expression, and regulation, yet few studies have experimentally explor...

Recent advances in the collection and interpretation of surface-sensitive vibrational spectroscopic measurements have made it possible to study the orientation of peptides and proteins in situ in a biologically relevant environment. However, interpretation of sum frequency generation (SFG) and attenuated total reflectance Fourier transform infrared (ATR-FTIR) vibrational spectroscopy is hindered by the fact that orientation cannot be inferred without some prior knowledge of the protein structure. In this work, molecular dynamics simulations were used to study the interfacial orientation and structural deformation of the short β-sheet peptide tachyplesin I at the polystyrene/water interface. By combining these results with ATR-FTIR and SFG measurements, reasonable agreement was found with the simulation results, suggesting that tachyplesin I lies parallel to the surface, although the simulation results imply a broader distribution of peptide twist angles than could be characterized using available experimental measurements. The interfacial structure was found to be deformable even when disulfide bonds were preserved, and these local deviations from a purely extended β-sheet conformation may be of importance to future developments in the interpretation of SFG and ATR-FTIR spectra.

Free-energy profiles describing the relative orientation of membrane proteins along predefined coordinates can be efficiently calculated by means of umbrella simulations. Such simulations generate reliable orientational distributions but are difficult to converge because of the very long equilibration times of the solvent and the lipid bilayer in explicit representation. Two implicit lipid membrane models are here applied in combination with the umbrella sampling strategy to the simulation of the transmembrane (TM) helical segment from virus protein U (Vpu). The models are used to study both orientation and energetics of this α-helical peptide as a function of hydrophobic mismatch. We observe that increasing the degree of positive hydrophobic mismatch increased the tilt angle of Vpu. These findings agree well with experimental data and as such validate the solvation models used in this study.

We combined molecular dynamics based free energy calculations with sum frequency generation (SFG) spectroscopy to study the orientational distribution of solvated peptides near hydrophobic surfaces. Using a simplified atomistic model of the polystyrene (PS) surface, molecular dynamics simulations have been applied to compute the orientational probability of an α-helical peptide, magainin 2, with respect to the PS/water interface. Free energy calculations revealed that the preferred (horizontal) peptide orientation was driven by the favorable interactions between the hydrophobic PS surface and the hydrophobic residues on the helix, and additional simulations examined the importance of small aggregate formation. Concentration-dependent measurements obtained via SFG vibrational spectroscopy suggest that, at very low peptide concentrations, magainin molecules tend to lie down at the PS/solution interface, which correlates well with the simulation results. When the concentration is increased, peptides exhibit behavior not captured by MD simulations using single helical peptides. A combination of simulations and experiments was shown to yield more reliable results with molecular-level insights into interaction between peptides and polymer surfaces.

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.

The mechanism of DNA compaction by dendrimers is key to the design of nanotechnologies that can deliver genetic material into cells. We present atomistic simulations, mesoscopic modeling and single-molecule pulling experiments describing DNA dendrimer interactions. All-atom molecular dynamics were used to characterize pulling-force-dependent interactions between DNA and generation-3 PAMAM amine-terminated dendrimers, and a free energy profile and mean forces along the interaction coordinate are calculated. The energy, force, and geometry parameters computed at the atomic level are input for a Monte Carlo model yielding mesoscopic force-extension curves. Actual experimental single-molecule curves obtained with optical tweezers are also presented, and they show remarkable agreement with the virtual curves from our model. The calculations reveal the microscopic origin of the hysteresis observed in the phase transition underlying compaction. A broad range of ionic and pulling parameters is sampled, and suggestions for windows of conditions to probe new single-molecule behavior are made.

We describe a strategy for constructing atomic resolution dynamical ensembles of RNA molecules, spanning up to millisecond timescales, that combines molecular dynamics (MD) simulations with NMR residual dipolar couplings (RDC) measured in elongated RNA. The ensembles are generated via a Monte Carlo procedure by selecting snap-shot from an MD trajectory that reproduce experimentally measured RDCs. Using this approach, we construct ensembles for two variants of the transactivation response element (TAR) containing three (HIV-1) and two (HIV-2) nucleotide bulges. The HIV-1 TAR ensemble reveals significant mobility in bulge residues C24 and U25 and to a lesser extent U23 and neighboring helical residue A22 that give rise to large amplitude spatially correlated twisting and bending helical motions. Omission of bulge residue C24 in HIV-2 TAR leads to a significant reduction in both the local mobility in and around the bulge and amplitude of inter-helical bending motions. In contrast, twisting motions of the helices remain comparable in amplitude to HIV-1 TAR and spatial correlations between them increase significantly. Comparison of the HIV-1 TAR dynamical ensemble and ligand bound TAR conformations reveals that several features of the binding pocket and global conformation are dynamically preformed, providing support for adaptive recognition via a ‘conformational selection’ type mechanism.

Third-generation (G3) poly(amidoamine) (PAMAM) dendrimers are simulated approaching 1,2-dimyristoyl- sn-glycero-3-phosphocholine (DMPC) bilayers with fully atomistic molecular dynamics, which enables the calculation of a free energy profile along the approach coordinate. Three different dendrimer terminations are examined: protonated primary amine, uncharged acetamide, and deprotonated carboxylic acid. As the dendrimer and lipids become closer, their attractive force increases (up to 240 pN) and the dendrimer becomes deformed as it interacts with the lipids. The total energy release upon binding of a G3-NH3+, G3-Ac, or G3-COO- dendrimer to a DMPC bilayer is, respectively, 36, 26, or 47 kcal/mol or, equivalently, 5.2, 3.2, or 4.7x10(-3) kcal/g. These results are analyzed in terms of the dendrimers' size, shape, and atomic distributions as well as proximity of individual lipid molecules and particular lipid atoms to the dendrimer. For example, an area of 9.6, 8.2, or 7.9 nm2 is covered on the bilayer for the G3-NH3+, G3-Ac, or G3-COO- dendrimers, respectively, while interacting strongly with 18-13 individual lipid molecules.

The molecular structures and enthalpy release during binding of poly(amidoamine) (PAMAM) dendrimers to 1,2-dimyristoyl-sn-glycero-3-phosphocholine (DMPC) bilayers were explored through atomistic molecular dynamics. Three PAMAM dendrimer terminations were examined: protonated primary amine, neutral acetamide, and deprotonated carboxylic acid. Fluid and gel lipid phases were examined to extract the effects of lipid tail mobility on the binding of generation-3 dendrimers, which are directly relevant to the nanoparticle interactions involving lipid rafts, endocytosis, lipid removal, and/or membrane pores. Upon binding to gel phase lipids, dendrimers remained spherical, had a constant radius of gyration, and approximately one-quarter of the terminal groups were in close proximity to the lipids. In contrast, upon binding to fluid phase bilayers, dendrimers flattened out with a large increase in their asphericity and radii of gyration. Although over twice as many dendrimer–lipid contacts were formed on fluid versus gel phase lipids, the dendrimer–lipid interaction energy was only 20% stronger. The greatest enthalpy release upon binding was between the charged dendrimers and the lipid bilayer. However, the stronger binding to fluid versus gel phase lipids was driven by the hydrophobic interactions between the inner dendrimer and lipid tails.

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.

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.

We investigate the electrostatics, energetics, and dynamics of dendrimer-DNA interactions that mimic protein-DNA complexes as a means to design facilitated mechanisms by which dendrimers can slide and search DNA for targets. By using all-atom molecular dynamics simulations, we calculated the free energy profiles of dendrimer-binding around the DNA via umbrella sampling. We also calculated electrostatic interaction maps in comparison to proteins, as well as the dynamical changes induced by DNA-dendrimer interactions via NMR-measurable order parameters. Our results show that for dendrimers to go around DNA, there is a free-energy barrier of 8.5 kcal/mol from the DNA major groove to DNA minor groove, with a minimum in the major groove. This barrier height makes it unlikely for an all-amine dendrimer to slide along DNA longitudinally, but following a helical path may be possible along the major groove. Comparison of the nonbonded interaction energy and the interaction free-energy profiles reveal a considerable entropic cost as the dendrimer binds to DNA. This is also supported by the mobility patterns obtained from NMR-measurable order parameter values, which show a decreased mobility of the dendrimer N-H bond vectors in the DNA-binding mode.

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 (MD) simulations, coarse-grained simulations and a statistical-mechanical model to study the initiation of DNA unzipping by an external force. The calculation of the potential of mean force profiles for the initial separation of the first few terminal base pairs in a DNA oligomer reveal that forces ranging between 130 and 230 pN are needed to disrupt the first base pair, values of an order of magnitude larger than those needed to disrupt base pairs in partially unzipped DNA. The force peak has an "echo," of approximately 50 pN, at the distance that unzips the second base pair. We show that the high peak needed to initiate unzipping derives from a free energy basin that is distinct from the basins of subsequent base pairs because of entropic contributions and we highlight the microscopic origin of the peak. Our results suggest a new window of exploration for single molecule experiments., Comment: 25 pages, 6 figures , Accepted for publication in Biophysical Journal

Thioredoxins (Trx) are a class of enzymes, which catalyze the reduction of disulphide bonds between two cysteine residues, commonly found in proteins. An experimental investigation into the reaction mechanisms employed by various species of Trx was carried out by Perez-Jimenez et al. using single molecule force-clamp spectroscopy. The experiment involved applying a pulling force along the disulfide bond of the protein substrate, and measuring the rate of the Trx-catalyzed reduction as a function of the pulling force. One interesting finding of the experiment was that some forms of thioredoxin exhibit a biphasic relationship for reduction rate as a function of force magnitude. For this project, a mathematical model of this system was created, which employs a Smoluchowski formalism in the vein of Agmon-Hopfied or Sumi-Markus models. The model describes the time evolution of the probability distribution function of the protein's configuration within a space defined as the internal protein coordinate, as it diffuses over a potential which distorts under the applied force, while losing probability density to the Bell model type “sink” term, which is also a function of the applied force, representing reactants going to products (disulphide bond cleavage). By numerically solving the Smoluchowski equation and integrating the resulting surface over both time and the protein coordinate to calculate lifetime for increasing values of applied force, the model successfully reproduced the experimentally observed values for disulphide bond reduction rate as a function of applied force. Parameterizing the Smoluchowski equation to fit the experimentally measured data points provided a means of drawing insights into a physical interpretation of the model including the relationship between degree of biphasic behavior and the distance along the reaction coordinate from the bottom of the reactant well to the top of the transition state.

Topoisomerases are enzymes of quintessence to the upkeep of superhelical DNA, and are vital for replication, transcription and recombination. An atomic-resolution model for human topoisomerase I in covalent complex with DNA is simulated using molecular dynamics with external potentials that mimic torque and bias the DNA duplex downstream of a single-strand cut to rotate around the intact strand, according to the prevailing enzymatic mechanism. The simulations reveal the first dynamical picture of how topoisomerase accommodates large-scale motion of DNA as it changes its supercoiling state, and indicate that relaxation of positive and negative supercoils are fundamentally different. To relax positive supercoils, two separate domains (the 'lips') of the protein open up by about 10-14 A, whereas to relax negative supercoils, a continuous loop connecting the upper and lower parts (and which was a hinge for opening the lips) stretches about 12 A while the lips remain unseparated. Normal mode analysis is additionally used to characterize the functional flexibility of the protein. Remarkably, the same combination of low-frequency eigenvectors exhibit the dominant contribution for both rotation mechanisms through a see-saw motion. The simulated mechanisms suggest mutations to control the relaxation of either type of supercoiling selectively and advance a hypothesis for the debated role of the N-terminal domain in supercoil relaxation.

Molecular dynamics simulations are presented for a Thermus aquaticus (Taq) DNA polymerase I complex (consisting of the protein, the primer-template DNA strands, and the incoming nucleotide) subjected to external forces. The results obtained with a force applied to the DNA template strand provide insights into the effect of the tension on the activity of the enzyme. At forces below 30pN a local model based on the parameters determined from the simulations, including the restricted motion of the DNA bases at the active site, yields a replication rate dependence on force in agreement with experiment. Simulations above 40pN reveal large conformational changes in the enzyme-bound DNA that may have a role in the force-induced exonucleolysis observed experimentally.