pyemma.msm.ui.MSM¶

class pyemma.msm.ui.MSM(T, dt='1 step')¶

Markov model with a given transition matrix

Parameters:

P (ndarray(n,n)) – transition matrix
dt (str, optional, default='1 step') –
Description of the physical time corresponding to the lag. May be used by analysis algorithms such as plotting tools to pretty-print the axes. By default ‘1 step’, i.e. there is no physical time unit. Specify by a number, whitespace and unit. Permitted units are (* is an arbitrary string):

‘fs’, ‘femtosecond*’

‘ps’, ‘picosecond*’

‘ns’, ‘nanosecond*’

‘us’, ‘microsecond*’

‘ms’, ‘millisecond*’

‘s’, ‘second*’

__init__(T, dt='1 step')¶

Methods

`__init__`(T[, dt])
`committor_backward`(A, B)	Backward committor from set A to set B
`committor_forward`(A, B)	Forward committor (also known as p_fold or splitting probability) from set A to set B
`correlation`(a[, b, maxtime, k, ncv])	Time-correlation for equilibrium experiment.
`eigenvalues`([k, ncv])	Compute the transition matrix eigenvalues
`eigenvectors_left`([k, ncv])	Compute the left transition matrix eigenvectors
`eigenvectors_right`([k, ncv])	Compute the right transition matrix eigenvectors
`expectation`(a)	Equilibrium expectation value of a given observable.
`fingerprint_correlation`(a[, b, k, ncv])	Dynamical fingerprint for equilibrium time-correlation experiment.
`fingerprint_relaxation`(p0, a[, k, ncv])	Dynamical fingerprint for perturbation/relaxation experiment.
`mfpt`(A, B)	Mean first passage times from set A to set B, in units of the input trajectory time step
`pcca`(m)	Runs PCCA++ [1]_ in order to compute a fuzzy metastable decomposition of MSM states
`relaxation`(p0, a[, maxtime, k, ncv])	Simulates a perturbation-relaxation experiment.
`timescales`([k, ncv])	The relaxation timescales corresponding to the eigenvalues

Attributes

`is_reversible`	Returns whether the MSM is reversible
`is_sparse`	Returns whether the MSM is sparse
`metastable_assignments`	Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_
`metastable_distributions`	Computes the probability distributions of active set states within each metastable set using the PCCA++ method [1]_ using Bayesian inversion as described in [2]_.
`metastable_memberships`	Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_.
`metastable_sets`	Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_
`nstates`	The active set of states on which all computations and estimations will be done
`stationary_distribution`	The stationary distribution, estimated on the active set.
`timestep`	Returns the physical time corresponding to one step of the transition matrix as string, e.g.
`transition_matrix`	The transition matrix, estimated on the active set.

committor_backward(A, B)¶

Backward committor from set A to set B

Parameters:	A (int or int array) – set of starting states B (int or int array) – set of target states

committor_forward(A, B)¶

Forward committor (also known as p_fold or splitting probability) from set A to set B

Parameters:	A (int or int array) – set of starting states B (int or int array) – set of target states

correlation(a, b=None, maxtime=None, k=None, ncv=None)¶

Time-correlation for equilibrium experiment.

In order to simulate a time-correlation experiment (e.g. fluorescence correlation spectroscopy [1]_, dynamical neutron scattering [2]_, ...), first compute the mean values of your experimental observable \(a\) by MSM state:

where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\). If a cross-correlation function is wanted, also apply the above computation to a second experimental observable \(b\).

Then the precise (i.e. without statistical error) autocorrelation function of \(f(x_t)\) given the Markov model is computed by correlation(a), and the precise cross-correlation function is computed by correlation(a,b). This is done by evaluating the equation

where \(acf\) stands for autocorrelation function and \(ccf\) stands for cross-correlation function, \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.

Note that instead of using this method you could generate a long synthetic trajectory from the MSM using generate_traj() and then estimating the time-correlation of your observable(s) directly from this trajectory. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely long synthetic trajectory the two results will agree exactly. The correlation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.

Parameters:

a ((M,) ndarray) – Observable, represented as vector on state space
maxtime (int or float) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\) By default (None), the maxtime will be set equal to the 3 times the slowest relaxation time of the MSM, because after this time the signal is constant.
b ((M,) ndarray (optional)) – Second observable, for cross-correlations
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power.
ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k

Returns:

times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the correlation has been computed
correlations (ndarray (N)) – Correlation values at given times

Examples

This example computes the autocorrelation function of a simple observable on a three-state Markov model and plots the result using matplotlib:

>>> import numpy as np
>>> import pyemma.msm as msm
>>>
>>> P = np.array([[0.99, 0.01, 0], [0.01, 0.9, 0.09], [0, 0.1, 0.9]])
>>> a = np.array([0.0, 0.5, 1.0])
>>> M = msm.markov_model(P)
>>> times, acf = M.correlation(a)
>>>
>>> import matplotlib.pylab as plt
>>> plt.plot(times, acf)

