The Fastest Simulation of Protein Folding Based on Torsion Angles


 Backgrounds:

Enormous number of possible conformations in the protein structure simulation have led molecular dynamics researchers to be frustrated until now. Some methods with defects ended their experiments into failure. This made them fail to determine the structure and function of folded protein in stable state with the lowest potential energy. This apparently exist in nature. The purpose of resolving a protein folding pathway that follows protein backbone residues torsional inertia was accomplished.
Results

A new method, torsion angle modeling, was adopted focused on the rotation of dihedral angles. The potential energy was calculated by rotating torsion angles of the peptide with 8 residues. It was found that when moving in the order of torsional inertia, 8 residues swivel in sequence. Six passes were repeated to find the lowest value.
Conclusion

The protein backbone torsion angle plays very important role in predicting protein structure. Actually it was thousand times faster or more than others to get the obvious pathway.

In nature protein folds faster and more naturally than this. Namely, in vivo folding is faster than in vitro. (Nicola et al., 1999;Kolb et al., 2000;Cabrita et al., 2010;O'Brien et al., 2011, Fedorov andBaldwin, 1999;Seckler et al., 1989) So Levinthal addressed paradoxically that there is the apparent pathway in the folding(J. T. P. DeBrunner and E. Munck, 1969). This is true. And Karplus suggested that the 'bias towards the native state' over much of the effective energy surface may govern the folding process(M. Karplus, 1997). We appreciate this is the pathway of folding that occurs in protein synthesis from mRNA transcript.
It is worth to look into cotranslational folding, because there are differences between simulated and experimental structures. We derived the structures from cotranslational and torsional algorithm.
Cotranslational folding simulation was managed before by using SAINT algorithm (Ellis et al., 2010). 3D molecular dynamics was used in this experiment. All movements in the Cartesian space were considered as possible in most molecular dynamics algorithms by most of the researchers. These employ strong covalent bond stretching and bond angle bending terms in their force elds. Though covalent bond stretching and bond angle bending are affected by the force eld, these model's description abandons the fact that only rotations about covalent single bond occurs.
We used, thus, backbone torsion angle as the only degree of freedom to describe the motion of protein folding. These models might hamper to interpret the mechanism of protein folding. And helpful manual modi cation is restricted from researchers. These di culties could be settled by torsion angle method. Also this is more realistic than other 3D models full of lattice approximations.

Methods
The polypeptide is free swivel chain and has levers which collide to water molecules. The collision in the protein folding in the cell follows the Brownian movement. It is much stronger. The calculation of the translational molecular kinetics is based on the translational enthalpy data. Thus, the kinetic energy of the water molecule is, And the velocity of water molecule gives, 1 calory is 4.1868 J (1 cal = 4.1868J), so Hence, velocity of water molecule is This velocity is fast enough to rotate backbone's residue with torques (Serway, R. A. and Jewett, Jr. J.W., 2003;F. H. Stillinger, 1975).
Torsional inertia in the swivels are different for each amino acid. In a backbone, only one residue rotates, the others do not. While the former has small inertia, the latter has larger one. The difference of torsional inertia comes from the type and number of atoms and the length of side chain, and the mass and length of the backbone on both sides of the residue. The middle residue has a small rotational inertia, because the rotational inertia is proportional to the square of the distance.
The following (Figure 1) represents only two residues among many in the backbone. Among these connected residues, only two residues with the least rotational inertia were shown. Between them, A has smaller rotational inertia than B. When water molecules collide, A rotates rst. B is a residue that stops while residue A rotates, and then returns as soon as A stops.
Residue A, has rotational momentum P a is, If B rotates, the rotation momentum P b of B is, When A does not rotates, the rotation momentum P b of B is If A rotates, the rotation momentum P b of B is If A stops rotation, P a becomes zero and P b is not reduced. So B can rotate more easily. If A rotates, P a becomes larger and P b is strongly reduced. Therefore, B cannot rotate easily.
In gure 2, the angle between the axel with A residue and B1 is 104.5°. Assuming that B rotates while A also rotates, residue A must handle this rotational movement of residue B. However, if it is 90°, because the radius of rotation(R = 1×R 0 ) is the longest, the rotational inertia of residue B is affected to residue A.
Since, residue A can hardly turn without the strongest force because rotational inertia(I = R 2 ×m) is the maximum.
To receive the half of the magnitude of the force mentioned above, it must be 135 °. The radius of B 2 is, And inertia I is, The bond angle between the two residues A and B is 104.5°. The inertia of B 2 at 104.5° is, This is almost the same as the value of 90° with the most di culty. In conclusion, residue B 2 cannot rotate with B 1 at the same time.
In short, all of the above are summarized as follows. If there are many residues of different sizes, they do not rotate at the same time. The residue with the smallest inertia rotates rst. As soon as it stops rotating, the residue with the next smallest inertia turns. When stops rotating, the next residue with more inertia than this rotates. In the case of multiple residues, the rotational inertia is ordered from the smallest.
The coordinates of short polypeptide with important points above was set. And it was compared with the structure from solution NMR. The cotranslational folding and torsional movement of atoms which are different from others was proved. The most representative model of each NMR assay was compared with the structure from the new folding algorithm. Typical structure alignment algorithms including TMalign(Y. Zhang and J. Skolnick, 2005) were not used in the comparison. These algorithms moves frame with insertions and deletions even in the case of the comparison of identical amino acid sequence. The result of the comparison was represented by logPr and RamRMSD(S. Jung, et al., 2011). And it was illustrated in the graph of torsion angle along the residue number( Figure 3). The change of the potential energy during the initial folding and following optimization is also displayed in the graph( Figure 4).

Data set
A single asymmetric chain structure was adopted in the PDB archive on the condition that the length of the chain was 8, the sequence identity was 90% or less, and that there were no heteroatoms. Structures which only contains protein without nucleic acids of DNA, RNA, or DNA-RNA hybrid was selected. Two structures of 1n9v and 1oeh was set. And 1oeh was abandoned because it was fragmented. Finally, 1n9v which is the angiotensin peptide was used in this work. 'DRVYIHPF' is the amino acid sequence of 1n9v.

Cotranslation Folding of Initial Structure
Cotranslational folding was performed with ProtTorter using torsion angles. Whenever a new amino acid was added, the potential energy was calculated, considering every conformation following the change of ϕ and ψ angles. As peptide bonds revolve around the backbone, we supposed they move 1 degree by 1 degree. So, 360 cases were observed. Local minima and the global minimum was gotten. These were to predict the initial structure of angiotensin.

Iterative Optimizations
This simulation became accurate by the peptide collision with water molecules, residues received larger torques, and planar amino acids to rotate faster. The length of the side chain was estimated as the maximum of bonds from the Cα atom. The order of calculation of the pair of bonds of the residue from the sum of the priority from torsional property.
Consequently, the 4 th peptide bond between Tyr(Y) and Ile(I) was rst calculated(1+5=6), and the bond between Tyr(Y) and Val(V)(1+8=9). The order of priority is 'D-(3)-R-(4)-V-(2)-Y-(1)-I-(6)-H-(7)-P-(5)-F'. The numbers in the round brackets designate the order of the linkage among residues. The dihedral angles anking peptide bond are φ and ψ angles. Between these two dihedral angles, we calculated the one close to the higher priority rst, and the other close to the lower priority later.
In this order, we iteratively optimized the initial structure from cotranslational folding. This was performed until the convergence of potential energy. The structure in which potential energy converged in six times from the initial structure was taken.

Results
Angiotensin was simulated using ProtTorter . While other simulation programs are invisible and untouchable, this can check each structure in every step. This can directly manipulate. Energetically stable local minima structure which is still smaller than the reference structure was observed.
The structure of 1n9v from the simulation with ProtTorter is shown in Fig. 5. Loop structure was found in most of the 8 residues except the rst aspartate residue in the angiotensin peptide. The dihedral angles(ϕ, ψ) were usually around the range of (25º-30º, 0º-5º). We can see the vivid difference between this and other simulations (Fig. 3). Experimental structure regularly oscillated above and below the 0º, while simulated angles were positive in the N-terminal region. These deviated far from the C-terminal region.
After simulating, we obtained angles and local energy minima. And all of these were arranged into Table  1 for clarity. The initial and six iterations were arrayed in the row of the table and seven bonds in the column. In the end, 49 arguments were placed in the table. The ψ angle is written on the top of the cell and the φ in the bottom. In each cell, angles and the number of local minima are written.  (4) 27 (9) 29 (9) 26 (8) 26 (8) 26 (8) 26 (8) 4th bond Ψ 1(2) 2(6) 0(8) 4(8) 3(7) 3 (7) 3 (9) Φ -150 (7) -39 (10 86 † initial structure generated from cotranslational folding with torsional energy calculation ‡ optimized structure following the folding path determined by the torsional propensity The structure was simulated from very large search space. As Table 1 shows, initial structure of cotranslational folding was from about 52 conformations. For this peptide composed with 8 residues, the most stable structure was from usually about 70-90 structures. It is very e cient to compare with typical molecular dynamics. logPr value (Table 2) signi es the difference of the two compared structures with more weight on the more closer similarity. There are eight dihedral angles each for residue from 1 to 8 in both reference and simulated structures. When the two angles of the same residue is very similar, those values were made to be equivalent. On the contrary, they were made to be different. RamRM SD 0.00 † initial structure generated from cotranslational torsional folding ‡ optimized structure following the sequence of folding based on the torsional propensity * experimentally determined structure The logPr values of simulated structures increased toward the later iterations implying the convergence in optimizations. The lowest logPr value of -15.18 was observed from the pair of 5th and 6th optimization. The fact that lower logPr values in the pairs of nearer iterations than the farther ones was found. The highest logPr value among the pairs of adjacent passes was − 4.62 in init. and opt. 1 (Table 2). There were difference between the cotranslational path and the torsional one.
RamRMSD (Table 2) is the RMS(root mean square) deviation between the positions of residues on the Ramachandran plot (Ramakrishinan and Ramachandan, 1965). RamRMSD is similar to logPr. This includes the growing similarity among later iterations. Pairs in closer passes had lower RamRMSD values than farther passes. The highest RamRMSD among the pairs of adjacent passes was 47.17. This was calculated from the pair between the initial and the rst optimization pass.
In Fig. 4, the change of energy in the folding of initial structure and in the optimization were illustrated.
The potential energy drastically uctuated in the simulation of initial structure in cotranslational folding. This partly indicates that addition of amino acid is either favorable or unfavorable in each different circumstance. This uctuation is different from following iterations. This re ects the strong effect of the change of con gurations. During the six passes of optimizations, the potential energy decreased saltatorily. This shows that there are a few critical bonds which strongly in uence the potential energy of the whole molecule. Demonstrating the fast convergence of the algorithm to the global energy minima, the potential energy remained as being conserved after three passes of optimizations.
Comparing this simulation and others' experimental structure of 1n9v with RamRMSD and logPr, this is more stable than others by global minimum of − 1.704(kcal/mol). The most correlation between each generated structure from its initial structure with adjacent passes have displayed. This increased for the later rounds of iterations. During the folding simulation, the energy dropped saltatorily (Fig. 4).
The structure from NMR spectroscopy was very different from this simulation. The average of all logPr values ranged from − 0.80 to -1.01. RamRMSD varied from 116.68 to 136.28. Although the simulated structure is somewhat different from the reference experiment, it is quite appreciable regarding the low and negative potential energy of -1.704(kcal/mol). This negative potential energy remark that this structure is stable in the vacuum environment. This structure is not only a low and stable energy structure but also a possible actual energy minimum because it is an energy minimum along the torsional propensity path. The torsional propensity path must be the path from Levinthal paradox.

Discussion
Results suggest ve parts to be discussed as follows. First, the difference of structure was due to the electrostatic interaction of atoms and the torsional barrier of rotatable bonds. Given motive force, a stronger turn is induced. This shortens the length of the loop structure. The structure of α-helix was observed from lattice model without the consideration of any detailed electrostatic or torsional potential energy (Leach, 2001). Thus, additional restraints of non-electrostatic interaction would induce the current loop structure into well-known helices.
Second, another reason for the difference is the utilized force eld. There was difference between NMR spectroscopy and this simulation. It is because that was conducted within an aqueous solution and this was performed under the vacuum environment. The difference between the simulated and the NMR structures was brought by the neglect of the interaction of solvents with the protein molecule. And hydrophobic effect and free energy from solvent accessible surface area could be obtained from experiments. This could be applied to structure simulation.
Third, it is very interesting for its fast convergence of the iterations. Although there is a false convergence, It is quite fast nding converging structure in 6 passes. Converging energy minima were quickly obtained following this method.
Fourth, simpli ed representation of the three dimensional structure of a protein in torsional system was applied. This regenerates the movements of atoms of polypeptide chain in the cellular environment. The fundamental characters of ribosome bound cotranslational folding could be generated. Three dimensional information can be transformed into one dimension by computing easily. This could be operated with sequence alignment algorithms in personal computer fast and correctly as BLAST (Altshul, 1990).
Fifth, ProtTorter adopted torsional representation of atomic movements . The results showed fast convergence to the stable form and which was negative and big in the potential energy. However, this path should be solidly validated referring longer polypeptide chains and larger numbers of test proteins.
The structure from this program is different from that of the representative NMR in torsion angle. This occurred in folding pathway or force eld.

Conclusion
Initial structure formation corresponded cotranslational protein folding. Torsion angels of the residue on backbone rotates one by one following torsional inertia. Optimization based on torsional inertia has reduced the number of candidate structures for the native structure. Stable conformation with low potential energy was obtained after input angiotensin sequence into ProtTorter which has the functions of torsion angle rotations and potential energy calculation. Coordinates of local minima to the stable structure were observed in very short time. Torque between two residues. The torsional inertia is determined from the length of the rotational radius and the mass of the object. The components of torque in the polypeptide chain are displayed.

Figure 2
Difference of torque among various bond angles. The magnitude of the torque varies according to the angles of bending of the rotating from stalled residues. The torque is maximized as it turns orthogonally while it vanishes when it turns in parallel.  Change of Potential Energy during the Initializations and Optimization. The change of potential energy during the initial structure generation and further iterative optimization is shown in blue line. The potential energy shows drastic uctuations during the initialization of cotranslational folding which partly indicates that addition of amino acid is not always either favorable or unfavorable. This uctuation is different from following iterative optimizations re ecting the strong effect of the change of con gurations. During the optimization processes, the potential energy decreased with saltatory tendency. After three iterations of optimizations, the potential energy remained as being rather conserved indicating the fast convergence of the algorithm to the local energy minima.

Figure 5
Structure of 1n9v. The 3D stick model of the structure of the angiotensin peptide(PDB entry 1n9v; "DRVYIHPF") is displayed. Loop structure was found among the most of the peptide(from 2nd to 8th residue) and drawn with blue cylinder. Carbon atoms are shown with gray color, nitrogen atoms with blue color, oxygen atoms with red color, and hydrogen atoms with white color. Some hydrogen atoms are shown as single spheres ignoring the bondage information.