Laplacian score and genetic algorithm based automatic feature selection for Markov State Models in adaptive sampling based molecular dynamics

Markov State Models (MSM)

MSMs are used to model randomly changing systems. They have been used to study time-series events in which the current state depends only on the previous state and not on any other states before it (Tang, Bevan & Grover, 2017). It has a wide variety of applications in the field of physics, chemistry, medicine, finance, management etc. It has been used for biological modelling such as the simulation of brain function (George & Hawkins, 2009), viral infection of single cells (Gupta & Rawlings, 2014), analysis of bacterial genome (Skewes & Welch, 2013), as well as MD simulations using adaptive sampling techniques (Doerr et al., 2016; Doerr & De Fabritiis, 2014) etc.

Figure 1 represents the workflow adopted towards the application of MSM to analyze MD trajectories of bio-molecules. The goal is to find meta-stable conformers, with a few variations as explained below. Prediction accuracy of MSM, as indicative by GMRQ score, can be improved by selecting optimal features from the simulation trajectory data.

MSM is represented by an Nx N matrix which gives information about the probability of transition between N states (Singhal, Snow & Pande, 2004; Noé & Fischer, 2008). In Fig. 2, x_ij is the average of transition count matrix and its transpose. The probability of each row in the average transition count matrix is calculated to obtain the MSM. From the MSM, the eigen-flux is calculated. Eigen flux (-n1 for s4, n2 for s3, n3 for s2 and n4 for s1) of each state shows the transition from source state to the sink state, and this is used to identify the slowest process (Beauchamp et al., 2011). The N states in the MSM represents the state decomposition into which the simulation data is clustered. The ideal MSM, or in other words, the ideal state decomposition reveals even the slowest dynamical process. In this specific use case, it is the conformational change occurring in the simulated system. The probability with which transition to a state j occurs is calculated as (1)${\displaystyle p_j\left(t+\tau \right)=\sum _{i=1}^N{p}_i\left(t\right)T_{ij}\left(\tau \right)}$ It can be generalized as (2)${\displaystyle p\left(t+\tau \right)=\sum _{i=1}^{N}p\left(t\right)T\left(\tau \right)}$ where p_j(t + τ) is a population vector after time τ; τ is lag time; p_i(t) is a population vector at time t; T_ij(τ) is the probability of transition from state i to state j; T(τ) is the probability transition matrix that represents the MSM.

The eigenvalues λ_i and eigenfunctions ψ_i of T(τ) can decompose the transition probability matrix (Shukla et al., 2015) as in Eq. (3). The eigenfunctions represents the slowest processes that occurs in the simulated system. (3)${\displaystyle T\left(\tau \right)⊛{\psi }_i={\lambda }_i{\psi }_i}$

Generalized Matrix Rayleigh Quotient (GMRQ) derived from the variational principle of conformation dynamics (Noé & Nuske, 2013) is a metric used to evaluate the efficiency of the MSM (McGibbon & Pande, 2015). MSMs with higher GMRQ scores provide higher discriminatory abilities towards identifying the slowest dynamical process occurring in the system (McGibbon & Pande, 2015). Variational principle is used in the calculus of variations, or changes to develop a function that maximizes the variables dependent on it. Based on variational principle, the eigen values and eigen functions are estimated for the MSM constructed using the simulated data. GMRQ score of the MSM is the sum of the first p eigenvalues of T(τ). They also identify the slow processes, p. The upper boundary for the GMRQ score is set by variational principle as shown below (Husic & Pande, 2017; Noé & Nuske, 2013): (4)${\displaystyle GMRQ\equiv \sum _{i=1}^p{\hat{\lambda }}_i\le \sum _{i=1}^p{\lambda }_i}$ where ${\hat{\lambda }}_i$; λ_i are estimated and real i-th eigenvalues respectively. The efficiency of MSM is improved by trying to reach the upper boundary set by the sum of real eigenvalues (λ_i). MSMs with higher GMRQ scores are better able to identify the slowest conformational change in the system (Noé & Nuske, 2013; McGibbon & Pande, 2015).

