Domain Motions and Functionally-Key Residues of l-Alanine Dehydrogenase Revealed by an Elastic Network Model

Mycobacterium tuberculosis l-alanine dehydrogenase (l-MtAlaDH) plays an important role in catalyzing l-alanine to ammonia and pyruvate, which has been considered to be a potential target for tuberculosis treatment. In the present work, the functional domain motions encoded in the structure of l-MtAlaDH were investigated by using the Gaussian network model (GNM) and the anisotropy network model (ANM). The slowest modes for the open-apo and closed-holo structures of the enzyme show that the domain motions have a common hinge axis centered in residues Met133 and Met301. Accompanying the conformational transition, both the 1,4-dihydronicotinamide adenine dinucleotide (NAD)-binding domain (NBD) and the substrate-binding domain (SBD) move in a highly coupled way. The first three slowest modes of ANM exhibit the open-closed, rotation and twist motions of l-MtAlaDH, respectively. The calculation of the fast modes reveals the residues responsible for the stability of the protein, and some of them are involved in the interaction with the ligand. Then, the functionally-important residues relevant to the binding of the ligand were identified by using a thermodynamic method. Our computational results are consistent with the experimental data, which will help us to understand the physical mechanism for the function of l-MtAlaDH.


Introduction
Tuberculosis is still a serious threat to human health worldwide, which is caused by the infection of the Mycobacterium tuberculosis (Mt) [1]. World Health Organization (WHO) has pointed out that at least one-third of the population in the world could be infected with latent Mt [2]. Up to now, the treatment of tuberculosis is still a difficult challenge due to the emergence of the drug-resistant strains of Mt [3]. At present, it is found that many anti-tubercular drugs become ineffective for the treatment of Mt, and thus, the development of new anti-tubercular drugs to treat this disease is urgent [4]. When the pathogens are under nutrient starvation regimes or anaerobic growth conditions, Rv2780 encoding L-alanine dehydrogenase (L-AlaDH) is overexpressed in persistent Mt [5]. The production and activity of L-AlaDH were improved when Mt was transferred from the aerobic to the anaerobic growth conditions [6]. L-AlaDH can catalyze the 1,4-dihydronicotinamide adenine dinucleotide (NAD)-dependent conversion of L-alanine to ammonia and pyruvate or convert pyruvate to L-alanine in the reverse direction. The catalytic mechanism can be expressed as [7]: L-alanine + NAD + + H 2 O é Pyruvate + NADH + H + + NH4 (1) Because the crystallographic structure of human L-AlaDH has not been determined, L-alanine dehydrogenase from Mycobacterium tuberculosis (L-MtAlaDH) is studied here. L-MtAlaDH is a secretory antigen, which is closely relevant to bacterial persistence during the persistent stage of infection [8]. Therefore, L-MtAlaDH is considered to be a potential target for the development of anti-tubercular drugs [9].
Two forms of structures for L-MtAlaDH have been determined in the Protein Data Bank (PDB), including the dimer ( Figure 1C) and the hexamer ( Figure 1D) states. The hexameric molecule is made of three-fold non-crystallographic axes, and two monomers nearby form the fold axis. As shown in Figure 1A, each monomer of L-MtAlaDH is composed of a substrate-binding domain (Residues 1-128 and 309-371) and a NAD-binding domain (Residues 129-308), which are connected by two α-helixes [9]. The NAD-binding domain includes a central seven-stranded β-sheet, as well as several surrounding α-helices. The substrate-binding domain consists of an eight-stranded β-sheet and several α-helices [9]. The crystal structure of the open-apo and closed-holo states of the protein have been determined by X-ray crystallography, which can be downloaded from the Protein Data Bank (PDB) with the access codes 2VOE and 2VHZ, respectively [8,9]. The experimental studies have demonstrated that upon the binding of the coenzyme, L-MtAlaDH undergoes a large-scale conformational transition from the open-apo to the closed-holo states, in which the substrate-binding domain (SBD) rotates by 16˝toward the NAD-binding domain (NBD) [9]. Considering that the conformational transitions between open and closed structures are essential for the function of L-MtAlaDH, many experimental and theoretical studies have focused on the investigation of the functional domain motions and identification of the related key residues of the protein. The open-closed domain motions of L-MtAlaDH upon the binding of the coenzyme have been studied with molecular dynamics (MD) simulations at the atomic level [10]. The simulation results have shown that the conformational transition of L-MtAlaDH is achieved through the open-closed and twisting motions of the NBD and the SBD. Based on the simulation results, it is proposed that two loops (Residues 94-99 and 238-251) are responsible for the binding of NAD, and another loop (Residues 267-293) is crucially responsible for the binding of substrate [10]. The experiments of Agren et al. have shown that L-MtAlaDH completely lost its catalytic activity when the conserved active site residues His96 and Asp270 were mutated. These two residues play an essential role in the enzymatic reaction at the active site and exhibit strong stability in biological evolutionary process [9].
In this work, the domain motions and the related functionally-key residues encoded in the structure of L-MtAlaDH were investigated by using the coarse-grained Gaussian network model (GNM) and the anisotropy elastic network model (ANM) [11][12][13]. The GNM and ANM are simple yet effective methods to explore the intrinsic dynamics of proteins [14][15][16]. It have been proven in numerous application studies that the low-frequency motion modes calculated by GNM and ANM represent the large-scale collective motions usually relevant to the functions of the protein [17][18][19]; whereas the high-frequency motion modes reflect the geometric irregularity in the protein structure, and the residues motivated in these modes are considered to be crucial for the protein stability [20,21]. In the present work, the functional domain motions of L-MtAlaDH were studied through analyzing the slow motion modes of GNM and ANM, and the key residues responsible for the stability of the protein were identified by calculating the fast modes of GNM and ANM. Besides that, the functionally-important residues responsible for substrate binding of the protein were revealed by using a thermodynamic method proposed by Su et al. [22,23]. The bound ligand is shown as a yellow ball-and-stick model in the closed-holo structure. The secondary substructures α-helices and β-sheets of the protein are displayed in red and blue colors, respectively. The protein is divided into the substrate-binding domain and the 1,4-dihydronicotinamide adenine dinucleotide (NAD)-binding domain, which is marked in the structure; (C) The dimeric L-MtAlaDH molecule. One monomer is shown in the solid ribbon model, and the other is shown in gray; (D) The hexameric L-MtAlaDH molecule. One monomer is shown in the solid ribbon model, and others are shown in gray.

