Learning Correlations between Internal Coordinates to improve 3D Cartesian Coordinates for Proteins

We consider a generic representation problem of internal coordinates (bond lengths, valence angles, and dihedral angles) and their transformation to 3-dimensional Cartesian coordinates of a biomolecule. We show that the internal-to-Cartesian process relies on correctly predicting chemically subtle correlations among the internal coordinates themselves, and learning these correlations increases the fidelity of the Cartesian representation. This general problem has been solved with machine learning for proteins, but with appropriately formulated data is extensible to any type of chain biomolecule including RNA, DNA, and lipids. We show that the internal-to-Cartesian process relies on correctly predicting chemically subtle correlations among the internal coordinates themselves, and learning these correlations increases the fidelity of the Cartesian representation. We developed a machine learning algorithm, Int2Cart, to predict bond lengths and bond angles from backbone torsion angles and residue types of a protein, and allows reconstruction of protein structures better than using fixed bond lengths and bond angles, or a static library method that relies on backbone torsion angles and residue types on a single residue. The Int2Cart algorithm has been implemented as an individual python package at https://github.com/THGLab/int2cart.

INTRODUCTION Biomolecular structures are described using two widely used mathematical representations: internal coordinates and Cartesian coordinates. The internal coordinate representation is defined by a set of bond lengths, bond angles, and dihedral or torsion angles, and provides a compact description in terms of the Z-matrix. In contrast, a Cartesian representation defines all of the atomic positions in Euclidean x,y,z coordinates and additionally captures the orientation of a molecule in space. Both representations are useful in certain contexts and applications. Internal coordinates can be beneficial for geometry optimizations 1 and are the preferred description for NMR structure determination and refinement in which the bond lengths and bond angles are generally taken as fixed 2 , while Cartesian coordinates are the preferred format of molecular dynamics simulations 3 and X-ray crystallography structures provided in the Protein Data Bank (PDB) repository 4 . Figure 1 considers the internal coordinates of a protein backbone that contains the three torsion angles φ (C − N − C α − C), ψ (N − C α − C − N ), and ω (C α − C − N − C α ), bond lengths N − C α (d 1 ), C α − C (d 2 ), and C − N (d 3 ), and bond angles N − C α − C (θ 1 ), C α − C − N (θ 2 ), and C − N − C α (θ 3 ); side chain information that may effect the backbone could also include C α − C β (r 1 ) and N − C α − C β (α 1 ) for example. When all of these quantities are specified exactly, the back-transformation from internal coordinates will also result in a perfect 3D Cartesian reconstruction of the protein backbone structure, using algorithms such as the natural extension reference frame (NeRF) 5 . However, the Cartesian reconstruction is almost universally defined by only the backbone torsions while holding the bond lengths and angles fixed at mean values to decrease the complexity of structure modelling. Sometimes the variations on the ω torsion angles are also ignored and taken as fixed values of 0°or 180°, due to the planar nature of the peptide bond. 6 One might assume that a protein structure can be reconstructed in Cartesian coordinates quite well utilizing fixed bond lengths and angles since they typically have quite small variations around their means. However, even small deviations from the mean of the stiff degrees of freedom can strongly influence the Cartesian reconstruction. The origin of this error arises especially clearly from the nature of chain molecules: as the protein chain gets longer, small errors in bond lengths and bond angles can quickly accumulate and result in significant differences in the final back-transformed structure. According to a study by Holmes et al. 7 on globular proteins, the RMSD errors incurred in the internal coordinate back-transformations to C α Cartesian positions under fixed bond lengths and angles is ∼6 Å for an average 150-amino-acid protein, and can be as high as 40 Å for larger proteins.
Alternatively, one could replace the assumption of fixed bonds and bond angles with a statistical approach that conditions on sequence or structural correlations in the PDB. Given the many types of correlations that exist between the internal coordinates of globular proteins, such as the φ and ψ torsion angles of the Ramanchandran plot 8 and the correlation between backbone and sidechain torsion angles used in the Dunbrack rotamer library 9 , the correlations among the stiff bond lengths and angles with the flexible torsions should not be surprising. Earlier studies on the relationship between bond angles and φ, ψ torsion angles or amino acid types were mostly focused on the N − C α − C bond angle, using both statistical methods and quantum mechanics calculation on model dipeptides. [10][11][12][13] . The work by Berkholz 14 found that by using a static library for backbone bond angles dependent on backbone φ, ψ angles and residue types, the median RMSDs of protein reconstruction normalized to 100 amino acids is 2.85 Å. Following this pioneering study, more recent work by Roberto and co-workers have extended the correlation to all bond lengths and bond angles centered on backbone N, C α and C atoms. 15,16 Similarly, Lundgren et al. studied the correlation between protein backbone angles, secondary structure, and sidechain orientations. 17 , and Ashraya et al. evaluated the steric-clash Ramanchandran maps conditioned on bond geometries. 18 However, none of these studies have considered the ω torsion angle dependence of the internal coordinates. This work provides a more comprehensive machine learning approach that both quantifies and learns internal coordinate correlations within a deeper amino acid sequence context, that in turn provides a more accurate prediction of the 3D Cartesian coordinates relative to the errors incurred under the standard assumptions of fixed bond lengths and angles. By capturing the subtle correlations observed among internal coordinates, the Int2Cart (Internal to Cartesian) algorithm reduces the reconstruction RMSD error to ∼2.3 Å for test proteins less than 150-amino-acids, and an average RMSD of ∼3.5 Å over the entire test set for globular proteins as large as 543 amino acids. In addition we consider the "out-of-distribution" performance of Int2Cart on two more difficult test systems, including 42 proteins obtained from the 12th Community Assessment of Structure Prediction (CASP-12) as well as accurate reconstructions for intrinsically disordered protein (IDP) ensembles that we find generates lowers errors and fewer undesirable steric clashes. The Int2Cart deep learning model developed for proteins should be general to other types of chain molecules including RNA, DNA, and lipid molecules if appropriate data sets are made available.

