
==== Front
J Phys Chem B
J Phys Chem B
jp
jpcbfk
The Journal of Physical Chemistry. B
1520-6106
1520-5207
American Chemical Society

39145603
10.1021/acs.jpcb.3c08213
Article
Kemeny Constant-Based Optimization of Network Clustering Using Graph Neural Networks
Martino Sam Alexander
https://orcid.org/0000-0002-5003-2233
Morado João
Li Chenghao
Lu Zhenghao
https://orcid.org/0000-0002-9823-4766
Rosta Edina *
Department of Physics and Astronomy, University College London, London WC1E 6BT, U.K.
* Email: e.rosta@ucl.ac.uk.
15 08 2024
29 08 2024
128 34 81038115
16 12 2023
08 07 2024
28 06 2024
© 2024 The Authors. Published by American Chemical Society
2024
The Authors
https://creativecommons.org/licenses/by/4.0/ Permits the broadest form of re-use including for commercial purposes, provided that author attribution and integrity are maintained (https://creativecommons.org/licenses/by/4.0/).

The recent trend in using network and graph structures to represent a variety of different data types has renewed interest in the graph partitioning (GP) problem. This interest stems from the need for general methods that can both efficiently identify network communities and reduce the dimensionality of large graphs while satisfying various application-specific criteria. Traditional clustering algorithms often struggle to capture the complex relationships within graphs and generalize to arbitrary clustering criteria. The emergence of graph neural networks (GNNs) as a powerful framework for learning representations of graph data provides new approaches to solving the problem. Previous work has shown GNNs to be capable of proposing partitionings using a variety of criteria. However, these approaches have not yet been extended to Markov chains or kinetic networks. These arise frequently in the study of molecular systems and are of particular interest to the biomolecular modeling community. In this work, we propose several GNN-based architectures to tackle the GP problem for Markov Chains described as kinetic networks. This approach aims to maximize the Kemeny constant, which is a variational quantity and it represents the sum of time scales of the system. We propose using an encoder-decoder architecture and show how simple GraphSAGE-based GNNs with linear layers can outperform much larger and more expressive attention-based models in this context. As a proof of concept, we first demonstrate the method’s ability to cluster randomly connected graphs. We also use a linear chain architecture corresponding to a 1D free energy profile as our kinetic network. Subsequently, we demonstrate the effectiveness of our method through experiments on a data set derived from molecular dynamics. We compare the performance of our method to other partitioning techniques, such as PCCA+. We explore the importance of feature and hyperparameter selection and propose a general strategy for large-scale parallel training of GNNs for discovering optimal graph partitionings.

H2020 European Research Council 10.13039/100010663 757850 Engineering and Physical Sciences Research Council 10.13039/501100000266 EP/R013012/1 document-id-old-9jp3c08213
document-id-new-14jp3c08213
ccc-price
Special Issue

Published as part of The Journal of Physical Chemistry Bvirtual special issue “Gregory A. Voth Festschrift”.
==== Body
pmc1 Introduction

Markov chains (MCs), and their extensions, have become a key tool across multiple domains to model dynamical systems.1,2 In their simplest form, MCs incorporate sequential data as a Markovian process into a stochastic model, where the probability of observing the next event depends solely on the preceding event. By eliminating the long-timescale dependencies that often occur in complex systems, MCs offer a robust and formal method to interpret high-dimensional data sets in the emergent era of big data. The simplicity, transferability, and capability of MCs reduce the difficulties associated with modeling and interpreting ordered dynamics and have led to a wide range of use cases. These applications encompass modeling character sequences in computer science and genomics,3,4 predicting web-user activity,5 and modeling time-dependent processes in physics and chemistry,6,7 to name just a few.

An application that has widely adopted these techniques is a biomolecular simulation, where the Markov state model (MSM) formalism has emerged as a strategy for interpreting the results of molecular dynamics (MD) simulations as transitions between aggregated conformational states.8−10 The surge in popularity of the MSM mirrors the increased ability to quickly generate large quantities of MD data, providing a useful tool to condense massive data sets into a small, intuitive probabilistic format. The MSM offers intuitive descriptions of complex dynamic systems as transitions between metastable minima and saddle points serving as transition states,11 helping to elucidate the underlying functional kinetics of interest.12,13 Over the past few decades, several domain-specific advances led to the development of techniques to automatically generate MSM representations of molecular systems. Simultaneously, these techniques have been employed to further sample the multitude of different configurations of biomolecular systems.14,15

Despite being smaller than the huge data sets from which they originate, many MSMs remain high-dimensional, potentially consisting of hundreds of states. This large number of states is often unmanageable for downstream tasks such as drug design and biomolecular engineering.16 As size often becomes the limiting factor to establishing an intuitive understanding of systems, many different approaches have been proposed to reduce the size of MCs and generate optimal ones.17−20 However, the ability to generate large data sets is only increasing, as modern coarse-graining protocols21 and new developments in machine learning22 allow for larger time scales to be probed. Traditional graph partitioning (GP) algorithms, such as spectral clustering, k-means, or multilevel methods, and their derivatives, have been widely employed to address this problem.23−25 However, these algorithms often face limitations in capturing the complex and intricate dependencies on the global structure in MCs. Additionally, the generality of these methods means that they often disregard the critical dynamic information held within MCs. Alternative Markov-Chain Monte Carlo (MCMC) approaches have been applied for dimensionality reduction on MCs, which can be used with most clustering criteria.26−28 Nevertheless, these approaches encounter difficulties in efficiently sampling the most relevant regions of partitioning space, resulting in impractical computation times and frequent overlooking of the most optimal solutions.

In this article, we introduce machine learning (ML) models based on graph neural networks (GNN) to optimize dynamical network clusterings. We build upon previous work29 for reducing the dimensionality of MCs, and propose ML architectures leveraging recent developments in using GNNs to solve the more general GP problem.30−32 Although several existing approaches use GNNs to minimize partitioning cuts33 or modularity,34 these applications do not extend to kinetic graphs or optimization based on dynamical criteria. Additionally, these approaches mostly focus on providing approximate partitionings for large graphs rather than exploring their potential for identifying optimal partitionings in smaller networks. We extend these methods by representing MCs as kinetic networks and show how performing gradient descent on a large number of randomly initialized small GNNs can yield optimal network partitionings by using elaborate optimization criteria. As a metric to quantitatively assess the degree to which a clustering preserves the underlying dynamics of the MSM, we use the change of the Kemeny constant (KC). KC is an inherently dynamical quantity derived from mean first passage times (MFPTs) associated with any kinetic network, including those used within biomolecular simulations.35 We also compare our clustering results with a commonly used kinetic network clustering method for MSMs, the robust Perron-cluster analysis (PCCA+), which does not use the KC as the optimization criterion.36

The proposed MC clustering method can produce a faithful low-dimensional replication of a system’s underlying dynamics while also providing a quantifiable measurement related to the amount of kinetic information lost. Along with outlining a general strategy for the large-scale training of GNNs for graph partitioning, we show our method’s performance when minimizing the change in KC on several different kinds of synthetic graphs and an example derived from MD simulations. We also tested more complex GNN architectures, showing how large-scale training of simple GNNs is the most effective for the task and hypothesizing why this is the case.

