Search

While the linear Pearson correlation coefficient represents a well-established normalized measure to quantify the interrelation of two stochastic variables $X$ and $Y$, it fails for multidimensional variables such as Cartesian coordinates. Avoiding any assumption about the underlying data, the mutual information $I(X, Y)$ does account for multidimensional correlations. However, unlike the normalized Pearson correlation, it has no upper bound ($I \in [0, \infty)$), i.e., it is not clear if say, $I = 0.4$ corresponds to a low or a high correlation. Moreover, the mutual information (MI) involves the estimation of high-dimensional probability densities (e.g., six-dimensional for Cartesian coordinates), which requires a k-nearest neighbor algorithm, such as the estimator by Kraskov et al. [Phys. Rev. E 69, 066138 (2004)]. As existing methods to normalize the MI cannot be used in connection with this estimator, a new approach is presented, which uses an entropy estimation method that is invariant under variable transformations. The algorithm is numerically efficient and does not require more effort than the calculation of the (un-normalized) MI. After validating the method by applying it to various toy models, the normalized MI between the $C_{\alpha}$ -coordinates of T4 lysozyme is considered and compared to a correlation analysis of inter-residue contacts.

A novel mixed quantum-classical approach to simulating nonadiabatic dynamics of molecules at metal surfaces is presented. The method combines the numerically exact hierarchical equations of motion approach for the quantum electronic degrees of freedom with Langevin dynamics for the classical degrees of freedom, namely, low-frequency vibrational modes within the molecule. The approach extends previous mixed quantum-classical methods based on Langevin equations to models containing strong electron-electron or quantum electronic-vibrational interactions, while maintaining a nonperturbative and non-Markovian treatment of the molecule-metal coupling. To demonstrate the approach, nonequilibrium transport observables are calculated for a molecular nanojunction containing strong interactions.

The time-dependent relaxation of a dynamical system may exhibit a power-law behavior that is superimposed by log-periodic oscillations. Sornette [Phys. Rep. 297, 239 (1998)] showed that this behavior can be explained by a discrete scale invariance of the system, which is associated with discrete and equidistant timescales on a logarithmic scale. Examples include such diverse fields as financial crashes, random diffusion, and quantum topological materials. Recent time-resolved experiments and molecular dynamics simulations suggest that discrete scale invariance may also apply to hierarchical dynamics in proteins, where several fast local conformational changes are a prerequisite for a slow global transition to occur. Employing entropy-based timescale analysis and Markov state modeling to a simple one-dimensional hierarchical model and biomolecular simulation data, it is found that hierarchical systems quite generally give rise to logarithmically spaced discrete timescales. By introducing a one-dimensional reaction coordinate that collectively accounts for the hierarchically coupled degrees of freedom, the free energy landscape exhibits a characteristic staircase shape with two metastable end states, which causes the log-periodic time evolution of the system. The period of the log-oscillations reflects the effective roughness of the energy landscape, and can in simple cases be interpreted in terms of the barriers of the staircase landscape.

To sample rare events, dissipation-corrected targeted molecular dynamics (dcTMD) applies a constant velocity constraint along a one-dimensional reaction coordinate $s$, which drives an atomistic system from an initial state into a target state. Employing a cumulant approximation of Jarzynski's identity, the free energy $\Delta G (s)$ is calculated from the mean external work and dissipated work of the process. By calculating the friction coefficient $\Gamma (s)$ from the dissipated work, in a second step the equilibrium dynamics of the process can be studied by propagating a Langevin equation. While so far dcTMD has been mostly applied to study the unbinding of protein-ligand complexes, here its applicability to rare conformational transitions within a protein and the prediction of their kinetics is investigated. As this typically requires the introduction of multiple collective variables $\{x_j\}= \vec{x}$, a theoretical framework is outlined to calculate the associated free energy $\Delta G (\vec{x})$ and friction $\matrix{\Gamma}(\vec{x})$ from dcTMD simulations along coordinate $s$. Adopting the $\alpha$-$\beta$ transition of alanine dipeptide as well as the open-closed transition of T4 lysozyme as representative examples, the virtues and shortcomings of dcTMD to predict protein conformational transitions and the related kinetics are studied., Comment: This preprint is the unedited version of a manuscript that has been sent to a scientific journal for consideration as a publication and can be downloaded for private use only. Copyright with the journal and its publisher after publication

