Filtering molecular dynamics trajectories to reveal low-frequency collective motions: phospholipase A2

A novel method for analysing molecular dynamics trajectories has been developed, which filters out high frequencies using digital signal processing techniques and facilitates focusing on the low-frequency collective motions of proteins. These motions involve low energy slow motions, which lead to important biological phenomena such as domain closure and allosteric effects in enzymes. The filtering method treats each of the atomic trajectories obtained from the molecular dynamics simulation as a "signal". The trajectories of each of the atoms in the system (or any subset of interest) are Fourier transformed to the frequency domain, a filtering function is applied and then an inverse transformation back to the time domain yields the filtered trajectory. The filtering method has been used to study the dynamics of the enzyme phospholipase A2. In the filtered trajectory, all the high frequency bond and valence angle vibrations were eliminated, leaving only low-frequency motion, mainly fluctuations in torsions and conformational transitions. Analysis of this trajectory revealed interesting motions of the protein, including concerted movements of helices, and changes in shape of the active site cavity. Unlike normal mode analysis, which has been used to study the motion of proteins, this method does not require converged minimizations or diagonalization of a matrix of second derivatives. In addition, anharmonicity, multiple minima and conformational transitions are treated explicitly. Thus, the filtering method avoids most of the approximations implicit in other investigations of the dynamic behaviour of large systems.     

Extraction of the energetics of selected types of motion from molecular dynamics trajectories by filtering

A novel method for analyzing molecular dynamics trajectories has been developed which enables the study of selected motions and the corresponding energetics. In particular, it is possible to filter out the high-frequency motions and focus on the structural and energetic features of low-frequency collective motions. The trajectories of the properties of interest are Fourier transformed to the frequency domain, a filtering function is applied, and then an inverse transformation back to the time domain yields the filtered trajectory. The method is demonstrated for harmonic fluctuations and conformational transitions of acetamide and N-acetylalanine N-methylamide, as models for peptides and proteins.     

A milestoning study of the kinetics of an allosteric transition: atomically detailed simulations of deoxy Scapharca hemoglobin

Atomically detailed simulations are used to compute the kinetics of the R-to-T transition in deoxy Scapharca hemoglobin. A computational approach called milestoning is utilized that combines 1), an efficient reaction path algorithm; and 2), a "fragment and glue" approach for classical trajectories. Milestoning computes the R-to-T transition kinetics on the microsecond timescale based on atomically detailed trajectories that rarely exceed a nanosecond. Eleven reference hypersurfaces (milestones) are constructed along the reaction coordinate, which is computed with a global path optimization algorithm. Two-hundred classical trajectories are calculated for each of the milestones to collect local distributions of first passage times. These local distributions are used in a non-Markovian theory to compute the overall timescale. Exponential enrichment of reactive trajectories, an important component of the milestoning approach, makes these calculations possible. The overall timescale of the reaction is estimated as 10 +/- 9 micros, in accord with available experimental data. The barrier is not sharp and is spread over four milestones. Even after the most significant structural changes are completed (phenylalanine F4 ring flips), highly collective and activated motions continue. The calculations suggest an additional late free energy barrier.     

Reduction of intrinsic kinetic and thermodynamic barriers for enzyme-catalysed proton transfers from carbon acid substrates

Many enzymes catalyse the heterolytic abstraction of the alpha-proton from a carbon acid substrate. Gerlt and Gassman have applied Marcus formalism to such proton transfer reactions to argue that transition states for concerted general acid-general base catalysed enolization at enzyme active sites occur late on the reaction coordinate (J. Am. Chem. Soc. 115 (1993) 11552). We postulate that as an enzyme evolves, it may decrease deltaG++ for a proton transfer step associated with substrate enolization by following the path of steepest descent on the two-dimensional surface corresponding to deltaG++, as defined by Marcus formalism. We show that for an enzyme that has decreased deltaG++ following the path of steepest descent, the values of the intrinsic kinetic (deltaG++(int,E)) and thermodynamic (deltaG(E)0) barriers for proton transfer reactions on the enzyme may be predicted from the known values of deltaG++(int,N) and deltaG(N)0 for the corresponding non-enzymic reaction and the free energy of activation on the enzyme (deltaG++(E)). In addition, the enzymic transition state will occur later on the reaction coordinate than the corresponding non-enzymic transition state (i.e. x++(E)>x++(N)) if the condition (6 - square root 2)/82deltaG++(int,N).     

Multiple time step diffusive Langevin dynamics for proteins

