The PTPS crystal structure of the P. falciparum was retrieved from the online Protein Data Bank (PDB ID: 1Y13). The protein functional unit consists of six identical chains, a total of 978 residues. This structure was used to construct an ENM on the C_β atomic coordinates of the protein, as found in the PDB file. The C_β atoms were selected because it has been shown that they provide a better representation of the side chains orientation. A harmonic potential within a cut-off distance (Rc) of 15 Å was used to account for the pairwise interactions between all of the C_β atoms.

For motif analysis, the protein sequences of PTPS from 20 different species, including 9 Plasmodium species, four mammalian species, four bacteria species and three fungi species (Supplementary Table S1) were retrieved using the P. falciparum PTPS sequence as a query in a BLAST search (Altschul et al., 1990) to identify other Plasmodium homologs in the PlasmoDB (Aurrecoechea et al., 2009) and mammalian, bacterial and fungi homologs in UniProt (Bateman, 2019). The BLAST search tool was used with default alignment parameters.

NMA was predominantly performed using the MODE-TASK tool suite (Ross et al., 2018) which employs the ANM (Atilgan et al., 2001) to construct an ENM of the protein structure. Although the ANM has been described previously (Atilgan et al., 2001), we provide a summary here. In the ANM, each residue of the P. falciparum PTPS structure was represented by a node placed at its C_β atom coordinate and all interactions between each pair of residues separated within a defined cutoff distance, Rc = 15 Å, was modeled as a set of springs with a uniform force constant γ. This yielded a network that contained N nodes and M springs representing the total number of interactions defined in the network, such that any given pair of nodes within Rc of each other will interact in accord with a conventional harmonic potential. The normal modes of this network are calculated in the absence of an external force, under an equilibrium condition. In the general case of N residues connected by M springs, the Hessian Matrix H is a 3N × 3N super-matrix that may be derived from the second derivatives of the overall potential V, with respect to the components of R_i, where R_j are the fluctuation vectors of the individual residues. Therefore, the Hessian matrix describes the force constant of the system. H is composed of N × N super-elements, i.e.,

where each super-element H_ij is a 3 × 3 matrix that holds the anisotropic information regarding the orientation of nodes i,j The ij^th super-element (i≠j) of H is defined as:

At equilibrium, the second derivatives may be calculated for the ANM using the β-carbon position vectors of PDB structures such that the elements of the off-diagonal H_ij are given by the equation:

And the elements of the diagonal super-elements H_ii are given by the equations:

For the diagonal elements of H_ii and

For the off-diagonal elements of H_ii

H can be decomposed into 3N-6 eigenvalues and 3N-6 eigenvectors that correspond to the respective frequencies and directions of the individual non-trivial modes. The modes with the lowest frequencies are termed the slowest modes and define the most collective, or global, motions of the protein. The highest frequency modes describe the more localized motions of the protein. The VMD program (Humphrey et al., 1996) was used to visualize the eigenvectors describing the structural change in each mode. Lastly, mechanical stiffness calculations (Eyal and Bahar, 2008) were performed on the ANM of the PTPS using the ProDy Python package (Bakan et al., 2011, 2014) and the obtained results were visualized using Matplotlib library and VMD.

In addition, the Gaussian Network Model (GNM) calculations were performed using the DynOmics portal (Li et al., 2017), specifically to identify hinge regions within the structure of PTPS. The GNM calculations were performed using a cutoff distance of 15 Å, with a spring constant scaling factor of 10 for contact distances ≤4.0 Å and a distance scaling exponent of 2. The Dynomics portal was further used to validate the high-frequency vibrating residues and to analyze the mechanical properties of PTPS, by obtaining a color-coded representation of the PTPS structure based on the mobility of the residues from the resultant low and high GNM modes.

