Let \( x_i=(r_{i1},r_{i2},…,r_{iN})\) denote the \( 3N\) -dimensional vector containing the three-dimensional Cartesian coordinates \( r_{i\alpha}\) of the \( N\) atoms in molecular structure \( i\) . Given a normalized atom weight vector \( \boldsymbol{w} = \{w_\alpha\}\;(\alpha = 1, …, N)\) with \( \sum_\alpha w_\alpha = 1\) and \( w_\alpha \ge 0\) , the weighted mean square deviation (MSD) between two structures after optimal rigid-body superposition is given by

Here \( R \in SE (3)\) denotes a rigid-body transformation. As a member of the special Euclidean group, \( R\) consists of a translation (of the entire structure) and a rotation (about its center, as defined by the weighted average of atom positions with weights \( \alpha\) ). \( |\cdot|\) represents the 3D Euclidean distance. The optimal \( R\) for a given weight vector \( \boldsymbol{w}\) is found by the Kabsch algorithm.[2, 3] The root-mean-square deviation is \( \text{RMSD} = \text{MSD}^{1/2}\) . Classical RMSD corresponds to the uniform prior \( w_\alpha = W_\alpha \equiv 1/N\) for all \( \alpha\) .

In a Bayesian formulation of RMSD structural alignment, we assume a priori that configurations \( x\) in \( 3N\) -dimensional configuration space are distributed as

where \( s\) is the reference structure and \( \sigma\) is the expectation of the root-mean-squared fluctuation per atom. As a Bayes prior for the weight vector, we assume that the optimal atom weight vector \( \boldsymbol{w}\) entering the \( MSD \) is distributed as

where \( \boldsymbol{W}\) is a vector of reference weights and \( \theta\) is a confidence parameter. Analogous to Bayesian ensemble refinement [12], we then use the Bayes relation to define a Bayes posterior for \( s\) and \( \boldsymbol{w}\) given configurations \( x_1,x_2,…\) ,

with an effective free energy

The optimal reference and weight vector are then defined by maximizing the Bayes posterior, or equivalently by minimizing the effective free energy,

Maximizing the posterior is unchanged if we multiply \( g\) by the constant \( 2\sigma^2\) . Setting the prior confidence to \( \theta = M/2\) fixes the entropy penalty at a per-frame weight \( \sigma^2\) , so its total weight becomes \( M\sigma^2\) . With this choice the confidence parameter is no longer free, and we minimize the effective free energy

with respect to \( s\) and \( \boldsymbol{w}\) . In the following, we use \( \sigma\) or

as fully equivalent parameters to control the trade-off between minimizing the \( MSD \) and the KL divergence. From here on \( \theta\) denotes \( M\sigma^2\) , not the prior confidence of eq 3. With this Bayesian perspective in mind, we refer to the RMSD alignment procedure using the optimal reference structure \( s\) and optimal weights \( \boldsymbol{w}\) as “Bayesian inference of RMSD” (BRMSD).

We note the close analogy to an expression derived originally in the “ensemble refinement of SAXS” (EROS) method [11] and then worked out in detail in ref [12]. Here, the MSD term replaces the \( \chi^2\) error of ensemble refinement. As in ensemble refinement, an expression identical to eq 8 is obtained by combining the usual metric of MSD with a relative entropy regularizer to penalize non-uniform atom weight distributions.[10, 11, 12]

To find the weighted average structure \( s\) of an ensemble of \( M\) structures \( \{x_i\}\) , we minimize \( G\) over the average structure \( s\) , the rigid-body superpositions \( \{R_i\}\) of \( x_i\) and \( s\) , the reference structure \( s\) , and the atom weights \( \boldsymbol{w}\) , subject to the normalization condition of the weights. The first term in \( G\) , as given in eq 8, is the aggregate mean-square deviation from the average structure. The second term represents the relative entropy of \( \boldsymbol{w}\) with respect to the uniform prior \( W_\alpha = 1/N\) , scaled by \( \theta = M\sigma^2\) . The minimization is achieved by setting derivatives of \( G\) with respect to \( s_\alpha\) and \( w_\alpha\) to zero. This yields the self-consistency equations

where \( R_i^{(n)}\) is the optimal rigid-body translation and rotation of \( x_i\) for weights \( \boldsymbol{w}^{(n)}\) with \( s^{n}\) as reference, \( R_i^{(n)}r_{i\alpha}\) its application to atom \( \alpha\) of \( x_i\) , and the weights \( w_\alpha^{(n)}\) are normalized after each iteration step \( n\) . Starting from \( s^{(1)} = x_1\) and \( w_\alpha^{(1)} = W_\alpha\) , eqs 11 and 10 are iterated towards convergence.

