MSM lag time cannot be used for variational model selection
Brooke E. Husic Vijay S. Pande








Department of Chemistry, Stanford University, Stanford CA 94305, USA





The variational principle for conformational dynamics has enabled the systematic construction of Markov state models through the optimization of hyperparameters by approximating the transfer operator.
In this note we discuss why lag time of the operator being approximated must be held constant in the variational approach.

MSM lag time cannot be used for variational model selection
    Vijay S. Pande
    December 30, 2023
===========================================================



Markov state models (MSMs) are a powerful master equation framework for the analysis of molecular dynamics (MD) datasets that involve a complete partition of the conformational space into disjoint states. <cit.> By representing each frame of a MD dataset as its state label, the populations of and conditional pairwise transition probabilities between the states can be counted, leading to thermodynamic and kinetic information about the system, respectively. This information is represented by a transition matrix, which contains all the information necessary to propagate the system forward in time.
The transition matrix is the discrete-time approximation to the transfer operator 𝒯(τ), which is characterized by its lag time τ. The transfer operator propagates the system, represented by a normalized probability density u_t(x), forward by a time step of τ, and admits a decomposition into eigenfunctions and eigenvalues (see Ref. Bowman_Book14, Ch. 3),



    𝒯(τ) ∘ u_t(x) = u_t+τ(x), 
       𝒯(τ) ∘ψ_i = λ_i ψ_i.



The eigenvalues λ_i are real and numbered in decreasing order. The unique highest eigenvalue λ_1=1 corresponds to the stationary distribution, and the subsequent eigenvalue/eigenfunction pairs represent dynamical processes in the time series. Importantly, the timescale of each process can only be retrieved with knowledge of the lag time at which the operator was defined using the equation,


    t_i = -τ/logλ_i.


Choosing a lag time at which the system is Markovian depends on what type of system is being modeled. At a long enough lag time for the system to be approximated as a Markov process, intrastate transitions occur much more quickly than interstate transitions. The appropriate lag time depends on the system of study: for protein folding, 50 ns might be appropriate; for electron dynamics, a suitable lag time might be on the order of femtoseconds.
If a system is Markovian at a lag time τ (if the intrastate transitions occur more quickly than τ), then the system will be Markovian at all lag times greater than τ and the timescales of the subprocesses will be constant for all Markovian lag times. This idea has motivated the use of implied timescales plots to choose a lag time. <cit.> Lag times after which the timescales “level out” are assumed to be Markovian, and usually the shortest such time is chosen for the most temporal resolution.

In practice, we usually do not know the true eigenfunctions ψ in (<ref>) and instead need to guess them. For a MSM, this means choosing how to divide phase space into disjoint states. Until recently, choosing how to define the states occupied by a dynamic system represented a bottleneck in MSM methods development, and heuristic, hand-selected states were common. However, the derivation of a variational principle for conformational dynamics by <cit.> in 2013 opened the door for a systematic approach to choosing the states of a system. 
Our guess, or ansatz eigenfunctions, ψ̂_i will admit corresponding eigenvalues λ̂_i. Using our ansatz, we can state the variational principle derived by Noé and Nüske, <cit.>


    GMRQ≡∑_i=1^mλ̂_̂î≤∑_i=1^mλ_i,


where GMRQ stands for generalized matrix Rayleigh quotient, which is the form of the approximator when the first m eigenfunctions are estimated simultaneously. By recalling the relation of the eigenvalues and operator lag time τ to the system timescales in (<ref>), we see that the variational principle establishes an upper bound on the timescales of the slowest m processes in the dynamical system. In practical cases, the variational bound can be exceeded due to statistically undersampled processes; therefore, the GMRQ must be evaluated under cross-validation as described in Ref. McGibbon_JChemPhys15a.

When we variationally choose a set of eigenfunctions, we can only compare them if we are trying to approximate the same transfer operator. Therefore, the lag time τ must not be changed when the ansatz is changed, and cannot be variationally optimized using the GMRQ—instead, it must be determined using such techniques as implied timescale plots.
In contrast, all transformation and dimensionality reduction choices leading up to the state decomposition are ideal hyperparameters to optimize using the GMRQ. This might include:


	
  * RMSD cutoffs for geometric clustering;
    
  * internal coordinate choices such as dihedral angles or contacts pairs, including which angles and pairs to include, and any transformations thereof; 
    
  * internal parameters for time-structure based independent component analysis (tICA) such as tICA lag time, number of components retained, and any transformations of these components;
    
  * clustering algorithm and number of clusters;

but, as discussed above, cannot include:


	
  * the operator lag time, or
    
  * the number of timescales scored.

These choices are illustrated in Fig. <ref>. For protein folding, we refer the reader to Ref. Husic_JChemPhys16 for a systematic study of these choices in the context of the VAC. 



In practice, we recommend starting with reasonable parameters for the system of study and choosing a valid lag time. Then, at the chosen lag time, perform a hyperparameter search for any of the state decomposition choices listed above. This two step process can then be repeated, alternating lag time validation using fixed hyperparameters with hyperparameter searches for a fixed lag time.

Perhaps the most natural way to understand