Adaptive sampling MD based on MSM and Optimized feature selection for MSM

We have used simulation trajectories that are produced by adaptive sampling MD, based on the MSM. In this protocol, the first step in the construction of MSM involves the extraction of feature, namely, dihedral angles, distance between amino-acid residues, root mean square deviation (RMSD) from the raw Cartesian coordination, etc. This step is called featurization. The selection of optimal features or collective variables have a critical impact on the thermodynamic and kinetic property calculations predicted by the model, and is explored in many studies (Shamsi, Cheng & Shukla, 2018; Mittal & Shukla, 2017; Sultan et al., 2014). Consequently, using the optimal feature set for the construction of MSM affects its prediction accuracy. Identification of the most critical dihedral or contact features assist in identifying the optimum conformers. These conformers are then used as starting structures for the next epoch of adaptive sampling in MSM based adaptive sampling MD. Appropriate feature selection circumvents issues due to energy barriers leading to convergence of thermodynamic properties (such as free energy of binding in the case of protein-ligand binding simulation), or folding of the unfolded protein (Doerr & De Fabritiis, 2014).

Feature selection is an unsupervised machine learning task. One method of categorizing them identifies three sub-classes, namely, filters (Devakumari & Thangavel, 2010; Mitra, Murthy & Pal, 2002; Tabakhi & Moradi, 2015), wrappers (Dy & Brodley, 2004; Dutta, Dutta & Sil, 2014; Breaban & Luchian, 2011), and hybrids (Solorio-Fernández, Carrasco-Ochoa & Martínez-Trinidad, 2016; Li, Lu & Wu, 2006). In filter methods, the importance of a feature is studied based on the intrinsic properties of the data. They are not selected by training on the model on which it is to be used. In wrapper-based methods, the features are first used to train the model. Then their impact on the performance of the model is evaluated. Filter methods are computationally less expensive than wrapper methods. Hybrid methods combine the advantages of both filter and wrapper methods (Li et al., 2008). In this method, the feature subset works best when the same model used for feature selection is used subsequently (Solorio-Fernández, Carrasco-Ochoa & Martínez-Trinidad, 2020). Laplacian-score based feature selection is a filter-based method. This method finds the features that have the power to preserve the clusters in the data. Features that can maintain the structure of the plotted nearest neighbor graph are selected.

Laplacian scores have been used to reveal hidden structures in the complex conformational space of the intrinsically disordered peptide, Aβ(1-42) (Sgourakis et al., 2011). Laplacian scores identified crucial interactions in these conformers, which provided new drug design hypothesis that can, potentially, be used to discover drugs that inhibit the oligomers and fibril formation critical to the progression of Alzheimer’s disease. In this study, inspired by the success of the study mentioned above, the initial population of a GA’s chromosomes have genes pre-filtered through Laplacian scores-based ordering. In bio-molecules, conformational changes happen in localized regions (Fan et al., 2015). Laplacian score ranks the features based on its ability to preserve the local structure of the graph constructed based on the k-nearest neighbor algorithm (He, Cai & Niyogi, 2006). Due to this property of the Laplacian score, it is used to identify a subset of the features. The feature subset selection is further optimized through the GA.

Datasets to evaluate the approach

Molecular dynamics trajectories of five systems are used to validate our approach by analyzing the kinetic and thermodynamic properties predicted by the MSMs. The features selected were derived from the GA run, whose initial population was based on our Laplacian score ranking scheme. The simulation dataset of the two protein-ligand complexes are Piperidine-Thrombin Doerr et al. (2016) and Benzamidine-Trypsin Scherer et al. (2015). The protein folding trajectories of Fs Peptide (Beauchamp et al., 2011), WW domain (Lindorff-Larsen et al., 2011) (generously shared by D.E Shaw Research), and Villin simulation Doerr et al. (2016) were also analyzed as a part of this study.