The inverse of H is equivalent to the covariance matrix C that is composed of N × N super-elements. Each off-diagonal, ijth, super-element of H^–1 contains the 3 × 3 matrix of correlations between the x-, y- and z-components of fluctuation vectors of residues i and j, while the ith super-element of H^–1 describes the self-correlation between the components of fluctuation vectors of residue i. For any given mode, the mean square fluctuations (MSF) of individual residues can be obtained by summing the components of fluctuation. The MSF profiles from individual normal modes were calculated using MODE-TASK (Ross et al., 2018).

The BIO3D R package for the exploratory analysis of the structure and sequence data (Grant et al., 2006) was used to compute the deformation energy and residue cross-correlation. The BIO3D cross-correlation S(i, j) is given by:

where Δr_i and Δr_j are the displacement vectors for atoms i and j, respectively. The elements S(i, j) are stored in matrix form and displayed as a three-dimensional dynamical cross-correlation map. If Sij = 1 the fluctuations of atoms i and j are completely correlated (same period and same phase), if Sij = -1 the fluctuations of atoms i and j are completely anti-correlated (same period and opposite phase), and if Sij = 0 the fluctuations of i and j are not correlated.

The deformation energy signifies the energy density due to deformation as a function of position (Hinsen, 1998) given by Eq. 7. It provides a measurement of the individual atom’s energy contribution toward the structural deformation. PTPS deformation energy was derived from the eigen energies and vectors of the first not trivial 20 normal modes according to the equation:

where r_i, r_j donate the displacement of atom i and j in the mode to be analyzed, Rij⁽⁰⁾ is the pair distance vector (Ri–Rj) in the input arrangement and K is the pair force constant.

Motif discovery was performed using Multiple Expectation Maximisation for Motif Elicitation (MEME) vs 4.11 (Bailey et al., 2015) to identify highly conserved motifs in the Plasmodium PTPS enzymes. A fasta file containing PTPS protein sequences was parsed to the MEME analysis software. The motifs search size was set to the range between 6 and 20 residues. The Motif alignment search tool (MAST) was then used to identify overlapping motifs (Bailey and Gribskov, 1998). A Python script was then used to analyze MAST files and MEME log files. Motif conservation was calculated as a number of sites per the total number of sequences, and the results were displayed as heatmaps. Motifs that were uniquely conserved in Plasmodium species were mapped onto the 3D structure of P. falciparum PTPS and visualized using the PyMOL Molecular Graphics System (DeLano, 2014).

In this study, the ANM was applied to the structure of PTPS to classify its collective motions that have the propensity to lead the enzyme from one conformational state to another, and thus modulate its function. Further investigation of the residue fluctuations within the modes was performed to identify highly active regions that may drive these motions. The nature of allosteric modulation requires a high degree of collectivity which is often well-described by the low-frequency modes (Bahar et al., 2010). The NMA of PTPS (a 3-fold symmetry structure) yielded a total of 2,934 modes. The first 20 slowest non-trivial modes were selected for the characterization of the global intrinsic motions of the protein structure.

Normal modes obtained from symmetrical structures are highly susceptible to degeneracy, thus producing degenerate and non-degenerate modes. The degenerate modes share the same frequency and consequently any orthogonal transformation. In contrast, non-degenerate modes characterize unique directions of motions that often capture global meaningful motions that account for large conformational changes. Previous studies showed that dominant conformational changes of complexes were captured within the slowest non-degenerate modes (Atilgan et al., 2001; Chennubhotla et al., 2005; Shrivastava and Bahar, 2006; Wako and Endo, 2011; Isin et al., 2012; Lee et al., 2017; Ross et al., 2018). Here, we identified eight non-degenerate modes of PTPS. These non-degenerate modes exhibited unique eigenvalues (i.e., frequencies) and displayed unique global motions. Table 1 shows the first 20 non-trivial normal modes of PTPS, their associated frequencies, degeneracy levels and the contribution of each mode to the overall motion of the protein. Although the contribution of the individual modes was low and did not reveal a single dominant motion, we suspect that this may be a consequence of the large size of the protein and the extensive degeneracy of the normal modes. Here we have analyzed the first 20 slowest modes which only represent 0.68% of the total modes, yet when combined they account for 11.08% of total motion of the protein. Similarly, the 8 non-degenerate modes only represent 0.27% of the total modes but account for 3.64% of total motion. The displacement vectors of the individual non-degenerate modes were studied as well as the MSF average over the 20 slowest and 20 fastest frequency normal modes.