\( G\) has a lower bound of zero since the MSD term is a sum of squared norms and \( S_\mathrm{KL}(\boldsymbol{w}\|\boldsymbol{W})\ge 0\) by Gibbs’ inequality[10]. Each alternating minimization step decreases \( G\) monotonically. The weight update minimizes \( G\) exactly at fixed \( s\) , and the structure update minimizes the MSD term at fixed \( \boldsymbol{w}\) . Such a decrease of a function bounded from below guarantees convergence to a stationary point. Since \( G\) is non-convex, this point may be a local, rather than a global minimum. Therefore, although for alignment and smoothing, a single initialization from \( s^{(1)} = x_1\) and \( w_\alpha^{(1)} = W_\alpha\) may suffice, for clustering and domain identification, we run \( n_r\) random restarts and retain the solution with the lowest \( G\) .

The sum in the exponent of eq 10 is the per-atom mean-square fluctuation (MSF) about the current average. Accordingly, atoms with large fluctuations receive exponentially suppressed weights, while structurally rigid atoms are up-weighted. Atoms whose summed squared deviation \( \sum_i |s_\alpha - R_i r_{i\alpha}|^2\) exceeds \( \theta = M\sigma^2\) are exponentially suppressed, whereas those below retain near-uniform weight. Equivalently, an atom is suppressed once its mean-square fluctuation averaged over the \( M\) frames exceeds \( \sigma^2\) , so the weights are set by the per-atom mean-square fluctuation rather than by the number of frames \( M\) used to estimate it, provided the sampling is converged. At \( \sigma \to \infty\) the entropy term dominates and weights stay uniform, recovering classical RMSD. The optimal \( \sigma\) depends on the module and on the fluctuation amplitude of the system.

The entropy-based effective atom count is the exponential of the Shannon entropy of the weight distribution (eq 12),

where, \( n_\text{eff}\) equals \( N\) for uniform weights and decreases monotonically as weight concentrates.

The BRMSD alignment concentrates weight on the most globally rigid regions of the protein, such as the TMD groove in SND3. However, the most functionally relevant regions may not always coincide with the rigid structural core. For instance, a flexible domain critical for substrate binding will warrant non-negligible weight alongside the globally rigid scaffold. Hence, depending on the problem at hand, one may wish to focus the alignment on a given domain of interest \( D\) while still including the remainder of the structure. We augment the free-energy functional of eq 8 with a concentration penalty on atoms in \( D\) (eq 13).

where \( n_D\) is the number of atoms in \( D\) and \( \mu > 0\) controls the strength of the concentration bias toward \( D\) . Taking the derivative with respect to \( w_\alpha\) yields two separate update rules.

Atoms outside \( D\) follow the standard alignment update (with equivalent eq 14 and eq 10), whereas atoms inside \( D\) experience a reduced effective concentration parameter \( \theta + \mu\) and an additional prior bias \( W_\alpha^{\theta/(\theta+\mu)}\) that amplifies their weight relative to atoms outside \( D\) . In the limit \( \mu \gg \theta\) , the atoms in \( D\) receive weight \( w_\alpha \propto n_D^{-1}\) (uniform within \( D\) ) regardless of their MSF. As before, all weights are renormalized after each update, and the structure update follows eq 11.

The dimensionless ratio \( \mu/\theta\) controls the degree of domain focusing. A practical criterion for selecting the operating point \( (\mu/\theta)_\text{op}\) is that \( n_\text{eff}\) should remain above \( \lvert D \rvert\) . If \( n_\text{eff} < \lvert D \rvert\) , the alignment is effectively performed on fewer atoms than the domain contains, which can cause numerical instability.

Structural alignment averages uniformly over all \( M\) frames. High-frequency positional noise from thermal fluctuations is thereby retained in each individual frame, even after alignment. To suppress this noise while preserving slow conformational dynamics, we replace the global average with a local one centered on each frame \( j\) . We introduce frame-dependent weights \( p(i|j) \ge 0\) , \( \sum_i p(i|j) = 1\) , that localize structure \( i\) relative to frame \( j\) . A triangular window of half-width \( W\) sets \( p(i|j) \propto \max(0, W - |i-j|)\) or a uniform window can be used. The smoothed-frame functional replaces eq 8 by

minimized independently for each frame \( j\) . The fixed-point equations are identical in form to eqs 11–10 with the uniform sum replaced by the windowed sum \( \sum_i p(i|j)\) . Because the window weights are normalized (\( \sum_i p(i|j) = 1\) ), this windowed sum is a per-frame average, so the entropy penalty here uses \( \theta = \sigma^2\) (the effective frame count is \( M = 1\) ). The converged \( s_j\) is the smoothed coordinate of frame \( j\) , and the collection \( \{s_j\}\) constitutes the smoothed trajectory.