METHODS
Dataset preparation. We have used SidechainNet 19 as a preprocessed dataset that uses clustering techniques to extract protein sequences and structures with defined similarity cutoffs, to reduce bias in the original PDB structures, and to prevent information leakage from the training set to the test set relevant to assessing the machine learning generalization. 19,20 The SidechainNet dataset represents each protein by its amino acid sequence, backbone and sidechain torsion angles (φ, ψ, ω, χ 1 , χ 2 , etc), backbone bond angles θ 1 -θ 3 , as well as the all-atom 3D coordinates. For this study we ignore the sidechain torsions as we are only reconstructing backbones, and also supplement the protein dataset with backbone bond lengths d 1 -d 3 calculated from the 3D coordinates for training, validation and test sets. We also identified some θ 2 and θ 3 bond angles that were incorrect due to missing atoms in the next residue, and they were masked out along with the residue at the end of the protein chain. We combined validation sets from 10% to 50% similarity cutoffs to serve as our primary test data. 20 When needed, we separated test set proteins at the broken chain positions and only retained chains longer than 50 consecutive amino acids. 20 Additionally we compiled two alternative and more difficult test data sets to validate the transferability of our model: the prediction targets from CASP-12 as well as reconstructions for IDP structural ensembles.
Our final training dataset contained 41,380 proteins with a minimum sequence length of 20 and maximum length of 4914 amino acids. The primary test set was comprised of 140 protein or protein fragments with sequence lengths between 23 and 543 amino acids. The challenging 42 proteins or protein fragments in the CASP12 test set ranged in length from 51 to 599 amino acids. The IDP ensemble comprised 1000 conformations for the N-terminal 92 residues of the Sic1 protein; 21,22 we extracted the backbone torsions and rebuilt the Cartesian structures for each conformer under different assumptions about the bond lengths and bond angles as reported in Results.
Neural network design. The structure of the deep neural network, Int2Cart, is depicted in Figure 2. We utilized 3 layers of stacked bidirectional gated recurrent units (GRU) as the central component, each of which contains a hidden state h t with its information updated by the reset and update mechanisms for each element in the input sequence through the following set of equations 23 : are the trainable parameters of the model, x t is the input to the cell at the current timestep, and r t and z t represent the reset and update gates, which are numbers between (0, 1) that control how much information to retain in the new update vector h t and how the new hidden state vector h t is composed from the update vectorh t and the old hidden state h t−1 . σ denotes the sigmoid function, and represents element-wise multiplication. Dropout was applied to the hidden states between layers, so that x is a Bernoulli random variable that zeros out elements in the hidden state vector with a probability defined by the dropout rate.
The inputs into the first layer of GRU cells are the φ, ψ and ω torsion angles and the amino acid types. Each torsion angle, a, was represented by a Gaussian smearing function discretized to a vector of length 180 to account for uncertainty in the data, denoted x ia where α i is the actual φ, ψ or ω angle andx a = (−180 + 2 * a) (both in degrees), and in this work ensures that the periodicity of the angles is taken into account. Each smeared torsion angle vector is further transformed through two fully-connected layers with 90 and 64 units each and Rectified Linear Units (ReLU) activation 24 to generate latent representations of the torsion angles. The residue types are encoded by a trainable embedding dictionary and formulated into latent vectors of length 64. The torsion angle latent vectors and the embedded residue types are then concatenated and transformed together through 2 fully-connected layers with 128 and 64 units and ReLU activation and constitute the inputs into the GRU cells.
The hidden state output from the last GRU layer is connected with multiple outputs to predict the backbone bond lengths and bond angles, or optionally sidechain bond lengths and bond angles as well. Each output is a fully-connected neural network with a hidden layer of size 100 and activation ReLU, and the output has size of 1 without any activation. The raw outputs are scaled by the standard deviation and translated by the mean value of that data type in the training dataset. The means and standard deviations we used are provided in Supplementary Table 1.
Training details of the Int2Cart machine learning method. The neural network was trained by minimizing the weighted mean square error loss function where w i controls the weighting for different data types in the loss function, y i are the predictions from the model andŷ i are the actual values from the data set. In practice we used the same weighting for all the data types. Missing data targets were masked out during the training. We used the Adam optimizer 25 with an initial learning rate of 0.001 and an exponentially decayed learning rate schedule, so lr i = exp(−i * α) where i is epoch number and α = 0.05 in our case. The model was trained for a total of 100 epochs using a batch size of 128.
Building all-atom Cartesian structures from internal angle model predictions. With the full profile of backbone torsion angles and predictions of bond lengths and bond angles from the model, the 3D Cartesian structure of the protein containing all backbone atoms is reconstructed using the SidechainNet package. 19 It utilizes the natural extension reference frame (NeRF) algorithm 5 to sequentially calculate the position of the next atom with the positions of three previous atoms and the new bond length, bond angle and torsion angle. The all-atom backbone Cartesian structures for all the protein fragments in our test data set are built from either the Int2Cart algorithm vs. a standard baseline of using fixed bond lengths and bond angles (Fixed) with no conditioning on torsions or amino acid type.

