Protein folding intermediates on the dimensionality reduced landscape with UMAP and native contact likelihood

In molecular dynamics (MD) simulations of proteins, it is essential to investigate the conformational landscapes with a high dimensionality, which include all information on protein stability, fluctuation, and transitions between different states.¹ We usually focus on thermodynamics and/or kinetics of proteins on the landscapes spanned with a few collective variables (CVs).² Since a protein consisting of N atoms has 3N − 6 internal degrees of freedom, a large reduction in the landscape dimensionality is required. Principal component analysis (PCA)^3–5 or time-lagged independent component analysis (tICA)⁶ have been often used to reduce the dimensionality of conformational landscapes and describe low-frequency modes of globular protein dynamics.^4,5,7,8

Due to an increase of computational power and advances of simulation algorithms, hundreds of μs or longer simulations have become more feasible.⁹ Then, MD simulations have allowed us to investigate not only small fluctuations around experimental structures but also much larger transitions.¹⁰ Folding/unfolding dynamics of small proteins have also been examined using atomistic or coarse-grained (CG) MD simulations.^11–13 For these long MD trajectories, PCA or tICA, which are categorized as “linear” dimensionality reduction methods, sometimes failed to describe essential dynamics of proteins.^14,15 In our recent MD simulations of spike protein on the surface of SARS-CoV-2, PCA with Cartesian coordinates could not visualize up-down motions of the receptor binding domain clearly.¹⁶ It was mainly because the low-frequency motions from PCA emphasized functionally less important loop fluctuations. Therefore, we used the CG beads, each of which represents a domain motion in spike protein as feature variable of PCA¹⁷ and distinguished up, down, and other intermediate structures on the dimensionality reduced landscape.

One of the simplest approach to overcome the problems in conventional PCA or tICA is applying “non-linear” dimensionality reduction methods, such as ISOMAP,^18,19 diffusion map,^20,21 sketch-map,^22,23 t-distributed stochastic neighbor embedding (t-SNE),^15,24 or uniform manifold approximation and projection (UMAP).^25,26 Based on various criteria for data point proximity, these methods map objects from high-dimensional space to lower dimensional one, without drastically changing distances between data points.²⁴ ISOMAP and diffusion map define a distance of two data points by computing the shortest path length that connects in-between data points,¹⁸ whereas diffusion map focuses on distance-dependent transition probability rather than the distance itself.²⁰ Sketch-map aims to reproduce the dissimilarity between any two data points in the high-dimensional space.²² t-SNE computes a weight for two data points based on the “direct” distance independent of in-between data points.²⁴ The strategy of UMAP is similar to t-SNE, but UMAP only focuses on the pairs of data points that are defined as “neighbors” in a high-dimensional space.²⁵ While UMAP shows similar or better performances compared to t-SNE in several applications,²⁷ the computational cost of UMAP is much smaller than t-SNE. UMAP were already applied to protein MD trajectories by Peng Tao’s group.^15,26 They performed atomistic MD simulations of several globular proteins for about 1 µs and examined thermodynamics on the conformational landscapes and kinetics with Markov State Modeling (MSM) after dimensionality reductions by PCA, tICA, t-SNE, and UMAP. They concluded that both t-SNE and UMAP outperform PCA and tICA in the structural classification on the landscapes and MSM for better description of thermodynamics and kinetics of proteins.

In both linear (PCA or tICA) and non-linear (t-SNE or UMAP) dimensionality reduction methods, the selection of feature variables to describe a macromolecular structure is essential.² Cartesian coordinates or dihedral angles of a target system are conventionally used in PCA or tICA. However, there are two well-known problems in analyses of protein folding simulations, when we use Cartesian coordinates in PCA. The first one is that the average conformation may lose its proper chemical structures and then low-frequency motions from PCA can be meaningless. The other one is that various disordered structures away from the native one are difficult to be distinguished in the Cartesian PCA. Gerhard Stock’s group used inter-residue contact distances as feature variables for the folding of HP35 to overcome the second problem.^28,29 They could not clearly distinguish the native and intermediate states using the feature variables, suggesting to use backbone dihedral angles.