To visualize the atomic displacements during the collective motions of PTPS, the respective eigenvectors of each of the eight non-degenerate modes were projected onto the structure of PTPS (Figure 2). Movies were constructed by projecting the eigenvectors onto the structure as a set of frames in which the vectors were added to the original atomic coordinates in increasing steps and then visualized using VMD. The atomic displacement during these modes was then examined to identify distinctive motions that are associated with the enzyme’s tunnel gating such as opening, closing or rotation as well as any other notable movements that resulted in structural changes around the tunnel or active site. The eight non-degenerate modes captured coupling movements between the tunnel gating, the N-terminal β-strands, and the central helices thus demonstrating that these regions primarily regulate the global dynamics of the protein. Furthermore, the observed structural changes of the tunnel were often accompanied by structural deformation around the active site region.

The modal analysis captured four distinctive motions that resulted in structural changes around the tunnel and the active site of PTPS. Mode 3 captured an asymmetric global twist of the terminal domains of the two PTPS trimers, characterized by their rotation in opposite directions (Supplementary Movie 1). From visual inspection, it appeared that the twisting of the terminal domains controlled the expansion and narrowing of the tunnel. Most notably, the diameter of the tunnel was smaller when the terminal domains followed a wringing motion with respect to the principal axis. We, therefore, suggest that the fluctuations of the terminal domains modulate the tunnel gating and consequently induces structural changes around the active site region. In mode 6, a notable tilt of the entire terminal domains was observed, leading to a side to side motion of the protein structure (Supplementary Movie 2). A stretching motion that resulted in the lateral expansion and contraction of the tunnel was captured in mode 7. From visual inspection, it appeared that this expansion and contraction was largely driven by the extensive movement of the terminal and central helices (Supplementary Movie 3). Mode 11 displayed a tilting motion of two central helices in an opposite direction to the third central helix (Supplementary Movie 4).

The identified motions within the tunnel were often accompanied by structural changes in the active site region, in which the tunnel acted as a connecting vessel to allow the entry of the substrate and provide the active site with flexibility. This further demonstrates how the tunnel gating mechanism exerts such control on the active site and illustrates the functional relevance of this cavity in the catalytic mechanism of the protein.

Visualization of the non-degenerate modes revealed three modes in which the wringing of the N-terminal domains promoted the surface exposure of buried regions within the tunnel. In particular, Mode 14 captured the outward and inward movement of the PTPS tunnel (likened to an engulfing movement) (Supplementary Movie 5). Mode 17 presented a prominent bending of the central helices from side to side, resulting in the protein tightening in the same direction (Supplementary Movie 6). Mode 18 featured a breathing motion that is characterized by the upward and downward movement of the central β-sheets (Supplementary Movie 7). Mode 19 displayed a clockwise rotation of the protein core that is associated with the terminal regions twisting in the opposite direction (Supplementary Movie 8). Within the eight identified non-degenerate modes, residue 76 to 114 of the N-terminal β-strands and residue 136 to 150 of the central helices exhibited the highest mobility (Figure 9A), which propelled the global motions of the protein. The expansion of the tunnel may promote the binding of PTPS to accessory proteins. Other proteins of this family are known to associate with other proteins, for example, the T-fold enzyme GCH1 is known to physically associate with another pentameric enzyme via the central opening of its tunnel in order to regulate its activity (Maita et al., 2004; Higgins and Gross, 2011).

