Substrate Effect on Catalytic Loop and Global Dynamics of Triosephosphate Isomerase

Abstract

The opening/closure of the catalytic loop 6 over the active site in apo triosephosphate isomerase (TIM) has been previously shown to be driven by the global motions of the enzyme, specifically the counter-clockwise rotation of the subunits. In this work, the effect of the substrate dihydroxyacetone phosphate (DHAP) on TIM dynamics is assessed using two apo and two DHAP-bound molecular dynamics (MD) trajectories (each 60 ns long). Multiple events of catalytic loop opening/closure take place during 60 ns runs for both apo TIM and its DHAP-complex. However, counter-clockwise rotation observed in apo TIM is suppressed and bending-type motions are linked to loop dynamics in the presence of DHAP. Bound DHAP molecules also reduce the overall mobility of the enzyme and change the pattern of orientational cross-correlations, mostly those within each subunit. The fluctuations of pseudodihedral angles of the loop 6 residues are enhanced towards the C-terminus, when DHAP is bound at the active site.

1. Introduction

Flexibility and dynamics are crucial for the biological function of proteins. In certain enzymes, conformational changes, such as transformation between open and closed structures, are linked with substrate binding and product release to/from the active site. In this work, we will focus on the widely-studied [1] triosephosphate isomerase (TIM), which is an essential enzyme in the glycolytic pathway.

TIM is a highly efficient enzyme, catalyzing the interconversion of dihydroxyacetone phosphate (DHAP) and D-glyceraldehyde 3-phosphate (GAP). It is a homodimer, with each subunit adopting the TIM-barrel fold. TIM is fully active as a dimer, even though there is no coupling between the two catalytic sites. In fact, each catalytic site is composed of residues N11, K13, H95 and E165 (numbered based on chicken TIM) that belong to a single subunit. Residue E165 is located on the so-called catalytic loop 6 (E165-A176). Closure of loop 6 over the active site is necessary during catalysis so that exposure of the active site to solvent and thereby formation of a potentially toxic byproduct are inhibited [2]. However, experiments have indicated that loop 6 closure is not ligand-gated, meaning that it is also observed in the apo state [3].

Alignment of the available crystal structures for TIM (both apo and ligand-bound) predominantly discloses the conformational changes of loop 6 between open and closed states. Thus, earlier computational studies have focused only on the loop 6 region flexibility and dynamics [4,5,6,7]. Through normal mode analysis based on elastic network model (ENM), global motions of the intact dimeric TIM [8] have been revealed for the first time. More importantly, the opening/closure of loop 6 has been found to be coupled with certain collective modes, both in ENM and molecular dynamics (MD) studies of apo TIM dimer [8,9].

More recently, several studies have appeared that emphasize such coupling between energetically-favored global modes of motion and functional loop dynamics [10,11,12] for a larger set of proteins. Specifically, for a set of 10 enzymes with different sizes and oligomerization states, it has been shown through ENM that the experimentally determined catalytic loop (10–20 residue-long) reconfigurations are driven, at least to a certain extent, by collective modes of the intact enzyme in the apo state [10]. It must be noted that ENM is a purely entropic model based on non-specific harmonic interactions between closely spaced residue pairs in the protein. Inasmuch as, the entropic effect seems significant for guiding the loop in the apo enzyme, the enthalpic effects due to any bound ligand(s) should also be assessed in terms of loop closure. For this aim, classical MD simulations seem to be appropriate for incorporating specific electrostatic interactions.

Having established the relation between global modes and functional loop dynamics, here we will elaborate on the dynamics of TIM with the substrate DHAP bound in one or both of the catalytic sites through MD simulations. The enzyme-DHAP complex is predominantly detected by experiments [13]. In previous ENM studies [8], an almost one-to-one correspondence has been observed between the collective modes of two different TIM crystal structures, one being the apo state with the catalytic loops open and the other ligand-bound with closed loops. In our current study, we performed 60-ns independent MD simulations on apo TIM and its complex with DHAP, so as to perform a detailed comparison of their structural dynamics with the inclusion of enthalpic effects.

2. Results and Discussion