RESULTS
Statistical analysis of the protein training set. Given the large collection of protein structures deposited into the PDB, we first consider a statistical analysis of protein bond lengths and bond angles when analyzed over the training set. Overall the distributions of these internal coordinate values are mostly Gaussian with relatively small standard deviations of ∼ 0.01 Å for bond lengths and ∼ 2.6°for bond angles. Figure 3 a-f depicts the deviations from the mean bond length and angle values for a given (φ, ψ) combination, and confirms the existence of strong correlations among the internal coordinates averaged over the training data set. Specifically, the θ 1 angle is larger for the right-hand and left-hand helix regions in the Ramanchandran plot, while the beta-sheet regions have more narrow θ 1 angles, with deviations from the mean as large as 7.5°. The θ 2 values are strongly correlated with the ψ torsion angle, with larger angles when ψ is between −100°and 0°, and smaller angles than the mean otherwise. The θ 3 values for nearly all of the (φ, ψ) combinations are larger than average, but have smaller angles for helix regions. Similarly, the d 1 and d 2 bond lengths show greater correlations with the φ torsion angle, with a preference for larger values when φ is between −50°and +50°, in which the bond lengths change by as much as 0.02 Å. Finally the correlation for the peptide d 3 bond with the backbone torsions is weak, consistent with its partial double bond character, except for a few hot spots where it can vary up to 0.04 Å. These correlations are statistically meaningful, because the standard deviations in each bin are smaller than the mean value differences (Supplementary Figure 1), which means the statistical bias is more significant than the variance. We have also considered the relationships between backbone ω torsion angles with bond angles (Figure 3 g-i) and bond lengths (Supplementary Figure 2), and found interesting correlations between internal coordinates and ω torsion angles. The majority of peptide bonds in proteins are in trans-conformation, with ω torsion angles close to 180°. However, cis-peptide bonds tend to associate with smaller θ 1 angles and larger θ 2 and θ 3 angles. This result also makes structural sense since cis-peptide bonds incur more steric repulsion between sidechains of two consecutive residues, and larger θ 2 , θ 3 and smaller θ 1 values allow the sidechains to be more separated. On the other hand, the correlation between bond lengths and ω torsion angles are not obvious (Supplementary Figure 2). These correlations dependent on ω are also important for accurately predicting internal coordinates from backbone torsion angles, as we will show later.
When we consider the observed distributions of all six internal coordinates conditioned on the residue type (Supplementary Figure 3), we find that the distributions are quite similar between amino acids with only subtle differences in the shape of the peaks, with the exception of glycine, which tends to have d 1 and d 2 values that are smaller, and θ 1 angles that are larger than other residues. Proline also defines an exception, with larger d 1 values due to the formation of the fivemembered ring that requires longer bond lengths. However the bond length and angle distributions when conditioned on backbone torsions and residue type exhibit notable variations across all twenty amino acids as seen in Figure 4 for the θ 1 bond angle, as well as for the other backbone bonds and angles provided in Supplementary Figures 4-8.
Machine learning of sequence and structural correlations. While the correlation graphs just described could serve as a source for bond lengths and angles when backbone torsion angles and residue types are provided, we are still missing the deeper sequence correlations that are buried beneath the statistics of the single residue results. Therefore, we trained a deep neural network on the same data in order to capture the more subtle correlations among the internal coordinates conditioned on amino acid sequence. After training, the Int2Cart neural network was used to predict the primary test set which has at most 50% sequence and structural similarity with the training proteins. The results are summarized in Table 1 in terms of root-mean-square error (RMSE) and Pearson correlation coefficients (R). We find that the bond length predictions are within computational uncertainty as measured by the variance in the data set, while predictions on the bond angles are more successful in terms of the RMSEs that are smaller than the dataset variance. Table 1: Int2Cart prediction accuracy on backbone bond lengths and bond angles. Accuracy is assessed in terms of root-mean-square error (RMSE) and Pearson correlation coefficients (R).