the separate treatment of the lag time is to consider its true role in kinetic model building.  While previous MSM approaches have treated it effectively like a hyperparameter (e.g., choosing a lag time based on flattening of implied timescales), in actuality, this approach is fundamentally philosophically incorrect. The lag time must be chosen a priori by the researcher, as it directly reflects the resolution of interest to study. Given a method which can directly identify the relevant degrees of freedom, choosing a lag time of picoseconds would bring water dynamics into the state space, vs. nanoseconds for backbone and side chain dynamics or microseconds for longer time slow scale rearrangements. For this reason, it simply doesn’t make sense to let the model choose the lag time and instead one must have the protocol choose the best model given a pre-chosen set lag time.

MSMs are just one example of the general set of models to which the VAC applies. The popular tICA framework <cit.> can also be variationally optimized. When using tICA as an intermediate step in MSM construction, the tICA lag time may be varied and optimized. However, in the case where the tICA model is the entity being evaluated, the tICA lag time and number of components scored must be held constant in order to ensure that the same operator is being approximated. Additional extensions of the VAC can be found in Ref. Wu_Arxiv17. We also refer the interested reader to Ref. McGibbon_JChemPhys15b, which presents continuous-time Markov processes that do not have lag times.
We would also like to note that the VAC is not a panacea: the slowest dynamical processes are often assumed, but not guaranteed, to be the processes of interest, and it is important to verify this for each analysis. We anticipate this note will help guide hyperparameter optimization when using VAC. The open-source software Osprey <cit.> has been designed for this type of optimization and is available on <msmbuilder.org>.


The authors are grateful to Matt Harrigan, Carlos Hernández, Jade Shi, Anton Sinitskiy, Nate Stanely, and Muneeb Sultan for discussion and manuscript feedback.
We acknowledge the National Institutes of Health under No. NIH R01-GM62868 for funding. V.S.P. is a consultant & SAB member of Schrodinger, LLC and Globavir, sits on the Board of Directors of Apeel Inc, Freenome Inc, Omada Health, Patient Ping, Rigetti Computing, and is a General Partner at Andreessen Horowitz. 











10
fxundefined [1]
 ifx#1

fnum [1]
 #1firstoftwo
 secondoftwo
 
fx [1]
 #1firstoftwo
 secondoftwo
 





noop [0]secondoftwo

ref[1]@startlink#1@href
href[1]#1@endlink
anitize@url [0]`
12`$12`&12`#12`1̂2`_12`%12
startlink[1]
endlink[0]

rl [1]href #1












@bib@innerbibempty

[Bowman, Pande, and Noé(2014)]Bowman_Book14
  author author G. R. Bowman, author V. S. Pande,   and author F. Noé  (publisher Springer, year 2014)NoStop
[Swope, Pitera, and Suits(2004)]Swope_JPCB04a
  author author W. C. Swope, author J. W. Pitera,   and author F. Suits, 10.1021/jp037421y journal journal J.
  Phys. Chem. B volume 108, pages 6571
  (year 2004)NoStop
[Noé and Nüske(2013)]Noe_MMS13
  author author F. Noé and author F. Nüske, 10.1137/110858616 journal
  journal Multiscale Model. Simul. volume
  11, pages 635 (year 2013)NoStop
[McGibbon and Pande(2015a)]McGibbon_JChemPhys15a
  author author R. T. McGibbon and author V. S. Pande, http://scitation.aip.org/content/aip/journal/jcp/142/12/10.1063/1.4916292
  journal journal J. Chem. Phys. volume 142, eid 124105 (year
  2015a)NoStop
[Husic et al.(2016)Husic,
  McGibbon, Sultan, and Pande]Husic_JChemPhys16
  author author B. E. Husic, author R. T. McGibbon,
  author M. M. Sultan,  and  author V. S. Pande, http://scitation.aip.org/content/aip/journal/jcp/145/19/10.1063/1.4967809;jsessionid=Td617tvPkXmj4MqYBi1JOe8C.x-aip-live-02
  journal journal J. Chem. Phys. volume 145, eid 194103 (year
  2016)NoStop
[Pérez-Hernández et al.(2013)Pérez-Hernández, Paul,
  Giorgino, De Fabritiis, and Noé]PerezHernandez_JChemPhys13
  author author G. Pérez-Hernández, author F. Paul, author T. Giorgino, author G. De Fabritiis,  and author F. Noé, http://dx.doi.org/10.1063/1.4811489 journal journal J. Chem. Phys. volume
  139, eid 015102 (year 2013)NoStop
[Schwantes and Pande(2015)]Schwantes_JCTC15
  author author C. R. Schwantes and author V. S. Pande, 10.1021/ct5007357 journal
  journal J. Chem. Theory Comput. volume
  11, pages 600 (year 2015)NoStop
[Wu and Noé(2017)]Wu_Arxiv17
  author author H. Wu and author F. Noé,  @noop  journal journal arXiv preprint
  arXiv:1707.04659  (year 2017)NoStop
[McGibbon and Pande(2015b)]McGibbon_JChemPhys15b
  author author R. T. McGibbon and author V. S. Pande, http://dx.doi.org/10.1063/1.4926516 journal journal J. Chem. Phys. volume
  143, eid 034109 (year
  2015b)NoStop
[McGibbon et al.(2016)McGibbon, Hernández, Harrigan,
  Kearnes, Sultan, Jastrzebski,
  Husic, and Pande]McGibbon_JOSS16
  author author R. T. McGibbon, author C. X. Hernández, author M. P. Harrigan, author S. Kearnes,
  author M. M. Sultan, author S. Jastrzebski, author
  B. E. Husic,  and author
  V. S. Pande, http://dx.doi.org/10.21105/joss.00034 journal journal J. Open Source Software volume 1
  (year 2016)NoStop