Structural alignment identifies a single average structure and a single weight vector for the full ensemble. When the ensemble spans multiple conformational states, a single average is insufficient. We extend the alignment framework to \( K\) cluster centers by replacing the uniform frame average with a soft assignment. Each frame \( i\) is assigned to cluster \( a\) with responsibility[25] \( q(a|i) \ge 0\) , \( \sum_a q(a|i) = 1\) . Each cluster carries its own center \( s_a\) and, crucially, its own atom weight vector \( \boldsymbol{w}_a\) , so the alignment metric adapts to the internal rigidity pattern of each conformational state rather than imposing a single global metric. The joint free-energy functional is given by eq 17.

Here, the third term penalizes deterministic assignments (large \( \tau\) encourages soft, diffuse responsibilities, small \( \tau\) recovers hard \( K\) -means). Setting \( \partial G/\partial q(a|i) = 0\) subject to \( \sum_a q(a|i) = 1\) via a Lagrange multiplier yields the softmax responsibility update (eq 18), where \( s_a\) is the center of cluster \( a\) , \( \theta\) controls atom weight concentration (eq 9), and \( \tau > 0\) (units Å\( ^2\) ) controls cluster assignment softness. The fixed-point iteration is

The cluster center \( s_{a\alpha}\) is the responsibility-weighted average structure and represents the fixed point that minimizes the weighted MSD across all frames. Each cluster has its own weight vector \( \boldsymbol{w}_a\) , so the alignment metric adapts to each cluster’s internal rigidity pattern.

Conformational clustering partitions frames into groups (clusters), whereas rigid domain identification partitions atoms into groups (domains). Both share the BRMSD free-energy structure, but the variable of interest shifts from frame assignments \( q(a|i)\) to atom weight vectors \( \boldsymbol{w}_k\) .

The identification of rigid domains is done by sequential peeling of atoms. A series of independent optimizations with \( K=1\) domain is performed on shrinking atom pools. At each step, the BRMSD functional (eq 8) is minimized on the current pool to obtain a single weight vector \( \boldsymbol{w}\) concentrated on the most rigid sub-structure within that pool. Atoms with \( w_\alpha > t \cdot \max_\beta w_\beta\) are assigned to the current domain and removed from the pool. The threshold \( t\) controls pool composition for subsequent steps. Small values over-claim atoms, depleting later pools, and large values underclaim, leaving rigid residues in the pool that can dominate subsequent steps. The procedure repeats on the remaining atoms until no atoms satisfy the threshold or \( K\) domains have been recovered. The peeling order is data-determined, such that the most globally rigid sub-structure is identified first. Hence, domains emerge in order of decreasing rigidity. At each step, \( \theta = M\sigma^2\) is held fixed and independent of the pool size, preserving the physical meaning of \( \sigma\) as the per-atom root-mean-squared fluctuation (RMSF) threshold across all steps. Domain detection quality is assessed by the Jaccard index \( J(A,B) = |A \cap B|/|A \cup B|\) between each recovered domain and the reference assignment.

An alternative formulation seeks \( K\) non-overlapping weight vectors simultaneously, enforcing mutual exclusivity through an overlap penalty on the dot products \( \boldsymbol{w}_j \cdot \boldsymbol{w}_k\) .

The weight update becomes

This update is obtained by differentiating \( G\) with respect to \( w_{k\alpha}\) , adding a Lagrange multiplier for the simplex constraint \( \sum_\alpha w_{k\alpha} = 1\) , and setting the derivative to zero. The overlap gradient \( \mu\sum_{j\ne k} (\boldsymbol{w}_j\cdot\boldsymbol{w}_k)w_{j\alpha}\) enters additively in the exponent, coupling domain \( k\) ’s weights to all other domains. Since the coupling is simultaneous, the weight updates for different \( k\) are not independent. In practice, simultaneous optimization can fail when the overlap-penalty gradient is too weak to prevent weight vectors from collapsing onto the globally most rigid region (see Results). Sequential peeling resolves this by construction.

Implementation details, including timing benchmarks for all modules (Table 1), are provided in the Supporting Information.

For completeness, we also sketch an extension to a Mahalonobis-like RMSD that accounts for anisotropic fluctuations in atom positions. Let us assume that we have an initial alignment, which gives us the aligned structures \( x_i\) . In a form of PCA,[26] we then construct the average \( s=\sum_{i=1}^Mx_i/M\) of the \( M\) structures and the covariance matrix \( C=\sum_{i=1}^M\delta x_i\delta x_i^T/M\) with \( \delta x_i=x_i-s\) and \( T\) the transpose. A spectral decomposition of \( C\) (or a singular value decomposition of the \( x_i\) ) then gives us \( n\) eigenvalues \( \lambda_i>0\) with orthonormal eigenvectors \( v_i\) , i.e., \( v_i^Tv_j=\delta_{ij}\) with \( \delta_{ij}\) the Kronecker delta. At least six of the eigenvalues will be zero, corresponding to rigid body translations and rotations, \( n\le 3N-6\) . We define the free energy by giving weights \( w_j\) to the PCA modes \( v_j\) as

