Estimating the probability of coincidental similarity between atomic displacement parameters with machine learning

The location of an atom is usually defined by its coordinate vector μ and its symmetric (ideally) semi-definite ADP tensor U once the resolution of the crystallographic data is better than 1.2 Å and the atoms appear as separate electron density peaks. Crystallographic refinement algorithms may or may not enforce meaningful U. For the purpose of this analysis, non-positive definite atom records are ignored. B_eq and B_est are then defined as [14]:

$$\begin{equation}{B_{{\text{eq}}}} = \frac{8}{3}{\pi ^2}\left( {{E_1} + {E_2} + {E_3}} \right),\end{equation} \tag{ 1 }$$

$$\begin{equation}{B_{{\text{est}}}} = \frac{8}{3}{\pi ^2}\sqrt {\frac{{{E_1} + {E_2} + {E_3}}}{{E_1^{ - 1} + E_2^{ - 1} + E_3^{ - 1}}}} ,\end{equation} \tag{ 2 }$$

where E₁, E₂ and E₃ are eigenvalues of U. U is equivalent to the covariance matrix of a multivariate normal distribution. Using the multivariate normal description of the atoms, the probability of their positions can be estimated in the crystal structure. The graphical representation of U is an ellipsoid. The eigenvectors and eigenvalues of the U tensor define the axis directions and their relative lengths, respectively. The surface of the ellipsoid is drawn at an arbitrary probability level of the multivariate normal distribution. The machine learning treatment of correlated diffraction intensity observations [2] and atom positions is very similar.

The Bayesian model (equations (3)–(6)) approximated U matrices as fixed Wishart distributions [8]. The Wishart distribution is frequently used as a conjugate prior to covariance matrices. The prior distributions are described in equations (3) and (5). Equation (4) shows a deterministic transformation and equation (6) defines the probability of the observation U. Equations (5) and (6) represent the Wishart probability distributions:

$$\begin{equation}p\left( {{\nu _0}} \right) = \sqrt {\frac{2}{{{\sigma ^2}\pi }}} {e^{ - \frac{1}{{2{\sigma ^2}}}{\nu _0}^2}} = \sqrt {\frac{2}{{9\pi }}} {e^{ - \frac{1}{{18}}{\nu _0}^2}};\,{\nu _0} \in \left[ {0,\,\infty } \right],\end{equation} \tag{ 3 }$$

wher σ is the scale parameter of the half-normal distribution (σ = 3 a priori),

$$\begin{equation}\nu = {\upsilon _0} + 5.\end{equation} \tag{ 4 }$$

$$\begin{equation}p\left( {\mathbf{V}} \right) = \frac{{{{\left| {\mathbf{V}} \right|}^{\frac{1}{2}}}}}{{{2^{\frac{{15}}{2}}}{{{\Gamma }}_p}\left( {\frac{5}{2}} \right)}}{e^{ - \frac{1}{2}tr\left( {\mathbf{V}} \right)}};{\mathbf{V}} \succ 0\end{equation} \tag{ 5 }$$

$$\begin{equation}p\left( {{\mathbf{U}}{\text{|}}\nu ,{\mathbf{V}}} \right) = \frac{{{{\left| {\mathbf{U}} \right|}^{\frac{{\nu - k - 1}}{2}}}}}{{{2^{\frac{{\nu k}}{2}}}{{\left| {\frac{{\mathbf{V}}}{\nu }} \right|}^{\frac{\nu }{2}}}{{{\Gamma }}_p}\left( {\frac{\nu }{2}} \right)}}{e^{ - \frac{1}{2}tr\left( {{{\left| {\frac{{\mathbf{V}}}{\nu }} \right|}^{ - 1}}{\mathbf{U}}} \right)}} = \frac{{{{\left| {\mathbf{U}} \right|}^{\frac{{{\nu _0} + 1}}{2}}}}}{{{2^{\frac{{3\left( {{\nu _0} + 5} \right)}}{2}}}{{\left| {\frac{{\mathbf{V}}}{\nu }} \right|}^{\frac{{{\nu _0} + 5}}{2}}}{{{\Gamma }}_p}\left( {\frac{{{\nu _0} + 5}}{2}} \right)}}{e^{ - \frac{1}{2}tr\left( {{{\left| {\frac{{\mathbf{V}}}{\nu }} \right|}^{ - 1}}{\mathbf{U}}} \right)}},\end{equation} \tag{ 6 }$$