2.1. RMSD and Residue MSF

Four independent 60 ns runs are utilized for sampling the conformational space of the apo and DHAP-bound enzyme more effectively. The root mean square deviations (RMSD) of the snapshots with respect to the initial (energetically minimized) structure are plotted in Figure S1 for each trajectory. The RMSD calculation is carried out based on C^α atoms after alignment of each snapshot onto the initial structure. RMSD profiles fall within reasonable limits (2–3 Å), implying the reliability of each trajectory. Based on these profiles, first 10 ns of each simulation is discarded as equilibration period and the remaining 10–60 ns range is considered in the analysis of time-averaged properties. In Complex1 simulation, both DHAP molecules stay bound to the catalytic sites located in each subunit. In Complex2, DHAP in monomer A detaches from the catalytic site at around 13^th ns (just after equilibration period), whereas DHAP in the catalytic site of monomer B stays bound throughout the same run. For this reason, averaged properties for the complex are obtained using Complex1 A and B and Complex2 B subunits. In the case of apo simulations, A and B subunits from both runs are considered in the averages.

Mean square fluctuations (MSF) of alpha-carbon atoms of each residue are calculated based on the average structure of each trajectory. Average profiles for the apo and complex are plotted in Figure 1. The overall mobility of the enzyme is constrained in the presence of the DHAP. Interestingly, this effect is not apparent in residues that are in close contact with DHAP (see inset in Figure 1). In contrast, mobility of relatively distant residues that include the interface between subunits and some outlying helices (shown in blue in the inset) are effectively reduced. Specifically, only three residues (shown in red) that are located on loop 6 exhibit higher MSF in the complex: P166 at the N-terminus (23% increase in MSF), T175 and A176 at the C-terminus (~45%).

Figure 1. Residue MSF for apo and complex simulations. Inset figure shows the complex residues with reduced (blue) and increased (red) MSF by 20% compared to apo. Green regions have similar MSF in both apo and complex simulations. DHAP is shown with stick representation.

Figure 1. Residue MSF for apo and complex simulations. Inset figure shows the complex residues with reduced (blue) and increased (red) MSF by 20% compared to apo. Green regions have similar MSF in both apo and complex simulations. DHAP is shown with stick representation.

During the simulations, DHAP molecules make H-bonding with certain residues: N11, K13, H95, E165, G171, S211, G232 and G233. In this list, the first four are the catalytic residues. G171 is located at the tip of loop 6, S211 is on loop 7 and residues G232 and G233 are on loop 8. All of these residues have been reported as H-bonding contributor residues to the active site geometry [1]. In the case of Complex2 monomer A, H-bonds between DHAP and residues N11, H95 are no longer observed after 9th ns as seen in

2b. This is followed by the breakage of H-bonds with E165, G233. Finally at 13 ns, H-bonding between DHAP and K13, S211 is broken, which leads to the detachment of the substrate from the active site. The number of H-bonds made between TIM residues and DHAP throughout each simulation is given in Figure S2.

2.2. Loop 6 Opening/Closure is Observed in the Presence of DHAP

The open/closed states of loop 6 can be detected looking at the C^α-C^α distances between loop 6 tip residues (I170, G171, T172) and residue Y208 [7,9]. These distances are plotted as a function of time for both chains in each simulation as in Figure 2, Figures S3–4. There is no correlation between the A and B subunits for the specific distances reported. Solid (lower) and dashed (upper) horizontal lines in Figure 2 refer to the I170-Y208 distances in the ligand-bound/closed (1TPH, smallest distance among the two chains) and apo/open (8TIM, largest distance) X-ray structures, respectively. In previous MD studies on apo TIM from chicken [7,9] and Trypanosoma cruzi [10], multiple opening/closure events of loop 6 have been observed during independent 60–100 ns runs.

Figure 2. The distance between C^α atoms of loop 6 tip residue (I170) and a reference residue (Y208) as a function of time. Left and right panels summarize the profiles for apo and complex simulations. Solid and dashed lines indicate the respective closed and open states of the loop based on crystal structures. Loop 6 samples open conformations even in the presence of DHAP. It opens/closes multiple times during apo and complex simulations.