Comparison of the Predicted Mean-Square Fluctuations with the Experimental Temperature Factors (B-Factors)
According to Equation (7), the B-factors of the residues in the open-apo and closed-holo structures are calculated. In Equation (7), the only adjustable parameter is γ, which is determined by normalizing the computational B-factors to the experimental data. The resulting γ  The bound ligand is shown as a yellow ball-and-stick model in the closed-holo structure. The secondary substructures α-helices and β-sheets of the protein are displayed in red and blue colors, respectively. The protein is divided into the substrate-binding domain and the 1,4-dihydronicotinamide adenine dinucleotide (NAD)-binding domain, which is marked in the structure; (C) The dimeric L-MtAlaDH molecule. One monomer is shown in the solid ribbon model, and the other is shown in gray; (D) The hexameric L-MtAlaDH molecule. One monomer is shown in the solid ribbon model, and others are shown in gray.

Comparison of the Predicted Mean-Square Fluctuations with the Experimental Temperature Factors (B-Factors)
According to Equation (7), the B-factors of the residues in the open-apo and closed-holo structures are calculated. In Equation (7), the only adjustable parameter is γ, which is determined by normalizing the computational B-factors to the experimental data. The resulting k B T{γ values used for the open-apo and closed-holo structures in the GNM are 3.68 and 6.90 Å 2 , respectively. The resulting k B T{γ values used for the open-apo and closed-holo structures in the ANM are 0.45 and 11.64 Å 2 , respectively. Figure 2 shows the comparison between the calculated and the experimental B-factors for the two studied systems. It is found that for GNM, the correlation coefficient between the calculated and experimental Temperature Factors (B-Factors) is 0.50 for the open-apo structure and 0.63 for the closed-holo structure, respectively. For ANM, the correlation coefficients are 0.49 and 0.51 for the open-apo and closed-holo structures, respectively. The correlation coefficient for the two proteins studied in this work is similar to that for other protein [12,24], which provides the validity for the study of MtAlaDH by using the GNM and the ANM.
Int. J. Mol. Sci. 2015, 16, page-page the two proteins studied in this work is similar to that for other protein [12,24], which provides the validity for the study of MtAlaDH by using the GNM and the ANM.

The Slow Motion Modes
In the GNM, the slow motion modes represent the large-scale collective motions, which are usually relevant to the functions of the protein [25]. Figure 3A shows the first slowest mode calculated by GNM for the open-apo and closed-holo structures of L-MtAlaDH. From Figure 3A, it is found that the first modes for the two studied proteins have a common hinge axes located in Met133 and Met301, with their fluctuation values equal to zero. The fluctuation profiles in Figure  3A also exhibit that the whole structure of the protein is divided into two distinct domains with large fluctuations. Comparing to the slowest mode of the open-apo structure, it can be found that the fluctuation of some residues around the residue His96 is reduced remarkably in the closed-holo structure. This is consistent with the simulation results that the root mean-square fluctuation (RMSF) of the residues around the residue His96 becomes small in the closed state of the protein

The Slow Motion Modes
In the GNM, the slow motion modes represent the large-scale collective motions, which are usually relevant to the functions of the protein [25]. Figure 3A shows the first slowest mode calculated by GNM for the open-apo and closed-holo structures of L-MtAlaDH. From Figure 3A, it is found that the first modes for the two studied proteins have a common hinge axes located in Met133 and Met301, with their fluctuation values equal to zero. The fluctuation profiles in Figure 3A also exhibit that the whole structure of the protein is divided into two distinct domains with large fluctuations.
Comparing to the slowest mode of the open-apo structure, it can be found that the fluctuation of some residues around the residue His96 is reduced remarkably in the closed-holo structure. This is consistent with the simulation results that the root mean-square fluctuation (RMSF) of the residues around the residue His96 becomes small in the closed state of the protein [10]. It is also found that the fluctuation value of the residues around Asp270 and Pro242 is also approximate to zero, although those residues do not belong to the hinge region (Residues 126-133, 304-320) [10]. These two residues are located at the ligand-binding cleft. This result indicates that the ligand-binding cleft become stable in the closed structure of the protein. Besides the two residues, the residues Met133, Thr178 and Met301 of the NAD-binding domain in the closed-holo structure are located in the minima of the fluctuation curve in Figure 3A. These residues are the ligand-binding sites and are locate in the jaws of the ligand-binding pocket [8]. The crystal structure of the closed-holo state of the protein exhibits that the nicotinamide ring of the ligand is bound in the residues Met133, Ser134, Ala137, Ile267, Asp270 and Met301, and the pyrophosphate moiety forms contacts with the main-chain atoms of residues 178 and 179 [9]. The decrease of the fluctuation of these residues indicates that they are more stable in the closed-holo structure than the open-apo structure. Figure 3B,C show the second and third slowest modes calculated by the GNM for the open-apo and closed-holo structures of L-MtAlaDH, respectively. The overlap between the open-apo and closed-holo structure of the second and third slowest modes is better than the first mode. In Figure 3B [10]. It is also found that the fluctuation value of the residues around Asp270 and Pro242 is also approximate to zero, although those residues do not belong to the hinge region (Residues 126-133, 304-320) [10]. These two residues are located at the ligand-binding cleft. This result indicates that the ligand-binding cleft become stable in the closed structure of the protein. Besides the two residues, the residues Met133, Thr178 and Met301 of the NAD-binding domain in the closed-holo structure are located in the minima of the fluctuation curve in Figure 3A. These residues are the ligand-binding sites and are locate in the jaws of the ligand-binding pocket [8]. The crystal structure of the closed-holo state of the protein exhibits that the nicotinamide ring of the ligand is bound in the residues Met133, Ser134, Ala137, Ile267, Asp270 and Met301, and the pyrophosphate moiety forms contacts with the main-chain atoms of residues 178 and 179 [9]. The decrease of the fluctuation of these residues indicates that they are more stable in the closed-holo structure than the open-apo structure. Figure 3B,C show the second and third slowest modes calculated by the GNM for the open-apo and closed-holo structures of L-MtAlaDH, respectively. The overlap between the open-apo and closed-holo structure of the second and third slowest modes is better than the first mode. In Figure 3B  It should be noted that besides the structures studied in this work, there are seven other structures of this protein that have been deposited in PDB. The first three slowest modes of these structures of this protein have also been calculated by using the GNM (data not shown). The mode shapes of these structures are found to be similar to those of the open-apo or closed-holo structures It should be noted that besides the structures studied in this work, there are seven other structures of this protein that have been deposited in PDB. The first three slowest modes of these structures of this protein have also been calculated by using the GNM (data not shown). The mode shapes of these structures are found to be similar to those of the open-apo or closed-holo structures studied in this work, which indicates that the functional motions are largely determined by the structural topology of the protein instead of the detailed local structures.
GNM can only describe the amplitude of residue fluctuations, but cannot give the direction of the fluctuations. In order to reveal the direction of the motions, ANM is used for the open-apo and closed-holo structures. The first three slowest normal modes calculated by ANM are shown in Figure 4A,B for the open-apo and the closed-holo structures of L-MtAlaDH, respectively. In Figure 4, the amplitude and direction of the motions are represented by the length and direction of the blue arrows, respectively. The red arrows and black dashed lines denote the motion direction and the motion axis, respectively. The first slowest mode ( Figure 4 (Mode 1)) of the two structures exhibits a twist motion of the two domains in the proteins. In this motion, the NAD-binding domain (NBD) and the substrate-binding domain (SBD) rotate reversely around their respective axes. The rotation axes are displayed in Figure 4 (Mode 1). This motion mode might be important for the two domains in the proteins to adjust their relative position to bind or dissociate the ligand. The long blue arrows represent the amplitude of the motions. This is consistent with the first slowest mode of GNM. The second slowest mode describes a hinge-bending motion of the NBD and SBD for the two studied proteins, as shown in Figure 4 (Mode 2). The α-helixes connecting the NBD and SBD serve as the hinge of the bending motion. The motion results in the opening and closing of the ligand-binding pocket. The studies by MD simulations also showed that the SBD rotates about 14.2˝around the hinge axis toward the NBD of the protein [10]. This is consistent with the second slowest motion mode of GNM, where the fluctuation of the C-terminus is larger than other residues. The third slowest mode ( Figure 4 (Mode 3)) represents a reverse rotation of the NBD and SBD around their respective axes. Additionally, the rotation axes of these two domains are marked by black dashed lines in this figure.
Int. J. Mol. Sci. 2015, 16, page-page studied in this work, which indicates that the functional motions are largely determined by the structural topology of the protein instead of the detailed local structures.
GNM can only describe the amplitude of residue fluctuations, but cannot give the direction of the fluctuations. In order to reveal the direction of the motions, ANM is used for the open-apo and closed-holo structures. The first three slowest normal modes calculated by ANM are shown in Figure 4A,B for the open-apo and the closed-holo structures of L-MtAlaDH, respectively. In Figure 4, the amplitude and direction of the motions are represented by the length and direction of the blue arrows, respectively. The red arrows and black dashed lines denote the motion direction and the motion axis, respectively. The first slowest mode ( Figure 4 (Mode 1)) of the two structures exhibits a twist motion of the two domains in the proteins. In this motion, the NAD-binding domain (NBD) and the substrate-binding domain (SBD) rotate reversely around their respective axes. The rotation axes are displayed in Figure 4 (Mode 1). This motion mode might be important for the two domains in the proteins to adjust their relative position to bind or dissociate the ligand. The long blue arrows represent the amplitude of the motions. This is consistent with the first slowest mode of GNM. The second slowest mode describes a hinge-bending motion of the NBD and SBD for the two studied proteins, as shown in Figure 4 (Mode 2). The α-helixes connecting the NBD and SBD serve as the hinge of the bending motion. The motion results in the opening and closing of the ligand-binding pocket. The studies by MD simulations also showed that the SBD rotates about 14.2° around the hinge axis toward the NBD of the protein [10]. This is consistent with the second slowest motion mode of GNM, where the fluctuation of the C-terminus is larger than other residues. The third slowest mode (Figure 4 (Mode 3)) represents a reverse rotation of the NBD and SBD around their respective axes. Additionally, the rotation axes of these two domains are marked by black dashed lines in this figure.  [26]. In this figure, the protein structures are displayed in gray tubes. The amplitude and direction of the motions are represented by the length and direction of the blue arrows, respectively. The red arrows and black dashed lines denote the motion direction and the motion axis, respectively.
In order to investigate whether the normal modes can reveal the conformational change observed by the experiment, the overlap between each mode and the experimental conformational change was calculated using Equation (17). Figure 5 shows the overlap for the first 50 modes. It is found that Mode 2 mostly contributes to the conformational change, where the value of the overlap is 0.658. In order to investigate whether the normal modes can reveal the conformational change observed by the experiment, the overlap between each mode and the experimental conformational change was calculated using Equation (17). Figure 5 shows the overlap for the first 50 modes. It is found that Mode 2 mostly contributes to the conformational change, where the value of the overlap is 0.658. In order to analyze the coupling between the residues motions, the cross-correlation between residue fluctuations is calculated using Equation (10). Because the low-frequency normal modes are usually relevant to protein functions, the first 30 slowest modes are considered in the calculation. The calculation results are displayed in Figure 6. The value of the cross-correlation ranges from −1 to 1, where the positive value represents the residues moving in the same direction and the negative value means the residues moving in the opposite direction. In Figure 6, the regions in orange-red color represent positive correlation and the blue regions denote negative correlation. It is found that the cross-correlation map, both for the open-apo and closed-holo structures, is divided into five orange regions and four blue regions separated by the residues Met133 and Met301, which means that these two residues serve as the hinges for the correlated motions. This is consistent with the result obtained by the analysis of the slowest normal mode, in which these two residues act as the motion hinges with almost zero fluctuations. In Figure 6, the blue regions correspond to the negative correlation between the motions of NBD and SBD, illustrating the open-closed conformational transition of the protein. The orange-red region in the center of the maps (Residues 133-301) represents the residues of the NBD moving as a whole, whereas the other four orange-red regions indicate the residues of the SBD also moving in a coupled way. This result indicates that the two domains keep their rigid structure during the conformational transitions. Comparing the closed-holo structure to the open-apo structure of the protein, more contacts are formed for the residues in the jaws of the ligand-binding pocket. These residue contacts may reduce the correlation between domain motions. Therefore, the blue and orange-red colors are weakened in the correlation map for the closed-holo structure compared to that for the open-apo structure, as shown in Figure 6. From Figure 6B, it is found that the fluctuation correlations among the residues Arg15, Loops L1 (Residues 94-99), L2 (Residues 238-247) and L3 (Residues 267-270) are significantly improved in the closed-holo structure because of their interactions with the ligand. In Figure 6, there are several regions with high correlation values, marked by the rectangles. These regions correspond to the highly coupled movement of residues 134-241 and 243-300 of NBD and residues 16-74 of SBD. In order to analyze the coupling between the residues motions, the cross-correlation between residue fluctuations is calculated using Equation (10). Because the low-frequency normal modes are usually relevant to protein functions, the first 30 slowest modes are considered in the calculation. The calculation results are displayed in Figure 6. The value of the cross-correlation ranges from´1 to 1, where the positive value represents the residues moving in the same direction and the negative value means the residues moving in the opposite direction. In Figure 6, the regions in orange-red color represent positive correlation and the blue regions denote negative correlation. It is found that the cross-correlation map, both for the open-apo and closed-holo structures, is divided into five orange regions and four blue regions separated by the residues Met133 and Met301, which means that these two residues serve as the hinges for the correlated motions. This is consistent with the result obtained by the analysis of the slowest normal mode, in which these two residues act as the motion hinges with almost zero fluctuations. In Figure 6, the blue regions correspond to the negative correlation between the motions of NBD and SBD, illustrating the open-closed conformational transition of the protein. The orange-red region in the center of the maps (Residues 133-301) represents the residues of the NBD moving as a whole, whereas the other four orange-red regions indicate the residues of the SBD also moving in a coupled way. This result indicates that the two domains keep their rigid structure during the conformational transitions. Comparing the closed-holo structure to the open-apo structure of the protein, more contacts are formed for the residues in the jaws of the ligand-binding pocket. These residue contacts may reduce the correlation between domain motions. Therefore, the blue and orange-red colors are weakened in the correlation map for the closed-holo structure compared to that for the open-apo structure, as shown in Figure 6. From Figure 6B, it is found that the fluctuation correlations among the residues Arg15, Loops L1 (Residues 94-99), L2 (Residues 238-247) and L3 (Residues 267-270) are significantly improved in the closed-holo structure because of their interactions with the ligand. In Figure 6, there are several regions with high correlation values, marked by the rectangles. These regions correspond to the highly coupled movement of residues 134-241 and 243-300 of NBD and residues 16-74 of SBD.

The Fast Motion Modes
In GNM, the fast motion modes reflect the geometric irregularity in the local structure, and the fluctuations in the fast motion modes are accompanied with the remarkable increases of the entropy. Then, the residues activated in the fast motion modes are regarded as the key residues that are important for the stability of the protein structure [21,25]. Figure 7 shows the eight fastest modes for the open-apo and closed-holo structures, respectively. It is found that the mode shapes for the open-apo and closed-holo structures are extremely similar, and the residues at the peaks of the fluctuation are almost the same, implying that the same set of residues is responsible for the stability of the open-apo and closed-holo structures of the protein. This result implies that the local structures are similar for these two proteins, and the domains of the protein keep their rigid structures during the open-closed motions. Figure 7 shows that the hot residues responsible for the stability of the proteins include Leu149, Val173, Asp198, Val263, Ile267, Cys274 and Val298 in NBD and Gln36, Val64, Thr112 and Thr345 in SBD. From the crystal structure of the proteins, it is found that these key residues are tightly packed, and most of them lie in the ligand-binding pocket, among which the residues Asp198 and Ile267 form polar interactions with the ligand NAD [8]. Our results imply that the residues crucial for protein stability also play important roles for the binding of the substrate to the protein.

The Fast Motion Modes
In GNM, the fast motion modes reflect the geometric irregularity in the local structure, and the fluctuations in the fast motion modes are accompanied with the remarkable increases of the entropy. Then, the residues activated in the fast motion modes are regarded as the key residues that are important for the stability of the protein structure [21,25]. Figure 7 shows the eight fastest modes for the open-apo and closed-holo structures, respectively. It is found that the mode shapes for the open-apo and closed-holo structures are extremely similar, and the residues at the peaks of the fluctuation are almost the same, implying that the same set of residues is responsible for the stability of the open-apo and closed-holo structures of the protein. This result implies that the local structures are similar for these two proteins, and the domains of the protein keep their rigid structures during the open-closed motions. Figure 7 shows that the hot residues responsible for the stability of the proteins include Leu149, Val173, Asp198, Val263, Ile267, Cys274 and Val298 in NBD and Gln36, Val64, Thr112 and Thr345 in SBD. From the crystal structure of the proteins, it is found that these key residues are tightly packed, and most of them lie in the ligand-binding pocket, among which the residues Asp198 and Ile267 form polar interactions with the ligand NAD [8]. Our results imply that the residues crucial for protein stability also play important roles for the binding of the substrate to the protein.

The Functionally-Important Residues Identified with the Thermodynamic Method
In our work, the functionally-key residues responsible for the ligand binding were identified by using a thermodynamic method proposed by Su et al. [22,23,27]. The closed-holo structure of the protein is adopted as the complex structure, and the corresponding GNM was constructed. Then, the ligand NAD was removed from the closed-holo structure, and the corresponding GNM was also constructed. Based on these two constructed models, each residue of the protein was perturbed, and the ΔΔG value in response to the perturbation was calculated according to the method proposed in the Methods Section. In our method, the free energy change ∆∆G for each residue is the sum of the free energy changes for all of the spring perturbations involved in this residue. Residues with relatively high ΔΔG values were regarded as the key residues involved in the ligand binding of the protein. The calculated result is shown in Figure 8. It was found that nine clusters of residues exhibit high ΔΔG values, which are centered at the residues His96, Glu118, Ala126, Met133, Thr178, Lys203, Pro242, Asp270 and Ala338, respectively.
From Figure 8, it is found that almost all of these key residues are located at the ligand-binding site in the NBD. Tripathi et al. have proposed that NAD interacts extensively with enzyme through the interactions with Met301, Asp270, Ile267, Thr178, Lys203, Asp198, Met133, His96, Lys75 and Arg15 [8]. Most of our predicted residues belong to the residues directly interacting with NAD. The analysis of the slowest normal mode discussed in the previous section has identified Met133 to be the motion hinge. This residue plays an important role for the opening and closing of the

The Functionally-Important Residues Identified with the Thermodynamic Method
In our work, the functionally-key residues responsible for the ligand binding were identified by using a thermodynamic method proposed by Su et al. [22,23,27]. The closed-holo structure of the protein is adopted as the complex structure, and the corresponding GNM was constructed. Then, the ligand NAD was removed from the closed-holo structure, and the corresponding GNM was also constructed. Based on these two constructed models, each residue of the protein was perturbed, and the ∆∆G value in response to the perturbation was calculated according to the method proposed in the Methods Section. In our method, the free energy change ∆∆G for each residue is the sum of the free energy changes for all of the spring perturbations involved in this residue. Residues with relatively high ∆∆G values were regarded as the key residues involved in the ligand binding of the protein. The calculated result is shown in Figure 8. It was found that nine clusters of residues exhibit high ∆∆G values, which are centered at the residues His96, Glu118, Ala126, Met133, Thr178, Lys203, Pro242, Asp270 and Ala338, respectively.
From Figure 8, it is found that almost all of these key residues are located at the ligand-binding site in the NBD. Tripathi et al. have proposed that NAD interacts extensively with enzyme through the interactions with Met301, Asp270, Ile267, Thr178, Lys203, Asp198, Met133, His96, Lys75 and Arg15 [8]. Most of our predicted residues belong to the residues directly interacting with NAD. The analysis of the slowest normal mode discussed in the previous section has identified Met133 to be the motion hinge. This residue plays an important role for the opening and closing of the ligand-binding pocket, and thus, it is crucial for the binding of the ligand to the protein. The residues His96 and Asp270 were identified as the conserved active site residues by Agren's group [9]. The residue Ala126 is also located in the hinge region between the two domains of the protein. In the slowest motion mode of the closed structure of the protein, the residues Thr178, Pro242 exhibit small fluctuations, which imply that these residues play important roles in the binding of the ligand, stabilizing the closed-holo structure of L-MtAlaDH. The functionally-key residues predicted by the thermodynamic method are consistent with the experimental observations and the results obtained by the analysis of the slow modes.
Int. J. Mol. Sci. 2015, 16, page-page 10 ligand-binding pocket, and thus, it is crucial for the binding of the ligand to the protein. The residues His96 and Asp270 were identified as the conserved active site residues by Agren's group [9]. The residue Ala126 is also located in the hinge region between the two domains of the protein.
In the slowest motion mode of the closed structure of the protein, the residues Thr178, Pro242 exhibit small fluctuations, which imply that these residues play important roles in the binding of the ligand, stabilizing the closed-holo structure of L-MtAlaDH. The functionally-key residues predicted by the thermodynamic method are consistent with the experimental observations and the results obtained by the analysis of the slow modes.

The Gaussian Network Model and the Anisotropy Elastic Network Model
In the Gaussian network model (GNM) [11,28], the tertiary structure of a protein is simplified as an elastic network. Each residue of the protein is represented by a point that is located at its Cα atom, and all of the pairwise residues within a cutoff distance (8 Å is adopted here) are connected by a harmonic springs with a uniform force constant. Then, the internal energy of the protein can be expressed as [15]: where γ is the spring constant; R Δ represents fluctuations of the residues; E is the identity matrix; ⊗ is the direct product of the matrix; T represents the transpose of the matrix; Γ is the Kirchhoff matrix with N N × elements that can be expressed as [11,15]: where ij R is the separation between the i-th and j-th C α atoms and c r is the cutoff distance.
The cross-correlation between the fluctuations of the i-th and j-th residues is written as:

The Gaussian Network Model and the Anisotropy Elastic Network Model
In the Gaussian network model (GNM) [11,28], the tertiary structure of a protein is simplified as an elastic network. Each residue of the protein is represented by a point that is located at its C α atom, and all of the pairwise residues within a cutoff distance (8 Å is adopted here) are connected by a harmonic springs with a uniform force constant. Then, the internal energy of the protein can be expressed as [15]: where γ is the spring constant; ∆R represents fluctuations of the residues; E is the identity matrix; b is the direct product of the matrix; T represents the transpose of the matrix; Γ is the Kirchhoff matrix with NˆN elements that can be expressed as [11,15]: where R ij is the separation between the i-th and j-th C α atoms and r c is the cutoff distance. The cross-correlation between the fluctuations of the i-th and j-th residues is written as: When i = j, the mean-square fluctuation of the i-th residue can be written as: where k B is the Boltzmann constant; T is absolute temperature.
The inverse of the Kirchhoff matrix can be calculated by diagonal decomposition as: where U is an orthogonal matrix whose columns represent the eigenvectors of Γ; Λ is a diagonal matrix, and the diagonal elements are the eigenvalues of Γ.
According to the Debye-Waller theory, the B-factor of the i-th residue can be calculated by: The mean-square fluctuation of the i-th residue in h the k-th mode is written by: The cross-correlation between the fluctuations of different residues is written by: In the GNM, the cross-correlation is normalized as: The elements of H are a submatrix with a size of 3ˆ3. The ij-th submatrix h ij is: When i " j, where the elements B 2 V Bx i By j , B 2 V Bx i By i of the submatrix h ij can be analytically expressed as: where x, y and z represent the coordinates of each atom.

Correlation Coefficient and Overlap
In our work, the correlation coefficient between the calculated and the experimental B-factors is written as: where x i and y i are the calculated and experimental values for the B-factor of the i-th C α atom, x and y are the mean values of x i and y i and N is the total number of C α atoms in the protein.
Because the open-apo and closed-holo structures of the protein have been determined by the X-ray crystallographic experiment, the conformational change can be calculated. Then, we investigated whether the slow motion modes obtained by ANM can account for the experimental conformational change, which can be evaluated by the overlap between the normal modes and the conformational change. The value of overlap can be calculated by: (   17) where Ñ η i is the i-th normal mode, ∆ Ñ R exp is the experimental conformational change andˇˇˇˇ∆ Ñ R expˇi s the norm of ∆ Ñ R exp .

Thermodynamic Cycle Method
Base on the elastic network model (ENM), a thermodynamic method was developed by our group to reveal the functionally-important residues involved in the binding of ligand [22,23,27]. According to the theory of statistical physics, the vibrational entropy of ENM can be expressed as [25]: where xHy is the average value of the Hamiltonian; F "´k B TInZ is the Helmholtz free energy; T is the temperature; N is the number of residues in the protein; k B is the Boltzmann constant; Z is the configurational integral part of the vibrational partition function that can be calculated by Z " ş expp´H{k B Tqd t∆Ru; t∆Ru represents the 3N-dimensional column vector of fluctuations. When a perturbation of the spring between the i-th and j-th residues is introduced, the vibrational entropy of the protein will be changed as: where ∆γ ij is the change of the force constant that accounts for the perturbation of the spring. According to Equation (18), the derivation of the entropy with respect to γ ij can be expressed as [22,23]: In this method, all of the residues in the protein were perturbed one by one, and then, the change of the free energy for the protein-ligand binding was calculated.
In Figure 9, ∆G 1 denotes the ligand-binding free energy without the introduction of the perturbation, and ∆G 2 denotes the value of the binding free energy with a residue perturbation. In this work, the binding ligand is 1,4-dihydronicotinamide adenine dinucleotide (NAD). The receptor* and complex* represent the systems where a residue perturbation is introduced. We want to calculate the change of the free energy resulting from the spring perturbation, i.e., When a perturbation of the spring between the i-th and j-th residues is introduced, the vibrational entropy of the protein will be changed as: where Δγij is the change of the force constant that accounts for the perturbation of the spring. According to Equation (18), the derivation of the entropy with respect to γij can be expressed as [22,23]: In this method, all of the residues in the protein were perturbed one by one, and then, the change of the free energy for the protein-ligand binding was calculated.
In Figure 9, ΔG1 denotes the ligand-binding free energy without the introduction of the perturbation, and ΔG2 denotes the value of the binding free energy with a residue perturbation. In this work, the binding ligand is 1,4-dihydronicotinamide adenine dinucleotide (NAD). The receptor* and complex* represent the systems where a residue perturbation is introduced. We want to calculate the change of the free energy resulting from the spring perturbation, i.e., 2 1 Figure 9. The thermodynamic cycle, ΔG1 denotes the ligand-binding free energy without the introduction of the perturbation, and ΔG2 denotes the value of the binding free energy with a residue perturbation. The receptor* and complex* represent the systems where a residue perturbation is introduced. Two nonphysical processes were constructed, i.e., Receptor → Receptor* and Complex → Complex*. The free energy changes associated with these two constructed processes are denoted as ΔG3 and ΔG4.
The residues with relatively large ΔΔG values are identified as the functionally-important residues, which are thermodynamically coupled with the binding of ligand to the protein. In order to calculate the free energy changes more easily, a thermodynamic cycle was proposed. Two nonphysical processes were constructed, i.e., Receptor → Receptor* and Complex → Complex*. The free energy changes associated with these two constructed processes are denoted as ΔG3 and ΔG4. Then, Equation (21) can be expressed as: When a perturbation was introduced to a residue of the protein, the force constant of the corresponding spring is reduced, and the perturbation of γij results in the same change of the potential energy of Receptor and Complex, i.e., ΔU3 = ΔU4. Then: In our method, when a perturbation was introduced to a residue, all of the springs connecting to the residue should be perturbed in the same way. In this way, the ΔΔG value resulting from the perturbation of a residue is calculated by the sum of the ΔΔG values of all of the perturbed springs involved in this residue. It should be pointed out that only the first 30 slowest normal modes of the Figure 9. The thermodynamic cycle, ∆G 1 denotes the ligand-binding free energy without the introduction of the perturbation, and ∆G 2 denotes the value of the binding free energy with a residue perturbation. The receptor* and complex* represent the systems where a residue perturbation is introduced. Two nonphysical processes were constructed, i.e., Receptor Ñ Receptor* and Complex Ñ Complex*. The free energy changes associated with these two constructed processes are denoted as ∆G 3 and ∆G 4 .
The residues with relatively large ∆∆G values are identified as the functionally-important residues, which are thermodynamically coupled with the binding of ligand to the protein. In order to calculate the free energy changes more easily, a thermodynamic cycle was proposed. Two nonphysical processes were constructed, i.e., Receptor Ñ Receptor* and Complex Ñ Complex*. The free energy changes associated with these two constructed processes are denoted as ∆G 3 and ∆G 4 . Then, Equation (21) can be expressed as: When a perturbation was introduced to a residue of the protein, the force constant of the corresponding spring is reduced, and the perturbation of γ ij results in the same change of the potential energy of Receptor and Complex, i.e., ∆U 3 = ∆U 4 . Then: ∆∆G " p∆U 4´T ∆S 4 q´p∆U 3´T ∆S 3 q "´Tp∆S 4´∆ S 3 q In our method, when a perturbation was introduced to a residue, all of the springs connecting to the residue should be perturbed in the same way. In this way, the ∆∆G value resulting from the perturbation of a residue is calculated by the sum of the ∆∆G values of all of the perturbed springs involved in this residue. It should be pointed out that only the first 30 slowest normal modes of the closed-holo structure were considered in our calculation, and for the protein-ligand complex, all of the modes were projected onto these slowest modes.

Conclusions
The open-closed conformational changes of L-MtAlaDH induced by the coenzyme are important for the catalytic function of the protein. In this work, the domain motions and the functionally-key residues encoded in the structure of the protein were investigated by using the coarse-grained Gaussian network model (GNM) and the anisotropy network model (ANM). The calculated B-factors are consistent with the experimental data from the X-ray crystallographic experiments, which provide the validity for the application of the GNM to the studied system L-MtAlaDH. The analysis of the slowest normal mode for the open-apo and the closed-holo structures shows that the two domains in the proteins undergo large-scale motions, where the motion axes are centered at Met133 and Met301. The ANM calculation results indicate that both for the open-apo and the closed-holo structures, the first three slow modes respectively represent the twist, hinge-bending and rotation motions of the proteins. These motions are responsible for the conformational transition of the proteins and contribute to the opening and closing of the ligand binding pocket. The key residues responsible for the stability of the proteins were successfully identified by analyzing the high-frequency normal modes. It is found that the key residues are tightly packed, and most of them lie in the ligand-binding pocket. The results of cross-correlation analysis implied that the residues in two domains of the protein move in a highly coupled way. In addition, the key residues responsible for the ligand binding were identified by using the thermodynamic cycle method. The predicted functionally-key residues are consistent with the experimental observations. Our method is a simple and coarse-grained model, which considered the connections between C α atoms only. Our methods should be combined with other information to determine the functionally-important residues more exactly and should be combined with more novel methods to calculate the interaction of the side chain. We will construct a more detailed model, in which the side chains of residues are explicitly considered, in the next step to improve our method.