2 Theory

2.1 Graphs

A graph G with N nodes can be described as a set of vertices V = {v1, ..., vN} and edges eij ∈ E, which are often represented as a binary adjacency matrix, A, where1

Each node, vi, can be associated with a dfeat dimensional feature vector, Fi, similarly admitting a matrix representation as a (N × dfeat) matrix, F. These features can be derived from the underlying data or generated at run time in the case of a featureless graph. This makes it possible to define a general graph structure as2

Therefore, in practice, a complete specification of the graph structure can be achieved solely by using the matrices A and F.

In what follows, the terms clustering and partitioning are used interchangeably to refer to the GP problem. The general partitioning problem consists of dividing the nodes of G into M discrete nonoverlapping subsets, L ∈ {1, ..., M}, such that specific criteria are met. These objectives can be general and applicable to nearly all graphs and include, e.g., optimizing connectivity, minimizing intercluster edges, or maximizing intracluster similarities. Furthermore, these objectives can also be domain-specific, based on the underlying model the network represents, such as how the KC represents the dynamical information on the system. The (N × M) partitioning assignment matrix, S, assigns the nodes of G into M distinct states, such that3

Since the discrete nature of the assignment matrix poses difficulties for gradient-based methods, it is common to relax this constraint and work with continuous probabilities instead. In practice, this means introducing a soft assignment matrix instead, where SiL now denotes the probability of node vi belonging to cluster L, subject to the conservation of probability:

where L runs over all M possible clusters.

Some of the most commonly used clustering criteria used for the GP problem are the modularity, minimum cut, and Davies–Bouldin index (DBI), which are defined as follows:

2.1.1 Modularity

The modularity, Q, measures how the connectivity of graph partitioning differs from what would be expected in a purely random graph with the same node degree distribution. Consequently, this metric provides a quantitative assessment of how statistically surprising a proposed partitioning is, with higher values implying the presence of underlying communities within the data.37 The modularity is calculated using the following equation:4

where m is the total number of edges in the graph, di and dj are the degrees of nodes vi and vj, respectively, Aij is given by eq 1, and SiL and SjL are given by eq 3. Networks characterized by high modularity exhibit dense intracluster node connections but sparse intercluster node connections. While there have been criticisms of the underlying assumptions for modularity, it is a widely used and useful partitioning criterion in many areas.38−40

2.1.2 Minimum Cut

The minimum cut, C, of a given partitioning is the number of intercluster edges, which can be calculated using5

where all indices and matrix values are given as previously defined. The simplicity of this criterion has led to its popularity across multiple domains, spurred on by its proven equivalence to the max-flow criterion.41 It stands as one of the most studied partitioning criteria,42−44 alongside its further generalizations commonly used to obtain more balanced clusters with minimal cuts.45,46 Intuitively, the minimum cut partitioning of a network corresponds to the partitioning with the least amount of intercluster edges.

2.1.3 Davies–Bouldin Index

The DBI quantifies the average similarity between clusters,47 calculated based on the distance between clusters relative to the size of the clusters themselves, such that the DBI formula reads6

where L runs over all possible M clusters, ΔKL denotes the distance between cluster centroids, i.e., ΔKL = ∥ BK – BL ∥, and σL represents the cluster diameter, calculated as the average distance between the feature vectors and the centroid of the respective cluster:7

where Fi represents the feature vector of node i, and the centroid BL is computed by averaging the feature vectors of all N nodes i within cluster L. A lower DBI value (bounded from below by 0) indicates better partitioning, as it suggests higher intercluster dissimilarity in comparison to intracluster similarity. We employed the DBI score as a metric for assessing the clustering performance.

2.2 Markov Chains

The dynamical clustering criterion based on KC we use in this study is derived from the kinetic properties of an MC. In previous work, we explored the relationship between the MFPTs and KC, providing thorough derivations of all of the quantities presented below. Here, we briefly summarize the key results required and direct the reader to our previous paper for a detailed reference.48

In a discrete-time, homogeneous MC comprising n states, the probability vector of a random walk occupying any of the n-states at time t, p(t), can be expressed using a stochastic rate matrix, K, and a set of initial occupation probabilities, p(0), such that8

where t = l τ, since in the discrete case, transitions between states occur at integer multiples of a given time interval τ. The quantity Q(τ) = eKτ is known as the Markov matrix and is commonly used as the fundamental representation of a Markov chain. When this representation is used, time propagation corresponds simply to repeatedly multiplying Q by the current state of the system. To ensure the conservation of probability, we have9

The matrix K, or similarly Q, can be calculated from experimental or simulated transition data, such as MD trajectories or a series of string-like sequences typically used in genomics. K fully defines the dynamic properties of the MC and can be used to derive several useful quantities related to its underlying dynamics. Assuming detailed balance, the spectral decomposition of K into left, ψil, and right, ϕil, eigenvectors and corresponding eigenvalues, λl, is given by10

where the index l orders the eigenvalues in descending order. These eigenvalues are upper-bound by 0 and are guaranteed to be real due to the detailed balanced assumption. The eigenvector corresponding to the largest right eigenvector, ψ(1), is known as the equilibrium probability and is commonly denoted by peq, corresponding to the distribution the system converges to as t → ∞. The other eigenvectors contain useful kinetic information corresponding to the relaxation and mixing times of the MC. These eigenvectors can be used to derive an expression for the average time it takes for a random walk starting from state i to first arrive at state j, i.e., the MFPTs of an MC, denoted as tji:11

which can be alternatively expressed in the matrix form as12

KC is the surprising result that the sum over all states of the product of the MFPTs with the equilibrium probabilities, pjeq, is a constant value for all i, i.e.:13

where the fact that ∑j ϕj(l) ψj(l) = 1, ∀ l and ∑j ψj(l) = δl1 was used.48 Intuitively, therefore, the Kemeny constant is an important quantity to measure the overall time scales present in a dynamical system, as it is equal to the sum of all relaxation times, 1/|λl|.

By considering a clustering with the assignment matrix, S, from eq 3, we can define the equilibrium probability of cluster L, PLeq, as simply the sum of the populations of the individual nodes:14

KC can then be redefined in terms of a partitioned MC’s KC, ζS = ∑LPLeqt̂LK, where t̂LK is the coarse-grained mean first passage time from state K to L. We can define t̂LK via:4815

which leads to a variational form compared with the unclustered KC:2916

Rearranging this equation yields our partitioning metric ΔK:17

From this equation, it is clear that for the trivial single-state partitioning M = 1, ΔK = ζ holds, and when M = N, the value of ΔK is 0. As KC represents the sum of relaxation times in a system, and clustering speeds up relaxation times by removing intracluster relaxation times, maximizing the KC of a proposed partitioning aims to preserve the original time scales as closely as possible. From eq 17 (with the last equality requiring a crisp clustering), we can see that maximizing the KC is identical to minimizing ΔK, so this provides a useful metric of how kinetically disconnected the clusters are in a form suitable for gradient-based optimization methods like neural networks.