Figure 2. The distance between C^α atoms of loop 6 tip residue (I170) and a reference residue (Y208) as a function of time. Left and right panels summarize the profiles for apo and complex simulations. Solid and dashed lines indicate the respective closed and open states of the loop based on crystal structures. Loop 6 samples open conformations even in the presence of DHAP. It opens/closes multiple times during apo and complex simulations.

It is known that for effective catalysis this loop should stay closed so that the active site is shielded from solvent [1,2]. In the presence of DHAP, the catalytic loop still opens and closes several times during our independent 60 ns runs. Concurrently, either one or both DHAP molecules stay bound at the active site. This interesting situation indicates that electrostatic interactions between substrate and loop 6 are not strong enough to keep the catalytic loop closed over the active site for extended periods. In fact, this could provide a means for effective product release (DHAP), which is not the rate-limiting step in the conversion of GAP to DHAP [1,14]. Also in conformity with our current findings, there are two crystal structures reported, in which the loop 6 is in open conformation despite the presence of a ligand [1,15,16].

Furthermore, the distance histograms for the three tip residues on loop 6 are calculated using all recorded snapshots for apo (left panels) and DHAP-bound (right panels) subunits in Figure 3. In fact, larger distances are sampled in the complex compared to apo simulations. Multiple peaks are observed around 11.5, 12.5 and 14 Å for I170-Y208, 15 and 17.5 Å for G171-Y208, 14.5, 17.5 and 18.5 Å for T172-Y208 distances in complex simulation, whereas the values for apo cases are around 12 Å for I170-Y208, 11.5 and 15 Å for G171-Y208 and 17 Å for T172-Y208 distances.

2.3. Essential Dynamics

In order to investigate the effects of DHAP on the enzyme dynamics, PCA is applied on each trajectory, as explained in Materials and Methods. Contribution of first five PCs to overall motion (percentage variance captured by PCs) is given in Table S1. Table 1 lists the overlap of loop 6 motion with the loop closure direction in the essential modes for both chains.

In apo simulations, PC1 corresponds to counter-clockwise rotation of subunits, which contributes 34% and 19% to overall motion in Apo1 and Apo2 simulations, respectively. In apo simulations, this motion is coupled with loop 6 opening/closing motion as seen in Figure 4a,b and Table 1 (more stressed in Apo1). In TIM-DHAP complexes, this counter-clockwise rotation does not appear in the essential modes for the subunits with DHAP. Still in some PCs, a similar type of rotation is detected as seen in the PC overlap figures (Figure S5) and Figure 4, but it lacks full coordination possibly due to the presence of DHAP molecules.

Table 1. Loop 6 overlap with closure direction (overall PCA).

PC	Apo1		Apo2		Complex1		Complex2
	A	B	A	B	A	B	A	B
1	0.76	0.59	0.37	0.44	0.58	0.68	0.75	0.18
2	0.42	0.04	0.24	0.14	0.76	0.61	0.74	0.10
3	0.14	0.08	0.06	0.58	0.17	0.75	0.22	0.00
4	0.78	0.81	0.49	0.26	0.28	0.49	0.81	0.54
5	0.41	0.16	0.24	0.48	0.55	0.11	0.25	0.54

Table 1. Loop 6 overlap with closure direction (overall PCA).

PC	Apo1		Apo2		Complex1		Complex2
	A	B	A	B	A	B	A	B
1	0.76	0.59	0.37	0.44	0.58	0.68	0.75	0.18
2	0.42	0.04	0.24	0.14	0.76	0.61	0.74	0.10
3	0.14	0.08	0.06	0.58	0.17	0.75	0.22	0.00
4	0.78	0.81	0.49	0.26	0.28	0.49	0.81	0.54
5	0.41	0.16	0.24	0.48	0.55	0.11	0.25	0.54

