Protein Conformational States—A First Principles Bayesian Method ^†

Abstract

Automated identification of protein conformational states from simulation of an ensemble of structures is a hard problem because it requires teaching a computer to recognize shapes. We adapt the naïve Bayes classifier from the machine learning community for use on atom-to-atom pairwise contacts. The result is an unsupervised learning algorithm that samples a ‘distribution’ over potential classification schemes. We apply the classifier to a series of test structures and one real protein, showing that it identifies the conformational transition with >95% accuracy in most cases. A nontrivial feature of our adaptation is a new connection to information entropy that allows us to vary the level of structural detail without spoiling the categorization. This is confirmed by comparing results as the number of atoms and time-samples are varied over 1.5 orders of magnitude. Further, the method’s derivation from Bayesian analysis on the set of inter-atomic contacts makes it easy to understand and extend to more complex cases.

1. Introduction

The conventional description of protein dynamics asserts that proteins posses intrinsic conformational states [1]. An enzyme may cycle between catalytic and open states [2]. An ion channel may open and close its central pore [3]. A chaperone protein assists transformation of large, hydrophobic proteins from initial, linear, to folded shapes [4]. X-ray and cryo electron-microscopy reveals conformations with small motions on the 1–5 Ångstrom level for those proteins that crystallize [5]. Neutron scattering and nuclear magnetic resonance structures of room temperature proteins show greater shape variability, but are usually able to classify structures into a few ‘canonical’ structures.

Advances in molecular modeling have made it possible to simulate the protein folding process, generating very large numbers of samples, free energy landscapes, and information on kinetics. Nevertheless, computational identification of distinct conformational states from molecular simulations has remained an active area of methodological research. Approaches fall roughly into two classes—linear and nonlinear. Linear methods work with the real $3N$-dimensional vector of atomic coordinates as a linear space, using pairwise alignments and Euclidean (RMSD) distances between structures. Nonlinear methods work with ‘features’ derived from nonlinear functions of the atomic coordinates—like pairwise distances between alpha-carbons. Well-known linear techniques include principle component analysis and linkage clustering based on pair-RMSDs. Nonlinear techniques include methods based on internal coordinates (relative distances and angles between groups of atoms) and T-distributed stochastic neighbor embedding [6].

Both linear and nonlinear techniques can be applied to subsets or groups of atoms to classify site-specific or inter-domain structures, respectively. This can be done trivially by handing off different point-sets to the analysis. Thus, most methods can be adapted to run with reasonable compute time and yield classifications most relevant to a motion under investigation. Automatically detecting protein secondary and tertiary structures, however, remains a challenging problem which we cannot address here.

One of the principle methods developed for visualizing domain motions in proteins is DynDom [7,8]. DynDom works from two input structures and determines relative domain rotations (using internal coordinates). This allows predicting transition motions from experimentally observed conformers, but not conformations from observed motions. On the other hand, many structure-to-structure similarity classification methods have emerged for this problem [9,10,11]. A recent review of such dimensionality reduction methods for protein conformational spaces noted that nonlinear methods are generally better than Cartesian or linear ones, but that the complexity of assumptions behind those models makes them difficult to work with and adapt [6,12].

This work presents a complete inference method derived from one single statistical hypothesis: that conformational states are defined by sets of contacting residues. Specifically, we hypothesize that the conformational state, k, uniquely determines which pairs of residues u, v, will be touching. Like a weighted coin flip, the contact probability is ${\mu }_{k,(u:v)}$—independently from all the other contacting pairs. Each conformational state is thus characterized by a vector, ${\mu }_k$, encoding the set of contacting pairs in state k.

The statistical model derived from this problem statement is termed a Bernoulli mixture model for binary feature classification [13]. The problem setup is similar to the Naive Bayes method [14]. However, because the categories are not known in advance, this is an unsupervised learning and classification problem.

