Correction: Path Similarity Analysis: A Method for Quantifying Macromolecular Pathways.

[This corrects the article DOI: 10.1371/journal.pcbi.1004568.].

A broader view on jamming: from spring networks to circle packings.

Jamming occurs when objects like grains are packed tightly together (e.g. grain silos). It is highly cooperative and can lead to phenomena like earthquakes, traffic jams, etc. In this paper we point out the paramount importance of the underlying contact network for jammed systems; the network must have one contact in excess of isostaticity and a finite bulk modulus. Isostatic means that the number of degrees of freedom is exactly balanced by the number of constraints. This defines a large class of networks that can be constructed without the necessity of packing particles together compressively (either in the lab or computationally). One such construction, which we explore here, involves setting up the Delaunay triangulation of a Poisson disk sampling and then removing edges to maximize the bulk modulus until the isostatic plus one edge is reached. This construction works in any dimensions and here we give results in 2D where we also show how such networks can be transformed into disk packs.

Ring correlations in random networks.

We examine the correlations between rings in random network glasses in two dimensions as a function of their separation. Initially, we use the topological separation (measured by the number of intervening rings), but this leads to pseudo-long-range correlations due to a lack of topological charge neutrality in the shells surrounding a central ring. This effect is associated with the noncircular nature of the shells. It is, therefore, necessary to use the geometrical distance between ring centers. Hence we find a generalization of the Aboav-Weaire law out to larger distances, with the correlations between rings decaying away when two rings are more than about three rings apart.

Anchored boundary conditions for locally isostatic networks.

Finite pieces of locally isostatic networks have a large number of floppy modes because of missing constraints at the surface. Here we show that by imposing suitable boundary conditions at the surface the network can be rendered effectively isostatic. We refer to these as anchored boundary conditions. An important example is formed by a two-dimensional network of corner sharing triangles, which is the focus of this paper. Another way of rendering such networks isostatic is by adding an external wire along which all unpinned vertices can slide (sliding boundary conditions). This approach also allows for the incorporation of boundaries associated with internal holes and complex sample geometries, which are illustrated with examples. The recent synthesis of bilayers of vitreous silica has provided impetus for this work. Experimental results from the imaging of finite pieces at the atomic level need such boundary conditions, if the observed structure is to be computer refined so that the interior atoms have the perception of being in an infinite isostatic environment.

Partial unfolding and refolding for structure refinement: A unified approach of geometric simulations and molecular dynamics.

The most successful protein structure prediction methods to date have been template-based modeling (TBM) or homology modeling, which predicts protein structure based on experimental structures. These high accuracy predictions sometimes retain structural errors due to incorrect templates or a lack of accurate templates in the case of low sequence similarity, making these structures inadequate in drug-design studies or molecular dynamics simulations. We have developed a new physics based approach to the protein refinement problem by mimicking the mechanism of chaperons that rehabilitate misfolded proteins. The template structure is unfolded by selectively (targeted) pulling on different portions of the protein using the geometric based technique FRODA, and then refolded using hierarchically restrained replica exchange molecular dynamics simulations (hr-REMD). FRODA unfolding is used to create a diverse set of topologies for surveying near native-like structures from a template and to provide a set of persistent contacts to be employed during re-folding. We have tested our approach on 13 previous CASP targets and observed that this method of folding an ensemble of partially unfolded structures, through the hierarchical addition of contact restraints (that is, first local and then nonlocal interactions), leads to a refolding of the structure along with refinement in most cases (12/13). Although this approach yields refined models through advancement in sampling, the task of blind selection of the best refined models still needs to be solved. Overall, the method can be useful for improved sampling for low resolution models where certain of the portions of the structure are incorrectly modeled.

Path Similarity Analysis: A Method for Quantifying Macromolecular Pathways.

