Coarse-Grained Description of Protein Internal Dynamics: An Optimal Strategy for Decomposing Proteins in Rigid Subunits

Covariance matrix and essential dynamical spaces

A natural formulation and motivation of the method discussed here is provided in terms of the few collective coordinates that best account for the structural fluctuations in a given protein (2). These generalized coordinates, also termed essential dynamical spaces or low-energy modes, are aptly identified as the eigenvectors associated to the largest eigenvalues of the covariance matrix, C,

where the brackets denote the canonical ensemble average over structural configurations and δr_i^μ indicates the μ^th Cartesian component of the vector displacement of the ith C_α from the average reference position. The covariance matrix is obtained either from atomistic molecular dynamics simulations or from exactly-solvable elastic network approaches (3,14,24,25). In the latter case, the canonical weight of a protein configuration is controlled by a model potential energy, F, that is quadratic in terms of the displacements of the amino acids from a preassigned reference structure, $F={\sum }_{\text{ij},\mu \nu }\delta r_{\text{i}}^{\mu }{M}_{\text{ij}}^{\mu \nu }\delta r_{\text{j}}^{\nu }$. The quadratic nature of F implies that the eigenvectors of C associated to the largest eigenvalues correspond to the lowest-energy modes of fluctuation of the system.

The elastic network model employed here is based on the β-Gaussian model of Micheletti et al. (25), where each amino acid is represented by two centroids (for the backbone and side chain). To avoid artifacts in the dynamical domain decomposition arising from overestimation of the mobility of exposed loops/termini, the range of the pairwise interaction of centroids at distance x is prolonged beyond the default value of 7.5 Å and weighted by an exponentially decreasing term, exp[−(x − 7.5 Å)/2 Å].

Rigid-block decomposition

The collective character of low-energy modes in proteins suggests that the internal dynamics of these biomolecules might be described in terms of the relative motion of a limited number of approximately rigid subunits. This picture holds if the distance fluctuation of amino-acid pairs within the blocks is appreciably smaller than for pairs in different blocks. In the following we shall discuss a general variational framework, apt for numerical implementation, for performing an optimal decomposition of a protein into a preassigned number of approximately-rigid groups.

A convenient criterion of optimality is given by the maximization, over the possible groupings, of the system mean-square fluctuation captured by the lowest-energy modes that is compatible with a nearly-rigid character of each group of amino acids.

For definiteness we discuss a subdivision of a protein's amino acids into Q putatively-rigid groups. We shall further indicate with ${\overrightarrow{v}}_{\text{l}}$, the l^th (nonzero) lowest-energy mode of the system $C{\overrightarrow{v}}_{\text{l}}={\lambda }_{\text{l}}{\overrightarrow{v}}_{\text{l}}$, with λ_l being the associated eigenvalue. We consider the following decomposition of the mode

where ${\overrightarrow{v}}_{\text{l}}^{\text{rb}}$ is the best rigid-body fit of the mode and $\Delta {\overrightarrow{v}}_{\text{l}}$ is the correction that accounts for the internal distortions within the putatively rigid blocks. Notice that the spatial contiguity of amino acids belonging to the same group is not enforced a priori. The first term of the right-hand side is entirely specified by 6Q parameters that correspond to the roto-translations of each rigid-body block. The modes are assumed to describe perturbative fluctuations of the reference structure (it should be borne in mind that this condition may be poorly met by highly mobile regions such as exposed loops or termini). The rigid-body motion of a group of residues consequently admits a linear parameterization in terms of the roto-translational degrees of freedom. In fact, the l^th modal displacement of the ith amino acid belonging, say, to the q^th group, is given by