Data type
0.007 Å 0.50 Cartesian coordinate reconstructions. Given the three torsion angles [φ, ψ, ω] for each residue over the entire protein sequence as input, we next consider how well the Cartesian coordinates are reconstructed based on whether bond and angle geometries are held fixed or learned from Int2Cart. Figure 5a reports the distribution of root-mean-square deviations (RMSD) for all backbone atoms with respect to the actual PDB structure for all proteins in the test set using Int2Cart and the Fixed method. In general, the reconstructed RMSD distributions for the Int2Cart structures are centered around lower RMSD values compared with the Fixed structures, with a median RMSD of 2.14 Å when normalized to 100 amino acids, and the average RMSD over the entire test set is 3.5 Å. By contrast the Fixed model yields a median RMSD of 3.24 Å when all proteins are normalized to 100 amino acids, and the average over the entire test set is 5.1 Å. The Int2Cart results are notably better than previous results, including the Berkholtz study 14 that conditioned on φ, ψ and amino acid type, showing that deeper sequence correlations are beneficial to learn. Figure 5b illustrates that the proteins reconstructed by assuming fixed bond lengths and bond angles have lost significant secondary structure integrity compared to the reference structures, whereas the Int2Cart structures retained a much higher proportion of secondary structural elements. The pairwise RMSDs for the test proteins reconstructed by Int2Cart are compared to those reconstructed using the Fixed approach in Figure 5c, and the RMSD difference between these two methods as a function of sequence length is shown in Figure 5d. It is evident that the vast majority of the test set proteins benefit from the machine learned bond lengths and bond angles, with an average improvement of 2 − 4 Å RMSD over using Fixed bond lengths and bond angles. There is no obvious correlation between the RMSD improvements made by Int2Cart over Fixed in regards sequence length, although the largest improvements occur in those proteins with longer amino acid sequences.
Beyond the anecdotal case in Figure 5b, we performed more extended analysis of Int2Cart and Fixed performance in regards the radius of gyration (Rg) and secondary structure recovery rate (SS-match) over the whole test dataset. Although we find that the Int2Cart Cartesian predictions have closer R g values to the ground truth structures, the Fixed Cartesian structures still yields a comparably good result as seen in Supplementary Figure 9. However Figure 6a shows that Int2Cart systematically improves upon Fixed in regards the SS-match values with a higher proportion of test set proteins that have SS-match values larger than 0.8 for Int2Cart predicted structures (Figure 6b), which translates to more than 80% of the residues having correct secondary structure assignments.
Ablation studies. To understand the importance of various inputs for prediction accuracy of Int2Cart and the importance of the bond length and bond angle targets in reconstructing Cartesian structures, we performed an ablation study by training separate deep learning models using subsets of the inputs, and reconstructing structures using only predicted bond lengths, only predicted bond angles, or using both. We have also trained models with additional inputs of χ 1 torsion angles, along with r 1 and α 1 sidechain bond lengths and bond angles as additional outputs, to evaluate how including sidechain information could improve prediction and reconstruction of backbone structures. All ablation trials are reported in Table 2.
The differences in predictions of the backbone bond lengths from different deep learning models are not significant, but prediction accuracy for backbone bond angles RMSE and reconstructed Cartesian structure RMSD are highly dependent on what information is available to the model. Specifically, a machine learning model that only knows about the residue types performs the worst with 5 Å in the reconstruction RMSD. Unsurprisingly based on statistical analysis of the PDB, backbone φ and ψ torsion angles provide more information than residue types alone, and allows the reconstruction RMSD to decrease to 4.24 Å on average. Including both φ, ψ and residue types further decreases the average reconstruction error to 4.05 Å. As expected from the correlation of bond lengths and bond angles with ω torsion angles as well, including exact values for ω torsion angles also significantly improves the model and allows the reconstructed structure RMSD to decrease further to 3.56 Å across the whole test set, and to 2.14 Å for proteins normalized to 100 amino acids. When we tested the inclusion of sidechain χ 1 torsion angles, we find that the 3D reconstruction model is even better, achieving an average reconstruction structure RMSD of 3.05 Å. That is probably because the χ 1 torsion angles are indicative of avoidable steric clashes between protein backbones and side chains to create more accurate descriptions of subsequent backbone bond geometries, even though side chain atoms are not explicitly treated during structure reconstruction. Table 3: Ablation studies of Cartesian coordinate reconstructions. Cartesian structure reconstruction RMSD using fixed bond lengths and angles, using predicted bond lengths and fixed angles, fixed bond lengths and predicted bond angles, and both predicted bond lengths and bond angles.

