Molecular Dynamics Simulation of Cholera Toxin A-1 Polypeptide

3.1 Total energies of the system

The energetically favorable system represents a system that depicted the dynamic and thermodynamic behavior of a normal protein as persisted in nature. Figure 2 showed the variation in total energy with respect to simulation time for CTA1 at 283, 310 and 323 K respectively. A decrease in the total energy occurs at the beginning of the simulation and is continued up to 100 picoseconds (ps) time span, then the total energy of cholera toxin A-1 polypeptide is almost constant for the rest of the simulation time. The mean total energy was found to be – 4.47 × 10⁵ kJ mol^-1 at 283 K, and at physiological temperature, the mean value of total energy was – 4.19 × 10⁵ kJ mol^-1. In simulation performed at physiological temperature an energy increase from – 4.3 to – 4.2 × 10⁵ kJ mol^-1 occurred from the start, which then became constant, however a small fluctuation was observed throughout the simulation time of 1000 ps. At 323 K, a shift in total energy was observed at the beginning and the mean total energy of – 4.075× 10⁵ kJ mol^-1 was recorded, and remained almost constant with minor variations as shown in Figure 2. These results also revealed that with a rise in simulation temperature the total energy of the protein in the box is found to increase, which causes the simulated system to become unstable and provided proof that with an increase in temperature the cholera toxin A-1 polypeptide became thermally unstable. For every 10-digit rise in temperature the total energy rise was about 0.1× 10⁵ kJ mol^-1. It has been confirmed from the above mentioned figure that the total energy values remained negative and CTA1 was thermally stable in the temperature range of 283 to 323 K.

Figure 1

Cholera toxin A1 in combination with adenosine ribosyl factor 6 from protein data bank file, entry number 2A5F [11].

Figure 2

Total energy fluctuations versus simulation time for CTA1 at 283 K (black), 310 K (red) and 323 K (green).

3.2 Root mean square deviation

One of the validity tests of standard MD studies is the determination of root mean square deviation (RMSD) [39].

We also performed this obligatory analysis to determine the stability of the molecular dynamics simulation trajectories. The time variable 185 α-Carbon (Cα) RMSD of CTA1 was calculated at the respective temperatures, taking the energetically unminimized crystal structures as a reference. Figure 3 shows a sudden jump of RMSD values at the very beginning and then making a plateau, thus convergence was observed. However, the results of Cα RMSD of CTA1 at three different temperatures have very small convergence, but the fluctuations were very prominent as shown in Figure 3. The MD simulation of CTA1 at 283 K showed Cα RMSD of 0.0267 nm at the beginning and then rose to 0.031 nm, subsequently after this a continuous fluctuation was present throughout the simulation time. The average Cα RMSD for this system was calculated to be 0.0308 nm at 283 K (Figure 3). Almost similar Cα RMSD results were obtained for MD simulation at 310 and 323 K with the same degree of fluctuations as in the case of 283 K. These high flexibilities in Cα RMSD of CTA1 might be due to the partially disordered secondary structure of cholera toxin A-1 polypeptide. These results were in agreement with the past analysis of MD simulation on various proteins [40].

Figure 3

Root mean square deviation of Cα atoms of CTA1 at 283 K (black), 310 K (red) and 323 K (green).

3.3 Root mean square fluctuation (RMSF)

The RMSF analysis is the most productive and fruitful character of MD simulation for enzymes and other biochemical compounds that are harmful to human life in one aspect or another. It is because of these observations that we determined the fluctuations of different residues of the compounds under study. In the literature, there are two opposite views of designing the inhibitory drugs for toxins, one is by targeting either the most rigid residues and secondly those that are highly flexible ones. But both the mobile and immobile residues and regions of an enzyme have specific functional rolesin different biological processes. For example molecular recognition and catalytic properties are due to the mobility possessed by definite and specific residues while the rigidity of residues produced different structures in beta-folds that are indispensable for substrate interaction and catalytic action [41]. The MD study performed by Grottesi and Sansom [42] on the venom of scorpion showed that the most immobilized residues, that have structural and functional role in the toxin, could be exploited for designing inhibitory drugs, while an opposite view has been given by Caceres et al. [43]. Based on MD investigations, they utilized the most flexible regions and loops as a substrate for designing the inhibitory drugs which act as ligands and could be readily accommodated inside those flexible pockets and regions. Similarly the different types of interactions between the enzyme residues and drugs are possible as the mobile residues have the ability to reduce the steric clashes and could be more receptive conformationally.