References

[1]	Noe, F., S. Doose, I. Daidone, M. Loellmann, J. D. Chodera, M. Sauer and J. C. Smith. 2011 Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. Proc. Natl. Acad. Sci. USA 108, 4822-4827.

[2]	Lindner, B., Z. Yi, J.-H. Prinz, J. C. Smith and F. Noe. 2013. Dynamic Neutron Scattering from Conformational Dynamics I: Theory and Markov models. J. Chem. Phys. 139, 175101.

eigenvalues(k=None, ncv=None)¶

Compute the transition matrix eigenvalues

Parameters:	k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1 ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns:	ts – transition matrix eigenvalues \(\lambda_i, i = 1,...,k\)., sorted by descending norm.
Return type:	ndarray(m)

eigenvectors_left(k=None, ncv=None)¶

Compute the left transition matrix eigenvectors

Parameters:	k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1 ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns:	L – left eigenvectors in a row matrix. l_ij is the j’th component of the i’th left eigenvector
Return type:	ndarray(k,n)

eigenvectors_right(k=None, ncv=None)¶

Compute the right transition matrix eigenvectors

Parameters:	k (int) – number of timescales to be computed. By default identical to the number of eigenvalues computed minus 1 ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns:	R – right eigenvectors in a column matrix. r_ij is the i’th component of the j’th right eigenvector
Return type:	ndarray(n,k)

expectation(a)¶

Equilibrium expectation value of a given observable.

Parameters:	a ((M,) ndarray) – Observable vector
Returns:	val – Equilibrium expectation value fo the given observable
Return type:	float

Notes

The equilibrium expectation value of an observable a is defined as follows

\[\mathbb{E}_{\mu}[a] = \sum_i \mu_i a_i\]

\(\mu=(\mu_i)\) is the stationary vector of the transition matrix \(T\).

fingerprint_correlation(a, b=None, k=None, ncv=None)¶

Dynamical fingerprint for equilibrium time-correlation experiment.

Parameters:

a ((M,) ndarray) – Observable, represented as vector on state space
b ((M,) ndarray (optional)) – Second observable, for cross-correlations
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation. This option is only relevant for sparse matrices and long times for which an eigenvalue decomposition will be done instead of using the matrix power
ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k

Returns:

timescales ((N,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
amplitudes ((N,) ndarray) – Amplitudes for the correlation experiment

References

[1]	Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.

fingerprint_relaxation(p0, a, k=None, ncv=None)¶

Dynamical fingerprint for perturbation/relaxation experiment.

Parameters:

p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
a ((n,) ndarray) – Observable, represented as vector on state space
lag (int or int array) – List of lag time or lag times (in units of the transition matrix lag time tau) at which to compute correlation
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation
ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k

Returns:

timescales ((N,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
amplitudes ((N,) ndarray) – Amplitudes for the relaxation experiment

References

[1]	Noe, F, S Doose, I Daidone, M Loellmann, M Sauer, J D Chodera and J Smith. 2010. Dynamical fingerprints for probing individual relaxation processes in biomolecular dynamics with simulations and kinetic experiments. PNAS 108 (12): 4822-4827.

is_reversible¶: Returns whether the MSM is reversible

is_sparse¶: Returns whether the MSM is sparse

metastable_assignments¶

Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_

pcca() needs to be called first before this attribute is available.

This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!

Returns:
Return type:	For each active set state, the metastable state it is located in.

See also

pcca: to compute the metastable decomposition

References

[1]	Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179

metastable_distributions¶

Computes the probability distributions of active set states within each metastable set using the PCCA++ method [1]_ using Bayesian inversion as described in [2]_.

pcca() needs to be called first before this attribute is available.

Returns:	p_out – A matrix containing the probability distribution of each active set state, given that we are in one of the m metastable sets. i.e. p(state \| metastable). The row sums of p_out are 1.
Return type:	ndarray((m,n))

See also

pcca: to compute the metastable decomposition

References

[1]	Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179

[2]	F. Noe, H. Wu, J.-H. Prinz and N. Plattner: Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114 (2013)

metastable_memberships¶

Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_.

pcca() needs to be called first before this attribute is available.

Returns:	M – A matrix containing the probability or membership of each state to be assigned to each metastable set. i.e. p(metastable \| state). The row sums of M are 1.
Return type:	ndarray((n,m))

See also

pcca: to compute the metastable decomposition

References

[1]	Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179

metastable_sets¶

Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_

pcca() needs to be called first before this attribute is available.

This is only recommended for visualization purposes. You cannot compute any actual quantity of the coarse-grained kinetics without employing the fuzzy memberships!

Returns:
Return type:	A list of length equal to metastable states. Each element is an array with microstate indexes contained in it

See also

pcca: to compute the metastable decomposition

References

[1]	Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179

mfpt(A, B)¶

Mean first passage times from set A to set B, in units of the input trajectory time step

Parameters:	A (int or int array) – set of starting states B (int or int array) – set of target states

nstates¶: The active set of states on which all computations and estimations will be done

pcca(m)¶

Runs PCCA++ [1]_ in order to compute a fuzzy metastable decomposition of MSM states

After calling this method you can access metastable_memberships(), metastable_distributions(), metastable_sets() and metastable_assignments()

Parameters:	m (int) – Number of metastable sets
Returns:	pcca_obj – An object containing all PCCA quantities. However, you can also ingore this return value and instead retrieve the quantities of your interest with the following MSM functions: `metastable_memberships()`, `metastable_distributions()`, `metastable_sets()` and `metastable_assignments()`.
Return type:	`PCCA`

References

[1]	Roeblitz, S and M Weber. 2013. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 7 (2): 147-179

relaxation(p0, a, maxtime=None, k=None, ncv=None)¶

Simulates a perturbation-relaxation experiment.

In perturbation-relaxation experiments such as temperature-jump, pH-jump, pressure jump or rapid mixing experiments, an ensemble of molecules is initially prepared in an off-equilibrium distribution and the expectation value of some experimental observable is then followed over time as the ensemble relaxes towards equilibrium.

In order to simulate such an experiment, first determine the distribution of states at which the experiment is started, \(p_0\) and compute the mean values of your experimental observable \(a\) by MSM state:

where \(S_i\) is the set of configurations belonging to MSM state \(i\) and \(f()\) is a function that computes the experimental observable of interest for configuration \(x_t\).

Then the precise (i.e. without statistical error) time-dependent expectation value of \(f(x_t)\) given the Markov model is computed by relaxation(p0, a). This is done by evaluating the equation

where \(E\) stands for the expectation value that relaxes to its equilibrium value that is identical to expectation(a), \(\mathbf{P(\tau)}\) is the transition matrix at lag time \(\tau\), \(\boldsymbol{\pi}\) is the equilibrium distribution of \(\mathbf{P}\), and \(k\) is the time index.

Note that instead of using this method you could generate many synthetic trajectory from the MSM using generate_traj() that with starting points drawn from the initial distribution and then estimating the time-dependent expectation value by an ensemble average. However, there is no reason to do this because the present method does that calculation without any sampling, and only in the limit of an infinitely many trajectories the two results will agree exactly. The relaxation function computed by the present method still has statistical uncertainty from the fact that the underlying MSM transition matrix has statistical uncertainty when being estimated from data, but there is no additional (and unnecessary) uncertainty due to synthetic trajectory generation.

Parameters:

p0 ((n,) ndarray) – Initial distribution for a relaxation experiment
a ((n,) ndarray) – Observable, represented as vector on state space
maxtime (int or float, optional, default = None) – Maximum time (in units of the input trajectory time step) until which the correlation function will be evaluated. Internally, the correlation function can only be computed in integer multiples of the Markov model lag time, and therefore the actual last time point will be computed at \(\mathrm{ceil}(\mathrm{maxtime} / \tau)\). By default (None), the maxtime will be set equal to the 3 times the slowest relaxation time of the MSM, because after this time the signal is constant.
k (int (optional)) – Number of eigenvalues and eigenvectors to use for computation
ncv (int (optional)) – Only relevant for sparse matrices and large lag times, where the relaxation will be computes using an eigenvalue decomposition. The number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k

Returns:

times (ndarray (N)) – Time points (in units of the input trajectory time step) at which the relaxation has been computed
res (ndarray) – Array of expectation value at given times

stationary_distribution¶

The stationary distribution, estimated on the active set.

For example, for connectivity=’largest’ it will be the transition matrix amongst the largest set of reversibly connected states

timescales(k=None, ncv=None)¶

The relaxation timescales corresponding to the eigenvalues

Parameters:	k (int) – number of timescales to be computed. As a result, k+1 eigenvalues will be computed ncv (int (optional)) – Relevant for eigenvalue decomposition of reversible transition matrices. ncv is the number of Lanczos vectors generated, ncv must be greater than k; it is recommended that ncv > 2*k
Returns:	ts – relaxation timescales in units of the input trajectory time step, defined by \(-tau / ln \| \lambda_i \|, i = 2,...,k+1\).
Return type:	ndarray(m)

timestep¶: Returns the physical time corresponding to one step of the transition matrix as string, e.g. ‘10 ps’

transition_matrix¶: The transition matrix, estimated on the active set. For example, for connectivity=’largest’ it will be the transition matrix amongst the largest set of reversibly connected states