We present an algorithm for simulating the long time scale dynamics of proteins and other macromolecules. Our method applies the concept of multiple time step integration to the diffusive Langevin equation, in which short time scale dynamics are replaced by friction and noise. The macromolecular force field is represented at atomic resolution. Slow motions are modeled by constrained Langevin dynamics with very large time steps, while faster degrees of freedom are kept in local thermal equilibrium. In the limit of a sufficiently large molecule, our algorithm is shown to reduce the CPU time required by two orders of magnitude. We test the algorithm on two systems, alanine dipeptide and bovine pancreatic trypsin inhibitor (BPTI), and find that it accurately calculates a variety of equilibrium and dynamical properties. In the case of BPTI, the CPU time required is reduced by nearly a factor of 60 compared to a conventional, unconstrained Langevin simulation using the same force field. Proteins 30:215–227, 1998. © 1998 Wiley-Liss, Inc.     

Multiple Time Step Brownian Dynamics for Long Time Simulation of Biomolecules

We report a multiple time step algorithm applied to an atomistic Brownian dynamics simulation for simulating the long time scale dynamics of biomolecules. The algorithm was based on the original multiple time step method; a short time step was used to keep faster motions in local equilibrium. When applied to a 28-mer # # ! folded peptide, the simulation gave stable trajectories and the computation time was reduced by a factor of 160 compared to a conventional molecular dynamics simulation using explicit water molecules. We applied it for the folding simulation of a 13-mer ! -helical peptide, giving a successful folding simulation. These results indicate that the Brownian dynamics with the multiple time step algorithm is useful for studies of biomolecular motions by long time simulation.     

Enzymatic circularization of a malto-octaose linear chain studied by stochastic reaction path calculations on cyclodextrin glycosyltransferase

Cyclodextrin glycosyltransferase (CGTase) is an enzyme belonging to the alpha-amylase family that forms cyclodextrins (circularly linked oligosaccharides) from starch. X-ray work has indicated that this cyclization reaction of CGTase involves a 23-A movement of the nonreducing end of a linear malto-oligosaccharide from a remote binding position into the enzyme acceptor site. We have studied the dynamics of this sugar chain circularization through reaction path calculations. We used the new method of the stochastic path, which is based on path integral theory, to compute an approximate molecular dynamics trajectory of the large (75-kDa) CGTase from Bacillus circulans strain 251 on a millisecond time scale. The result was checked for consistency with site-directed mutagenesis data. The combined data show how aromatic residues and a hydrophobic cavity at the surface of CGTase actively catalyze the sugar chain movement. Therefore, by using approximate trajectories, reaction path calculations can give a unique insight into the dynamics of complex enzyme reactions.     

Picosecond timescale rigid-helix and side-chain motions in deoxymyoglobin

The contribution of rigidbody motions to the atomic trajectories in a 100 ps molecular dynamics simulation of deoxymyoglobin is examined. Two typesof rigid-body motions are considered: one in which the helices are rigid units and one in which the side-chains are rigid units. Using a quaternionbased algorithm, fits of the rigid reference structures are made to each time frame of the simulation to derive trajectories of the rigid-body motions. The fitted trajectories are analysed in terms of atomic position fluctuations, mean-square displacements as a function of time, velocity autocorrelation functions and densities of states. The results are compared with the corresponding quantities calculated from the full trajectory. The relative contribution of the rigid helix motions to the helix atom dynamics depends on which quantity is examined and on which subset of atoms is chosen: rigid-helix motions contribute 86% of the rms helix backbone atomic position fluctuations, but 30% of the helix,: atom (backbone and side-chain) mean square displacements and only 1.1% of total kinetic energy. Only very low-frequency motions contribute to the rigid-helix dynamics; the rigid-body analysis allows characteristic rigid-helix vibrations to be identified and described. Treating the side-chains as rigid bodies is foundto be an excellent approximation to both their diffusive and vibrationalmean-square displacements: 96% of side-chain atom mean-square displacements originate from rigid side-Chain motions. However, the errors in theside-chain atomic positional fits are not always small. An analysis is madeof factors contributing to the positional error for different types of side-chain. © Wiley-Liss, Inc.     

Leveraging nonlinear dynamic models to predict progression of neuroimaging biomarkers