Based on our modal analysis, we suggest that the flexibility of the tunnel is essential for the efficiency of the active site. The tunnel appears to modulate entry of the substrate via expansion and contraction which is primarily driven by the rotation of the N-terminal regions, thus allowing the substrate to pass through the tunnel into the active site for catalysis. Furthermore, we hypothesize that the twist and shear motions propel the substrate through the tunnel. Taken together, these results further the understanding of the functional dynamics of the PTPS enzyme, which will serve as significant guidelines toward the design of allosteric modulators that target structural regions that are pertinent to the function of the enzyme, and therefore the life cycle of the malarial parasite.

The atomic MSF profiles of the eight non-degenerate low-frequency modes were calculated and plotted (Figure 3). In each of the individual profiles, the N-terminal β-strands and the central helices have shown a notable fluctuation when compared to the protein core. Sharp peaks on the MSF profile distinguished these regions (Figure 3). The calculated MSF of the slowest 20 normal modes and the experimental B-factor were compared to validate the fluctuations that were predicted in the normal modes. Figure 4B shows the atomic fluctuation over the first 20 non-trivial low-frequency modes mapped onto the PTPS structure. The combination of the slowest 20 modes identifies the PTPS terminal domains as the most mobile (red) regions of the structure (Figure 4B). The simplified models of the ANM were shown to agree with the experimental B-factor results, as presented in Figure 4.

Deformation is defined as any structural change caused by either an external force or a change in temperature (Bao, 2002). When the applied force is sufficient, a notable deformation can be observed; otherwise, the structure will resist the applied force and revert to its original state. After obtaining the normal modes, an analysis of the deformation energy distributed across the structure was performed to observe which regions of the structure were deformed during these global motions. The obtained results showed that the extensive fluctuations of the N-terminal domains, as seen in the MSF profiles, caused a significant build-up of deformation energy in the tunnel region of the enzyme. This demonstrates that the fluctuation of N-terminal domains has a long-distance effect on the tunnel core, which is located ∼30 Å away. The deformation energy values were calculated as the sum of the contributions of the first 20 non-trivial modes. The derived energy values were mapped onto the PTPS structure (Figure 5). We observed low deformation energy in the fluctuating terminal domains (blue), which was accompanied by high deformation energy within the tunnel region (red).

A heatmap showing the cross-correlation of the C_β atom pairs is shown in Figure 6. Only motions in the first 20 non-trivial modes were selected in the calculation to highlight the residues involved in the protein collective motions. The residue cross-correlation analysis highlighted regions moving in a concerted manner, which therefore illustrates their involvement in the structural dynamics of the enzyme. The off-diagonal elements presented a positive correlation coefficient within the same chain, more specifically the residues of the terminal β-strands displayed correlation motion. Given the fact that PTPS is a multimeric protein, chains within one trimer exhibited a similar direction of motion. Anticorrelated motions were identified across the two dimers which can be understood due to the twisting or wringing of the helices across the trimers. The observed anticorrelation across the two trimers designates the opposing twisting or wringing of the terminals, therefore allowing the tunnel to open and close accordingly.

Protein motifs are evolutionarily conserved sequence patterns that might unearth a biological function. Identifying the motifs and their structural locations can, therefore, provide insight into regions that modulate protein function. It can further reveal key residues that are unique and typically required to retain the protein function and structural stability (Mackenzie and Grigoryan, 2017; Ross et al., 2017; Zheng and Grigoryan, 2017).

Following the analysis of the dynamics of the enzyme, we performed motif discovery to identify highly conserved motifs in P. falciparum PTPS and residues that are potentially involved in maintaining the secondary structure of the enzyme. The motif discovery was performed on the PTPS homologue sequences from four mammalian species including human, four bacterial species, three fungal species, and nine other Plasmodium species. A total of 27 motifs were identified. The analysis revealed four motifs (motifs 5–8) that were conserved in the Plasmodium PTPS enzyme sequences but were not detected in any of the mammalian species. Motifs 6–8 were uniquely conserved in the Plasmodium species, while motif 5 was also found among the bacterial and fungal PTPS enzymes. A heat map illustrating the occurrence of the motifs is shown in Figure 7.

