MCPath: Monte Carlo path generation approach to predict likely allosteric pathways and functional residues.

Allosteric mechanism of proteins is essential in biomolecular signaling. An important aspect underlying this mechanism is the communication pathways connecting functional residues. Here, a Monte Carlo (MC) path generation approach is proposed and implemented to define likely allosteric pathways through generating an ensemble of maximum probability paths. The protein structure is considered as a network of amino acid residues, and inter-residue interactions are described by an atomistic potential function. PDZ domain structures are presented as case studies. The analysis for bovine rhodopsin and three myosin structures are also provided as supplementary case studies. The suggested pathways and the residues constituting the pathways are maximally probable and mostly agree with the previous studies. Overall, it is demonstrated that the communication pathways could be multiple and intrinsically disposed, and the MC path generation approach provides an effective tool for the prediction of key residues that mediate the allosteric communication in an ensemble of pathways and functionally plausible residues. The MCPath server is available at http://safir.prc.boun.edu.tr/clbet_server.

INTRODUCTION

Allostery is a key concept for regulation of protein activity. Important processes such as binding or catalytic activity can be modulated by perturbations at distal sites such as ligand binding, chemical modification or changes in the environment. The classical models in the ‘old view’ require conformational changes (1–3), yet the ‘new view’ emphasizes the pre-existence of conformational states and dynamics in the allostery (4). Also, the perspective of single common route has been shifted to multiple pathways in allosteric communication between functional sites (5). A perturbation creates a signal that propagates through non-covalent interactions in multiple pathways from the allosteric to the active/binding sites, which is an important aspect underlying the allosteric mechanism and which still remains elusive.

The nuclear magnetic resonance (NMR) is a powerful tool to investigate allosteric pathways experimentally (6–9). The detection of dynamics of allosteric movement provides significant information for the investigation of signal propagation pathways (7). Site-directed mutagenesis has also been useful for the determination of allosteric sites and key residues that play a role in signaling (10,11). Complementing the experimental efforts, various computational approaches have been developed toward understanding allosteric signaling mechanisms (5). In those approaches, energetic connectivity from evolutionary data (12), anisotropic thermal diffusion as a form of energy propagation (13), optimal path generation on residue networks with knowledge-based potentials (14), perturbation response scanning (15), information propagation based on a Markov process (16), protein dynamics from modeled energy landscape (17), geometrical interpretation of protein structures (18), interaction energy networks (19), information theory (20) have been used. Nevertheless, to our knowledge, a few web servers and databases are available for determination of signaling residues such as ALLOSMOD (17) and ASD (21).

Here, we propose and implement a Monte Carlo (MC) path generation method to define all likely allosteric communication pathways by generating an ensemble of maximum probability paths and also an infinitely long path for plausible functional residues based on the graph centrality measures (22,23). The pathways are the assembly of consecutive and probable paths of inter-residue interactions between two allosteric sites. The interaction weights of each residue with the others based on a form of atomistic potential function and the MC sampling provides an ensemble of stochastic pathways rather than a single pathway without a complete enumeration of all paths. This should provide maximally probable routes for a perturbation to be dissipated in the presence of conformational and dynamics states. The proposed method is elaborated mainly on PDZ domain representative proteins, and also tested on bovine rhodopsin and three myosin structures of different conformations that are provided as supplementary cases.

MATERIALS AND METHODS

Energy of inter-residue interactions

A protein can be considered as a network of residues, where the residue–residue interactions are described by a form of the Lennard-Jones 12-6 potential that provides the interaction probability of each residue with its neighbors. The Lennard-Jones 12-6 potential was modified as follows (24): where, E(r) is the Lennard-Jones potential function of inter-atomic interaction with minimum energy ε between two atoms of r distance and collision diameter σ, u is the unit step function, rmin () is the radius of the minimum energy and below which the energy is considered as Emin [] and rcut is the maximum radius (5.5 Å) considered for an interaction above which the energy is taken as 0. The repulsion part in this form of the potential function is modified to avoid the over punishment of the accidental get-togethers in crystal structures. Experimentally obtained van der Waals parameters were used in the calculations.

MC path generation

In the MC path generation, each step of the path is created based on the weighted inter-residue interactions. A step is defined as the direct propagation of a signal between two interacting residues. For a protein with M atoms and N residues, all atom–atom interaction energies E(r) have been calculated based on a given structure and stored in a square energy matrix with M × M dimensions. Then the energy between each atom in every residue pair is summed up to obtain a single energy value E′(r) between those two residues. The atomistic energy matrix thus reduces to the residue energy matrix with dimensions N × N.

The occurrence of an interaction between residue pairs i and j is proportional to the Boltzmann weight and normalized to define a probability measure as follows:

The probability matrix of dimensions N × N comprises the normalized probabilities for the interaction of N residues with the remaining N − 1 residues. The diagonal elements have 0 values. The premise here is that the inter-residue interaction describes the highest probably flow path from the subset of contacts in a given structure for the propagation of a perturbation. With the MC path generation, the random numbers are used to mediate the flow on a subset of steps accessible to each residue according to the weighted inter-residue interactions using Equation (2). For the determination of functional residues, the graph centrality parameters such as closeness and betweenness are calculated. To obtain non-zero values for all residues and to have a proper sampling of paths, infinitely long pathways are generated. The closeness is the inverse of the average of the shortest path between a residue and all other residues given as follows:

The betweenness is the reflector of the frequency of a residue that is a part in all shortest path pairs within the complete structure formulated as follows: where gjk denotes the shortest path between j and k, and gjik denotes a shortest path between j and k and passing through vertex i.

WEB SERVICE

The MCPath server (http://safir.prc.boun.edu.tr/clbet_server) is based on MC path generation runs with a PDB ID or an uploaded structure. The user has three different uses of the web service:

(A) Generating paths starting from a given residue and number of steps: User needs to give the initial residue index and chain ID, and the path length The output lists all the pathways ranked with their probabilities and populated pathways. Populated pathway is referred to a pathway that occurs more than once.

(B) Generating paths between two given residues: User needs to specify the initial and final residue indexes and chain ID. The output lists all the pathways ranked with probabilities and populated pathways.

In the generation of allosteric pathways (A and B), the output also includes the top three populated pathways shown on the structure’s ribbon diagram. If all the pathways are distinct, the three pathways with the highest probabilities are demonstrated. With increasing sizes of proteins, as more number of steps needs to link the two allosteric points, observing a pathway repeating exactly with the same steps of the path becomes rarer. This could be resolved by clustering the pathways with a similarity measure.

(C) Generating a long path with an infinite number of steps for the calculation of graph centrality measures such as closeness and betweenness: User needs to give the path length. The output lists and plots the closeness and betweenness values with the suggested functional residues, which are identified from the maxima of the curves. The path length might depend on the protein’s size. Longer path lengths are suggested for larger protein structures.

The flowchart of the process is given in Figure 1.

Figure 1.

The flow chart of the MCPath server.

RESULTS AND DISCUSSION

PDZ domains are small globular modular proteins that mediate protein–protein interactions (25). These proteins assume intra-domain allosteric activity and extensively been studied by experimental and computation means. Here, the PDZ representative proteins, third PDZ (PDZ3) domain from Post-Synaptic Density Protein 95 (PSD-95) and PDZ2 domain from Human Phosphatase (HPTP1E) are mainly focused. Intra-molecular communication pathways and functional sites are searched by the MCPath Server. The predictions on bovine rhodopsin and three myosin structures of different conformations are also presented as supplementary case studies (Supplementary Figures S1 and S2).

PDZ Domain from PSD-95

The PDZ3 domain from synaptic protein PSD-95 in complex with a C-terminal peptide is selected for the analysis (PDB ID: 1BE9) (26). The populated pathways between the active site His372 at the ligand-binding pocket and the allosteric site Leu353 (12) by the MC path generation using the MCPath server are mapped on the structure in Figure 2. The two distinct pathways that any perturbation would follow to propagate are visualized with lines. The interaction transmits through β2 in the most populated pathway, and it goes through the peptide and then into the core of the domain via β2 and anti-parallel β3 in the third populated pathway. The second populated pathway is a subset of the third pathway. The residues by site-directed mutagenesis (27) highly agree with the MCPath’s results (Figure 2). The present findings also agree with the results of the previous computational methods, energetic connectivity from the evolutionary data (12) and anisotropic thermal diffusion (13). Ile341 previously suggested (13,27) to be involved in the allosteric interaction indeed participate not in the top three but next populated pathways. The most important result here is that the energetic coupling suggests a pathway through the peptide and the others do not, yet the MCPath’s most populated pathways suggest both ways.

Figure 2.

Two possible communication pathways between the ligand-binding site His372 (red) and the allosteric site Leu353 (green) in the PDZ3 domain from PSD-95 (PDB ID: 1BE9), with the C-terminal ligand (yellow). The residues in the allosteric pathway are colored in gray. The most populated pathway (blue) and second populated pathway (magenta) are shown. The second most populated pathway is a subset of third populated pathway. The table lists the three populated pathways and the other computational (12,13) and experimental (27) results.

To identify plausible functional residues, the graph centrality measures are obtained by the MCPath server. The closeness values of the network based on an infinitely long pathway are plotted in Figure 3. The residues at the maxima of the closeness value curve are listed. These residues are highly conserved with the conservation scores 6–9 (28), except Tyr397. The latter residues include the active site residues Leu379 and Phe325 (29) and also the residues observed in the communication pathways in Figure 2. Apparently, the residues with high closeness values are predisposed to be on an allosteric pathway or for another functional role.

Figure 3.

The closeness values of the PDZ domain from PSD-95. The plausible functional residues suggested by the closeness values (maxima of the curve) are listed.

Different runs (of the same or varying lengths) might result in a slight shift in the maxima values, nevertheless the overall shape of the closeness/betweenness curves give consistently the same regions of the structure for high centrality.

PDZ2 domain from human phosphatase

The PDZ2 domain from human phosphatase HPTP1E in complex with a C-terminal peptide is also studied here (PDB ID: 1D5G) (30). Although the structures of PSD-95 and hPTP1E are similar, the analysis may result in different pathways. The analysis was performed on the first NMR model structure. The populated pathways between the peptide-binding site His71 and the allosteric site Val85 (12) are investigated by the MCPath server. His71 is a characteristic binding residue for PDZ2 domains (30). The most popular communication pathways by the MCPath server are mapped on the structure in Figure 4a. One of the most populated pathways is a relatively short pathway passing through α2.The second populated pathway is longer and passes through the core of protein structure via β2 in the N-terminal region. When an ensemble of pathways of 10 steps long from His71 is generated without the end (target) point specificity, the resulting populated pathways are moving through either α2 or the protein core via a group of β turns in the C-terminal region to β2 (Figure 4b). The pathways passing through protein core nevertheless need more steps to reach the allosteric site Val85, which is targeted in the guided pathway generations. On the other hand, when an ensemble of pathways with eight steps long are generated starting from Ile20, another key peptide-binding site (29), two populated pathways reaching to the two distal sites, Val40 and Thr81, are obtained (Figure 4c).

Figure 4.

(a) Two possible communication pathways between the active site His71 (red) and the allosteric site Val85 (yellow) (12) in the PDZ2 domain from HPTP1e (PDB ID: 1D5G): the most populated (magenta) and the second populated (blue) pathways. (b) Two populated pathways from an ensemble of 10-step paths generated starting from His71 (red): the most populated (magenta) and the second populated (blue) pathways. (c) Two populated pathways from an ensemble of 8-step paths starting from Ile20 (red): the most populated pathway (magenta) that reaches to the distal site Thr81 and the second populated pathway (blue) ending the distal site Val40. (d) Experimentally verified energetically coupled residues (magenta) (12) and the residues suggested for the allosteric signaling based on the methyl-containing side chain dynamics (blue) (9). The binding site (His71) is in red. C-terminal ligand is shown in yellow stick.

Figure 4d shows the experimentally identified allosteric pathways: changes in the side chain methyl dynamics with the peptide binding (31) and energetic connectivity by using evolutionary data that were experimentally validated (12). The generated pathways by the MCPath server agree with the latter studies; pathways starting from the two key peptide-binding sites encapsulate the allosteric interactions described in those studies. Two distal allosteric sites, Val40 and Thr81, are identified as in the previous NMR studies (9,10,20,32). In this analysis, the calculations were performed without the peptide because the high affinity of signal transmission through peptide may block the prediction of other signaling pathways.

The closeness values by the centrality analyses are plotted and the functional residues are listed in Figure 5. The results are agreed by the NMR study (31), where Val40 and Leu66 were suggested as distal residues affected by the peptide binding. Leu18 and Ile20 are the binding sites. Ile35 was identified as interacting with both distal sites using information theoretical method (20). Also, the residues suggested by the closeness values highly cover the residues in the allosteric pathways presented in Figure 4. The suggested functional residues have high conservation scores of 5–9 (28).

Figure 5.

The closeness values for the PDZ2 domain from HPTP1e. The functional residues suggested by the closeness values (maxima of the curve) are listed.

CONCLUSION

The MC Path generation creates a new network on a given protein structure composed of the paths created based on the weighted inter-residue interactions. Here we demonstrate that the maximum probability pathways with the stochastic nature of the MCPath generation provide possible propagation routes for a perturbation. The analysis is based on a static structure, yet it may still be used for the cases involving large conformational changes by providing a set of most probable routes with commonalities. Populated pathways imply the predisposition of multiple pathways, which could be used in one or in different functional states of the structure. The centrality parameters also provided by the MCPath server labels the plausible functional regions of the structure, such as active sites or catalytic sites and also the regions with a role in allosteric communication.

SUPPLEMENTARY DATA

Supplementary Data are available at NAR Online: Supplementary Figures 1–2 and Supplementary References [33-40].

FUNDING

TUBITAK Project 110T088 and Betil Fund. Funding for open access charge: Betil fund.

Conflict of interest statement. None declared.