Nuria Cirauqui Diaz

Elisa Frezza

Juliette Martin
email: juliette.martin@ibcp.fr

Nuria Cirauqui Diaz

Using Normal Mode Analysis on protein structural models. How far can we go on our predictions?

Keywords: Normal mode analysis, internal coordinates, coarse-grained models, protein flexibility, protein motions, conformational changes, model refinement, structural models

Using normal mode analysis on protein structural models. How far can we go on our predictions?

INTRODUCTION

Normal Mode Analysis (NMA) is a computationally inexpensive method extensively used to predict large amplitude motions in proteins, which often relate to biologically relevant functions [START_REF] Bahar | Global dynamics of proteins: bridging between structure and function[END_REF][START_REF] Kim | Robust elastic network model: A general modeling for precise understanding of protein dynamics[END_REF][START_REF] Rueda | Thorough Validation of Protein Normal Mode Analysis: A Comparative Study with Essential Dynamics[END_REF][START_REF] Yang | How well can we understand large-scale protein motions using normal modes of elastic network models?[END_REF] . In many cases, a few low-energy normal modes account for most of the structural differences between two conformational states [START_REF] Hinsen | Analysis of domain motions in large proteins[END_REF][START_REF] Krebs | Normal mode analysis of macromolecular motions in a database framework: developing mode concentration as a useful classifying statistic[END_REF][START_REF] Marques | Hinge-bending motion in citrate synthase arising from normal mode calculations[END_REF][START_REF] Tama | Conformational change of proteins arising from normal mode calculations[END_REF] . Moreover, the atomic fluctuations obtained by normal modes match experimental B-factors [START_REF] Bahar | Vibrational Dynamics of Folded Proteins: Significance of Slow and Fast Motions in Relation to Function and Stability[END_REF][START_REF] Keskin | Proteins with similar architecture exhibit similar large-scale dynamic behavior[END_REF] , as well as those calculated based on molecular dynamics simulations [START_REF] Rueda | Thorough Validation of Protein Normal Mode Analysis: A Comparative Study with Essential Dynamics[END_REF][START_REF] Micheletti | Accurate and efficient description of protein vibrational dynamics: comparing molecular dynamics and Gaussian models[END_REF][START_REF] Pontiggia | Anharmonicity and self-similarity of the free energy landscape of protein G[END_REF][START_REF] Skjaerven | Principal component and normal mode analysis of proteins; a quantitative comparison using the GroEL subunit[END_REF][START_REF] Yang | Close correspondence between the motions from principal component analysis of multiple HIV-1 protease structures and elastic network modes[END_REF] . For these reasons, NMA is widely applied to study important transition pathways of biomolecules [START_REF] Dobbins | Insights into protein flexibility: The relationship between normal modes and conformational change upon protein-protein docking[END_REF][START_REF] Zheng | Protein conformational transitions explored by mixed elastic network models[END_REF] , so as to get insights into allosteric pathways [START_REF] Fanelli | Structure network analysis to gain insights into GPCR function[END_REF][START_REF] Rodgers | Modulation of global low-frequency motions underlies allosteric regulation: demonstration in CRP/FNR family transcription factors[END_REF][START_REF] Seckler | The interplay of structure and dynamics: insights from a survey of HIV-1 reverse transcriptase crystal structures[END_REF] . It can also be applied to generate input conformations for molecular docking [START_REF] Dobbins | Insights into protein flexibility: The relationship between normal modes and conformational change upon protein-protein docking[END_REF][START_REF] Cavasotto | Normal mode-based approaches in receptor ensemble docking[END_REF][START_REF] Cavasotto | Representing receptor flexibility in ligand docking through relevant normal modes[END_REF][START_REF] Dietzen | On the applicability of elastic network normal modes in small-molecule docking[END_REF][START_REF] Jiménez-García | LightDock: a new multi-scale approach to protein-protein docking[END_REF][START_REF] Lindahl | Refinement of docked protein-ligand and protein-DNA structures using low frequency normal mode amplitude optimization[END_REF] , and even for structural refinement of X-ray diffraction data [START_REF] Kundu | Dynamics of proteins in crystals: comparison of experiment with simple models[END_REF][START_REF] Poon | Normal mode refinement of anisotropic thermal parameters for a supramolecular complex at 3.42-A crystallographic resolution[END_REF] , small angle X-ray scattering (SAXS) [START_REF] Panjkovich | Deciphering conformational transitions of proteins by small angle X-ray scattering and normal mode analysis[END_REF] or cryo-electron microscopy (cryo-EM) [START_REF] Brink | Experimental verification of conformational variation of human fatty acid synthase as predicted by normal mode analysis[END_REF] .

NMA can be performed following different approaches and algorithms [START_REF] Dykeman | Normal mode analysis and applications in biological physics[END_REF] . There are two main strategies for computing normal modes, depending on the description of the degrees of freedom: Cartesian Coordinate Space (CCS) and Internal Coordinate Space (ICS). The former is the most popular approach and is computationally simpler [START_REF] Bahar | Coarse-grained normal mode analysis in structural biology[END_REF][START_REF] Bahar | Normal mode analysis of biomolecular structures: functional mechanisms of membrane proteins[END_REF][START_REF] Stein | A Systematic Study of the Energetics Involved in Structural Changes upon Association and Connectivity in Protein Interaction Networks[END_REF][START_REF] Wako | Normal mode analysis as a method to derive protein dynamics information from the Protein Data Bank[END_REF] . In this approach, the Cartesian coordinates of all the atoms or a subset of them are used as variables. However, several works showed the advantages of computing NMA in ICS, usually using dihedral angles as variables, though any other internal property (i.e. bonds and valence angles) could be also considered [START_REF] Wako | Normal mode analysis as a method to derive protein dynamics information from the Protein Data Bank[END_REF][START_REF] Bray | Optimized torsion-angle normal modes reproduce conformational changes more accurately than cartesian modes[END_REF][START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF][START_REF] Hoffmann | NOLB: Nonlinear Rigid Block Normal-Mode Analysis Method[END_REF][START_REF] López-Blanco | iMODS: internal coordinates normal mode analysis server[END_REF][START_REF] Mendez | Torsional network model: normal modes in torsion angle space better correlate with conformation changes in proteins[END_REF] .

First, normal mode analysis in ICS allows to extend the validity of the harmonic approximation of the conformational energy hypersurface [START_REF] Kitao | Comparison of normal mode analyses on a small globular protein in dihedral angle space and Cartesian coordinate space[END_REF] . Second, this approach offers another advantage as it includes a reduced number of variables, which results in a reduced number of modes. In this way, while for NMA in CCS most of the relevant protein motions are located within the first non-trivial 10-20 modes [START_REF] Stein | A Systematic Study of the Energetics Involved in Structural Changes upon Association and Connectivity in Protein Interaction Networks[END_REF] , this number can be reduced to 5 in NMA in ICS [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] and fewer modes contribute to the conformational change [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF][START_REF] Mendez | Torsional network model: normal modes in torsion angle space better correlate with conformation changes in proteins[END_REF] . Finally, this method provides a third, intrinsic advantage: the atom connectivity is conserved, allowing to model larger conformational changes. Indeed, our previous studies showed the better performance of NMA in ICS in generating protein conformations of the bound state, when starting from the unbound structure [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] .