where \( \boldsymbol{R}=I_N\otimes R\) applies the translation and rotation atom-wise to the configurations. Minimization with respect to \( s\) gives the average position, \( s=\sum_{i=1}^M\boldsymbol{R}x_i/M\) , as before. Setting the derivative with respect to the normalized weights \( w_j\) to zero gives us an update rule for the mode weights,

However, we have not currently implemented this procedure because the rigid-body alignment code to optimize \( R\) will have to be adapted to the more general metric. After iteration to self-consistency, one would then interpolate between an alignment according to the Mahalanobis distance (\( \theta=0\) ) and the usual Euclidean distance (\( \theta\to\infty\) ).

We benchmarked BRMSD on two systems that represent complementary challenges. For alignment and smoothing, we used SND3, which is a 191-residue eukaryotic translocon-associated protein recently resolved by cryo-EM at 3.1 Å resolution (PDB: 9I78).[19] Three transmembrane helices (TMD-2, residues 28–60; TMD-3, residues 86–115; TMD-4, residues 119–134) form a buried TM groove that constitutes the rigid functional core. An additional, shorter transmembrane helix (TMD-1, residues 5–21) is embedded in the membrane and is less rigid than the other three TMDs (it serves as the target domain for the focused alignment benchmark). Two cytosolic claw domains (N-claw, residues 61–85; C-claw, residues 136–173) and a C-terminal tail (residues 174–191) are substantially more mobile. We analyzed a 300 ns all-atom MD trajectory of the isolated protein chain, sampled at 100 ps intervals. The SND3 trajectory samples equilibrium fluctuations around a single mean structure without discrete conformational states.

The conformational clustering and domain identification were benchmarked on adenylate kinase (ADK), a phosphotransferase that undergoes a large-scale open-to-closed conformational transition upon substrate binding [21]. Crystallography and mutational studies have established three structurally distinct regions in ADK to single-residue precision, namely the CORE (residues 1–29, 60–121, and 160–214), NMPbind (residues 30–59), and LID (residues 122–159) domains.[21] The two endpoint crystal structures, open (PDB ID: 4AKE) and closed (PDB ID: 1AKE), are characterized by NMPbin–LID distances of 32.1 and 19.6 Å, and LID–CORE distances of 40.0 and 30.0 Å, respectively. To evaluate clustering performance, we utilized 200 independent Dynamic Importance Sampling (DIMS) trajectories [27] (with a total of 19,691 frames) that span this open-to-closed transition. This ensemble represents a non-equilibrium transition rather than a continuous equilibrium trajectory. For domain detection validation, we employed 28,282 Framework Rigidity Optimized Dynamics Algorithm (FRODA) frames [28]. FRODA propagates rigid-body constraints on covalent geometry, so intra-domain RMSF is substantially lower than inter-domain RMSF by construction.

We scanned \( \sigma\) from 0.1 to 10 Å on the SND3 trajectory (Figure 1). As \( \sigma\) decreases, the per-residue weights \( w_{\rm {residue}}\) grow monotonically in the TM groove (TMD-2, TMD-3, TMD-4) while collapsing toward zero in the more mobile regions (Figure 1A). The per-residue RMSF profile is mostly insensitive to \( \sigma\) (Figure 1B), confirming that the weight distribution robustly identifies the same rigid core regardless of regularization strength.

The role of \( \sigma\) is not to alter which residues are mobile, but to control the sharpness of the core-periphery boundary used for alignment and domain decomposition. Larger \( \sigma\) yields a diffuse weighting that spreads weight across many residues, while smaller \( \sigma\) sharpens the boundary, concentrating weight on the most rigid atoms. However, excessively small \( \sigma\) over-regularizes the distribution, collapsing the effective alignment set to too few atoms and risking sensitivity to local structural noise rather than global domain motion.

Pan/Zoom

The optimal \( \sigma\) is chosen from two diagnostics plotted in Figure 2A, namely the ratio \( w_\text{groove}/w_\text{rest}\) with \( w_\text{D} = \sum_{\alpha\in{\rm {D}}}w_\alpha\;(D={\rm groove, rest})\) (shown in red) and the effective atom count \( n_\text{eff}\) (eq 12) (shown in green). We define \( \sigma_\text{op}\) as the smallest \( \sigma\) at which \( n_\text{eff} \geq 0.20\,N\) . At this value, at least one-fifth of all atoms contribute to the alignment and weight concentration on the rigid core is near its practical maximum. At \( \sigma_\text{op} = 0.43\)  Å, the ratio reaches 6.5 and \( n_\text{eff} = 604\) , so the groove residues carry the dominant weight while a statistically meaningful fraction of the chain (604 of 3013 atoms) still contributes to the alignment. The cumulative weight-MSF Lorenz curve at \( \sigma_\text{op}\) quantifies this discrimination. The 50% most-rigid atoms carry approximately 90% of the total weight (Figure 6).