In Complex1, where both DHAP molecules stay bound throughout the simulation, bending of the two subunits with respect to the interface is emphasized in the lower indexed PCs, namely PC1 and PC2. These modes, shown in Figure 4 and Figure 5, are coupled to loop opening/closing (see Table 1). Complex2 with one DHAP (bound to subunit B) represents an intermediate case between apo and Complex1 and this is reflected in its essential modes. Specifically, PC1 and PC2 of the A subunit (without DHAP) exhibit rotation and bending motion, respectively, as in the case of Apo1. In contrast, the PC2 and PC3 for the bound subunit are bending, where loop opening/closure is not observed. PC4 for the same simulation is the combination of these two motions as seen in Figure 5.

Figure 3. Distance histograms for I170-Y208, G171-Y208 and T172-Y208 for (a, c, e) apo and (b, d, f) complex simulations, respectively. Distance distribution is obtained by dividing the number of occurrence of the distance to the total number of snapshots for 10–60 ns range.

Figure 3. Distance histograms for I170-Y208, G171-Y208 and T172-Y208 for (a, c, e) apo and (b, d, f) complex simulations, respectively. Distance distribution is obtained by dividing the number of occurrence of the distance to the total number of snapshots for 10–60 ns range.

Moreover, PC2 of Apo1 simulation corresponds to bending type motion for both subunits and this motion is also observed in PC3 of Apo2 simulation for monomer A only and PC4 of the same simulation for both subunits. From Table 1, coupling of loop opening/closing motion with bending type motion is not stressed in apo simulations. In summary, the presence of DHAP hinders the coupling of counter-clockwise rotation with loop opening/closing. But loop opening/closing can still be observed in the essential modes, mostly coupled with bending in the presence of two DHAP molecules.

Figure 4. Most dominant dynamics from PC1 of (a) Apo1, (b) Apo2, (c) Complex1 and (d) Complex2. Subunit A (orange), subunit B (blue) with loop 6 (red) are shown in ribbon and DHAP in magenta-orange stick representation. In (a) and (b), PC1 corresponds to counter-clockwise rotation (CCR) of two subunits, coupled with loop 6 opening/closing. In (c), PC1 of Complex1 simulation is a bending type motion rather than counter-clockwise rotation but loop 6 opening/closing is observed in this mode. In (d), subunit A (without DHAP) rotates like in the Apo1 PC1 coupled with loop opening/closing, whereas a mixed-type of motion is observed in subunit B (with DHAP), where the helix 5 undergoes bending without loop 6 opening/closing.

Figure 4. Most dominant dynamics from PC1 of (a) Apo1, (b) Apo2, (c) Complex1 and (d) Complex2. Subunit A (orange), subunit B (blue) with loop 6 (red) are shown in ribbon and DHAP in magenta-orange stick representation. In (a) and (b), PC1 corresponds to counter-clockwise rotation (CCR) of two subunits, coupled with loop 6 opening/closing. In (c), PC1 of Complex1 simulation is a bending type motion rather than counter-clockwise rotation but loop 6 opening/closing is observed in this mode. In (d), subunit A (without DHAP) rotates like in the Apo1 PC1 coupled with loop opening/closing, whereas a mixed-type of motion is observed in subunit B (with DHAP), where the helix 5 undergoes bending without loop 6 opening/closing.

Figure 5. Other PCs for Apo1, Apo2, Complex1 and Complex2 simulations.

Figure 5. Other PCs for Apo1, Apo2, Complex1 and Complex2 simulations.

Normalized orientational cross-correlations between residue fluctuations (Equation 1) are calculated for apo and complex simulations to observe whether or not DHAP affects the correlations within and between the subunits. The results are averaged over the subunits for the symmetric cases (Apo1-Apo2 and Complex1) and given as intra and inter-subunit maps in Figure 6. Corresponding maps for each run without averaging over the subunits are given as supplementary information in Figure S6. Even though there are some differences, the overall features of Apo1 and Apo2 seem consistent.

Figure 6. Normalized orientational cross-correlation maps for apo and complex1. The results for apo are averaged over all the subunits in Apo1 and Apo2 simulations and for complex, over Complex1 A and B subunits. The correlated regions with loop 6 (L6, represented by the average of three tip residues I170, G171, T172) are given in cartoon representations between intra and inter-subunit maps and colored using blue (negative correlation)-white (no correlation)-red (positive correlation) coloring scheme. In the complex intra-subunit map, two blocks- Regions 1 and 2- are negatively correlated with each other. Region 1 is also positively correlated with the same region of the other subunit (named as Region 3) in the inter-subunit map.