wher U are the ADP observations, V and ν are the positive definite-scale matrix and the degrees of freedom of the Wishart distribution, respectively. The rank of the matrix is k = 3. Γ_p is the multivariate gamma function (p = 3).

The relationship between the model parameters is illustrated in figure 1. In order to generate the ν₀ parameter a half-normal prior distribution was assumed with a σ parameter of 3. The prior scale matrix V is another Wishart prior distribution with an identity scaling matrix and ν parameter set to 5. The samples of the matrix V were generated by Bartlett decomposition [9] in order to maintain the symmetry and semi-definite positivity of its posterior samples. The ν₀ parameter of the observations was incremented by a fixed value of 5 to ensure that the Wishart likelihood function is always defined. The Wishart likelihood function in equation (6) has a degree of freedom ν and a scale parameter V/ν. The training set consisted of 150 protein atoms randomly selected from the crystal structure of a 0.93 Å resolution complex of bovine trypsin in a complex with a canonical inhibitor (PDB entry 2xtt) [15]. The remaining protein atoms formed the test set.

Zoom In Zoom Out Reset image size

The model was implemented with the Python library PyMC3 [16], and the posterior distribution was estimated with the No-U-Turn Sampling (NUTS) algorithm [17] using the default parameters. The tuning was performed for 2000 Markov Chain Monte Carlo (MCMC) iterations, and the joint posterior distribution was represented by the last 1000 iterations. Four MCMC chains were generated to evaluate the reproducibility of convergence. The ADPs of the randomly selected atoms were assigned to separate 3 × 3 matrices, as vectorization of covariance matrices was not implemented in version 3.6 of PyMC3. This likely limited the performance of the algorithm and increased the memory usage. A total of 150 observations were already sufficient for reproducible convergence to very similar posterior distributions even if a different set of atoms was selected. Another limitation of this implementation was that the MCMC algorithm could only use one processor core, and calculating four MCMC chains took 190 s on a Linux workstation (i7 3970X CPU at 3.50 GHz clock frequency). This runtime corresponds to 12 000 iterations and excludes initialization/compilation time.

2.3.1. Distance calculations

For the distance calculations, the ADPs were converted to dimensionless matrices by dividing the tensor element magnitudes by 1 Ų. The atomic positions were converted into dimensionless vectors by dividing the elements by 1 Å. The metrics between the dimensionless U_A and U_B were defined as:

• (a)  
  Euclidean distance is Frobenius's norm of the difference matrix between U_A and U_B :

  $$\begin{equation}{d_{\text{E}}} = \|{{\mathbf{U}}_{\mathbf{A}}} - {{\mathbf{U}}_{\mathbf{B}}}\|_F,\end{equation} \tag{ 7 }$$
• (b)  
  Riemann distance:

  $$\begin{equation}{d_{\text{R}}} = \sqrt {\sum\limits_i {\text{log}}{{\left( {{{\boldsymbol{\lambda}}_{\mathbf{i}}}} \right)}^2}} ,\end{equation} \tag{ 8 }$$

  wher λ_i are the joint eigenvalues of U_A and U_B (assuming common eigenvectors).
• (c)  
  Log Euclidean distance is Frobenius's norm of the difference matrix between log(U_A ) and log(U_B ):

  $$\begin{equation}{d_{lE}} = \|\log \left( {{{\mathbf{U}}_{\mathbf{A}}}} \right) - {\text{log}}{\left( {{{\mathbf{U}}_{\mathbf{B}}}} \right)\|_F},\end{equation} \tag{ 9 }$$
