Aromatic side-chain flipping acts a probe to understand conformational dynamics of biomolecules. In this post, we will see how to predict transition kinetics associated with different side-chain rotameric conformations using Markov State Modelling (MSM). MSM allows one to predict kinetics and equilibrium population associated with conformational dynamics. In this post we will use BACE1 as a model system. BACE1 is a target for Alzheimer’s disease. Conserved tyrosine present in flap of BACE1 governs conformational dynamics as highlighted by Bhakat and Soderhjelm (Bhakat, S.; Soderhjelm, P. Flap Dynamics in Pepsin-Like Aspartic Proteases: A Computational Perspective Using Plasmepsin-II and BACE-1 as Model Systems. Journal of Chemical Information and Modeling 2022, 62, 914-926).
We have generated ~16 microsecond molecular dynamics simulation to sample N to F transition and vice-versa.
In order to sample transition we will first featurize our molecular dynamics data onto sin/cos transformed χ1 and χ2 angles of Tyrosine and use the features to do clustering and generate MSM. The following script does just that.
#msmbuilder imports from msmbuilder.featurizer import DihedralFeaturizer from msmbuilder.utils import verbosedump,verboseload from msmbuilder.utils import load,dump #other imports import os,glob,shutil import numpy as np import mdtraj as md import pandas as pd import pickle #prettier plots #Loading the trajectory ref = md.load('prot.pdb') # Zero indexed. This selectes Tyr73 from PDB a = ref.top.select("resid 72") # Uploading PBC fixed XTC file ds = dataset("*.xtc", topology="prot.pdb", atom_indices=a, stride=20) #Featurization, sin/cos transformed which uses Chi1 and Chi2 angle of Tyr73 featurizer = DihedralFeaturizer(types=['chi1', 'chi2']) dump(featurizer,"transformed_raw_featurizer.pkl") # Sin/cos untransformed featurizer f=DihedralFeaturizer(types=['chi1', 'chi2'], sincos=False) dump(f,"raw_featurizer.pkl") # Fitting the featurizer to the .xtc file diheds = featurizer.fit_transform(ds) dump(diheds, "features.pkl") dihedsraw = f.fit_transform(ds) dump(dihedsraw, "raw-features.pkl") print(diheds.shape) #Robust scaling, scaling featurization from msmbuilder.preprocessing import RobustScaler scaler = RobustScaler() scaled_diheds = scaler.fit_transform(diheds) dump(scaled_diheds, "scaled-transformed-features.pkl")
Next we will use features.pkl to do clustering using micro state K-means algorithm with K=500 (see the notebook here: https://github.com/sbhakat/MSM-ring-flipping/blob/main/pcca-msm.ipynb)
We will then perform an implied timescale test. We can clearly see from Figure 2 that the a lag time 200 is a good choice for MSM.
MSM with lag time 200 allowed us to predict equilibrium weighted free energy surface along χ1 and χ2 angles of Tyrosine (Figure 3).
Finally we will generate coarse grained representation of MSM using PCCA. PCCA managed to separate different rotational states of Tyr as shown in Figure 4.
Finally we will calculate kinetics of ring flipping using msm.mfpt() module of PyEMMA.
Here I show how to predict kinetics associated with ring flipping which is fast flipping event (< 1 ms). One can use NMR 13C relaxation dispersion experiment to validate our MSM predicted values.