Pan/Zoom

The weight optimization automatically concentrates weight on the TM groove without manual residue selection. The resulting BRMSD (\( 0.50\pm0.08\)  Å) at \( \sigma_\text{op}=0.43\)  Å is \( \sim\) 3 times smaller than the standard RMSD (\( 1.64\pm0.25\)  Å, uniform weights) (Figure 7). The domain-decomposed RMSD time series at \( \sigma_\text{op}\) (Figure 2B) shows that the TM-groove RMSD is substantially lower and less variable than both the non-groove domain (Rest) and the standard all-atom RMSD, confirming that weight concentration on rigid atoms suppresses the dominant source of conformational noise. Per-residue RMSF and atom weights at \( \sigma_\text{op}\) , shown in Figure 8, demonstrate the exponential suppression of flexible domains. The BRMSD-aligned trajectory is shown in Movie S1 alongside the raw trajectory.

To demonstrate the focused alignment module, we selected TMD-1 in SND3. Although it is membrane-embedded and structurally ordered, it contributes marginally to the global alignment (\( w_\text{TMD1} = 0.002\) at \( \sigma_\text{op}\) ) because its fluctuations are larger than those of TMD-2/3/4. Focused alignment recovers the internal dynamics of TMD-1 by explicitly biasing the weight optimization toward the TMD-1 atoms.

We scanned \( \mu/\theta\) ratio from 0 to 0.8 at \( \sigma = 0.43\)  Å (Figure 9). Figure 3A shows the per-residue normalized weight profiles for global alignment at \( \sigma_\text{op}\) (blue) and focused alignment at \( \mu/\theta = 0.53\) (orange). The latter raises the per-residue weights sharply in TMD-1 while the TM-groove weight decreases and the non-groove rest remains suppressed. As \( \mu/\theta\) increases, \( w_\text{TMD1}\) fraction grows from 0.002 to 0.54 while the \( w_\text{groove}\) fraction decreases from 0.87 to 0.44 and \( w_\text{rest}\) fraction falls to 0.02 (Figure 3B). The TMD-1 mean RMSF decreases from 1.49 Å at global alignment (\( \mu/\theta = 0\) ) to 1.14 Å with focused alignment at \( \mu/\theta = 0.53\) (Figure 3C). This 24% reduction reflects the removal of rigid-body translations and rotations of TMD-1 relative to the rest of the protein. The per-residue weight and RMSF profiles for all \( \mu/\theta\) values are shown in Figure 9A.

Pan/Zoom

The operating point \( (\mu/\theta)_\text{op} = 0.53\) is identified from the \( n_\text{eff}\) criterion. At this value \( n_\text{eff} = 303\) , which just exceeds the TMD-1 atom count \( \lvert D \rvert = 281\) (Figure 3B). At \( \mu/\theta = 0.55\) , \( n_\text{eff}\) drops below \( \lvert D \rvert\) (to 125), marking the onset of numerical over-concentration. Beyond the crossover at \( \mu/\theta \approx 0.53\) , the focusing penalty overwhelms the entropy term and almost all weight collapses onto TMD-1 (\( w_\text{TMD1}\to1.0\) , \( w_\text{groove}\to0.002\) ; Figure 3B). The alignment is then performed on TMD-1 alone, so the now-unweighted TM groove and rest of the protein are no longer superimposed and their RMSF rises sharply (Figure 3C). This over-concentration, where \( n_\text{eff}<\lvert D\rvert\) , marks the upper usable limit of \( \mu/\theta\) . The per-\( \mu/\theta\) RMSF profiles in Figure 9B confirm that the groove and rest RMSF values are largely insensitive to the focusing parameter, isolating the TMD-1 RMSF reduction as a genuine alignment effect rather than an artifact.

We smoothed the SND3 trajectory using a triangular window kernel, sweeping the half-width \( W\) over 5, 10, 20, 30, and 50 frames (full window 1–10.0 ns) at \( \sigma_\text{op} = 0.43\)  Å (Figure 10). We quantify the smoothing quality by the mean BRMSD of each raw frame from its locally windowed average. This quantity increases with window width because wider windows produce local averages further from any individual frame, ranging from \( 0.30\)  Å at \( W = 5\) frames to \( 0.37\)  Å at \( W = 50\) frames. We adopt \( W_\text{op} = 20\) frames (4.0 ns full window). The 4.0 ns full window also sits below the slow TM-groove modulation timescale, so it suppresses sub-nanosecond fluctuations without averaging out the slower collective motion. The mean BRMSD from the windowed average at \( W_\text{op}\) is \( 0.35\)  Å.