Figure 6. Normalized orientational cross-correlation maps for apo and complex1. The results for apo are averaged over all the subunits in Apo1 and Apo2 simulations and for complex, over Complex1 A and B subunits. The correlated regions with loop 6 (L6, represented by the average of three tip residues I170, G171, T172) are given in cartoon representations between intra and inter-subunit maps and colored using blue (negative correlation)-white (no correlation)-red (positive correlation) coloring scheme. In the complex intra-subunit map, two blocks- Regions 1 and 2- are negatively correlated with each other. Region 1 is also positively correlated with the same region of the other subunit (named as Region 3) in the inter-subunit map.

In general, strong positive and negative correlations are observed both in apo and complex cases. However, presence of DHAP leads to a change in the regions correlated with loop 6, which are shown in blue-white-red color code on the cartoon representations in the middle panels of Figure 6. Blue and red colors indicate negative and positively correlated regions with loop tip residues, respectively. In apo, loop 6 is correlated with the regions including helix 1 and helix 8, whereas it is anti-correlated with the same regions in the complex. The correlation of loop 6 with helix 6 C-terminus is negative in the absence of DHAP and turns into slightly positive in the complex. These loop 6 correlations are shown in more detail in Figures S7-S9 for the N-terminus (P166-W168), tip residues (I170-T172) and C-terminus (K174-A176), respectively.

Moreover, a “block” type behavior becomes distinct in the intra-correlations of the complex, where the blocks labeled as “Region 1” and “Region 2” are shown on the cartoon representation of the complex. Region 1 and 2 exhibit clear negative correlations in the intra-subunit map. And in the inter-subunit map, Region 1 shows strong positive correlations with Region 3 (corresponding to Region 1 but on the other subunit). In overall, DHAP presence affects both intra- and inter-subunit correlations in the enzyme.

2.4. Loop Dynamics and Flexibility

To assess the effect of DHAP on loop 6 dynamics, PCA is applied on the loop 6 region only, extracted independently for chains A and B. From Table 2, loop opening/closing motion is the dominant motion (PC1) in Apo1 (59% and 58% for Apo1 A-B), Complex1 (45%, 58% for A and B, respectively) and Complex2 for A chain (60%). Same motion is observed in the PC2 of Apo2 for both chains (22%, 37% for A and B, respectively) and Complex2 B chain (15%). Lateral movement of the loop is the second dominant motion in Apo1 and Complex1 simulations, and it is PC1 of Apo2 A-B and of the Complex2 B. These two dominant motions of the catalytic loop 6, shown in Figure 7, do not seem affected from the presence of the DHAP. However there are certain differences in the flexibility of the catalytic loop when DHAP is bound to the catalytic site.

Table 2. Percentage variance captured by PC modes of loop 6 and overlap with closure direction

	Apo1		Apo2		Complex1		Complex2
PC	A	B	A	B	A	B	A	B
Percentage Variance Captured
1	59	58	45	37	45	58	60	58
2	21	17	22	29	19	18	16	15
Overlap with Closure
1	0.79	0.71	0.09	0.12	0.68	0.73	0.84	0.03
2	0.03	0.49	0.83	0.75	0.51	0.06	0.15	0.75

Table 2. Percentage variance captured by PC modes of loop 6 and overlap with closure direction

	Apo1		Apo2		Complex1		Complex2
PC	A	B	A	B	A	B	A	B
Percentage Variance Captured
1	59	58	45	37	45	58	60	58
2	21	17	22	29	19	18	16	15
Overlap with Closure
1	0.79	0.71	0.09	0.12	0.68	0.73	0.84	0.03
2	0.03	0.49	0.83	0.75	0.51	0.06	0.15	0.75