Diverse classes of proteins function through large-scale conformational changes and various sophisticated computational algorithms have been proposed to enhance sampling of these macromolecular transition paths. Because such paths are curves in a high-dimensional space, it has been difficult to quantitatively compare multiple paths, a necessary prerequisite to, for instance, assess the quality of different algorithms. We introduce a method named Path Similarity Analysis (PSA) that enables us to quantify the similarity between two arbitrary paths and extract the atomic-scale determinants responsible for their differences. PSA utilizes the full information available in 3N-dimensional configuration space trajectories by employing the Hausdorff or Fréchet metrics (adopted from computational geometry) to quantify the degree of similarity between piecewise-linear curves. It thus completely avoids relying on projections into low dimensional spaces, as used in traditional approaches. To elucidate the principles of PSA, we quantified the effect of path roughness induced by thermal fluctuations using a toy model system. Using, as an example, the closed-to-open transitions of the enzyme adenylate kinase (AdK) in its substrate-free form, we compared a range of protein transition path-generating algorithms. Molecular dynamics-based dynamic importance sampling (DIMS) MD and targeted MD (TMD) and the purely geometric FRODA (Framework Rigidity Optimized Dynamics Algorithm) were tested along with seven other methods publicly available on servers, including several based on the popular elastic network model (ENM). PSA with clustering revealed that paths produced by a given method are more similar to each other than to those from another method and, for instance, that the ENM-based methods produced relatively similar paths. PSA applied to ensembles of DIMS MD and FRODA trajectories of the conformational transition of diphtheria toxin, a particularly challenging example. For the AdK transition, the new concept of a Hausdorff-pair map enabled us to extract the molecular structural determinants responsible for differences in pathways, namely a set of conserved salt bridges whose charge-charge interactions are fully modelled in DIMS MD but not in FRODA. PSA has the potential to enhance our understanding of transition path sampling methods, validate them, and to provide a new approach to analyzing conformational transitions.

Rigidity loss in disordered systems: three scenarios.

We reveal significant qualitative differences in the rigidity transition of three types of disordered network materials: randomly diluted spring networks, jammed sphere packings, and stress-relieved networks that are diluted using a protocol that avoids the appearance of floppy regions. The marginal state of jammed and stress-relieved networks are globally isostatic, while marginal randomly diluted networks show both overconstrained and underconstrained regions. When a single bond is added to or removed from these isostatic systems, jammed networks become globally overconstrained or floppy, whereas the effect on stress-relieved networks is more local and limited. These differences are also reflected in the linear elastic properties and point to the highly effective and unusual role of global self-organization in jammed sphere packings.

Ring statistics of silica bilayers.

The recent synthesis and characterisation of bilayers of vitreous silica has produced valuable new information on ring sizes and distributions. In this paper, we compare the ring statistics of experimental samples with computer generated samples. The average ring size is fixed at six by topology, but the width, skewness and other moments of the distribution of ring edges are characteristics of particular samples. We examine the Aboav-Weaire law that quantifies the propensity of smaller rings to be adjacent to larger rings, and find similar results for available experimental samples which however differ somewhat from computer-generated bilayers. We introduce a new law for the areas of rings of various sizes.

Two exactly soluble models of rigidity percolation.

We summarize results for two exactly soluble classes of bond-diluted models for rigidity percolation, which can serve as a benchmark for numerical and approximate methods. For bond dilution problems involving rigidity, the number of floppy modes F plays the role of a free energy. Both models involve pathological lattices with two-dimensional vector displacements. The first model involves hierarchical lattices where renormalization group calculations can be used to give exact solutions. Algebraic scaling transformations produce a transition of the second order, with an unstable critical point and associated scaling laws at a mean coordination =4.41, which is above the 'mean field' value =4 predicted by Maxwell constraint counting. The order parameter exponent associated with the spanning rigid cluster geometry is ß=0.0775 and that associated with the divergence of the correlation length and the anomalous lattice dimension d is dν=3.533. The second model involves Bethe lattices where the rigidity transition is massively first order by a mean coordination =3.94 slightly below that predicted by Maxwell constraint counting. We show how a Maxwell equal area construction can be used to locate the first-order transition and how this result agrees with simulation results on larger random-bond lattices using the pebble game algorithm.

Fast calculation of the infrared spectra of large biomolecules.

Vibrational spectra of proteins potentially give insight into biologically significant molecular motion and the proportions of different types of secondary structure. Vibrational spectra can be calculated either from normal modes obtained by diagonalizing the mass-weighted Hessian or from the time autocorrelation function derived from molecular dynamics trajectories. The Hessian matrix is calculated from force fields because it is not practical to calculate the Hessian from quantum mechanics for large molecules. As an alternative to molecular dynamics the spectral response can be calculated from a time autocorrelation derived from numerical solution of the harmonic equations of motion, resulting in calculations at least 4 times faster. Because the calculation also scales linearly with number of atoms, N, it is faster than normal-mode calculations that scale as N (3) for proteins with more then 4,700 atoms. Using this method it is practical to perform all-atom calculations for large biological systems, for example viral capsids, with the order of 10(5) atoms.

Bond percolation in higher dimensions.