Our approach

Figure 1 provides a detailed view of the protocol we have implemented, while the principal choices for each component of the workflow are shown in Table 1. The hyperparameters for the protocol were chosen based on the published research, HTMD and Msmbuilder tutorial (Chen et al., 2018; Doerr et al., 2016; Beauchamp et al., 2011). Scikit-learn (Pedregosa et al., 2011) was used for feature preprocessing, and to calculate the Laplacian scores, IPython (Pérez & Granger, 2007) for the execution of the python scripts in the form of a pipeline. MDTraj (McGibbon et al., 2015) was used to analyze the simulation trajectories, and HTMD (Doerr et al., 2016) was used for the construction of free energy heatmaps, standard free energy plots and cktest graphs for the MSMs. These graphs were constructed for selected residues to implement our approach, and for all residues to re-implement previous work (for comparison). MSMBuilder (Beauchamp et al., 2011) was used for MSM construction and GMRQ scoring.

• Featurization: Featurization is performed to transform rotational and transitional motion from the simulation data to a vector. Dihedral angles and distance between CA atoms of amino-acid residues, from the Cartesian coordinates were extracted as features from the simulation data (for the 5 bio-systems). Dihedral featurization comprises backbone and side chain dihedral angles, along with sine and cosine of these angles (Beauchamp et al., 2011). All of the above are calculated for each frame in the simulation trajectory. In our approach, we used the sine and cosine of back bone dihedral angles for dihedral features. It has been proved to be one of the important metrics that is able to identify different conformers in an MD trajectory (Cossio, Laio & Pietrucci, 2011). Contact featurization calculates the distance between each pair of amino acid residue, or the distance between the specified pair of residues for each frame (Beauchamp et al., 2011). In this study, we implemented contact featurization as a vector of distances between CA atom of each amino acid residue and ligand for the protein-ligand systems, and the distance between the CA atoms of amino acid residues for the protein-folding simulation data.

• Laplacian scoring: An average, across all trajectories, of the Laplacian scores obtained for each feature is calculated. The features were then sorted based on this average Laplacian score (of each feature). This ordered list is used to build the initial population for the GA. Laplacian Score (LS) is fundamentally based on Laplacian Eigenmaps and Locality Preserving Projection. The algorithm for calculating the Laplacians (He, Cai & Niyogi, 2006) is given below. Let L_r denote the Laplacian score of the rth feature. Let f_ri denote the ith sample of the rth feature; where i = 1, …, m.

  1. Construct a nearest neighbor graph—G with m nodes—. The ith node corresponds to x_i. We connect the edge between nodes i and j, if x_i and x_j are close, i.e., x_i is among k nearest neighbors of x_j, or x_j is among k nearest neighbors of x_i.

  2. The weight matrix S of the graph, models the local structure of the data space. ${\text{S}}_{ij}={\text{e}}^{\frac{{\left(x_i{x}_j\right)}^2}{t}}$, when nodes i and j are connected, and S_ij = 0 otherwise.

  3. For the r^th feature, we define: f_r = [f_r1, f_r2, …, f_rm]^T, D = diag(S1), 1 = [1, …, 1]^T, L = DS where the matrix L is called graph Laplacian (Chung, 1997). Let (5)${\displaystyle \widetilde{f_r}=f_r-\frac{{f_r}^{T}D\mathit{1}}{{\mathit{1}}^{T}D\mathit{1}}\mathit{1}}$

  4. The Laplacian score of the r^th feature is: (6)${\displaystyle L_r=\frac{f_r^{-T}L\widetilde{f_r}}{f_r^{-T}D\widetilde{f_r}}}$

