Keywords: Ultrametricity; Energy landscapes; Random graphs; Kramers kinetics; Mahalanobis metric

1 Introduction↩︎

The hierarchical organization of energy landscapes is one of the central concepts of modern physics of disordered systems, spin glass theory, and macromolecular biophysics. Mathematically, such an organization is described by ultrametric spaces – metric spaces in which, for any three points \(x\), \(y\), \(z\), the metric \(d\) satisfies the condition \[d(x,y)\le\max\{d(x,z),d(z,y)\},\] which is called the strong triangle inequality. The fundamental properties of ultrametrics, including their connection to rooted trees and hierarchical clustering, are detailed in the classic review [1] and monographs (see, e.g., [2], [3]). These ideas have found applications in a wide variety of fields – from phylogenetic analysis [4] and taxonomy theory to statistical physics and \(p\)-adic analysis [5]–[7]. A key event that stimulated interest in ultrametricity in physics was Parisi’s discovery of the ultrametric organization of state space in the Sherrington–Kirkpatrick spin glass model [8], [9]. Further research [10], [11] and rigorous mathematical proofs [12]–[14] confirmed that ultrametricity is not an artificial construct but reflects the real geometry of landscapes with many competing minima separated by barriers of increasing height. Almost simultaneously with these works, Frauenfelder and colleagues, while investigating the kinetics of CO binding to myoglobin, formulated the hypothesis of a hierarchical organization of the conformational space of proteins [15]–[17]. According to this hypothesis, protein dynamics represent a random walk on an ultrametric tree of metastable substates grouped into clusters separated by barriers. The development of these concepts led to the formation of a whole research area at the intersection of \(p\)-adic analysis and molecular biophysics [18]–[21], within which specific models of \(p\)-adic diffusion on protein energy landscapes were constructed.

The motivation for this work comes from recent results obtained by the author in analyzing the ultrametricity of a kinetic metric on the space of RNA secondary structures [22], as well as subsequent unpublished numerical studies of protein molecules. These studies demonstrated that the kinetic metric on the space of quasi-equilibrium states of biopolymers possesses an unexpectedly high degree of nontrivial ultrametricity. This observation raised a fundamental question – is the discovered hierarchy a unique consequence of biological evolution optimizing the energy landscape of a protein molecule, or does it represent a universal property of any sparse, locally tree-like networks with a rugged energy relief?

To answer this question, in this work we apply the mathematical apparatus developed in [22] to a model of a random energy landscape defined on Erdős–Rényi graphs. We generate an ensemble of sparse random graphs, assign a random energy from a uniform distribution over a given interval to each vertex, and apply the algorithms for basin finding, spectral decomposition, and kinetic metric computation used in previous works. The principal difference between such a system and biopolymer systems is the complete absence of physico-chemical constraints on the graph topology and energy distribution, which allows us to investigate the ultrametricity of the system in its purest, abstract form.

The central result of our work is a mathematical proof of the fact that as the width of the energy distribution tends to infinity, the kinetic Mahalanobis metric, constructed from the spectral decomposition of the Kramers transition rate matrix, converges pointwise to the classical single-linkage ultrametric. This limit theorem not only explains the monotonic increase in the degree of nontrivial ultrametricity with increasing free energy distribution interval observed in computational experiments but also establishes a direct mathematical link between two fundamental constructions – the kinetic metric, which accounts for all possible transition paths, and the minimax metric, defined by a single optimal path.

The paper is organized as follows. Section 2 describes the model of a random Erdős–Rényi graph and phase transitions in such networks. Section 3 is devoted to the construction of the Kramers transition rate matrix and the proof of its spectral properties. Section 4 describes the gradient descent algorithm for identifying macrostates, the spectral decomposition description, and the construction of the kinetic Mahalanobis metric. Section 5 contains the definition of the formal ultrametricity criterion. Section 6 presents the computational experiment results and their analysis. Section 7 contains the formulation and proof of the limit theorem on the convergence of the kinetic metric to the single-linkage metric. Section 8 discusses the physical meaning of the obtained results, provides a comparison with \(p\)-adic models, and outlines directions for further research.

2 Model of a Random Energy Landscape on Sparse Graphs↩︎

As a basic topological structure, we use the classical random Erdős–Rényi graph model \(G(V,E)\) [23], [24]. The graph is constructed on a set of vertices \(\mathcal{M}=\{1,2,\dots,V\}\), which correspond to the Markov microstates of the system. The edges of the graph \(\mathcal{E}\) are formed by randomly selecting \(E\) unordered pairs of vertices with equal probability without replacement, which guarantees the absence of multiple edges.

The key parameter determining the global topology of an Erdős–Rényi graph is the average vertex degree: \[\langle k\rangle=\frac{2E}{V}.\]

The theory of random graphs establishes the existence of a phase transition at \(\langle k\rangle=1\). Namely, the following statement holds [24]. Consider a random graph with \(V\) vertices and \(E\) edges. Let \(V\to\infty\), while the average vertex degree \(\langle k\rangle=2E/V\) remains constant. Then, for \(\langle k\rangle<1\), all connected components of the graph are almost surely trees or contain exactly one cycle, and the size of the largest connected component is of order \(O(\log V)\). At \(\langle k\rangle=1\), the size of the largest connected component is of order \(O(V^{2/3})\). For \(\langle k\rangle>1\), a unique giant connected component almost surely emerges, containing a positive fraction \(S>0\) of all vertices, where \(S\) is determined as the unique positive solution of the equation \(S=1-\exp(-\langle k\rangle S).\) All other connected components have size \(O(\log V)\).

The most important property of sparse Erdős–Rényi graphs is their local tree-likeness. A graph is called locally tree-like if, for any vertex \(v\), its \(r\)-neighborhood \(N_{r}(v)=\{u:d_{\mathrm{graph}}(u,v)\le r\}\) with high probability contains no cycles for \(r=O(\log\log V)\). For an Erdős–Rényi graph with \(\langle k\rangle=O(1)\), the probability of having a cycle of length \(l\) in the neighborhood of an arbitrary vertex is estimated as \(O(\langle k\rangle^{l}/V)\), which tends to zero as \(V\to\infty\). This means that locally the graph is isomorphic to a Cayley tree (regular Bethe tree) with branching degree \(\langle k\rangle\). This property is of crucial importance for kinetics: the absence of local cycles means that between any two distant vertices, there exists, essentially, a unique trunk path. Alternative paths require going to higher hierarchical levels of the tree, which forms a natural barrier hierarchy.

Consider some realization of an Erdős–Rényi graph \(G=(\mathcal{M},\mathcal{E})\). Assign to each vertex \(i\in\mathcal{M}\) of the graph a random free energy value \(F_{i}\), where the vertex energies are chosen independently from a uniform distribution over a given interval \((F_{\min},F_{\max})\). Next, we consider a Markov random walk process on the graph \(G=(\mathcal{M},\mathcal{E})\). To describe the evolution of the system, we use the master equation: \[\frac{df_{i}(t)}{dt}=\sum_{j:\{i,j\}\in\mathcal{E}}\left(k_{j\to i}f_{j}(t)-k_{i\to j}f_{i}(t)\right),\] where \(f_{i}(t)\) is the probability of being in microstate (vertex) \(i\), \(k_{i\to j}\) are the rate constants for transitions between microstates, and the summation is carried out only over neighboring vertices, i.e., over vertices connected by edges.

3 The Kramers Transition Rate Matrix and Its Spectral Properties↩︎

According to the theory of chemical reaction rates [25], the transition rate constant is determined by the height of the energy barrier. In the absence of information about the exact geometry of saddle points, a conservative lower estimate of the barrier through the free energies of the initial and final states is used: for two neighboring microstates \(i\) and \(j\), the rate constant of the transition \(i\to j\) has the form \[k_{i\to j}=\nu_{0}\exp\left(-\frac{\max\{F_{i},F_{j}\}-F_{i}}{RT}\right),\label{eq95kramers95rate95asymmetric}\tag{1}\] where \(\nu_{0}>0\) is the frequency factor (pre-exponential factor) with dimensions \([\text{time}]^{-1}\), \(T\) is the temperature, \(R\) is the universal gas constant, \(F_{i}\) is the free energy expressed in molar units, and the quantity \(\max\{F_{i},F_{j}\}\) is used as a lower estimate of the free energy of the transition state. This estimate (1 ) is standard for models where the exact geometry of the energy surface is unavailable [26]. The constants \(k_{i\to j}\) satisfy detailed balance: \[k_{i\to j}\cdot w_{i}=k_{j\to i}\cdot w_{j},\] where \(w_{i}=\exp(-F_{i}/RT)\) is the Boltzmann weight of microstate \(i\).

Define the matrix \(K\) of size \(V\times V\), where \(V=|\mathcal{M}|\), as follows: \[K_{ij}=\begin{cases} \nu_{0}\exp\left(-\dfrac{\max\{F_{i},F_{j}\}-F_{j}}{RT}\right), & \text{for }\{i,j\}\in\mathcal{E},\\[10pt] 0, & \text{for }\{i,j\}\notin\mathcal{E},\;i\neq j, \end{cases}\] \[K_{ii}=-\sum_{j\neq i}K_{ji}.\] Note that in the convention adopted here, the element \(K_{ij}\) (\(i\neq j\)) represents the transition rate from state \(j\) to state \(i\), as a result of which the sum of the elements of each column of the matrix \(K\) is zero (\(\sum_{i}K_{ij}=0\)), and the generator acts on the column vector of probabilities \(\mathbf{f}\) from the right according to the equation \(\dot{\mathbf{f}}=K\mathbf{f}\). The matrix \(K\) is the generator of a continuous-time Markov process. It is non-symmetric, which complicates the numerical analysis of its spectrum, defined by the standard eigenvalue equation \(K\phi=\lambda\phi\). A standard technique is to transition to a symmetric matrix \(S\), which has the same eigenvalues as the matrix \(K\). Define \(S\) through the similarity transformation \(S=W^{-1/2}KW^{1/2}\), where \(W=\mathrm{diag}(w_{1},\dots,w_{V})\). For \(i\neq j\), the elements of the matrix \(S\) are: \[S_{ij}=\frac{1}{\sqrt{w_{i}}}K_{ij}\sqrt{w_{j}}=\nu_{0}\exp\left(-\frac{|F_{i}-F_{j}|}{2RT}\right)\;\text{for }\{i,j\}\in\mathcal{E},\label{eq95S95offdiag}\tag{2}\] and \(S_{ij}=0\) for non-adjacent vertices. The diagonal elements under this transformation are preserved: \[S_{ii}=K_{ii}=-\sum_{j\neq i}K_{ji}=-\sum_{j\neq i}\nu_{0}\exp\left(-\frac{\max\{F_{i},F_{j}\}-F_{i}}{RT}\right).\label{eq95S95diag}\tag{3}\] The matrix \(S\) constructed in this way is symmetric and has a spectrum that exactly coincides with the spectrum of \(K\). However, \(S\) is not a generator of a Markov process in the strict sense, since the sums of elements in its columns are not zero. The stationary distribution of the original process (eigenvalue \(\lambda_{0}=0\)) corresponds to an eigenvector of the matrix \(S\) of the form \(\boldsymbol{\psi}_{0}=\frac{1}{\sqrt{Z}}W^{1/2}\mathbf{1}\), where \(\mathbf{1}\) is a vector of ones, and \(Z=\sum_{i=1}^{V}w_{i}\) is the total partition function of the system.