2.3 Graph Neural Networks

To describe GNNs in this work we adopt a message-passing formalism49 where a GNN-based encoder maps a node vi in G to a context-aware embedding using T separate layers via the iterative process:18

where t ∈ {1, ..., T} denotes the current network depth, is the local neighborhood of all nodes connected to node vi by edges eij, and hit represents the hidden representation of node vi at depth t. For consistency, the initial node feature vectors Fi are represented by hi0.

The functions Mt and Ut are commonly implemented as neural networks, known as the message-passing function and the vertex update function, respectively. Their exact formula depends on the type of graph layer used.

Two different graph layers are tested in our methodology: the GraphSAGE50 and the GATv251 layers. The former was chosen to provide a simpler, lower parameter function, which requires fewer resources to train, while the latter provides a more complex layer with the capability to better attune to all nodes. GraphSAGE layers use a simple message-passing function consisting of either the sum or mean, and a vertex update function that concatenates these expressions and feeds them into a feed-forward neural network:19

20

where [·||·] denotes concatenation and σ is any nonlinear function.

GATv2 layers replace this with a formulation based on scaled-dot-product attention seen in transformer layers. This differs by performing the aggregation step individually for each node in a neighborhood to calculate a context-aware attention mask, which is then used to perform a weighted vertex update:21

22

which can be generalized across multiple attention heads.

3 Methodology

The proposed method assumes that a given neural network architecture, Φ(θ), is sufficiently expressive such that the mappings it explores, from the space of dfeat dimensional node features, , to the space of possible M-state partitionings, , includes the optimal solution, Sopt, to the GP problem given our clustering criteria, along with a region of parameter space that can be optimized via gradient descent to yield it. Taking inspiration from Monte Carlo approaches to similar problems, we look to determine Sopt by initializing Φ with several different sets of parameters θ and optimizing them in parallel using modern ML tools. In this framework, GNNs are an effective Ansatz to begin partitioning graphs using complex criteria. They also provide a generic and intuitive optimization framework with gradient descent. As Φ is a general model, operating in without any problem-specific restrictions on its values or gradients, it can yield trivial solutions to eq 17, such as the aforementioned M = N case. To overcome this, we take inspiration from similar work on producing balanced clusters33,34 and added a penalty function to ΔK, which drives the network away from trivial solutions and gives the final form of the loss function, :23

Optimal KC partitionings can have unbalanced subsets, with clusters that contain only a few nodes often being the most interesting transition states.11 Therefore, the form of the penalty function is designed to penalize any invalid clustering (where any cluster has zero nodes), while also allowing the network to converge significantly unbalanced partitions with at least one node. The final penalty function is based on a Gaussian centered at 0, parametrized in this work using σ = 1.0 and η = 15. This means that the penalty increases as any cluster approaches 0 nodes and becomes negligible for clusters with 1 or more nodes.

As MCs do not fit into the graph structure outlined in eq 2, simple methods were derived for interpreting node connectivity using an adjacency matrix, A, and a node feature matrix, F. The connectivity of the graph was decided using a cutoff method, taking all node pairs (vi, vj) with (Qij + Qji)/2 > c, where c is a constant. In practice, we found c = 0.015 to be sufficient for all graphs we worked with, although in other applications, this could be changed depending on the density of the resulting network.

A systematic exploration of various node features was undertaken. We tested using the columns of various matrices that can be identified with nodes, such as the adjacency matrix, the rate matrix, the Markov matrix, and the MFPTs, as node features. For each node, we represented its connectivity, transition rates/probabilities to other nodes, or average time required to reach each node by splitting these matrices into vectors. Alongside these, principal component analysis (PCA), implemented as the singular value decomposition (SVD) of each of these matrices, of all previously mentioned matrices, and the eigenvectors of the rate or Markov matrix were also tested. Finally, we also tested both a simple linear layer, which combined the adjacency, rate/Markov, and MFPT columns for each node, and a trainable embedding layer, similar to those often used as dictionaries in large language models.52 More details regarding the features used are given in SI(A).

All neural network architectures tested used a GNN-based encoder, mapping the dfeat dimensional features into an embedding dimension dembed via a hidden representation of the size dhidden. Mean pooling was used for the aggregation function in the GraphSAGE layers as it provided better stability during training. The encoder was then followed by either a series of linear layers or dot-product attention-based transformer layers of dimension dhidden as a decoder.53 When using attention, the output of each transformer layer was normalized and added back to the input before being processed by a feed-forward layer. Our implementation used two transformer layers with a dropout layer between them. The output of a single trainable graph layer was used as positional encoding before the decoder. This transformer-like, attention-based decoder can process a whole cluster at once as a sequence of node embeddings. It was chosen to provide a comparison between a simpler model and one that is larger, more expressive, and better suited to handle highly contextual sequences such as the node-cluster assignment matrix. All neural networks used are described in detail in Figure 1.

Figure 1 Details of all four GNN architectures used, along with a description of the preprocessing required to begin training, together with an overview of the training loop process for a single gradient optimization step. Implementations of the GraphSAGE and GATv2 GNN layers used in the encoder are provided along with a schematic description of the transformer decoder (TD) and linear decoder (LD) architectures. Activation functions for the various layers are not included for the sake of brevity.

Building on the simplest model with GraphSAGE layers and the linear decoder (GraphSAGE + LD), we validated this model and extended it in two directions to test whether more complex architectures could improve the performance. Besides the GraphSAGE encoder, we also implemented the GATv2 encoder. In addition to the linear decoder, we also implemented a transformer decoder (TD). This results in four models: (i) GraphSAGE layers and the linear decoder (GraphSAGE + LD), (ii) linear decoder with GATv2 layers (GATv2 + LD), (iii) GraphSAGE with a transformer decoder (GraphSAGE + TD), and finally (iv) GATv2 with a transformer decoder (GraphSAGE + TD).

To enable fast and efficient training of the networks, a general parallel methodology for training GNNs for partitioning was established, as outlined in Figure 2. This requires three additional functions to be defined as well as a valid loss function and the number of processes P:

Figure 2 Schematic description of the algorithm used to train multiple networks in parallel, and the use of the general functions , , and .

3.1 Generator Function

: This function returns valid parameters θ for a given neural network Φ. This provides the initial starting point in parameter space for optimization, and complex generators can be used to provide a quicker convergence. In this work, however, simple initialization using both Kaiming and Glorot algorithms was tested. We used the Glorot method due to its slight advantage observed in early tests.54,55 Additionally, we explored bespoke initializations based on converged models, where a loss function was used that corresponded to PCCA+ clustering.

3.2 Culling Function

: This function determines whether or not a worker thread has finished training a network and is called after a predetermined number of optimization steps has finished. Observing individual runs, it is common to see neural networks that initialize in local minima of parameter space and are unable to be trained. An effective culling function allows for these types of runs to be terminated early, allowing for more time to be spent in interesting regions of the parameter space. To show the validity of the method, all experiments in this study terminate after a predetermined number of training steps.