• Genetic Algorithm (GA):

  • Initial population - Features are shortlisted based on Laplacian score. To avoid over-fitting or under-fitting, we use the thumb rule of selecting $\frac{1}{10}$th of the number of samples (Sánchez & García, 2018) as the target number of features. This subset is then used as the initial population for GA. Under-fitting occurs when the number of features that represent the model are significantly lesser than the number of data points available. In MD centered protein systems, the number of features or the dimensionality is high due to the nature of the problem, and thus probability of under-fitting is negligible.

  • Fitness calculation - The fitness of the feature subset is scored based on GMRQ score of the respective MSM. Dimensionality reduction is performed using TICA (Pérez-Hernández et al., 2013). Dimensionality reduction produces a linear combination of the features, and the top components (called time structure-based independent components or tICs) capture the most prominent processes in the simulation. Further, clustering of frames is performed using MiniBatchKmeans clustering (McGibbon & Pande, 2015). The MSM is constructed using the clusters identified and the transition probability between these clusters. The number of clusters, lag time of TICA and MSM, number of components used for TICA and n_timescales which represents the timescale of the n slow processes is presented in Table 1. K-fold cross-validation (K = 5) is used to avoid over-fitting. The value of GMRQ score for the MSM is calculated based on the k-fold cross validation.

  • Selection - Selection of parents is performed for the production of off-springs for the next generation. This is performed based on the fitness of the feature subset. This helps in the selection of ideal traits, which in this case is the important CVs.

  • Crossover - The crossover probability decides the number of parents taking part and, in effect, the number of off-springs generated. We have used a single point crossover.

  • Mutation - It is performed to increase the gene variations in each of the generations. Single point mutation is carried out, in which a randomly selected gene is replaced with the highest ranked feature based on its Laplacian score, which is not already present in the chromosome.

  • Repeat the above four steps for the specified number of generations

• MSM construction and GMRQ scoring: The MSM is constructed using the most optimal CVs obtained after the GA. This is then compared with MSM constructed using all the features in terms of the GMRQ score, and also compared using the implied timescale. The methodology is summarized in the Fig. 1

Details on Bio-systems studied

Five systems were analyzed in this work (The original complexes have been uploaded in figshare).

• Benzamidine-Trypsin (Doerr et al., 2016) : This work studies the binding of serine protease beta-trypsin to inhibitor benzamidine. It comprises 10 microseconds of simulation data, with component trajectories of 200 nanoseconds each.

• Piperidine-Thrombin (Scherer et al., 2015) : Thrombin is a serine protease that acts in the coagulation pathway to convert factor XI. The simulation dataset of Piperidine-Thrombin comprises 810 trajectories of 200 nanoseconds each, resulting in a cumulative simulation length of 162 microseconds.

• Fs Peptide (Beauchamp et al., 2011) : Fs peptide is a system widely used to study intricacies of protein folding. The simulations were carried out using the AMBER99SB-ILDN force field with OpenMM 6.0.1. The Fs peptide dataset includes 28 trajectories of 500 nanoseconds each, totalling to 14 microseconds of data.

• Villin(Doerr et al., 2016) : Villin is an actin binding protein. The Villin dataset consists of 1,374 trajectories of 100 nanoseconds each and total simulation data of length 137.4 microseconds.

• WW domain (Lindorff-Larsen et al., 2011) : WW domain is a protein domain that plays an important role in the interaction between protein ligands. The WW-domain dataset consists of 325 trajectories with a total of 1,137 microseconds.

Comparison of GMRQ scores

In this section, we have compared the GMRQ scores of MSM models obtained when the full set of dihedral angle-based features and the full set of features calculated using all contacts are considered, against the ones where we use only selected features using our approach. The summary of these results are shown in Table 2. For the Benzamidine-Trypsin system, the model with select features has a higher GMRQ score of 3.21 and 2.59. The GMRQ scores of MSMs constructed for Villin and Fs Peptide systems with reduced set of dihedral and contact features are both higher by 19.15%, 19.95%, 17.44%, 11.25% respectively, as compared to when constructed with all the features. This indicates that the features identified by our method is able to capture the most important features and reduce the noise by avoiding unwanted features in certain cases. However, in the case of Piperidine-Thrombin and WW-domain systems, the GMRQ score of the MSM constructed with reduced set of selected features is less than the score obtained with all the features selected. This may be due the fact that the interactions captured by the reduced set of features is not enough to capture the slowest process that occurs in the respective simulation data.

