htmd.model module#

Markov state models are a statistical tool for analysing molecular simulations which has met with lots of success. The Model class here, encapsulates all functionallity for the calculation of Markov models while hiding unnecessary details under the hood. It uses DeepTime [1] internally to calculate Markov models.

References

class htmd.model.Model(data: MetricData | None = None, file: str | None = None)#

Bases: object

Constructor for the Model class.

Parameters:

data (MetricData object, optional) – A MetricData object containing the discretized trajectories
file (str, optional) – Path to a previously saved Model file to load. Provide either data or file.

Examples

>>> model = Model(data)

Methods

Attributes

property P: ndarray#: The transition probability matrix

Conducts a Chapman-Kolmogorov test.

The Chapman-Kolmogorov test validates the model estimated at the production lag time tau (the lag passed to markovModel()) by checking the identity T(k*tau) = T(tau)^k. When lags is not given, the test is anchored on the production lag and the test lag times are the integer multiples tau, 2*tau, 3*tau, .... The production-lag model is used as the reference, so the model you actually built is the one being validated.

Parameters:

lags (list or range or np.ndarray, optional) – Explicit list of lag times to test. If given, the smallest lag is used as the reference model, so it should be the production lag.
maxlag (float, optional) – Largest lag time to test. Used if lags is None. By default it will use the mode of the trajectory lengths. The number of multiples is chosen so that k*tau does not exceed maxlag.
numlags (int, optional) – Upper bound on the number of multiples of the production lag to test. Used if lags is None.
units (str, optional) – Units of the lag times. By default it will use frames.
plot (bool) – If the method should display the plot of the CK test
save (str) – Path of the file in which to save the figure
errors (int, optional) – Number of Bayesian samples to use for the error bars. If set to None, it will not plot error bars and use a Maximum Likelihood MSM.
estimator (str) – Which MSM estimator to use for each lag time. 'mle' (default) or 'oom'. 'oom' cannot be combined with errors (Bayesian sampling).

property cluster_ofmicro: ndarray#

Mapping of microstates to clusters

Numpy array which at index i has the index of the cluster corresponding to microstate i.

copy() → Model#

Produces a deep copy of the object

Returns:: data – A copy of the current object
Return type:: Model object

Examples

>>> model2 = model.copy()

createCoreSetModel(threshold: float = 0.5) → tuple#

Given an MSM this function detects the states belonging to a core set and returns a new model consisting only of these states.

Parameters:

threshold (float) – Membership probability threshold over which microstates belong to the core of a macrostate

Returns:

newmodel (Model object) – A new model object
frames (list) – A list of the frames that were kept in the new model

createState(microstates: list | ndarray | None = None, indexpairs: list | None = None)#

Creates a new state. Works both for new clusters and macrostates.

If creating a new cluster, it just reassigns the given frames to the new cluster. If creating a new macrostate, it removes the given microstates from their previous macrostate, creates a new one and assigns them to it.

Parameters:

microstates (list or np.ndarray) – The microstates to split out into a new macrostate.
indexpairs (list) – List of lists. Each row is a simulation index-frame pair which should be added to a new cluster.

eqDistribution(plot: bool = True, save: str | None = None) → ndarray#

Obtain and plot the equilibrium probabilities of each macrostate

Parameters:

plot (bool, optional) – Disable plotting of the equilibrium distribution by setting it to False
save (str) – Path of the file in which to save the figure

Returns:

eq – An array of equilibrium probabilities of the macrostates

Return type:

np.ndarray

Examples

>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> model.eqDistribution()

Get samples of MSM states in Molecule classes

Parameters:

states (list or np.ndarray, optional) – A list of states to visualize
statetype (str, optional) – The type of state to visualize. Can be ‘macro’, ‘micro’ or ‘cluster’.
wrapsel (str or np.ndarray, optional) – Atom selection used for wrapping, as an atom selection string, a boolean mask, or an integer index array (see Molecule.atomselect). Set to None to disable wrapping. See more here
alignsel (str or np.ndarray, optional) – Atom selection used for aligning all frames, as an atom selection string, a boolean mask, or an integer index array. Set to None to disable aligning. See more here
alignmol (Molecule object) – A reference molecule onto which to align all others
samplemode (str, optional) – How to obtain the samples from the states. Can be ‘weighted’ or ‘random’.
numsamples (int) – Number of samples (conformations) for each state.
simlist (numpy.ndarray of Sim objects) – Optionally pass a different (but matching, i.e. filtered) simlist for creating the Molecules.

Returns:

mols – A list of Molecule objects containing the samples of each state

Return type:

list of Molecule objects

Examples

>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> mols = model.getStates()
>>> for m in mols:
>>>     m.view()

load(filename: str)#

Load a MetricData object from disk

Parameters:: filename (str) – Path to the saved MetricData object

Examples

>>> model = Model()
>>> model.load('./model.dat')

property macro_ofcluster: ndarray#

Mapping of clusters to macrostates