The connection between the spectra of matrices \(K\) and \(S\) is established directly through the similarity transformation. Substituting the eigenvector change \(\phi=W^{1/2}\psi\) into the original equation \(K\phi=\lambda\phi\) and multiplying on the left by \(W^{-1/2}\) leads to the standard symmetric eigenvalue problem: \[S\psi=\lambda\psi.\label{eq95standard95eigenvalue}\tag{4}\] Consequently, all eigenvalues \(\lambda_{k}\), and hence the relaxation times \(\tau_{k}=-1/\lambda_{k}\), are preserved upon transition from the original process to the symmetrized description. The replacement of matrix \(K\) with matrix \(S\) is solely due to the computational efficiency of diagonalizing symmetric matrices. The eigenvalues and, consequently, the weights \(1/|\lambda_{k}|\) in the metric are preserved. The eigenvectors \(\psi_{k}\) of the symmetrized matrix \(S\) are orthonormal in the standard Euclidean inner product.

It is easy to see that the matrix \(S\) is negative semi-definite. Indeed, consider an arbitrary vector \(x\in\mathbb{R}^{m}\) and expand the quadratic form for the symmetric matrix \(S\), using the definitions of its elements (2 )–(3 ): \[x^{T}Sx=\sum_{i}S_{ii}x^{2}_{i}+\sum_{i\neq j}S_{ij}x_{i}x_{j}\] \[=-\sum_{i}\left(\sum_{j\neq i}K_{ji}\right)x^{2}_{i}+\sum_{\{i,j\}\in\mathcal{E}}2S_{ij}x_{i}x_{j}.\label{xSx}\tag{5}\] Using the connection between the elements \(S\) and \(K\), write \(S_{ij}=\frac{1}{\sqrt{w_{i}}}K_{ij}\sqrt{w_{j}}=\frac{1}{\sqrt{w_{j}}}K_{ji}\sqrt{w_{i}}\), whence \(K_{ji}=S_{ij}\sqrt{w_{j}/w_{i}}\). Substituting this expression into (5 ), after elementary transformations we obtain: \[x^{T}Sx=-\sum_{\{i,j\}\in\mathcal{E}}S_{ij}\left(\sqrt{\frac{w_{j}}{w_{i}}}x^{2}_{i}+\sqrt{\frac{w_{i}}{w_{j}}}x^{2}_{j}-2x_{i}x_{j}\right)\] \[=-\sum_{\{i,j\}\in\mathcal{E}}S_{ij}\left(\sqrt[4]{\frac{w_{j}}{w_{i}}}x_{i}-\sqrt[4]{\frac{w_{i}}{w_{j}}}x_{j}\right)^{2}\le0.\label{eq95negative95semidef}\tag{6}\] Inequality (6 ) holds due to the non-negativity of all off-diagonal elements \(S_{ij}\) (\(i\neq j\)). Equality to zero is achieved if and only if \(x_{i}/\sqrt{w_{i}}=x_{j}/\sqrt{w_{j}}\) for all pairs \(\{i,j\}\in\mathcal{E}\), that is, when the vector \(x\) is proportional to \(\sqrt{w}\) on each connected component of the graph \(G\).

From negative semi-definiteness, it follows that \(S\) has a full set of eigenvalues \(\lambda_{0},\lambda_{1},\dots,\lambda_{V-1}\in\mathbb{R}\), satisfying the condition \(0=\lambda_{0}\ge\lambda_{1}\ge\lambda_{2}\ge\dots\ge\lambda_{V-1}\). The quantities \(\tau_{k}=-1/\lambda_{k}\) (\(k\ge1\)) have the dimension of time and are relaxation times.

4 Basins of Attraction and the Kinetic Metric↩︎

Let a parameter \(\varepsilon_{\mathrm{eq}}>0\) be given, called the energy equivalence threshold. An attraction point is defined as such a microstate \(i\) for which all neighbors (in the sense of graph \(G\)) have a free energy not less than \(F_{i}-\varepsilon_{\mathrm{eq}}\): \[\forall j\in\mathcal{N}(i):F_{j}\ge F_{i}-\varepsilon_{\mathrm{eq}},\] where \(\mathcal{N}(i)\) is the set of neighbors of vertex \(i\) in graph \(G\). Attraction points connected by a path in graph \(G\), all edges of which connect microstates with a free energy difference not exceeding \(\varepsilon_{\mathrm{eq}}\), are merged into one attraction point (this procedure allows merging degenerate minima into gentle plateaus). Each attraction point (possibly after merging) becomes the center of a basin of attraction.

The basin of attraction \(\mathcal{B}_{a}\), corresponding to a given attraction point \(i_{a}\), is the set of all microstates \(i\) for which the gradient descent trajectory ends at \(i_{a}\). The gradient descent procedure is defined as follows. For each microstate \(i\) not yet assigned to any basin, among all its neighbors \(j\) having a free energy strictly less than \(F_{i}-\varepsilon_{\mathrm{eq}}\), the neighbor with the minimum energy is chosen, and microstate \(i\) is assigned to the same basin as this neighbor. If a microstate has no neighbors with energy strictly less than \(F_{i}-\varepsilon_{\mathrm{eq}}\), then it itself is an attraction point and forms a new basin. The procedure is repeated for all microstates in ascending order of free energy. As a result, a family of basins of attraction \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K}\}\) is formed, where the number \(K\) is determined by the graph topology and the value of \(\varepsilon_{\mathrm{eq}}\). Each basin \(\mathcal{B}_{a}\) represents a subset of microstates from which gradient descent trajectories lead to the same attraction point.

An important aspect of constructing basins of attraction is the problem of microstate graph connectivity. Due to the stochastic nature, graph \(G\) may break up into several connected components. To correctly account for this factor, we will use an analysis by connected components, completely analogous to that implemented in [22]. Let graph \(G\) break up into \(C\) connected components \(G_{1},\dots,G_{C}\). For each component \(G_{c}\) containing at least \(K_{\min}\) basins of attraction (where \(K_{\min}\) is a given threshold), a symmetrized transition rate matrix \(S_{c}\) is independently constructed (as a submatrix of the original matrix on the indices of microstates belonging to \(G_{c}\)). To assess the degree of component isolation, a fragmentation index is introduced, defined as \[f_{\mathrm{inter}}=1-\frac{\sum_{c}\binom{K_{c}}{3}}{\binom{K_{\mathrm{total}}}{3}},\] where \(K_{c}\) is the number of basins in component \(c\), \(K_{\mathrm{total}}=\sum_{c}K_{c}\) is the total number of basins in significant components, and \(\binom{K}{3}\) is the number of unordered triples. The index \(f_{\mathrm{inter}}\) takes values from 0 (all basins belong to one component) to 1 (each basin is isolated). A high value of \(f_{\mathrm{inter}}\) indicates that the significant basins of attraction are distributed across several mutually disconnected graph components. In this case, the observed ultrametric structure characterizes the mutual arrangement of basins within each individual component but does not reflect the global organization of the entire energy landscape, since transitions between different connected components are forbidden.

Let us proceed to the construction of the kinetic metric. For each basin of attraction \(\mathcal{B}_{a}\), its conditional equilibrium distribution vector \(\mathbf{p}_{a}\in\mathbb{R}^{V}\) is defined with components \[p_{a}(i)=\begin{cases} w_{i}/W_{a}, & \text{if }i\in\mathcal{B}_{a},\\ 0, & \text{otherwise}, \end{cases}\] where \(W_{a}=\sum_{j\in\mathcal{B}_{a}}w_{j}\) is the partition function of the basin.

Next, a fundamental question arises: how to correctly construct a kinetic metric using the eigenvectors of the symmetric matrix \(S\), if the original physical process is described by the non-symmetric matrix \(K\)? Direct projection of \(\mathbf{p}_{a}\) onto the eigenvectors \(\boldsymbol{\psi}_{k}\) of matrix \(S\) would be erroneous, since \(\boldsymbol{\psi}_{k}\) are not eigenvectors of \(K\). Below we provide a justification of exactly how projections should be computed to obtain a metric corresponding to the original process.

The original matrix \(K\) satisfies the detailed balance condition: \(K_{ij}w_{j}=K_{ji}w_{i}\), where \(w_{i}=\exp(-F_{i}/RT)\). Introduce the diagonal matrix of stationary weights \(W=\mathrm{diag}(w_{1},\dots,w_{V})\). The detailed balance condition is equivalent to the symmetry of the matrix \(KW\): \((KW)^{T}=KW\), where \(W=\mathrm{diag}(w_{1},\dots,w_{V})\) is the diagonal matrix of stationary weights. This means that the matrix \(K\) is similar to the symmetric matrix \(S=W^{-1/2}KW^{1/2}\) and, consequently, is diagonalizable with a real spectrum.