At \( W_\text{op}\) , the mean BRMSD of the smoothed trajectory from the global average drops from \( 0.50 \pm 0.08\)  Å (raw) to \( 0.35 \pm 0.04\)  Å, a 30% reduction (Figure 11). The smoothed trajectory retains the slow TM groove modulations visible at timescales above 4 ns while suppressing sub-nanosecond fluctuations. Movie S2 shows the raw and smoothed trajectories side by side. Per-residue RMSF in Figure 12 shows that smoothing reduces fast fluctuations uniformly across the chain, while the slow, large-amplitude movements are retained.

Unlike alignment and domain identification, conformational clustering of the open-to-closed transition requires a broad weight distribution rather than a concentrated one. The signal that separates the two states of ADK resides in the mobile NMPbind and LID domains, so concentrating weight on the rigid CORE (small \( \sigma\) ) suppresses precisely those displacements and collapses the cluster centres onto a single point. A sweep of the alignment \( \sigma\) shows that the Calinski–Harabasz (CH) index (eq 26) is almost zero for \( \sigma \lesssim 2\)  Å, where the cluster centers merge. CH index rises sharply between \( \sigma=3\) and \( 4\)  Å and plateaus for larger \( \sigma\) . We adopt \( \sigma=4.0\)  Å as the onset of this stable, high-CH regime (Figure 13).

To select \( \tau\) , we swept over the values 1, 2, 5, 10, 20, and 50 Å\( ^2\) at \( K = 2\) , evaluating the CH index (eq 26), minimum centre separation (eq 27), and population balance (Figure 14). The CH index and centre separation remain high for \( \tau \le 5\)  Å\( ^2\) and collapse to zero at \( \tau \ge 10\)  Å\( ^2\) , where the responsibilities become uniform and the two centres merge. We adopt the largest \( \tau\) that preserves discrimination. At \( \tau_\text{op} = 5.0\)  Å\( ^2\) the two centres are well separated (centre separation \( = 3.20\)  Å) and the partition is most balanced (\( |\bar{q}_1 - 0.5| = 0.02\) ), whereas smaller \( \tau\) over-sharpens toward hard \( K\) -means without improving the partition.

\( K\) -screening at \( \tau_\text{op}\) identifies \( K = 2\) as optimal. At \( K = 3\) , the minimum inter-centre separation collapses to 0.0 Å across restarts, indicating fully degenerate centres. No third distinct conformational state was detectable by BRMSD clustering at these parameters.

The two clusters correspond to open-like (soft population 51.8%, hard count 10,610 frames) and closed-like (soft population 48.2%, hard count 9,081 frames) conformations (Figure 4). The corresponding cluster centres (responsibility-weighted average structures, eq 20) are located at \( d_{\rm {NMP}-LID}\) = 26.2 Å, \( d_{\rm {LID}-CORE}\) = 37.9 Å (open) and \( d_{\rm {NMP}-LID}\) = 18.8 Å, \( d_{\rm {LID}-CORE}\) = 33.7 Å (closed). They are denoted by the black circles in Figure 4A. The near equal populations reflect the pathway-sampling design of the DIMS ensemble and do not correspond to thermodynamic equilibrium free energies.[29] DIMS trajectories sample structural coverage of the open-to-closed transition, not a Boltzmann distribution, so population fractions cannot be interpreted as equilibrium state probabilities.

Pan/Zoom

Plotting each frame by its weighted BRMSD to the two cluster centres (Figure 4B), each distance evaluated in that cluster’s own weight vector, separates the ensemble into two off-diagonal arcs that meet only along the \( d_\text{closed} = d_\text{open}\) diagonal. Because every axis uses a cluster-adapted metric, a frame of one state measured against the other centre is driven to large distance by precisely the mobile LID and NMPbind atoms that distinguish the two states, giving a cleaner separation than the geometric order parameters of Figure 4A. The diagonal crossing marks the transition dividing surface (\( q \approx 0.5\) ), and its low point density indicates rapid transit through the barrier region. The outward curvature of each arc reflects the arced, hinge-bending transition path between the open and closed endpoints rather than a straight interpolation between the two centres.

The responsibilities (\( q(a|i)\) , where \( a\) is the cluster and \( i\) is the structure) correlate strongly with both order parameters. Pearson correlation coefficient of \( q(\text{cluster 1}|i)\) with \( d_\text{NMP-LID}\) and \( d_\text{LID-CORE}\) are \( -0.88\) and \( -0.94\) , respectively (Figure 15). The median of the ratio \( d_\text{between}/d_\text{within}\) (between-cluster and within-cluster distances) is 2.48 and all frames have \( d_\text{between} > d_\text{within}\) . This confirms that the two conformational states are genuinely distinct in BRMSD space (Figure 16). Per-residue backbone displacement between the two cluster centres (Figure 17) reveals the hierarchical displacement order. LID (mean 5.68 Å) \( >\) NMPbind (4.73 Å) \( \gg\) CORE (1.11 Å), consistent with the ground-truth domain assignments.