In order to investigate the effect of DHAP molecules on catalytic loop flexibility, the pseudodihedral angle of residue, φ_i, is calculated using C^α coordinates of residues i−1, i, i+1 and i+2. RMS fluctuations of φ_i for the catalytic loop residues are first calculated for A and B subunits in each run. Then average pseudodihedral RMSF for different subunits and runs are calculated and plotted in Figure 8. In both apo and complex simulations, residues W168-A176 on the loop 6 are more flexible than the others, which is consistent with the previous findings [7,9,17]. However, when the complex and apo cases are compared with each other, the flexibility increases in the residues P166, A169 and G171-T177 more than 35% in the presence of DHAP.

Figure 7. Apo1 and Complex1 most dominant PCs from loop 6 A PCA. PC1 corresponds to loop opening/closing and PC2 lateral movement of the loop for both simulations.

Figure 7. Apo1 and Complex1 most dominant PCs from loop 6 A PCA. PC1 corresponds to loop opening/closing and PC2 lateral movement of the loop for both simulations.

Figure 8. Average pseudodihedral RMS fluctuations for loop 6 residues in apo and complex simulations.

Figure 8. Average pseudodihedral RMS fluctuations for loop 6 residues in apo and complex simulations.

3. Materials and Methods

3.1. MD Simulations

MD simulations are carried out for apo and chicken TIM-DHAP complex using AMBER 11 simulation package [18] with ff03 [19] force field. The analysis is performed on 10–60 ns time range of the following four independent 60 ns MD runs of chicken TIM: Apo1 [9] and Apo2 for apo TIM, Complex1 (extended simulation of [20]) and Complex 2 for TIM with DHAP.

X-ray crystal structures for TIM are extracted from Protein Data Bank (PDB [21]). The crystal structure of chicken TIM (PDB id: 8TIM) at 2.50 Å resolution is used as the initial structure for apo chicken TIM simulations. The initial structure for complex TIM simulations is chicken TIM crystal structure (PDB id: 1TPH [22]) at 1.80 Å resolution. Two phosphoglycolohydroxamate (PGH) ligands in the original complex are replaced with DHAP, since PGH is an intermediate analog of DHAP, having similar structure and orientation.

The same simulation protocol as in our previous work on chicken TIM [9] is adopted here. NPT simulations at 300 K and 1 atm are carried out using a truncated octahedron periodic box filled with TIP3P water molecules [23] and neutralized with ions (see Table 3). The simulations are carried out using a time step of 2 fs and SHAKE algorithm [24]. Long-range electrostatic interactions are handled using Ewald summation technique with the particle-mesh method [25], using a cutoff distance of 9 Å. In TIM-DHAP simulations, the same DHAP parameters are used as in Alakent et al.’s work [20]. The simulation details for all runs are summarized in Table 3.

Table 3. Simulation system details.

	Apo1	Apo2	Complex1	Complex2
Total no. of atoms	57532	53395	60667	59185
No. of water molecules	16674	15295	17729	17235
No. of Na+/Cl− ions	8 Cl−	8 Cl−	2 Na+	2 Na+
Box dimensions (Å)	107	106	108	109

Table 3. Simulation system details.

	Apo1	Apo2	Complex1	Complex2
Total no. of atoms	57532	53395	60667	59185
No. of water molecules	16674	15295	17729	17235
No. of Na+/Cl− ions	8 Cl−	8 Cl−	2 Na+	2 Na+
Box dimensions (Å)	107	106	108	109

3.2. PCA

3.2.1. Essential Dynamics

PCA is applied on the MD trajectory of apo and complex enzyme in order to determine functionally-relevant, collective dynamics of the whole protein during the simulation [26]. For each simulation, the trajectory is first superimposed on the initial structure using ptraj tool of AMBER 11 [18]. After the superimposition, mean position for the alpha-carbon atoms is calculated to obtain the average structure for the simulation. The covariance matrix C of size 3N × 3N is constructed using the deviation of alpha-carbon atoms from the average structure. C is diagonalized in order to obtain eigenvalue decomposition as C = P L P^T. The columns of P matrix are the eigenvectors of unit length, corresponding to the principal components of the trajectory and diagonal entries of L matrix contains the eigenvalues in decreasing order. In this text PC1, which is the first column of P, refers to the linear combination of alpha-carbon displacements with maximal variance, and it is succeeded by PC2, PC3, and so on.