Consider for the matrix \(K\) the right eigenvectors \(\boldsymbol{\phi}_{k}\), satisfying \(K\boldsymbol{\phi}_{k}=\lambda_{k}\boldsymbol{\phi}_{k}\), and the left eigenvectors \(\tilde{\boldsymbol{\phi}}_{k}\), satisfying \(\tilde{\boldsymbol{\phi}}^{T}_{k}K=\lambda_{k}\tilde{\boldsymbol{\phi}}^{T}_{k}\). From the detailed balance condition, it follows that the left and right vectors are related by a simple algebraic relation: \[\tilde{\boldsymbol{\phi}}_{k}=W^{-1}\boldsymbol{\phi}_{k}.\] Indeed, transposing the equality \(K\boldsymbol{\phi}_{k}=\lambda_{k}\boldsymbol{\phi}_{k}\), we obtain \(\boldsymbol{\phi}^{T}_{k}K^{T}=\lambda_{k}\boldsymbol{\phi}^{T}_{k}\). Multiplying by \(W^{-1}\) and using \(K^{T}W^{-1}=W^{-1}K\) (which is equivalent to detailed balance), we arrive at \(\boldsymbol{\phi}^{T}_{k}W^{-1}K=\lambda_{k}\boldsymbol{\phi}^{T}_{k}W^{-1}\), i.e., \(\tilde{\boldsymbol{\phi}}^{T}_{k}K=\lambda_{k}\tilde{\boldsymbol{\phi}}^{T}_{k}\) with \(\tilde{\boldsymbol{\phi}}_{k}=W^{-1}\boldsymbol{\phi}_{k}\). The right and left eigenvectors form a biorthogonal system: \[\langle\tilde{\boldsymbol{\phi}}_{k},\boldsymbol{\phi}_{\ell}\rangle=\delta_{k\ell}.\] An arbitrary initial probability distribution \(\mathbf{f}(0)\) is expanded over the right eigenvectors with coefficients computed via the inner product with the left vectors: \[\mathbf{f}(0)=\sum_{k}c_{k}\boldsymbol{\phi}_{k},\;c_{k}=\langle\tilde{\boldsymbol{\phi}}_{k},\mathbf{f}(0)\rangle.\] This is the standard biorthogonal expansion. The physical projection of the macrostate \(\mathbf{p}_{a}\) onto the \(k\)-th relaxation mode is computed via the inner product with the left eigenvector: \[\pi^{(\mathrm{phys})}_{ak}=\langle\tilde{\boldsymbol{\phi}}_{k},\mathbf{p}_{a}\rangle=\sum_{i}w^{-1}_{i}\phi_{k}(i)\frac{w_{i}}{W_{a}}=\frac{1}{W_{a}}\sum_{i\in\mathcal{B}_{a}}\phi_{k}(i).\label{eq95phys95projection}\tag{7}\] It is these coefficients that are the physical projections and must be used in the formula for the Mahalanobis distance.

Computing the eigenvectors of the non-symmetric matrix \(K\) for large sparse systems is numerically unstable and computationally expensive. Therefore, we use the symmetric matrix \(S=W^{-1/2}KW^{1/2}\), whose eigenvectors \(\boldsymbol{\psi}_{k}\) are related to the right eigenvectors of \(K\) by a similarity transformation: \[\boldsymbol{\psi}_{k}=W^{-1/2}\boldsymbol{\phi}_{k}\;\Longleftrightarrow\;\boldsymbol{\phi}_{k}=W^{1/2}\boldsymbol{\psi}_{k}.\] Substituting this relation into formula (7 ), we obtain \(\pi^{(\mathrm{phys})}_{ak}=\frac{1}{W_{a}}\sum_{i\in\mathcal{B}_{a}}w^{1/2}_{i}\psi_{k}(i).\) Therefore, the physical projection is expressed through the standard Euclidean inner product in the symmetric basis: \[\pi^{(\mathrm{phys})}_{ak}=\langle\mathbf{u}_{a},\boldsymbol{\psi}_{k}\rangle,\label{eq95final95projection}\tag{8}\] where \(\mathbf{u}_{a}\) is a vector with components \(u_{a}(i)=w^{1/2}_{i}/W_{a}\) for \(i\in\mathcal{B}_{a}\) and \(0\) otherwise. Thus, using vectors \(\mathbf{u}_{a}\) together with eigenvectors \(\boldsymbol{\psi}_{k}\) of matrix \(S\) yields exactly the same numerical projection values as using the left/right vector pair of the original matrix \(K\).

The square of the Mahalanobis distance [27] between basins \(\mathcal{B}_{a}\) and \(\mathcal{B}_{b}\) is defined as the weighted sum of squared differences of the physical projections: \[D^{2}(\mathcal{B}_{a},\mathcal{B}_{b})=\sum^{n_{\mathrm{phys}}}_{k=1}\frac{1}{|\lambda_{k}|}\,\left(\pi^{(\mathrm{phys})}_{ak}-\pi^{(\mathrm{phys})}_{bk}\right)^{2}.\label{eq95mahalanobis95macro}\tag{9}\] Since the eigenvalues \(\lambda_{k}\) of matrices \(K\) and \(S\) coincide (they are invariant under the similarity transformation), and the projections \(\pi^{(\mathrm{phys})}_{ak}\) are computed correctly via (8 ), this metric is the exact kinetic metric of the original Markov process. The replacement of the non-symmetric problem with a symmetric one is purely computational in nature and does not change the physical content of the metric.

The weights \(1/|\lambda_{k}|\) in formula (9 ) have the following physical interpretation. The quantity \(|\lambda_{k}|\) characterizes the rate of the \(k\)-th relaxation process – the slower the process (the smaller \(|\lambda_{k}|\)), the greater the contribution made by the difference in projections onto the corresponding mode to the distance between macrostates. This reflects the fact that transitions between macrostates having different projections onto slow modes are hindered by high energy barriers.

The coefficients \(\pi^{(\mathrm{phys})}_{ak}\) are the coordinates of the vector \(\mathbf{u}_{a}\) in the basis of eigenvectors \(\{\boldsymbol{\psi}_{k}\}^{V-1}_{k=1}\) of matrix \(S\). Physically, \(\pi^{(\mathrm{phys})}_{ak}\) determines the amplitude with which the \(k\)-th relaxation mode is represented in the probability distribution concentrated on basin \(\mathcal{B}_{a}\): the larger \(|\pi^{(\mathrm{phys})}_{ak}|\), the greater the contribution of the \(k\)-th mode to the dynamics of the transition from this basin. The quantity \(|\lambda_{k}|^{-1}\) is the relaxation time of the \(k\)-th mode. The distance \(D(\mathcal{B}_{a},\mathcal{B}_{b})\) is computed through the differences of projections \(\pi^{(\mathrm{phys})}_{ak}-\pi^{(\mathrm{phys})}_{bk}\), since it is the difference in mode amplitudes between two basins that determines the kinetic non-equivalence of their relaxation properties: the more the projections onto slow modes differ, the greater the effective transition time between basins.

The sum of squared projections \(\sum^{V-1}_{k=1}(\pi^{(\mathrm{phys})}_{ak})^{2}\) is not equal to one and depends on the basin, since the vectors \(\mathbf{u}_{a}\) are not normalized: \(\|\mathbf{u}_{a}\|^{2}=1/W_{a}\). From the orthonormality of the basis \(\{\boldsymbol{\psi}_{k}\}^{V-1}_{k=0}\) and the fact that the normalized zero eigenvector is \(\boldsymbol{\psi}_{0}=W^{1/2}\mathbf{1}/\sqrt{Z}\), the projection onto the zero mode is \(\langle\mathbf{u}_{a},\boldsymbol{\psi}_{0}\rangle=1/\sqrt{Z}\). According to Parseval’s identity: \[\|\mathbf{u}_{a}\|^{2}=\langle\mathbf{u}_{a},\boldsymbol{\psi}_{0}\rangle^{2}+\sum^{V-1}_{k=1}(\pi^{(\mathrm{phys})}_{ak})^{2}\;\implies\;\frac{1}{W_{a}}=\frac{1}{Z}+\sum^{V-1}_{k=1}(\pi^{(\mathrm{phys})}_{ak})^{2},\] from which it follows that \[\sum^{V-1}_{k=1}(\pi^{(\mathrm{phys})}_{ak})^{2}=\frac{1}{W_{a}}-\frac{1}{Z}.\] Different values of this sum for different basins reflect their unequal kinetic prominence. The smaller the statistical weight \(W_{a}\) (i.e., the narrower and deeper the basin), the larger \(\sum_{k}(\pi^{(\mathrm{phys})}_{ak})^{2}\), which means the dominance of slow modes in the dynamics of the transition from such a basin. Conversely, basins with large \(W_{a}\) (wide and gentle) have a small sum of squared projections, as their dynamics are concentrated in fast modes. The complete invariant, independent of the basin, is \[W_{a}\left(\langle\mathbf{u}_{a},\boldsymbol{\psi}_{0}\rangle^{2}+\sum^{V-1}_{k=1}(\pi^{(\mathrm{phys})}_{ak})^{2}\right)=W_{a}\left(\frac{1}{Z}+\frac{1}{W_{a}}-\frac{1}{Z}\right)=1.\label{eq95inv}\tag{10}\] Relation (10 ) means that upon multiplication by the statistical weight, the sum of all squared projections (including the zero mode) is normalized to one for any basin.

Introduce on the set of found basins of attraction \(\left\{ \mathcal{B}_{a}\right\}\) an equivalence relation \(\sim_{m}\), defined by the equality of projections onto the subspace \(\mathcal{V}\): \[\mathcal{B}_{a}\sim_{m}\mathcal{B}_{b}\;\Longleftrightarrow\;\pi^{(\mathrm{phys})}_{a}=\pi^{(\mathrm{phys})}_{b}\;\Longleftrightarrow\;D(\mathcal{B}_{a},\mathcal{B}_{b})=0.\] Denote by \(\left\{ \mathcal{B}_{a}\right\} /\sim_{m}\) the factor set of equivalence classes under this relation. Then it can be proved (see [22] for details) that the function \(D\) induces a metric on the factor set \(\left\{ \mathcal{B}_{a}\right\} /\!\sim_{m}\). We emphasize that the metric \(D\) is not automatically an ultrametric, since the Euclidean metric in the weighted projection space does not necessarily satisfy the strong triangle inequality. Therefore, checking the ultrametricity of \(D\) is a substantive task, the result of which is determined by the real structure of the energy landscape.

5 Ultrametricity Criterion↩︎

Let \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K_{c}}\}\) be the set of basins in connected component \(c\). To check the ultrametricity of a given component, we will analyze all unordered triples of distinct elements from \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K_{c}}\}\).

Let \((\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime},\mathcal{B}^{\prime\prime\prime})\) be a triple of distinct elements from \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K_{c}}\}\). Sort the three pairwise distances in ascending order: \(d_{\min}\le d_{\text{med}}\le d_{\max}\). Fix two real parameters \(\varepsilon\) and \(\delta\), satisfying the condition \(0<\varepsilon<\delta\).