Adopting a $300 \, \mu$s-long molecular dynamics (MD) trajectory of the reversible folding of villin headpiece (HP35) published by D. E. Shaw Research, we recently constructed a Markov state model (MSM) of the folding process based on interresidue contacts [J. Chem. Theory Comput. 2023, ${\bf {19}}$, 3391]. The model reproduces the MD folding times of the system and predicts that both the native basin and the unfolded region of the free energy landscape are partitioned into several metastable substates that are structurally well characterized. Recognizing the need to establish well-defined but nontrivial benchmark problems, in this Perspective we study to what extent and in what sense this MSM may be employed as a reference model. To this end, we test the robustness of the MSM by comparing it to models that use alternative combinations of features, dimensionality reduction methods and clustering schemes. The study suggests some main characteristics of the folding of HP35, which should be reproduced by any other competitive model of the system. Moreover, the discussion reveals which parts of the MSM workflow matter most for the considered problem, and illustrates the promises and possible pitfalls of state-based models for the interpretation of biomolecular simulations., Comment: arXiv admin note: text overlap with arXiv:2303.03814

Markov state models represent a popular means to interpret molecular dynamics trajectories in terms of memoryless transitions between metastable conformational states. To provide a mechanistic understanding of the considered biomolecular process, these states should reflect structurally distinct conformations and ensure a timescale separation between fast intrastate and slow interstate dynamics. Adopting the folding of villin headpiece (HP35) as a well-established model problem, here we discuss the selection of suitable input coordinates or `features', such as backbone dihedral angles and interresidue distances. We show that dihedral angles account accurately for the structure of the native energy basin of HP35, while the unfolded region of the free energy landscape and the folding process are best described by tertiary contacts of the protein. To construct a contact-based model, we consider various ways to define and select contact distances, and introduce a low-pass filtering of the feature trajectory as well as a correlation-based characterization of states. Relying on input data that faithfully account for the mechanistic origin of the studied process, the states of the resulting Markov model are clearly discriminated by the features, describe consistently the hierarchical structure of the free energy landscape, and$\unicode{0x2014}$as a consequence$\unicode{0x2014}$correctly reproduce the slow timescales of the process.

Protein-ligand (un)binding simulations are a recent focus of biased molecular dynamics simulations. Such binding and unbinding can occur via different pathways in and out of a binding site. We here present a theoretical framework how to compute kinetics along separate paths and to combine the path-specific rates into global binding and unbinding rates for comparison with experiment. Using dissipation-corrected targeted molecular dynamics in combination with temperature-boosted Langevin equation simulations [Nat. Commun. \textbf{11}, 2918 (2020)] applied to a two-dimensional model and the trypsin-benzamidine complex as test systems, we assess the robustness of the procedure and discuss aspects of its practical applicability to predict multisecond kinetics of complex biomolecular systems., Comment: This preprint is the unedited version of a manuscript that has been published in J. Chen. Phys. and can be downloaded for private use only. Copyright with AIP and and the authors

While allostery is of paramount importance for protein signaling and regulation, the underlying dynamical process of allosteric communication is not well understood. PDZ3 domain represents a prime example of an allosteric single-domain protein, as it features a well-established long-range coupling between the C-terminal $\alpha_3$-helix and ligand binding. In an intriguing experiment, Hamm and coworkers employed photoswitching of the $\alpha_3$-helix to initiate a conformational change of PDZ3 that propagates from the C-terminus to the bound ligand within 200 ns. Performing extensive nonequilibrium molecular dynamics simulations, the modeling of the experiment reproduces the measured timescales and reveals a detailed picture of the allosteric communication in PDZ3. In particular, a correlation analysis identifies a network of contacts connecting the $\alpha_3$-helix and the core of the protein, which move in a concerted manner. Representing a one-step process and involving direct $\alpha_3$-ligand contacts, this cooperative transition is considered as elementary step in the propagation of conformational change.

Photoproteins such as bacteriorhodopsin (bR) and rhodopsin (Rho) need to effectively dissipate photoinduced excess energy to prevent their damage. Another well-studied G protein-coupled receptor (GPCR) containing 7 transmembrane (TM) helices is the B2 adrenergic receptor (B2AR), for which energy dissipation paths have been linked with allosteric communication. To study the vibrational energy transport in the active and inactive states of these GPCRs, a master equation approach [J. Chem. Phys. 152, 045103 (2020)] is employed, which uses scaling rules that allow to calculate energy transport rates solely based on the protein structure. Despite their structural similarities, the three GPCRs reveal quite different strategies to redistribute excess energy. While bR quickly removes the energy using the TM7 helix as a "lightning rod", Rho exhibits a rather poor energy dissipation, which might eventually require the hydrolysis of the Schiff base between the protein and the retinal chromophore to prevent overheating. Heating the ligand adrenaline of B2AR, the resulting energy transport network of the protein is found to change significantly upon switching from the active to the inactive state. While the energy flow may highlight aspects of the interresidue couplings of B2AR, it seems not particularly suited to explain allosteric phenomena.

To interpret molecular dynamics simulations of biomolecular systems, systematic dimensionality reduction methods are commonly employed. Among others, this includes principal component analysis (PCA) and time-lagged independent component analysis (TICA), which aim to maximize the variance and the timescale of the first components, respectively. A crucial first step of such an analysis is the identification of suitable and relevant input coordinates (the so-called features), such as backbone dihedral angles and interresidue distances. As typically only a small subset of those coordinates is involved in a specific biomolecular process, it is important to discard the remaining uncorrelated motions or weakly correlated noise coordinates. This is because they may exhibit large amplitudes or long timescales and therefore will be erroneously be considered important by PCA and TICA, respectively. To discriminate collective motions underlying functional dynamics from uncorrelated motions, the correlation matrix of the input coordinates is block-diagonalized by a clustering method. This strategy avoids possible bias due to presumed functional observables and conformational states or variation principles that maximize variance or timescales. Considering several linear and nonlinear correlation measures and various clustering algorithms, it is shown that the combination of linear correlation and the Leiden community detection algorithm yields excellent results for all considered model systems. These include the functional motion of T4 lysozyme to demonstrate the successful identification of collective motion, as well as the folding of villin headpiece to highlight the physical interpretation of the correlated motions in terms of a functional mechanism.

The friction coefficient of fluids may become a function of the velocity at increased external driving. This non-Newtonian behavior is of general theoretical interest as well as of great practical importance, e.g., for the design of lubricants. While the effect has been observed in large-scale atomistic simulations of bulk liquids, its theoretical formulation and microscopic origin is not well understood. Here we use dissipation-corrected targeted molecular dynamics, which pulls apart two tagged liquid molecules in the presence of surrounding molecules and analyzes this nonequilibrium process via a generalized Langevin equation. The approach is based on a second-order cumulant expansion of Jarzynski's identity, which is shown to be valid for fluids and therefore allows for an exact computation of the friction profile as well of the underlying memory kernel. We show that velocity-dependent friction in fluids results from an intricate interplay of near-order structural effects and the non-Markovian behavior of the friction memory kernel. For complex fluids such as the model lubricant \alkane, the memory kernel exhibits a stretched-exponential long-time decay, which reflects the multitude of timescales of the system., Comment: This preprint is the unedited version of a manuscript that has been submitted to a scientific journal for consideration as publication and may be downloaded for personal use only. Copyright will be with the publisher after publication

Given nonstationary data from molecular dynamics simulations, a Markovian Langevin model is constructed that aims to reproduce the time evolution of the underlying process. While at equilibrium the free energy landscape is sampled, nonequilibrium processes can be associated with a biased energy landscape, which accounts for finite sampling effects and external driving. Extending the data-driven Langevin equation (dLE) approach [Phys.\ Rev.\ Lett.\ {\bf 115}, 050602 (2015)] to the modeling of nonequilibrium processes, an efficient way to calculate multidimensional Langevin fields is outlined. The dLE is shown to correctly account for various nonequilibrium processes, including the enforced dissociation of sodium chloride in water, the pressure-jump induced nucleation of a liquid of hard spheres, and the conformational dynamics of a helical peptide sampled from nonstationary short trajectories., Comment: J. Phys. Chem. B 2021. This preprint is the unformatted and early version of a manuscript that has been published in an updated version in the Journal of Physical Chemistry B and may be downloaded for private use only. Copyright will with the journal / ACS

While allostery is of paramount importance for protein regulation, the underlying dynamical process of ligand (un)binding at one site, resulting time evolution of the protein structure, and change of the binding affinity at a remote site is not well understood. Here the ligand-induced conformational transition in a widely studied model system of allostery, the PDZ2 domain, is investigated by transient infrared spectroscopy accompanied by molecular dynamics simulations. To this end, an azobenzene derived photoswitch is linked to a peptide ligand in a way that its binding affinity to the PDZ2 domain changes upon switching, thus initiating an allosteric transition in the PDZ2 domain protein. The subsequent response of the protein, covering four decades of time ranging from $\sim$1~ns to $\sim$10~$\mu$s, can be rationalize by a remodelling of its rugged free energy landscape, with ver subtle shifts in the populations of a small number of structurally well defined states. It is proposed that structurally and dynamically driven allostery, often discussed as limiting scenarios of allosteric communication, actually go hand-in-hand, allowing the protein to adapt its free energy landscape to incoming signals., Comment: This unedited earlier version of the manuscript may be downloaded for personal use only. Any other use requires prior permission of the author and the National Academy of Sciences USA. The final manuscript was published in Proceedings of the National Academy of Sciences USA 117, 26031-26039 (2020) and can be found under https://www.pnas.org/content/117/42/26031.short

Coarse-graining of fully atomistic molecular dynamics simulations is a long-standing goal in order to allow the description of processes occurring on biologically relevant timescales. For example, the prediction of pathways, rates and rate-limiting steps in protein-ligand unbinding is crucial for modern drug discovery. To achieve the enhanced sampling, we first perform dissipation-corrected targeted molecular dynamics simulations, which yield free energy and friction profiles of the molecular process under consideration. In a second step, we use these fields to perform temperature-boosted Langevin simulations which account for the desired molecular kinetics occurring on multisecond timescales and beyond. Adopting the dissociation of solvated sodium chloride as well as trypsin-benzamidine and Hsp90-inhibitor protein-ligand complexes as test problems, we are able to reproduce rates from molecular dynamics simulation and experiments within a factor of 2-20, and dissociation constants within a factor of 1-4. Analysis of the friction profiles reveals that binding and unbinding dynamics are mediated by changes of the surrounding hydration shells in all investigated systems., Comment: This unedited earlier version of the manuscript may be downloaded for personal use only. The final manuscript was published in Nature Communications 11, 2918 (2020) as open access publication and is available at https://www.nature.com/articles/s41467-020-16655-1

Principal component analysis (PCA) represents a standard approach to identify collective variables $\{x_i\}\!=\!\boldsymbol{x}$, which can be used to construct the free energy landscape $\Delta G(\boldsymbol{x})$ of a molecular system. While PCA is routinely applied to equilibrium molecular dynamics (MD) simulations, it is less obvious how to extend the approach to nonequilibrium simulation techniques. This includes, e.g., the definition of the statistical averages employed in PCA, as well as the relation between the equilibrium free energy landscape $\Delta G(\boldsymbol{x})$ and energy landscapes $\Delta{\cal G} (\boldsymbol{x})$ obtained from nonequilibrium MD. As an example for a nonequilibrium method, `targeted MD' is considered which employs a moving distance constraint to enforce rare transitions along some biasing coordinate $s$. The introduced bias can be described by a weighting function $P(s)$, which provides a direct relation between equilibrium and nonequilibrium data, and thus establishes a well-defined way to perform PCA on nonequilibrium data. While the resulting distribution ${\cal P}(\boldsymbol{x})$ and energy $\Delta{\cal G} \propto \ln {\cal P}$ will not reflect the equilibrium state of the system, the nonequilibrium energy landscape $\Delta{\cal G} (\boldsymbol{x})$ may directly reveal the molecular reaction mechanism. Applied to targeted MD simulations of the unfolding of decaalanine, for example, a PCA performed on backbone dihedral angles is shown to discriminate several unfolding pathways. Although the formulation is in principle exact, its practical use depends critically on the choice of the biasing coordinate $s$, which should account for a naturally occurring motion between two well-defined end-states of the system., Comment: This article may be downloaded for personal use only. Any other use requires prior permission of the author and AIP Publishing. This article appeared in J. Chem. Phys., 150(20), 204110 and may be found at https://aip.scitation.org/doi/10.1063/1.5089636