We first applied the domain identification module to SND3 to assess whether any sub-domain structure is resolvable. Sequential peeling recovers the domains when the alignment \( \sigma\) is set to the rigidity scale of TMD-1. A \( \sigma\) sweep locates the cleanest separation at \( \sigma = 1.0\)  Å (Figure 18), where round 1 recovers the TM groove and round 3 recovers TMD-1 as its own domain (Figure 19) with a Jaccard index of 0.44 to the TMD-1 reference (Figure 20). Round 2 falls on the cytosolic claws but does not form a clean domain (\( J = 0.28\) ). The separation is sharp in \( \sigma\) . It peaks at 1.0 Å and collapses for \( \sigma \geq 1.6\)  Å, where TMD-1 merges into the groove. This value matches the internal RMSF of TMD-1 measured by focused alignment (1.14 Å). Only the stiffest half of TMD-1 is recovered in round 3 because it is a partially rigid membrane-embedded helix rather than a fully independent domain.

For ADK, we applied sequential peeling to the 28,282-frame FRODA backbone ensemble. FRODA generates motions by propagating rigid-body constraints on covalent geometry, so its domain boundaries are structurally unambiguous by construction. Hence the domains are recovered more cleanly in ADK than SND3 (300 ns equilibrium run). However, simultaneous \( K = 3\) optimization fails to detect the three different rigid domains in ADK. This happens because all three weight vectors collapse onto CORE regardless of the overlap penalty \( \mu\) . Every domain begins with uniform weights, so the first iteration concentrates all vectors onto the globally lowest-MSF region (CORE), from which the overlap-penalty gradient is too weak to escape. Sequential peeling resolves this by running three independent \( K = 1\) optimizations on shrinking atom pools. The peeling order follows global rigidity. The most rigid domain is always found first, so the order (CORE \( \to\) LID \( \to\) NMPbind) is data-determined rather than user-specified.

Pan/Zoom

A \( \sigma\) sweep at \( K = 1\) on the full 856-atom pool identifies \( \sigma_\text{op} = 3.0\)  Å as maximizing \( J_\text{CORE}\) (Figure 21). Here \( J_X\) is the Jaccard score for domain \( X\) , which measures the overlap between the atoms recovered by the BRMSD module and the experimental ground truth, so this criterion requires ground-truth domain labels. Figure 5 reports the recovered domains through the per-residue weight (the mean BRMSD weight of each residue’s backbone atoms). Each round is normalized by its own maximum, so that the most rigid residue in a round equals one. Round 1 (856 atoms, 5 restarts) recovers CORE with \( J_\text{CORE} = 0.94\) and claims 616 atoms above the \( t = 0.50\) threshold. Round 2 (240 non-CORE atoms, 5 restarts) identifies the LID domain (\( J_\text{LID} = 0.99\) ), and Round 3 (96 remaining atoms, 5 restarts) recovers NMPbind with \( J_\text{NMPbind} = 1.00\) . The Jaccard scores for all three domains are given in Figure 22. At \( t = 0.50\) the CORE mask absorbs 26 of the 120 NMPbind backbone atoms at the CORE–NMPbind boundary, but the remaining NMPbind and LID cores stay cleanly separated, so both flexible domains are recovered essentially perfectly.

A sensitivity sweep over \( t \in [0.40, 0.95]\) shows that \( J_\text{CORE}\) is flat (the step-1 optimization is independent of \( t\) ), while \( J_\text{LID}\) and \( J_\text{NMPbind}\) degrade monotonically as \( t\) increases. For \( t \leq 0.60\) , the CORE mask slightly overclaims (\( N_\text{core} = 580\) –636 vs. the ground-truth 584), but the non-CORE pool still separates cleanly into LID and NMPbind, giving \( J_\text{LID} \geq 0.94\) and \( J_\text{NMPbind} = 1.00\) . As \( t\) increases beyond 0.70, CORE is progressively underclaimed (\( N_\text{core} = 413\) at \( t = 0.80\) , 201 at \( t = 0.90\) ), leaving rigid CORE hinge residues in the pool that dominate the later steps and collapse both \( J_\text{LID}\) and \( J_\text{NMPbind}\) toward zero.