We collect results for bond percolation on various lattices from two to fourteen dimensions that, in the limit of large dimension d or number of neighbors z, smoothly approach a randomly diluted Erdos-Rényi graph. We include results on bond-diluted hypersphere packs in up to nine dimensions, which show the mean coordination, excess kurtosis, and skewness evolving smoothly with dimension towards the Erdos-Rényi limit.

Amorphous graphene: a realization of Zachariasen's glass.

Amorphous graphene is a realization of a two-dimensional Zachariasen glass as first proposed 80 years ago. Planar continuous random networks of this archetypal two-dimensional network are generated by two complementary simulation methods. In the first, a Monte Carlo bond switching algorithm is employed to systematically amorphize a crystalline graphene sheet. In the second, molecular dynamics simulations are utilized to quench from the high temperature liquid state. The two approaches lead to similar results as detailed here, through the pair distribution function and the associated diffraction pattern. Details of the structure, including ring statistics and angular distortions, are shown to be sensitive to preparation conditions, and await experimental confirmation.

Collective dynamics differentiates functional divergence in protein evolution.

Protein evolution is most commonly studied by analyzing related protein sequences and generating ancestral sequences through Bayesian and Maximum Likelihood methods, and/or by resurrecting ancestral proteins in the lab and performing ligand binding studies to determine function. Structural and dynamic evolution have largely been left out of molecular evolution studies. Here we incorporate both structure and dynamics to elucidate the molecular principles behind the divergence in the evolutionary path of the steroid receptor proteins. We determine the likely structure of three evolutionarily diverged ancestral steroid receptor proteins using the Zipping and Assembly Method with FRODA (ZAMF). Our predictions are within ~2.7 Å all-atom RMSD of the respective crystal structures of the ancestral steroid receptors. Beyond static structure prediction, a particular feature of ZAMF is that it generates protein dynamics information. We investigate the differences in conformational dynamics of diverged proteins by obtaining the most collective motion through essential dynamics. Strikingly, our analysis shows that evolutionarily diverged proteins of the same family do not share the same dynamic subspace, while those sharing the same function are simultaneously clustered together and distant from those, that have functionally diverged. Dynamic analysis also enables those mutations that most affect dynamics to be identified. It correctly predicts all mutations (functional and permissive) necessary to evolve new function and ~60% of permissive mutations necessary to recover ancestral function.

DNA Bending through Large Angles Is Aided by Ionic Screening.

We used adaptive umbrella sampling on a modified version of the roll angle to simulate the bending of DNA dodecamers. Simulations were carried out with the AMBER and CHARMM force fields for 10 sequences in which the central base pair step was varied. On long length scales, the DNA behavior was found to be consistent with the worm-like chain model. Persistence lengths calculated directly from the simulated structures and indirectly through the use of sequence-independent coarse-grained models based on simulation data were similar to literature values. On short length scales, the free energy cost of bending DNA was found to be consistent with the worm-like chain model for small and intermediate bending angles. At large angles, the bending free energy as a function of the roll angle became linear, suggesting a relative increase in flexibility at larger roll angles. Counterions congregated on the concave side of the highly bent DNA and screened the repulsion of the phosphate groups, facilitating the bending.

Protein unfolding under force: crack propagation in a network.

The mechanical unfolding of a set of 12 proteins with diverse topologies is investigated using an all-atom constraint-based model. Proteins are represented as polypeptides cross-linked by hydrogen bonds, salt bridges, and hydrophobic contacts, each modeled as a harmonic inequality constraint capable of supporting a finite load before breaking. Stereochemically acceptable unfolding pathways are generated by minimally overloading the network in an iterative fashion, analogous to crack propagation in solids. By comparing the pathways to those from molecular dynamics simulations and intermediates identified from experiment, it is demonstrated that the dominant unfolding pathways for 9 of the 12 proteins studied are well described by crack propagation in a network.

Comparison of pathways from the geometric targeting method and targeted molecular dynamics in nitrogen regulatory protein C.