Bernoulli mixture models have been applied extensively in the field of text subject analysis [15], optical character recognition [14], and image feature classification [13]. Essentially all of these applications have been successful at building extremely accurate classification models. The latter work also presents a thorough summary of sampling methods.

However, there remain difficulties sampling the distribution over categories, $\mu$, especially when the number of categories and reference classifications are not known in advance. The well-known expectation-maximization algorithm (EM) [16] is available in principle, but is not a replacement for sampling. Theoretical work on the EM method [17] shows that redundant categories will result in many circumstances. In this work, we have introduced a prior that eliminates redundant categorizations.

This work is structured as follows. Section 2 presents the underlying probability distribution of categories, then outlines a novel method for quickly sampling parameter space—achieving category inference. Full technical details are present in Appendix A.1 and Appendix A.2. Section 3 describes test problems on which the method is demonstrated. Section 4 presents results demonstrating that the method creates structurally meaningful categories with >90% accuracy. Although the potential application space is vast, this work focuses on proving method robustness using well-defined synthetic test problems. Each follows a time sequence mimicking domain motions in proteins—so that the classification accuracy can be judged by correctly assigning categories in time-order. For the practitioner interested in trying the method directly, full source code and scripts reproducing the test cases in this work are available under the GPLv3 license (Supplementary Materials or [18]).

2. Theory

A naïve Bayes model (for bit-strings) assumes that structural input samples, $i=1,\dots ,N$ are generated by first selecting a conformational state, $z_i\in \{1,\dots ,K\}$, with probability ${\pi }_{z_i}$, and then independently deciding whether each point-to-point contact ($x_{ij},j=1,\dots ,M$), is made with probability ${\mu }_{z_{i}j}$. If contact j is made in sample number i, then $x_{ij}=1$. The bit-vector, $x_i$, is said to posess feature j. Otherwise $x_{ij}=0$, and the feature is absent. The model parameters are thus $\theta =(K,\pi ,\mu )$.

It leads to a sampling distribution,

$$\begin{array}{}\mathrm{(1)}& P\left(xz\right|\theta )& =\prod _i{\pi }_{z_i}\prod _{j=1}^F{\mu }_{z_{i}j}^{x_{ij}}{(1-{\mu }_{z_{i}j})}^{1-x_{ij}}\mathrm{(2)}& & =\prod _{k=1}^K{\pi }_k^{N_k}\prod _{j=1}^M{\mu }_{kj}^{N_{kj}}{(1-{\mu }_{kj})}^{N_k-N_{kj}}.\end{array}$$

The second line above notes that, once the categorization, z, is known, the sampling distribution is easy to express in terms of feature counts in $\{i:z_i=k\}$—the set of samples assigned to category k,

$$\begin{array}{cccc}\hfill N_k& =\left|\right\{i:z_i=k\left\}\right|,\hfill & \hfill N_{kj}& =\sum _{\{i:z_i=k\}}{x}_{ij}.\hfill \end{array}$$

(3)

The first is the number of samples in set k, and the second is the number of times each contact is seen in that set.

According to Bayes’ theorem [19], we can turn this around to predict two important things—the probability that sample i belongs to category k, (read z given x and $\theta$),

$$P(z_i=k|x_i\theta )\propto \prod _{j=1}^M{\mu }_{kj}^{x_{ij}}{(1-{\mu }_{kj})}^{1-x_{ij}},$$

(4)

and also the probability distribution over all possible parameters,

$$P\left(\theta z\right|xI)=C(x\left)P\right(xz\left|\theta \right)P\left(\theta \right|I),$$

(5)

where $C\left(x\right)$ is an x-dependent normalization constant. Sampling this distribution provides everything—the categorizations, z, the conformational states, $\pi ,\mu$, and even a predicted number of categories, K.