3.3 Termination Function

: This function determines when to stop training and, when combined with effective choices for and , can allow for the implementation of complex methodologies for determining the convergence of an experiment.

This methodology aims to be as general as possible to make implementation of future applications for different criteria straightforward.

To test the validity and effectiveness of using GNNs to cluster Markov chains, the method was first tested on synthetic kinetic graphs. First, a simple stochastic block model (SBM)56 was employed to generate random networks with four neighborhoods. The initial set of networks used sampled intercluster connection probabilities from a Gaussian distribution scaled to be in the range [0, 0.25] and set intracluster probabilities to a constant of 0.5. This led to networks with dense clusters, where optimal partitioning always successfully separated them. We also designed SBM networks where the optimal clustering was less clear, sampling all connection probabilities from the same Gaussian as before. For both sets of examples, SBM networks with disconnected nodes were discarded, and the random-walk normalized Laplacian matrix, L, was used to generate rate matrices via the following equation:24

where I is the identity matrix, D is a matrix with the degrees of each node along the leading diagonal, and A is the adjacency matrix given by eq 1. Results on these synthetic graphs were compared with clusterings, and respective ΔK values were obtained from the PCCA+ clustering. We note that as PCCA+ uses a different kinetics-based objective function29 that reduces the clustering problem to M2 dimensions, the optimization is readily solvable using Schur decomposition, therefore only a single optimal cluster is obtained using PCCA+. Each SBM experiment trained 100 independently initialized models for 250 steps of gradient descent on 500 randomly generated graphs, ranging from N = 150 to N = 300 nodes in increments of 50.

To establish a baseline performance for all models on a single graph, the network was tested on a graph corresponding to a rate matrix derived from an analytical 1D potential simulation.29 The system consisted of 100 nodes distributed along the potential, which comprised four distinct wells: three with 20 nodes and one with 40 nodes. The optimal KC partitioning successfully separated all 4 wells into individual partitions. This example of a linear chain network architecture was chosen to test the network in situations where connectivity is limited. All four models were tested 1000 times on this example for 500 gradient descent steps each. A feature comparison was also performed on this system, where 1000 GraphSAGE + LD models were trained for 500 steps using each set of features described above.

Ultimately, the method was tested on a network derived from MD data. We used the MSM benchmark system of PyEMMA57 corresponding to 500 ns long implicit water MD simulations of a small 5-amino-acid system, here referred to as the pentapeptide system. The simulation initially underwent clustering into 250 states using K-means on tICA projected data, followed by the calculation of a Markov matrix,58,59 which served as the input to our method. Each experiment consisted of training 1000 models for 500 gradient descent steps to cluster the 250-state MSM into 5 states. Given the relevance of this example to biomolecular simulations, an additional hyperparameter search was conducted to identify the optimal network architecture. For this purpose, we again trained 1000 GraphSAGE + LD models for 500 steps while varying the number of layers in both the decoder and the encoder, along with the dimensions dhidden and dembed. Furthermore, we performed a feature comparison following the same methodology as previously described. We also compared the results of this example to the PCCA+ and parallel tempering variational clustering (PTVC), our previous MCMC approach.29 Additionally, we used the DBI score with the Markov matrix as the node features vector as a metric to assess the quality of the partitionings produced.

The neural networks and relevant code were implemented in Python 3.10.6 using PyTorch 2.060 or Tensorflow.61 Network diagrams were generated using the Netgraph library.62 To perform gradient descent steps, the Adam optimizer63 was employed with a learning rate of 1 × 10–4. All networks utilized dembed = 32 and dhidden = 64, with 4 GNN layers in the encoder and 2 layers in the decoder, unless otherwise specified. Due to the numerical sensitivity and dependence of the MFPT calculation on accurate eigenvalues, 64-bit floats were used for all calculations, which were run on the CPU. Each experiment involved the simultaneous execution of 90 processes in parallel with each process running on a single core. To compute the DBI scores, the scikit-learn implementation was utilized.64

4 Results and Discussion

Our GNN approach successfully found the optimal clustering for all SBM networks in the first set (e.g., Figure 3a), where the clusters were well-defined in 100 training runs or less. The optimal cluster matched the PCCA+ solution, confirming the viability of KC as the optimization criterion. Figure 3b,c details the performance on the second set of SBM graphs featuring a more diverse range of networks. Our implementation identified partitionings that better optimized ΔK than the PCCA+ (Figure S1). However, the improvement diminished as the number of nodes in the network decreased, as often the clustering became clearer. The stochastic nature of the GNN-based method, however, lends itself more to these systems where several viable partitionings have a low ΔK value. In these examples, the optimal clustering is close to many other similar solutions; therefore, it is helpful to evaluate many initializations to evaluate most of the very low ΔK partitionings. Violin plot distributions of optimal KC values found are also shown in Figure S1.

Figure 3 Performance of the SBM network clustering for randomly sampled cluster connection probabilities. (a) Example of partitioning produced by the GNN optimization. (b, c) Quantitative assessment of the method’s performance over 100 runs for the 500 different networks. (b) Proportion of runs where the method yields a partitioning with a lower ΔK value than PCCA+. (c) Average time needed in seconds for one training run when limited to run on one 2.20 GHz CPU core.

The 1D potential (Figure 4a) was used to test our approach on sparse graphs where connectivity is limited to a linear chain and the majority of effective ΔK partitionings lie in a small region of partitioning space. In this regime, gradient descent steps can become less effective as the majority of initial partitions are far from the region of interest and can be unable to escape their local minima. Fortunately, simply training more neural networks can alleviate the problem by providing a wider range of starting points for optimization. In future applications, an effective choice of the initialization function could speed up convergence, whereby can learn and then avoid initializing in uninteresting regions of parameter space. To further test whether different node features could improve the frequency of successful optimization outcomes, we tested several additional node features, sampling the outcome 1000 times each. In total, we tested nine different node feature types: (i) adjacency matrix, (ii) the PCA vectors created from the adjacency matrix, (iii) rate matrix, (iv) the PCA vectors created from the rate matrix, (v) MFPTs (default in other examples), (vi) the PCA vectors created from the MFPT matrix, (vii) eigenvectors of the rate matrix, (viii) trainable embeddings, and (ix) a neural network. The feature comparison results (Figure 4b) show that it is possible to use improved node features, such as the adjacency matrix or a trainable embedding in this case, to achieve a much more successful optimization for the 1D model potential-derived network.

Figure 4 (a) Plot showing the 1D potential used for an example of a linear chain network. Nodes are colored according to the most optimal 4-state clustering for the system. (b) Violin plots for feature comparisons on the 1D potential network. Each experiment was run 1000 times with each set of features. The corresponding distribution of final ΔK/KC values is shown, along with a dotted blue line signifying the ΔK/KC value of the optimal solution for each problem.