In the current study, we investigate the applicability of UMAP for describing folding/unfolding dynamics of a small globular protein, B1 domain of protein G (hereafter, we call it protein G) (Fig. 1). Due to slow folding/unfolding dynamics, a CG model of the protein is used in our MD simulation. One of the key differences between Peng Tao’s works and ours is the feature variables in UMAP to describe dimensionality reduced landscapes. Here, we propose native contact likelihood as feature variables and examine the performances of PCA, tICA, and UMAP. This work demonstrates usefulness of UMAP together with proper feature variables in thermodynamic analysis on the dimensionality reduced landscape and kinetic analysis with Markov state modeling (MSM).^30–34 This approach would be applicable to large-scale conformational transitions of multi-domain proteins and other biomolecular complexes.

Here, we briefly describe theoretical framework of the non-linear dimensionality reduction by UMAP.²⁵ UMAP consists of the following three steps: (i) definition of the “neighborhoods” of each data point, (ii) construction of a fuzzy topological structure, and (iii) a nonlinear projection from original space to dimensionality reduced one, preserving the data-point proximity. In linear methods, such as PCA, we can evaluate the reduced dimensionality in terms of contributions of selected components. In contrast, we cannot use such a criterion in UMAP, and the dimension of the reduced space in UMAP is arbitrary. In this study, we set it to 2, in order to understand the reduced conformational landscape intuitively.

