core set Markov state models

[anyschaefer]

Hey,

I would like to have an overview of the current software that is able to generate core set MSMs with milestoning.
Specifically, I am looking for a set up that creates the transition matrix with the time-lagged correlation between forward $m_i^+(t)$ and backward $m_i^-(t)$ milestoning processes, with time $t$ and lagtime $\tau$ .
Basically, I need two main steps:
the projection of discrete trajectory of length $N$ into the milestoning processes

$$m_i^-(t) = \cases{1 & \text{if process is in core set i}\cr 1 & \text{if process came from core set i}\cr 0 & \text{otherwise}}$$

$$m_i^+(t) = \cases{1 & \text{if process is in core set i}\cr 1 & \text{if process go core set i}\cr 0 & \text{otherwise}}$$

and the estimation of the count matrix, with elements

$$\frac{1}{N-\tau} \sum_t m_i^-(t)m_j^+(t+\tau)$$

I know that pyemma has implemented the milestoning method, but if I understand it correctly it is only based on the backward process, right?
Have you integrated the method in deeptime?

All the best,
Joana

[thempel]

Hi Joana,
Yes, the pyemma implementation is based on the backward process only. If I remember correctly, this required some discrete trajectory filtering but otherwise used the same ML-MSM estimator for obtaining the transition matrix from the that. Not sure if there were some unresolved questions regarding the combination of the reversible MSM estimator with the forward-backward formalism that you are describing above. In any case, I don't think it was ever ported to deeptime, even though in my opinion, core-set MSMs are a super useful tool especially for large systems / noisy data.
Cheers,
Tim

[anyschaefer]

Hey Tim,
okay thanks!
We have developed a routine for core set MSMs (kellerlab) in recent years. I would be interested in transferring this to deeptime and directly integrating the reweighting. I'm so curious to see how this works with larger systems.
I'm thinking of extending count_matrix with a new function, and I need to implement the trajectory projection. Can you guide me where to find your discretization methods?
All the best,
joana

[thempel]

Sure! I am not sure if I understand everything correctly though. What you describe as a count matrix in your original post looks like it is normalized -- is it a count matrix or an ML transition matrix estimator? And: What kind of trajectory projection are you doing, i.e., is it a continuous to core-set map or do you start from raw discrete trajectories and discard some states and lump others into core-sets?

Based on my current understanding, this could be a possible implementation:

• Implement functions that extract the forward and backward trajectories from the raw discrete trajectories likely belong into deeptime/markov/_discrete_states.py (I am not sure what this class here is used for, but the place seems right.)
• EITHER A): Extend TransitionCountEstimator's count method (deeptime/markov/_transition_counting.py#L589) with a new kwarg, maybe milestoning_method. Functions from 1) can be used to combine resulting forward and backward trajectories into a count matrix here. Extend the API to treat milestone trajectories and allow access to this new kwarg from top level API of MSM methods.
• OR B): Create a new class, MilestoningTransitionCountEstimator or so, that implements your method. This should inherit most things from the other default estimator. IMO this could be a cleaner solution, but would require some way of accessing different estimators from the top level API. Maybe @clonker has an opinion on that.

Finally some amazing code from pyemma -- maybe some of it can be recycled.

• Possibly efficient implementation of assigning unassigned states to cores:
  https://github.com/markovmodel/PyEMMA/blob/3327f28b49e388e1ce4a6a83ab2f0c0ac7ca5050/pyemma/msm/estimators/_dtraj_stats.py#L131-L138
• Function that assigns microstates to core-sets: https://github.com/markovmodel/PyEMMA/blob/3327f28b49e388e1ce4a6a83ab2f0c0ac7ca5050/pyemma/util/discrete_trajectories.py#L227
• Some unit test cases: https://github.com/markovmodel/PyEMMA/blob/3327f28b49e388e1ce4a6a83ab2f0c0ac7ca5050/pyemma/util/tests/test_discrete_trajectories.py#L199

[anyschaefer]

Hey Tim,

thanks for your tips, they were helpful. I have worked out the rough framework for the implementation, see below.