where ${\overrightarrow{r}}_{\text{q}}^{\text{cm}}$ denotes the center of mass of the q^th group and ${\overrightarrow{r}}_{\text{i}}$ is the position of the ith C_α. Notice that the same translation and rotation parameters (${\overrightarrow{t}}_{\text{l}}\text{and}{\overrightarrow{\omega }}_{\text{l}}$, respectively) are used for each residue in a given group. The optimal rigid-body approximation to ${\overrightarrow{v}}_{\text{l}}$ is obtained by maximizing the norm of ${\overrightarrow{v}}_{\text{l}}^{\text{rb}}$ over the 6Q-dimensional parameter space. The accuracy of the fit is readily obtained by computing the fraction of the norm of the essential dynamical space captured by the rigid-body approximation:

Considering the space of the top n essential modes, the above expression generalizes to

Unless otherwise stated, considerations will be limited, as customary, to the n = 10 lowest-energy modes, which are usually sufficient to capture most of the fluctuations in a given system (alternatively n could be chosen so to capture a preassigned fraction of the overall internal fluctuations). The goal of the procedure is to identify, within the possible partitioning of amino acids into Q dynamical domains, the one yielding the largest possible value of f.

Elastic deformation energy

Because the volume of configuration space (number of distinct amino acids groupings) is large, the optimization of f in Eq. 5 requires the stochastic exploration of tens of thousands of possible amino-acid groupings. Even if optimized algorithms exist for the calculation of the rigid fits (26), the number of repeated matrix operations involved is so large that it is not computationally convenient to extremize directly the quantity f. A more effective strategy is to perform a preliminary exploration of configuration space by optimizing a simple objective function to efficiently identify candidate subdivisions over which f is finally evaluated and maximized.

The objective function that we consider is

where σ_i = 1 ··· Q denotes the group to which amino acid i belongs, δ is the Krönecker delta, ${\overrightarrow{d}}_{\text{ij}}^0$ is the distance vector of amino acids i and j in the reference conformation, R_c is an interaction cutoff distance set equal to 7 Å, and n = 10. The sought optimal grouping of amino acids is the one that minimizes F, similarly to the spirit of graph clustering methods (27). For systems consisting of truly-rigid subparts, this will be analogous to maximizing f. The first term in the sum represents the cost of the average elastic energy associated to the internal deformation of the molecule. This term penalizes fluctuations in the distance of any two points belonging to the same putatively-rigid group, consistently with the definition of rigid bodies. The second term introduces a penalty, controlled by the parameter α ≥ 0, for dynamical domains consisting of regions that are disconnected in space. Upon increasing α, in fact, the term disfavors the number of pairs of neighboring amino acids (those closer than the cutoff distance R_c = 7 Å) that belong to different groups. The optimization of F, therefore, leads to group assignments that minimize the interface area between the groups, while not strictly enforcing the spatial compactness of the domains. The minimization of F is straightforwardly done within a simulated annealing protocol, with elementary moves corresponding to changes of the group assignment of individual amino acids. The corresponding changes of F only require the summation of N precalculated quantities, N being the number of amino acids in the protein. Proteins of 200 residues can thus be subdivided into, e.g., 10 dynamical domains in ∼1 min, on present-day personal computers.

The search for the optimal solution is carried out separately for increasing values of α. Notice that for sufficiently large α, the minimization of F eventually leads to solutions having fewer groups than Q, which are hence not considered (this is intuitively expected, as in the limit α → ∞ the presence of boundaries is forbidden, and a single dynamical domain is returned by the minimization of F). The solution corresponding to the largest value of f is taken as providing the best subdivision.

Adenylate kinase

Adenylate kinase is a phosphotransferase regulating the relative abundance of AMP, ADP, and ATP within the cell. The enzyme is composed by a central core and two domains, the ATP binding domain (Lid) and the AMP binding one, which are highly mobile. In the available “closed” crystallographic state (Protein Data Bank (PDB) code: 1ake), they are displaced toward the core by >7 Å, with respect to the “open” crystal structure (PDB code: 4ake).

Recent experiments have indicated that the enzyme is spontaneously capable of interconverting between the closed and open forms even in the absence of ligands (8,9). This points to the predisposition of AKE's internal dynamics to bridging the open/closed conformations and to the absence of large free-energy barriers separating the two reference states, consistently with indications from atomistic molecular dynamics simulations of AKE (10,29,30).