Source of bond lengths & bond angles
Reconstructed structure RMSD (Å) Fixed bond lengths and angles 5.07±3.47 Predicted bond lengths and fixed bond angles 5.06±3.46 Fixed bond lengths and predicted bond angles 3.51±2.65 Predicted bond lengths and bond angles 3.51±2.65 Table 3 shows the Cartesian structure RMSDs for protein structures reconstructed with predicted bond lengths, predicted bond angles and both, and compared with using fixed bond lengths and bond angles. Interestingly, the reconstruction quality does not depend on the direct prediction of bond lengths, as it essentially has no effect on the reconstructed structures, which my have been anticipated from the fact that bond length errors are on par with the variance. But this final ablation study provides direct evidence that accurate predictions of bond angles are of primary importance for the quality of the reconstructed Cartesian structures.
Evaluation of the method on more challenging test datasets. Finally we consider two independent test data sets that are quite different from the originally defined test set from SidechainNet, performed to evaluate the transferability of the Int2Cart model. Table 4 reports the prediction RMSE for all bond length and bond angles targets of the CASP-12 protein data set and the mean Cartesian reconstruction error. Compared with the previous test dataset results, the predictions on θ 3 becomes slightly worse in terms of RMSE and correlation coefficients, but all other prediction targets are very close to the original results described in Table 1. The reconstructed structure RMSD on the CASP-12 test dataset is on average 4.5 Å, but this shift to higher RMSD values relative to the original test set is due to the higher proportion of longer proteins. This is supported by the fact that the reconstructed RMSD as a function of sequence length for the CASP-12 data are indistinguishable from the primary test dataset as illustrated in Supplementary Figure 10. Additionally, the reconstructed RMSD for proteins normalized to 100 amino acids in the CASP12 test set is 2.06 Å, in close agreement with the results of the original test set. Table 4: Int2Cart prediction accuracy and reconstructed Cartesian structure evaluated on the more challenging CASP12 dataset. Accuracy is assessed in terms of root-mean-square error (RMSE) and Pearson correlation coefficients (R) for bond and angle predictions, and Cartesian reconstruction evaluated as root-mean-square-deviation (RMSD) with respect to the reference PDB structure.
Reconstruction RMSD all CASP12 proteins 4.49 Å CASP12 proteins (normalized to 100 aa) 2.06 Å In the second stress test we show that our Int2Cart method can improve upon the Cartesian reconstruction of an ensemble of structures of a disordered protein compared to Fixed bond lengths and angles. Figure 7 compares the Cartesian reconstruction RMSD distributions for Int2Cart and Fixed for the Sic1 IDP ensemble, in which we find that the Int2Cart method is overall closer to the original ensemble, with a 3.1 Å average RMSD compared to the Fixed method that has a mean RMSD of 3.4 Å. We have also checked the number of steric clashes in the structures generated from these two methods. A steric clash is defined as two atoms in the structure that are closer to 0.6 times the sum of the van der Waals radii of the two atoms. 28 Out of the 1000 conformations, 73 structures generated from Int2Cart contained steric clashes, which means 92.7% of the structures are clash-free. By comparison, 102 structures generated using fixed bond lengths and bond angles contained steric clashes, which translates to 89.8% of clash-free structures. A higher proportion of clash-free structures is meaningful because typically structures containing clashes are discarded, and a method with higher proportion of clash-free structures wastes less computational resources, and supports the application of the Int2Cart algorithm to disordered protein ensembles.