Numpy array which at index i has the index of the macrostate corresponding to cluster i. Clusters which were not connected and thus are not in the model have a macrostate value of -1.

property macro_ofmicro: ndarray#

Mapping of microstates to macrostates

Numpy array which at index i has the index of the macrostate corresponding to microstate i.

property macronum: int#: Number of macrostates

markovModel(lag: float, macronum: int, units: str = 'frames', estimator: str = 'mle')#

Build a Markov model at a given lag time and calculate metastable states

Parameters:

lag (float) – The lag time at which to calculate the Markov state model. The units are specified with the units argument.
macronum (int) – The number of macrostates (metastable states) to produce
units (str) – The units of lag. Can be ‘frames’ or any time unit given as a string.
estimator (str) – Which MSM estimator to use. 'mle' (default) builds a reversible maximum-likelihood MSM. 'oom' uses deeptime’s OOM-reweighted estimator, which corrects for the effect of starting out of equilibrium and is suited to short, adaptively sampled trajectories. Note: 'oom' can occasionally yield negative stationary probabilities on ill-conditioned or sparsely sampled data.

Examples

>>> model = Model(data)
>>> model.markovModel(150, 4)  # 150 frames lag, 4 macrostates

maxConnectedLag(lags: list | ndarray, njobs: int = 1)#

Heuristic for getting the lagtime before a timescale drops.

It calculates the last lagtime before a drop occurs in the first implied timescale due to disconnected states. If the top timescale is closer to the second top timescale at the previous lagtime than to itself at the previous lagtime it means that a drop occured. The lagtime before the drop is returned.

Parameters:

lags (list or np.ndarray) – A list of lag times for which to calculate the implied timescales
njobs (int) – Number of parallel jobs (currently unused).

Returns:

ml – The maximum lagtime before a drop occurs in the top timescale

Return type:

int

Examples

>>> model = Model(data)
>>> model.maxConnectedLag(list(range(1, 100, 5)))

property micro_ofcluster: ndarray#

Mapping of clusters to microstates

Numpy array which at index i has the index of the microstate corresponding to cluster i. Clusters which were not connected and thus are not in the model have a microstate value of -1.

property micronum: int#: Number of microstates

plotFES(dimX: int, dimY: int, temperature: float, states: bool = False, s: float = 10, cmap=None, fescmap=None, statescmap=None, plot: bool = True, save: str | None = None, data: MetricData | None = None, levels: int = 7)#

Plots the free energy surface on any given two dimensions. Can also plot positions of states on top.

Parameters:

dimX (int) – Index of projected dimension to use for the X axis.
dimY (int) – Index of projected dimension to use for the Y axis.
temperature (float) – Simulation temperature.
states (bool) – If True, will plot scatter plot of microstates coloured by macro state on top of FES.
s (float) – Marker size for states.
cmap (str or matplotlib colormap) – Sets the Matplotlib colormap for both fescmap and statescmap
fescmap (str or matplotlib colormap) – Matplotlib colormap for the free energy surface
statescmap (str or matplotlib colormap) – Matplotlib colormap for the states
plot (bool) – If the method should display the plot of the FES. If both plot=False and save=None, the method will return the figure and axes.
save (str) – Path of the file in which to save the figure. If both plot=False and save=None, the method will return the figure and axes.
data (MetricData object) – Optionally you can pass a different MetricData object than the one used to build the model. For example if the user wants to build a model on distances but wants to plot the FES on top of RMSD values. The MetricData object needs to have the same simlist as the Model.
levels (int) – Number of contour levels to use for the free energy surface.

Returns:

f (matplotlib.figure.Figure) – The matplotlib figure object
ax (matplotlib.axes.Axes) – The matplotlib axes object

Examples

>>> import matplotlib as plt
>>> model.plotFES(0, 1, 300)
>>> model.plotFES(2, 3, 300, data=otherdata, states=True, fescmap=plt.cm.gray)

Plot the implied timescales of MSMs of various lag times

Parameters:

lags (list or range or np.ndarray) – Specify explicitly at which lag times to compute the timescales.
minlag (float) – The minimum lag time for the timescales. Used in combination with maxlag and numlags.
maxlag (float) – The maximum lag time for the timescales. If None will default to the mode length of the trajectories.
numlags (int) – The number of points to place between minlag and maxlag.
units (str) – The units of lag. Can be ‘frames’ or any time unit given as a string.
errors (int) – Number of Bayesian samples used to calculate errors with Bayesian MSMs. If None, no errors are calculated.
nits (int) – Number of implied timescales to calculate. If None, all are calculated.
results (bool) – If the method should return the calculated implied timescales
plot (bool) – If the method should display the plot of implied timescales
save (str) – Path of the file in which to save the figure
njobs (int) – DEPRECATED: Number of parallel jobs to spawn for calculation of timescales. Negative numbers are used for spawning jobs as many as CPU threads. -1: for all CPUs -2: for all except one etc.
ylog (bool) – Set to False to get linear y axis instead of logarithmic
estimator (str) – Which MSM estimator to use for each lag time. 'mle' (default) or 'oom'. 'oom' cannot be combined with errors (Bayesian sampling).