In order to know the flexible and rigid behavior of residues of the toxin in a dynamic state, residue wise Cα RMSF were computed at three different temperatures. In MD simulation of CTA1 at 283 K, isoleucine 83 showed the most rigid behavior with RMSF of 0.0129 nm, as shown in Figure 4. Similarly isoleucine 14, 122; leucine 39, valine 58, aspargine 140, tyrosine 126 and tryptophan 172 were also inflexible. On the other hand, glycine 184 was the most fluctuating having Cα-RMSF of 0.0276 nm and aspargine 1 has CαRMSF of 0.025 nm, the reason was that these residues are present at the C and N-termini having no interaction with other residues. Similarly, glycine 77, alanine 101, 157 and proline 116 are also highly flexible. The average Cα RMSF of CTA1 polypeptide at 283 K was 0.0175 nm. At the physiological temperature, the CTA1 residues fluctuate differently. Glutamine 121 was the most rigid throughout the simulation having Cα RMSF of 0.0124 nm, while valine 85 and 111, tyrosine 28, 40, 126 and 148; aspartic acid 57, isoleucine 178, aspargine 140 and histidine 138 are also rigid. The high fluctuating residue is glycine 184 at the carboxyl terminus having CαRMSF of 0.028 nm. The other residues are glycine 18, 49, 71, 77, 116 and alanine 101 are also highly flexible. The average Cα RMSF at 310 K was 0.0175 nm and at 323 K aspargine 41 was highly conserved with RMSF of 0.0125 nm. Similarly valine 58, leucine 3 and 99, glutamine 121, isoleucine 122, and histidine 179 are also rigid as shown in Figure 4. While the highly flexible residues were glycine 18, 49, 71, 116 and 131 but glycine 184 has the highest Cα RMSF of 0.0259 nm. The average RMSF at 323 K was found to be 0.0175 nm. Based on these RMSF results, we can suggest to use the molecular modeling and docking methods to explore the whole enzymatic portion of CTA1 for possible binding pockets. These molecular docking studies will be helpful in designing drugs that can bind to the active sites and other pockets in CTA1. Based on the available methods of molecular docking, both rigid and flexible residues need to be targeted to discover suitable small molecule ligands that acts as inhibitors of CTA1 enzymatic activity.

Figure 4

RMSF for each Cα atom of CTA1 at 283 K (black), 310 K (red) and 323 K (green).

3.4 Radius of gyration

In order to assess the folding and unfolding pathway of CTA1, we calculated the radius of gyration (Rg) at different temperatures. The result varies for CTA1 in the dynamic phase for 1000 ps. A continuous fluctuation in Rg values at 283 K was oberserved, which exhibited the folding, unfolding of protein as the time passed as shown in Figure 5. The mean Rg value was found to be 1.734 nm, with a minimum at 1.732 nm at 740 ps, while the maximum was 1.736 nm at 180 ps and at 550 ps. The MD simulation run at physiological temperature showed the mean Rg value of 1.734 nm. In the current study the highest value recorded at this temperature was 1.736 nm near 170 ps while the smallest one was 1.717 nm at 500 ps, similarly, the mean Rg value of CTA1 was observed to be 1.734 nm at 310 K (Figure 5). The radius of gyration results of CTA1 simulation at 323 K looked different from the first two-simulations of 283 and 310 K. One can observe a continuous fluctuation in the radius of gyration values, but not so prominent when compared to the earlier two MD simulations at 283 and 310 K (Figure 5). Both the minimum and maximum Rg values of 1.742 at 190 ps and 1.737 nm at 810 ps and 990 ps were observed respectively. Thus, at high temperatures the cholera toxin A1 polypeptide was expanded due to increased motion of different groups of the polypeptide, as a result it follows a stable pathway for its unfolding.

Figure 5

Radius of gyration of Cholera toxin at temperature of 283 K (black), 310 K (red) and 323 K (green).

3.5 Solvent accessible surface area (SASA)