The conserved Plasmodium motifs were mapped onto the crystal structure to identify their location (Figure 8). Motif 5 was located in the central β-sheet strands forming the tunnel cavity, motif 6 and 7 were located in the N-terminal antiparallel β-strands, as well as a loop linking the β-sheet strands and motif 8 was located in the central α-helices region. Motif 7 and 8 were shorter and more conserved relative to other motifs, which further suggests functional importance. As shown by the NMA, these regions were responsible for modulating the PTPS conformational transitions of the tunnel. The results presented significant structural and sequence differences between the Plasmodium PTPS enzyme and the human PTPS. As all of the identified motifs were not found in the human PTPS, this vital difference can be exploited for the attainment of drug selectivity and future antimalarial drug design.

The low-frequency normal modes are often associated with large amplitude conformational changes which are essential for function (Mahajan and Sanejouand, 2015). The first 20 modes illustrate the substantial contribution of the PTPS helices to the protein’s global motion, in which they regularly control the tunnel movement via wringing or bending motions. In the MSF profile, the most substantial contributions to the atomic fluctuations emerged from the terminal and central helices (Figure 9A). The residues G35, K97, N107, and S145 showed notable fluctuations in the MSF profile. This demonstrates their key role in the global dynamics of the enzyme. Notably, residue N107 was found in motif 7 (position 9) and S145 was found in motif 8 (position 10) (Figure 8), further highlighting residues that are highly conserved in the Plasmodium species and which were active in modes that captured notable structural changes in the enzyme.

As the frequency increases, the modes become more localized and are accompanied by fast vibration of individual residues. It has previously been reported that high-frequency vibrating residues overlap with residues that are highly conserved and have an essential role in stabilizing the protein, and thus maintaining its function (Haliloglu et al., 2005). High-frequency vibrating residues in PTPS were identified from the average MSF calculated over the 20 highest frequency modes (Figure 9B). These included the active site residues H29, H41, and H43 as well as other important residues for the enzyme catalysis. The residues L78, I83, V159, E161, and A167 showed high fluctuation in the MSF profile of the fastest modes. Motif analysis showed that these residues do in fact overlap with sites of high conservation, with L78 and I83 found in motif 6 at (position 3 and 8), V159, E161, and A167 in motif 5 (position 6, 8, and 14), respectively (Figure 8).

Previous studies established that a cluster of hydrophobic residues surrounds the active site pocket and interacts with the substrate ring (Bürgisser et al., 1995; Colloc’h et al., 2000). Here we show that several of these hydrophobic residues also fluctuated in the high-frequency modes (Figure 10B). Residues E161 and T127, which are located at the bottom of the active site pocket (Figure 10A), showed significant fluctuation of the fast frequency modes, with E161 displaying the highest residue fluctuation (Figure 9B). These two residues were previously reported for their key role in substrate recognition and binding, as they both act as proton donors and acceptors during catalysis (Nar et al., 1994; Bürgisser et al., 1995; Ploom et al., 1999; Nar, 2011). Furthermore, Nar and colleagues reported that T105, T106, and E107, located around the active site pocket, constitute an acceptor site for the substrate ring during catalysis in the Rattus norvegicus. PTPS structure (PDB ID: 1B6Z). In the P. falciparum PTPS structure these residues are equivalent to S126, T127, E129, all three of which corresponded to high-frequency vibrating residues. Overall, the combined results obtained from our NMA and motif analysis have revealed conserved residues in the Plasmodium species that have key structural significance, with strong potential to modulate the stability and function of the enzyme. Thus, the characterized regions can then be proposed as alternatives that can be targeted and kept in consideration for future drug design efforts.