Using biomarkers to model disease course effectively and make early prediction is a challenging but critical path to improving diagnostic accuracy and designing preventive trials for neurological disorders. Leveraging the domain knowledge that certain neuroimaging biomarkers may reflect the disease pathology, we propose a model inspired by the neural mass model from cognitive neuroscience to jointly model nonlinear dynamic trajectories of the biomarkers. Under a nonlinear mixed‐effects model framework, we introduce subject‐ and biomarker‐specific random inflection points to characterize the critical time of underlying disease progression as reflected in the biomarkers. A latent liability score is shared across biomarkers to pool information. Our model allows assessing how the underlying disease progression will affect the trajectories of the biomarkers, and, thus, is potentially useful for individual disease management or preventive therapeutics. We propose an EM algorithm for maximum likelihood estimation, where in the E step, a normal approximation is used to facilitate numerical integration. We perform extensive simulation studies and apply the method to analyze data from a large multisite natural history study of Huntington's Disease (HD). The results show that some neuroimaging biomarker inflection points are early signs of the HD onset. Finally, we develop an online tool to provide the individual prediction of the biomarker trajectories given the medical history and baseline measurements.     

Approximate quantum trajectory dynamics for reactive processes in condensed phase

A method of molecular dynamics with quantum corrections, practical for studies of large molecular systems, is reviewed. The approach is based on the Bohmian formulation of the time-dependent Schrödinger equation in which a wavefunction is represented by an ensemble of interdependent trajectories. The quantum effects come from the quantum potential acting on trajectories on par with the usual classical potential. The quantum potential is determined from the evolving nuclear wavefunction, i.e. from the quantum trajectory (QT) ensemble itself. For practical and conceptual reasons the quantum potential and corresponding quantum nuclear effect are computed only for the selected light nuclei. For studies of reactive chemical processes, the classical potential is computed on-the-fly using the density functional tight binding method of electronic structure. A massively parallel implementation, based on the message passing interface allows for efficient simulations of ensembles of thousands of trajectories describing systems of up to 200 atoms. As a biochemical application, the approximate QT approach is used to model the tunnelling-dominated proton transfer in soybean-lipoxygenase-1. A materials science application is represented by a study of the nuclear quantum effect on adsorption of hydrogen and deuterium on a C₃₇H₁₅ molecule, which is a model ‘flake’ of graphene.     

Real-time interactive frequency filtering of molecular dynamics trajectories

Molecular dynamics simulations of atomic motion in protein and nucleic acid molecules must be done on a femtosecond time-scale. Much of this rapid motion is unimportant for the slower changes that are most relevant to biological function (conformational changes, substrate binding, protein folding). The high-frequency motion makes simulations computationally expensive. More importantly, the high frequencies obscure visualization of the relevant dynamics processes. Sessions, Dauber-Osguthorpe and Osguthorpe presented a method for removing high-frequency motions from atomic co-ordinates of trajectories generated by simulation. While that study used fast Fourier methods and emphasized the use of filtering for analysis of trajectories, this communication describes a new method that makes it much easier to use frequency filtering in programs that display trajectories as a sequence of moving images. Tests of the method on systems extending from pure water to proteins and nucleic acid molecules in vacuo and in solution have demonstrated its general utility. Impressed with the power and simplicity of the new method, we wish to present it in sufficient detail to allow others to implement it themselves.     

Two-arm trajectory planning in a manipulation task

This paper addresses the trajectory planning problem for a task which requires positioning and orienting an object firmly grasped by two hands at a visually specified goal configuration in the horizontal plane. The motor task involves three degrees of freedom (two translational and one rotational), and the motions of the arms are constrained by the physical coupling through the held object. Experimentally measured trajectories of two arms in the coordinated positioning/orienting task are presented. The hypothesis that the rotational and translational components of motions are decoupled and independently planned is tested. Two explicit mathematical models to account for the kinematic features of the two-arm motions are formulated, and the predictions of the models are compared with the experimental data. Both models extend the minimum-jerk model to the two-arm coordinated motions case. The trajectories predicted by the models were found to be in qualitative agreement with the experimental data. However, neither model could account for the observed configuration dependence of the motions, nor for some of the properties of the measured velocity components of the motions. Our findings support the idea that the rotational and translational components of two-arm motions in the positioning/orienting task are independently planned in extra-personal space, and are further combined in a hierarchical fashion to produce the observed motions. The tested models may serve as a basis for further investigations of issues pertinent to the generation of two-arm trajectories. Received: 27 March 1995 / Accepted in revised form: 17 September 1996     