Geometric targeting (GT) is a recently introduced method for rapidly generating all-atom pathways from one protein state to another, based on geometric rather than energetic considerations. To generate pathways, a bias is applied that gradually moves atoms toward a target structure, while a set of geometric constraints between atoms is enforced to keep the structure stereochemically acceptable. In this work, we compare conformational pathways generated from GT to pathways from the much more computationally intensive and commonly used targeted molecular dynamics (TMD) technique, for a complicated conformational change in the signaling protein nitrogen regulatory protein C. We show that the all-atom pathways from GT are similar to previously reported TMD pathways for this protein, by comparing motion along six progress variables that describe the various structural changes. The results suggest that for nitrogen regulatory protein C, finding an all-atom pathway is primarily a problem of geometry, and that a detailed force field in this case constitutes an unnecessary extra layer of detail. We also show that the pathway snapshots from GT have good structure quality, by measuring various structure quality metrics. Transient hydrogen bonds detected by the two methods show some similarities but also some differences. The results justify the usage of GT as a rapid, approximate alternative to TMD for generating stereochemically acceptable all-atom pathways in highly constrained protein systems.

Generating stereochemically acceptable protein pathways.

We describe a new method for rapidly generating stereochemically acceptable pathways in proteins. The method, called geometric targeting, is publicly available at the webserver http://pathways.asu.edu, and includes tools for visualization of the pathway and creating movie files for use in presentations. The user submits an initial structure and a target structure, and a pathway between the two input states is generated automatically. Besides visualization, the structural quality of the pathways makes them useful as input pathways into pathway refinement techniques and further computations. The approach in geometric targeting is to gradually change the system's RMSD relative to the target structure while enforcing a set of geometric constraints. The generated pathways are not minimum free energy pathways, but they are geometrically plausible pathways that maintain good covalent bond distances and angles, keep backbone dihedral angles in allowed Ramachandran regions, avoid eclipsed side-chain torsion angles, avoid non-bonded overlap, and maintain a set of hydrogen bonds and hydrophobic contacts. Resulting pathways for over 20 proteins featuring a wide variety of conformational changes are reported here, including the very large GroEL complex.

Flexibility of ideal zeolite frameworks.

We explore the flexibility windows of the 194 presently-known zeolite frameworks. The flexibility window represents a range of densities within which an ideal zeolite framework is stress-free. Here, we consider the ideal zeolite to be an assembly of rigid corner-sharing perfect tetrahedra. The corner linkages between tetrahedra are hard-sphere oxygen atoms, which are presumed to act as freely-rotating, force-free, spherical joints. All other inter-tetrahedral forces, such as coulomb interactions, are ignored. Thus, the flexibility window represents the null-space of the kinematic matrix that governs the allowable internal motions of the ideal zeolite framework. We show that almost all of the known aluminosilicate or aluminophosphate zeolites exhibit a flexibility window. Consequently, the presence of flexibility in a hypothetical framework topology promises to be a valuable indicator of synthetic feasibility. We describe computational methods for exploring the flexibility window, and discuss some of the exceptions to this flexibility rule.

The long-wavelength limit of the structure factor of amorphous silicon and vitreous silica.

Liquids are in thermal equilibrium and have a non-zero structure factor S(Q --> 0) = [-(2)]/ = rho(0)k(B)Tchi(T) in the long-wavelength limit where rho(0) is the number density, T is the temperature, Q is the scattering vector and chi(T) is the isothermal compressibility. The first part of this result involving the number N (or density) fluctuations is a purely geometrical result and does not involve any assumptions about thermal equilibrium or ergodicity, so is obeyed by all materials. From a large computer model of amorphous silicon, local number fluctuations extrapolate to give S(0) = 0.035 +/- 0.001. The same computation on a large model of vitreous silica using only the silicon atoms and rescaling the distances gives S(0) = 0.039 +/- 0.001, which suggests that this numerical result is robust and perhaps similar for all amorphous tetrahedral networks. For vitreous silica, it is found that S(0) = 0.116 +/- 0.003, close to the experimental value of S(0) = 0.0900 +/- 0.0048 obtained recently by small-angle neutron scattering. Further experimental and modeling studies are needed to determine the relationship between the fictive temperature and structure.

Hierarchical plasticity from pair distance fluctuations.

Observations, experiments and simulations often generate large numbers of snapshots of configurations of complex many-body systems. It is important to find methods of extracting useful information from these ensembles of snapshots in order to document the motion as the system evolves in time. Some of the most interesting information is contained in the relative motion of individual constituents, rather than their absolute motion. We present a novel statistical method for identifying hierarchies of plastically connected objects in a system from a series of two or more snapshot configurations. These plastic clusters are distinctive in that although their members tend to remain loosely connected, the clusters may be deformed plastically. This method is demonstrated for a number of systems, including an exactly soluble freely jointed polymer chain model, a two-dimensional simulation of two species of interacting bodies and a protein. These concepts are implemented as TIMME, the Tool for Identifying Mobility in Macromolecular Ensembles.