The solvent accessible surface area of CTA1 polypeptide was considered in order to find the role of different hydrophobic and hydrophilic forces and residues exposed to the solvent at different temperatures. The hydrophobic residues present at the end of the polypeptide chain are proline and tryptophan but the CTA1 is also strongly hydrophilic which has been proved from water of crystallization present in the crystal structure [7, 11]. The cholera toxin simulation at 283, 310 and 323 K showed a mean total SASA of 111 nm², which was exposed to the solvent in the system. The hydrophilic SASA was greater than the hydrophobic SASA as shown in the supplementary figures 1 A-C. The mean values of hydrophilic and hydrophobic SASA were observed to be 58.5 nm² and 52.5 nm² respectively for the CTA1; thus, there was a difference of 6 nm². The SASA calculated for CTA1 was not only equal, but nearly remained constant with a fluctuation of about 2 nm² at all the simulation temperatures during the time span. Thus, this small fluctuation can be assigned to the folding and unfolding of toxin, which has been confirmed from the radius of gyration studies. It was further confirmed from SASA studies that CTA1 exhibits nearly 50% hydrophobic and 50% hydrophilic SASA. Due to this partial hydrophobic nature, it can easily move through the plasma membrane of the cell and the membrane of the endoplasmic reticulum while the partial hydrophilic nature provides it the opportunity to make interaction in the hydrophilic cytosol.

According to Ampapathi et al. [16] cholera toxin is hydrophobic in nature. They divided the CTA1 polypeptide into three regions as CTA1₁, CTA1₂, and CTA1₃, where the 1^st region is globular, 2^nd forms an extended bridge while a 3^rd is highly hydrophobic due to the hydrophobic amino acids present in these regions. Contrary to all this, Zhang et al. [2,3] and O’Neal et al. [11] obtained the crystals of cholera toxin with X-ray diffraction analysis which confirmed its hydrophilic nature because most of the crystals retained up to eighty water of crystallization molecules and depicted its hydrophilic nature. Our results are consistent with those of Zhang et al. [2,3] and O’Neal et al. [11] but are opposite to those of Ampapathi et al. [16]. From the above discussion, it is concluded that cholera toxin can spread through aerosols and it is a potent and easy to carry biological weapon. Therefore in order to prevent its misuse by rogue countries and terrorist groups as a biological weapon, inhibitory drugs and methods should be developed for its control.

3.6 Secondary structure analysis