• (d)  
  Log det distance:

  $$\begin{equation}{d_{lD}} = \sqrt {\log \left( {\det \left( {\frac{{{{\mathbf{U}}_{\mathbf{A}}} + {{\mathbf{U}}_{\mathbf{B}}}}}{2}} \right)} \right) - \frac{1}{2} \times \log \left( {\det \left( {{{\mathbf{U}}_{\mathbf{A}}}} \right){\text{det}}\left( {{{\mathbf{U}}_{\mathbf{B}}}} \right)} \right)} ,\end{equation} \tag{ 10 }$$
• (e)  
  Wasserstein distance [18]:

  $$\begin{equation}{d_{\text{W}}} = {\left( {tr\left( {{{\mathbf{U}}_{\mathbf{A}}} + {{\mathbf{U}}_{\mathbf{B}}} - 2{{\left( {{{\mathbf{U}}_{\mathbf{A}}}^{1/2}{{\mathbf{U}}_{\mathbf{B}}}{{\mathbf{U}}_{\mathbf{A}}}^{1/2}} \right)}^{1/2}}} \right)} \right)^{1/2}},\end{equation} \tag{ 11 }$$
• (f)  
  Symmetric Kullback–Leibler distance [19]:

  $$\begin{equation}{d_{{\text{KL}}}} = \frac{1}{2}\left( {tr\left( {{{\mathbf{U}}_{\mathbf{B}}}^{ - 1}{{\mathbf{U}}_{\mathbf{A}}}} \right) + tr\left( {{{\mathbf{U}}_{\mathbf{A}}}^{ - 1}{{\mathbf{U}}_{\mathbf{B}}}} \right) - 6} \right),\end{equation} \tag{ 12 }$$
• (g)  
  B_eq distance:

  $$\begin{equation}{d_{\text{B}}} = \frac{8}{3}{\pi ^2}\left| {tr\left( {{{\mathbf{U}}_{\mathbf{A}}}} \right) - tr\left( {{{\mathbf{U}}_{\mathbf{B}}}} \right)} \right|,\end{equation} \tag{ 13 }$$
• (h)  
  Eigenvalue distance:

$$\begin{equation}{d_\lambda } = \|{\boldsymbol{\lambda} _{\mathbf{A}}} - {\boldsymbol{\lambda} _{\mathbf{B}}}\|,\end{equation} \tag{ 14 }$$

where λ_A and λ_B are the eigenvalues of U_A and U_B , respectively. The distance is defined as the vector norm of the differences between λ_A and λ_B .

In addition, the distance between atom positions was defined as:

$$\begin{equation}{d_{\text{r}}} = \|{\boldsymbol{\mu} _{\mathbf{A}}} - {\boldsymbol{\mu} _{\mathbf{B}}}\|,\end{equation} \tag{ 15 }$$

where μ_A and μ_B are the coordinates of atoms A and B. The distance is defined as the vector norm of the differences between μ_A and μ_B .

The matrix logarithm of U is defined as:

$$\begin{equation}{e^{{\text{log}}\left( {\mathbf{U}} \right)}} = {\mathbf{U}}.\end{equation} \tag{ 16 }$$

The metrics were implemented using the Python libraries NumPy [20, 21], SciPy [22] and pyRiemann [23]. The Python libraries Matplotlib [24], Graphviz [25] and Seaborn [24] were used for data visualization. Validation of the method using synthetic data is detailed in the supporting online information (figures S1–S4 and tables S1 and S2 (available online at stacks.iop.org/MLST/2/035033/mmedia)).