Now I have run into a problem. In all existing deeptime models, you assume that the matrix on which we build the MSM is the pure transition matrix, which should be positive in all entries and row-normalised. This assumption is backed up in higher level ( in MarkovStateModel in self.update_transition_matrix(transition_matrix)) to deep inside to the code in the functions defining the stationary distribution (in deeptime/markov/tools/analysis and so on).

I would say the condition of positive entries does not apply to the csMSM, since the matrix passed into the deeptime/markov/msm/MarkovStateModel is the dot product of transition matrix ($\mathbf{P}$) and inverse of the ‘mass matrix’ ($\mathbf{M}^{-1} $)
$\mathbf{c}^T\mathbf{P}(\tau)\mathbf{M}^{-1} = \lambda\mathbf{c}^T,$

with
$P_{ij}(\tau)=\frac{1}{N_T-n_T}\sum^{N_T-n_T}_{t=1} m_i^-(t)m_j^+(t+n_T)$

$M_{ij}=\frac{1}{N_T}\sum^{N_T}_{t=1}m_i^-(t) m_j^+(t)$

So my problem is that whenever I have a negative entry in $\mathbf{P}(\tau)\mathbf{M}^{-1}$, the deeptime MSM model refuses to analyse it.

I now see two possibilities to solve this problem:

• First, to touch the MarkovStateModel in update_transition_matrix(transition_matrix) and lower level functions:
  T_nonrev = transition_matrix_non_reversible(C)
from ...analysis import stationary_distribution
pi = stationary_distribution(T_nonrev, transition_matrix_negative) # ADD HERE something to allow negative matrix elements
  (although I don't think you can simply use a modified form update_transition_matrix(transition_matrix) because all basic functions for stationary distribution, eigenvalue and eigenvector calculation test the negativity of the input matrix)

• or create a 'new MarkovStateModel' for the core set MSM
  there I would only evaluate eigenvectors and eigenvalues because I can't/wouldn't code all functionalities of deeptime MSM model again

Or do you think there is another way to allow a negative input to MarkovStateModel?

all the best
Joana

---

framework for the implementation

As you suggested, I have added in:

• deeptime/markov/_milestoning_transition_counting.py
  MilestoningModel(TransitionCountModel)
MilestoningEstimator(TransitionCountEstimator)
  -> extends the class for the overlap matrix needed in a core set MSM:

• deeptime/markov/tools/estimation/api.py
  milestoning_count_matrix
transition_matrix
  -> basic function to implement transition matrix ($\mathbf{P}$) and inverse of the ‘mass matrix’ ($\mathbf{M}^{-1} $)

• deeptime/markov/tools/estimation/dense/sparse
  transition_matrix_mass_matrix_weighted
  -> calculate dot product of transition matrix ($\mathbf{P}$) and inverse of the ‘mass matrix’ ($\mathbf{M}^{-1} $)

• deeptime/markov/_discrete_states.py
  MilestoningProcesses()
  -> to project the clustered MD output to the milestoning processes $m_i^-(t)$ and $m_j^+(t)$$

deeptime/markov/
-> more changes in utils, transition matrix and types moduels

• deeptime/markov/msm/
  CoreSetMSM(MaximumLikelihoodMSM)
  -> Core Set estimator for MSMs (:class:MarkovStateModel <deeptime.markov.msm.MarkovStateModel>) given milestoning processe.

[thempel]

Hi Joana, sorry for the delay, unfortunately I didn't come up with great ideas in the meantime.

In my memory, most of the methods that we inherited from pyemma are defined only for stochastic matrices, i.e., all transition matrix elements need to be >=0. If I recall correctly, for non-positive transition matrices matrices, the eigenvalue decomposition generally yields results that are not following the assumption of downstream methods anymore (afair eigenvalues do not need to be real valued eigenvalues with largest eigenvalue 1, etc). Thus, apart from the changes to the code, we would need to make sure that all methods that have this requirement only get matrices that are valid for them, or that we only remove the assertions when this is actually allowed. What analyses are you planning to do with the csMSM transition matrix, and are they still defined for your case?