3.2.2. Loop 6 Dynamics

In order to observe the dominant motions of the catalytic loop 6 (on either A or B chain) during the simulations, PCA is based on the (3 × 12) coordinates of each catalytic loop (12-residue long loop: E165-A176), which are extracted from the superimposed trajectory of each simulation.

3.3. Overlap Calculations

3.3.1. Overall Overlap of PCs

After aligning each trajectory onto the Apo1 mean structure [27], the inner product of first 10 modes is calculated for the combinations of trajectories (Apo1-Apo2, Apo1-Complex1, etc.). The overall overlap figures for the simulations are given in Figure S5.

3.3.2. Overlap of Loop 6 Motion with Loop 6 Closing Direction

Loop 6 closing direction is the normalized displacement vector from the open (PDB id: 8TIM) to the closed (PDB id: 1TPH) X-ray structure (aligned on each trajectory), calculated for E165-A176 residues of each subunit A and B. The vector for loop 6 in each principal component is obtained by extracting the coordinates of E165-A176 in both chains. After normalizing all the vectors, the dot product with displacement vector is calculated and absolute value of the overlaps is reported in Table 1 and Table 2.

3.4. Normalized Orientational Cross-Correlations Between Residue Fluctuations

In order to determine negative and positively correlated regions within the subunits (intra) and between the subunits (inter), normalized orientational cross-correlations O(i,j) between residue fluctuations are calculated based on reconstructed trajectories using the first five PCs as below

$$\text{O}\left(\text{i},\text{j}\right)=<\Delta {\mathbf{\text{R}}}_{\text{i}}·\Delta {\mathbf{\text{R}}}_{\text{j}}>/{\left[<\Delta {\mathbf{\text{R}}}_{\text{i}}·\Delta {\mathbf{\text{R}}}_{\text{i}}><\Delta {\mathbf{\text{R}}}_{\text{j}}·\Delta {\mathbf{\text{R}}}_{\text{j}}>\right]}^{1/2}$$

(1)

Here ∆R_i is the fluctuation of the alpha-carbon of the residue i in the position vector R_i. The interval for O(i,j) is [−1,1] with −1 and 1 corresponding to strong negative and positively correlated residue pairs, respectively. O(i,j) = 0 means the residue fluctuations are uncorrelated in terms of orientation.

4. Conclusions

In this study, we investigated the effect of DHAP’s presence on TIM’s catalytic loop motions and overall enzyme dynamics by means of MD simulations. In all simulations, multiple opening/closure events of loop 6 take place even in the presence of DHAP. The distances between the catalytic loop tip residues I170, G171, T172 and residue Y208 reach higher values when DHAP is bound, meaning that the loop adopts more open positions compared to the apo enzyme. Moreover, the average pseudodihedral RMSF shows that the conformational flexibility of the loop 6 is increased in the complex, especially toward the C-terminus.

DHAP presence also affects the overall mobility of the enzyme, reflected by the reduced MSF of most residues including the interface. Only three residues located on the catalytic loop (P166 at the N-terminus, T175 and A176 at the C-terminus) show higher mobility in the complex. PCA reveals that in apo enzyme, counter-clockwise rotation -the dominant mode- drives the opening/closure of the catalytic loop, in conformity with our previous findings [8,9,10]. In contrast, this mode is suppressed in the complex and the bending of the subunits appears in the dominant PCs coupled to loop opening/closure. Furthermore, the normalized orientational cross-correlations are altered into two anti-correlated blocks, as a result of ligand binding. Moreover, the correlations of the loop with the rest of the enzyme are also affected to a great extent.

In summary, the specific interactions of DHAP modify the global and loop dynamics to a certain extent, which have not been observed in previous studies based on purely entropic ENM [8]. In another MD study, the frequencies of intra-minimum motions based on short MD runs have been shown to be altered in the presence of DHAP [20]. Still our simulations indicate that multiple loop opening/closure events are present and coupled to the global motions of the enzyme, both in apo and complex form. Therefore, entropic effects seem to have a significant role in the mechanisms of DHAP binding/release in TIM.