In Bayesian probability, a prior distribution has to be assumed by the researcher. The prior characterizes the parameter space, independently from any sampled data. Our prior distribution over parameters, introduced below, is $P\left(\theta \right|I)$. Since the parameters directly determine the sampling distribution, the prior does not affect it [$P\left(xz\right|\theta )=P(xz\left|\theta I\right)$]. Note that this work juggles between two different priors, I and U, because the inference problem is simpler using $P\left(\theta \right|U)$, but $P\left(\theta \right|I)$ eliminates redundant solutions.

We choose a prior probability,

$$\begin{array}{cc}\hfill P\left(\theta \right|I)& \propto P\left(\theta \right|U)\prod _{k<l}\left(1-B({\mu }_k,{\mu }_l)\right),\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill P\left(\theta \right|U)& =\frac{\mathsf{\Gamma }\left(K\alpha \right)}{\mathsf{\Gamma }{\left(\alpha \right)}^K}\prod _{k=1}^K{\pi }_k^{\alpha -1}\hfill \end{array}$$

(7)

that enforces uniqueness of the categories. Here $B({\mu }_k,{\mu }_l)=B_{kl}$ is the Bhattacharyya similarity between distributions p and q,

$$B(p,q)=\prod _{j=1}^M\sqrt{p_j{q}_j}+\sqrt{(1-p_j)(1-q_j)}.$$

(8)

If two categories share the same distribution, then $p=q$, and $B(p,q)=1$. This forces our estimator to return zero likelihood that ${\mu }_k={\mu }_l$ for any $k\ne l$.

The second part of Equation (7), $P\left(\theta \right|U)$, is a conventional prior used for Bernoulli mixture inference in the literature. We use $\alpha =M+1$ throughout this work. Appendix A.1 contains a detailed justification for this choice. Essentially, it forces the likelihood for dividing a category into two parts to be asymptotically insensitive to the number of features, M. A proof of this fact, as well as a useful connection to the information entropy of compression is also present in Appendix A.1.

Sampling Method

Our sampling method is traditional Markov Chain Monte Carlo using four types of moves: a recategorization move, where categories, z, are assigned according to $P\left(z\right|\theta xI)$, a reclassification move, where $\theta$ is sampled from $P\left(\theta \right|zxU)$ and accepted with probability ${\prod }_{k<l}(1-B_{kl})$, and one split and one join rule. The function, $P\left(\theta \right|zxU)$, referred to here is just $P\left(\theta \right|zxI)$ without the Bhattacharyya distance terms, and with a different constant prefactor.

For the split trial move, one of the categories, k, is split into two new categories. Every member of k is re-assigned into one of two new sets, labaled L or R. Join moves are the opposite of split moves. This re-categorization changes the set labels, z. Specifically, $z_i$ goes from k to L or R for every i formerly at $z_i=k$.

Generating split or join trial moves was done by randomly choosing either one category to split or choosing two categories to join. For splits, member i of category k is moved to set L with probability ${\eta }^{x_{ij}}{(1-\eta )}^{1-x_{ij}}$. We used $\eta =0.9$, but any $\eta \in (0.5,1.0)$ should work in principle. If all elements end up in L or R, the partitioning is re-done. To concentrate splitting on productive cases, we did not attempt to split categories with $N_{kj}=0$ or $N_{kj}=N_k$. Immediately after splitting or joining categories, a reclassification move (re-assigning $\theta$) was performed. Category split moves were accepted using the Metropolis criterion, which is the smaller of $P_{\mathrm{gen},\mathrm{join}}{P}_{\mathrm{split}}/P_{\mathrm{gen},\mathrm{split}}{P}_{\mathrm{joined}}$ or one. Explicit formulas for the move generation probabilities ($P_{\mathrm{gen},\mathrm{split}}$, etc.) are provided in Appendix A.2.

Figure 1 provides a graphical summary of this inference scheme. Each conformational sample is mapped to a bit-string, which is used as the basis for inferring $\mu$. Inference proceeds by sampling potential parameters until a good explanation for the data is found. Trial moves that re-categorize and update $\mu$ look horizontally to find better category prototypes. Trial moves that split categories based on presence or absence of some features allow us to traverse category space vertically.