When using the rate matrix as the node feature in this example, the GNN fails to find the correct clustering even once. This may be due to the difficulty in discriminating nodes at high energies (with low populations) solely on the basis of the rate matrix, as these nodes have similar transitions. The adjacency matrix, however, provides more information about the node location in the global structure of the network, leading to a noticeable advantage. Our GNN implementation allows for simple changes like this to enable defining the node features on a problem-specific basis, therefore giving the method additional viability across a wide range of MCs that vary greatly in structure.

We also tested the four GNN architectures on the 1D potential, training 1000 models of each variety using the adjacency matrix as node features. Results are shown in Figure 5b and described in Table 1. We found that simpler models perform better at finding optimal partitionings over more complex models in this framework. Despite the better average performance the GATv2 and TD models provided, no models, but the simplest GNN was able to find the optimal clustering over 1000 runs.

Figure 5 (a) Clustering corresponding to the lowest ΔK/KC values found by all GNN models other than GraphSAGE + LD. (b) Distribution of ΔK/KC values from training 1000 copies of each neural network architecture on the 1D potential derived network. A logarithmic scale is used to highlight the large differences between the different models. The blue dashed line indicates the value associated with optimal clustering for this system.

Table 1 Network Architecture Comparison for the 1D Potential Derived Networka

model name	mean	minimum	global minimum accuracy	
GraphSAGE + LD	5884.74 (±8372.62)	169.69	5.6%	
GATv2 + LD	1988.50 (±2696.11)	806.48	0.0%	
GraphSAGE + TD	2201.59 (±1796.45)	806.48	0.0%	
GATv2 + TD	1508.19 (±2728.35)	806.48	0.0%	
PCCA+	169.69 (±0.0)	169.69	N/A	
a The four network architectures are compared in terms of the mean (standard deviation) and minimum for 1000 experiments along with the percentage of times the optimal solution was found.

The wide variance of the observed results showcases the poor ability of all models to optimize initial partitionings on this network when started from many regions of parameter space; however, in the smaller models, this is mitigated by the smaller solution space explored and the training time required. The large magnitude of some solutions is indicative of a partitioning problem where there is only a small region of interest in partitioning space, leading larger models like the GATv2 and TD models to struggle to find this region, possibly due to the more complex and unclear gradients calculated when the node embeddings are more codependent and additional resources are required to train the extra parameters. There is a sharp cutoff in the performance of all larger models, with the best value attained being consistently around 806 and the corresponding clustering shown in Figure 5a. We were unable to pinpoint the exact cause of this lower bound. We also note that in current MSM building strategies, only metastable clusters are identified, unless discretized collective variables are introduced. Typically, a k-mean clustering is performed that assigns very rarely sampled transition regions to metastable states. Therefore, the optimal clustering even from more complex GNNs (Figure 5a) is able to separate the three main metastable states, which is in line with current practical MSM applications. This problem is particularly challenging for GP objective functions that are not derived from kinetics-based measures, and such measures produce qualitatively incorrect optimal clustering, such as modularity.29 Accordingly, the nonoptimal solution of the GNNs represented in Figure 5a corresponds to a lower DBI value than that of Figure 4a, which suggests that the DBI measure is also inadequate for this system (Figure S2). The optimal ΔK clustering, however, is identical to the solution found by PCCA+ in this example and is consistently found by the PTVC method due to the small size of this network.

Examples derived from numerical simulations of biomolecular systems often result in noisy derived MCs with dense connections such as in the pentapeptide case. Here, the optimal PCCA+ clustering does not correspond to the optimal ΔK solution, with four nodes assigned to different clusters, which are highlighted in Figure 6a and contribute 7.4% of the overall equilibrium probabilities (Figure S3). As shown in Figure 7b and Table 2, again, the simplest GraphSAGE and linear decoder method outperforms all other GNN architectures, similarly to the 1D example above. The PTVC method does not identify the optimal ΔK clustering in general when using a random initialization but has been observed to find the correct solution when given the PCCA+ as an initial clustering. Here, we also compare the optimal solutions from different methods in Figure 7a by evaluating the DBI score. Interestingly, our optimal KC clustering also corresponds to the lowest DBI score, indicating a better clustering using this measure as well.

Figure 6 (a) Comparison of the PCCA+ solution with the optimal ΔK clustering for the pentapeptide system, which has DBI scores of 0.823 and 0.758, respectively. Networks are displayed using a force-directed layout weighted with the MFPTs. Nodes that differ between the two clusters are highlighted with a yellow border. (b) Violin plot distributions for 1000 experiments with each set of nine features. The corresponding distribution of final ΔK/KC values is shown as a violin plot, along with a blue dashed line for the optimal ΔK/KC value.

Figure 7 (a) Bar plot showing the DBI values associated with the optimal ΔK partitionings found by each method used. (b) Violin distribution plot of ΔK/K values acquired from training 1000 copies of each neural network architecture on the pentapeptide system. The blue dashed line indicates the value associated with optimal clustering.

Table 2 Network Architecture Comparison for the Pentapeptide Systema

model name	mean	minimum	global minimum accuracy	
GraphSAGE + LD	271.03 (±4.80)	258.67	2.8%	
GATv2 + LD	273.28 (±5.71)	262.81	0.0%	
GraphSAGE + TD	266.66 (±5.94)	259.098	0.0%	
GATv2 + TD	274.92 (±5.81)	262.74	0.0%	
PCCA+	258.86 (±0.0)	258.86	N/A	
a The four network architectures are compared in terms of the mean (standard deviation) and minimum for 1000 experiments, along with the percentage of times the optimal solution was found.

As the pentapeptide system is relatively small, significant structural differences cannot be identified within the small microstates. Correctly identifying these microstates is an open problem, with significant implications for the dynamics of the MSM.65 We expect more significant differences between states for larger networks that correspond to larger systems, going beyond prototype peptide simulations.

In this example, our method finds optimal clustering less frequently than in the synthetic examples of SBM networks. This could be due to the increased number of low ΔK partitionings available, which provide many local minima that are more difficult to escape. However, this is not a problem unique to the GNN method, as the PTVC method encounters the same issue, having difficulty finding the relevant areas of partitioning space, and is unable to find the optimal clustering in this example using random initialization. The dense and interconnected nature of the nodes means it is difficult for PCCA+-based methods to identify the most optimal clustering as well, as more optimal clusters may not differ much in terms of meta-stability. The larger models suffer from a vanishing gradient due to the heavily conditional nature of optimizing all node partition probabilities at the same time, resulting in an inability to find the optimal clustering even once while nevertheless providing a better result on average than the GraphSAGE and LD model. The best clusterings found by the GraphSAGE and TD models are shown in Figure S4 for reference. Still, they are outclassed by the PCCA+, and all but the GraphSAGE with linear decoder model provide worse performance than PTVC, highlighting the efficiency of simpler networks for this kind of task.

While further analysis of learning rate adjustments could provide an effective solution to train larger and more complex models on this problem, it is clear that training simple networks is both quicker and more reliable in finding the optimal clustering. Feature comparisons on the GraphSAGE + LD example show the success of the method in finding the optimal partitioning over 1000 individual runs, with all features tested except for using the eigenvectors (Figure 6b). MFPTs show the best performance here, finding the optimal solution the greatest number of times. Most likely, this is due to the computational simplicity in finding effective mappings from the MFPTs to optimized partitionings, as the value of eq 23 directly depends upon them. The poor performance of the eigenvectors is most likely down to numerical stability, as some of them are very small 64-bit floats.