For example, in a recent computational/theoretical study carried out by some of us (30), the dynamical evolution of the free E. coli adenylate kinase was followed from two starting structures for as long as 100 ns. Over this extensive time-span, which is nevertheless much smaller than the experimental interconversion time, the molecule was found to populate a fair number of structurally distinct substates. Most of the structural fluctuations within and across the substates were described by very few low-energy collective modes entailing the independent motion of the Lid and AMP-bd subdomains with respect to the core (see Figs. 6 and 9 in (30)).

The rigid-block decomposition scheme here applied to AKE lends naturally to assessing if, and to what extent, the molecule's internal dynamics can be described in terms of a few parts that move as nearly-rigid units. We begin by considering the fluctuations within the substate where the 50-ns-long trajectory started from the open structure, 4ake, dwelled for ∼10 ns (30). The reference structure for the substate, which is the most populated of the MD trajectory, is provided in Fig. 1 a, along with the representation of the lowest energy mode. The mobility the Lid and of the AMP-binding subdomains, corresponding to regions 117–164 and 30–64, respectively, is evident.

The n = 10 lowest energy modes within the substate were used to subdivide the enzyme into Q = 2, 3…10 dynamical domains. A representation of the subdivisions into three and four groups is provided in Fig. 1, b and c. The fraction of essential dynamics motion (see Eq. 5) captured by the various subdivisions is shown in Fig. 1 d.

The graph indicates that a very limited number of dynamical domains is already sufficient to account for most of the essential dynamics. In fact, subdivisions into Q = 2, 3, and 4 blocks capture as much as 52%, 77%, and 83% of the fluctuations entailed by the n = 10 essential modes (which account for 80% of the overall mobility).

The subdivision for Q = 2 identifies region 122–156 as an approximately rigid, but highly mobile, unit. The region overlaps well with the Lid indicated before. The less mobile AMP-binding domain is identified as a distinct unit when using Q = 3. In fact, for Q = 3, the regions corresponding to the two mobile nearly-rigid subdomains are 122–158 and 32–59, and are compatible with the customary tripartite subdomain division of AKE. If the entire 50-ns-long trajectory is used rather than the most populated substate, it is found that the boundary of the AMP-binding domain is virtually unaltered (sequence interval 32–60). The larger configurational space spanned by the more mobile Lid domain instead reflects into an extension of the both the left and right subdomain boundaries by ∼10 residues, thus covering the interval 112–167.

In summary, for Q = 3, the three units cover five sequence intervals: one for each of the two mobile domains and three for the nearly-fixed core. It is interesting to compare this dynamics-based subdivision with the one provided by the TLS analysis of crystallographic data (22), which enforces the sequence continuity of each rigid block. The TLS decomposition of 4ake into five intervals (as many as those found with Q = 3) returns the following segments: 1–27, 78–116, and 171–214, identifiable with the core, and 28–77 and 117–170, compatible with the AMP-bd and Lid subdomains, respectively. With the exception of one of the AMP-bd/core boundaries, the TLS subdivisions and those of our analysis of the full MD trajectory are mismatched by only approximately five residues. They are hence generally consistent, despite the differences not only in method but also for the nature of the input data (crystallographic B-factors for TLS and MD data for our method). An important distinction between the two results is, however, that the three segments constituting the core regions are encompassed in a single rigid unit by this variational method, while the others are treated as independent ones within the TLS scheme.

The optimal subdivision was compared also with the one returned by the DynDom server (16), which requires the input of two structures representing the conformational variability of the molecule of interest. Accordingly, from the set of MD-sampled conformers we selected the pair with the largest root mean-square deviation. DynDom returned a subdivision in two domains, the smallest corresponding to the Lid (sequence interval 110–169) plus a small loop (residues 6–12) and the other to the core plus the AMP-binding domain. Interestingly, this latter subdomain is recognized as a separate dynamical domain if the open and closed crystallographic conformers of AKE (1ake, 4ake) are used as input structures.