Triples of elements are classified according to the following rules. A triple \((\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime},\mathcal{B}^{\prime\prime\prime})\) is called trivially ultrametric if the condition \[\frac{d_{\max}-d_{\min}}{d_{\min}}\le\varepsilon.\label{eq95trivial}\tag{11}\] is satisfied. In the case \(d_{\min}=0\), the triple is considered trivially ultrametric only if \(d_{\max}=0\). Trivial ultrametricity means that all three distances are approximately equal. Such a situation can occur, for example, when points are randomly chosen in a high-dimensional space [28], [29] and does not necessarily indicate the presence of a hierarchical structure. A triple \((\mathcal{B}^{\prime},\mathcal{B}^{\prime\prime},\mathcal{B}^{\prime\prime\prime})\) is called nontrivially ultrametric if two conditions are satisfied simultaneously: \[\frac{d_{\max}-d_{\text{med}}}{d_{\text{med}}}\le\varepsilon,\;\frac{d_{\text{med}}-d_{\min}}{d_{\text{med}}}>\delta.\label{eq95nontrivial}\tag{12}\] Nontrivial ultrametricity corresponds to the case where the two largest distances are approximately equal, and the third is significantly smaller. It is precisely this structure of triples that is typical for ultrametric spaces isomorphic to the set of leaves of a multilevel rooted tree. A triple that satisfies neither condition (11 ) nor condition (12 ) is called non-ultrametric.

Let \(\mathcal{T}\) be the set of all unordered triples of distinct elements from \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K_{c}}\}\). Its cardinality is \(\binom{K_{c}}{3}\). Let \(\mathcal{T}_{\text{nt}}\subset\mathcal{T}\) be the subset of nontrivially ultrametric triples. The degree of nontrivial ultrametricity of the component space \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K_{c}}\}\) is defined as \[u_{\text{nt}}=\frac{|\mathcal{T}_{\text{nt}}|}{|\mathcal{T}|}\times100\%.\label{eq95u95nt}\tag{13}\] Similarly, the degree of trivial ultrametricity \(u_{\text{tr}}\) and the degree of non-ultrametricity \(u_{\text{non}}\) are defined, with \(u_{\text{nt}}+u_{\text{tr}}+u_{\text{non}}=100\%\). The degrees of trivial ultrametricity, nontrivial ultrametricity, and the degree of non-ultrametricity for the entire set of basins are defined as weighted sums of the corresponding ultrametricity degrees of its connected components.

Note that in continuous metrics, exact equalities required by conditions (11 ) and (12 ) for \(\delta=0\), \(\varepsilon=0\) have zero probability. Therefore, in computational experiments, approximate inequalities with specific values of \(\varepsilon\) and \(\delta\) are used, chosen based on the required classification strictness.

6 Computational Results↩︎

To perform the calculations, we developed a specialized program in the Python interpreter, the source code of which is publicly available [30]. Using this program, a series of computational experiments was conducted to test the hypothesis of the ultrametric organization of random energy landscapes with the following fixed parameters: number of vertices \(V=5000\), number of edges \(E=5000\) (corresponding to an average degree \(\langle k\rangle=2.0\) and a theoretical giant component fraction \(S\approx0.7968\)); Kramers temperature \(T=300\) K (\(RT\approx2.49\) kJ/mol); energy equivalence threshold \(\varepsilon_{\mathrm{eq}}=3.0\) kJ/mol; minimum basin size \(s_{\min}=10\) vertices; minimum number of basins in a connected component \(K_{\min}=5\); ultrametric triple classification parameters \(\varepsilon=0.05\), \(\delta=0.10\); number of requested eigenmodes \(n_{\mathrm{modes}}=50\); spectral gap threshold \(\theta=10^{6}\); number of independent realizations for each parameter set \(N=10\). The variable parameter was the width of the uniform free energy distribution interval \(\Delta F=F_{\max}-F_{\min}\) (with fixed \(F_{\min}=0\)). Seven values of the free energy distribution interval were investigated: \(\Delta F=10\), \(20\), \(50\), \(100\), \(200\), \(500\), and \(1000\) kJ/mol.

Analysis of connected components for all realizations showed results fully consistent with the theory of Erdős–Rényi random graphs. The average size of the giant component was \(4000.5\pm34.8\) vertices with an average number of edges inside it of \(4801.3\pm31.9\). The fraction of vertices in the giant component \(S_{\mathrm{obs}}\approx0.800\) practically coincides with the theoretical value \(S=0.7968\). The remaining vertices broke up into many small components, the size of which did not exceed a few tens of vertices and which were automatically discarded by the algorithm as not containing a sufficient number of basins. The fragmentation index \(f_{\mathrm{inter}}\) in all successful realizations was equal to zero, indicating that all significant basins belonged to a single giant connected component.

The main results, averaged over ten independent realizations for each value of \(\Delta F\), are presented in Table 1.

As seen from Table 1, a clear monotonic dependence of the degree of nontrivial ultrametricity on the energy distribution width is observed. For small values of \(\Delta F=10\) and \(20\) kJ/mol, the kinetic Mahalanobis metric reveals moderate ultrametricity: about \(42\%\) and \(61\%\) of basin triples are classified as nontrivially ultrametric. At \(\Delta F=50\) kJ/mol (which corresponds to \(\sim20RT\)), the kinetic Mahalanobis metric reveals a high degree of ultrametricity: about \(81\%\) of triples are classified as nontrivially ultrametric. As \(\Delta F\) increases to \(100\) kJ/mol, the fraction of nontrivially ultrametric triples reaches a plateau at about \(\approx83\%\). At \(\Delta F=200\) kJ/mol, another increase occurs: \(u_{\mathrm{nt}}\) rises to \(86.5\%\), while the fraction of non-ultrametric triples drops to \(13\%\). A further increase in \(\Delta F\) to \(500\) kJ/mol and \(1000\) kJ/mol is accompanied by a continuation of the monotonic growth of \(u_{\mathrm{nt}}\) to \(90.6\%\) and \(95.9\%\), respectively. The fraction of non-ultrametric triples at \(\Delta F=1000\) kJ/mol decreases to \(3.6\%\), indicating that the system has entered the asymptotic regime.

The average number of basins of attraction \(\langle K\rangle\) also demonstrates an increasing trend with growing \(\Delta F\) (from \(24.5\) to \(30.3\)); however, this increase is relatively weak compared to the change in \(u_{\mathrm{nt}}\). This observation indicates that the increase in the degree of ultrametricity is not a simple consequence of the landscape fragmenting into a larger number of basins, but reflects a fundamental change in the structure of the kinetic metric.

To illustrate the nature of nontrivially ultrametric triples, Table 2 shows characteristic examples of distance triples for \(\Delta F=50\), \(\Delta F=200\), and \(\Delta F=1000\) kJ/mol (taken from realizations with the maximum \(u_{\mathrm{nt}}\)).

One notes the qualitative difference in the structure of distances between the regimes with small and large \(\Delta F\). At \(\Delta F=50\) kJ/mol, the absolute distance values are of order \(10^{11}\)–\(10^{18}\), and the ratios \(d_{\max}/d_{\min}\) reach several orders of magnitude. At \(\Delta F=200\) kJ/mol, the distances increase to order \(10^{13}\)–\(10^{21}\), with the two largest distances \(d_{\max}\) and \(d_{\mathrm{med}}\) coinciding with high accuracy, while \(d_{\min}\) is smaller by two to three orders of magnitude. At \(\Delta F=1000\) kJ/mol, distances reach astronomical values of order \(10^{63}\)–\(10^{79}\), while \(d_{\max}\) and \(d_{\mathrm{med}}\) are practically indistinguishable, and \(d_{\min}\) is smaller by tens of orders of magnitude. Such a structure is a direct sign that for large \(\Delta F\), the Mahalanobis metric approximates a strict ultrametric ever more accurately, in which for any nontrivial triple, the two largest distances are exactly equal.

The reader may be surprised that the absolute values of the Mahalanobis distances \(D\) reach extremely large magnitudes (of order \(10^{18}\) and even \(10^{63}\)). It is important to emphasize that this is not a numerical artifact or a computational error, but represents a fundamental physical consequence of the Arrhenius law and Kramers theory. The kinetic distance \(D\) is inversely proportional to the square root of the modulus of the transition rate matrix eigenvalue (\(D\sim1/\sqrt{|\lambda|}\)), which physically characterizes the relaxation time (or mean committor time) between macrostates. According to Kramers theory, the transition rate through a minimax barrier of height \(U^{*}\) is exponentially suppressed: \(|\lambda|\sim\exp(-U^{*}/RT)\). Consequently, the kinetic distance grows exponentially with the barrier height and also accounts for entropic factors (statistical weights of basins): \[D\sim\exp\left(\frac{U^{*}}{2RT}\right).\] In our computational experiments with a distribution width \(\Delta F=50\) kJ/mol, the distances between basins already reach values on the order of \(10^{18}\). When the interval increases to \(\Delta F=1000\) kJ/mol and temperature \(T=300\) K (\(RT\approx2.49\) kJ/mol), the typical heights of separating barriers \(U^{*}\) amount to hundreds of kJ/mol. Thus, a barrier of height \(U^{*}\approx730\) kJ/mol gives an asymptotic estimate \(D\sim\exp(730/4.98)\sim10^{63}\), which exactly matches the observed values. The waiting times for such transitions at room temperature exceed the age of the Universe by many orders of magnitude, meaning absolute kinetic isolation of the basins. Thus, the giant values of \(D\) are a direct numerical confirmation of the Limit Theorem formulated in Section 7: the logarithmic asymptotics of the metric \(\frac{1}{\beta}\ln D\) exactly corresponds to the minimax barrier height \(\frac{1}{2}\rho^{(0)}\). For practical analysis and visualization of such landscapes (e.g., when constructing dendrograms), a standard and mathematically correct procedure is the use of a logarithmic distance scale (\(\ln D\)), which is fully consistent with the traditions of spin glass physics and \(p\)-adic analysis.

7 Limit Theorem on the Convergence of the Kinetic Metric to the Single-Linkage Metric↩︎

Let us formulate and prove the main theoretical result of the work. Fix a connected graph \(G=(\mathcal{M},\mathcal{E})\) with a set of vertices \(\mathcal{M}=\{1,\dots,V\}\). Assign to each vertex \(i\) a non-negative number \(F_{i}\ge0\), which we will call the base free energy. Without loss of generality, we assume all quantities \(F_{i}\) are pairwise distinct. Introduce a scale parameter \(\beta>0\) (playing the role of inverse temperature \(1/RT\)) and consider the family of generators \(K(\beta)\).

