## 2 Overview of the Workflow  
In our approach, the identification of a common set of conformational motions among different proteins is based on the analysis of their dynamics in a CG representation; from here, a representative set of normal modes is identified through a dynamics-based clustering of the proteins comprising the initial dataset. The selected, representative modes are then orthonormalized and ordered, so as to obtain the final basis set. An overview of the workflow is given in Figure 1 and explained in detail in the following paragraphs. The starting point is the identification of a set of proteins (Figure 1a). The choice of this dataset is arbitrary and independent of the pipeline; however, the number of proteins that the dataset contains is supposed to be large enough so as to be representative of the families or superfamilies that are included, meaning that the more distant are the members in terms of homology, the larger should be the dataset. This is necessary to ensure the sufficient generality of the resulting basis set of conformational motions. The selected set of structures is used to run pairwise dynamics-based protein alignments with the ALADYN software developed by some of us [51] (Figure 1b). ALADYN takes two input structures and performs the maximization of a score function that takes into account the spatial superposition of protein regions that have similar motion. The dynamical information is encoded in the low-energy (large-amplitude) eigenvectors obtained from the diagonalization of the interaction matrix \(M_{ij}\) of the Hamiltonian function of the \(\beta\)-Gaussian network model:  
\[H=\frac{1}{2}\sum_{ij}\delta\vec{x}_{i}M_{ij}\delta\vec{x}_{j} \tag{1}\]  
where \(\delta\vec{x}_{i}\) is the displacement vector of the i-th bead with respect to the equilibrium configuration. Once the eigenvectors have been obtained, the extent and consistency of the alignment are quantified through the root-mean-squared inner product (RMSIP) between the spaces given by the first 10 modes of each aligned protein. If we call \(N_{i}\) and \(N_{j}\) the total number of residues in the chains of the two aligned proteins, the RMSIP calculation is limited to a subset \(q<N_{i},N_{j}\) of marked C\({}_{\alpha}\). These subsets of amino acids are chosen by firstly grouping the amino acids into groups of 10 subsequent ones, then maximizing a single scoring parameter via the standard Metropolis criterion over the space of possible pairs of groups among the two proteins' sequences, as exhaustively explained in [51]. Specifically, the RMSIP is defined as:  
\[\text{RMSIP}(\{\{\vec{v}_{i}^{k}\}_{i},\{\vec{w}_{m}^{k}\}_{j}\})=\text{RMSIP}_{ ij}:=\sqrt{\frac{1}{10}\sum_{l,m=1}^{10}\left|\sum_{k=1}^{q}\vec{v}_{l}^{k}\cdot \vec{w}_{m}^{k}\right|^{2}} \tag{2}\]  
The \(\text{RMSIP}\in[0,1]\) takes on the value of 1 in the case of the perfect correspondence of the spaces and 0 in the case of their complete orthogonality. The quantity (1.0-RMSIP), which still takes values in the interval \([0,1]\), is therefore suitable to define a distance in dynamics between two proteins after alignment. The statistical significance of the alignment, quantified by means of a \(z\)-score, is taken into account by weighting the RMSIP by the hyperbolic tangent of the module of the \(z\)-score, so as to give more importance to the most reliable results. The distance in dynamics between two aligned proteins \(i\) and \(j\) is therefore defined as:  
\[d_{ij}=1.0-\left(\text{RMSIP}_{ij}\cdot\text{tanh}|z_{ij}|\right) \tag{3}\]  
After all the pairwise alignments between the elements of the dataset are performed, a distance matrix that expresses differences in the large-scale dynamics is obtained (Figure 1c); then, the dataset undergoes hierarchical clustering [52] based on this distance matrix, in  
Figure 1: Schematic representation of the workflow proposed.