We conclude by commenting on the subdivision in Q = 4 dynamical units of the most populated MD substate. With respect to the Q = 3 case, the boundaries of the two mobile domains are only slightly adjusted to 35–60 and 118–159, respectively. However, a new domain, comprising several sequence segments 7–25, 108–117, 160–174, and 195–214, is identified at the interface between the core and the Lid. This group of hinge residues have consistently been shown, by independent methods (28,30), to be subject to a significative strain during the free enzyme dynamical evolution.

HIV-1 protease

As a further example of the dynamics-based decomposition we consider the HIV-1 protease dimer complexed with a peptide substrate. The internal dynamics of the enzyme has a distinctive collective character that has been extensively investigated in past years by means of both atomistic MD simulations as well as coarse-grained models (see (5,7,33–35) and references therein). The enzyme flaps, in fact, carry out particular large-scale movements, resulting in detectable anticorrelated displacements of the flap's tip, which contacts the bound (inhibiting) peptide and the distal region of the flaps, where several mutations causing drug resistance are located (5).

To illustrate the applicability of the method in the absence of data from atomistic simulations, we obtained the essential dynamical spaces from the βGM elastic network model of Micheletti et al. (25), modified as described in the Methods section.

The complex shown in Fig. 1 e, which corresponds to the equilibrium structure of the MD study of Piana et al. (5), was again subdivided from 2 to 10 domains (see Fig. 1, f–h). The fraction of internal dynamics captured by the various decompositions is shown in Fig. 1 i. The curve has a slightly slower increasing trend compared to AKE (Fig. 1 d). In fact, when Q = 3, 4 domains are used, ∼72% and 79%, respectively, of the essential dynamical fluctuations is captured for the HIV-1 protease/substrate complex.

The subdivisions into Q = 2, 3, 4 approximately-rigid units are represented in Fig. 1, f–h. As for AKE, the units are compactly organized in space but do not cover a single stretch of the primary sequence. The sequence-disconnected nature of the domains does not lead simply to a detailed comparison with the TLS decomposition. We shall therefore restrict ourselves to considering the primary hinge-points, represented by amino acids 20, 35, 57, and 70, which emerge from the precalculated subdivision offered by the TLS webserver of the HIV-1 PR monomer (PDB structure 1t3r) in five-to-seven intervals. The first three hinges fall within three amino acids (along the primary sequence) from boundaries indentified for the optimal subdivision in two primary domains, Q = 2 (see the Supporting Material), suggestive of good consistency.

As a further comparison we considered the precalculated subdivision of the HIV-1 PR monomer offered by the DynDom server (based on structures 1aid and 1hsg). The returned subdivision consisted of two domains, the smaller one comprising segments 32–60 and 75–77, and broadly corresponding to the monomer flap. Though it should be borne in mind that the subdivision might depend on the fact that only one monomer is considered (multimers are not accepted by the DynDom server), the identified modular nature of the flaps is compatible with salient aspects of the TLS and variational decomposition.

The inspection of the optimal subdivisions in Fig. 1, f–h, prompts two considerations. The motion of the flaps is largely consistent with a coordinated rotatory movement around the central fulcrum regions. It is evident from the Q = 3, 4 cases that within the nearly-rigid parts comprised by the flaps, the points at the two extremes are displaced in opposite directions. On the one hand, this feature illustrates that the motion of nearly-rigid units in proteins can be sufficiently general to allow for the presence of anticorrelated motion within its constituents parts. On the other hand, the analysis supports the qualitative description of the flap motion first given by Piana et al. (5), based on the visual inspection of the first essential mode of a multi-nanosecond MD simulation (see Fig. 6 b in (5)). Indeed, if the first mode only is used for decomposing the proteins in to Q = 2 blocks, it is found that each entire flap is identified as a nearly-rigid unit (see the Supporting Material).