DISCUSSION AND CONCLUSION
In this work we have developed a new machine learning approach to the generic representation problem of internal coordinates (bond lengths, valence angles, and dihedral angles) and how to increase the fidelity of the back-transformation to 3D Cartesian coordinates. The Int2Cart algorithm utilizes a gated recurrent unit neural network to predict real-valued backbone bond lengths and bond angles for each residue of a complete protein sequence given its torsion angle profile. In summary Int2Cart can reconstruct the Cartesian structure of 95% of test proteins with RMSDs that are significant improvements over the fixed backbone bond lengths and bond angles that are the standard practice in a large variety of protein modelling approaches. The success of our algorithm on the more dissimilar CASP-12 prediction targets and across IDP ensembles further validates that the Int2Cart algorithm is transferable among different types of proteins, and can consistently improve the quality of Cartesian structure reconstruction.
In its current form the Int2Cart algorithm only generates backbone structures for the target proteins, although we can improve Cartesian reconstruction performance with the inclusion the χ 1 torsion and predicting r 1 and α 1 . Theoretical approaches such as the Monte Carlo Side Chain Ensemble (MC-SCE) approach can utilize the backbone from Int2Cart to calculate side chain ensembles in order to complete the full structure. 29 It is also clear that there is still room for improvement in the Cartesian reconstruction of larger proteins, and the inherent scaling of error with respect to sequence length is inevitable for a deep learning model that predicts internal coordinates in a sequential manner (i.e. a GRU model). Therefore, it may be possible to improve the quality of Cartesian structure reconstruction with a distance-based neural network model, i.e., by representing the 3D coordinates of the structure directly. Nevertheless, the model in its current form already provides a useful computational tool to greatly improve the quality of protein structures reconstructed from backbone torsion angles only, regardless of globular folded proteins or disordered protein ensembles. We envision Int2Cart should see broad use in structure refinement and validation 30,31 and development of protein force fields that could benefit from more accurate valence models of backbone bond lengths and bond angles conditioned on other geometrical or sequence features. 32 Finally, the Int2Cart GRU neural network model could also be useful for other chain molecules, only requiring retraining with new data if available for systems such as nucleic acids and lipids.       Supplementary Figure 10: Comparison for reconstructed structure RMSD as a function of sequence length for primary test datset and CASP12 test dataset. The higher reconstructed structure RMSD for structures in CASP12 test dataset is most likely due to higher proportion of longer proteins compared with results of the primary test dataset.