One advantage of the PTVC method, which uses an MCMC-based algorithm, over our method is the ability to be initialized from any arbitrary partitioning. The results shown are for the PTVC when given no initial clustering information; however, when the method is initialized from the PCCA+ solution it consistently finds the optimal solution on all occasions. In the GNN models, we also changed our initialization function by pretraining the models to replicate the PCCA+ solution using mean-squared error loss that is minimal for the optimal PCCA+ solution. Once the GNN is trained to this specific loss function and identifies the PCCA+ clustering, it is then further trained using the Kemeny loss in eq 23. This pretraining, however, did not improve the performance of the GNN optimization (data not shown).

Figure 8 shows the results of a hyperparameter search using the best-performing simple GNN model to identify the optimal parameters for the network in the pentapeptide example. It highlights the problem that the network’s effectiveness is highly dependent upon the chosen hyper-parameters of the neural network and the size of the network. The best parameter space for a given problem must be sufficiently large to be able to express the wide multitude of possible S matrices while also being small enough to allow for effective optimization in a feasible time. This is highlighted by the poorer performance of most networks featuring four linear layers in the decoder, where most of the time the method either fails to find the optimal solution or does so rarely it is classified as an outlier. Future work can further verify whether a varying or learnable adjustment to the learning rate could provide a reasonable improvement, as this would allow networks to more efficiently probe the relevant regions of parameter space. Alternatively, further work could evaluate if it is possible to quantify the relationship between the size of the graph being clustered and the most optimal network size.

Figure 8 Parameter grid search on the pentapeptide example; 1000 networks were initialized for each set of parameters, and the resulting ΔK/KC distribution is shown. The gray background indicates the parameters identified as optimal.

5 Conclusions

We describe a GNN-based network clustering optimization method that uses the computing infrastructure and tooling that have risen around GNNs to provide a working basis for difficult optimization problems on graphs. The benefits are numerous as GNNs provide an ample Ansatz for beginning optimization for partitioning phenomena in situations where otherwise finding an efficient representation for the task proves to be difficult.

We provide a novel implementation of a GP approach to cluster graphs based on kinetic optimization criteria. Here, we use KC to obtain metastable kinetically optimal clusters. We outline a general methodology for training GNNs to optimize graph partitionings and compare the performance of four different architectures across three different kinetic network types using several different node features. Our results show that smaller and simpler networks often show superior performance in our test cases, as larger networks become more difficult to optimize and often cannot find the optimal solutions. However, as we are interested only in the best solution with the maximal Kemeny constant, multiple evaluations can lead to the correct global solution, even if only obtained less frequently in individual initializations. GraphSAGE layers combined with a linear decoder achieved both greater consistency and accuracy than any of the tested attention-based alternatives.

To show that our proposed method provides a viable alternative to existing methods in terms of the coarse-grained systems, we compare optimal clusterings obtained with those identified using the PCCA+ clustering. It is important to note, however, that as PCCA+ uses a different kinetics-based objective function, the task of identifying the optimal clustering using the KC-based measure is independent of this comparison. Testing the method on SBM graphs generated to have a clear optimal partitioning, we can attain results similar to those of PCCA+, but with more optimal values of the Kemeny constant. A linear chain corresponding to an analytic 1D free energy profile also poses a difficult task for GNNs due to low connectivity, while PCCA+ in this case can produce the optimal solution. We show how we can adapt the node features the GNN uses and the number of runs to alleviate this problem and improve the performance. However, optimal clustering is only found in about 5% of the attempts in only one GNN architecture, and therefore, improving the reliability of the method still requires future work. The outline of our general methodology, introducing the strategy-defining generator , culling , and termination functions, provides an ample framework for future work in this direction.

Our MD-simulation-derived MSM for a pentapeptide is the most complex example we considered here. It has a slightly more optimal ΔK clustering that seems to better capture metastable states than the PCCA+-derived clusters, as also seen by comparing the corresponding DBI scores. Only one of the four architectures can find this optimal solution in less than 3% of the time, highlighting the problem of GNNs, which further work on functions that can hopefully provide a solution. Additionally, GNNs do not provide a straightforward way to initialize the search from a close-to-optimal solution, which makes the PTVC algorithm more efficient.

In this work, we explored a number of node features that were used by the GNNs to enable efficient optimization. The MFPT-based features were some of the most successful ones, as these directly relate to the definition of ΔK. However, for examples derived from biomolecular simulations, features that are specifically relevant to physical and structural observables (angles, distances, energy terms, etc.) could outperform the general features used here. This is a promising avenue for future directions.

In summary, we have shown the validity of using GNNs as an effective method to solve the GP problem applied to MCs, alongside describing a general framework outlining the strategy for optimizing partitionings using GNNs in parallel on graph-like structures. We propose using such a framework to solve the many different and difficult nonconvex optimization tasks that are commonplace in computational fields using graph structures, taking full advantage of the speed, availability, and efficiency of modern ML tooling and associated hardware. The general training paradigm outlined provides a working basis for generalizing this strategy to different problems, and the results highlight the usefulness of simpler GNN architectures when solving such optimization problems. It remains to be seen how this can be improved by implementing more complex parameter space sampling strategies to aid in the identification of more optimal solutions, taking inspiration from similar techniques used for enhanced sampling in MD.

Supporting Information Available

The Supporting Information is available free of charge at https://pubs.acs.org/doi/10.1021/acs.jpcb.3c08213.Comparison of the ΔK values of the partitionings found by PCCA+ (the dotted black line) and the GraphSAGE + LD GNN method on four SBM networks; DBI values corresponding to clusterings found for the 1D example; histogram of the equilibrium probability distributions of each node in the pentapeptide example; and comparison of the pentapeptide network partitionings with the lowest ΔK value found by each method used (PDF)

Supplementary Material

jp3c08213_si_001.pdf

The authors declare no competing financial interest.

Acknowledgments

