pyemma.msm.MSM¶
-
class
pyemma.msm.
MSM
(P, pi=None, reversible=None, dt_model='1 step', neig=None, ncv=None)¶ Markov model with a given transition matrix
-
__init__
(P, pi=None, reversible=None, dt_model='1 step', neig=None, ncv=None)¶ Markov model with a given transition matrix
Parameters: - P (ndarray(n,n)) – transition matrix
- pi (ndarray(n), optional, default=None) – stationary distribution. Can be optionally given in case if it was already computed, e.g. by the estimator.
- reversible (bool, optional, default=None) – whether P is reversible with respect to its stationary distribution. If None (default), will be determined from P
- dt_model (str, optional, default='1 step') –
Description of the physical time corresponding to one time step of the MSM (aka lag time). 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*’
- neig (int or None) – The number of eigenvalues / eigenvectors to be kept. If set to None, defaults will be used. For a dense MSM the default is all eigenvalues. For a sparse MSM the default is 10.
- 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.
Methods
__init__
(P[, pi, reversible, dt_model, ...])Markov model with a given transition matrix 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])Compute the transition matrix eigenvalues eigenvectors_left
([k])Compute the left transition matrix eigenvectors eigenvectors_right
([k])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. get_model_params
([deep])Get parameters for this model. 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]_ to compute a metastable decomposition of MSM states propagate
(p0, k)Propagates the initial distribution p0 k times relaxation
(p0, a[, maxtime, k, ncv])Simulates a perturbation-relaxation experiment. set_model_params
(P[, pi, reversible, ...])Call to set all basic model parameters. simulate
(N[, start, stop, dt])Generates a realization of the Markov Model timescales
([k])The relaxation timescales corresponding to the eigenvalues update_model_params
(**params)Update given model parameter if they are set to specific values Attributes
P
The transition matrix on the active set. dt_model
Description of the physical time corresponding to the lag. is_reversible
Returns whether the MSM is reversible is_sparse
Returns whether the MSM is sparse metastable_assignments
Assignment of states to metastable sets using PCCA++ metastable_distributions
Probability of metastable states to visit an MSM state by PCCA++ metastable_memberships
Probabilities of MSM states to belong to a metastable state by PCCA++ metastable_sets
Metastable sets using PCCA++ neig
number of eigenvalues to compute. nstates
Number of active states on which all computations and estimations are done pi
The stationary distribution on the MSM states reversible
Returns whether the MSM is reversible sparse
Returns whether the MSM is sparse stationary_distribution
The stationary distribution on the MSM states timestep_model
Physical time corresponding to one transition matrix step, e.g. transition_matrix
The transition matrix on the active set. -
P
¶ The transition matrix 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 [NDD11], dynamical neutron scattering [LYP13], ...), first compute the mean values of your experimental observable \(a\) by MSM state:
\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]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 accurate (i.e. without statistical error) autocorrelation function of \(f(x_t)\) given the Markov model is computed by correlation(a), and the accurate cross-correlation function is computed by correlation(a,b). This is done by evaluating the equation
\[\begin{split}acf_a(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{a} \\ ccf_{a,b}(k\tau) &= \mathbf{a}^\top \mathrm{diag}(\boldsymbol{\pi}) \mathbf{P(\tau)}^k \mathbf{b}\end{split}\]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 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 ((n,) 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 5 times the slowest relaxation time of the MSM, because after this time the signal is almost constant.
- b ((n,) 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 computed 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
[NDD11] 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. [LYP13] 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.
-
dt_model
¶ Description of the physical time corresponding to the lag.
-
eigenvalues
(k=None)¶ Compute the transition matrix eigenvalues
Parameters: k (int) – number of eigenvalues to be returned. By default will return all available eigenvalues Returns: ts – transition matrix eigenvalues \(\lambda_i, i = 1, ..., k\)., sorted by descending norm. Return type: ndarray(k,)
-
eigenvectors_left
(k=None)¶ Compute the left transition matrix eigenvectors
Parameters: k (int) – number of eigenvectors to be returned. By default all available eigenvectors. 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)¶ Compute the right transition matrix eigenvectors
Parameters: k (int) – number of eigenvectors to be computed. By default all available eigenvectors. 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 ((n,) ndarray) – Observable vector on the MSM state space 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 \pi_i a_i\]\(\pi=(\pi_i)\) is the stationary vector of the transition matrix \(P\).
-
fingerprint_correlation
(a, b=None, k=None, ncv=None)¶ Dynamical fingerprint for equilibrium time-correlation experiment.
Parameters: - a ((n,) ndarray) – Observable, represented as vector on MSM state space
- b ((n,) 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 computed 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 ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
- amplitudes ((k,) ndarray) – Amplitudes for the correlation experiment
References
Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].
[NDD11] 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
- 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 ((k,) ndarray) – Time-scales (in units of the input trajectory time step) of the transition matrix
- amplitudes ((k,) ndarray) – Amplitudes for the relaxation experiment
References
Spectral densities are commonly used in spectroscopy. Dynamical fingerprints are a useful representation for computational spectroscopy results and have been introduced in [NDD11].
[NDD11] 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.
-
get_model_params
(deep=True)¶ Get parameters for this model.
Parameters: deep (boolean, optional) – If True, will return the parameters for this estimator and contained subobjects that are estimators. Returns: params – Parameter names mapped to their values. Return type: mapping of string to any
-
is_reversible
¶ Returns whether the MSM is reversible
-
is_sparse
¶ Returns whether the MSM is sparse
-
metastable_assignments
¶ Assignment of states to metastable sets using PCCA++
Computes the assignment to metastable sets for active set states using the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute 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: assignments – For each MSM state, the metastable state it is located in. Return type: ndarray (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
-
metastable_distributions
¶ Probability of metastable states to visit an MSM state by PCCA++
Computes the probability distributions of active set states within each metastable set by combining the the PCCA++ method [1]_ with Bayesian inversion as described in [2].
pcca()
needs to be called first to make this attribute 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, 147-179. [2] F. Noe, H. Wu, J.-H. Prinz and N. Plattner. 2013. Projected and hidden Markov models for calculating kinetics and metastable states of complex molecules J. Chem. Phys. 139, 184114.
-
metastable_memberships
¶ Probabilities of MSM states to belong to a metastable state by PCCA++
Computes the memberships of active set states to metastable sets with the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute 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
¶ Metastable sets using PCCA++
Computes the metastable sets of active set states within each metastable set using the PCCA++ method [1]_.
pcca()
needs to be called first to make this attribute 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: sets – A list of length equal to metastable states. Each element is an array with microstate indexes contained in it Return type: list of ndarray 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
-
neig
¶ number of eigenvalues to compute.
-
nstates
¶ Number of active states on which all computations and estimations are done
-
pcca
(m)¶ Runs PCCA++ [1]_ to compute a metastable decomposition of MSM states
After calling this method you can access
metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
.Parameters: m (int) – Number of metastable sets Returns: pcca_obj – An object containing all PCCA quantities. However, you can also ignore this return value and instead retrieve the quantities of your interest with the following MSM functions: metastable_memberships()
,metastable_distributions()
,metastable_sets()
andmetastable_assignments()
.Return type: PCCA
Notes
If you coarse grain with PCCA++, the order of the obtained memberships might not be preserved. This also applies for
metastable_memberships()
,metastable_distributions()
,metastable_sets()
,metastable_assignments()
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
-
pi
¶ The stationary distribution on the MSM states
-
propagate
(p0, k)¶ Propagates the initial distribution p0 k times
Computes the product
\[p_k = p_0^T P^k\]If the lag time of transition matrix \(P\) is \(\tau\), this will provide the probability distribution at time \(k \tau\).
Parameters: - p0 (ndarray(n,)) – Initial distribution. Vector of size of the active set.
- k (int) – Number of time steps
Returns: pk – Distribution after k steps. Vector of size of the active set.
Return type: ndarray(n,)
-
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:
\[a_i = \frac{1}{N_i} \sum_{x_t \in S_i} f(x_t)\]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 accurate (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
\[E_a(k\tau) = \mathbf{p_0}^{\top} \mathbf{P(\tau)}^k \mathbf{a}\]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 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) – 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 5 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
-
reversible
¶ Returns whether the MSM is reversible
-
set_model_params
(P, pi=None, reversible=None, dt_model='1 step', neig=None)¶ Call to set all basic model parameters.
Sets or updates given model parameters. This argument list of this method must contain the full list of essential, or independent model parameters. It can additionally contain derived parameters, e.g. in order to save computational costs of re-computing them.
Parameters: - P (ndarray(n,n)) – transition matrix
- pi (ndarray(n), optional, default=None) – stationary distribution. Can be optionally given in case if it was already computed, e.g. by the estimator.
- reversible (bool, optional, default=None) – whether P is reversible with respect to its stationary distribution. If None (default), will be determined from P
- dt_model (str, optional, default='1 step') –
Description of the physical time corresponding to the model time step. 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*’
- neig (int or None) – The number of eigenvalues / eigenvectors to be kept. If set to None, defaults will be used. For a dense MSM the default is all eigenvalues. For a sparse MSM the default is 10.
Notes
Explicitly define all independent model parameters in the argument list of this function (by mandatory or keyword arguments)
-
simulate
(N, start=None, stop=None, dt=1)¶ Generates a realization of the Markov Model
Parameters: - N (int) – trajectory length in steps of the lag time
- start (int, optional, default = None) – starting hidden state. If not given, will sample from the stationary distribution of the hidden transition matrix.
- stop (int or int-array-like, optional, default = None) – stopping hidden set. If given, the trajectory will be stopped before N steps once a hidden state of the stop set is reached
- dt (int) – trajectory will be saved every dt time steps. Internally, the dt’th power of P is taken to ensure a more efficient simulation.
Returns: htraj – The state trajectory with length N/dt
Return type: (N/dt, ) ndarray
-
sparse
¶ Returns whether the MSM is sparse
-
stationary_distribution
¶ The stationary distribution on the MSM states
-
timescales
(k=None)¶ The relaxation timescales corresponding to the eigenvalues
Parameters: k (int) – number of timescales to be returned. By default all available eigenvalues, minus 1. 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_model
¶ Physical time corresponding to one transition matrix step, e.g. ‘10 ps’
-
transition_matrix
¶ The transition matrix on the active set.
-
update_model_params
(**params)¶ Update given model parameter if they are set to specific values
-