3. Test Systems

The ability of our sampling procedure to predict categories was tested on three geometrical systems: ‘chomp’ (a closing angle), ‘heli’ (a rotating line), and ‘glob’ (three rotating spheres). Each system was generated as a time-series of 1000 frames for approximately P total particles in 2 or 3-dimensional space. After generation, Gaussian random noise of width $\sigma =0.1$ was added to every degree of freedom.

Figure 2 shows images of these three test trajectories. A complete description of the coordinate generation methods is present in Appendix A.3. All three systems were processed into binary feature data by calculating pairwise distances between all points. Pairs of points within 2 distance units were translated to 1 (representing contact). When forming the feature vectors, x, we removed features, j, for which every sample showed the same result (all contacting or all disconnected). It is important to note that this removal changes M seen by the algorithm—usually decreasing it well below $P(P-1)/2$.

Critically, we repeated these classifications for a range of material points, P. This tested robustness of the unsupervised classification problem with respect to the amount of features available. Adding more points without changing the geometry of the problem should not change the number of categories detected.

For each run, five independent MCMC chains were started from an assignment of all points to a single category. Each chain ran for 1000 steps. Samples were collected every 10 steps—starting from step 500. The acceptance probabilities for category split/join Monte Carlo moves varied around 10–13%.

4. Results

We analyzed the results of MC in two different ways. First, the categories assigned were tested for grouping in time. Since the contact lists (on which the categorization is based) varied slowly over time, we expect categories to come in ‘runs’. Second, we computed histograms over the number of categories, K. This is a strong test of the method’s sensitivity to the number of material points, P.

As expected, we found a high degree of correlation between categories and time for every case. Similar time-points were grouped into similar categories. To quantify these results, we counted transitions between category indices in time-order. For a perfect categorization, the number of transitions should equal the number of categories minus one. We computed the categorization accuracy in two ways. For each system, the left columns of

Table 1. Categorization accuracy for each system type and number of points, P. The left column for each system is the percent of category assignments identical to their previous time (not including $K-1$ required). The right column for each is the percent of category assignments that do not switch between runs.

P	Chomp		Heli		Glob
18	89.2	98.0	89.9	98.9	91.1	96.7
42	94.7	98.9	93.7	99.0	93.5	97.3
78	96.3	99.3	95.3	99.4	95.5	98.0
114	96.8	99.4	96.5	99.5	96.4	98.4
150	97.2	99.4	96.7	99.4	97.0	98.5
222	97.8	99.5	97.0	99.7	98.2	99.1
258	98.4	99.7	97.7	99.5	98.6	99.4
330	98.4	99.6	97.8	99.7	98.7	99.4
438	98.5	99.7	98.1	99.5	98.8	99.4

and

Table 2. Categorization accuracy for each test system as the sample number, N, varies. The number of points remains fixed at $P=222$. Column labels are as in Table 1.

N	Chomp		Heli		Glob
125	99.3	99.2	99.8	99.2	99.5	99.3
250	99.4	99.6	99.3	99.6	99.2	99.4
500	98.8	99.6	97.7	99.7	98.6	99.3
1000	97.7	99.6	96.9	99.6	98.2	99.2
2000	96.9	99.4	96.5	99.6	98.0	99.1

are 100 minus the percent of mis-categorized frames. For lossy categorization at time-boundaries, we expect oscillation between two values. We quantified this by forming a transition matrix between categories, and removing transitions along the ‘most likely path’. Excluding this boundary oscillation we found nearly 100% accuracy for the categorizations. Those are shown in the right columns of Table 1 and Table 2.