The original Kramers Markov generator \(K(\beta)\) is defined by its off-diagonal elements \[K_{ij}(\beta)=\exp\!\bigl(-\beta(\max(F_{i},F_{j})-F_{j})\bigr),\;\{i,j\}\in\mathcal{E},\] and diagonal elements \(K_{ii}(\beta)=-\sum_{j\neq i}K_{ji}(\beta)\) (here the pre-exponential factor \(\nu_{0}\) is omitted, as it does not affect the logarithmic asymptotics). The matrix \(K(\beta)\) satisfies the detailed balance condition \(K_{ij}w_{j}=K_{ji}w_{i}\) with stationary Boltzmann weights \(w_{i}(\beta)=\exp(-\beta F_{i})\).

For the correct spectral decomposition of the non-symmetric matrix \(K\), we use the biorthogonal system of left and right eigenvectors. The right eigenvectors \(\boldsymbol{\phi}_{k}\) satisfy the equation \(K\boldsymbol{\phi}_{k}=\lambda_{k}\boldsymbol{\phi}_{k}\), and the left eigenvectors \(\tilde{\boldsymbol{\phi}}_{k}\) satisfy the equation \(\tilde{\boldsymbol{\phi}}^{T}_{k}K=\lambda_{k}\tilde{\boldsymbol{\phi}}^{T}_{k}\). From the detailed balance condition, it follows that the left and right vectors are related by the relation \(\tilde{\boldsymbol{\phi}}_{k}=W^{-1}\boldsymbol{\phi}_{k}\), where \(W=\mathrm{diag}(w_{1},\dots,w_{V})\). They form a biorthogonal basis: \(\langle\tilde{\boldsymbol{\phi}}_{k},\boldsymbol{\phi}_{m}\rangle=\delta_{km}\).

The partition of the set of vertices into basins of attraction does not depend on \(\beta\). Fix this partition \(\{\mathcal{B}_{1},\dots,\mathcal{B}_{K}\}\). For each basin \(\mathcal{B}_{a}\), define the macrostate through the conditional equilibrium probability distribution \(\mathbf{p}_{a}\) with components \[(\mathbf{p}_{a})_{i}=p_{a,i}=\begin{cases} w_{i}/W_{a}, & i\in\mathcal{B}_{a},\\ 0, & i\notin\mathcal{B}_{a}, \end{cases}\label{eq95p95a95i}\tag{14}\] where \(W_{a}=\sum_{j\in\mathcal{B}_{a}}w_{j}\) is the partition function of the basin. The physical projection of the macrostate \(\mathbf{p}_{a}\) onto the \(k\)-th relaxation mode is computed via the inner product with the left eigenvector: \[c_{ak}=\langle\tilde{\boldsymbol{\phi}}_{k},\mathbf{p}_{a}\rangle=\sum_{i\in\mathcal{B}_{a}}w^{-1}_{i}\phi_{k,i}\frac{w_{i}}{W_{a}}=\frac{1}{W_{a}}\sum_{i\in\mathcal{B}_{a}}\phi_{k,i}.\label{eq95proj95def95final}\tag{15}\] Then the kinetic Mahalanobis metric on the set of basins is defined as \[D^{2}_{\beta}(\mathcal{B}_{a},\mathcal{B}_{b})=\sum^{V-1}_{k=1}\frac{1}{|\lambda_{k}(\beta)|}\,\bigl(c_{ak}-c_{bk}\bigr)^{2}.\label{eq95D95def95final}\tag{16}\] The summation starts from \(k=1\), since \(\lambda_{0}=0\) corresponds to the stationary distribution and is excluded from the metric definition.

Let us define the single-linkage metric. Without loss of generality, we will assume that \(\min_{i\in\mathcal{M}}F_{i}=0\). Indeed, the generator matrix \(K\) is invariant under a simultaneous shift of all free energies by an arbitrary constant \(C\), since the transition rates depend only on energy differences. The square of the kinetic metric \(D^{2}_{\beta}\) under such a shift is multiplied by a global factor \(e^{\beta C}\) (and the metric \(D_{\beta}\) itself is multiplied by \(e^{\beta C/2}\)), which leads to an additive shift of the logarithmic asymptotics by \(C/2\). Since the single-linkage distance \(\rho^{(0)}\) also increases by \(C\) under a shift of energies by \(C\), both sides of the theorem formulated below shift by the same amount \(C/2\). Therefore, the choice of normalization \(\min F_{i}=0\) does not limit the generality of the result.

Next, for each edge \(\{i,j\}\in\mathcal{E}\), define its weight \[U_{ij}=\max(F_{i},F_{j}).\] For an arbitrary path \(\gamma=(v_{0},v_{1},\dots,v_{L})\) in the graph \(G\), its height is defined as \[h(\gamma)=\max_{0\le\ell<L}U_{v_{\ell}v_{\ell+1}}.\] The single-linkage distance between vertices \(p\) and \(q\) is defined as \[\rho^{(0)}(p,q)=\min_{\gamma:p\to q}h(\gamma).\] On the set of basins, the distance is given by the formula \[\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})=\min_{p\in\mathcal{B}_{a},\,q\in\mathcal{B}_{b}}\rho^{(0)}(p,q)\] for \(a\neq b\) and \(\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{a})=0\). The function \(\rho^{(0)}\) is an ultrametric: symmetry is obvious, and the strong triangle inequality \[\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})\le\max\!\bigl(\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{c}),\rho^{(0)}(\mathcal{B}_{c},\mathcal{B}_{b})\bigr)\] follows from the fact that the concatenation of optimal paths \(\mathcal{B}_{a}\to\mathcal{B}_{c}\) and \(\mathcal{B}_{c}\to\mathcal{B}_{b}\) (connected inside basin \(\mathcal{B}_{c}\) by a path, all edges of which, by the construction of basins, have weight strictly less than the height of any separating barrier) gives a path from \(\mathcal{B}_{a}\) to \(\mathcal{B}_{b}\), whose height does not exceed the maximum of the heights of the two original paths.

Also define the bottom energy of basin \(\mathcal{B}_{a}\) as the minimum value of the base free energy over the microstates of this basin: \[F_{\min}(\mathcal{B}_{a})=\min_{i\in\mathcal{B}_{a}}F_{i}.\]

Theorem 1. For any two basins \(\mathcal{B}_{a}\) and \(\mathcal{B}_{b}\), the following limit relation holds: \[\lim_{\beta\to\infty}\frac{1}{\beta}\ln D_{\beta}(\mathcal{B}_{a},\mathcal{B}_{b})=\frac{1}{2}\,\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b}).\]

Note that although the physical principle of the dominance of minimax barriers is well known in the theory of large deviations for continuous diffusion processes (Freidlin–Wentzell theory [31]), a rigorous derivation of this asymptotics directly from the spectral definition of the kinetic Mahalanobis metric on a discrete graph is a non-trivial mathematical problem, not reducible to classical results. The proof of the formulated limit theorem relies on the apparatus of linear algebra, graph theory, and variational calculus. We will sequentially transition from the spectral definition of the metric to a variational representation, and then obtain asymptotic lower and upper bounds.

Lemma 1 (Variational representation of the kinetic metric). The square of the kinetic Mahalanobis metric \(D^{2}_{\beta}(\mathcal{B}_{a},\mathcal{B}_{b})\) admits the following representation: \[D^{2}_{\beta}(\mathcal{B}_{a},\mathcal{B}_{b})=\max_{\mathbf{h}\in\mathbb{R}^{V}:\,\langle\mathbf{w},\mathbf{h}\rangle=0}\frac{\bigl(\langle\mathbf{p}_{a},\mathbf{h}\rangle-\langle\mathbf{p}_{b},\mathbf{h}\rangle\bigr)^{2}}{\mathcal{E}(\mathbf{h},\mathbf{h})},\] where \(\mathbf{w}=(w_{1},\dots,w_{V})^{T}\), \(w_{i}=\exp(-\beta F_{i})\), the vectors \(\mathbf{p}_{a}\) and \(\mathbf{p}_{b}\) are given by their components \((p_{a})_{i}=w_{i}/W_{a}\) for \(i\in\mathcal{B}_{a}\) (and zeros outside \(\mathcal{B}_{a}\)), and \(\mathcal{E}(\mathbf{h},\mathbf{h})\) is the Dirichlet quadratic form \[\mathcal{E}(\mathbf{h},\mathbf{h})=\sum_{\{i,j\}\in\mathcal{E}}C_{ij}(h_{i}-h_{j})^{2},\;C_{ij}=\exp\!\bigl(-\beta\max(F_{i},F_{j})\bigr).\]

Proof of Lemma 1. Introduce the symmetric matrix \(L=-S\). From the negative semidefiniteness of \(S\), proved in Section 3, it follows that \(L\) is a positive semidefinite matrix. The kernel of matrix \(L\) is one-dimensional and spanned by the normalized vector \(\boldsymbol{\psi}_{0}=\frac{1}{\sqrt{Z}}W^{1/2}\mathbf{1}\) with components \(\frac{\sqrt{w_{i}}}{\sqrt{Z}}\), where \(Z=\sum^{V}_{i=1}w_{i}\) is the total partition function. According to the standard theory of spectral decomposition of symmetric matrices (see, e.g., [32]), the Moore–Penrose pseudoinverse matrix \(L^{+}\) is defined as \(L^{+}=\sum^{V-1}_{k=1}\frac{1}{\mu_{k}}\boldsymbol{\psi}_{k}\boldsymbol{\psi}^{T}_{k}\), where \(\mu_{k}=-\lambda_{k}>0\) are the nonzero eigenvalues of \(L\), and \(\boldsymbol{\psi}_{k}\) (\(k\ge1\)) are the corresponding orthonormal eigenvectors.