We determine the distance between two samples from the same Wishart random distribution repeatedly and evaluate what distance corresponds to a certain level of cumulative probability. In order to compare covariance matrices, there are many metrics are available. All definitions assume that for identical covariance matrices the distance is 0 and most definitions are commutative, i.e. the distance is the same in both directions. Perhaps it is most computationally efficient to calculate the Euclidean distance (d_E). A general case of the Euclidean distance is the Riemann distance (d_R) defined between two points in a Riemann manifold. Locally, the geometry of a smooth Riemann manifold approaches Euclidean geometry, which makes the two interchangeable when very similar ADPs are compared. Unfortunately, d_E and d_R started to diverge rapidly. Atom pairs with similar d_E can have different d_R and vice versa. These metrics do not have to be particularly large to observe this tendency. Figure 6 shows a comparison of different covariance distances between pairs of randomly selected atoms in the protein structure and between two samples from the simulated Wishart distribution, respectively. For d_E, d_R, d_lE, d_W and d_lD, the two distance probability distributions differ to a large extent. Our null hypothesis is that there is no relationship between ADPs of atoms and their similarity can be predicted by our statistical random model. Table 1 compares distances at different significance levels, i.e. at which distance the pairwise distances derived from our random model reach 1% and 5% cumulative density. When protein atoms are compared, the probability of finding pairs with very short distances is much higher. For example, 1% cumulative density is reached at distance d_E of 0.052 in the random model, while 5% cumulative density of protein comparisons has already reached d_E of 0.034. For any of these 5% atom pairs we can reject the null hypothesis. By putting cumulative density into context, we can focus on pairwise comparisons to a single, generic atom from all other atoms in the structure. The crystal structure contains 1876 protein atoms. If we accept 1% false positive pairwise relationships, 27.7% (520 atoms) of the protein atoms are significantly more similar to a single atom than expected by chance. This number of atoms is much larger than the first and second coordination shell around a typical atom and a substantial fraction of similarities develops over a relatively long distance. Clearly, a protein is more inhomogeneous than this simple model suggests. While some atoms have almost no significantly similar partners, other atoms belong to an even larger network of significantly similar atoms.

Zoom In Zoom Out Reset image size

The cumulative density of d_B at short distances is very similar when the experimental data and the random model are compared. It may be counterintuitive to realize that samples pulled from the same distribution do not have a maximum around zero for their pairwise distances and this is a characteristic of multivariate distributions. For the univariate B_eq distance (d_B) the most frequent pairwise distances are close to zero both for comparison of the random andexperimental data. This means that the high similarity of B_eq is not surprising and expected by chance, i.e. the similarity does not carry substantial information. The metric d_B is not good for detecting similarities. The distance distribution differences between the random and experimental model do diverge when distances are long, but in this study, we focus on the coincidental similarities not on coincidental differences.

Figure 7 shows the correlation between ADP metrics, and the distance between atom coordinates is also included (d_r). The diagonal plots represent the kernel density estimates of each distance distribution. The d_lD, d_lE, d_R and d_KL metrics are highly correlated with one another. Therefore, only one of them, d_R, was retained in the final analysis. The kernel density estimate of d_r represents the pair distribution (p(r)) function, which provides information about the shape of the molecule. When d_r is plotted against different ADP metrics, it is apparent that atoms can be far apart and yet have nearly identical ADPs, and they are much more similar than expected by chance. The converse does not have to be true. When atoms are close together, their ADPs do have increased similarity. This can be observed due to the relative lack of points under the diagonal. There are almost no atom pairs within bonding distances (1–4 Å), which have highly dissimilar ADPs. This local similarity could originate from crystallographic restraints, rigid body motions and other types of local constraints. There is also a low frequency of short distances in the p(r) function. Notably, metrics that are less correlated with d_B, such as d_R show more diverse atom pairs at closer interatomic distances. B_eq is a pure measure of displacement amplitude without any directional information. This indicates that the local ADP similarity is mostly restricted to the displacement amplitudes rather than to the displacement directions. This is in line with the general observation that the crystallographic B-factors tend to increase from the core towards the surface of proteins.

Zoom In Zoom Out Reset image size

The isotropic displacement amplitude distance d_B and the metrics based on the shape parameters (eigenvalues) of U (d_λ ) have poor correlation with other metrics, especially when the distances are short. This would normally be a favorable, complementary property, but figure 6 and table 1 show that short distances of d_λ and d_B occur easily by chance and do not carry substantial information.