Across many fields, a problem of interest is to predict the transition rates between nodes of a network, given limited stationary state and dynamical information. We give a solution using the principle of Maximum Caliber. We find the transition rate matrix by maximizing the path entropy of a random walker on the network constrained to reproducing a stationary distribution and a few dynamical averages. A main finding here is that when constrained only by the mean jump rate, the rate matrix is given by a square-root dependence of the rate, $\omega_{ab} \propto \sqrt{p_b/p_a}$, on $p_a$ and $p_b$, the stationary state populations at nodes a and b. We give two examples of our approach. First, we show that this method correctly predicts the correlated rates in a biochemical network of two genes, where we know the exact results from prior simulation. Second, we show that it correctly predicts rates of peptide conformational transitions, when compared to molecular dynamics simulations. This method can be used to infer large numbers of rates on known networks where smaller numbers of steady-state node populations are known.

A computational approach to describe the energy relaxation of a high-frequency vibrational mode in a fluctuating heterogeneous environment is outlined. Extending previous work [H. Fujisaki, Y. Zhang, and J.E. Straub, J. Chem. Phys. {\bf 124}, 144910 (2006)], second-order time-dependent perturbation theory is employed which includes the fluctuations of the parameters in the Hamiltonian within the vibrational adiabatic approximation. This means that the time-dependent vibrational frequencies along an MD trajectory are obtained via a partial geometry optimization of the solute with fixed solvent and a subsequent normal mode calculation. Adopting the amide I mode of N-methylacetamide in heavy water as a test problem, it is shown that the inclusion of dynamic fluctuations may significantly change the vibrational energy relaxation. In particular, it is found that relaxation occurs in two phases, because for short times ($\lesssim$ 200 fs) the spectral density appears continuous due to the frequency-time uncertainty relation, while at longer times the discrete nature of the bath becomes apparent. Considering the excellent agreement between theory and experiment, it is speculated if this behavior can explain the experimentally obtained biphasic relaxation the amide I mode of N-methylacetamide., Comment: 24 pages, 7 figures, submitted to J. Chem. Phys

Employing extensive quantum-chemical calculations at the DFT/B3LYP and MP2 level, a quartic force field of isolated N-methylacetamide is constructed. Taking into account 24 vibrational degrees of freedom, the model is employed to perform numerically exact vibrational configuration interaction calculations of the vibrational energy relaxation of the amide I mode. It is found that the energy transfer pathways may sensitively depend on details of the theoretical description. Moreover, the exact reference calculations were used to study the applicability and accuracy of (i) the quasiclassical trajectory method, (ii) time-dependent second-order perturbation theory, and (iii) the instantaneous normal mode description of frequency fluctuations. Based on the results, several strategies to describe vibrational energy relaxation in biomolecular systems are discussed., Comment: 18 pages, 6 figures, submitted to J. Phys. Chem. A