Set \(\mathbf{p}_{a}=(p_{a,1},\dots,p_{a,V})^{T}\) (and analogously \(\mathbf{p}_{b}\)), where \(p_{a,i}\) are defined in (14 ). Then the vector \(\mathbf{u}_{a}=W^{-1/2}\mathbf{p}_{a}\) has components \(u_{a,i}=w^{1/2}_{i}/W_{a}\) for \(i\in\mathcal{B}_{a}\) and zeros outside \(\mathcal{B}_{a}\), which exactly coincides with the definition of the vector \(\mathbf{u}_{a}\) from Section 4. According to (9 ) and (8 ) \(D^{2}_{\beta}=(\mathbf{u}_{a}-\mathbf{u}_{b})^{T}L^{+}(\mathbf{u}_{a}-\mathbf{u}_{b})\) or \[D^{2}_{\beta}=\mathbf{y}^{T}L^{+}\mathbf{y},\;\text{where }\mathbf{y}=W^{-1/2}\mathbf{p}_{a}-W^{-1/2}\mathbf{p}_{b}.\] The vector \(\mathbf{y}\) is orthogonal to the kernel of \(L\), since \(\langle\mathbf{y},\boldsymbol{\psi}_{0}\rangle=\frac{1}{\sqrt{Z}}\left(\sum_{i\in\mathcal{B}_{a}}\frac{w_{i}}{W_{a}}-\sum_{j\in\mathcal{B}_{b}}\frac{w_{j}}{W_{b}}\right)=\frac{1}{\sqrt{Z}}(1-1)=0\). Consequently, \(\mathbf{y}\) belongs to the image of matrix \(L\) (\(\mathbf{y}\in\mathrm{Im}(L)\)).

Let us prove the variational principle: for any \(\mathbf{y}\in\mathrm{Im}(L)\), \(\mathbf{y}\neq\mathbf{0}\), \[\mathbf{y}^{T}L^{+}\mathbf{y}=\max_{\mathbf{x}\in\mathrm{Im}(L),\,\mathbf{x}\neq\mathbf{0}}\frac{(\mathbf{y}^{T}\mathbf{x})^{2}}{\mathbf{x}^{T}L\mathbf{x}}.\label{eq95var95principle}\tag{17}\] (If \(\mathbf{y}=\mathbf{0}\), then both sides are identically zero, and the statement is trivial).

We obtain an upper bound for (17 ). Expand the vectors \(\mathbf{x},\mathbf{y}\in\mathrm{Im}(L)\) in the basis of eigenvectors \(\{\boldsymbol{\psi}_{k}\}^{V-1}_{k=1}\): \(\mathbf{y}=\sum^{V-1}_{k=1}d_{k}\boldsymbol{\psi}_{k}\) and \(\mathbf{x}=\sum^{V-1}_{k=1}c_{k}\boldsymbol{\psi}_{k}\). Then the numerator and denominator of the fraction take the form \[(\mathbf{y}^{T}\mathbf{x})^{2}=\left(\sum^{V-1}_{k=1}c_{k}d_{k}\right)^{2},\;\mathbf{x}^{T}L\mathbf{x}=\sum^{V-1}_{k=1}\mu_{k}c^{2}_{k}.\] The desired quadratic form is \(\mathbf{y}^{T}L^{+}\mathbf{y}=\sum^{V-1}_{k=1}\frac{d^{2}_{k}}{\mu_{k}}\). Apply the classical Cauchy–Bunyakovsky inequality to the vectors \(\mathbf{u}=(c_{1}\sqrt{\mu_{1}},\dots,c_{V-1}\sqrt{\mu_{V-1}})^{T}\) and \(\mathbf{v}=(d_{1}/\sqrt{\mu_{1}},\dots,d_{V-1}/\sqrt{\mu_{V-1}})^{T}\): \[(\mathbf{u}^{T}\mathbf{v})^{2}\le\|\mathbf{u}\|^{2}\|\mathbf{v}\|^{2}\implies\left(\sum^{V-1}_{k=1}c_{k}d_{k}\right)^{2}\le\left(\sum^{V-1}_{k=1}\mu_{k}c^{2}_{k}\right)\left(\sum^{V-1}_{k=1}\frac{d^{2}_{k}}{\mu_{k}}\right).\] Substituting here the expressions for the numerator, denominator, and quadratic form, we obtain \[(\mathbf{y}^{T}\mathbf{x})^{2}\le(\mathbf{x}^{T}L\mathbf{x})(\mathbf{y}^{T}L^{+}\mathbf{y}).\] Since \(\mathbf{x}\in\mathrm{Im}(L)\) and \(\mathbf{x}\neq\mathbf{0}\), we have \(\mathbf{x}^{T}L\mathbf{x}>0\) and arrive at the upper bound \[\frac{(\mathbf{y}^{T}\mathbf{x})^{2}}{\mathbf{x}^{T}L\mathbf{x}}\le\mathbf{y}^{T}L^{+}\mathbf{y}.\label{eq95up}\tag{18}\]

We prove that the found upper bound (18 ) is the exact maximum. To do this, consider the vector \(\mathbf{x}^{*}\in\mathrm{Im}(L)\) on which it is attained. Set \(\mathbf{x}^{*}=L^{+}\mathbf{y}\). Since \(\mathbf{y}\in\mathrm{Im}(L)\) and \(\mathbf{y}\neq\mathbf{0}\), the vector \(\mathbf{x}^{*}\neq\mathbf{0}\) and, by symmetry of \(L\), \(\mathbf{x}^{*}\in\mathrm{Im}(L^{+})=\mathrm{Im}(L)\). Compute the value of the functional \(\max_{\mathbf{x}\in\mathrm{Im}(L),\,\mathbf{x}\neq\mathbf{0}}\frac{(\mathbf{y}^{T}\mathbf{x})^{2}}{\mathbf{x}^{T}L\mathbf{x}}\) on the vector \(\mathbf{x}^{*}\). The numerator and denominator respectively are \[(\mathbf{y}^{T}\mathbf{x}^{*})^{2}=(\mathbf{y}^{T}L^{+}\mathbf{y})^{2},\] \[(\mathbf{x}^{*})^{T}L\mathbf{x}^{*}=(L^{+}\mathbf{y})^{T}L(L^{+}\mathbf{y})=\mathbf{y}^{T}(L^{+})^{T}LL^{+}\mathbf{y}.\] Since \(L\) is symmetric, \(L^{+}\) is also symmetric, i.e., \((L^{+})^{T}=L^{+}\). For the symmetric Moore–Penrose pseudoinverse, the equality \(L^{+}LL^{+}=L^{+}\) holds. Consequently, the denominator is exactly equal to \[\mathbf{y}^{T}L^{+}LL^{+}\mathbf{y}=\mathbf{y}^{T}L^{+}\mathbf{y}.\] Since \(\mathbf{y}\in\mathrm{Im}(L)\) and \(\mathbf{y}\neq\mathbf{0}\), the quantity \(\mathbf{y}^{T}L^{+}\mathbf{y}>0\). Then the value of the fraction on the vector \(\mathbf{x}^{*}\) is \[\frac{(\mathbf{y}^{T}\mathbf{x}^{*})^{2}}{(\mathbf{x}^{*})^{T}L\mathbf{x}^{*}}=\frac{(\mathbf{y}^{T}L^{+}\mathbf{y})^{2}}{\mathbf{y}^{T}L^{+}\mathbf{y}}=\mathbf{y}^{T}L^{+}\mathbf{y}.\] Thus, the upper bound is attained on the vector \(\mathbf{x}^{*}=L^{+}\mathbf{y}\), which proves equality (17 ).

Next, make the change of variables \(x_{i}=\sqrt{w_{i}}h_{i}\). The condition \(\mathbf{x}\in\mathrm{Im}(L)\) is equivalent to \(\langle\mathbf{x},\boldsymbol{\psi}_{0}\rangle=0\), which, taking into account the normalization of \(\boldsymbol{\psi}_{0}\), takes the form \(\frac{1}{\sqrt{Z}}\sum_{i}w_{i}h_{i}=0\), or simply \(\sum_{i}w_{i}h_{i}=\langle\mathbf{w},\mathbf{h}\rangle=0\). The numerator in (17 ) transforms as follows: \[\mathbf{y}^{T}\mathbf{x}=\sum_{i\in\mathcal{B}_{a}}\frac{w_{i}}{W_{a}}h_{i}-\sum_{i\in\mathcal{B}_{b}}\frac{w_{i}}{W_{b}}h_{i}=\langle\mathbf{p}_{a},\mathbf{h}\rangle-\langle\mathbf{p}_{b},\mathbf{h}\rangle.\] For the denominator, we use the definition of the elements of matrix \(S\) from Section 3: \[\mathbf{x}^{T}L\mathbf{x}=\mathbf{x}^{T}(-S)\mathbf{x}=\sum_{\{i,j\}\in\mathcal{E}}S_{ij}\left(\sqrt[4]{\frac{w_{j}}{w_{i}}}x_{i}-\sqrt[4]{\frac{w_{i}}{w_{j}}}x_{j}\right)^{2}.\] Substituting \(x_{i}=\sqrt{w_{i}}h_{i}\), we obtain \(\sqrt[4]{w_{j}/w_{i}}\sqrt{w_{i}}h_{i}=w^{1/4}_{j}w^{1/4}_{i}h_{i}\). The difference in parentheses takes the form \(w^{1/4}_{i}w^{1/4}_{j}(h_{i}-h_{j})\). Squaring and multiplying by \(S_{ij}=\exp(-\beta\left|F_{i}-F_{j}\right|/2)\), we find: \[S_{ij}\sqrt{w_{i}w_{j}}=\exp\!\Bigl(-\beta\frac{\left|F_{i}-F_{j}\right|+F_{i}+F_{j}}{2}\Bigr)=\exp\!\bigl(-\beta\max(F_{i},F_{j})\bigr)=C_{ij}.\] Thus, \(\mathbf{x}^{T}L\mathbf{x}=\mathcal{E}(\mathbf{h},\mathbf{h})\). Substituting these expressions into (17 ) completes the proof of Lemma 1.

Lemma 2 (Topological properties of the single-linkage distance). Let \(\rho=\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})\). Consider the subgraph \(G_{<\rho}\), containing all vertices \(\mathcal{M}\) and only those edges \(\{i,j\}\in\mathcal{E}\) for which \(U_{ij}=\max(F_{i},F_{j})<\rho\). Then the sets \(\mathcal{B}_{a}\) and \(\mathcal{B}_{b}\) belong to different connected components of the subgraph \(G_{<\rho}\).

Proof of Lemma 2. Suppose the contrary: let \(\mathcal{B}_{a}\) and \(\mathcal{B}_{b}\) lie in one connected component of \(G_{<\rho}\). Then there exists a path \(\gamma\) between some vertices \(p\in\mathcal{B}_{a}\) and \(q\in\mathcal{B}_{b}\), all edges of which have weight strictly less than \(\rho\). Consequently, the height of this path \(h(\gamma)<\rho\). By the definition of the single-linkage distance, \(\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})\le h(\gamma)<\rho\), which contradicts the equality \(\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})=\rho\).