The Table 3 shows that our approach is able to mitigate model inaccuracies caused by the over-fitting that occurs when the length of simulation (and hence, number of observations) is relatively shorter, while the number of features are higher. The variance in the GMRQ score for the model with selected residues is 5 to 10 times lesser than when all residues are used for MSM construction. This shows that the MSM constructed using selected residues is able to generalize better, and overcome over-fitting issues. The WW-domain is not given in the Table 3 as simulation length of 1,137 microseconds is adequate, and hence, over-fitting does not occur.

Some studies use the Chapman–Kolmogorov property (Noé et al., 2009; Prinz et al., 2011; Bowman, Pande & Noé, 2013) as a measure of self-consistency of individual MSMs. Nonetheless, since our goal was to find the optimum MSM that identifies the slowest implied timescale, we have used GMRQ scores for the primary analysis. The analysis of Chapman–Kolmogorov test is given for the three systems that produced better results (see Table 2 for the ranking) through our approach is provided in form of figures in the (Figs. S1 through S3).

Standard free energy prediction

The free energy surface is indicative of the thermodynamic and kinetic properties of the system. The free energy surface and the standard free energy plots of the meta-stable states, identified as macro-states in MSMs, of the Villin system are given in Fig. 4. One of the main advantages of an MSM is its ability to predict the thermodynamic and kinetic properties of the system. The meta-stable states identified by the MSM model map to different conformers of the bio-molecule in the simulation. These are identified as different points on the free energy surface. The stable, lowest energy state is the bonded state in protein-ligand binding simulation, and folded state in protein folding simulation. For the sake of brevity, the same for the other two high scoring MSM models, namely, Benzamidine-Trypsin system, Fs Peptide and Villin (dihedral features only) are given in Figs. S4 to S6.

The free energy surfaces of the Villin system for full-set and selected of features are given in Figs. 4A and 4C respectively. The standard free energy of each of the meta-stable states identified in the corresponding free energy heat map are shown for full-set and selected of features in Figs. 4A and 4D. Figure 4A shows four meta-stable states on the free energy surface that was plotted using all the contact features. The standard free energy of each of the meta-stable states identified in Fig. 4A) is shown in Fig. 4B). The lowest energy of −2.29 kcal/mol, observed for the meta-stable state 2, is identified as the sink state by the MSM model. Similarly, four meta-stable states are identified (Fig. 4C) for the MSM model that used contact features between selected residues that were identified by our approach. The standard free energy of each of the meta-stable states is plotted in Fig. 4D, and lowest energy of −2.31 kcal/mol is observed for the meta-stable state 2, identified as the sink state by the MSM model. The lower standard energy stable state (−2.31 kcal/mol) is identified by the MSM model constructed using selected residues identified by our approach.

Implied timescale comparisons

Implied timescale values refer to the time taken by the slowest processes captured by MSMs. Table 4 shows that the MSMs constructed using the selected dihedral and contact features of Benzamidine-Trypsin, Villin and Fs Peptide are able to capture slowest processes. Since our goal is to capture the slowest timescale that represents the folded or the protein-ligand bonded state, the other four slower timescales is given in the Fig. S7, for the interested reader. Five timescales are compared for systems with better MSMs as identified by our approach, namely, Benzamidine-Trypsin, Fs Peptide and Villin.

The MSM constructed using the full-set of all the contacts with the ligand, in the Piperidine-Thrombin system, is able to capture the slower process. Nonetheless, it can be seen that, in this case, GMRQ score for the model from the reduced set is not significantly higher against the one with the full-set of contact-based features.