choderalab/msms

Markov state models

MSMs

Markov state models

Current Models: Src https://github.com/choderalab/MSMs/blob/master/src/README.md

msms's Issues

Plotting other structures onto MSM

I guess we need to check this a bit more thoroughly in general re: ensembler models, but currently have a basic question:

If I plot the first frame of run0-clone0.h5 is that equivalent to plotting the starting structure?

Import of matplotlib.pyplot and numpy frequently missing in these

Maybe there is a reason I am overlooking...

Use contact map features?

John relayed Frank's suggestion to use a subset of inter-residue contacts as input features for MSM-building: we would like to use just "the subset of contacts that change during a simulation." I've been looking into this a bit for modeling the data from 10468 (Abl).

tl;dr -- promising!

My biggest question is:

How should we choose the contact threshold?

To determine whether a given inter-residue contact "changes," we binarize the inter-residue distance as either > or <= a threshold distance, then select any residue pairs that aren't all 0s or all 1s throughout our simulations. PyEMMA has a default setting for this of 5 angstroms. However, the subset of contacts that will appear to be "changing" clearly depends on what threshold you use:

Note that the set of selected contacts at a threshold x+epsilon is not simply a superset of the contacts selected at a lower threshold x:

Here I've just set the threshold=1.0nm since that induces an acceptably small number of features (~4000) that seem to represent the extent of the protein, but I don't know how to do this in a principled way.

Comparing feature sets:

A feature set that does a good job at linearly representing the slow relaxation processes of a system will have a high GMRQ, i.e. a large sum of leading eigenvalues in the tICA solution.

I compared 3 feature sets in terms of their GMRQ for fixed-rank approximation (m=100) of the lag_time=100 propagator. I also computed the maximum cumulative sum of leading eigenvalues.

Feature set 1: Backbone torsions (all available: ‘phi’, ‘psi’, ‘omega’, ‘chi1’, ‘chi2’, ‘chi3’, ‘chi4’)

The sum of the first 100 eigenvalues = 55.7.

The sum of the first 542 eigenvalues = 76.3.

Feature set 2: Inter-residue distances between contacts that change (threshold=1.0nm)
(evaluated on a subset of 10 trajectories)
The sum of the first 100 eigenvalues = 61.9.

The sum of the first 1999 eigenvalues = 92.2.

Feature set 3: Kabsch-Sander hydrogen-bond energies
The sum of the first 100 eigenvalues = 40.4.

The sum of the first 468 eigenvalues = 48.5.

Visualizing the results

From the eigenvalue plots above we can see that the first 2 tICs represent only a tiny fraction of the kinetic behavior of the kinase, for any of these choices of input features-- projecting onto the first 2 tICs is a very lossy representation of the tICA model. Can we create a more trustworthy 2D picture?

Here I'm using the first 100 tICs as inputs into t-SNE, which produces a nonlinear 2D embedding that preserves local neighborhoods at the cost of introducing global distortions. This is useful if we want a picture of the cluster structure of a dataset, but means we can't read into the x and y axes as "components." Each dot represents a simulation frame, dots are colored by the order in which they were collected, and time-adjacent frames are connected by a thin grey line. Areas that are visited more than once should have a mix of colors in them. Areas that are color-homogenuous appear to have been visited during just one contiguous chunk of trajectory.

Backbone torsions

Contact maps

Note that the above pictures were computed using just 10 out of the ~500 trajectories available.

If we recompute the "Contact map --> top-100 tICs --> t-SNE" projection for all available simulation data, then we get the following picture:

This picture suggests that the majority of the available data samples transitions between large relevant states (the dense multi-colored regions in the map), and there may be a large number of negligibly occupied microstates that are visited ~once (the little "sprinkles" on the periphery). Caveat: It's possible that t-SNE is over-partitioning the dataset (i.e. hallucinating clusters)-- although I ran a quick sanity check that the k-nearest neighbor graph (for k=30) is mostly identical in the 2D embedding and in the 100D tICA model, there are some additional sanity checks I can run before we interpret these pictures too much.

Other comments

MDTraj offers several schemes for computing inter-residue distances. I've been using the scheme='ca' (which returns the distance between their alpha carbons), since it's significantly cheaper than scheme='closest-heavy' (which I presume is still cheaper than scheme='closest'), but I don't know which of the available schemes is preferred here.

Test Robustness of tICA with atomPairs

So I did simple experiment where I grabbed random selections of atom pairs, then calculated the tica eigenvectors. They're pretty robust even with a small subset of atom pairs. This suggests that we could probably do something very affordable here and get very reproducible results.

Here is the driver code:

import itertools
import numpy as np
import mdtraj as md
import msmbuilder as msmb

trj0 = md.load("./system.subset.pdb")
trajectories = [md.load("./Trajectories/trj%d.h5" % i) for i in range(5)]

atom_indices = np.arange(trj0.n_atoms)
pair_indices = np.array(list(itertools.combinations(atom_indices, 2)))

n_choose = 1000
n_trials = 10
for i in range(n_trials):
    pair_indices_subset = pair_indices[np.random.choice(len(pair_indices), n_choose, replace=False)]
    metric = msmb.metrics.AtomPairs(atom_pairs=pair_indices_subset)
    tica = msmb.reduce.tICA(1, prep_metric=metric)
    for trj in trajectories:
        tica.train(trajectory=trj)
    tica.solve()
    print tica.vals[0:5]