Integrating the posterior probability (Equation (2) times Equation (7)) over $\pi$ leads to factors like $\mathsf{\Gamma }{(N+\alpha K)}^{-1}$, which seem prohibitively costly as N increases. We therefore wanted to check that the number of categories does not decrease as features or samples are added. Figure 3 shows the sampled probability distributions over K, the number of categories, for increasing P (Figure 3a) and N (Figure 3b). Interestingly, for every system, about five conformational states were deduced at $N=1000$. As N increased, however, more categories were deduced and the distribution spread to higher numbers. This is probably reasonable, since more values of the ‘time’ coordinate generated a more fine-grained motion.

Adenylate Kinase Open/Closed Transition

Finally, we tested the classification method against a protein with a well-characterized conformational transition. Adenylate kinase (ADK) converts ATP, ADP, and AMP by closing around substrate molecules [20]. The transition from closed to open was simulated in Reference [5] using steered molecular dynamics on a reaction coordinate interpolating between the electron density maps of PDB IDs 1AKE (Reference [20], closed) and 4AKE (Reference [2], open). The simulation data we used did not contain ligands, but did contain water and ions. Our analysis only made use of the alpha carbon (C${}_{\alpha }$) positions.

Features were calculated for each of $N=3900$ equally spaced frames during steered dynamics by testing whether C${}_{\alpha }$ to C${}_{\alpha }$ distances were less than 5 Å. These structures contain $P=214$ C${}_{\alpha }$-s. Sampling was carried out as described in Section 3, but 8 independent MC chains were sampled for 1250 steps (instead of 5 for 1000). The acceptance probability of split/join moves was 17%. During sampling, we saved the parameters, $\theta$, possessing maximum likelihood values at each K.

Our implementation (Reference [18]) is parallelized so that each thread carries out an independent Monte Carlo chain. The run-time for the ADK example was less than 30 min. Larger runs on proteins up to 400 residues using 75,000 sampled protein conformations have been carried out using this method in under an hour. The implementation also contains a series of vizualization tools to highlight the protein regions responsible for detected conformation-to-conformation differences.

We then extracted conformational states with the highest probability for landing in each category as representative points for that category. Since the reference open and closed PDB conformations formed extreme points, our representative structures approached them more nearly as K increased. Figure 4 shows that the two end-point conformations ended up very close to the open and closed states from the PDB.

5. Conclusions

The method developed here is ideally suited for the unsupervised structural classification problem. It has been derived from a first-principles Bayesian analysis of the set of atoms which interact within a structure. This work solved the central problem of defining an appropriate prior distribution over parameter space and implementing an efficient sampling method.

Tests on sample conformational transitions identified more categories than naïvely expected because it generated milestones along the motion’s time-coordinate. However, all categorizations were shown to have excellent accuracy as judged by picking out the correct time-sequence. Finally, a test on Adenylate kinase verified that these conclusions easily generalize to protein motions.

The central results of this work showed that the method behaves well under large variations in the number of features and samples. These observations validate our choice for the prior distribution, since changes in $P\left(\theta \right|I)$ will have large effects on the distribution over number of categories, $P\left(K\right|xI)$. They also show robustness of the MC sampling method itself, since relatively high acceptance rates were achieved.

Many future applications and the development of this method are possible. Changes in the pair contacts between states could be analyzed more thoroughly, as in the DynDom method [7]. Probabilities of assignment to each conformational state can be used as reaction coordinates. We are presently applying the method to classify chemical compound space using binary MACCS fingerprints [21], and to characterize conformational space of SARS-CoV-2 proteins simulated using replica exchange molecular dynamics [22].

Now that the principle has been demonstrated, more informative classification schemes can be devised. Adding information like hydrogen-bonding, salt bridge formation, and secondary structure annotation will allow the Bayesian framework able to recognize categories more like a biochemist would. The method can also be focused on active domains and binding sites by adding more points and shorter distance cutoffs for key residues. The insensitivity to M shown in this work provides a high degree of confidence that any amount of additional data will improve the overall categorization without spoiling the classification already achieved using coarser-level data.