We thank Dr. Alessia Annibale for helpful suggestions, Madhav Sharma for his feedback on the manuscript, and Bin Liang for his initial contributions to pretraining. The authors acknowledge funding from the ERC Starting Grant (Project 757850 BioNet) and EPSRC (grant no. EP/R013012/1). This project made use of time on ARCHER granted via the UK High-End Computing Consortium for Biomolecular Simulation, HECBioSim (http://hecbiosim.ac.uk). The authors acknowledge the use of the UCL Myriad High Performance Computing Facility (Myriad@UCL), and associated support services, in the completion of this work.
==== Refs
References

Okamoto H. ; Qiu X. Detecting hierarchical organization of pervasive communities by modular decomposition of Markov chain. Sci. Rep. 2022, 12 , 20211 10.1038/s41598-022-24567-x.36418410
Von Hilgers P. ; Langville A. N. The Five Greatest Applications of Markov chains. In Proceedings of the Markov Anniversary meeting, 2006; pp 155–158.
Sanderson C. ; Guenter S. Short text authorship attribution via sequence kernels, Markov chains and author unmasking: An investigation. In Proceedings of the 2006 Conference on Empirical Methods in Natural Language Processing, 2006; pp 482–491.
Schbath S. An overview on the distribution of word counts in Markov chains. Journal of Computational Biology 2000, 7 , 193–201. 10.1089/10665270050081469.10890396
Bianchini M. ; Gori M. ; Scarselli F. Inside pagerank. ACM Transactions on Internet Technology (TOIT) 2005, 5 , 92–128. 10.1145/1052934.1052938.
Berthiaux H. ; Mizonov V. Applications of Markov chains in particulate process engineering: a review. Canadian journal of chemical engineering 2004, 82 , 1143–1168. 10.1002/cjce.5450820602.
Shukla D. ; Hernández C. X. ; Weber J. K. ; Pande V. S. Markov State Models Provide Insights into Dynamic Modulation of Protein Function. Acc. Chem. Res. 2015, 48 , 414–422. 10.1021/ar5002999.25625937
PMID: 25625937

Swope W. C. ; Pitera J. W. ; Suits F. Describing Protein Folding Kinetics by Molecular Dynamics Simulations. 1. Theory. J. Phys. Chem. B 2004, 108 , 6571–6581. 10.1021/jp037421y.
Singhal N. ; Snow C. D. ; Pande V. S. Using path sampling to build better Markovian state models: Predicting the folding rate and mechanism of a tryptophan zipper beta hairpin. J. Chem. Phys. 2004, 121 , 415–425. 10.1063/1.1738647.15260562
Pan A. C. ; Roux B. Building Markov state models along pathways to determine free energies and rates of transitions. J. Chem. Phys. 2008, 129 , 064107 10.1063/1.2959573.18715051
Martini L. ; Kells A. ; Covino R. ; Hummer G. ; Buchete N.-V. ; Rosta E. Variational Identification of Markovian Transition States. Phys. Rev. X 2017, 7 , 031060 10.1103/PhysRevX.7.031060.
Lane T. J. ; Bowman G. R. ; Beauchamp K. ; Voelz V. A. ; Pande V. S. Markov state model reveals folding and functional dynamics in ultra-long MD trajectories. J. Am. Chem. Soc. 2011, 133 , 18413–18419. 10.1021/ja207470h.21988563
Bowman G. R. ; Bolin E. R. ; Hart K. M. ; Maguire B. C. ; Marqusee S. Discovery of multiple hidden allosteric sites by combining Markov state models and experiments. Proc. Natl. Acad. Sci. U. S. A. 2015, 112 , 2734–2739. 10.1073/pnas.1417811112.25730859
Mardt A. ; Pasquali L. ; Wu H. ; Noé F. Author Correction: VAMPnets for deep learning of molecular kinetics. Nat. Commun. 2018, 9 , 4443 10.1038/s41467-018-06999-0.30349135
Jo S. ; Suh D. ; He Z. ; Chipot C. ; Roux B. Leveraging the information from Markov state models to improve the convergence of umbrella sampling simulations. J. Phys. Chem. B 2016, 120 , 8733–8742. 10.1021/acs.jpcb.6b05125.27409349
Husic B. E. ; Pande V. S. Markov State Models: From an Art to a Science. J. Am. Chem. Soc. 2018, 140 , 2386–2396. 10.1021/jacs.7b12191.29323881
PMID: 29323881

Deuflhard P. ; Weber M. Robust Perron cluster analysis in conformation dynamics. Linear algebra and its applications 2005, 398 , 161–184. 10.1016/j.laa.2004.10.026.
Hummer G. ; Szabo A. Optimal dimensionality reduction of multistate kinetic and Markov-state models. J. Phys. Chem. B 2015, 119 , 9029–9037. 10.1021/jp508375q.25296279
Buchete N.-V. ; Hummer G. Coarse master equations for peptide folding dynamics. J. Phys. Chem. B 2008, 112 , 6057–6069. 10.1021/jp0761665.18232681
Reuter B. ; Fackeldey K. ; Weber M. Generalized Markov modeling of nonreversible molecular kinetics. J. Chem. Phys. 2019, 150 , 174103 10.1063/1.5064530.31067901
Peng Y. ; Pak A. J. ; Durumeric A. E. ; Sahrmann P. G. ; Mani S. ; Jin J. ; Loose T. D. ; Beiter J. ; Voth G. A. OpenMSCG: A Software Tool for Bottom-Up Coarse-Graining. J. Phys. Chem. B 2023, 127 , 8537–8550. 10.1021/acs.jpcb.3c04473.37791670
Wang J. ; Olsson S. ; Wehmeyer C. ; Pérez A. ; Charron N. E. ; De Fabritiis G. ; Noé F. ; Clementi C. Machine learning of coarse-grained molecular dynamics force fields. ACS central science 2019, 5 , 755–767. 10.1021/acscentsci.8b00913.31139712
Li Y. ; Dong Z. Effect of clustering algorithm on establishing Markov state model for molecular dynamics simulations. J. Chem. Inf. Model. 2016, 56 , 1205–1215. 10.1021/acs.jcim.6b00181.27249546
Nagel D. ; Sartore S. ; Stock G. Toward a Benchmark for Markov State Models: The Folding of HP35. J. Phys. Chem. Lett. 2023, 14 , 6956–6967. 10.1021/acs.jpclett.3c01561.37504674
Röblitz S. ; Weber M. Fuzzy spectral clustering by PCCA+: application to Markov state models and data classification. Advances in Data Analysis and Classification 2013, 7 , 147–179. 10.1007/s11634-013-0134-6.
Zhu B. Zhu Graph partition by Swendsen-Wang cuts. In Proceedings Ninth IEEE International Conference on Computer Vision, 2003; pp 320–327.
Anderson D. F. ; Higham D. J. Multilevel Monte Carlo for continuous time Markov chains, with applications in biochemical kinetics. Multiscale Modeling & Simulation 2012, 10 , 146–179. 10.1137/110840546.
Luo Z. T. ; Sang H. ; Mallick B. A Bayesian contiguous partitioning method for learning clustered latent variables. J. Mach. Learn. Res. 2021, 22 , 1748–1799.
Koskin V. ; Kells A. ; Clayton J. ; Hartmann A. K. ; Annibale A. ; Rosta E. Variational kinetic clustering of complex networks. J. Chem. Phys. 2023, 158 , 104112 10.1063/5.0105099.36922127
Shaham U. ; Stanton K. ; Li H. ; Nadler B. ; Basri R. ; Kluger Y. SpectralNet: Spectral clustering using deep neural networks. In 6th International Conference on Learning Representations, ICLR 2018, 2018.
Kawamoto T. ; Tsubaki M. ; Obuchi T. Mean-field theory of graph neural networks in graph partitioning. J. Stat. Mech.: Theory Exp. 2018, 31 , 124007 10.1088/1742-5468/ab3456.
Liao R. ; Brockschmidt M. ; Tarlow D. ; Gaunt A. L. ; Urtasun R. ; Zemel R. Graph partition neural networks for semi-supervised classification, 2018; arXiv preprint arXiv:1803.06272 (accessed December 10, 2023).
Nazi A. ; Hang W. ; Goldie A. ; Ravi S. ; Mirhoseini A. A Deep Learning Framework for Graph Partitioning; ICLR, 2019.
Tsitsulin A. ; Palowitch J. ; Perozzi B. ; Müller E. Graph clustering with graph neural networks. J. Mach. Learn. Res. 2023, 24 , 1–21.
Singhal N. ; Pande V. S. Error analysis and efficient sampling in Markovian state models for molecular dynamics. J. Chem. Phys. 2005, 123 , 204909 10.1063/1.2116947.16351319
Kube S. ; Weber M. A coarse graining method for the identification of transition rates between molecular conformations. J. Chem. Phys. 2007, 126 , 024103 10.1063/1.2404953.17228939
Newman M. E. Modularity and community structure in networks. Proc. Natl. Acad. Sci. U. S. A. 2006, 103 , 8577–8582. 10.1073/pnas.0601602103.16723398
Guimera R. ; Sales-Pardo M. ; Amaral L. A. N. Modularity from fluctuations in random graphs and complex networks. Phys. Rev. E 2004, 70 , 025101 10.1103/PhysRevE.70.025101.
Good B. H. ; De Montjoye Y.-A. ; Clauset A. Performance of modularity maximization in practical contexts. Phys. Rev. E 2010, 81 , 046106 10.1103/PhysRevE.81.046106.
Agarwal G. ; Kempe D. Modularity-maximizing graph communities via mathematical programming. European Physical Journal B 2008, 66 , 409–418. 10.1140/epjb/e2008-00425-1.
Leighton T. ; Rao S. Multicommodity max-flow min-cut theorems and their use in designing approximation algorithms. Journal of the ACM (JACM) 1999, 46 , 787–832. 10.1145/331524.331526.
Ding C. H. ; He X. ; Zha H. ; Gu M. ; Simon H. D. A min-max cut algorithm for graph partitioning and data clustering. In Proceedings 2001 IEEE international conference on data mining, 2001; pp 107–114.
Hagen L. ; Kahng A. Fast spectral methods for ratio cut partitioning and clustering. In 1991 IEEE international conference on computer-aided design digest of technical papers, 1991; pp 10–11.
Even G. ; Naor J. ; Rao S. ; Schieber B. Fast approximate graph partitioning algorithms. SIAM Journal on Computing 1999, 28 , 2187–2214. 10.1137/S0097539796308217.
Andreev K. ; Räcke H. Balanced graph partitioning. In Proceedings of the sixteenth annual ACM symposium on Parallelism in algorithms and architectures, 2004; pp 120–124.
Krauthgamer R. ; Naor J. ; Schwartz R. Partitioning graphs into balanced components. In Proceedings of the twentieth annual ACM-SIAM symposium on Discrete algorithms, 2009; pp 942–949.
Davies D. L. ; Bouldin D. W. A Cluster Separation Measure. IEEE Transactions on Pattern Analysis and Machine Intelligence 1979, PAMI-1 , 224–227. 10.1109/TPAMI.1979.4766909.
Kells A. ; Koskin V. ; Rosta E. ; Annibale A. Correlation functions, mean first passage times, and the Kemeny constant. J. Chem. Phys. 2020, 152 , 104108 10.1063/1.5143504.32171226
Gilmer J. ; Schoenholz S. S. ; Riley P. F. ; Vinyals O. ; Dahl G. E. Neural message passing for quantum chemistry. Int. Conf. Mach. Learn. 2017, 1263–1272.
Hamilton W. ; Ying Z. ; Leskovec J. Inductive representation learning on large graphs. In Advances in neural information processing systems, 2017; Vol. 30 .
Brody S. ; Alon U. ; Yahav E. How Attentive are Graph Attention Networks? In International Conference on Learning Representations, 2021.
Mikolov T. ; Sutskever I. ; Chen K. ; Corrado G. S. ; Dean J. .Distributed representations of words and phrases and their compositionality. In Advances in neural information processing systems, 2013; Vol. 26 .
Klambauer G. ; Unterthiner T. ; Mayr A. ; Hochreiter S. Self-normalizing neural networks. In Advances in neural information processing systems, 2017; Vol. 30 .
He K. ; Zhang X. ; Ren S. ; Sun J. Delving deep into rectifiers: Surpassing human-level performance on imagenet classification. In Proceedings of the IEEE international conference on computer vision, 2015; pp 1026–1034.
Glorot X. ; Bengio Y. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, 2010; pp 249–256.
Abbe E. Community detection and stochastic block models: recent developments. J. Mach. Learn. Res. 2018, 18 , 1–86. 10.1561/9781680834772.
Scherer M. K. ; Trendelkamp-Schroer B. ; Paul F. ; Pérez-Hernández G. ; Hoffmann M. ; Plattner N. ; Wehmeyer C. ; Prinz J.-H. ; Noé F. PyEMMA 2: A software package for estimation, validation, and analysis of Markov models. J. Chem. Theory Comput. 2015, 11 , 5525–5542. 10.1021/acs.jctc.5b00743.26574340
Schwantes C. R. ; Pande V. S. Improvements in Markov state model construction reveal many non-native interactions in the folding of NTL9. J. Chem. Theory Comput. 2013, 9 , 2000–2009. 10.1021/ct300878a.23750122
Pérez-Hernández G. ; Paul F. ; Giorgino T. ; De Fabritiis G. ; Noé F. Identification of slow molecular order parameters for Markov model construction. J. Chem. Phys. 2013, 139 , 015102 10.1063/1.4811489.23822324
Paszke A. ; Gross S. ; Massa F. ; Lerer A. ; Bradbury J. ; Chanan G. ; Killeen T. ; Lin Z. ; Gimelshein N. ; Antiga L. Pytorch: An imperative style, high-performance deep learning library. In Advances in neural information processing systems, 2019; Vol. 32 .
Abadi M. TensorFlow: Large-Scale Machine Learning on Heterogeneous Systems, 2015; https://www.tensorflow.org/, Software available from tensorflow.org (accessed April 3, 2024).
Brodersen P. J. N. Netgraph: Publication-quality Network Visualisations in Python. Journal of Open Source Software 2023, 8 , 5372 10.21105/joss.05372.
Kingma D. P. ; Ba J. Adam: A Method for Stochastic Optimization. In ICLR (Poster), 2015.
Pedregosa F. ; Varoquaux G. ; Gramfort A. ; Michel V. ; Thirion B. ; Grisel O. ; Blondel M. ; Prettenhofer P. ; Weiss R. ; Dubourg V. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 2011, 12 , 2825–2830.
Arbon R. E. ; Zhu Y. ; Mey A. S. Markov state models: to optimize or not to optimize. J. Chem. Theory Comput. 2024, 20 , 977–988. 10.1021/acs.jctc.3c01134.38163961