The secondary structure of cholera toxin is of great interest because of its unusual entry mechanism of translocation inside the cell cytosol [44-46]. It utilizes the retrotranslocation mechanism and deceives the proteolytic enzymes of endoplasmic reticulum and also uses the process of unfolding and refolding mechanism [15,17,20]. The secondary structure of CTA1 was calculated through the timeline plugin module of VMD having structural identification (STRIDE) of protein algorithm. The STRIDE algorithm computes the secondary structure of the given protein as a function of time based on the changes in the hydrogen bond analysis, main chain dihedral angles and hydrogen bond distances of amino acid residues. Supplementary figure 2 A-C showed strips representing the secondary structure of CTA1 at respective temperatures as a function of time. The secondary structures contain the regular conformation units of α-helices, 3-10 helices, β-sheets and β-turns while the irregular conformations are loops or coils. In the CTA1 simulation at 283 K, the secondary structure from the N-terminal, aspartic acid one made a coil continuously for 1000 ps time span. Lysine 2 showed a transition of β-sheet and coil, but the β-sheet form was dominant. The amino acid residues from leucine 3 up to aspartic acid 7 form a continuous extended bridge of β-sheet for 1000 ps simulation time. The three residues in CTA1 sequence from serine 8 to proline 10 were in irregular conformation of the coil form. Proline 11 to serine 17 was in a magenta color strip, that represented the α-helix at certain positions, transition in this secondary structure to β-turn occurred which were denoted by green patches in supplementary figure 2-A. Serine 18, which was dominant in coil form, transferred several times to β-turn form during the course of simulation time. After this, glycine 19 and leucine 20 were in β-sheet. Methionine 21 to proline 22 is in coil structure, but transition has occurred to β-turn due to methionine. Then arginine 23 to arginine 31 was in continuous α-helix. But glycine 32 to methionine 35 again was in coil conformation. The transition of tan color to white color that was representative of the transtition of the isolated bridge form to the coil form was present throughout the simulation time span. Isoleucine 377 and aspargine 38 were in coil secondary structure continuously. For the third time α-helix was formed by residues from leucine 39 to arginine 44. Glycine 45 again forms a coil continuously. A continuous β-turn was formed by threonine 46 to glycine 49. A coil was formed by phenylalanine 50 and valine 51. The β-turn represented by green color was formed by arginine 52 to aspartic acid 54. The transition of β-turn into isolated bridge and vice versa, was formed by aspartic acid 55. Glycine 56 was a part of a coil. β-Sheet was then formed by tyrosine 57 and valine 58 continuously throughout the simulation. Serine 59, threonine 60 presented a transition of coil and β-sheets several times. Serine 61 formed the isolated bridge, but transition to β-sheet occurs during the simulation. Coil form was again observed due to isoleucine 62 and serine 63. The transition of isolated bridges and β-turn into one another was produced by threonine 73 to leucine 75, but leucine 64 and leucine 75 were in isolated bridge form of secondary structure. The β-turn was formed again by serine 76 to serine 79. Threonine 80 to isoleucine 86 was in the β-sheet form but threonine 80 transfers to isolated bridge form for certain time periods of simulation. Alanine 87 and threonine 88-formed coil but alanine 87 has a transition to β-sheet. The β-turn was shown by residues of alanine 89 to aspargine 91. Methionine 92 to aspargine 94 formed β-sheet while valine 95 to leucine 99 formed an α-helix. Glycine 100 to tyrosine 102 showed a transition of 3-10 helix shown by blue color band and its transition to β-turn several times during the simulation as reflected in supplementray figure 2-A. Serine 103 to histidine 5 formed β-turn. Proline 106 to aspartic acid 108 also showed transition from 3-10 helix to β-turn. Glutamine 109 produced a transition of coil and β-turn during the simulation. Aspartic acid 110 had transition of β-sheet and coil form of secondary structure several times during the simulation. The β-sheet was again represented by valine 111 to leucine 114. The two glycine residues of 115 and 116 formed β-turn. Then again isoleucine 117 and proline 118 formed β-sheet continuously. Another transition of secondary structure of β-turn into 3-10 helix and vice versa was shown by tyrosine 119 to glutamine 121. Isoleucine 122 to histidine 129 formed the β-sheet but transition to coil was shown by arginine 127 and to the isolated bridge by valine 128 and histidine 129 as observed in supplementary figure 2-A. An uninterrupted β-turn was produced by phenylalanine 130 and glycine 131. Valine 132 was switched to β-sheet and isolated bridge. The leucine 132 residue also showed similar behavior as above but coil form replaced the isolated bridge structure. Aspartic acid 134 to leucine 137 were in coil form but for a small time it shifted to β-sheet. The β-sheet was again formed by histidine 138 and arginine 139, and then β-turn was formed by aspargine 140 to tyrosine 143. Arginine 144 formed a coil but turned into β-turn. Another alteration of 3-10 helix into β-turn was produced by aspartic acid 145 to tyrosine 147. Next β-turn was represented by tyrosine 148 to leucine 151. Residues from aspartic acid 152 to proline 155 were in coil form. The 3-10 helix was regularly presented by alanine 156 to aspartic acid 158. Glycine 159 to glycine 164 were mostly in β-turn form, but changed to coil during the course of simulation by glycine 159 and alanine 163 and changed to 3-10 helix form and coil form by tyrosine 160 to leucine 162. A complete coil was shown by phenylalanine 165 and a full β-turn was presented by proline 166 to histidine 169. Another evolution of secondary structure of 3-10 helix to coil form was reflected in the supplementary figure 2-A by arginine 170. Alanine 171 to arginine 173 formed 3-10 helix but shifted to β-turn. Another β-turn was formed by glutamic acid 174 to proline 176. Tryptophan 177 to histidine 179 also reflects a transition of 3-10 helix into β-turn, however histidine 180 also showed coils form for several simulation times. The last conversion was produced in secondary structure of CTA1 by alanine 181 of coil form and β-turn while proline 182 to serine 185 were continuously in the coil form. From the results it was clear that the secondary structure graph produced through STRIDE algorithm showed several transitions of the secondary structure at the 283 K simulation. Further regular units were interrupted by irregular confirmation of secondary structure in this simulated system. The simulation performed at 310 K and 320 K for CTA1 (Supplementary figure 2-B, C) showed almost the same amount of secondary structure but the time and amount of transition were different. From these results we concluded that CTA1 has partially disordered secondary structure and can be placed in the intrinsically disordered (ID) protein class. The behavior of cholera toxin A1 polypeptide of regular units of helices and β-sheets interrupted by coils were in close agreement to the work of Ampapathi et al. [16] who observed the secondary structure through circular dichroism spectroscopy and 2D nuclear magnetic resonance spectroscopy and that of others [13-16].