Returns:

If given results=True this method will return the following data
its (np.ndarray) – The calculated implied timescales. 2D array with dimensions (len(lags), nits)
lags (np.ndarray) – A list of the lag times that were used to calculate the implied timescales

Examples

>>> model = Model(data)
>>> model.plotTimescales()
>>> model.plotTimescales(lags=list(range(1,100,5)))
>>> model.plotTimescales(minlag=0.1, maxlag=20, numlags=25, units='ns')

Samples frames from a set of states

Parameters:

states (list or int) – A list of state indexes from which we want to sample
frames (int or list) – An integer with the number of frames we want to sample per state or a list of same length as states which contains the number of frames we want from each of the states. If set to None it will return all frames of the states.
statetype (str) – The type of state we want to sample from. Can be ‘micro’, ‘macro’ or ‘cluster’.
replacement (bool) – If we want to sample with or without replacement.
samplemode (str) – What sampling strategy to use. Can be ‘random’, ‘even’ or ‘weighted’. For statetype == ‘macro’ this can be set to ‘even’ to sample evenly from all microstates in the macrostate or to ‘weighted’ to sample proportional to the equilibrium probability of each microstate in the macrostate.
allframes (bool) – Deprecated. Use frames=None instead.

Returns:

absframes (numpy.ndarray) – An array which contains for each state an array containing absolute trajectory frames
relframes (numpy.ndarray) – An array which contains for each state a 2D array containing the trajectory ID and frame number for each of the sampled frames

Examples

>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> model.sampleStates(range(5), [10, 3, 2, 50, 1])  # Sample from all 5 macrostates
>>> model.sampleStates(range(model.micronum), samplesnum, statetype='micro')  # Sample from all microstates

save(filename: str)#

Save a Model object to disk

Parameters:: filename (str) – Path of the file in which to save the object

Examples

>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> model.save('./model.dat')

Visualize macro/micro/cluster states in VMD

Parameters:

states (list or np.ndarray, optional) – A list of states to visualize
statetype (str, optional) – The type of state to visualize. Can be ‘macro’, ‘micro’ or ‘cluster’.
protein (bool, optional) – Set to True to enable pure protein system visualization
ligand (str or np.ndarray, optional) – Atom selection for the ligand, as an atom selection string, a boolean mask, or an integer index array (see Molecule.atomselect). See more here
mols (list, optional) – A list of Molecule objects to visualize
numsamples (int) – Number of samples (conformations) for each state.
wrapsel (str or np.ndarray, optional) – Atom selection used for wrapping, as an atom selection string, a boolean mask, or an integer index array (see Molecule.atomselect). Set to None to disable wrapping. See more here
alignsel (str or np.ndarray, optional) – Atom selection used for aligning all frames, as an atom selection string, a boolean mask, or an integer index array. Set to None to disable aligning. See more here
gui (bool) – Set to True to enable the GUI in the NGL viewer.
simlist (numpy.ndarray of Sim objects) – Optionally pass a different (but matching, i.e. filtered) simlist for visualizing the states.

Examples

>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> model.viewStates(protein=True)

>>> model.viewStates(ligand='resname MOL')

htmd.model.getStateStatistic(reference: Model | MetricData, data: MetricData, states: list | range | ~numpy.ndarray, statetype: str = 'macro', weighted: bool = False, method=<function mean>, axis: int | None = 0) → list#

Calculates properties of the states.

Calculates properties of data corresponding to states. Can calculate for example the mean distances of atoms in a state, or the standard deviation of the RMSD in a state.

Parameters:

reference (Model object or MetricData object) – A model containing the state definitions or a MetricData object containing cluster definitions
data (MetricData object) – A projection corresponding to the conformations in the states of the model. The data and the model need to share the same simlist
states (list or np.ndarray) – A list of the states for which we want to calculate the properties
statetype (str) – The state type. Can be ‘macro’, ‘micro’ or ‘cluster’.
weighted (bool) – If the properties of the macrostates should be calculated as the weighted average of their microstates, where the weights are the equilibrium probabilities of the microstates
method – A function pointer for the method that calculates our desired property. e.g. np.mean, np.std
axis (int) – On which axis of the data the method should operate on

Returns:

statistic – A list which contains in each element the desired statistic of a specific state specifies in states

Return type:

list

Examples

>>> data = MetricDistance.project(sims, 'protein and name CA', 'resname MOL')
>>> model = Model(data)
>>> model.markovModel(100, 5)
>>> # Get the standard deviation of distances in all macrostates
>>> getStateStatistic(model, data, list(range(5)), method=np.std)

htmd.model.macroAccumulate(model: Model, microvalue: ndarray) → ndarray#

Accumulate values of macrostates from a microstate array

Parameters:

model (Model object) – The model which to use to accumulate
microvalue (np.ndarray) – An array of values corresponding to the microstates of the model

Returns:

macrovalue – An array of values corresponding to the macrostates of the model

Return type:

np.ndarray