In Fig. 2, we assume N data points (⁠$xi,i=1,2,…,N$ in original space. Each data point corresponds to a snapshot structure in MD simulations. In UMAP, the number of neighborhood points (NN) from each point is fixed as a hyperparameter. Neighborhood points of x_i for NN = 2 are x_j and x_k, whose distances from x_i are the two smallest ones. A topological network in original space is then determined by connecting all the neighborhood points (Fig. 2, left top). To preserve relative data positions in the topological network before and after nonlinear projection, we utilize a fuzzy topological structure characterized with a set of weights $pij|i,j=1,2,…$ of the connection between x_i and x_j. This value, p_ij, is the largest (=1) if the distance between x_i and x_j is the smallest in all the neighborhood points. The weight becomes smaller if the distance is greater. The fuzzy topological structure in dimensionality reduced space is optimized to mimic that in original space as much as possible (Fig. 2, right top). Mathematically, UMAP minimizes a cross entropy H defined as

by changing the coordinates in the dimensionality reduced space $yi|i=1,2,…$⁠.

The default value of NN in UMAP is 15. If we select a larger value of NN, the global manifold in the original space is well preserved in the dimensionality reduced one, while the computational costs increase. In contrast, a smaller value of NN gives us a better description of local structures in the topological network, while leading to a failure to preserve the global manifold. Practically, we have to test several values of NN in each application.

In this study, UMAP was applied to reduce the dimensionality of conformational landscape of protein G. UMAP calculations were performed by using a Python library implemented by the original developers.³⁵ For comparison, PCA and tICA were carried out using python libraries in PyEMMA 2, a useful library for the analysis of MD trajectories.³⁶ In addition, we prepared and used wrapper scripts for UMAP and PyEMMA 2 to facilitate UMAP analysis of MD trajectories. We set the lag time in tICA to 2.5 ns.

In the current study, we use quantities derived from native contact pairs in the folded structure. It is defined as

for one native contact pair by atoms i, j. d_ij,0 is the native contact distance of this pair, and d_ij is a distance of the same atom pair for a given structure. σ_ij is the variance of the native contact distance, and ideally, this parameter should reflect an interaction strength of corresponding contact pair. The quantity defined in Eq. (2) is able to be interpreted as a likelihood of the contact formation between atoms i and j, given a pairwise distance d_ij. Therefore, we call this quantity by “contact likelihood.” The average contact likelihood for protein-G is designated to be equivalent to q-value.³⁷ This contact likelihood becomes 1 when the native contact is formed and becomes 0 when atoms i, j are far away. The unfolded structures distribute around the origin $0,0,0,…,0T$ in contact likelihood space, while the folded and native-like structures distribute around $1,1,1,…,1T$⁠. Intermediate structures are expected to distribute between the unfolded and folded structures.

In this study, σ_ij in Eq. (2) is determined using the crystallographic B-factors of protein G. There is a relationship between the B-factor and the mean square displacement ⟨u²⟩ of the atom that B = 8π²⟨u²⟩. We can calculate σ_ij for the Cα-based CG model using the reproductive property of Gaussian

where $Bni$ represents the B-factor of n_ith atom in the ith residue and N_i is the number of atoms in the residue. Details of the determination of σ_ij values from the B-factor of protein G is explained in the supplementary material, Fig. S1.

We performed CG MD simulation of protein G, which contains 56 amino-acid residues. It consists of N-terminal β-hairpin by its two β-strands (β₁₂), one α-helix, and C-terminal β-hairpin by the other two strands (β₃₄) (Fig. 1). β₁₂ and β₃₄ contact with β₁ and β₄ strands, and consequently these four strands form an anti-parallel β-sheet. The KB-Go model provided to MMTSB web service was applied to protein G.^38,39 This model takes into account of only Cα atoms in a protein. The crystal structure [Protein Data Bank (PDB) ID: 1PGB⁴⁰] was used to set up the KB-Go model and to define σ_ij in the native contact likelihood as defined in Eq. (2). The number of native contact pairs in protein G (1PGB) is 84.

Both the native and non-native contact interactions were calculated without any truncation. 5.4 µs MD simulation with the KB-Go model was performed using the GENESIS program.^41,42 In the simulation, system temperature was kept constant at 330 K in NVT ensemble with Langevin thermostat⁴³ of which friction coefficient was 0.01 ps⁻¹. The time step was set to 20 fs with SHAKE constraints⁴⁴ for all bonds, and the simulation snapshot was saved in a trajectory at every 5 ps.

To examine dynamical information from MD simulations, MSM has become one of the most popular methods.^45,46 Suppose that structures in given simulation trajectories are discretized into n states with a clustering method. The construction of MSM is equivalent to calculate a n × n transition probability matrix $Pτ$⁠. The element, $Pijτ$ (i, j = 1, 2, …, n), is a transition probability defined as

with a certain lag time τ. x denotes a conformation of the target system, and n is the number of states. i and j are the state indices. The time evolution of a probability vector, $pt=p1t,p2t,…,pntT$⁠, each of which gives a probability that the system visits a state at time t, can be written using the matrix, such as

If we assume thermodynamic equilibrium, which imposes a symmetrical transition matrix, we obtain more detailed dynamical information by solving the eigenvalue problem of $Pτ$⁠. The first eigenvalue should be 1, and its eigenvector represents the equilibrium distribution. All the other eigenvectors correspond to certain transition modes. The magnitude of each eigenvalue indicates the contribution of each mode to the conformational transition. Then, plotting eigenvectors into space on which MSM is constructed gives us an intuitive view for the conformational transitions.

Practically, the state discretization or clustering is carried out using the essential feature variables or slow coordinates instead of Cartesian coordinates of biomolecules. Linear (PCA or tICA) or non-linear (UMAP and others) dimensionality reduction methods are applied to get the slow coordinates. Frequently, tICA is used to obtain a set of slow coordinates.^47–49 In this study, we compare the performance of MSM using UMAP with those using PCA and tICA.

To determine the number of discretized states n in MSM, we compute a generalized matrix Rayleigh quotient (GMRQ),^50,51 defined as

Here, v_i is the ith eigenvector of $Pτ$⁠, and T denotes the transposition of a vector or a matrix. The computation of GMRQ requires two independent datasets: one for training data and another for validation. For this purpose, we randomly split MD trajectories into the two sets. If we set the number of states too low, the GMRQ value becomes incorrectly small due to discretization errors. Increasing the number of states can improve the GMRQ value, while discretization with too many states declines the GMRQ value, resulting from the overfitting to the training MSM. In this study, state discretization was done by k-means clustering method,⁵² and all processes accompanying MSM analysis, including k-means clustering, were carried out using PyEMMA 2.³⁶ 

We performed 5.4 µs CG MD simulation of protein G. Total number of structures sampled in the MD simulation reached more than one million. In the trajectories, protein G frequently stayed at the folded or unfolded states for several hundred ns. In the current analysis, we reduced the number of snapshots in MD simulation trajectories. First, we computed root mean square deviations (RMSD) of all snapshots against the native structure defined as

Here, N is the number of Cα atoms. r_i and $ri0$ are coordinates of the ith Cα atom in MD snapshot and the native structure, respectively. Based on the RMSD values, we define the folded and unfolded states of protein G if the RMSD values are <5 and >10 Å, respectively. Then, we divided whole trajectories into 27 subsets, each of which includes 40 000 structures. We extracted only 10 subsets that involve the folding/unfolding transitions (totally 400 000 structures) in terms of the RMSD changes. By discarding the other 17 non-transition trajectories, the landscape in the current analysis does not show correct free-energy changes while it keeps the global shape of the folding/unfolding transition landscape. Hereafter, we call it as just a “conformational landscape.” Figure 3 shows a time course of RMSD value in the concatenation of extracted trajectories, and folding/unfolding transitions occurred more than twenty times. The RMSD values of protein G are ranging from 0 to 30 Å, indicating that various structures are sampled in our CG MD simulation.

We compared several feature variables for the dimensionality reduction of conformational landscape. In Fig. 4, we plot histograms of the RMSD [Fig. 4(a)] and averaged contact likelihood [Fig. 4(b)]. The folded and unfolded states are well separated in both cases, but the unfolded state shows different distributions in these two histograms. In the RMSD histogram, the unfolded state is widely distributed compared to the folded one, suggesting that conformational heterogeneity of the unfolded state observed in the CG MD simulations. In the averaged contact likelihood histogram, the folded and unfolded states show similar deviations from their respective averages, reflecting that this feature variable categorizes some of different unfolded structures into one class. As a result, the unfolded state has a compact distribution in the contact likelihood space. To examine the relationship between this feature variable and conformational classifications, we calculate the silhouette score for the folded and unfolded states.⁵³ As mentioned, the folded and unfolded states are defined in terms of RMSD values. Silhouette score s_i for the ith MD structure x_i is calculated with these formulas,

C₁ is a state to which x_i belongs, and C₂ is the other. N₁ and N₂ are numbers of data points in C₁ and C₂, respectively. ‖⋯‖ denotes a norm. a_i is a measure of the compactness of C₁ in data space, and b_i indicates how well x_i is separated from C₂. s_i can take on any value between −1 < s_i < 1, but should be larger for better separation. If s_i takes a negative value, the classification may be inappropriate. We average silhouette scores for the folded and unfolded states to obtain s_fold and s_unfold, respectively.

The silhouette scores using Cartesian coordinates, torsion angles, and contact likelihood as feature variables are compared in Table I. s_fold with Cartesian coordinates has a high value close to 1, whereas s_unfold is almost zero. This means that unfolded structures distribute broadly around one peak for the folded state, as clearly visualized in the supplementary material, Fig. S2a. s_fold and s_unfold with torsion angles as feature variables show a similar tendency, while s_fold is much smaller. On torsion angle space, structures in the folded state are scattered more than those with Cartesian coordinates (supplementary material, Fig. S2b). In contrast to these feature variables, both s_fold and s_unfold are positive when we use the contact likelihood as defined in Eq. (2), suggesting that the contact likelihood can describe conformational heterogeneity in the folded and unfolded states almost equally. Hereafter, we use the contact likelihood as feature variables in the dimensionality reduction of conformational landscapes.

Next, we examined two-dimensional conformational landscapes of protein G obtained from PCA, tICA, and UMAP. Regarding UMAP, NN is a user-dependent hyperparameter and the users should select a value of NN that produce a reasonable landscape. In this study, we tested five different NN values around its default value of 15 and, consequently, selected the default value, considering continuity of a dimensionality reduced landscape. In Fig. 5(a), the UMAP landscape consists of two major basins (Major1 and Major2), which are connected with two local basins, Bridge1 and Bridge2. In addition, there are two bulges, which are named as Sub1 and Sub2. In Fig. 5(b), six snapshots are taken from Major1, Major2, Bridge1, Bridge2, Sub1, and Sub2, representing key structural features of individual states. As all the structures in Major1 do not form β₁₂ and β₃₄ at the same time (Fig. 6), Major1 is regarded as the unfolded state. In contrast, structures in Major2 keep all the secondary structures observed in the folded state. This is consistent with the scatter plot colored with the RMSD values against the native structure (supplementary material, Fig. S3a). The structures in bridge regions represent unfolded or partial unfolded ones. In Bridge1, neither β₁₂ nor β₃₄ is formed, while the N- and C-terminal regions are close to each other. In Bridge2, only β₃₄ is formed, keeping the interaction with the N-terminal region. These structures appear to be transient structures in the folding/unfolding transitions. The RMSD values from the native structure are less than 5 Å in Sub1 and Sub2, while they are slightly different from the native one. β₃₄ in structures belongs to Sub1 is dissociated from α-helix and there is an insertion of the C-terminal region between β₁₂ and α-helix in Sub2.

The landscapes obtained from PCA and tICA are shown in Figs. 5(c) and 5(d), respectively. Total contributions of two dominant components in PCA and tICA are 72.6% and 84.5%, respectively. We obtained an optimal lag time of 2.5 ns in tICA. Both landscapes involve Major1 and Major2 as two most stable basins (supplementary material, Figs. S3b and S3c). “Bridge” structures are distributed between the two major basins as we see in the UMAP landscape. Major1 and Major2 are mainly distinguished on the first independent component (IC) of tICA, indicating that the folding/unfolding transition is the slowest mode. Furthermore, the second IC resolves another metastable state, Major3 [Fig. 5(d)]. According to the supplementary material, Fig. S3c, Major3 includes almost folded structures considering the RMSD values in that state, but IC2 obviously indicates that there are transitions between Major3 and the other states including folded Major2 with a certain relaxation time. Therefore, Major3 should be distinguished from Major2, even though they are not distinguished in the PCA landscape or in terms of RMSD. In both PCA and tICA, we can interpret the obtained axes (components) as modes, which can be interpreted as structural features. For instance, PC1 and IC1 obviously correspond to the formation of native contacts. On the other hand, we cannot simply interpret Dimension1 in the UMAP landscape as native contact formation, because there are little changes of contact formation in Major1, even though the distribution of Major1 occupies almost half the size of the whole landscape along with Dimension1. This is one of the biggest differences between linear and nonlinear dimensionality reduction methods. Folding transition mechanisms are understood in UMAP landscape using representative structures (see Fig. 5) and projected structural features (see Fig. 6).

Folding/unfolding trajectories in MD simulations are reasonably mapped along the landscapes in all three cases (supplementary material, Figs. S4 and S5). Then, how the dimensionality reduced landscapes from UMAP, PCA, and tICA provide molecular mechanisms for folding dynamics of protein G? According to the previous protein G simulations based on the same CG model,³⁸ the β₃₄ formation occurred in advance of β₁₂. Figure 6 clearly shows that the UMAP landscape describes that the transition from Bridge1 to Bridge2 corresponds to the β₃₄ formation, and subsequently β₁₂ is formed toward the native structure. The β₃₄ formation seems to happen earlier than the β₁₂ formation in the PCA landscape. In contrast. The tICA landscape does not distinct the order of β₁₂ and β₃₄ formations (Fig. 6 right), suggesting that the first IC might not describe the folding process as good as UMAP and PCA. This would be reasonable if we consider rapid folding events of protein G in CG MD simulations (typically happened within 0.05–0.25 ns, supplementary material, Figs. S4 and S5). Each folding event occurs too rapid compared to the current lag time (2.5 ns) used in the current tICA. A proper balance between the lag time and transition times is required for obtaining meaningful conformational landscape reduced with tICA.

In the UMAP landscape, Sub1 and Sub2 are separated from the folded state, even though the RMSD values of the structures are less than 5 Å. In Fig. 7, we project snapshots belong to Sub1 and Sub2 in the PCA or tICA landscape and found that Major 3 in the tICA landscape is almost equivalent to Sub1 in the UMAP landscape, suggesting that UMAP is able to distinct metastable states that are kinetically trapped. Instead, Sub1 is merged into the folded state in the PCA landscape. Sub2 is not resolved as an isolated cluster just using the first and second components in PCA and tICA. With the third component of PCA and tICA, Sub2 is successfully isolated as shown in the supplementary material, Fig. S6. These results suggest that the UMAP landscape successfully isolates the folded, unfolded, and essential intermediate states more clearly on and off the folding pathways.

Next, we carried out MSM after reducing the dimensionality of conformational space of protein G with PCA, tICA, and UMAP. First, we compute GMRQ as a function of the number of clusters to determine the common number of clusters in MSM with different dimensionality reduction methods. A value of lag time in GMRQ analysis was fixed to 2.5 ns according to the result of tICA. In the supplementary material, Fig. S7a, we can see that the optimum number of states is 80 in case of UMAP. In contrast, there are not clear tendency in the GMRQ plots of PCA and tICA clustering. We, therefore, decided to use the same number of states 80 for the three MSMs. Note that this number of clusters is much larger than that in conventional MSM analyses. This large number of clusters is just a sufficient condition because this value was determined according to only the GMRQ convergence in the microstate MSM. An essential number of clusters for this system should be investigated through lumping microstates into a small number of macrostates.⁴⁶ We also analyze characteristics of clustering with different methods, by comparing the probability distributions of intra- and inter-cluster deviations of contact likelihood defined as

(supplementary material, Fig. S7b). Here, CL_im is contact likelihood for the mth contact in the ith snapshot. M is the number of contacts in the native structure, and N is the number of possible pairs for i and j. Snapshots i and j are in the same cluster for intra-cluster deviation, and not for inter-cluster deviation. PCA and tICA show similar probability distributions of intra-cluster deviations with a single peak around 0.2 although the distribution of PCA slightly shifts smaller than that of tICA. Unlike these two methods, probability distribution in UMAP clustering has two peaks with a split around 0.2. The supplementary material, Fig. S7c, shows PCA/tICA/UMAP landscapes with the cluster centers colored by average intra-cluster deviation values. The UMAP landscape has so many clusters, each of which intra-cluster deviation is around 0.1, resulting in the bimodal distribution. Regarding the inter-cluster deviation, there is no difference among three methods, and note that these tendencies of intra- and inter-cluster deviations also appear in the case of different number of states (supplementary material, Fig. S8).

Based on the state discretization in Fig. S7c, we construct MSMs and validate them using implied time scale (ITS) defined as

Here, $λiτ$ is the ith eigenvalue of $Pτ$⁠. If the constructed MSM is Markovian, ITS would be a constant value with the lag time longer than a certain value τ_T (τ_T is called Markovian lag time). In other words, by using ITS analysis, we can identify significant transition modes of which ITS value converged with a lag time longer than τ_T. Figure 8(a) shows the dependency of ITS on a lag time, and all three MSMs are converged with τ_T = 2.5 ns. Using this τ_T, we construct MSMs with UMAP, PCA, and tICA. In each MSM, the stationary distribution (first eigenvector) is similar to the original conformational landscape (supplementary material, Fig. S9). Other eigenvectors representing slow dynamics of the system can resolve three transition modes in the UMAP MSM [left in Fig. 8(a)], while there are only one or two transition modes in the PCA/tICA MSM [middle and right of Fig. 8(a)].

We also plot the landscapes colored by eigenvector values [Fig. 8(b)]. The second eigenvector modes in these MSMs are in common: they describe folding/unfolding transition of protein G. The third eigenvector mode in the UMAP and tICA MSMs represents a transition to the metastable states, Sub1 or Major3. The fourth eigenvector mode, which is meaningful only in the UMAP MSM, describes the transition to Sub2. These transition modes are consistent with the results discussed in Sec. III C. UMAP is, thus, applicable to the preprocessing of MSM and can resolve slow dynamics of a target biomolecule with fewer dimensions than tICA.

In this study, we proposed a new approach to reduce the dimensionality of conformational landscape with UMAP and the contact likelihood. The method was applied to CG MD simulation trajectories of folding dynamics of protein G. By the comparison between UMAP and PCA/tICA, we have shown that (i) UMAP can distinct essential intermediate structures as well as the folded/unfolded states, (ii) folding mechanisms of protein G are well explained in the dimensionality reduced space, and (iii) MSM based on UMAP is better than MSM with tICA or PCA. Interestingly, UMAP is able to resolve the metastable states, which might be distinguished only by tICA with dynamical information. This suggests potential applications of this method for investigating kinetics of protein motions. Regardless of dimensionality reduction methods, the selection of feature variables is found to be essential. Here, we used native contact likelihood as feature variables of folding/unfolding transition of protein G. Considering that the native contact is a key quantity to describe the nature of protein folding funnels,^11,54,55 the feature variables would be applicable to folding/unfolding dynamics of other globular proteins. Another advantage of UMAP for the dimensionality reduction of conformational landscape is that it does not require dynamical information. This opens a new possibility to apply it to MD simulation trajectories with biased samplings or enhanced conformational sampling algorithms.^56–58 Since these methods allow us to explore a wider conformational space of a protein or other biomolecules, dimensionality reduced landscape via UMAP and appropriate feature variables would provide essential information on thermodynamics and kinetics of the molecules more easily than the conventional methods.

See the supplementary material for the methodological explanations on contact likelihood, Figs. S1–S9 for the UMAP, PCA, and tICA landscapes.

We used the computational resources provided by the HPCI System Research Project (Project Nos. hp210107 and hp210177) and those in RIKEN Advanced Center for Computing and Communication (Hokusai BigWaterfall). This work was supported, in part, by MEXT/JSPS KAKENHI [Grant Nos. 19H05645 and 21H05249 (to Y.S.)], RIKEN pioneering projects in “Biology of Intracellular Environments,” “Glycolipidlogue” (to Y.S.), and MEXT “Program for Promoting Research on the Supercomputer Fugaku (Biomolecular dynamics in a living cell/MD-driven Precision Medicine).” M.O. received financial support from RIKEN special post-doctoral fellowship.

The authors have no conflicts to disclose.

Mao Oide: Conceptualization (lead); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Investigation (lead); Methodology (lead); Software (lead); Visualization (lead); Writing – original draft (equal); Writing – review & editing (supporting). Yuji Sugita: Conceptualization (supporting); Funding acquisition (lead); Supervision (lead); Validation (lead); Writing – original draft (equal); Writing – review & editing (lead).

The data that support the findings of this study are available from the corresponding author upon reasonable request.