Analysis of kinematic invariances of multijoint reaching movement

 There is a no unique relationship between the trajectory of the hand, represented in cartesian or extrinsic space, and its trajectory in joint angle or intrinsic space in the general condition of joint redundancy. The goal of this work is to analyze the relation between planning the trajectory of a multijoint movement in these two coordinate systems. We show that the cartesian trajectory can be planned based on the task parameters (target coordinates, etc.) prior to and independently of angular trajectories. Angular time profiles are calculated from the cartesian trajectory to serve as a basis for muscle control commands. A unified differential equation that allows planning trajectories in cartesian and angular spaces simultaneously is proposed. Due to joint redundancy, each cartesian trajectory corresponds to a family of angular trajectories which can account for the substantial variability of the latter. A set of strategies for multijoint motor control following from this model is considered; one of them coincides with the frog wiping reflex model and resolves the kinematic inverse problem without inversion. The model trajectories exhibit certain properties observed in human multijoint reaching movements such as movement equifinality, straight end-point paths, bell-shaped tangential velocity profiles, speed-sensitive and speed-insensitive movement strategies, peculiarities of the response to double-step targets, and variations of angular trajectory without variations of the limb end-point trajectory in cartesian space. In humans, those properties are almost independent of limb configuration, target location, movement duration, and load. In the model, these properties are invariant to an affine transform of cartesian space. This implies that these properties are not a special goal of the motor control system but emerge from movement kinematics that reflect limb geometry, dynamics, and elementary principles of motor control used in planning. All the results are given analytically and, in order to compare the model with experimental results, by computer simulations. Received: 6 April 1994/Accepted in revised form: 25 April 1995     

Mechanisms underlying the generation of averaged modified trajectories

In this work we have studied what mechanisms might possibly underlie arm trajectory modification when reaching toward visual targets. The double-step target displacement paradigm was used with inter-stimulus intervals (ISIs) in the range of 10-300 ms. For short ISIs, a high percentage of the movements were found to be initially directed in between the first and second target locations (averaged trajectories). The initial direction of motion was found to depend on the target configuration, and on : the time difference between the presentation of the second stimulus and movement onset. To account for the kinematic features of the averaged trajectories two modification schemes were compared: the superposition scheme and the abort-replan scheme. According to the superposition scheme, the modified trajectories result from the vectorial addition of two elemental motions: one for moving between the initial hand position and an intermediate location, and a second one for moving between that intermediate location and the final target. According to the abort-replan scheme, the initial plan for moving toward the intermediate location is aborted and smoothly replaced by a new plan for moving from the hand position at the time the trajectory is modified to the final target location. In both tested schemes we hypothesized that due to the quick displacement of the stimulus, the internally specified intermediate goal might be influenced by both stimuli and may be different from the location of the first stimulus. It was found that the statistically most successful model in accounting for the measured data is based on the superposition scheme. It is suggested that superposition of simple independent elemental motions might be a general principle for the generation of modified motions, which allows for efficient, parallel planning. For increasing values of the inferred locations of the intermediate targets were found to gradually shift from the first toward the second target locations along a path that curved toward the initial hand position. These inferred locations show a strong resemblance to the intermediate locations of saccades generated in a similar double-step paradigm. These similarities in the specification of target locations used in the generation of eye and hand movements may serve to simplify visuomotor integration. Received: 22 June 1994 / Accepted in revised form: 15 September 1994     

Locally accessible conformations of proteins: multiple molecular dynamics simulations of crambin.

Multiple molecular dynamics (MD) simulations of crambin with different initial atomic velocities are used to sample conformations in the vicinity of the native structure. Individual trajectories of length up to 5 ns sample only a fraction of the conformational distribution generated by ten independent 120 ps trajectories at 300 K. The backbone atom conformational space distribution is analyzed using principal components analysis (PCA). Four different major conformational regions are found. In general, a trajectory samples only one region and few transitions between the regions are observed. Consequently, the averages of structural and dynamic properties over the ten trajectories differ significantly from those obtained from individual trajectories. The nature of the conformational sampling has important consequences for the utilization of MD simulations for a wide range of problems, such as comparisons with X-ray or NMR data. The overall average structure is significantly closer to the X-ray structure than any of the individual trajectory average structures. The high frequency (less than 10 ps) atomic fluctuations from the ten trajectories tend to be similar, but the lower frequency (100 ps) motions are different. To improve conformational sampling in molecular dynamics simulations of proteins, as in nucleic acids, multiple trajectories with different initial conditions should be used rather than a single long trajectory.     

Sequential estimation for prescribed statistical accuracy in stochastic simulation of biological systems