Our earlier analysis identified the reproducible similarity of Cβ atoms in the amino acid residues His-57 and Ser-195 using a hierarchical clustering approach in trypsin in a complex with a small molecular inhibitor benzamidine [7]. We can now verify the robustness of the claim in a different crystal structure, by looking at a different protein-inhibitor complex refined by different crystallographic software, and put the results into a statistical context. In our earlier analysis, we only used a variant of Euclidean distance and one type of clustering algorithm. Hierarchical clustering tends to commit an atom early to a cluster, making it unavailable for other plausible groupings. This is especially problematic when the ends of the branches (small clusters) contain the most relevant information. Small changes in the metrics and algorithms used for clustering can result in a different tree. Thus, hierarchical clustering is excellent for capturing some patterns, but often misses links that are revealed in a more systematic analysis. In table 2, the similarities are ranked after increasing d_R. At this high level of similarity, d_E and d_W do not diverge strongly, but if the same atoms are sorted according to other metrics, the ranking changes slightly. Distances of these magnitudes are well below the 1% significance level (table 1). The metric d_λ and especially d_B correlate poorly with d_R, indicating that the displacement amplitudes alone do not contain sufficient information. The interatomic distance d_r varies greatly in this group of highly similar atoms to Cβ of His-57. At rank 1 and 2, we find atoms 16.2 and 24.9 Å apart, respectively. Only at rank 3, do we find a covalently bonded atom (Cα of His-57). All other similarities arise in the absence of covalent bonding and van der Waals interactions. We find Cβ of Ser-195 from our previous study [7] at rank 23 in the company of other carbon atoms, and in particular Cβ atoms. Recurring atoms are observed from the amino acid residues Ser-195, Ser-214, Lys-87 and Ile-89. Residue Ser-195 is the most important member of the catalytic triad as it performs the first nucleophilic attack on the substrate peptide or ester and forms covalent tetrahedral and acyl-enzyme intermediates during the catalytic cycle. Ser-214 is often attributed to the 'catalytic tetrad' [26] because its mutagenesis drastically hampers serine protease catalysis. Amide nitrogen of Ser-214 forms a hydrogen bond of variable strength to the P1 residue of the substrate [27]. Its side-chain hydroxyl group (Oγ) maintains an unconventional CH–O hydrogen bond to His-57; thus, Ser-214 maintains a link between the substrate and His-57 [28].

Table 2. The shortest ADP distances to Cβ of His-57. The color scale red-yellow-green illustrates the range between the lowest and highest values in the columns.

We established a statistical framework for comparing ADPs and revealed that ADP similarities in bovine trypsin do not arise by coincidence. We have still not given a causal explanation. Unfortunately, we can only speculate about its origin if we consider the spatial distribution alone. We discuss here two kinematic processes: constrained displacements and coherence of motion. Both types of motion could result in identical positional distributions. We illustrate the difference with a theoretical experiment, where we assume that an atom is part of a planar chemical group. Suppose that this atom is constrained in its freedom. It cannot move in the plane direction and it displaces only in- and out-of-plane directions. Suppose that we have a number of identical chemical groups separated by long distances so that no direct contact between them is possible. If there is a yet unexplained process, which aligns the chemical groups parallel to one another, then the displacement of atoms will occur in similar directions. They do not have to move in phase; the constraints provided by the planar chemical group explain the uniform directional distribution. One problem with this model is that we do not have a universally accepted mechanism for protein folding, which explains why chemical groups should orient themselves in similar directions. This model also requires very strong constraints, for example, the out-of-plane displacement has to be perfectly perpendicular to the planar chemical group. Moreover, empirically we can observe that chemical groups are not perfectly aligned in the protein folded structure, and they appear to be affected by steric hindrances that prevent crystalline order. However, we cannot dismiss the presence of atoms with very similar displacements.

We prefer an alternative explanation because it does not require a speculative mechanism that orients the chemical groups over a long distance. Instead, we can assume that the discussed atom in the planar chemical group displaces coherently (in phase) with the same type of atoms in other groups, at least occasionally. Then, the planar constraints on the atoms would then work the other way around. The parallel orientation of the planar chemical groups becomes a consequence of the displacement alignment of a particular atom type.