Proof of the Theorem. Denote \(\rho=\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b})\). Let us prove that \(\lim_{\beta\to\infty}\frac{1}{\beta}\ln D^{2}_{\beta}=\rho\). The proof will be broken down into 3 steps.

Step 1 (finding the lower bound). Let \(i_{a}=\arg\min_{i\in\mathcal{B}_{a}}F_{i}\) be the vertex with minimum energy in basin \(\mathcal{B}_{a}\). Let \(\mathcal{C}_{a}\) be the connected component of the subgraph \(G_{<\rho}\) containing \(i_{a}\). According to Lemma 2, \(\mathcal{B}_{b}\cap\mathcal{C}_{a}=\emptyset\). Define the trial vector \(\mathbf{h}^{*}\in\mathbb{R}^{V}\) as the indicator of the set \(\mathcal{C}_{a}\): \(h^{*}_{i}=1\) for \(i\in\mathcal{C}_{a}\) and \(h^{*}_{i}=0\) otherwise. For any edge \(\{i,j\}\in\mathcal{E}\), the difference \(h^{*}_{i}-h^{*}_{j}\) is nonzero only if the edge connects a vertex from \(\mathcal{C}_{a}\) with a vertex outside \(\mathcal{C}_{a}\). By definition of a connected component, such an edge is absent in \(G_{<\rho}\), therefore \(U_{ij}\ge\rho\). Compute the quadratic form \(\mathcal{E}(\mathbf{h}^{*},\mathbf{h}^{*})\): \[\mathcal{E}(\mathbf{h}^{*},\mathbf{h}^{*})=\sum_{\{i,j\}\in\mathcal{E}:\,h^{*}_{i}\neq h^{*}_{j}}\exp(-\beta U_{ij})\le\sum_{\{i,j\}\in\mathcal{E}}\exp(-\beta\rho)=|\mathcal{E}|\exp(-\beta\rho).\] Consider the intersection \(\mathcal{B}_{a}\cap\mathcal{C}_{a}\) and the difference \(\mathcal{B}_{a}\setminus\mathcal{C}_{a}\). Let us show that for any vertex \(j\in\mathcal{B}_{a}\setminus\mathcal{C}_{a}\), \(F_{j}\ge\rho\) holds. Indeed, if \(F_{j}<\rho\), then, since the gradient descent procedure from \(j\) to \(i_{a}\) strictly decreases in energy, all vertices and edges on this path would have weight strictly less than \(\rho\). Consequently, the entire path would lie in \(G_{<\rho}\), and \(j\) would belong to the same connected component as \(i_{a}\) (i.e., \(j\in\mathcal{C}_{a}\)), which contradicts the condition \(j\notin\mathcal{C}_{a}\). Thus, for all \(j\in\mathcal{B}_{a}\setminus\mathcal{C}_{a}\), the Boltzmann weight \(w_{j}\le\exp(-\beta\rho)\). The total weight of these vertices is bounded: \(\sum_{j\in\mathcal{B}_{a}\setminus\mathcal{C}_{a}}w_{j}\le V\exp(-\beta\rho)\). Since \(\mathcal{B}_{b}\cap\mathcal{C}_{a}=\emptyset\), the difference of projections (the expression in brackets in the numerator of the functional from Lemma 1) is bounded below as: \[\langle\mathbf{p}_{a},\mathbf{h}^{*}\rangle-\langle\mathbf{p}_{b},\mathbf{h}^{*}\rangle=\frac{1}{W_{a}}\sum_{j\in\mathcal{B}_{a}\cap\mathcal{C}_{a}}w_{j}-0=1-\frac{1}{W_{a}}\sum_{j\in\mathcal{B}_{a}\setminus\mathcal{C}_{a}}w_{j}\ge1-\frac{V\exp(-\beta\rho)}{W_{a}}.\] Since \(W_{a}\ge w_{i_{a}}=\exp(-\beta F_{i_{a}})\) and \(\rho>F_{i_{a}}\) (the barrier between basins is strictly higher than their minima), the fraction \(\frac{V\exp(-\beta\rho)}{W_{a}}\le V\exp(-\beta(\rho-F_{i_{a}}))\) tends exponentially to zero as \(\beta\to\infty\). Denote this correction as \(\delta(\beta)=o(1)\). Then the numerator of the variational functional (the square of this difference) for the trial vector \(\mathbf{h}^{*}\) is bounded below as \((1-\delta(\beta))^{2}=1-o(1)\). The vector \(\mathbf{h}^{*}\) may not satisfy the condition \(\langle\mathbf{w},\mathbf{h}^{*}\rangle=0\). Introduce a shift \(\tilde{\mathbf{h}}=\mathbf{h}^{*}-c\mathbf{1}\), where \(c=\langle\mathbf{w},\mathbf{h}^{*}\rangle/\langle\mathbf{w},\mathbf{1}\rangle\). The numerator is invariant with respect to a shift by a constant, since \(\langle\mathbf{p}_{a}-\mathbf{p}_{b},\mathbf{1}\rangle=1-1=0\). The quadratic form is also invariant: \(\mathcal{E}(\tilde{\mathbf{h}},\tilde{\mathbf{h}})=\mathcal{E}(\mathbf{h}^{*},\mathbf{h}^{*})\). Substituting \(\tilde{\mathbf{h}}\) into Lemma 1, we obtain: \[D^{2}_{\beta}\ge\frac{1-o(1)}{|\mathcal{E}|\exp(-\beta\rho)}=\frac{1-o(1)}{|\mathcal{E}|}\exp(\beta\rho).\] Consequently, \(\liminf_{\beta\to\infty}\frac{1}{\beta}\ln D^{2}_{\beta}\ge\rho\).

Step 2 (finding the upper bound). Let \(\mathbf{h}\in\mathbb{R}^{V}\) be an arbitrary vector satisfying the condition \(\langle\mathbf{w},\mathbf{h}\rangle=0\). Denote \(\Delta=\langle\mathbf{p}_{a},\mathbf{h}\rangle-\langle\mathbf{p}_{b},\mathbf{h}\rangle\). Let \(i_{a}\in\mathcal{B}_{a}\) and \(i_{b}\in\mathcal{B}_{b}\) be the vertices with minimum \(F\) values in their basins. Split the sum \(\langle\mathbf{p}_{a},\mathbf{h}\rangle=\sum_{j\in\mathcal{B}_{a}}(p_{a})_{j}h_{j}\) into two parts: over the set \(\mathcal{B}_{a}\cap\mathcal{C}_{a}\) and over the set \(\mathcal{B}_{a}\setminus\mathcal{C}_{a}\). For \(j\in\mathcal{B}_{a}\cap\mathcal{C}_{a}\), there exists a path \(\gamma_{j}\) to \(i_{a}\) inside \(\mathcal{C}_{a}\) of length at most \(V\). The maximum edge weight inside the finite set \(\mathcal{C}_{a}\) is strictly less than \(\rho\). Denote \(\rho_{a}=\max_{\{u,v\}\in\mathcal{E}(\mathcal{C}_{a})}U_{uv}<\rho\). By the Cauchy–Bunyakovsky inequality, we have \[|h_{j}-h_{i_{a}}|\le\sqrt{V\exp(\beta\rho_{a})\mathcal{E}(\mathbf{h},\mathbf{h})}.\] For \(j\in\mathcal{B}_{a}\setminus\mathcal{C}_{a}\), as proven in Step 1, \(F_{j}\ge\rho\). The gradient descent path to \(i_{a}\) has maximum weight \(F_{j}\). Consequently, \(|h_{j}-h_{i_{a}}|\le\sqrt{V\exp(\beta F_{j})\mathcal{E}(\mathbf{h},\mathbf{h})}\). The contribution of these vertices to the mean value is estimated as \[(p_{a})_{j}|h_{j}-h_{i_{a}}|=\frac{\exp(-\beta F_{j})}{W_{a}}\sqrt{V\exp(\beta F_{j})\mathcal{E}(\mathbf{h},\mathbf{h})}=\frac{\sqrt{V\mathcal{E}(\mathbf{h},\mathbf{h})}}{W_{a}}\exp(-\beta F_{j}/2).\] Summing over all \(j\) and taking into account that \(W_{a}\ge\exp(-\beta F_{i_{a}})\), we obtain that the contribution of the set \(\mathcal{B}_{a}\setminus\mathcal{C}_{a}\) is bounded by \(\frac{V^{3/2}\sqrt{\mathcal{E}(\mathbf{h},\mathbf{h})}}{W_{a}}\exp(-\beta\rho/2)\le V^{3/2}\exp(\beta(F_{i_{a}}-\rho/2))\sqrt{\mathcal{E}(\mathbf{h},\mathbf{h})}\). Since \(\rho_{a}<\rho\) and \(F_{i_{a}}-\rho/2<\rho/2\), both terms are exponentially smaller than the main term \(\exp(\beta\rho/2)\sqrt{\mathcal{E}(\mathbf{h},\mathbf{h})}\). Consequently: \[\langle\mathbf{p}_{a},\mathbf{h}\rangle=h_{i_{a}}+\varepsilon_{a},\;\text{where }|\varepsilon_{a}|=o\!\left(\exp(\beta\rho/2)\sqrt{\mathcal{E}(\mathbf{h},\mathbf{h})}\right).\] Similarly, for \(\mathcal{B}_{b}\), we have \(\langle\mathbf{p}_{b},\mathbf{h}\rangle=h_{i_{b}}+\varepsilon_{b}\) with an exponentially suppressed error \(\varepsilon_{b}\).

Now estimate the difference \(h_{i_{a}}-h_{i_{b}}\). Construct a composite path from \(i_{a}\) to \(i_{b}\) by joining:

(1) The path reverse to the gradient descent path from some vertex \(p\in\mathcal{B}_{a}\) to \(i_{a}\). Since energies decrease along the descent, the maximum edge weight on this path is \(F_{p}\).

(2) The optimal path \(\gamma^{*}\) from \(p\) to \(q\in\mathcal{B}_{b}\), whose height by definition is \(\rho\).

(3) The gradient descent path from \(q\) to \(i_{b}\) with maximum weight \(F_{q}\).

