article: Large macromolecular assemblies are often important for biological processes in cells. Allosteric communications between different parts of these molecular machines play critical roles in cellular signaling. Although studies of the topology and fluctuation dynamics of coarse-grained residue networks can yield important insights, they do not provide characterization of the time-dependent dynamic behavior of these macromolecular assemblies. Here we develop a novel approach called Perturbation-based Markovian Transmission (PMT) model to study globally the dynamic responses of the macromolecular assemblies. By monitoring simultaneous responses of all residues (>8,000) across many (>6) decades of time spanning from the initial perturbation until reaching equilibrium using a Krylov subspace projection method, we show that this approach can yield rich information. With criteria based on quantitative measurements of relaxation half-time, flow amplitude change, and oscillation dynamics, this approach can identify pivot residues that are important for macromolecular movement, messenger residues that are key to signal mediating, and anchor residues important for binding interactions. Based on a detailed analysis of the GroEL-GroES chaperone system, we found that our predictions have an accuracy of 71–84% judged by independent experimental studies reported in the literature. This approach is general and can be applied to other large macromolecular machineries such as the virus capsid and ribosomal complex.
  In the cellular environment, biological processes often involve large assemblies of macromolecules, such as ribosome complex in protein synthesis, assembly of virus capsid, and chaperone systems assisting protein folding. Understanding how signaling communications are transmitted and the dynamics of such allosteric communications in these biological nanomachines at multiple time scale is an important problem. To explore the roles of large-scale dynamics fluctuations of proteins, coarse grained models such as Gaussian network models (GNM) and elastic network models (ENM) have been developed and are now widely used [1]–[13]. These models are based on a simplified representation of protein structures. The eigen modes of these models are then used to explore the fluctuation dynamics of proteins around their native conformations. Chennubhotla and Bahar have further extended the network model by introducing a Markovian process to describe how signal propagates in proteins [8]. In another study, by clustering residues into groups based on simplifications of the Markovian transition matrix, models of allosteric communication paths can be constructed [8]. By examining pair correlation between residues at slow eigen modes, clusters of coupled residues in potassium channel that are important for inter-subunit cooperativity are identified [14]. Recent work based on topological analysis of shortest paths of network models showed that residues mediating signaling in proteins and modules of protein architecture can be identified [15], [16]. A drawback of current coarse grain network model based studies is that they do not characterize the time-dependent behavior of dynamic changes in protein assemblies. Furthermore, they often require additional mode analysis to assess the relevance of fluctuations associated with individual eigen mode [8]. Although the dominant modes most important for biological functions are often understood to be the slow modes [7], [8], their precise identification often cannot be determined a priori and requires significant amount of analysis with additional experimental information [10]. Frequently, the most important mode is different for different proteins. A recent study by Zheng, Brooks, and Thirumalai showed that by locating robustly conserved residues, evolutionary analysis of proteins through multiple sequence alignment can be helpful in identifying important dynamic modes [10]. In this study, we postulate that rich information about allosteric dynamics relevant for biological function can be directly obtained from native structures without any mode analysis or multiple sequence alignment. Our approach builds upon the Markovian stochastic model introduced by Chennubhotla and Bahar [8], and our earlier preliminary work [17]. Chennubhotla and Bahar have firmly established the connection between the Markovian transmission model and the physical elastic network model through the mapping between Markovian hitting time and the equilibrium fluctuation dynamics of residues [18]. Our aim here is different, and it is to study how different parts of a macromolecular machinery respond to signal perturbation. Here the perturbation is given as the initial condition, and our conjecture is that the dynamic response of residues can reveal biologically important information. Our approach is termed as the Perturbation-based Markovian Transmission (PMT) model. The main task of our study is to characterize the temporal dynamic response to signal perturbation across all residues in a large macromolecular assembly at all time scales, i. e. , from the beginning perturbation until the equilibrium state is reached. This is a challenging task, as it is difficult to follow the precise time trajectories of several thousand or more residues simultaneously over or more decades of time span. For this purpose, we show that using a matrix subspace projection method, the underlying large scale master equation governing the dynamics of Markovian signal propagation can be solved accurately. We also show that simple perturbation can reveal intrinsic dynamic properties of proteins. In addition, we present new results on the characterization of functionally important residues and the mechanism of signal transduction in the large macromolecular assembly of the GroEL-GroES system, which has been studied previously using Gaussian network model [8] and Brownian dynamics simulations [10].
  As an example, we use the PMT model to study the dynamic response of the chaperone complex in Escherichia coli. The function of this ATP regulated chaperone complex is to assist the folding of unfolded and misfolded substrate proteins while changing conformations among the T, R, and R″ states (pdb, and, respectively [23]–[25]). The structures of the GroEL-GroES chaperone are shown in Fig 1. We briefly review our current understanding of the allosteric communications and functional activity of the GroEL-GroES chaperone complex in E. coli. The cycle of conformational change starts from the T state, which contains neither the ADP nor the ATP molecules. In addition, the GroES subunit is not yet bound to the GroEL subunit (Fig 2). In step 1 during the transition from the T state to the R state, one of the trans-rings takes seven ATP molecules and interact with the unfolded or misfolded protein substrate. This is followed by conformational changes of the trans-ring to become a cis-ring, which is observed in the R state structure. The ATP molecules in the R state are then hydrolyzed to become ADP molecules. This hydrolysis triggers further conformational change, and the complex transits from the R state to the R″ state in step 2. In the meantime, the GroES subunit binds to the top two helices (H and I) of the cis-ring chains of the GroEL subunit in the R″ state, and the protein substrate is released into the central cavity of the chaperone complex. The last step of the process occurs after a new unfolded protein ligand bind to the original trans-ring chains of the GroEL subunit in the R″ state [26]. The effect is that the cis-ring in the R″ state transits and take the conformation of the trans-ring in the T state, releasing the the GroES subunit, the protein substrate, and the seven ADP molecules during the process. Now the complex goes back to step 1 and ready for another cycle [26], [27]. To understand the dynamic response of the chaperone GroEL-GroES complex, we analyze details of the time trajectories of the residues upon perturbation. We are interested in the response dynamics of residues in the chaperone complex. We describe a few examples of residues experimentally known to be functionally important. Pivot, messenger, and anchor residues have distinct patterns in dynamic responses to perturbation. These patterns are different and are mutually exclusive. A residue can be either a pivot, a messenger, or an anchor residue in a specific conformational state, but the same residue cannot have two or more different patterns of dynamic responses simultaneously in the same state. We postulate that additional functionally important residues can be identified through analysis of their dynamic responses upon perturbation. Below we describe the prediction criteria and additional residues predicted to play functionally important roles. In the majority of the cases, our predictions are supported by experimental studies. A complete list of identified pivot, messenger, and anchor residues can be found in Table 1, Table 2, and Table 3, respectively, along with corresponding experimental evidence supporting their roles in the function of GroEL-GroES. Details can be found in Supporting Information (Text S1). In this study, we have developed a new approach, called the Perturbation-based Markovian Transmission (PMT) model, to study the dynamic properties of large macromolecular assemblies. In this model, we probe the macromolecular assembly with an initial perturbation, and assume perturbations are transmitted by a Markovian processes through the connected network defined by the contacts between atoms and residues. This method enables us to monitor the time-dependent responses of all residues simultaneously from the moment of perturbation until all residues reach the equilibrium state. Our method can also be modified to study signal transmission pathways by applying initial perturbation to restricted local region such as binding surfaces. It is effective for very large macromolecular assemblies, requiring relatively small amount of computational time. For example, the computation of probability flows for all 8,000 residues in the GroEL-GroES system takes about 5 minutes on a machine with 2 dual core AMD Opteron processor (1. 8 GHz) and 4G memory. We have applied the PMT model to study the GroEL-GroES protein chaperone system. With an initial in silico uniform perturbation to all residues in the T, R, and R″ conformational states, we follow the time-dependent response of 8,000 residues simultaneous until the stationary state is reached after 5 decades in time. Our analysis of the responses of experimentally identified residues known to be functionally important show characteristic patterns. We have further characterized and categorized specific patterns for pivot residues, messenger residues, and anchor residues important for signal propagation in the GroEL-GroES system. Based on these patterns, we have further predicted additional residues that are likely to play important functional roles. Our predictions are largely consistent with experimental results on their roles during the allosteric conformational changes between states. In addition, we have also made a number of novel predictions of functionally important residues, which can be verified by experimental tests. A model of signal transmission pathway has also been developed by applying the intial perturbation to residues in local binding surface. The successful application of the PMT model to the GroEL-GroES chaperone system for studying its dynamic properties shows that this technique works well. It can be applied to other large macromolecular machineries for detailed computational studies. We envision that the PMT method can be applied to virus capsid, ribosomal complex, and other large allosteric multi-protein assemblies.
  The structure of the GroEL-GroES chaperone complex in the T, R, and R″ states are obtained from the Protein Data Bank (). The pdb IDs are, and, respectively [23]–[25]). The GroEL-GroES complex of the R″ state consists of 14 chains for the homo-oligemeric GroEL subunits, and 7 chains for the GroES subunits, with a total of 21 chains and 8,015 residues. The GroEL subunit structure contains 7 of the 14 chains (3,829 residues) in the T state, and 14 chains (7,658 residues) in the R state (Fig 1). Since there are multiple symmetrically related identical chains in each ring, we take the average values of the time-dependent probability flow of equivalent residues from different chains located in the same ring. We obtain the contact network of a protein molecule by identify contacts between different residues. Two atoms from different residues are regarded to be in contact if the Euclidean distance between them is Å. A residue is represented as a node, and an edge between residues exists if they have atomic contacts. The perturbations are signals and can be transmitted from a residue to its neighboring residues in contact with a probability flow following a Markovian processes. When the perturbation is transmitted from residue to a neighboring residue, the probability flow in our PMT model is defined as: if there is no atom-atom contacts between residue and, and if there is at least one atom-atom contacts between residue and. Here is the number of atom-atom contacts between residue and. If, we set equal to the number of atoms without contact with any atoms from other residues. The Markovian transition matrix is a stochastic matrix because each column sums up to 1: . In our Markovian model, all of the atom-atom contacts of residue with other residue (s) will transmit the same amount of probability flow. This is different from the model used in reference [8], where the probability flow is not equally distributed among different atom-atom contacts, as the probability flow in that model is weighted by both the size of the source residue and the size of the target residue: , where is the number of atoms in residue. For a Markovian process, the final stationary state distribution of the signal is independent of the initial perturbation, and has the property: . This unique stationary distribution can be written out in closed-form [40]. For our network model, we have: where. The probability at the stationary state for residue is the same as the fraction that the number of atomic contacts makes with other residues, plus the number of its non-interacting atoms, with the sum of all atomic contacts plus all non-interacting lone atoms in the whole structure (Supporting Information Text S1). To analytically obtain the full dynamics response of residues in the PMT model, we discuss the Krylov subspace matrix reduction method for solving the master equation. Noting the sparsity of the Markovian matrix, we follow the approach of Sidje [22]. Based on the analytical solution of matrix exponential shown in Eq. (3), one can expand in the Krylov subspace, which is defined as: (4) Denoting as the 2-norm of a vector or matrix, our approximation then becomes, where is the first unit basis vector, is a matrix formed by the orthonormal basis of the Krylov subspace, and the upper Heisenberg matrix, both computed from an Arnoldi algorithm [41], [42]. The error can be bounded by. We now only need to compute explicitly the matrix exponential. Because is much smaller than the total number of states, this is a simpler problem. A special form of the Padé rational of polynomials instead of Taylor expansion is used [22]: , where and. All these calculations are done using Sidje' s expokit software package [22]. We select [22]. We use sensitivity, specificity, and accuracy to assess the performance of our predictions in identifying functionally important residues. Sensitivity is defined as, specificity as, and accuracy as, where, , , , and are the numbers of true positive, false positive, true negative, false negative, and total number of residues, respectively.
summary: Biological processes in a cell often require complex molecular machineries with large macromolecular assemblies as components. An example is the chaperone system in the bacterium E. coli, which helps proteins to fold correctly. In these macromolecular machineries, signals are transmitted dynamically in order for biological functions to be carried out. Studying the dynamic process of signal transmission helps us to identify key elements of the macromolecular assemblies that are pivots for dynamic motions, communicators for interfacing with other molecules, and anchors that are key for signal transmission. In this study, we describe a novel computational method that can globally survey the dynamic responses of the macromolecular machinery to perturbation over the full time course by monitoring simultaneously all the elements at the amino acid residue level and at multiple time spans, from the initial perturbation until the system reaches equilibrium. We show that the key residues predicted by our computational method in the chaperone system of E. coli to a large extent are correct, as they often coincide with the ones identified by experimental studies. We also show that this computational method can make novel predictions about the importance of additional amino acid residues previously uncharacterized, which can be further tested in experimental studies. This approach can be applied to study other large macromolecular assemblies such as the virus capsid and ribosomal complex.
keywords: computational biology
title: Perturbation-based Markovian Transmission Model for Probing Allosteric Dynamics of Large Macromolecular Assembling: A Study of GroEL-GroES