The second observation regards the location of the catalytic site of HIV-1 protease with respect to the “primary dynamical boundaries.” By the latter, we mean the boundaries separating the most prominent rigidlike regions in a protein (i.e., when using Q = 2 or Q = 3). By inspecting Fig. 1, f–h, it is seen that the highlighted catalytic aspartic dyad, which has low mobility, straddles the rigid-domains interface for Q = 2 and is close to one or more domain boundaries for Q = 3 and 4.

This suggests that the proximity of the catalytic amino acids to the primary dynamical boundary is instrumental for accompanying the limited mobility of the aspartic dyad with a functionally-oriented modulation of the bound peptide substrate (38).

Catalytic site location

It is interesting to note that modulations of the active site analogous to HIV-1 protease have been found for several other proteolytic enzymes differing by catalytic chemistry and structural architecture (31,32,39). The existence of such common large-scale movements is taken as the starting point for examining what relationship, if any, exists between the location of the cleavage sites and the proximity of primary dynamical boundaries for other enzymes belonging to hydrolases. Is the cleavage site commonly located near primary dynamical boundary, as for HIV-1 protease?

The question appears particularly appealing also in view of recent considerations made by del Sol et al. (11) that functional sites in proteins with allosteric behavior are preferentially located at the boundary between regions that are modular in terms of contacting amino acids.

The question will be formulated within a rather comprehensive framework, where proteolytic enzymes are considered along with other enzymes members of class 3 of EC. The enzymes were taken from the list of 76 representatives of the main EC and CATH groups singled out in Zen et al. (39). The list, restricted for simplicity to monomeric enzymes (following the indication in annotated UNIPROT (40) entries), is reported in Table 1 along with the EC and CATH code and with the indication of the amino acids constituting the catalytic site.

For each entry, the boundary between the two primary dynamical domains was identified from the Q = 2 subdivision. To measure the separation of a catalytic residue from the primary dynamical boundary, we considered the distance of its C_α from the nearest C_α belonging to the other dynamical domain. The normalized distribution of these boundary distances is shown with a thick line in Fig. 2 along with the reference distribution (dashed line) of the boundary distance of every amino acid in the 15 proteins. The two distributions present appreciable differences as the catalytic residues are preferentially closer to the interface than other amino acids in the proteins.

The overall indication of catalytic-site/boundary proximity in hydrolases conveyed by Fig. 2 was complemented by a case-by-case analysis of the 15 enzymes in Table 1. This detailed investigation is necessary in view of the fact that the cumulated data in Fig. 2 reflect properties of a group of enzymes with a certain heterogeneity in length, structural architecture, and number of catalytic sites.

The Q = 2 subdivisions of the 15 enzymes are provided as Supporting Material and consistently reveal the good proximity of the cleavage sites with the boundaries between the dynamical domains. Here we limit the discussion to three enzymes whose dynamical role in the functional cleavage of peptides or nucleic acids has been previously considered (41–45), namely: exonuclease III (PDB 1ako); human adenovirus proteinase (PDB 1avp); and endo-1,3-1,4-β-D-glucan 4-glucanohydrolase (PDB 2ayh). Their Q = 2 subdivision is represented in Fig. 3, a–c.

Exonuclease III and adenovirus proteinase bind DNA in double- and single-stranded forms, respectively. In Zen et al. (39), a dynamics-based connection between them was established, which is particularly interesting as they are not evolutionarily related and are characterized by two different architectures, 4-Layer Sandwich (CATH: 3.60.10.10) and 3-Layer (αβα) Sandwich (CATH: 3.40.395.10), respectively. In both cases the catalytic residues are found to be located at the primary boundary. As visible in Fig. 3, the low-energy modes have a common character as they entail an outward/inward concerted movement between the two blocks in the surroundings of the catalytic sites, with the latter at the center. The analysis carried out on endo-1,3-1,4-β-D-glucan 4-glucanohydrolase also shows that the two catalytic residues of the enzyme are located in proximity to the interface between the two primary dynamical domains, both surrounded by loops forming a groove that can arguably accommodate the corresponding ligand.