Under the \( \theta = M\sigma^2\) convention, each atom’s converged weight depends only on its own mean-square fluctuation relative to \( \sigma\) , independent of the number of frames \( M\) and of the number of atoms \( N\) in the selection. We verified both invariances on the SND3 alignment at \( \sigma_\text{op} = 0.43\)  Å (Figure 23). We quantify the difference between two normalized weight vectors \( \boldsymbol{w}\) and \( \boldsymbol{w}'\) by the Jensen–Shannon distance, the square root of the Jensen–Shannon divergence (a symmetrized form of the KL divergence).

\( D_{\rm {JS}}\) is bounded in \( [0,1]\) and is well defined even when some weights vanish. \( D_{\rm {JS}} = 0\) for identical weight vectors and \( D_{\rm {JS}} = 1\) for distributions with disjoint support.

Recomputing the weights from alternate frames (1501 of 3001 frames) leaves the per-residue weights unchanged (\( D_\text{JS} = 0.012\) ), confirming \( M\) -invariance (Figure 23A). Removing an entire mobile region from the selection and re-optimizing on the retained atoms likewise leaves their weights unchanged. Dropping TMD-1 (\( D_\text{JS} = 0.014\) ) (Figure 23B) or the cytosolic claws (\( D_\text{JS} = 0.074\) ) (Figure 23D) has almost no effect, whereas removing the rigid TM groove, which carries the dominant weight, forces a different solution, yielding a higher \( D_\text{JS} = 0.827\) (Figure 23C). The weights are thus invariant to random changes in the ensemble but responsive to targeted changes.

The weight update is a nonlinear fixed-point map, so the free energy \( G\) can in principle support several self-consistent solutions. Seeding the SND3 alignment from 40 spatially localized weight vectors (Gaussian bumps of half-width 10 Å centered on random atoms) and clustering the converged solutions by \( D_\text{JS}\) shows that, at the operating point \( \sigma_\text{op}\) , all initializations converge to a single solution (Figure 24). As \( \sigma\) decreases below \( \sigma_\text{op}\) , the landscape fragments into multiple competing minima (Figure 24A), each concentrating weight on a different local rigid sub-structure (Figure 24C), and all with higher \( G\) than the global minimum (Figure 24B). This justifies the single initialization used for the alignment and smoothing modules at \( \sigma_\text{op}\) . It also explains why clustering and domain identification use random restarts.

Calinski–Harabasz (CH) Index

Minimum Inter-centre Separation

Population Balance

Jaccard Index

All five modules are implemented in the open-source Python package BRMSD. The package uses NumPy for numerical kernels, MDAnalysis[22, 23] for trajectory I/O and atom selection, and joblib for optional multi-threaded parallelism over restarts. Each module exposes a self-contained function (BRMSD_iterate, BRMSD_focus_iterate, BRMSD_smooth, soft_kmeans, identify_domains_sequential) with a common interface: coordinate array, \( \sigma\) value, and module-specific hyperparameters (\( W\) , \( \mu/\theta\) , \( \tau\) , \( K\) , \( \mu\) , \( n_r\) ). Python 3.10+ and MDAnalysis 2.x are required.

For the alignment and smoothing modules, a single initialization from the first frame (\( s^{(1)} = x_1\) , \( w_\alpha^{(1)} = W_\alpha\) ) is used. Convergence is declared when the per-iteration changes in both the average structure and the atom weights fall below the tolerance tol (default \( 10^{-3}\) ), that is \( \Delta_\text{struct} < \texttt{tol}\) and \( \Delta_\text{weights} < \texttt{tol}\) . For clustering and domain identification, \( n_r\) restarts are run in parallel using joblib, each initialized with random center selection. The restart with the lowest converged \( G\) is retained. For the benchmarks reported here, \( n_r = 3\) restarts are used for clustering and \( n_r = 5\) for domain identification. In both cases, the selected solution was reproduced in at least two of the \( n_r\) runs, confirming reproducibility of the reported optima.

The per-iteration computational cost scales as \( \mathcal{O}(MN)\) for the weight update and \( \mathcal{O}(N^3)\) for the Kabsch rotation, where \( M\) is the number of frames and \( N\) the number of atoms in the selection. Practical runtimes are summarized in Table 1. Alignment of the SND3 trajectory (\( M = 3001\) , \( N = 3013\) all-atom) converges in 31 iterations. Smoothing adds a multiplicative factor of \( M\) over alignment but benefits from multi-process parallelism over frames via joblib. Focused alignment converges similarly to standard alignment at each \( \mu/\theta\) value. Clustering and domain detection use joblib parallelism over restarts. Systems with substantially larger frame counts (\( M \gtrsim 10^5\) ) or atom selections (\( N \gtrsim 10^4\) ) would benefit from GPU-accelerated matrix operations, which the current implementation does not support.

Movie S1. Raw (left) and BRMSD-aligned (right) SND3 trajectories are shown side by side (\( \sigma_\text{op} = 0.43\)  Å). The lighter colors represent the weight-concentrated transmembrane groove.

Movie S2. Raw (left) and smoothed (right) SND3 trajectories are shown side by side (\( W_\text{op} = 20\) frames, 4 ns full triangular window, \( \sigma_\text{op} = 0.43\)  Å). Sub-nanosecond fluctuations are suppressed while slow movements are preserved.

Pan/Zoom