It is important to note that \(F_{p}\le\rho\) and \(F_{q}\le\rho\) (otherwise any path from \(p\) or \(q\) would have height strictly greater than \(\rho\), which contradicts the existence of \(\gamma^{*}\)). Consequently, the maximum edge weight on this entire composite path is \(\max(F_{p},\rho,F_{q})=\rho\). Since this composite path may contain self-intersections, extract from it a simple path \(\gamma_{\mathrm{simple}}\) from \(i_{a}\) to \(i_{b}\) by removing all cycles. The length of this simple path \(l\) is strictly less than the number of vertices in the graph: \(l\le V-1<V\). The maximum edge weight on \(\gamma_{\mathrm{simple}}\) does not exceed the maximum edge weight on the original composite path, i.e., it does not exceed \(\rho\). Since the simple path does not contain multiple edges, applying the Cauchy–Bunyakovsky inequality and taking into account that \(C_{ij}\ge\exp(-\beta\rho)\) for all edges of the path, we obtain a rigorous estimate \[(h_{i_{a}}-h_{i_{b}})^{2}\le l\sum_{\{i,j\}\in\gamma_{\mathrm{simple}}}(h_{i}-h_{j})^{2}\le V\exp(\beta\rho)\sum_{\{i,j\}\in\gamma_{\mathrm{simple}}}C_{ij}(h_{i}-h_{j})^{2}.\] Since each edge in \(\gamma_{\mathrm{simple}}\) occurs exactly once, the last sum is strictly majorized by the global Dirichlet quadratic form: \[(h_{i_{a}}-h_{i_{b}})^{2}\le V\exp(\beta\rho)\mathcal{E}(\mathbf{h},\mathbf{h}).\] Now write \(\Delta=\langle\mathbf{p}_{a},\mathbf{h}\rangle-\langle\mathbf{p}_{b},\mathbf{h}\rangle=(h_{i_{a}}-h_{i_{b}})+\varepsilon_{a}-\varepsilon_{b}\). Using the obtained estimate for \((h_{i_{a}}-h_{i_{b}})\) and the exponential smallness of the errors \(\varepsilon_{a},\varepsilon_{b}\), we find: \[\Delta^{2}\le3V\exp(\beta\rho)\mathcal{E}(\mathbf{h},\mathbf{h})\cdot(1+o(1)).\] Since this holds for any admissible \(\mathbf{h}\), the maximum in Lemma 1 satisfies the estimate \[D^{2}_{\beta}\le3V\exp(\beta\rho)\bigl(1+o(1)\bigr).\] Consequently, \(\limsup_{\beta\to\infty}\frac{1}{\beta}\ln D^{2}_{\beta}\le\rho\) (the factor \(3V\) disappears after dividing by \(\beta\) and taking the logarithm as \(\beta\to\infty\)).

Step 3 (completion of the proof). From the obtained bounds it follows that \(\lim_{\beta\to\infty}\frac{1}{\beta}\ln D^{2}_{\beta}=\rho\). Since \(D_{\beta}>0\) and \(\ln D_{\beta}=\frac{1}{2}\ln D^{2}_{\beta}\), we finally obtain: \[\lim_{\beta\to\infty}\frac{1}{\beta}\ln D_{\beta}(\mathcal{B}_{a},\mathcal{B}_{b})=\frac{1}{2}\,\rho^{(0)}(\mathcal{B}_{a},\mathcal{B}_{b}),\] which proves the Theorem.

8 Discussion↩︎

The theorem formulated in Section 7 allows us to explain the dependence of \(u_{\mathrm{nt}}\) on \(\Delta F\) observed in the computational experiment, presented in Table 1. For small \(\Delta F\) (e.g., \(\Delta F=10\) and \(20\) kJ/mol), the energy differences between neighboring vertices are comparable to the thermal energy \(RT\). In this regime, the Kramers matrix is not sharply inhomogeneous: transition probabilities through different edges differ by no more than a factor of \(\exp(8)\sim10^{3}\). Although this creates a certain hierarchy, it is not absolutely rigid. In particular, alternative paths between two basins passing through different saddle points can have comparable probabilities, and the resulting kinetic distance is determined by the superposition of contributions from several paths. Moreover, in this regime, finite-size effects and differences in basin depths give corrections of the same order of magnitude as the barrier term \(\frac{1}{2}\rho^{(0)}\), further violating ultrametric relations. The combination of these factors leads to a reduction in the fraction of nontrivially ultrametric triples (\(\approx42\%\) and \(\approx61\%\), respectively). As \(\Delta F\) increases, the ratio of the characteristic energy difference to \(RT\) increases. At \(\Delta F=50\) kJ/mol (\(\approx20\,RT\)), the difference in transition rates through different edges can reach \(\exp(20)\sim10^{9}\) times, and \(u_{\mathrm{nt}}\) increases to \(\approx81\%\). At \(\Delta F=200\) kJ/mol (\(\approx80\,RT\)), the difference in transition rates through different edges can reach \(\exp(80)\sim10^{35}\) times. In this regime, the transition probability through the highest edge on the path exponentially dominates the probabilities of all other transitions, and the kinetics begins to be determined by a single optimal path – the one that minimizes the maximum energy along the path. This corresponds to the construction of the single-linkage metric \(\rho^{(0)}\), which by construction is an ultrametric. Simultaneously, for large \(\Delta F\), the characteristic energy amplitude significantly exceeds the heights of suboptimal barriers, so that corrections associated with alternative paths become negligibly small compared to the barrier term \(\frac{1}{2}\rho^{(0)}\). In this regime, the kinetic metric asymptotically approaches a strict ultrametric. With a further increase in \(\Delta F\) to \(500\) and \(1000\) kJ/mol (\(\approx200\,RT\) and \(\approx400\,RT\)), the dominance of the optimal path and the suppression of contributions from suboptimal paths become even more pronounced, and \(u_{\mathrm{nt}}\) continues to grow monotonically, reaching \(90.6\%\) and \(95.9\%\), respectively.

Thus, the physical mechanism for the increase in ultrametricity consists of the joint action of two factors: the exponential suppression of contributions from suboptimal paths to the effective transition rate between basins, and the tendency of the system towards a limit regime where kinetics are fully determined by minimax barriers. As \(\Delta F\to\infty\), only the optimal path (in the minimax sense) gives a non-zero contribution, and the kinetic metric asymptotically converges to the exact single-linkage ultrametric. For finite \(\Delta F\), both suboptimal paths and finite-size effects (entropic pre-exponential factors) contribute, violating exact ultrametricity. However, for a graph of fixed size \(V\), the magnitude of these corrections to the minimax barrier is exponentially suppressed with increasing ratio \(\Delta F/(RT)\), as it is determined by the Arrhenius factor for the difference in heights between the optimal and suboptimal barriers \(\sim\exp(-\delta\cdot\Delta F/RT)\). Note that since the parameters \(\Delta F\) and \(T\) enter the theory only through the combination \(\Delta F/(RT)\), all the above conclusions about the increase in the degree of ultrametricity with increasing \(\Delta F\) at a fixed temperature apply equally to the low-temperature limit \(T\to0\) at fixed \(\Delta F\).

It is important to emphasize the role of graph topology in the described mechanism. A sparse Erdős–Rényi graph with \(\langle k\rangle=2.0\) is locally tree-like: for any vertex, its \(r\)-neighborhood for \(r=O(\log\log V)\), with probability tending to one as \(V\to\infty\), contains no cycles. As a result, the probability of the existence of several substantially different simple paths between two distant vertices is low. On an ideal tree, there is exactly one simple path between any two vertices, and the single-linkage metric automatically becomes an ultrametric. The presence of rare cycles in an Erdős–Rényi graph creates alternative paths, which, for finite \(\Delta F\), can give a noticeable contribution and thereby reduce the observed degree of ultrametricity. However, as \(\Delta F\) increases, even these alternative paths become separated by characteristic barrier heights, and the hierarchy is restored. Thus, the very topology of the Erdős–Rényi graph induces a hierarchical structure, and the large energy spread makes this structure kinetically explicit, while simultaneously suppressing the influence of suboptimal paths.

It should be emphasized that the limit theorem proved in Section 7 does not use assumptions about the average graph degree or its local tree-likeness and is formally valid for any finite connected graph, including a complete graph. However, for \(\langle k\rangle\gg1\), there exists a large number of alternative paths between vertices, and for finite \(\Delta F\), the contributions of many of them to the kinetic metric may be comparable. In this case, the convergence of \(u_{\mathrm{nt}}\) to unity with increasing \(\Delta F\) may be achieved at significantly larger values of \(\Delta F/(RT)\) than in the sparse case. Thus, local tree-likeness, characteristic of an Erdős–Rényi graph with \(\langle k\rangle=O(1)\), is not a necessary condition for asymptotic ultrametricity, but ensures its manifestation at moderate values of \(\Delta F\) accessible for direct numerical observation.

The obtained results establish a connection between the kinetic Mahalanobis metric and the \(p\)-adic parametrization of ultrametric spaces used in [18]–[21]. In \(p\)-adic models, the state space is assumed to be isomorphic to the set of leaves of a regular \(p\)-adic tree, and the distance between states is given by an explicit ultrametric. The limit theorem proved in Section 7 shows that for a locally tree-like graph with an exponentially wide energy distribution, the kinetic metric asymptotically converges to the single-linkage ultrametric without a priori postulating a tree structure. Deviations from strict ultrametricity for finite \(\Delta F\) are associated with the presence of suboptimal paths, the contribution of which is exponentially suppressed as \(\Delta F/(RT)\) increases. This result indicates that \(p\)-adic diffusion equations may be applicable to systems with locally tree-like topology and strong energy inhomogeneity, going beyond the biopolymer applications considered in [18]–[21].

Note that the obtained results also define a number of directions for further research. First, a quantitative estimate of the contribution of suboptimal paths to \(u_{\mathrm{non}}\) at finite \(\Delta F\) and its dependence on the parameters \(\Delta F\), \(T\), and the average graph degree \(\langle k\rangle\) is of interest. Second, the model used in this work assumes a uniform distribution of vertex free energies; generalization to other distributions (e.g., normal or exponential) will allow establishing whether the form of the distribution affects the rate of convergence of \(u_{\mathrm{nt}}\) to unity with increasing \(\Delta F\). Third, it is of interest to apply the developed formalism to random graphs with a higher average degree (\(\langle k\rangle>2\)), as well as to graphs containing a pronounced modular structure, which will allow estimating the influence of the number of cycles and topological inhomogeneity on the degree of ultrametricity for fixed energy landscape parameters.

References↩︎