The MSM analysis workflow#

HTMD’s analysis side is built around a fixed five-step pipeline. Each step has a dedicated class, and the output of one step is the input of the next. Knowing the shape of the pipeline makes the API obvious in retrospect.

        flowchart LR
    A[SimList] --> B["Projection<br/>(MetricDistance, MetricRmsd, ...)"]
    B --> C["Dimensionality reduction<br/>(TICA / GWPCA)"]
    C --> D["Clustering<br/>(MiniBatchKMeans, KCenters, RegCluster)"]
    D --> E["Markov state model<br/>(Model)"]
    E --> F["Kinetics<br/>(rates, MFPTs)"]

1. SimList: enumerate trajectories#

simlist() takes an explicit list of trajectory directories and a matching list of topology files (or one shared topology), and returns a list of Sim objects, each pointing at one trajectory together with its topology. It’s the moral equivalent of glob plus topology matching (you provide the globs), and it’s the only object you carry through the rest of the pipeline.

from htmd.ui import *
sims = simlist(glob("data/*/"), glob("input/*/structure.pdb"))

2. Projection: per-frame features#

A projection maps each frame of each simulation to a feature vector. Projections come from moleculekit’s metric* family (MetricDistance, MetricRmsd, MetricDihedral, MetricSecondaryStructure, and so on), configured once and applied across the whole SimList:

metr = Metric(sims)
metr.set(MetricDistance("name CA", "resname BEN", periodic="selections"))
data = metr.project()

The result is a MetricData object: a ragged list of per-trajectory (n_frames_i, n_features) arrays (different sims can have different lengths) plus the mapping back to individual frames.

3. Dimensionality reduction: keep what matters#

For all but the simplest systems, the raw projection has too many features for clustering to behave well. HTMD provides:

  • TICA: time-lagged independent component analysis, the default choice for MD.

  • GWPCA: globally-weighted PCA.

Both classes take a MetricData and a lag time (in frames or absolute time units) and return a new MetricData in the reduced space. TICA can also wrap a Metric directly for on-the-fly projection without materialising the full feature array first.

4. Clustering: discretise the state space#

Clustering assigns each frame to one of K microstates. HTMD wraps several clustering schemes:

  • MiniBatchKMeans from scikit-learn (the usual default).

  • KCenter: good for low-population states.

  • RegCluster: radius-based assignment around precomputed centres.

The result is a per-frame integer label stored on the same MetricData object.

5. Markov state model and kinetics#

Model builds the Markov state model from the clustered data: counts transitions at a chosen lag time, computes the transition matrix, identifies macrostates by PCCA+, and exposes equilibrium populations and timescales.

Kinetics consumes the Model and computes mean first-passage times, rate constants, and free energies between user-specified source and sink states.

data_clust.fstep = 0.1                          # ns per frame; Kinetics needs this
model = Model(data_clust)
model.markovModel(lag=20, macronum=4)
kin = Kinetics(model, temperature=300, concentration=1e-3)
kin.getRates()

Kinetics raises if model.data.fstep is None; it needs the per-frame timestep to convert rates into wall-clock units.

What this means in practice#

Each tutorial in Analysing trajectories is one trip through this pipeline on a different system. If a tutorial mentions a class you don’t recognise, find the box it sits in above and you’ll know what its inputs and outputs look like.