The consistent indication of Fig. 3 is that the catalytic site is located close to the primary interface. This fact appears particularly interesting when considering how the primary boundary is modulated by the lowest-energy modes of the domains (which are compatible with the opening/closing of the catalytic cleft (41,44,45)). By comparison to noninterfacing amino acids, it is found that interface residues cover a fairly large range of values both for overall mobility and for the degree of distortion of the local structural environment (see the Supporting Material). Interestingly, the catalytic site is accommodated at, or close to, an interface subregion having both low mobility and low-structural deformation. While these properties are consistent with the expected rigidity of the catalytic region, it is interesting that they can be realized in proximity of the primary dynamical boundary, where appreciable elastic strain can be built up due to the relative motion of the dynamical domains.

Approximately-rigid units: connectedness in sequence and space

We conclude the analysis with a systematic evaluation of the extent to which the subdivision of a protein into a limited number of approximately-rigid units results in dynamical domains that are compact in space and/or cover uninterrupted regions of the primary sequence. The interest in this question is twofold. On the one hand, it can provide indications on the viability, for computational efficiency, of enforcing a priori the proximity in sequence or space of the amino acids belonging to the same group. On the other hand, it can shed some light on the existence of consistent modular organizations of proteins at the level of sequence, structure, and dynamics.

An interesting general context where these questions can be posed is provided by multidomain proteins. We considered a data set of 90 protein monomers, with overall sequence identity <90% and constituted by three-to-six CATH domains each consisting of a single sequence interval (see the Supporting Material). Each protein was subdivided into a number of dynamical domains equal to the number of CATH domains. To analyze robust aspects of the sequence-integrity of the dynamical domains, the rigid-units subdivisions were postprocessed to eliminate domain fragments covering excessively short sequence intervals. Specifically, fragments smaller than 1/20th of the protein length (and in any case no longer than 10 amino acids) were removed. The amino acids in these fragments are reassigned to the nearest flanking unit. The resulting dynamical domains subdivisions (along the primary sequence) were compared with the ones provided by CATH.

It was found that only for 30 proteins out of 90, did the number of domain boundaries along the sequence match. Therefore, in two-thirds of the cases, the dynamical domain subdivisions gathered regions that were disconnected along the primary sequence, at variance with structural subdivisions employed in domain identifications. Notably, for the corresponding 30 cases, the dynamical subdivisions were very well consistent with the CATH ones. In fact, out of the 91 boundaries occurring in the 30 proteins, as many as 71 occurred at a separation of <10 residues of the CATH ones. By commonly-employed criteria (46) this reflects a very strong agreement of the subdivisions. It is worth noting that also for the 60 nonmatching proteins, most of the CATH subdivisions fall within 10 residues from the dynamical ones, which are, however, more numerous.

For all the 90 proteins we checked the extent to which the non-postprocessed dynamical domains, despite possibly comprising segments that are not contiguous in sequence, occupy compact regions in space. The compactness of a domain was ascertained by measuring the diameter of the graph given by the contact map of the residues (with a contact cutoff distance of 7.5 Å). A finite value of the diameter, which measures the minimum number of graph edges that need to be traversed for connecting any two nodes in the graph, indicates the spatial compactness of the domain. It was found that <5 dynamical domains out of 308 comprised disconnected, though nearby, regions. This provides an a posteriori indication of the fact that rigidlike units comprise amino acids that occupy spatially connected regions.

Finally, despite the differences in sequence integrity of the dynamical and CATH subdivisions, we performed a test to quantify the overlap between these two decompositions. We found that the mutual one-to-one overlap of the CATH domains and rigid units was, on average, 80%, which underscores a nontrivial, albeit not perfect, consistency between the two subdivision criteria. The degree of overlap, specialized for the two principal CATH codes (class and architectures) is provided as Supporting Material.