{'fraction_non_alphanumeric': 0.027668179128033144, 'fraction_numerical': 0.006904714933912014, 'mean_word_length': 4.908216783216783, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improved coarse-graining of Markov state models via explicit consideration of statistical uncertainty

Markov state models (MSMs)---or discrete-time master equation models---are a powerful way of modeling the structure and function of molecular systems like proteins. Unfortunately, MSMs with sufficiently many states to make a quantitative connection with experiments (often tens of thousands of states even for small systems) are generally too complicated to understand. Here, I present a Bayesian agglomerative clustering engine (BACE) for coarse-graining such Markov models, thereby reducing their complexity and making them more comprehensible. An important feature of this algorithm is its ability to explicitly account for statistical uncertainty in model parameters that arises from finite sampling. This advance builds on a number of recent works highlighting the importance of accounting for uncertainty in the analysis of MSMs and provides significant advantages over existing methods for coarse-graining Markov state models. The closed-form expression I derive here for determining which states to merge is equivalent to the generalized Jensen-Shannon divergence, an important measure from information theory that is related to the relative entropy. Therefore, the method has an appealing information theoretic interpretation in terms of minimizing information loss. The bottom-up nature of the algorithm likely makes it particularly well suited for constructing mesoscale models. I also present an extremely efficient expression for Bayesian model comparison that can be used to identify the most meaningful levels of the hierarchy of models from BACE.


I. INTRODUCTION
Markov state models (MSMs) are a powerful means of understanding dynamic processes on the molecular scale, like protein folding and function [1][2][3]. These discrete-time master equation models consist of a set of states-akin to local minima in the system's free energy landscape-and a matrix of transition probabilities between them. Typically, the states are identified via a kinetic clustering of molecular dynamics simulations.
Unfortunately, building MSMs and extracting understanding from them is still a challenging task. Ideally, MSMs would be constructed using a purely kinetic clustering of a simulation data set. Calculating the transition rate between two conformations is an unsolved problem though, so, a number of alternative methods for building MSMs have been developed [4][5][6][7][8][9]. Many of these approaches have converged on a two-stage process. First, the conformations sampled are clustered into microstates based on geometric criteria such that the degree of geometric similarity between conformations in the same state implies a kinetic similarity. Such models are excellent for making a quantitative connection with experiments because of their high temporal and spatial resolution. However, it is difficult to examine such models to gain an intuition for a system because the rugged nature of most biomolecules free energy landscapes requires that the initial microstate model have tens of thousands of states. Therefore, in a second stage, the initial state space is coarse-grained by lumping rapidly interconverting-or kinetically close-microstates together into macrostates to obtain a more compact and comprehensible model. Reasonable methods are now available for the first stage of this procedure [4][5][6][7][8][9], but there is still a need for more efficient and accurate methods for coarse-graining MSMs.
Most all of these methods operate on the maximum-likelihood estimate of the transition probability matrix and do not account for statistical uncertainty in these parameters due to finite sampling. For example, PCCA(+) uses the eigenspectrum of the transition matrix to find the partitioning that best captures the slowest transitions. Such methods are well suited to data-rich situations but often fail when poorly sampled transitions are present [13].
For example, Fig. 1 shows a case where PCCA fails due to a few poorly sampled transitions.  thin lines represent 1 transition. Therefore, the best coarse-graining into three states is to merge states 0-2, 3-5, and 6-8 because transitions within these groups are fast compared to transitions between the groups. BACE correctly identifies this optimal coarse-graining into three macrostates.
Eigenspectrum-based methods also have trouble creating mesoscale models-models with a large number of macrostates that are still quantitatively predictive yet are significantly more compact than the original microstate model-due to issues like propagating error.
Here, I present a Bayesian agglomerative clustering engine (BACE) for coarse-graining MSMs in a manner that accounts for model uncertainty and can easily create mesoscale models. Bayesian methods have found wide applications in the physical sciences, and in MSMs in particular [17][18][19][20], for their ability to deal with uncertainty. Inspired by the hierarchical nature of biomolecules' free energy landscapes, BACE performs an agglomerative clustering of microstates into macrostates by iteratively lumping together the most kinetically similar states-i.e. the most rapidly mixing states. The key equation derived here is a closed-form expression for a Bayes factor that quantifies how likely two states are to be kinetically identical. This expression is related to the relative entropy [21], an information theoretic measure that has found numerous applications in the physical sciences [22][23][24].
Indeed, the expression is actually equivalent to the generalized Jensen-Shannon divergence [25], which can be interpreted as the average amount of information a single measurement gives about which of two possible distributions is being sampled. I also present an approximate expression for model comparison that allows one to identify the most informative levels of the hierarchy of models generated with BACE. These methods could be applied directly to other Markov processes and could also be extended to other probabilistic models.

II. BACE ALGORITHM
The hierarchical structure of biomolecules' free energy landscapes naturally suggests a hierarchical approach to model construction. The free energy landscapes of almost all biomolecules are extremely rugged, having numerous local minima separated by barriers of different heights. Put another way, free energy basins in this landscape can typically be subdivided into smaller local minima, giving rise to a hierarchy of minima. Transitions across low barriers occur exponentially more often than those across higher barriers. Groups of local minima separated by low barriers will mix rapidly and, therefore, appear as a single larger state to other minima separated from them by larger barriers. Thus, these groups can satisfy a requirement for coarse-graining models called lumpability [27]. A microstate MSM is considered lumpable with respect to some set of macrostates if and only if, for every pair of macrostates M 1 and M 2 and any pair of microstates i and j in M 1 k∈M 2 p ik = k∈M 2 p jk where p ij is the probability that the system will transition to state j given that it is currently in state i.
We can exploit the concept of lumpability to construct coarse-grained models by progressively lumping together the most kinetically similar states-i.e. those with similar transition probabilities. Physically, this is equivalent to merging states that mix rapidly because they are only separated by a low free energy barrier. One might be tempted to use an L1 or L2 norm between the transition probabilities out of each pair of states to determine which are most similar. However, such an approach would ignore the fact that some states and transitions are better sampled than others and, therefore, would be susceptible to the same pitfalls as eigenspectrum-based methods.
I propose a Bayesian method for determining which states to lump together. Specifically, I propose to employ a Bayes factor comparing how likely the data observed for a pair of states is to have come from either different (P (different|C)) or the same (P (same|C)) underlying distribution of transition probabilities where C is the matrix of transition counts observed between all pairs of states. Bayes factors compare the evidence (or marginal likelihood, P (Model|Data)) for two different models. In Finally, it is valuable to derive an approximate expression for the BACE Bayes factor. One could evaluate the Bayes factor by sampling from the posterior distribution for each state.
However, doing so would require a number of calculations for every comparison of a pair of states. A single, closed-form expression-like the one derived in the next section-is significantly more efficient.
The BACE algorithm is then 4. Repeat steps 2 and 3 until only two states remain.
One could also stop the algorithm when the BACE Bayes factor reaches a certain threshold. For example, a log 10 (Bayes factor) of 1 indicates that the model in the numerator is significantly more likely (over ten times more likely) than the one in the denominator.
Therefore, if the minimum BACE Bayes factor between any pair of states reaches 1, then one could infer the any further merging of states would greatly reduce the quantitative accuracy of the model and stop the algorithm. However, if one's objective is to understand a system then continuing to merge states may be of great value. The resulting models will only be qualitatively correct, at best. However, their simplicity may allow more insight.
Hypotheses generated with these simple, qualitative models can then be tested with more complex, quantitative models and, ultimately, with experiments.
Future improvements to the algorithm could be made by including more complex moves.
For example, the current algorithm is greedy and, therefore, can never recover if two states are mistakenly merged. One could correct such mistakes by allowing microstates to move to more appropriate macrostates or by iteratively breaking macrostates apart and rebuilding them, as in Ref. [4]. Such moves are not undertaken here as they would reduce the efficiency of the method. Moreover, the present greedy algorithm performs quite well compared to other methods, as discussed in the Results section.

III. BACE BAYES FACTOR
To make the BACE algorithm computationally efficient, I derive a closed-form expression for the log of the BACE Bayes factor from Eq. 1. The final expression for the BACE Bayes factor is where C is the transition count matrix,Ĉ i is the number of transitions observed from state i, D(p i q) = k p ik log p ik q k is the relative entropy between probability distribution p i and q, p i is a vector of maximum likelihood transition probabilities from state i, and q =Ĉ i p i +Ĉ j p ĵ C i +Ĉ j is the vector of expected transition probabilities from combining states i and j. Note that this expression includes a comparison between p ij and q j that helps prevent the merger of disconnected states. For example, consider the simple model A ↔ B ↔ C (A and C are disconnected). If the BACE Bayes factor only compared transition probabilities to states other than the two being considered, then one could easily obtain lumpings like {A, C} and {B}. However, comparing the self-transition probabilities and exchange probabilities between the states being compared helps avoid these pathological situations (e.g. in this case p AA > 0 while p CA = 0, so these states are unlikely to appear kinetically close).
This expression is equivalent to the generalized Jesnen-Shannon divergence [25] and, therefore, has an appealing information theoretic interpretation. Given a sample drawn from one of two probability distributions, the Jensen-Shannon divergence is the average information that sample provides about the identity of the distribution it was drawn from [28]. The result is zero if the two distributions are equivalent and reaches its maximal value if the distributions are non-overlapping and a single data point, therefore, uniquely specifies which distribution it was drawn from. In this case, the larger the Bayes factor is, the more likely the data for each state are to have come from different underlying distributions. By iteratively merging the most kinetically similar states, BACE retains the most divergent states, which can be interpreted as keeping the states with the most information content.
To derive Eq. 2, we first recognize that every possible set of transition probabilities out of some initial state that satisfies 0 ≤p ij ≤ 1 and jp ij = 1 has some probability of generating the observed transitions out of that state. From Bayes rule, the posterior probability of some distribution (p i ) being the true underlying distribution given a set of observed transitions is where C i is a vector of transition counts out of state i and α i will be discussed shortly.
Assuming the transition probabilities for each state are independent, we can use a multinomial distribution for the likelihood A Dirichlet prior (D) is chosen as it is conjugate to the multinomial likelihood. That is, if the prior is a Dirichlet then the posterior is also a Dirichlet. The prior is then where α i is a vector of pseudocounts giving the expected number of transitions before any data is observed. We choose α ik = 1/n where n is the number of states because for a state to exist we must have observed at least one transition originating from that state and, prior to observing any data, the chance that that transition is to any particular state is equal [17,23].
Combining the expressions for the likelihood and prior, the posterior distribution from We can now calculate the log of the evidence for a particular model (M) log P (C i |M ) = log where H(p i ) = − k p ik log p ik is the entropy of p i and we have made the substitutionŝ C i = k C ik , p ik = C ik /Ĉ i (the maximum likelihood estimate of the transition probability), Γ(C ik + 1/n) ≈ Γ(C ik + 1) = C ik !, Γ(1/n) ≈ n, and Stirling's approximation. Note that the approximations made between Eqs. 4 and 5 break down for small sample sizes but this can be ignored safely. One could calculate the evidence more accurately by directly evaluating Eq. 4. However, this could lead to numerical errors as the Γ function tends to blow-up for the large inputs one is likely to encounter in real-world applications of this method.
Moreover, the closed-form expression for the Bayes factor based on Eq. 5 performs quite well in practice, as discussed in the Results section.
The BACE Bayes factor given in Eq. 2 is then the ratio of the evidence for the transition counts from states i and j coming from two different distributions versus a single distribution (log P (different|C) P (same|C) = log P (C|different)P (different)

P (C|same)P (same)
) where we assume the prior probabilities for the two models are equal and drop terms depending only on n as they simply introduce a constant that has no effect on the relative ordering of Bayes factors comparing various states.

IV. APPROXIMATE BAYESIAN MODEL COMPARISON
Bayesian model comparison is a powerful way to determine which of two models best explains a set of observations. Such methods are of great value here as they can be used to compare the results of BACE to other coarse-graining methods. Moreover, they can be used to decide which levels of the hierarchy of models from BACE are most deserving of further analysis. However, current methods [20] are too computationally demanding for this second task.
Using similar mathematical machinery to that employed in the derivation of BACE and paralleling the derivation in Ref. [20], we can also derive a closed-form expression for the log of the Bayes factor comparing two coarse-grainings-or lumpings-of an MSM, L 1 and where B and C are the transition count matrices at the macrostate and microstate levels, BACE is much better at dealing with statistical uncertainty in model parameters than current eigenspectrum-based methods. For example, it is able to correctly identify the three macrostates in the simple model shown in Fig. 1 even in the presence of the poorly sampled transitions that confound eigenspectrum-based methods. BACE also naturally lumps states with few samples into larger ones, whereas eigenspectrum-based methods tend to make such states into singleton macrostates. With BACE, a significantly better sampled state will dominate the Bayes factor when compared to a poorly sampled state, leading to a high likelihood that the poorly sampled state will be absorbed into its better sampled neighbor.
Methods like super-level-set hierarchical clustering (SHC) [13] and the most probable path algorithm [16] also deal with poorly sampled states by merging them into larger states.
However, in these methods, a state is considered poorly sampled if its population is below some user-defined cutoff. BACE offers the advantage of naturally identifying poorly sampled states without any reliance on user-defined input.
Beyond this qualitative improvement, a quantitative measure of model validity shows that coarse-grainings from BACE capture both the thermodynamics and kinetics of molecular systems better than existing eigenspectrum-based methods (Table I) Table I shows that BACE is typically many orders of magnitude better than eigenspectrum-based methods by this metric. Such quantitative comparisons are crucial because the complexity of most real-world MSMs renders a qualitative assessment of a coarse-graining's validity impossible.
Another advantage of BACE is that it generates an entire hierarchy of models. Having this hierarchy makes it possible to look for general properties that are robust to the degree of coarse-graining and, therefore, may be important properties of the system being investigated.
In addition, having this hierarchy allows the user to determine how many macrostates are appropriate to use. In theory, one could employ the Bayesian model comparison method  Fig. 1 to a full protein, the villin headpiece. The same number of states is used for each method. The numbers reported are the log 10 of the Bayes factor comparing how likely the coarse-graining from BACE is to have generated a given data set to how likely the model from PCCA (or PCCA+) is to have generated the data. Numbers greater than one suggest the model from BACE is significantly better at explaining the observed data than the models from other methods whereas numbers less than -1 would suggest the other methods are significantly better than BACE. These values were calculated with the model comparison method from Ref b Model with 181 microstates and 6 macrostates from Ref. [26].
accounting for reversibility from Ref. [20] to decide which levels of the hierarchy are most deserving of further analysis but, in practice, this would be impractical due to the time requirements of that method. However, both the BACE Bayes factor and the approximate model comparison method presented here correlate well with the reversible method ( Fig.   2) and, therefore, can be used to guide which levels of the hierarchy are pursued further.
Each Bayes factor changes more rapidly when more distinct states are lumped together, so models immediately preceding these dramatic jumps are ideal for further analysis. The BACE Bayes factor can even be used to visualize the hierarchical nature of a system's free energy landscape and choose appropriate levels for further analysis (Fig. 3). One could also combine the methods by using the approximate expression to guide the application of the

VII. SIMULATION DETAILS
The alanine dipeptide data used for Table I was taken from Ref. [26]. One hundred simulations were performed with Gromacs 4.5 [30] using the Amber96 force field [31] and the OBC GBSA implicit solvent [32]. Each trajectory is 500ps long, with conformations stored every 1ps.
The villin data used for Table I was taken from Ref. [33] and the macrostate definitions were taken from Ref. [29]. Five hundred simulations were performed with Gromacs deployed on the Folding@home distributed computing environment [33,34]. The Amber03 force field [35] and Tip3p explicit solvent were used. Each trajectory is up to 2µs long, with conformations stored every 50ps.