Stochastic simulation of biological systems proceeds by repeatedly generating sample paths or trajectories of the underlying stochastic process, from which many relevant and important system properties can be obtained. While a great deal of research is targeted towards accelerated trajectory generation, issues concerned with the variability across trajectories are often neglected. Advanced methods for properly quantifying the statistical accuracy and determining a reasonable number of trajectories are hardly addressed formally in the context of biological system simulation, though mathematical statistics provides a large body of powerful theory. We invoke this theory and show how mathematically well-founded sequential estimation approaches serve for systematically generating enough but not too many trajectories for achieving a certain prescribed accuracy. The practical applicability is demonstrated and illustrated by numerical examples through simulation studies of an immigration-death process and a gene regulatory network.     

An analysis of the influence of protein intrinsic dynamical properties on its thermal unfolding behavior

The influence of the protein topology-encoded dynamical properties on its thermal unfolding motions was studied in the present work. The intrinsic dynamics of protein topology was obtained by the anisotropic network model (ANM). The ANM has been largely used to investigate protein collective functional motions, but it is not well elucidated if this model can also reveal the preferred large-scale motions during protein unfolding. A small protein barnase is used as a typical case study to explore the relationship between protein topology-encoded dynamics and its unfolding motions. Three thermal unfolding simulations at 500 K were performed for barnase and the entire unfolding trajectories were sampled and partitioned into several windows. For each window, the preferred unfolding motions were investigated by essential dynamics analysis, and then associated with the intrinsic dynamical properties of the starting conformation in this window, which is detected by ANM. The results show that only a few slow normal modes imposed by protein structure are sufficient to give a significant overlap with the preferred unfolding motions. Especially, the large amplitude unfolding movements, which imply that the protein jumps out of a local energy basin, can be well described by a single or several ANM slow modes. Besides the global motions, it is also found that the local residual fluctuations encoded in protein structure are highly correlated with those in the protein unfolding process. Furthermore, we also investigated the relationship between protein intrinsic flexibility and its unfolding events. The results show that the intrinsic flexible regions tend to unfold early. Several early unfolding events can be predicted by analysis of protein structural flexibility. These results imply that protein structure-encoded dynamical properties have significant influences on protein unfolding motions.     

TransPath: a computational method for locating ion transit pathways through membrane proteins

The finely tuned structures of membrane channel proteins allow selective passage of ions through the available aqueous pores. To understand channel function, it is crucial to locate the pores and study their physical and chemical properties. Here, we propose a new pore-searching algorithm (TransPath), which uses the Configurational Bias Monte Carlo (CBMC) method to generate transmembrane trajectories driven by both geometric and electrostatic features. The trajectories are binned into groups determined by a vector distance criterion. From each group, a representative trajectory is selected based on the Rosenbluth weight, and the geometrically optimal path is obtained by simulated annealing. Candidate ion pathways can then be determined by analysis of the radius and potential profiles. The proposed method and its implementation are illustrated using the bacterial KcsA potassium channel as an example. The procedure is then applied to the more complex structures of the bacterial E. coli chloride channel homolog and a homology model of the ClC-0 channel.     

Motions and structural variability within toxins: implication for their use as scaffolds for protein engineering

Animal toxins are small proteins built on the basis of a few disulfide bonded frameworks. Because of their high variability in sequence and biologic function, these proteins are now used as templates for protein engineering. Here we report the extensive characterization of the structure and dynamics of two toxin folds, the "three-finger" fold and the short alpha/beta scorpion fold found in snake and scorpion venoms, respectively. These two folds have a very different architecture; the short alpha/beta scorpion fold is highly compact, whereas the "three-finger" fold is a beta structure presenting large flexible loops. First, the crystal structure of the snake toxin alpha was solved at 1.8-A resolution. Then, long molecular dynamics simulations (10 ns) in water boxes of the snake toxin alpha and the scorpion charybdotoxin were performed, starting either from the crystal or the solution structure. For both proteins, the crystal structure is stabilized by more hydrogen bonds than the solution structure, and the trajectory starting from the X-ray structure is more stable than the trajectory started from the NMR structure. The trajectories started from the X-ray structure are in agreement with the experimental NMR and X-ray data about the protein dynamics. Both proteins exhibit fast motions with an amplitude correlated to their secondary structure. In contrast, slower motions are essentially only observed in toxin alpha. The regions submitted to rare motions during the simulations are those that exhibit millisecond time-scale motions. Lastly, the structural variations within each fold family are described. The localization and the amplitude of these variations suggest that the regions presenting large-scale motions should be those tolerant to large insertions or deletions.