NMA is typically used in conjunction with coarse-grained protein models in which the pseudoatoms are connected by springs if their distance is closer than a chosen cutoff distance. This strategy allows reducing the computational cost and using the starting structure as a reference for calculations [START_REF] Wako | Normal mode analysis as a method to derive protein dynamics information from the Protein Data Bank[END_REF][START_REF] Atilgan | Anisotropy of fluctuation dynamics of proteins with an elastic network model[END_REF] . Coarse-grained models represent proteins in a broad range of particle size, from a single particle per protein [START_REF] Hills | Insights from Coarse-Grained Gō Models for Protein Folding and Dynamics[END_REF][START_REF] Peeters | Group theory of icosahedral virus capsid vibrations: a top-down approach[END_REF] , to several particles per residue [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Gopal | PRIMO/PRIMONA: A coarse-grained model for proteins and nucleic acids that preserves near-atomistic accuracy[END_REF][START_REF] Monticelli | The MARTINI Coarse-Grained Force Field: Extension to Proteins[END_REF][START_REF] Sterpone | The OPEP protein model: from single molecules, amyloid formation, crowding and hydrodynamics to DNA/RNA systems[END_REF][START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF] ; but the most common representation considers merely the alpha carbons [START_REF] Wako | Normal mode analysis as a method to derive protein dynamics information from the Protein Data Bank[END_REF] . Nonetheless, all-atom calculations are also being applied, which often allow for more accurate frequency calculation [START_REF] Dykeman | Normal mode analysis and applications in biological physics[END_REF] .

It has already been proven that both the protein representation and the variable system influence the number of low-energy modes needed to describe a relevant conformational change [START_REF] Dykeman | Normal mode analysis and applications in biological physics[END_REF][START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF] .

Several studies have demonstrated the conservation of large amplitude motions (and therefore low-energy modes) among homologous proteins [START_REF] Cirauqui | Structural, physicochemical and dynamic features conserved within the aerolysin pore-forming toxin family[END_REF][START_REF] Kolan | Elastic network normal mode dynamics reveal the GPCR activation mechanism[END_REF][START_REF] Laberge | Common dynamics of globin family proteins[END_REF][START_REF] Lai | A phylogenetic analysis of normal modes evolution in enzymes and its relationship to enzyme function[END_REF][START_REF] Leo-Macias | An analysis of core deformations in protein superfamilies[END_REF] , and even within a set of proteins that possess the same fold despite low sequence identity [START_REF] Cirauqui | Structural, physicochemical and dynamic features conserved within the aerolysin pore-forming toxin family[END_REF][START_REF] Hollup | Exploring the factors determining the dynamics of different protein folds[END_REF][START_REF] Zheng | A comparative study of motor-protein motions by using a simple elastic-network model[END_REF] . Low energy modes have thus proven to be robust to sequence variations [START_REF] Yang | Close correspondence between the motions from principal component analysis of multiple HIV-1 protease structures and elastic network modes[END_REF][START_REF] Hollup | Exploring the factors determining the dynamics of different protein folds[END_REF][START_REF] Zheng | Low-frequency normal modes that describe allosteric transitions in biological nanomachines are robust to sequence variations[END_REF] . Along with this observation comes the question of whether or not normal mode analysis could be applied on protein structural models.

Due to the lack of experimental data, it often happens that no structure is available for the protein of interest, and molecular modeling is the only available tool for predicting it. When working with molecular dynamics simulations, protein models usually deviate from the starting configuration, and therefore from the native structure, partly due to the limited accuracy of force fields [START_REF] Chen | Can molecular dynamics simulations provide high-resolution refinement of protein structure?[END_REF][START_REF] Heo | What makes it difficult to refine protein models further via molecular dynamics simulations?[END_REF] . On the other hand, normal mode calculations were presented as more robust [START_REF] Hollup | Exploring the factors determining the dynamics of different protein folds[END_REF][START_REF] Tiwari | Conservation of intrinsic dynamics in proteins -what have computational models taught us?[END_REF] , although no extensive analysis has been done in this regard. Hollup and co-worker showed, using computer-generated models of three proteins, that intrinsic dynamics are maintained when the main protein architecture is conserved [START_REF] Hollup | Exploring the factors determining the dynamics of different protein folds[END_REF] . Nonetheless, to our knowledge, no study at large scale has been conducted till date to assess the reliability of mode calculation when we are working with models instead of with the native structure.

To fill this gap, we herein thoroughly analyze the reliability of normal mode calculations performed on computational models. Due to its proven better performance, we used NMA in internal coordinate space to carry out our study, using the same protein representation as in our previous works [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . Nonetheless, to compare with standard approaches, we also computed NMA in CCS considering only Cα atoms. We performed this work on a large scale, taking advantage of the available models of the Critical Assessment of protein Structure Prediction (CASP) experiment [START_REF] Kinch | CASP 11 target classification[END_REF][START_REF] Moult | Critical assessment of methods of protein structure prediction (CASP)-Round XII[END_REF] , selecting 99 native structures and 420 models with values of Root Mean Square Deviation (RMSD) to the native structure up to 5 Å. We assessed the accuracy of normal modes by analyzing the mode overlap, the global flexibility and residue fluctuations. Moreover, we checked whether or not we could use NMA as a tool to improve protein models. Our results show that, for models with initial RMSD to the native structure smaller than 3 Å, results are reliable in terms of both motion prediction and flexibility, while the reliability decreases for less accurate models. Nonetheless, up to 5 Å, mode prediction could be still used to get a qualitative idea of the native motions, with most Root Weighted Square Inner Product (RWSIP) values native/model higher than 0.5. The prediction of flexible regions is accurate for models until 4 Å RMSD. Finally, we observed that NMA in internal coordinates was able to better refine models than NMA in CCS, while preserving structure valence.

METHODS

Data set

The two last experiments available when this work started, CASP11 and CASP12, were used as a source of structural models [START_REF] Moult | Critical assessment of methods of protein structure prediction (CASP)-Round XII[END_REF][START_REF] Moult | Critical assessment of methods of protein structure prediction: Progress and new directions in round XI[END_REF] . In CASP, proteins are divided by domains, defined as "structurally compact evolutionary modules in proteins that serve as the basic units of folding" [START_REF] Kinch | CASP 11 target classification[END_REF] . On multi-domain structures, the lowest modes usually correspond to motions between domains. We thus excluded multi-domain targets from our analysis, in order to focus on intra-domain motions. Incomplete models (i.e. not presenting all residues resolved in the native structure) were discarded. We selected only models with Cα RMSD to the native structure lower than 5 Å, because above that value comparisons will tend not to be meaningful. After those filters, we finished with a very unequal number of models per target. In order not to enrich our results by a particular protein; we selected a maximum of 5 models per target. For that, we performed clustering with the gromos method of the g_cluster program of the gromacs v5.1 package [START_REF] Berendsen | GROMACS: A message-passing parallel molecular dynamics implementation[END_REF] , so as to obtain 5 clusters at the maximum and select the center of each cluster. The final benchmark comprises 99 native structures and a total of 420 models. They are here named according to their code in the respective CASP experiment.

Classification of native structures by properties

The classification by oligomerization state was based on the biological subunit of the PDB entry, indicated on the CASP database. When this entry was not noted in CASP, we looked for it based on the protein sequence, and checked that both structures were similar, though some minimal differences in loop refinement were found. Secondary structure was assigned using DSSP [START_REF] Kabsch | Dictionary of protein secondary structure: Pattern recognition of hydrogen-bonded and geometrical features[END_REF][START_REF] Touw | A series of PDB-related databanks for everyday needs[END_REF] .

Residues with no DSSP assignment were classified loops. Proteins with more loop residues than the 3 rd quartile of the distribution were classified as loop-rich. Proteins with more than 5 times α-helix residues than β-sheet residues were classified as α -rich. Proteins with more than 5 times β-sheet residues than α-helix residues were classified as β -rich.

Internal Normal Mode Calculation

Normal mode analysis (NMA) is a technique to investigate the vibrational motion of a harmonic oscillating system around a minimum energy conformation. NMA are usually computed using Cartesian coordinate space [START_REF] Goldstein | Classical mechanics[END_REF][START_REF] Ma | Usefulness and limitations of normal mode analysis in modeling dynamics of biomolecular complexes[END_REF] . In our study, we computed NMA in internal coordinate space (iNMA) since this approach allows to extend the validity of the harmonic approximation of the conformational energy hypersurface [START_REF] Kitao | Comparison of normal mode analyses on a small globular protein in dihedral angle space and Cartesian coordinate space[END_REF] . Under a harmonic approximation in iNMA, Hessian and kinetic energy matrices can be diagonalized, allowing the analytic solution of the Lagrangian equations of motion. The eigenvectors of this matrix are the normal modes, and the eigenvalues are related to the squares of the associated frequencies. The protein main conformational changes are usually represented by a combination of some of the lowest frequency modes, but the amplitude of this motion is not defined because it depends on several conditions such as temperature.

Calculation protocol

The calculation of internal normal modes was performed using a home-made program, which has been described and evaluated in our previous studies [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . Briefly, an anisotropic elastic network model was used, with a modified Zacharias representation of the proteins [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF][START_REF] Zacharias | Protein-protein docking with a reduced protein model accounting for side-chain flexibility[END_REF] . In this coarse-grained model, N, C′, and C α are considered explicitly (rather than using a C α pseudoatom), whereas the side chains (with the exception of Gly) are represented by either one or two pseudoatoms (SC1, SC2), each representing groups of atoms within the side chain [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . As internal variables we considered merely the backbone dihedral angles φ and ψ, as they usually change significantly more than the peptide group dihedral ω. NMA in Cartesian coordinates was performed using the Bio3D program [START_REF] Grant | Bio3d: an R package for the comparative analysis of protein structures[END_REF] , with Cα coordinates as variables (cNMA). For sake of comparison, we used similar parameters as those used for internal NMA (iNMA). In other words, an anisotropic elastic network model (ENM) [START_REF] Atilgan | Anisotropy of fluctuation dynamics of proteins with an elastic network model[END_REF] was built considering a fixed force of 0.6 kcal mol -1 for atoms under a cutoff of 15 Å.

Conversion from internal coordinates to Cartesian ones

To allow a better comparison with cNMA, normal modes calculated in iNMA were converted to the Cartesian space. This conversion was performed considering all modes, using the Taylor expansion of the Cartesian coordinates to calculate atom displacements [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF] . Afterwards, the orthogonalization of the normal modes was achieved with the Gram-Schmidt process [START_REF] Arfken | Mathematical methods for physicists[END_REF] , via the gramSchmidt function of the pracma v1.9.9 R package ( https://www.rdocumentation.org/packages/pracma ).

Evaluation of similarity of Normal Mode calculations 2.3.1 Prediction of protein flexibility

The intrinsic protein flexibility (IFlex) was predicted using the following equation [START_REF] Dobbins | Insights into protein flexibility: The relationship between normal modes and conformational change upon protein-protein docking[END_REF] :

F lex I = √ 1 N ∑ N i=1 1 w i 2 (1)
where N is the number of non-trivial modes (3 n -6 for cNMA and 2 n -1 for iNMA, for a protein of n residues) and w i is the eigenvalue associated to mode i.

The weighted difference in IFlex values between a structural model and the corresponding native structure (ΔIFlex) was calculated using the following equation:

IF lex 100 Δ = ( IF lex j IF lex -IF lex i j ) * (2)
where IFlex i and IFlex j are the IFlex values of model i and native structure j , respectively.

Frequency units are expressed in 10THz.

Root Mean Square Inner Product (RMSIP) and Root Weighted Square Inner Product (RWSIP)

The RMSIP [START_REF] Fuglebakk | Measuring and comparing structural fluctuation patterns in large protein datasets[END_REF] is a measure of the similarity between two sets of normal modes or principal components of the covariance matrix, defined as:

M SIP R = √ N ∑ N i=1 ∑ N j=1 U •V ( i j ) 2 (3)
where U i and V j are the set of eigenvectors of the NMA of model ( i ) and native structure ( j ), respectively.

The RWSIP 63 is a similar measure, which also takes into account the eigenvalues. It is defined as:

W SIP R = √ v ∑ N i=1 u i j ∑ N i=1 ∑ N j=1 u v U •V i j ( i j ) 2 (4)
where u i and v i are the eigenvalues of the covariance matrix corresponding to the eigenvectors U i and V j , for model ( i ) and native structure ( j ), respectively.

Unless otherwise stated, according to common practice, we considered the first 10 non-trivial modes ( N =10) for RMSIP calculation. For iNMA, the RMSIP/RWSIP was computed after conversion to the Cartesian Coordinate Space.

Overlap

The overlap ( O i ) between a given eigenvector i and the difference vector between a model and its corresponding native structure was calculated using the inner product:

V • U O i = i (5)
where V is the difference vector between the initial structures native/model, and U i is a particular eigenvector of the model ( i ).

Root Mean Square Fluctuation (RMSF)

The mean square fluctuation of the atomic positions (RMSF) along the internal normal modes was performed as described previously, after conversion to the Cartesian Coordinate space [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF] .

For cNMA, the calculation was performed with the fluct.nma tool of the Bio3D package [START_REF] Grant | Bio3d: an R package for the comparative analysis of protein structures[END_REF] . From RMSF calculations, temperature factors (B-factors) can be obtained by the following equation and compared with experimental ones [START_REF] Kuzmanic | Determination of Ensemble-Average Pairwise Root Mean-Square Deviation from Experimental B-Factors[END_REF] :

RM SF B = 3 8π 2 2 (6)
In X-ray experiments, B-factors (B) are usually defined as a measure of the spatial fluctuations of atoms around their average position, and their motion is described as an isotropic Gaussian distribution of displacements about this average position [START_REF] Willis | Thermal Vibrations in Crystallography[END_REF] .

Generation of modified conformations by NMA

We generated modified conformations of native structures and structural models by applying a given amplitude and phase to them using normal modes. In order to study larger conformational changes, we used mode amplitudes beyond those corresponding to room temperature (termed scaled amplitudes), as already performed by others and us [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF][START_REF] Suhre | ElNemo: a normal mode web server for protein movement analysis and the generation of templates for molecular replacement[END_REF] . To note, the use of scaled amplitudes has shown to relate better with experimental B-factors [START_REF] Suhre | ElNemo: a normal mode web server for protein movement analysis and the generation of templates for molecular replacement[END_REF] .

In iNMA, we used separately the 5 first lowest frequency modes in each direction for generating sets of modified conformations [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . By increasing the amplitude, we obtained conformations differing from the initial structure with increasing RMSD values (0.5 to 5 Å) (see Figure S4). For cNMA, we used the Bio3D procedure (magnitude value of 50, step of 1) [START_REF] Grant | Bio3d: an R package for the comparative analysis of protein structures[END_REF] .

Structures with steric clashes (as defined later on) were removed. For the analysis of model refinement using NMA, we selected, from the pool of generated conformations, the one closest to the native structure.

Clash parameter

The increase of the modified amplitude may, at some point, result in the creation of clashes, or steric overlaps, on the generated conformations. The clash parameter between a pair of pseudoatoms i and j is defined as [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] :

r r d clash ij = p,i + p,j -ij -c (7)
where d ij is the distance between atoms i and j , r p,i and r p,j are the radii for the pseudoatoms i and j , respectively, and c is a cutoff for the Van der Waals overlap. For iNMA, backbone atoms C α , N and C were considered, and a value of 1.5 Å was used for c . For cNMA, as we had only C α coordinates, c was set to a value of 1 Å. We selected these values to still allow for some soft clashes but ensuring models that could be corrected by minimization. Positive values indicate the presence of strong steric overlaps, and therefore conformations presenting positive clash ij values were discarded.

Structural deformation by bond stretching

As it is well-known that cNMA deforms structures, for example by an overstretching of backbone bonds, we calculated the RMSD of the virtual Cα i -Cα i+1 bonds as [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] :

RM SD Cα -Cα i i+1 = √ N d ∑ N d i=1 (d -d ) Cα -Cα i i+1 0 Cα -Cα i i+1 2 (8)
where N d represents the number of C α i -C α i+1 virtual bonds and d and d 0 are, respectively, the length of each virtual bond in the generated conformation by NMA and in the initial structure.

Root Mean Square Deviation (RMSD) calculation

All structural fit and all RMSD calculations between two conformations were performed using the McLachlan algorithm [START_REF] Mclachlan | Rapid comparison of protein structures[END_REF] as implemented in the program ProFit v3.1 89 , considering only the backbone atoms Cα, C and N. For modified conformations generated by cNMA, only Cα atoms were used.

Evaluation of Model Quality

As an alternative to RMSD, which requires the knowledge of the native structure, model quality was assessed with QMEAN [START_REF] Benkert | Toward the estimation of the absolute quality of individual protein structure models[END_REF] , a well-established approach to estimate model quality in the absence of the native structure, based on a linear combination of six structural descriptors. Four descriptors are statistical and derived from known structures (solvent accessibility, backbone geometry, inter-atomic packing), and two descriptors evaluate the consistency of structural features with sequence-based predictions. Qmean6 score was computed by interactive access to the QMEAN server [START_REF] Benkert | QMEAN server for protein model quality estimation[END_REF] . The score ranges from 0 (low quality) to 1 (high quality).

Statistical approach

All statistical treatment was performed using the R statistical environment [START_REF] Core | R: A Language and Environment for Statistical Computing[END_REF] . The Spearman correlation coefficient (ρ) was used for data with non gaussian distribution [START_REF] Rees | Essential Statistics[END_REF] . Distributions were compared using the Wilcoxon test (two samples) [START_REF] Wilcoxon | Individual Comparisons by Ranking Methods[END_REF] and Kruskal-Wallis test (more than two samples) [START_REF] Kruskal | Use of Ranks in One-Criterion Variance Analysis[END_REF] . The area under the ROC curve (AUC) was computed with the pROC package [START_REF] Robin | pROC: an open-source package for R and S+ to analyze and compare ROC curves[END_REF] . The AUC were used for (i) quantifying the separation of IFlex distributions between two groups of proteins with different characteristics and (ii) quantifying the quality of the prediction of flexible regions using NMA-derived RMSF. In the latter case, in each native structure, residues with high experimental B-factor (value higher than the average + 1.5 standard deviation) were classified as flexible residues and NMA-derived RMSF were used to predict these flexible residues. AUC values vary between 0 (wrong prediction for all residues) and 1 (perfect prediction), where 0.5 denotes a random prediction. The rationale to use AUC here instead of Spearman correlation is that the goal is to identify peaks in flexibility profiles, whereas correlation is influenced by the precise ranking of all residues, within and outside of the peaks.

RESULTS

To build up our dataset, we selected 99 targets of the Critical Assessment of protein Structure Prediction (CASP) experiment, with its 99 native structures and a total of 420 models (see section Methods for more details) [START_REF] Moult | A large-scale experiment to assess protein structure prediction methods[END_REF] . To clarify the naming scheme, a target means a particular protein, the native structure (NS) is the experimental structure, and the model is a prediction of its tridimensional structure. CASP models have the advantage to cover a broad range of applied methodologies, from de novo to comparative modeling, from server to human generation.

Moreover, the native structure is available and the quality of the models has been evaluated using standardized measurements [START_REF] Kryshtafovych | CASP11 statistics and the prediction center evaluation system[END_REF] . Also, CASP targets tend to be selected due to their difficulty, mostly in terms of absence of near homologues, making it an interesting source for studying normal mode applications outside straightforward cases [START_REF] Kinch | CASP 11 target classification[END_REF] . Finally, they include all kinds of proteins, with different secondary structures, either in a monomeric form or in complex with another biomolecule, representing an unbiased set for very general applications. We did not consider models with initial RMSD to their native structure higher than 5 Å, as we wanted to focus in cases with no strong mistakes in the modeling. We recall that NMA was computed using internal coordinates (hereby called iNMA), which preserves protein internal geometry, and in Cartesian coordinates using a Cα model (hereby called cNMA) for comparison, as this is still the most popular approach.

Motions predicted by NMA are not strongly affected by modeling errors

According to Fuglebakk and co-workers , the Root Mean Square Inner Product (RMSIP) of the normal modes, defined in Equation ( 3) in subsection 2.3.2, so as the similar Root Weighted Square Inner Product (RWSIP), defined in Equation ( 4) in subsection 2.3.2, are some of the best known metrics to compare sets of protein motions [START_REF] Fuglebakk | Measuring and comparing structural fluctuation patterns in large protein datasets[END_REF] . Here, we used these quantities to evaluate whether or not motions calculated from structural models resembled those predicted for their native structures. It is important to stress that, in the case of iNMA, as several sets of torsions can lead to the same final conformation, and, moreover, modes are not strictly orthogonal in the internal space, it is mandatory to convert results from the internal to the Cartesian space [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF] before calculating the RMSIP/RWSIP. In Figure 1, we present the results for RWSIP defined as the sum of the overlaps between the first 10 modes, weighted by the eigenvalues (see the section Methods for more details). We used this number of modes since it is well known that they capture most of the relevant motions in both iNMA and cNMA [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . Globally, NMA is able to capture the main protein motions using structural models, in both ICS and CCS: average RWSIP is equal to 0.69 for iNMA and 0.73 for cNMA.

As shown in Figure 1, there is a clear trend of decreasing RWSIP values as RMSD increases. In other words, very accurate models, meaning those having less than 2 Å RMSD to their native structure, always presented high RWSIP values (RWSIP > 0.65 and > 0.80 for iNMA and cNMA respectively). This trend is less clear for iNMA (Figure 1A) than for cNMA (Figure 1B), suggesting that an accurate model in terms of RMSD to the native does not ensure a correct (i.e. similar to the native) calculation of protein motions in the internal space. These results indicate that iNMA can be more sensitive to local details, as already observed in our previous studies [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . Interestingly, for models of moderate quality (RMSD greater than 2.5 Å), iNMA produces fewer cases with low RWSIP values (<0.5) than cMNA. This suggests that iNMA could be more tolerant than cNMA to model inaccuracy for medium quality models.

It is well-known that for iNMA only few modes are sufficient to model protein motions with respect to cNMA, as already observed for the prediction of large conformational changes in our previous studies [START_REF] Frezza | Internal Normal Mode Analysis (iNMA) Applied to Protein Conformational Flexibility[END_REF][START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] . In fact, when only the first 5 lowest modes are considered, a more pronounced decrease in RMSIP/RWSIP values is observed for cNMA than iNMA, with a larger number of outliers (see Figure S1). Finally, our results were not influenced by the existence of gaps in the native structure and in models ( i.e. missing residue(s) inside a given chain of the structure), which occured in 25% of the proteins under investigation (see Figure S2).

RWSIP results indicate that protein motions predicted from models are similar to those predicted from the native structures. We wanted to check if this means that they are sampling a similar conformational space, that is, if we could obtain similar conformations by applying normal modes on either a native structure or a model. Indeed, for rigid proteins which at physiological conditions would not sample more than 1 Å RMSD, we cannot expect to obtain similar and meaningful conformations working with models with initial RMSD to the native higher than this value.

Then, we selected four of the most flexible proteins in the dataset (based on the IFlex value as defined in the section Methods and analyzed below) and created normal mode modified conformations for native structures and models with RMSD about 3 Å, as explained in the section Methods. Figure 2 shows these examples of conformational sampling along the first normal mode obtained by iNMA and cNMA. For these proteins, both iNMA and cNMA predict similar motions either from the native or from the model. These results encourage using protein models in the case of flexible proteins for generating conformations, for example for protein docking. This is important since the conformational space of very flexible proteins is more difficult to explore by other approaches, such as standard molecular dynamics. We would like to stress that unrealistic overstretching of the flexible protein regions was observed for cNMA (for example targets T0783-D1 and T0805-D1 in Figure 2).

Prediction of protein intrinsic flexibility is not affected by NMA approaches

Among the several methods published to estimate protein flexibility from NMA [START_REF] Dobbins | Insights into protein flexibility: The relationship between normal modes and conformational change upon protein-protein docking[END_REF][START_REF] Camps | FlexServ: an integrated tool for the analysis of protein flexibility[END_REF][START_REF] Kovacs | Predictions of protein flexibility: First-order measures[END_REF] , Dobbins and co-workers proposed a function to predict the intrinsic flexibility of a protein based on the eigenvalues (herein termed IFlex, see Equation (1) in the subsection 2.3.1), and they observed an agreement with the observed RMSD between two functional protein states (bound-unbound) [START_REF] Dobbins | Insights into protein flexibility: The relationship between normal modes and conformational change upon protein-protein docking[END_REF] .

Here, we used the same approach to compare the prediction of protein intrinsic flexibility obtained using iNMA (IFlex i ) and cNMA (IFlex c ) for native structures (NS). Figure 3 summarizes our results. Proteins are similarly ranked by predicted flexibility independently of the approach used for NMA calculation (Spearman correlation coefficient ρ=0.86 for IFlex values calculated with either iNMA or cNMA).

Some proteins have much higher values of IFlex compared to others. The most outlier target, T0865-D1, corresponds to the C-terminal coiled-coil domain of the SH3 domain-containing kinase-binding protein 1 from human (PDB code 2N64) [START_REF] Kühn | The adaptor protein CIN85 assembles intracellular signaling clusters for B cell activation[END_REF] . Therefore, its high IFlex value is expected, as coiled-coils are known for their high flexibility. Protein target T0805-D1 is a bacterial nitroreductase, likely oligomeric (see Table S1). Protein T0865-D1 is also oligomeric (see Table S1). As we performed NMA on a single protomer, and not on the oligomer, the missing interchain interactions and spatial constraints of the second protomer can easily justify the predicted high flexibility. Target T0820-D2 is an uncharacterized protein identified in marine metagenomic data, whose oligomer status is unknown.

To better explore the obtained differences in predicted flexibility, we classified our proteins based on some selected characteristics (presence of loops, missing regions, secondary structure and oligomeric state of the native structure; Table S1). No differences were observed in IFlex values related to the presence of gaps (i.e. missing regions) or secondary structure content (Table S2). However, a significant increase in IFlex values was obtained for oligomeric proteins and, when the calculations are done with iNMA, also for loop-rich proteins. To evaluate the separation between groups, we used AUC (area under the ROC curve) values, as explained in the section Methods. In all cases, the separation among groups was more marked for iNMA. The same differences among subgroups were obtained when the flexibility was predicted using the models instead of the native structures (Table S3), with greater significance in the case of iNMA.

Therefore, iNMA is able to better capture flexibility differences related to structural aspects, both from the native structure and from protein models.

Protein intrinsic flexibility can be well captured by structural models

To evaluate our capability to predict protein intrinsic flexibility using structural models, we used

IFlex to analyze the differences between intrinsic flexibility predicted for a given native structure and its models. To do so, we computed the weighted difference in IFlex values, calculated as the difference between IFlex computed for each model and their native structure divided by the value obtained for the native structure and multiplied by 100 (∆IFlex) (see Equation (2) in the subsection 2. Even more relevant is the fact that for models presenting up to 3 Å RMSD to the native structure, 90% of the points are within the interval [-10,10] of ∆IFlex with both approaches, meaning an error lower than 10% compared to the native structure IFlex (see Figure 4). This

shows that the calculation of IFlex is reliable in the case of small modelling errors. Even for models with RMSD as high as 5 Å, 90% of the data points are within the [-19,26] interval of ∆IFlex for iNMA and [-20,38] for cNMA, showing that the calculation on structural models approximates well the flexibility of the native structure. To note, for models with RMSD higher than 3 Å, we observe that the ∆IFlex values are less spread (i.e. narrower empirical 90% interval) in iNMA than cNMA, suggesting that iNMA is more tolerant to modelling errors in the case of models with medium quality.

Noteworthy, there is a good correlation between values of ΔIFlex calculated either by cNMA or iNMA (ρ = 0.77, data not shown). In addition, there is also a good correlation between IFlex predicted from native structures and IFlex predicted from models, with ρ equal to 0.78 for iNMA and 0.80 for cNMA (see Figure S3). The above results suggest that structural errors have a similar impact on the flexibility prediction in spite of the NMA approach used.

Flexible regions can be predicted from models

The atomic displacements predicted by NMA have been shown to well correlate with experimental B-factors [START_REF] Atilgan | Anisotropy of fluctuation dynamics of proteins with an elastic network model[END_REF][START_REF] Yang | Comparisons of experimental and computed protein anisotropic temperature factors[END_REF] . In this study, we also verified if this correlation is still valid for protein models. To do so, we computed the Root Mean Square Fluctuation (RMSF) along normal modes for each model, and compared it to the experimental B-factors of its respective native structure. We used AUC computation, as explained in the section Methods, to predict regions with high experimental B-factors using NMA-derived RMSF.

The results are shown in Figure 5. Panel 5A illustrates the three comparisons that have been carried out: B-factors versus RMSF of native structures (dashed red line), B-factors versus RMSF of models (dashed blue line), and RMSF of native structures versus RMSF of models (dashed green line). As shown in Figure 5A, flexible regions are predicted with high accuracy from NMA-derived RMSF computed on native structures: average AUC is equal to 0.83 for both iNMA and cNMA. The small discrepancy observed could moreover be explained by the fact that B-factors are heavily influenced by non-thermal contributions [START_REF] Hinsen | Structural flexibility in proteins: impact of the crystal environment[END_REF][START_REF] Soheilifard | Critical evaluation of simple network models of protein dynamics and their comparison with crystallographic B-factors[END_REF] and do not reflect only the structural flexibility of proteins [START_REF] Fuglebakk | Measuring and comparing structural fluctuation patterns in large protein datasets[END_REF] . When structural models are used instead of native structures, the prediction remains accurate: average AUC is equal to 0.79 for iNMA and 0.80 for cNMA. If we now compare the RMSF values between native structures and models, we observe very high prediction success: average AUC is equal to 0.92 for iNMA and 0.93 for cNMA.

As expected, the prediction of flexible regions from models (dashed blue line in Figure 5A) is impacted by the model quality, see Figure 5B. We observe a clear decrease in prediction accuracy for models with RMSD greater than 4 Å, using both NMA approaches. Panel 5C shows an example of a pair native/model, where we can observe how flexible regions are well identified using NMA, despite the medium accuracy of the model (RMSD to the native greater than 4 Å).

Our results suggest that NMA on structural models is capable of capturing the flexible regions of the native structure independently of the approach used.

iNMA allows model refinement preserving the protein valence structure

Although NMA has been mostly used for getting insights into the main protein motions, it is becoming more and more relevant in the generation of different conformations of a protein [START_REF] Zhang | Monte Carlo replica-exchange based ensemble docking of protein conformations[END_REF] , and it has also been used for improving protein models [START_REF] Rai | Modeling G protein-coupled receptors for structure-based drug discovery using low-frequency normal modes for refinement of homology models: Application to H3 antagonists[END_REF][START_REF] Warszycki | From Homology Models to a Set of Predictive Binding Pockets-a 5-HT1A Receptor Case Study[END_REF] . Here, we aimed to analyze how close we could move a model to its native structure by applying normal modes. To this aim, for each model, we generated a set of conformations by applying a subset of normal modes computed by iNMA and cNMA, using different amplitudes (see the section Methods for more details). The modified conformation presenting the minimum RMSD to the native was selected (RMSD NMA ), and compared to the initial RMSD between the model and its native structure (see Figure 6A-B). In Figure 6, to highlight the change in protein valence structure, the dots are colored based on the RMSD of C α i -C α i+1 virtual bonds between the n ew conformation and the initial model (see the section Methods for details). For very accurate models (low RMSD values), the improvement in terms of RMSD to the native is very small for both iNMA and cNMA. On the contrary, as the initial RMSD model/native increases, iNMA outperforms cNMA, with lower RMSD NMA values while preserving protein connectivity.

To quantify the best improvement in terms of native-likeness using a statistical approach, for each model we also computed the fractional RMSD, ∆RMSD%, defined as ∆RMSD% = 100(RMSD-RMSD NMA )/RMSD (Figure 6C-D). Large values of fractional RMSD mean that the model approaches the native structure by applying a specific normal mode and amplitude. By comparing the histograms in Figures 6C and6D, one can see that iNMA represents a promising approach for model refinement. In fact, the mean ∆RMSD% value is 9.4 and 5.6 for iNMA and cNMA, respectively while the median is equal to 7.0 and 3.9, respectively. Thus, iNMA is capable of refining models by reducing the initial RMSD value to the native by 10% on an average (red line in Figure 6C). It is worth noticing that these results were obtained by considering only one mode, and therefore expected to improve when two or more modes are combined [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] .

Not only does iNMA improve the predictions compared to cNMA, but the latter can dramatically overextend flexible regions [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] , by increasing Cα i -Cα i+1 virtual bonds to unrealistic values. It is consequently understandable that when structures derived from cNMA modes are used in other computational approaches or in combination with experimental data [START_REF] Gorba | Normal-Mode Flexible Fitting of High-Resolution Structure of Biological Molecules toward One-Dimensional Low-Resolution Data[END_REF][START_REF] Pérard | Structure of the full-length HCV IRES in solution[END_REF][START_REF] Schindler | Protein-protein and peptide-protein docking and refinement using ATTRACT in CAPRI[END_REF] , further structural refinement steps are necessary that can also lose most of the conformational change obtained by cNMA. On the contrary, iNMA allows to preserve the protein valence structure and it is a better candidate for integrative modelling and for improving computational models.

DISCUSSION

Internal normal mode analysis is a fast computational approach to predict large conformational changes without modifying protein connectivity. In this study, we investigated the capability of iNMA to predict protein motion, its intrinsic flexibility, atomic displacements and conformational changes using protein models instead of native structures. This step is necessary if we would like to use structural models in computational studies when its native structure is unknown. To assess this approach, we compared the results obtained with either iNMA or cNMA.

Our results show that, when working with models of medium quality (i.e. expected RMSD to the native structure lower than 5 Å), normal mode predictions do not diverge substantially with those of the native structure. The predicted motions do resemble (RWSIP above 0.5) and normal modes are able to locate the most flexible regions of the protein.

Therefore, in most cases, we could use NMA to generate modified conformations of computational models. Nonetheless, calculations are only strictly reliable for very accurate models (RMSD to the native up to 2-3 Å). Above that cutoff, the biggest errors can be obtained in the calculation of intrinsic flexibility, mostly for very flexible proteins.

As expected, we observed that iNMA is more sensitive to structural characteristics related to flexibility, ranking as more flexible those proteins whose native structure corresponds to an oligomeric state, and also those presenting very long loops. Moreover, it is shown as a better approach for creating modified conformations and for refining protein models while preserving the protein valence structure. Also, when applied for model refinement, it was able to get closer to the native structure than cNMA using only a single mode. These results are very promising and we expect that modelling refinement will be further improved by combining at least 2 modes, as already observed in our previous study on protein conformational changes [START_REF] Frezza | Internal Coordinate Normal Mode Analysis: A Strategy To Predict Protein Conformational Transitions[END_REF] .

Another important point is how to evaluate the expected accuracy of a model. There are several ways to assess it. For example, sequence identity template/target for homology modelling, or confidence scores in de novo modeling, which depend on the significance of threading template alignments and the convergence parameters. Moreover, some approaches are being developed with the aim of predicting model accuracy based on quality assessment tools, such as QMEANDisCo [START_REF] Studer | QMEANDisCo-distance constraints applied on model quality estimation[END_REF] . While those methods are promising, they are still under active development. Here, we observed no clear relationship between the models' QMEAN scores and RWSIP values obtained from the comparison of NMA computed on models and native structures (see the section Methods for the details and Figure S5). This means that a good score for the nativeness of a model does not ensure a native-like calculation of normal modes. Nonetheless, taking into account the expected RMSD of a model and the results presented here, it is possible to get an idea of the reliability of NMA performed on that model, and therefore the usage that we can make of any information coming from them.

From another perspective, our study also contributes to a better understanding of protein dynamics. Earlier studies based on NMA computed on either proteins with similar fold or protein models (Tiwari and Reuter 2018) helped to better understand how protein intrinsic dynamics is guided by its global architecture, but it is also fine-tuned for a specific function. Our results also support this picture: small modelling errors (up to 2-3 Å) do not strongly affect either protein flexibility or its dynamics, although some differences can be observed, which might be relevant in structures differing before the application of NMA by more than 3 Å.

In summary, the outcome of our study suggests that iNMA is a more suitable tool than cNMA for improving structural models, and for integrating them with either experimental data (i.e. SAXs or Cryo-EM) or other computational techniques, such as protein docking.

Nonetheless, doing so without prior knowledge of the native structure and/or its expected conformational changes is still challenging. Some advances have been made in this field, such as the reduction in the number of modes that needs to be considered (Frezza and Lavery, 2019). To further enhance iNMA as a predictive tool, we plan to conduct an analysis on a range of biological and physicochemical properties, searching for potential indicators pointing to relevant modes and movements. S1). (D) Sames as (C) but for modes calculated with cNMA.

3 . 1 )

 31 . Positive values of ΔIFlex denote models which are more flexible than the corresponding native structure, while negative values denote the opposite. The results are shown in Figure 4 as a function of the initial RMSD model/native. The median values of ΔIFlex are 0.86 for iNMA and 0.31 for cNMA, and the means are 3.70 and 1.87 respectively. Those positive values indicate a slight tendency to create models with higher flexibility than the native protein structure. Outliers with negative values (i.e. under-estimated flexibility) correspond to very flexible proteins highlighted in Figure 3.

Figure 1 .

 1 Figure 1. Motions predicted by NMA are not strongly affected by modeling errors. (A)

Figure 2 .

 2 Figure 2 . iNMA allows native-like conformational sampling of protein models with

Figure 3 .

 3 Figure 3. Prediction of protein intrinsic flexibility is not affected by NMA approaches . In

Figure 4 .

 4 Figure 4. Protein intrinsic flexibility is well captured by structural models. The weighted

Figure 5 .

 5 Figure 5 . Flexible regions can be predicted from models . (A) Overview of the comparison

Figure 6 .

 6 Figure 6. iNMA allows model refinement preserving the protein valence structure . (A)

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2 . iNMA allows native-like conformational sampling of protein models with structure

Figure 3 .

 3 Figure 3. Prediction of protein intrinsic flexibility is not affected by NMA approaches .

Figure 4 .

 4 Figure 4. Protein intrinsic flexibility is well captured by structural models.

Figure 5 .

 5 Figure 5 . Flexible regions can be predicted from models .

Figure 6 .

 6 Figure 6. iNMA allows model refinement preserving the protein valence structure .

Figure S1 .

 S1 Figure S1. Results of motion comparison are not impacted by considering either RMSIP or RWSIP, or by the

Figure S2 .

 S2 Figure S2. Gaps (i.e. missing regions) in structures have no impact on the mode reproduction . (A) For

Figure S3 .

 S3 Figure S3. There is a good correlation between IFlex predicted from native structures and IFlex predicted

Table S2 . Analysis of the predicted flexibility for different subgroups of native structures, classified based on structural features (related to Table S1

 S2 ). Results for iNMA (IFlex I ) and cNMA (IFlex C ). Bold fonts indicate significant differences at the 5% threshold (p-value of the Mann-Whitney or Wilcoxon test).

	Feature	Classification	Number of	Mean IFlex I	Mean IFlex C
			structures in	value	value
			dataset	(± standard	(± standard
				deviation)	deviation)
	Secondary	Both	56	0.97 ± 0.48	1.39 ± 0.51
	structure
		Beta	27	0.93 ± 0.15	1.35 ± 0.22
		Alpha	16	1.15 ± 0.66	1.92 ± 1.71
	p-value		0.0846300	0.1314000
	AUC		0.611	0.582
	Oligomerization	Monomer	56	0.

89 ± 0.14 1.3 ± 0.21 Oligomer 43 1.12 ± 0.65 1.67 ± 1.17

		p-value		0.0040640	0.0109300
		AUC		0.670	0.650
	Loops	No	74	0.96 ± 0.43	1.47 ± 0.90
		Yes	25	1.06 ± 0.52	1.45 ± 0.41

Table S3 . Analysis of the predicted flexibility for different subgroups of models, classified based on the structural features or their native structures (related to Tables S1 and S2).

 S3 Results for iNMA (IFlex I ) and cNMA (IFlex C ). Bold fonts indicate significant differences at the 5% threshold (p-value of the Mann-Whitney or Wilcoxon test).

	Feature	Classification	Number of	Mean IFlex I	Mean IFlex C
			models in	value	value
			dataset	(± standard	(± standard
				deviation)	deviation)
	Secondary	Both	241	0.96 ± 0.34	1.44 ± 0.62
	structure
		Beta	108	0.95 ± 0.28	1.41 ± 0.35
		Alpha	70	1.16 ± 0.60	1.95 ± 1.73
	p-value		0.000984	0.001655
	AUC		0.600	0.594
	Oligomerization	Monomer	234	0.

91 ± 0.24 1.37 ± 0.35

		Oligomer	185	1.08 ± 0.51	1.71 ± 1.25
		p-value		5e-07	0.0015520
		AUC		0.643	0.590
	Loops	No	316	2.21 ± 1.39	1.52 ± 0.98

This project/research has received funding from the European Union's Horizon 2020 Framework Programme for Research and Innovation under the Specific Grant Agreement No. 785907 (Human Brain Project SGA2). We would like to thank Guillaume Launay for useful discussions during the course of this study.

SUPPORTING INFORMATION

Figure S4. Scheme of the methodology for creation of generated conformations by iNMA. For each original structure, 10 conformations were created along the two directions of the first 5 normal modes, till an RMSD of 0.5 Å between the new conformation and the original one was reached (set shadowed in green). Then, the amplitude was increased to generate similarly 10 conformations but with an RMSD of 1 Å to the original structure (set shadowed in orange), and subsequently till an RMSD of 5 Å was obtained (not represented). The modes are indicated by the numbers 1-5, and the directions by a sign (+,-). The conformations with clashes, which were removed, are cancelled by a cross.


Didier Dubois
email: dubois@irit.fi

Henri Prade
email: prade@irit.fi

Refining aggregation operations in finite ordinal scales

Keywords: Aggregation, Qualitative scale, Decision

We discuss how to aggregate multiple criteria evaluations belonging to qualitative, linearly ordered scales. Qualitative aggregation operations such as min or max can be refined by discrimin and leximin orderings, in agreement with the Pareto ordering of vector evaluations. Further refinements of discrimin orderings, as well as the generalization of discrimin and lexirnin t o functions other than min or max, are presented. Lastly, it is pointed out that a generalized leximin is not sufficient for describing any aggregation structure. However, the definition of aggregation structures in qualitative scales amounts, in practice, to the specification of a small number of positionings of aggregation results.

Introduction

The fuzzy set framework provides a panoply of aggregation connectives for combining multiplecriteria evaluations, under a commensurability assumption: the merits of actions or items according to each criterion must be graded on the same scale, but it is usually on a numerical scale, namely [0,1]. In this paper, we investigate the possibility of not using numbers for modeling preference so as to avoid arbitrariness, improve robustness, and be more compatible with logical representations. Then the question is: can we solve multiple-criteria decision problems with a finite ordinal scale under a commensurability assumption? This might be of interest in information systems where we have to handle qualitative descriptions of preferences of users making queries.

Let be X = { a l < a2 < ... < a p ) be a finite ordinal scale. Basic aggregation operations on ordinal scales are minimum, maximum, order-statistics (i.e. the kth smallest element in a multiset), and the alpha-median med(x, y, a). Many wellbehaved aggregation operations on finite ordinal scales, such as associative smooth operations o odor, 2000), seem to be constant on significant subsets of their domains, which make these aggregations not so attractive in practice.

Besides, minimum or maximum, by retaining only the smallest or the largest evaluations, are not always sufficiently discriminating for rankordering possible choices in practice. [START_REF] Dubois | Refinem-ents of the maximin approach to decision-making in fuzzy environment[END_REF] have proposed to use two refinements of the strict min (resp. max) ordering, in agreement with the strict Pareto dominance, namely the discrimin (resp. discrimax) and the leximin (resp. leximax) ordering. Let i and y E X" be two eva-luation vectors. Let a(i, y) = {i, xi*yi) be the discriminating set of components for vectors x and -Y ---

Then define x >discrimin y

The leximin ordering is similarly defined once the two vectors have been increasingly reordered.

The extension of discrimin and leximin orderings to the comparison of vectors of different lengths has been discussed by Dubois and Fortemps (1 999).

In the following we more particularly study further refinements of discrimin on the one hand, and generalizations of discrimin and leximin t o monotone functions on the other hand, following ideas first suggested in [START_REF] Prade | Refinements of minimumbased ordering in between discrimin and leximin[END_REF] and [START_REF] Dubois | Pareto-optimality and qualita-tive aggregation structures[END_REF] respectively. Then we discuss how more general qualitative aggregation structures can be specified through the positioning of a few tuples.

Orderings between discrimin and leximin

Classical discrimin is based on the elimination of identical singletons at the same places in the comparison process of two vectors. Similarly, we can work with 2 elements subsets which are identical and pertain to the same pair of criteria.

However, simplifications can take place only one time. Thus, if the vectors are of the form ii = (u, v, u, s) and = (v, u, v, t) (with min (u, v) s min (s, t) in order to have the two vectors minequivalent), we may either delete components of ranks 1 and 2, or of ranks 2 and 3, leading in both cases to compare (u, s) and (v, t), and to consider the first vector as smaller as soon as u < min (v, s, t).

We can now introduce the definition of the (order) 2-discrimin. Let us build a set a 2 ( % $ ) as any maximal set {(i, j), disjoint pairs such that ui = vj and uj = vi) U {k, uk = vk). Then the 2discrimin is just the minimum-based ordering once components corresponding to singletons and pairs in a2(Y,$) are deleted. Note that a 2 ( k $ ) is not always unique as shown by the above example. However this does not affect the result of the comparison of the vectors after the deletion of the components in B 2 ( & ,$) as it can be checked from the above formal example, since the minimum aggregation is not sensitive to the place of the terms. The 2-discrimin also includes the deletion of identical components as in the ordinary discrimin, since it would be strange to delete pairs of identical values in the comparison but not single identical values (which may blur the comparison).

The 2-discrimin refines the discrimin ordering. It deletes pairs of values which play a neutral role in the aggregation, but which may lead to ties if these values are not ignored.

Clearly this idea can be extended to 3 elementssets as well or more generally with k elementssets. However the definition of the 3-discrimin should be further refined by choosing the permutation, which lead to a discriminant situation if possible, as shown by the following example. Let

There exist two overlapping permutations : where a < b c c c d < e c f < g is assumed.

Ignoring the values involved in the first permutation, we obtain v cmin u since b c c, while using the second permutation, we get v u since both vectors of remaining values lead to the same minimum a. Note that the discriminant permutation should involve the minimum value of the vectors, to be of interest. However they cannot exist two overlapping permutations which, after ignoring their components in the min-based comparison, would lead to opposite orderings (namely u > v and v > u). This holds for the k-dicrimin as well. This can be seen by considering vectors of the form where y and z are such that u =,inv. Simplifying by the first four components leads to compare min(x, t, y) with min(s, x, z), while deleting (s, x, t) with (t, s, x) leads to the comparison of min(t, a, b, y) with min(b, s, a, z). Clearly the two comparisons cannot disagree with each other (i.e. cannot lead to u > v and v > u respectively). This remark could be further generalized by replacing some of the above vector components by subvectors. In case of several overlaps, this analysis can be iterated.

So the generalized procedure for applying kdiscrimin is to look for permutations of orders 1, 2, . . ., k and to explore the different possibilities in case of overlapping permutations until a discriminating one is found. Note that 'kdiscrimin' amounts to a limited leximin on k-long subsequences. Thus 'k-discrimin' provides orderings in between discrimin and leximin. However the n-discrimin ordering may remain less discriminating than the leximin ordering.

The 'k-discrimin', defined as such, is not fully compatible with transitivity, since it can be anyway. This is particular for the k-discrimin (k 2 2), since the discrimin ordering >discrimin is indeed transitive.

Generalized discri-f and lexi-f

The discrimin and leximin ideas can be applied to functions other than min (or max). Given a finite, totally ordered set (X, z) with p elements top 1 and bottom 0, consider an aggregation function f : xn -, X, which, by definition, is increasing in the wide sense and such that f(1, 1

... 1) = 1 and f(0, ... 0) = 0.

Clearly the discriminating power of qualitative aggregation operations is bound to be very weak since f classifies xn into p ranked classes, and thus no aggregation function can be strictly increasing from xn to X. A consequence is in general the existence for any aggregation function f of a maximizing set S such that some f-optimal solutions in S are not Pareto optimal. One way out of this difficulty may be to use functions from xn to a bigger finite scale Y. However, this idea is not satisfactory from a practical point of view since the combinatorics of functions from xn to Y become rapidly prohibitive as Y is bigger, and are thus much higher than those of functions from xn to X.

When f = min, the natural way to tackle the problem has been to introduce relations that naturally refine the min-ordering, and restore the Pareto optimality of the selected maximal solutions. Such relations are the discrimin ordering and the leximin ordering. The discrimination power of the latter is maximal, i.e. it is equal to that of the most discriminating symmetric aggregation operations.

We try here to generalize this refinement technique to more general families of aggregation operation. We restrict to the case of symmetric functions. Consider a family { P ) of symmetric functions xP -+ X. For any positive integer p, P is supposed to be -) extensively pr$ferentially consistent with F1:

-) globally strictly monotone (assuming p 2 n): if xi > yi V i =l,p, then P(F) > P ( y ) where x and y E xn.

These conditions look natural in the scope of applications. The first condition is a weak form of preferential independence. They are satisfied by the minimum, the maximum (but not other orderstatistics). By convention f' is the identity function on X. Call {# ... f"...} a qualitative aggregation structure. The notation f is used when the number of arguments is not emphasized. The generalization of discrimin and leximin to such aggregation structures is as follows:

where M(x) denotes the multi-set induced by vector 2 (the same grade may appear several times in x), and where it can be easily checked that ~( x )

-~( y ) and ~( i ) -~( k ) have the same number of elements. It can be shown that under mild conditions such as global monotonicity and extensive preferential consistency, lexi-f and discri-f maximal solutions are indeed Pareto-optimal, and that the corresponding ordering of solutions is quite discriminant.

Specifying qualitative aggregation structures

Globally strictly monotone aggregation functions on X are easily proved to be idempotent on finite scales, since if 0 = xl < X* < .. .< X, = 1 it follows that f(xi ... xi) = xi. It rules out the Archimedean t-norms and conorm-like operations on finite sets [START_REF] Fodor | Smooth associative operations on finite ordinal scales[END_REF]. Moreover, the only associative idempotent aggregation operations different from min and max, the a-medians (f(x, y) = median(a, x, y)), are generally not globally strictly monotone, since they are constant on large subsets of x". The above extensions of leximin and discrimin orderings to these operations thus do not possess enough discrimination power; they can only be refined by directly adopting Paretoordering on the ranges where these aggregation functions are constant. The simplest non trivial example of finite totally ordered scale is X = { a , < a2 < a3}, that we shall write it as {1,2,3} for short. Adopting the lexi-f2 ordering for X does not leave many degrees o f freedom: one must indicate the relative position of f2(1, 3) and ? [START_REF] Dubois | Refinem-ents of the maximin approach to decision-making in fuzzy environment[END_REF][START_REF] Dubois | Refinem-ents of the maximin approach to decision-making in fuzzy environment[END_REF]. If ?(I, 3) > f2(2, 2), this is the leximax ordering. If ?(I, 3) < ?(2, 2), this is the lexirnin ordering. If ?(I, 3) = ?(2, 2), this is a kind of ordinal average (which is less discriminant). With three arguments, the ordering of 3-tuples is fixed by further positioning (2, 2, 2) with respect to (1, 1, 3) and (1, 3, 3) (note that (1, 3, 3) >lexi-f (1, 1, 3) in any case). The lexi-f positions of (1, 2, 2) with respect to (1, 1, 3), (1, 2, 3) with respect to (2, 2, 2), and (2, 2, 3) with respect to (1, 3, 3), are enforced by the position of (2, 2) with respect to (1, 3). If ?(I, 3) > ?(2, 2), then ?(I, 3, 3) 2 ?(2, 2, 2) and only the relative position of (2, 2, 2) and (1, 1, 3) is left open. 1f ?(2, 2) > ?(I, 3), then ?(2, 2, 2) 2 ?(I, 1, 3) and only the relative position of (2, 2, 2) and (1, [START_REF] Dubois | Computing im-proved optimal solutions to max-min flexible cons-traint satisfaction problems[END_REF][START_REF] Dubois | Computing im-proved optimal solutions to max-min flexible cons-traint satisfaction problems[END_REF] is left open.

The generation of complete preorderings o f tuples of elements from a finite ordered scale in agreement with Pareto-dominance and symmetry has been considered in Moura- [START_REF] Pires | Specifying fuzzy constraints interactions without using aggre-gation operators[END_REF] in the scope of fuzzy constraint satisfaction problems. A natural question is whether any such complete preorderings of xn can be obtained as a lexi-f ordering for some 1 qualitative aggregation structure {f . . . fn), and more generally, can be generated by a small number of extra constraints on the relative positioning of a few tuples.

Unfortunately the answer for lexi-f ordering is negative. A counterexample is obtained using a 4-element scale X = (0, 1, 2, 3 ) . Then ? is characterized by the relative positionings of (0, 2) w.r.t. (1, l), (1, 3) w.r.t. (2, 2), and (0, 3) w.r.t. (1, 2 1) and (2, 2). However, using f : X -, X, we have P(0, 3) E {?(o, 01, ?(I, 11, ?(2, 21, ? [START_REF] Dubois | Computing im-proved optimal solutions to max-min flexible cons-traint satisfaction problems[END_REF][START_REF] Dubois | Computing im-proved optimal solutions to max-min flexible cons-traint satisfaction problems[END_REF]). Neither the discri-f nor the lexi-f extension, nor even supplemented by Pareto-dominance itself, can generate the complete preorders such that (1, 1) < (0, 3) < (2, 2). There are 12 total orderings which are Pareto-compatible and respect symmetry in this example, and only 8 of them can be generated as a lexi-f ordering via an aggregation structure xn -, X. Generating the 2 other total orderings requires a function X -, Y where Y has 5 levels.

More generally the problem is how to efficiently describe Pareto-compatible rankings using a small number of conditions on the relative positioning of a few tuples. This includes the open question: how to characterize the descriptive power of lexi-f. An example of property which may simplify the study of qualitative aggregation structures is the following regularity : let (i, j) E x2 = (0, 1,2, ..., n), if ?(i, j) > ?(i+l, j-1), then ?(i+l, j+l) > ?(i+2, j>. The combinatorics of such regular aggregations functions look moderate and deserve further exploration. Another issue is the expansion of a

Concluding remarks

Numerical scales such as [0, 11 are not always meaningful in practice when users are only able to discriminate between a small number of satisfaction levels. Qualitative aggregation of fuzzy sets with membership grades belonging to such qualitative scales have been discussed. This can be relevant for modelling users' preferences and comparing solutions of constraint satisfaction problems.

The above study has exhibited some intrinsic limitations of the otherwise appealing finite setting for criteria aggregation using a single finite scale, whereby concise representations and functions having good algebraic properties turn out to lack expressivity, even under natural lexicographic-like extensions.

(

  a, c, d) with (d, e, c) and (d, a, c) with (c, d, e).

  established that x >k.discrimin y and y >k-discrimin -z only entail x zk-discrimin 2. This is due to the possible overlapping of k elements-sets which are deleted in the comparison, as shown by the following example x = (a c e), $ = (c a d), and z -----(b d a) where >2-discrimin Y, Y >2-discrimin > ---and neither x >2-discrimin z nor z >2-discrimin x with a c b c c c d c e. However 2 >,,,hin z