Hinge regions are often found between domains in a protein to allow flexibility, in which they permit the domains to move relative to one another or clamp down on a substrate (Towler et al., 2004; Yang and Bahar, 2005; Amusengeri and Tastan Bishop, 2019). The hinge residues, act as anchoring points and are involved in the propagation of large-scale conformational changes. The ANM provides a 3D description of motions identified in the modes, whereas in GNM the motion is projected to a mode space of N dimensions and therefore provides a description of atoms’ mean squared displacements. GNM calculations are considered to be the preferred method for predicting the magnitude of motions at the cost of losing directions (Atilgan et al., 2001) and thus easily allow the identification of hinge sites. The DynOmics portal was used to recognize hinge residues (Li et al., 2017). The server provided information about key sites involved in the collective mechanics and allostery of the protein by mapping of the coarse-grained conformations driven by collective modes to their full-atomic representations. The domain separations analysis based on modes of the GNM disclosed hinge sites (residues that exhibit minimal displacements in the softest two modes). The labeled residues in Figure 11 correspond to the interfacial residues whose neighboring residues have a different sign in the eigenvector relative to themselves. Thus, these interfacial residues act as hinges about which their neighboring residues are displaced in opposite directions. Specifically, we identified hinge sites that overlapped with residues of the active site: H29, H41, and H43 as well as E161. The functional importance of these residues was supported by previous studies in which the mutation of these residues resulted in either a complete or a dramatic loss of enzymatic activity (Bürgisser et al., 1995; Nar, 2011). Furthermore, the active residues were responsible for the coordination of the catalytically important metal ion (Khairallah et al., 2020). The residues Y45, L49, T127, S148, V159, E161 (as well as their neighboring residues) also corresponded to high frequency vibrating residues that were identified in the ANM (Figure 9B), suggesting that high-frequency vibrating residues correspond to rigid points about which significant conformational changes occur.

To further validate the high-frequency vibrating residues that were identified in the ANM, we used the DynOmics portal to analyze the mechanical properties of PTPS to identify the most rigid residues in the GNM (Figure 12). The results obtained from the DynOmics portal were in agreement with our findings from the ANM. Firstly, the N-terminal and central helices were predicted to have high mobility in the low-frequency modes, while the most rigid residues were located in the active site region (Figure 12A). Secondly, the GNM calculations performed using DynOmics located high energy hotspots in the active site region which directly overlapped with the high-frequency vibrating residues that were identified in the ANM (Figure 12B).

The ProDy Python package for protein structural dynamics analysis (Bakan et al., 2011, 2014) was also used, in particular, the ProDy MechStiff function to evaluate the mechanical stiffness of the PTPS protein (Eyal and Bahar, 2008). A mechanical stiffness matrix was then produced by calculating the effective force constants in response to uniaxial extensional forces exerted at each pair of residues (Eyal and Bahar, 2008; Figure 12C). The obtained results illustrated that pairs located in the active site and the surrounding region exhibited strong pairs of interactions. This suggests that the residues of the active site H29, H41 and H43 as well as its surrounding catalytic residues including T127 and E161 bear a relatively strong resistance to deformation. Figure 12E shows the results averaged over all pairs for each residue, which provides a profile of the mechanical resistance of individual residues to deformation. Some residues, especially those located in the active site and surrounding area were more rigid than others and more resistant to deformation, as indicated by this profile. More specifically, these residues involved the active site residues H29, H41, and H43 as well as S17, L78, I83, F93, T127, G131, and E161 which were also among the highly conserved, high frequency vibrating and hinge residues. A cartoon representation of a single chain of the P. falciparum PTPS is shown in Figure 12D. The figure shows pairs of residues with the strongest effective force constant and their location within the protein structure. The identified pairs of T127, E16, G131 were located around the active site area and were also among the high vibrating and hinge residues.