S E Hoff

F E Thomasen

K Lindorff-Larsen

M Bonomi
email: mbonomi@pasteur.fr

Accurate model and ensemble refinement using cryo-electron microscopy maps and Bayesian inference

Converting cryo-electron microscopy (cryo-EM) data into high-quality structural models is a challenging problem of outstanding importance. Current refinement methods often generate unbalanced models in which physico-chemical quality is sacrificed for excellent fit to the data. Furthermore, these techniques struggle to represent the conformational heterogeneity averaged out in low-resolution regions of density maps. Here we introduce EMMIVox, a Bayesian inference approach to determine single-structure models as well as structural ensembles from cryo-EM maps. EMMIVox automatically balances experimental information with accurate physico-chemical models of the system and the surrounding environment, including waters, lipids, and ions. Explicit treatment of data correlation and noise as well as inference of accurate B-factors enable determination of structural models and ensembles with both excellent fit to the data and high stereochemical quality, thus outperforming state-of-the-art refinement techniques. EMMIVox represents a flexible approach to determine high-quality structural models that will contribute to advancing our understanding of the molecular mechanisms underlying biological functions.

Introduction

Cryo-electron microscopy (cryo-EM) has become a powerful tool to determine the structures of complex biological systems. Recent advances in physical instrumentation and image processing algorithms have pushed the boundaries of what can be achieved with this technique in terms of resolution and coverage across a wide spectrum of system size and complexity [START_REF] Nogales | The development of cryo-EM into a mainstream structural biology technique[END_REF] . While most single-particle cryo-EM density maps still have a resolution between 3 Å and 4 Å (Fig. S1A), the number of systems determined at atomic resolution is steadily increasing, with the current record set by Apoferritin, resolved in 2020 at ~1.2 Å 2, 3 .

At the same time, tremendous progress in cryo-electron tomography (cryo-ET) is enabling the determination of the structure of complex biological systems in situ at sub-nanometer resolution [START_REF] Turk | The promise and the challenges of cryo-electron tomography[END_REF] (Fig. S1B). It is also becoming increasingly clear that single-structure models are not always a faithful representation of the three-dimensional (3D) cryo-EM density maps [START_REF] Bonomi | Determination of protein structural ensembles using cryoelectron microscopy[END_REF] . Small-scale continuous dynamics of functionally important flexible regions within biomolecules are typically averaged out during the reconstruction of 3D maps resulting in fuzzy density regions with reduced resolution. To accurately interpret these regions, we should therefore move away from single-structure models towards ensembles of conformations [START_REF] Tang | Conformational heterogeneity and probability distributions from single-particle cryo-electron microscopy[END_REF][START_REF] Doerr | A dynamic direction for cryo-EM[END_REF] . In all these situations, it is of paramount importance to convert cryo-EM data into highquality structural models, for example for in silico structure-based drug design [START_REF] Scapin | Cryo-EM for Small Molecules Discovery, Design, Understanding, and Application[END_REF] or training deep learning approaches on accurate structural data [START_REF] Jumper | Highly accurate protein structure prediction with AlphaFold[END_REF] .

Over the years, a variety of different metrics have been proposed to evaluate the quality of structural models obtained from cryo-EM maps [START_REF] Lawson | Cryo-EM model validation recommendations based on outcomes of the 2019 EMDataResource challenge[END_REF][START_REF] Farabella | TEMPy: a Python library for assessment of three-dimensional electron microscopy density fits[END_REF][START_REF] Pintilie | Measurement of atom resolvability in cryo-EM maps with Q-scores[END_REF][START_REF] Terashi | Residue-wise local quality estimation for protein models from cryo-EM maps[END_REF] . Generally speaking, these metrics fall into two categories: how well a model (or set of models) explains the observed map or directly the single-particle images (fit to the data) and how good the model is in terms of basic stereochemical parameters, such as length of chemical bonds, backbone and sidechain dihedral distributions, and clashes between close atoms. To evaluate the fit to the data, the most common approach is to compare the experimental map with the map predicted from a model, for example by calculating the cross correlation between the two maps over the entire 3D space or in the proximity of the structural model. In low-resolution areas, the density can be predicted either as an average over multiple conformations or from a single-structure model by introducing temperature factors (Bfactors), which implicitly account for intrinsic dynamics and other reasons behind the observed fuzzy density, such as errors in the image alignment or local damage caused by the electron beam. The main challenge in determining a high-quality structural model is to obtain a good balance between the fit to the observed density map and the overall stereochemical quality of the model. Several modelling approaches have been developed to refine structural models into cryo-EM maps [START_REF] Malhotra | Modelling structures in cryo-EM maps[END_REF] . These techniques rely on various approaches including homology modelling [START_REF] Sali | Comparative Protein Modeling by Satisfaction of Spatial Restraints[END_REF] , rigid-body fitting [START_REF] Wang | Automated structure refinement of macromolecular assemblies from cryo-EM maps using Rosetta[END_REF][START_REF] Pettersen | UCSF chimera -A visualization system for exploratory research and analysis[END_REF][START_REF] Joseph | Improved metrics for comparing structures of macromolecular assemblies determined by 3D electron-microscopy[END_REF][START_REF] Kawabata | Multiple Subunit Fitting into a Low-Resolution Density Map of a Macromolecular Complex Using a Gaussian Mixture Model[END_REF] , flexible refinement [START_REF] Igaev | Automated cryo-EM structure refinement using correlation-driven molecular dynamics[END_REF][START_REF] Trabuco | Flexible fitting of atomic structures into electron microscopy maps using molecular dynamics[END_REF][START_REF] Cragnolini | Cryo-EM structure and B-factor refinement with ensemble representation[END_REF][START_REF] Afonine | Real-space refinement in PHENIX for cryo-EM and crystallography[END_REF][START_REF] Blau | Gentle and fast all-atom model refinement to cryo-EM densities via a maximum likelihood approach[END_REF][START_REF] Topf | Protein structure fitting and refinement guided by cryo-EM density[END_REF][START_REF] Croll | ISOLDE: a physically realistic environment for model building into low-resolution electron-density maps[END_REF][START_REF] Emsley | Features and development of Coot[END_REF][START_REF] Scheres | Implementation of a Bayesian approach to cryo-EM structure determination[END_REF] , integrative modelling [START_REF] Bonomi | Bayesian Weighing of Electron Cryo-Microscopy Data for Integrative Structural Modeling[END_REF][START_REF] Van Zundert | Integrative Modeling of Biomolecular Complexes: HADDOCKing with Cryo-Electron Microscopy Data[END_REF] , and machine learning to automate model fitting to 3D density maps [START_REF] Giri | Deep learning for reconstructing protein structures from cryo-EM density maps: Recent advances and future directions[END_REF] . Most of these methods do not optimize B-factors, and therefore, while they are extremely useful to create structural models that occupy the space defined by the cryo-EM density, a quantitative evaluation of the model fit to the experimental map with the metrics commonly used for PDB validation is challenging.

On the other hand, a handful of methods have been proposed to model structural ensembles from either 2D single-particle images [START_REF] Cossio | Bayesian analysis of individual electron microscopy images: Towards structures of dynamic and heterogeneous biomolecular assemblies[END_REF][START_REF] Velazquez-Muriel | Assembly of macromolecular complexes by satisfaction of spatial restraints from electron microscopy images[END_REF][START_REF] Tang | Ensemble Reweighting Using Cryo-EM Particle Images[END_REF] or 3D maps [START_REF] Cragnolini | Cryo-EM structure and B-factor refinement with ensemble representation[END_REF][START_REF] Bonomi | Simultaneous Determination of Protein Structure and Dynamics Using Cryo-Electron Microscopy[END_REF][START_REF] Riley | qFit 3: Protein and ligand multiconformer modeling for X-ray crystallographic and single-particle cryo-EM density maps[END_REF] . As of today, one of the most popular refinement software commonly used prior to depositing models in the PDB database is PHENIX [START_REF] Afonine | Real-space refinement in PHENIX for cryo-EM and crystallography[END_REF] , which enables real-space refinement, optimization of B-factors, and modelling residues in alternative conformations. The approach implemented in PHENIX relies on an empirical scoring function to maximize the correlation between the cryo-EM map predicted from a model and the experimental map while trying to preserve stereochemical properties. This approach typically leads to an excellent fit to the data (Fig. S2AB) but often at the expense of physico-chemical properties, especially in terms of clashes between atoms (Fig. S2CD).

Here we present EMMIVox, a computational approach to determine single-structure models as well as conformational ensembles using cryo-EM maps. EMMIVox is based on a Bayesian inference framework [START_REF] Rieping | Inferential structure determination[END_REF] to balance automatically the experimental information with state-of-the-art physico-chemical models of the system and the surrounding environment. Explicit treatment of data correlation and uncertainty as well as accurate inference of B-factors contribute to the determination of structural models with excellent fit to the data without sacrificing the stereochemical quality of the models. We benchmarked our approach on nine complex biological systems and demonstrated that EMMIVox models outperformed those obtained with state-of-the-art refinement techniques and deposited in the PDB database in terms of several quality metrics.

We also illustrate how EMMIVox can be used in combination with medium-low resolution cryo-EM maps to refine coarse-grained models of large protein complexes and to determine conformational ensembles describing the structural heterogeneity hidden in low-resolution areas of atomistic cryo-EM maps.

EMMIVox is implemented in the open-source, freely available PLUMED library [START_REF] Tribello | PLUMED 2: New feathers for an old bird[END_REF][START_REF] Bonomi | Promoting transparency and reproducibility in enhanced molecular simulations[END_REF] (www.plumed.org) and aims at setting a new standard for single-structure and ensemble refinement by optimally integrating cryo-EM maps with accurate atomistic and coarse-grained physico-chemical models of the system as well as other experimental data, when available.

Results

This section is organized as follows. We first provide a general overview of the EMMIVox approach and illustrate its accuracy in refining both atomistic and coarse-grained single-structure models using cryo-EM maps at various resolutions. We then focus on ensemble refinement and present the results of our benchmark as well as one case study: the GPT type 1a tau filament from progressive supranuclear palsy neurodegenerative disease [START_REF] Shi | Structure-based classification of tauopathies[END_REF] .

Overview of the EMMIVox approach

EMMIVox is based on a Bayesian inference framework to generate a hybrid energy function that combines the molecular mechanics force fields used in classical Molecular Dynamics (MD) simulations with spatial restraints to enforce the agreement of a structural model with the density observed in the voxels of a cryo-EM map. To balance automatically stereochemical quality of the models and fit to the data, EMMIVox: i) pre-filters the voxels of a cryo-EM map to reduce the correlation between experimental data points and therefore help avoid data overweighting; ii) uses prior models of uncertainty in the experimental data obtained from independent 3D reconstructions (half maps) while allowing for the presence of random and systematic errors in the map; and iii) builds spatial restraints weighted by the estimated accuracy in each voxel. Furthermore, EMMIVox exploits modern atomistic and coarse-grained force fields to accurately describe the physico-chemical properties of a biological system as well as its environment, including water molecules, ions, lipids, and small molecules. A Monte Carlo optimization of residue-level B-factors coupled with structural refinement guided by the EMMIVox hybrid energy function enables the determination of single-structure models that fit the observed cryo-EM map while preserving a high stereochemical quality.

Instead of implicitly modelling local dynamics with B-factors, EMMIVox can be combined with metainference [START_REF] Bonomi | Metainference: A Bayesian inference method for heterogeneous systems[END_REF] to obtain structural ensembles that explicitly represent the continuous dynamics of highly flexible regions of the system. A detailed description of EMMIVox is provided in Materials and Methods.

Accurate atomistic single-structure refinement

We first benchmarked the accuracy of EMMIVox in refining single-structure models using as test systems the GPT type 1a tau filament (1.90 Å, PDB 7p6a) [START_REF] Shi | Structure-based classification of tauopathies[END_REF] , the ChRmine channelrhodopsin (2.02 Å, PDB 7w9w) [START_REF] Kishi | Structural basis for channel conduction in the pump-like channelrhodopsin ChRmine[END_REF] , the Anaplastic lymphoma kinase extracellular domain fragment in complex with an activating ligand (2.27 Å, PDB 7n00) [START_REF] Reshetnyak | Mechanism for the activation of the anaplastic lymphoma kinase receptor[END_REF] , a dimeric unphosporylated Pediculus humanus protein kinase (2.35 Å, PDB 7t4n) [START_REF] Gan | Activation mechanism of PINK1[END_REF] , the human CDK-activating kinase bound to the inhibitor ICEC0942 (2.50 Å, PDB 7b5o) [START_REF] Greber | 2.5 A-resolution structure of human CDK-activating kinase bound to the clinical inhibitor ICEC0942[END_REF] , the major facilitator superfamily domain containing 2A in complex with LPC-18:3 (3.03 Å, PDB 7mjs) [START_REF] Cater | Structural basis of omega-3 fatty acid transport across the blood-brain barrier[END_REF] , the Escherichia coli PBP1b (3.28 Å, PDB 7lq6) [START_REF] Caveney | CryoEM structure of the antibacterial target PBP1b at 3.3 A resolution[END_REF] , the mammalian peptide transporter PepT2 (3.50 Å, PDB 7nqk) [START_REF] Parker | Cryo-EM structure of PepT2 reveals structural basis for proton-coupled peptide and prodrug transport in mammals[END_REF] , and the SPP1 bacteriophage tail tube (4.00 Å, PDB 6yeg) [START_REF] Zinke | Architecture of the flexible tail tube of bacteriophage SPP1[END_REF] . These examples include soluble and membrane proteins, protein complexes, as well as systems containing ordered waters, small molecules, and lipids.

The quality of the EMMIVox models was evaluated and compared to the deposited PDBs using four key metrics. The physico-chemical quality was quantified by the clashscore [START_REF] Chen | MolProbity: all-atom structure validation for macromolecular crystallography[END_REF] and MolProbity score [START_REF] Chen | MolProbity: all-atom structure validation for macromolecular crystallography[END_REF] , while the fit to the data using the correlation coefficient between predicted and observed maps in the proximity of the model (𝐶𝐶 !"#$ ) [START_REF] Afonine | New tools for the analysis and validation of cryo-EM maps and atomic models[END_REF] and the EMRinger score [START_REF] Barad | EMRinger: side chain-directed model and map validation for 3D cryo-electron microscopy[END_REF] . EMMIVox models generally showed improved quality in all four metrics (Fig. 1). The improvements are particularly pronounced for clashscore that indicates the absence of serious clashes between atoms in the EMMIVox models, and the MolProbity score, suggesting a superior overall stereochemical quality. At the same time, EMMIVox models often improved the 𝐶𝐶 !"#$ indicating a better fit to the experimental cryo-EM map, and the EMRinger scores, particularly in maps with resolution better than 4 Å. Overall, the four metrics indicate that EMMIVox generates single-structure models that improve upon current models deposited in the PDB.

In the following sections, we will focus on the different EMMIVox components that contribute to the observed high quality of the models. Three key aspects affect the balance between stereochemical quality and fit to the data: i) the number of experimental data points (voxels) that are fit using spatial restraints; ii) the strength of these spatial restraints; iii) the accuracy of the cryo-EM map predictor from a model.

Removal of correlated data leads to balanced model refinement.

Neighboring voxels of a cryo-EM map contain correlated information, with the strength of correlation depending on the voxel size and the map resolution. Ignoring such correlation leads to overcounting the number of (independent) data points available and ultimately biasing the refinement towards (over)fitting the data at the expense of the stereochemical quality of the models. To address this point, we developed a pre-filtering procedure to subsample the set of cryo-EM voxels and reduce data correlation (Materials and Methods). We first tested our procedure using over 2400 cryo-EM density maps deposited in the EMDB with resolution ranging from 1.78 Å to 4 Å. Upon removal of correlated data, the median number of voxels per model atom decreased as the resolution of the cryo-EM map worsened and became independent of the voxel size (Fig. S3). This indicates, as expected, that the amount of spatial information provided by high-resolution maps is greater than in medium-low resolution maps.

We then benchmarked the quality of the EMMIVox single-structure refinement as a function of the amount of data removed. Removing a large portion of correlated data led to better stereochemical models compared to utilizing all available voxels (Fig. S4 andS5) at the expense of the quality of the fit to the entire cryo-EM map (Fig. S6 andS7). As fewer points were removed, the balance shifted towards the fit to the data.

We therefore identified a correlation coefficient equal to 0.8 as the optimal threshold for data removal in EMMIVox refinement to guarantee a good balance between stereochemical quality of the model and fit to the data. In our benchmark set, this threshold corresponded on average to removing 45% (± 30%) of the voxels in the proximity of the deposited PDB, resulting in 13 ± 5 voxels per atom.

Bayesian noise models reduce data overfitting. In EMMIVox, the strength of the spatial restraints used to enforce the model agreement with the cryo-EM map reflects the estimated accuracy (noise level) of the density in each voxel. Our Bayesian inference framework enables us to estimate the noise level on the fly based on the consistency between cryo-EM data, physico-chemical prior and potentially additional experimental data (Materials and Methods), and ultimately to downweigh voxels that are considered outlier data points during refinement. To guide noise inference, we developed priors based on the density variability in each voxel calculated from two independent 3D reconstructions, or half maps. In regions of the map where large variations are observed, our approach helps to avoid overfitting the data.

To exemplify this point, we examined the relation between per-residue fit to the map (local 𝐶𝐶 !"#$ ) and median noise level of the voxels around a given residue. This analysis illustrates that in the EMMIVox models residues in low-noise regions sometimes fit the map even better than in the deposited PDB, while in regions with high noise level the local fit to the map is often much worse (Fig. 2A andS8). In these regions, during refinement EMMIVox correctly reduced the weight of the experimental data in favor of the molecular mechanics forcefield, ultimately increasing the overall stereochemical quality of the model. To determine which specific physico-chemical properties of the system were better represented in our models, we computed the total number of atom clashes, hydrogen bonds, and salt bridges in the EMMIVox models and deposited PDBs across our entire benchmark (Fig. 2B). EMMIVox dramatically reduced the number of serious clashes observed in the deposited PDBs and optimized the positions of backbone and sidechain atoms so as to better describe hydrogen bond and salt bridge geometries (Fig. 2C). The improvement in these three physico-chemical fingerprints is more accentuated for residues in regions of the cryo-EM map with high noise, supporting the ability of our noise models to shift, when needed, the balance of the refinement towards the accurate molecular mechanics forcefield used by EMMIVox.

Accurate cryo-EM map prediction with Bayesian inference of B-factors. During refinement and validation,

a predictor of cryo-EM density from a 3D model (forward model) is required to measure the agreement of a model with the experimental map. Per-residue B-factors need to be determined to smoothen the prediction obtained with classical forward models derived from Gaussian fits of electron scattering factors and thus to explain fuzzy densities with an individual conformation. To reduce data overfitting, PHENIX adds restraints to avoid B-factors of residues close in space to deviate too much from each other. In presence of sharp transitions between ordered and disordered regions, these additional restraints lead to under (over) estimating the B-factors corresponding to flexible (rigid) residues, effectively compressing the space that B-factors can sample. Inspired by PHENIX, we developed a Bayesian inference approach to determine Bfactors that allows these restraints to be violated, especially in regions corresponding to order-disorder transitions (Materials and Methods; Fig. S9). The B-factors inferred by EMMIVox contribute to a more accurate density prediction across the whole system and ultimately to increase the model fit to the data.

Coarse-grained single-structure refinement

For prospective use in single-structure refinement using medium-low resolution cryo-EM and cryo-ET data, we developed a forward model to predict density maps from coarse-grained Martini 3 models, in which each amino acid is represented by a few beads [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] (Materials and Methods). To evaluate the accuracy of our coarse-grained forward model, we calculated density maps from the Martini representations of a set of 1909 cryo-EM structures and quantified the agreement between the predicted and experimental cryo-EM data.

For comparison, we also calculated density maps from the all-atom structures using the atomistic forward model. Comparing the 𝐶𝐶 !"#$ with the experimental maps given by the coarse-grained and atomistic forward models as a function of experimental resolution revealed that they perform equally well at lower resolutions (>4.0 Å), while the atomistic forward model is more accurate at higher resolutions (<4.0 Å) (Fig. 3). These results confirm the accuracy of our coarse-grained forward model and suggest that Martini models may be useful for single-structure refinement of medium-low resolution cryo-EM and cryo-ET data.

Additionally, our new forward model will be useful for validating and biasing Martini simulations using cryo-EM and cryo-ET data.

Ensemble refinement

When coupled with metainference 41 , EMMIVox can be used to model structural ensembles by interpreting low resolution areas of a cryo-EM map in terms of a mixture of conformational heterogeneity and noise (Materials and Methods). To determine the accuracy of this combined approach, we determined structural ensembles for all the systems in our single-structure refinement benchmark set and compared their quality of fit to the experimental map to the EMMIVox single-structure model (Fig. 4A andS10). The 𝐶𝐶 !"#$ scores for structural ensembles had a median increase of 13% compared to single-structure models, reaching up to a 𝐶𝐶 !"#$ of 0.95 for the ensemble of Escherichia coli PBP1b (3.28 Å, PDB 7lq6). The largest improvement was observed for the Bacteriophage SPP1 (4 Å, PDB 6yeg), whose 𝐶𝐶 !"#$ score increased by 28.4% from 0.60 observed in the single-structure model to 0.77 for the ensemble.

The superior model-data fit is not necessarily due only to the larger number of parameters used in ensemble refinement, i.e. multiple conformations, because single-structure models can implicitly account for dynamic effects using residue-level B-factors to describe fluctuations around an average conformation. The improved ensemble 𝐶𝐶 !"#$ suggests therefore that the dynamics of flexible regions might not be well represented by such Gaussian fluctuations. To investigate this point further, we examined in more detail the nature of the structural diversity observed in our EMMIVox ensembles. Across our entire benchmark set, 48% of residues showed a dynamic behavior that cannot be described by a unimodal distribution (Fig. 4B), neither at the backbone nor side-chain level. When restricting our analysis to medium-to-large scale dynamics, defined by residue Root Mean Square Fluctuations (RMSF) greater than 6 Å, the majority of residues (61%) sampled multimodal conformational distributions impacting a larger portion of backbone atoms as dynamics increases (Fig. 4C). These results, combined with the higher ensemble 𝐶𝐶 !"#$ , indicate that unimodal, Gaussian-like distributions centered around a single-structure model with fluctuations proportional to the B-factors cannot accurately describe the conformational heterogeneity of the backbone and side chains that are averaged out in cryo-EM maps.

To determine whether the EMMIVox ensembles overinterpret experimental noise as conformational heterogeneity, we compared for each system the per-residue B-factors of the single-structure model to the residue MSF within the EMMIVox ensemble (Fig. S11). For all systems we observed a correlation between these two quantities indicating that the majority of the (multi-modal) heterogeneity observed in our ensembles corresponds to an increased B-factor in the single-structure model. However, most systems contained several residues with large B-factors that were not matched by proportionately large MSFs. This suggests that fuzzy densities around these residues correspond to experimental noise rather than structural dynamics. One intriguing example is ChRmine channelrhodopsin (2 Å, PDB 7w9w). In the EMMIVox single-structure model we observed disproportionately high B-factors for residues 191-214 and 269-279 (Fig. 4D, dashed box), which are located on the intracellular side (Fig. 4E). Density around these residues is extremely fuzzy, most likely due to the residual density from the antibody introduced to facilitate image alignment and/or radiation damage as this region is unprotected by the micelle during data collection (Fig. 4E). EMMIVox correctly did not overinterpret this fuzzy density as an extremely dynamic region but generated an ensemble with reduced conformational heterogeneity compared to the single-structure Bfactors (Fig. 4F). This analysis, along with the previous observation that structural variety is often described by multimodal distributions, show that it is possible to get additional information about conformational heterogeneity beyond simply representing the continuous dynamics averaged in cryo-EM maps as Gaussian-like fluctuations proportional to the B-factors.

Structural dynamics from atomistic cryo-EM maps: the GPT type 1a tau filament

As cryo-EM approaches atomistic resolution the question arises as to whether ensemble descriptions of density maps are still necessary [START_REF] Fraser | What Will Computational Modeling Approaches Have to Say in the Era of Atomistic Cryo-EM Data?[END_REF] . We argue that ensemble representations are essential to extract all the information contained in these maps, such as the presence and location of semi-ordered water and lipids around proteins or local side chain flexibility, which cannot be captured by single-structure models nor by 2D particle classification. We demonstrate this concept using as example the GPT type 1a tau filament, which was recently determined at 1.9 Å resolution (PDB 7p6a, EMDB 13223) [START_REF] Shi | Structure-based classification of tauopathies[END_REF] .

Tau 1a proteins, which are thought to play a role in various neurodegenerative diseases, organize into tensof-nm-long filaments with fold-dependent repeating structures. With EMMIVox, we determined the structural ensemble of the tau 1a filament from an atomistic cryo-EM map and observed several interesting features. First, the location of ordered water molecules found in the deposited PDB was correctly identified by the water density within the EMMIVox ensemble without the use of symmetry constraints (Fig. 5AB). Notably, we also observed additional water density close to the ordered waters, which was not present in the deposited PDB. This density corresponds to a second hydration shell composed of semi-ordered waters in exchange with bulk solution. Second, it is known that the aromatic residues at the surfaces of fibrils may be important to fibril growth [START_REF] Daskalov | Contribution of Specific Residues of the beta-Solenoid Fold to HET-s Prion Function, Amyloid Structure and Stability[END_REF] . Recent NMR studies revealed that the dynamics of these aromatic rings vary significantly depending on the amino acid location, with surface exposed amino acids populating a larger number of conformations in rapid exchange compared to aromatics within the fibril core [START_REF] Becker | The Rigid Core and Flexible Surface of Amyloid Fibrils Probed by Magic-Angle-Spinning NMR Spectroscopy of Aromatic Residues[END_REF] . This behavior was also observed in our EMMIVox ensemble (Fig. 5CD). This raises the question as to whether other non-aromatic surface exposed amino acids occupy multiple conformations, and if this information can be extracted from atomistic cryo-EM maps. We calculated the population of the conformers of two surface-exposed residues from the EMMIVox ensemble: K343, which in the deposited PDB was modeled in two alternative conformations with equal occupancy, and K347. We confirm that K343 occupies two distinct conformations with equal probability of 50% (Fig. 5E). Interestingly, we also observed a previously unmodeled minor conformation for K347 with an occupancy of 30% compared to 70% for the major conformation present in the deposited PDB (Fig. 5F).

Discussion

Here we presented EMMIVox, a tool to determine accurate single-structure models as well as structural ensembles by combining MD simulations driven by state-of-the-art molecular mechanics force fields and cryo-EM density maps. EMMIVox incorporates cryo-EM data in the form of voxels using a Bayesian framework that explicitly accounts for data correlation and noise. The method automatically balances model fit to experimental data with stereochemical quality and produces single-structure models that outperform the structures deposited in the PDB in terms of multiple key metrics. EMMIVox systematically reduced the number of serious clashes observed in deposited structures and optimized the positions of backbone and side-chain atoms to improve the description of physico-chemical interactions fundamental for protein stability, such as hydrogen bonds and salt bridges. These single-structure models can be leveraged for all those applications that require high-quality data, such as in silico structure-based drug design or training machine learning approaches.

Even more importantly, EMMIVox can be used to determine structural ensembles that reveal the conformational heterogeneity hidden in low-resolution areas of cryo-EM maps. While 2D classification, manifold embedding and other advanced image processing techniques [START_REF] Scheres | Processing of Structurally Heterogeneous Cryo-EM Data in RELION[END_REF][START_REF] Punjani | cryoSPARC: algorithms for rapid unsupervised cryo-EM structure determination[END_REF][START_REF] Punjani | 3DFlex: determining structure and motion of flexible proteins from cryo-EM[END_REF][START_REF] Herreros | Estimating conformational landscapes from Cryo-EM particles by 3D Zernike polynomials[END_REF][START_REF] Kinman | Uncovering structural ensembles from single-particle cryo-EM data using cryoDRGN[END_REF][START_REF] Chen | Deep learning-based mixed-dimensional Gaussian mixture model for characterizing variability in cryo-EM[END_REF][START_REF] Frank | Continuous changes in structure mapped by manifold embedding of single-particle data in cryo-EM[END_REF] have the potential to identify distinct conformational states directly from the single-particle data, we have shown here that even atomistic maps still present regions in which dynamics of flexible regions is averaged out. EMMIVox relies on both the cryo-EM data and physico-chemical knowledge to determine structural ensembles that capture key information that is lost in single-structure models, such as the presence and population of minor sidechain conformational states, semi-ordered waters, lipids and ligands. Our analysis indicated that this structural variability often corresponds to multimodal conformational distributions that cannot be well described by Gaussian-like fluctuations around an average model and proportional to single-structure B-factors. While these might appear as minor details, we have shown in the past that dynamics of flexible regions extracted from cryo-EM maps often play a crucial role in biological function, for example for ligand recognition in the ASCT2 transporter [START_REF] Garibsingh | Rational design of ASCT2 inhibitors using an integrated experimentalcomputational approach[END_REF] and to understand the effect of post-translational modifications in microtubules [START_REF] Eshun-Wilson | Effects of alpha-tubulin acetylation on microtubule structure and stability[END_REF] .

EMMIVox presents multiple advantages with respect to our previous approach for structure determination from cryo-EM maps (EMMI) [START_REF] Bonomi | Simultaneous Determination of Protein Structure and Dynamics Using Cryo-Electron Microscopy[END_REF] . EMMIVox directly utilizes the density voxels of 3D cryo-EM density maps, replacing the previous representation method using Gaussian Mixture Models (GMM). While representing the cryo-EM density as voxels requires careful removal of correlated datapoints, it significantly increases usability, especially with atomistic cryo-EM density maps, for which determining accurate GMMs is prohibitively expensive. The direct use of voxels also enables estimation of experimental errors from deposited half-maps and the use of noise models for each density voxel, ensuring that various sources of error are not interpreted as structural dynamics when modeling structural ensembles.

Additionally, our novel Bayesian inference approach to determine B-factors, which is unique to EMMIVox, enables accurate single-structure refinement and quantitative comparison of the models to the deposited PDBs using standard metrics of quality.

In spite of significant improvements, EMMIVox still presents some limitations. First, the computational cost of EMMIVox is higher compared to the real-space refinement approach implemented in PHENIX. The main reason is the use of accurate models of the environment, i.e. explicit water molecules, lipid bilayer, and solution ions, and of the intermolecular forces, especially long-range interactions. Furthermore, while single-structure refinement can be performed on a GPU-enabled Desktop computer in 10 to 24 hours, the determination of structural ensembles requires simultaneous use of multiple computer nodes, which may be out of reach for some users. Second, EMMIVox was designed to capture the underlying conformational distribution of the selected single-particle images used to reconstruct the 3D density map. As this data is a subset of the entire particle set, it is important to note that the resulting ensemble may represent local dynamics around a stable macrostate and not the entire conformational landscape. Furthermore, it should be kept in mind that the cryo-EM cooling process might perturb the room temperature ensemble by reducing thermal motion and enabling equilibration into lower free-energy conformations [START_REF] Bock | Effects of cryo-EM cooling on structural ensembles[END_REF] . In EMMIVox, these issues are partially mitigated by using the room temperature ensemble provided by the molecular mechanics force field as a prior. Third, while EMMIVox has been shown to accurately and rapidly capture conformations present in available cryo-EM density maps, including those which are difficult to observe on short time scales in standard MD simulations (Fig. S12), EMMIVox may have trouble sampling conformations separated by large free-energy barriers. In this case, sampling issues can be alleviated by combining EMMIVox with various enhanced sampling techniques, such as metadynamics as recently proposed in the MEMMI approach [START_REF] Brotzakis | Determination of the Structure and Dynamics of the Fuzzy Coat of an Amyloid Fibril of IAPP Using Cryo-Electron Microscopy[END_REF] . Finally, because of these sampling limitations, an initial model that already fits the experimental density to a good extent is required as starting point. AlphaFold2 9 , ab initio machine learning and other automated model building techniques [START_REF] Giri | Deep learning for reconstructing protein structures from cryo-EM density maps: Recent advances and future directions[END_REF][START_REF] Jamali | Automated model building and protein identification in cryo-EM maps[END_REF] , or flexible fitting methods such as MDFF [START_REF] Trabuco | Flexible fitting of atomic structures into electron microscopy maps using molecular dynamics[END_REF] , Tempy-REFF [START_REF] Cragnolini | Cryo-EM structure and B-factor refinement with ensemble representation[END_REF] , and the novel maximum likelihood approach implemented in GROMACS [START_REF] Blau | Gentle and fast all-atom model refinement to cryo-EM densities via a maximum likelihood approach[END_REF] , can provide excellent starting models for accurate single-structure and/or ensemble refinement with EMMIVox. Despite these limitations, EMMIVox enables determining accurate structural and dynamic models using in vitro data, and holds promise as a valuable tool for future efforts focused on structure determination in situ.

As the resolution of cryo-ET improves (Fig. S1B), high-quality in cell data will become readily available and integrative approaches that combine various sources of experimental and computational data will be required to obtain accurate structural models. The coarse-grained forward model implemented in EMMIVox will enable refining models of large macromolecular architectures from sub-nanometer subtomogram averaging data. Additionally, the Bayesian framework on which EMMIVox is built enables automatic weighting of multiple sources of in silico and experimental data, making it a perfect framework for integrative structure and dynamic determination in cell. This integration of a diverse set of data is further facilitated by the implementation of EMMIVox in the PLUMED-ISDB 69 module of PLUMED [START_REF] Tribello | PLUMED 2: New feathers for an old bird[END_REF] , which makes it possible to combine a wide range of different types of experimental data.

In summary, EMMIVox represents a flexible instrument to convert cryo-EM data into high-quality structural models. Our approach can be used to reveal dynamic properties of proteins, lipids, ligands, waters, and ions by extracting information from cryo-EM density maps that would otherwise be lost. These models can advance our understanding of the molecular mechanisms that drive biological functions, and provide valuable information for structure-based drug design and training of novel machine learning approaches. Looking ahead, we envision that EMMIVox will be an indispensable part of integrative structural biology pipelines and will contribute to obtaining a more complete picture of highly intricate and dynamic systems in biologically relevant environments.

Materials and Methods

Theory of EMMIVox for single-structure refinement

General overview. EMMIVox is based on a Bayesian inference framework 37 that estimates the probability of a model 𝑀 given the information available about the system, including prior physico-chemical knowledge and newly acquired experimental data. The posterior probability 𝑝(𝑀|𝐷) of model 𝑀, which is defined in terms of its structure 𝑋 and other parameters, given data 𝐷 and prior knowledge is:

𝑝(𝑀|𝐷) ∝ 𝑝(𝐷|𝑀) • 𝑝(𝑀) (1)
where the likelihood function 𝑝(𝐷|𝑀) is the probability of observing data 𝐷 given 𝑀 and the prior 𝑝(𝑀)

is the probability of model 𝑀 given the prior information. At variance with the previous EMMI approach [START_REF] Bonomi | Bayesian Weighing of Electron Cryo-Microscopy Data for Integrative Structural Modeling[END_REF][START_REF] Bonomi | Simultaneous Determination of Protein Structure and Dynamics Using Cryo-Electron Microscopy[END_REF] , here we define the experimental data as a set of density values 𝐷 = {𝑑 % } observed in the voxels of a cryo-EM three-dimensional (3D) map. To define the likelihood function, one needs i) a forward model 𝑓 % (𝑋) to predict the density that would be observed in voxel i for structure 𝑋 in absence of noise, and ii) a noise model that defines the distribution of deviations between observed and predicted data.

Atomistic forward model.

To predict the density in voxel i of a cryo-EM map from an atomistic model, we used the fast Fourier transform of a 5-Gaussian fit of the electron scattering factors [START_REF] Peng | Robust parameterization of elastic and absorptive electron atomic scattering factors[END_REF] :

𝑓 % (𝑋, 𝑏 & ) = / / 𝐴 ',$ 1 𝜋 𝐵 ',$ + 𝑏 ' & /4 7 ) * exp 1- 𝜋 * <𝑋 ' -𝑉 % < * 𝐵 ',$ + 𝑏 ' & /4 7 + $,- . ',- (2)
where 𝑋 ' and 𝑉 % are the coordinates of atom j and the center of the voxel i, respectively. The Gaussians parameters 𝐴 $ and 𝐵 $ depend on the atom type (Tab. S1). The B-factors 𝑏 & = 2𝑏 𝑗 𝑓 3, which are defined for each atom j but are identical within individual residues, enable smoothing the forward model prediction to describe low-resolution regions with fuzzy density using a single-structure model. The external sum runs over all the non-hydrogen atoms, with exclusion of the carboxylate oxygens of glutamic and aspartic acid, as these groups are often damaged by the electron beam. To reduce the computational cost, we further restricted this summation to a neighbor list of atoms with distance cutoff from the voxel center of 1.0 nm.

Noise model. We use a Gaussian noise model for the density 𝑑 % observed in voxel i.:

𝑝(𝑑 % |𝑋, 𝑏 & , 𝜎 % ) = 1 𝜎 % √2𝜋 exp C- [𝑑 % -𝛼 ⋅ 𝑓 % (𝑋, 𝑏 & )] * 2 𝜎 % * H (3)
where the uncertainty parameter 𝜎 % * quantifies the level of experimental noise as well as errors in the forward model. These variable parameters allow to dynamically adjust the intensity of the cryo-EM spatial restraints during model refinement based on the accuracy in the data, which is inferred by the consistency between cryo-EM data, physico-chemical prior and potentially additional experimental data. High-noise voxels (outliers) are therefore automatically identified and downweighed in the refinement of the structural model in favor of the prior information. 𝛼 is a constant scaling factor between the entire set of voxels of the predicted and experimental maps, and it is optimized prior to production runs (Scale factor optimization).

Data correlation and likelihood function.

The data likelihood 𝑝(𝐷|𝑀) is typically expressed as a product of individual likelihoods 𝑝(𝑑 % |𝑀), one for each experimental data point, under the assumption that these points are independent. In case of cryo-EM data, this assumption does not hold because neighboring voxels of a 3D map contain correlated information. Ignoring such correlation would lead to overcounting the number of independent data points and ultimately biasing the refinement towards overfitting the data at the expenses of the stereochemical quality of the models. To reduce data correlation, we developed the following pre-filtering procedure:

1. We first selected all the voxels with non-negative density around 3.5 Å of the structure resolved in the deposited PDB. This includes protein atoms as well as ordered waters, ions, small-molecules and lipids. We sorted the selected voxels in descending order based on the value of the density: this constitutes our initial pool of data points 𝐷 = {𝑑 % };

2. We started from the first voxel 𝑑 -(highest density) in pool 𝐷 and calculated the density autocorrelation function for displacements up to a few voxels in the x, y, and z direction starting from this location. For each displacement, the autocorrelation function is calculated by averaging over multiple starting voxels in a cubic minibox of side equal to 4 Å and centered on 𝑑 -.

Equivalently, the autocorrelation can be expressed as the Pearson correlation coefficient of a series of density values and its space-lagged version;

3. We removed from pool 𝐷 all the voxels within the cubic minibox with Pearson correlation coefficient greater than a predefined threshold, with the exception of 𝑑 -.

4. We moved to the next voxel 𝑑 * in pool 𝐷 in descending order of density and re-applied the filtering procedure of step 2 and 3;

5. We iterated until reaching the voxel in pool 𝐷 with lowest density.

It should be noted that sub-sampling procedures are commonly used in structural modelling with experimental data, for example to remove correlation between points in Small-Angle X-ray scattering profiles [START_REF] Paissoni | Determination of Protein Structural Ensembles by Hybrid-Resolution SAXS Restrained Molecular Dynamics[END_REF] . PHENIX itself utilizes only the density in the positions occupied by the atoms during each step of refinement, which is obtained by interpolating the density in the voxels surrounding each atom [START_REF] Afonine | Real-space refinement in PHENIX for cryo-EM and crystallography[END_REF] . After pre-filtering the map with the procedure describe above, we assume that the 𝑁 1 selected voxels can be considered as independent data points and write the total likelihood function as:

𝑝(𝐷|𝑀) = K 𝑝(𝑑 % |𝑋, 𝑏 & , 𝜎 % ) . ! %,- (4)
Priors. Experimental cryo-EM density maps may contain errors from various sources including systematic errors during 3D reconstruction, noise from the use of a limited number of images taken at various sample orientations, and radiation damage from exposure to the electron beam. Furthermore, the forward model used to predict a map from a structural model is intrinsically inaccurate. While these errors might be difficult to measure, it is critically important that they are not misinterpreted as structural dynamics during the generation of structural ensembles. To guide the inference of the uncertainty parameters 𝜎 % that quantify the noise level in the map, we defined a lower bound 𝜎 % !%2 for each voxel by calculating the voxel-by-voxel variation between two independent reconstructions, or half maps. While this quantity does not describe all the sources of error at play, we expect that the total error in a given voxel cannot be lower than 𝜎 % !%2 . We imposed this lower bound using the following prior for the uncertainty parameters 𝜎 % 72 :

𝑝(𝜎 % ) = C 0 𝑓𝑜𝑟 𝜎 % ≤ 𝜎 % !%2 1/𝜎 % 𝑓𝑜𝑟 𝜎 % > 𝜎 % !%2 (5)
where 1/𝜎 % is a typical Jeffreys prior.

Inspired by PHENIX, we also added a restraint to guide B-factors inference and avoid that residues close in space have B-factors significantly different from each other. These restraints are enforced by:

𝑝(𝑏 & , 𝜎 & ) = K 1 𝜎 '$ √2𝜋 exp R- S𝑏 ' & -𝑏 $ & T * 2 𝜎 '$ * U (6) ',$
where the product is over all pairs of B-factors corresponding to residues for which at least a pair of atoms is closer than 5 Å. At variance with the B-factor restraint implemented in PHENIX, our Bayesian approach tolerates outliers, i.e. pair of close residues with significantly different B-factors, via the introduction of the uncertainty parameters 𝜎 & = W𝜎 '$ X. This approach allows us to better model during single-structure refinement those regions of the system that undergo a sudden transition between order and disorder.

Finally, as structural prior 𝑝(𝑋) we used state-of-the-art atomistic force fields 𝐸 33 (𝑋), which enabled us to accurately model the system as well as the environment, including explicit water molecules, ions, smallmolecules, and lipids:

𝑝(𝑋) = exp[-𝐸 33 (𝑋)/𝑘 4 𝑇] (7)
where 𝑘 4 is the Boltzmann constant and 𝑇 the temperature of the system.

Marginalization.

To avoid sampling the uncertainty parameters 𝜎 % and 𝜎 & , we marginalized the corresponding distributions. The resulting marginal data likelihood is:

𝑝(𝑑 % |𝑋, 𝑏 & ) = ] 𝑝(𝑑 % |𝑋, 𝑏 & , 𝜎 % ) 1 𝜎 % 𝑑𝜎 % 56 7 " #"$ = 0.5 𝑑 % -𝛼 ⋅ 𝑓 % (𝑋, 𝑏 & ) erf a 𝑑 % -𝛼 ⋅ 𝑓 % (𝑋, 𝑏 & ) √2 𝜎 % !%2 b (8)
while the B-factors prior becomes upon introduction of a Jeffreys prior 1/𝜎 𝑗𝑘 :

𝑝(𝑏 & ) = K ] 1 𝜎 '$ √2𝜋 exp R- S𝑏 ' & -𝑏 $ & T * 2 𝜎 '$ * U 1 𝜎 '$ 𝑑𝜎 '$ = K 0.5 𝑏 ' & -𝑏 $ & erf 1 𝑏 ' & -𝑏 $ & √2 𝜎 9& !%2 7 ',$ 56
7 %& #"$ ',$ (9)
where 𝜎 9& !%2 is a lower bound for the B-factor uncertainty parameters set equal to 0.1 nm 2 .

EMMIVox hybrid energy function. After defining all the component of our approach and marginalizing the uncertainty parameters, we obtain the final EMMIVox posterior distribution: (10) where the marginal data likelihood is given by Eq. 8, the B-factors prior by Eq. 9, and the structural prior by Eq. 7. To sample the posterior, we define the associated hybrid energy function as:

𝑝(𝑋, 𝑏 & e𝐷) ∝ K 𝑝(𝑑 % |𝑋, 𝑏 & ) . ! %,- ⋅ 𝑝(𝑏 & ) ⋅ 𝑝(𝑋)
𝐸(𝑋, 𝑏 & ) = -𝑘 4 𝑇 log 𝑝(𝑋, 𝑏 & e𝐷) = 𝐸 33 (𝑋) + 𝐸 :;<=>?@ (𝑋, 𝑏 & ) + 𝐸 9& (𝑏 & ) ( 11
)
The EMMIVox hybrid energy is therefore decomposed into: i) the molecular mechanics force field 𝐸 33 , ii)

the spatial restraints 𝐸 :;<=>?@ to enforce the model agreement with the cryo-EM map, and iii) the restraints 𝐸 9& to guide B-factors determination. A Gibbs sampling scheme is used to sample model coordinates with Molecular Dynamics (MD) and B-factors parameters with Monte Carlo (MC).

Development of the coarse-grained forward model

Parameterization. To enable the use of EMMIVox with coarse-grained representations, we developed a forward model to calculate density maps from Martini 3 protein structures [START_REF] Souza | Martini 3: a general purpose force field for coarse-grained molecular dynamics[END_REF] . We parameterized the model by fitting to predictions of the atomistic forward model. In the atomistic forward model, the atomic electron scattering factors are given by the Fourier transform of a 5-gaussian mixture (see section Atomistic forward model and Eq. 2). We used the same framework for the Martini forward model, but with bead electron scattering factors given by a single Gaussian, i.e. only two parameters Ak and Bk for each Martini bead type k. We determined an individual set of Ak and Bk for each bead position in each type of amino acid for a total of 52 different bead types using the following procedure. To capture the conformational variation in the atoms mapping to each bead type, we used a set of 2906 structures to fit the parameters. These were all the structures obtained with single-particle cryo-EM in the period 2020-2023, with resolution ranging from 1.2 Å to 11 Å. We divided the structures into groups of atoms, each corresponding to one Martini bead type, skipping any beads with missing atoms. For each group, we defined a cubic box of voxels centered on each atom and with side equal to 6 Å (for a total of 6859 evenly spaced voxels). In this set of voxels, we calculated a density map from the group of atoms using the atomistic forward model and from the corresponding Martini bead using the coarse-grained forward model with a grid scan of Ak and Bk. We scanned values of Ak from 0 to 8.0 in steps of 0.1 and Bk from 0 Å 2 to 40.0 Å 2 in steps of 0.5 Å 2 . We evaluated the agreement between the two forward models as a function of Ak and Bk using the squared error summed over all the voxels:

𝑀𝑆𝐸(𝐴 $ , 𝐵 $ ) = 1 𝑁 1 ?@𝑓 % @";A%2% (𝐴 B , 𝐵 $ ) -𝑓 % "A=!%#A%: B * (12)
. !

%

where NV is the number of voxels, 𝑓 % @";A%2% and 𝑓 % "A=!%#A%: are the voxel densities predicted by the Martini and atomistic forward models, respectively. For each of the 52 Martini bead types, we averaged the MSE grid scan over the bead type's occurrences in all 2906 structures and selected the set of Ak and Bk that minimized the average MSE (Tab. S2).

Validation.

To validate the Martini forward model and compare its accuracy with that of the atomistic forward model, we evaluated the agreement between predicted and experimental density maps for 1909 cryo-EM structures. These were selected from our initial set of 2906 structures as those that could be easily mapped to the Martini 3 representation using the martinize2 python script without any additional modifications to the structures [START_REF] Kroon | Martinize2 and Vermouth: Unified Framework for Topology Generation[END_REF] . We calculated density maps from the atomistic and Martini representations using the respective forward models and calculated the 𝐶𝐶 !"#$ with the experimental density maps. We used the same protocol described below (in Details of the single-structure refinement benchmark) to preprocess the experimental maps and to optimize the B-factors and scale factors of the predicted maps to maximize 𝐶𝐶 !"#$ . The B-factors were sampled for 1000 MC steps and the scale factor was scanned from 0.7 to 1.3 in steps of 0.05. No data correlation filtering was used for the experimental maps. To evaluate the difference between the forward models as a function of the experimental resolution, we binned the structures by resolution (bin size equal to 1.0 Å) and calculated for each bin the average difference in 𝐶𝐶 !"#$ between the atomistic and Martini forward models:

Δ𝐶𝐶 "A=!%#A%:>@";A%2% = 1 𝑁 C ?@𝐶𝐶 % "A=!%#A%: -𝐶𝐶 % @";A%2% B (

. '

%

where NS is the number of structures in the bin.

Theory of EMMIVox for ensemble refinement

To model structural ensembles fitting a cryo-EM density map, we coupled EMMIVox with metainference [START_REF] Bonomi | Metainference: A Bayesian inference method for heterogeneous systems[END_REF] .

Metainference is a general Bayesian inference approach that enables modelling structural ensembles using any kind of ensemble-averaged experimental data as well as prior physico-chemical information. Inspired by the maximum entropy/replica-averaged approach [START_REF] Cavalli | Molecular dynamics simulations with replicaaveraged structural restraints generate structural ensembles according to the maximum entropy principle[END_REF] , the metainference posterior is expressed in terms of several copies of the system (or replicas), which represent the conformational heterogeneity of the system.

With a Gaussian data likelihood, as in the EMMIVox case, the metainference posterior is:

𝑝 @D @{𝑋 ; }, F𝜎 ;,% G H𝐷B ∝ J 𝑝(𝑋 ; ) ⋅ J 1 √2𝜋𝜎 ;,% exp Q- [𝑑 % -𝑓 % ({𝑋 ; } )] * 2𝜎 ;,% * T %,- . ( ;,- 𝑝(𝜎 ;,% ) (14)
where {𝑋 ; } are the structures of the 𝑁 E replicas of the system and F𝜎 ;,% G the error parameters, one per replica 𝑟 and experimental data point 𝑖. The predicted experimental observable 𝑓 % ({𝑋 ; }) is calculated as the average of 𝑓 % over the 𝑁 E replicas. Sampling of the metainference posterior is performed by a multi-replica MD/MC simulation driven by the hybrid energy function associated to Eq. 14. This approach has been extensively used to determine structural ensembles of highly dynamic and disordered systems using NMR spectroscopy [START_REF] Heller | Small-molecule sequestration of amyloid-beta as a drug discovery strategy for Alzheimer's disease[END_REF][START_REF] Heller | Sequence Specificity in the Entropy-Driven Binding of a Small Molecule and a Disordered Peptide[END_REF] as well as SAXS/SANS 77, 78 data. Metainference was also used to determine structural ensembles from cryo-EM data in our previous EMMI approach based on a Gaussian Mixture Models representation of the cryo-EM map [START_REF] Bonomi | Bayesian Weighing of Electron Cryo-Microscopy Data for Integrative Structural Modeling[END_REF][START_REF] Bonomi | Simultaneous Determination of Protein Structure and Dynamics Using Cryo-Electron Microscopy[END_REF] . In combination with EMMIVox and the more accurate noise models developed here, metainference enables ensemble-refinement by interpreting low resolution areas of a cryo-EM map in terms of a mixture of conformational heterogeneity and noise.

Details of the single-structure refinement benchmark

Details of the systems. Nine systems determined at a resolution ranging from 1.9 Å to 4.0 Å were chosen to benchmark the EMMIVox single-structure refinement (Tab. S3): the GPT type 1a tau filament, the ChRmine channelrhodopsin, the Anaplastic lymphoma kinase extracellular domain fragment in complex with an activating ligand, a mutant of dimeric unphosporylated Pediculus humanus protein kinase, the human CDK-activating kinase bound to the clinical inhibitor ICEC0942, the major facilitator superfamily domain containing 2A in complex with LPC-18:3, the Escherichia coli PBP1b, the mammalian peptide transporter PepT2, and the SPP1 bacteriophage tail tube.

Setup and general MD details.

Missing residues in the deposited PDB were modelled using GalaxyFill [START_REF] Coutsias | A kinematic view of loop closure[END_REF] or, in case missing residues could not be placed, with Modeller 15 v. 10.1. The resulting model was then processed using the CHARMM-GUI 80 server. Membrane proteins (PDB ids 7w9w, 7mjs, and 7nqk) were inserted in a homogeneous POPC lipid bilayer. Each system was solvated in a triclinic box with dimensions chosen in such a way that the edge of the box was 1.0 nm away from the closest model atom. K+ and Clwere added to ensure charge neutrality at concentration equal to 0.15 M. The CHARMM36m 81 forcefield was used for proteins and lipids, CgenFF for small molecules [START_REF] Vanommeslaeghe | CHARMM General Force Field: A Force Field for Drug-Like Molecules Compatible with the CHARMM All-Atom Additive Biological Force Fields[END_REF] , and the mTIP3P 83 model for water molecules. In all simulations the equations of motion were integrated by a leap-frog algorithm with timestep equal to 2 fs. The smooth particle mesh Ewald 84 method was used to calculate electrostatic interactions with a cutoff equal to 1.2 nm. Van der Waals interactions were gradually switched off at 1.0 nm and cut off at 1.2 nm. All simulations were carried out using GROMACS 85 v. 2020.5 equipped with the development version of PLUMED 38 . To optimize performances, EMMIVox is implemented in PLUMED with libtorch [START_REF] Paszke | PyTorch: An Imperative Style, High-Performance Deep Learning Library[END_REF] to efficiently calculate the cryo-EM forward model and the hybrid energy function on the GPU.

Cryo-EM map preprocessing.

For each system, we downloaded the cryo-EM map as well as the two half maps from the EMDB database. We applied our pre-processing procedure to select the voxels within 3.5 Å of the model atoms and filter them to reduce data correlation, with a threshold of 0.7, 0.8, 0.9, and 1.0 (no filtering). The two experimental half maps were used to calculate the lower bound for the density uncertainty parameters. Since CHARMM-GUI processing of the deposited PDB often translates and rotates the input conformation, we calculated the transformation that aligns initial and final models and applied it to all the voxels selected for refinement. Finally, a single map data file was created with the list of voxels to be used by PLUMED for refinement and, for each voxel, the following information: (transformed) coordinates of the voxel center, density value, and uncertainty lower bound.

Equilibration. Energy minimization was performed on each system using the steepest decent approach. A 1 ns-long NPT equilibration was then performed with the Bussi-Donadio-Parrinello thermostat [START_REF] Bussi | Canonical sampling through velocity rescaling[END_REF] and the Berendsen barostat 88 , set at 300 K and 1 atm respectively. In systems containing a lipid bilayer, the pressure coupling type was semi-isotropic to allow deformations in the xy plane independent from the z-axis. Next, NVT simulations were carried out for 2 ns using the Bussi-Donadio-Parrinello thermostat at 300 K. During the last two steps of equilibration, positional restraints were applied to all the heavy atoms of the protein and any other components contributing to the predicted cryo-EM density map.

Modelling of ordered waters. In presence of ordered waters in the deposited PDB, a special treatment was required. First, resolved waters as well as all the water molecules within 3.5 Å in the initial CHARMM-GUI solvated model (buffer waters) contributed to the map prediction via our forward model. Second, the positions of ordered waters were restrained during equilibration; the buffer waters were restrained to stay within 8 Å from a reference protein atom, defined as the closest to each water molecule in the initial model.

Scale factor optimization.

To reduce the number of free parameters to sample during the production run, we determined an optimal scaling factor between predicted and observed cryo-EM maps. We analyzed with the PLUMED driver tool the trajectory obtained with positional restraints during the NVT equilibration for different values of the scaling factors in the range from 0.5 to 1.5 at intervals of 0.05. For each value of the scaling factor, we sampled with MC the B-factors while re-reading the trajectory, with a maximum MC move per B-factor equal to 0.05 nm 2 . To optimize sampling, B-factors were initialized using an empirical relationship between map resolution (res, in nm) and average B-factor values (in nm 2 ) determined from 8000 deposited PDBs obtained from cryo-EM maps with resolution less than 5 Å:

𝑏 & = 6.95408 * 𝑟𝑒𝑠 * -0.01 * 2.45697 (15) The scaling factor that resulted in the lowest EMMIVox hybrid energy along the entire trajectory was selected to be used in the production simulations described in the following section.

Production. All single-structure refinement production simulations were performed in the NVT ensemble with Bussi-Donadio-Parrinello thermostat at 300 K for 20 ns (NPT ensemble with Parrinello-Rahman barostat 89 at 1 atm for membrane proteins). Before calculating the EMMIVox hybrid energy, structural discontinuities due to Periodic Boundary Conditions (PBC) were fixed on-the-fly by PLUMED.

Furthermore, to optimize performances, we: i) updated the forward model neighbor list every 50 MD steps;

ii) sampled the B-factors every 500 MD steps with maximum MC move equal to equal to 0.05 nm 2 ; iii) used a multiple-time step algorithm for PBC reconstruction and EMMIVox hybrid energy calculation with stride equal to 4 MD steps [START_REF] Ferrarotti | Accurate Multiple Time Step in Biased Molecular Simulations[END_REF] . The system trajectory was saved to file every 10 ps for subsequent analysis.

Ordered and buffer waters that contributed to the map prediction were restrained to stay within 8 Å from a reference protein atom, defined as the closest to each water molecule in the initial model. These restraints allowed: i) exchanges between the ordered molecules resolved in the PDB and their surrounding buffer waters; ii) to identify additional sites of (semi)ordered molecules in the proximity of the resolved waters, for example a second coordination shell.

Analysis.

To generate the final single-structure refined model, we first extracted the conformation with lowest EMMIVox hybrid energy from the production run as well as the associated B-factors. We then performed an (EMMIVox hybrid) energy minimization using the steepest decent approach. During minimization, B-factors were sampled starting from those found in the initial model using a modified MC sampler in which only downhill moves were accepted. B-factors were sampled every 100 steps with a maximum MC move equal to 0.1 nm 2 . The forward model neighbor list was updated at every step. Next, we selected the conformation at the end of the minimization and fixed discontinuities in the structure due to PBC with the PLUMED driver tool. We then generated a PDB file containing only the heavy atoms of the systems that were used to predict the cryo-EM density map during production. This model was realigned to the original cryo-EM map downloaded from the EMDB using the inverse transformation computed in the section Cryo-EM map preprocessing. Finally, we added to the PDB the B-factors obtained at the end of the final minimization. MolProbity was used to compute the clashscore and MolProbity score, while PHENIX v. 1.15.2 was used to evaluate the fit to the experimental map with the EMRinger score.

We implemented the calculation of 𝐶𝐶 !"#$ 51 in a GPU-enabled python script and generalized it to compute this score from a structural ensemble. At variance with PHENIX, we did not optimize an isotropic B-factor before the 𝐶𝐶 !"#$ calculation. Despite this difference, the 𝐶𝐶 !"#$ calculated by PHENIX and by our python script are strongly correlated (Fig. S13). For all the fit-to-data calculations, the original map as downloaded from the EMDB was used, regardless of the data correlation cutoff and the total number of voxels used in the refinement. For the analysis of physico-chemical fingerprints (Fig. 2), we used MDAnalysis 91 v. 2.0 to calculate the number of hydrogen bonds from a single-structure model with donoracceptor distance and angle cutoff equal to 3.0 Å and 150 degrees, respectively. Hydrogen bonds between anionic carboxylate (RCOO -) of either aspartic acid or glutamic acid and the cationic ammonium (RNH3 + ) from lysine or arginine were classified as salt bridges. Clashes between atoms with overlap between van der Waals radii greater than 0.4 Å were calculated with PHENIX.clashscore.

Details of the ensemble refinement benchmark

Setup. We determined EMMIVox structural ensembles fitting the cryo-EM map for all the systems selected for the single-structure refinement benchmark. To prepare the ensemble simulations, we first extracted the conformation with lowest EMMIVox hybrid energy from the single-structure refinement production run and identified the minimum B-factor 𝑏 !%2 & across all residues of this conformation. We then extracted 16 frames from the second half of the single-structure refinement production runs equally distributed in time.

These conformations will be used as starting points of the production run described in the following section.

Production. EMMIVox ensemble simulations were performed with the same settings used in singlestructure refinement, except for the fact that during ensemble refinement B-factors were not sampled but kept constant and all equal to the 𝑏 !%2 & identified in the previous step. This value of B-factors can be considered as either a baseline dynamic of the most rigid residue of the system or a minimum noise level of the best resolved residue. In both cases, using such constant B-factor during ensemble refinement contributes to avoiding overinterpreting noise in terms of conformational heterogeneity. 16 metainference replicas were used, each one simulated for 20 ns, resulting in 320 ns of ensemble trajectory for each system.

Analysis. After the EMMIVox ensemble simulations were completed, the trajectories of all the metainference replicas were concatenated, resulting in the complete structural ensemble of the system. After fixing with PLUMED the structural discontinuities due to PBC, we re-aligned the EMMIVox ensemble to the original cryo-EM map downloaded from the EMDB using the inverse transformation computed in the section Cryo-EM map preprocessing. We then evaluated the fit to the data by calculating the 𝐶𝐶 !"#$ of the average cryo-EM map from the EMMIVox ensemble and the experimental map. To make a fair comparison, we selected the voxels for the 𝐶𝐶 !"#$ calculation based only on the single-structure model using the standard convention [START_REF] Afonine | New tools for the analysis and validation of cryo-EM maps and atomic models[END_REF] . In this way, the calculation of 𝐶𝐶 !"#$ for single-structure and ensemble models was performed on the same set of voxels. The analysis of multimodality of the conformational distribution of each individual residue (Fig. 4) was performed using the Folding Test of Unimodality implemented in libfolding [START_REF] Siffer | Are your data gathered? The Folding Test of Unimodality[END_REF] (https://github.com/asiffer/python3-libfolding).

Software and data availability

EMMIVox is implemented as a part of the PLUMED-ISDB

List of Figures

Fig. 1 .Fig. 2 .

 12 Fig. 1. EMMIVox single-structure refinement benchmark. a-i) Assessment of the quality of the EMMIVox single-structure model (colored bars) and deposited PDB (grey bars) for each of the nine systems in our benchmark set. Models are evaluated based on their stereochemical quality and fit with the experimental cryo-EM map. The stereochemical quality is measured by two metrics: i) clashscore, which corresponds to the number of serious clashes per 1000 atoms, and ii) MolProbity score, which is a global measure of quality that combines clashscore, percentage of Ramachandran dihedrals in non-favored regionsand percentage of bad sidechain rotamers. For both metrics, lower values correspond to higher quality models. The model fit to the data is measure by i) 𝐶𝐶 !"#$ , which is the cross-correlation between experimental map and map predicted from the model in the proximity of the structural model, and ii) EMRinger, which measures the precise fitting of an atomic model into the map. For both metrics, higher values correspond to models that better fit the data.

Fig. 3 Fig. 4 .Fig. 5 .

 345 Fig. 3 EMMIVox coarse-grained refinement. a) Difference between atomistic and coarse-grained (Martini 3) forward model in the agreement with 1909 experimental cryo-EM density maps. The plot shows the average difference in 𝐶𝐶 !"#$ (Δ𝐶𝐶) given by the two forward models as a function of the experimental resolution. Error bars show the standard deviation. b-c) Example structures and cryo-EM density maps at medium (b: PDB 7nnh 93 , 4.0 Å) and high (c: PDB 7o6q 94 , 1.88 Å) resolution. From left to right: experimental density maps and density maps predicted by the atomistic and coarse-grained forward models.The predicted maps are colored by local cross correlation to the experimental map (low: red, high: blue) and the global 𝐶𝐶 !"#$ with the experimental map is shown.

  [START_REF] Bonomi | Integrative structural and dynamical biology with PLUMED-ISDB[END_REF] module in the development version (GitHub master branch) of PLUMED[START_REF] Tribello | PLUMED 2: New feathers for an old bird[END_REF] (https://github.com/plumed/plumed.github.io). The GROMACS topologies and PLUMED input files used in our benchmark as well as the EMMIVox refined single-structure models

	are available in PLUMED-NEST (www.plumed-nest.org), the public repository of the PLUMED
	consortium 39 , as plumID:23.041. Scripts to prepare and analyze EMMIVox simulations as well as complete
	tutorials for single-structure and ensemble refinement are available on GitHub
	(https://github.com/maxbonomi/EMMIVox).

Acknowledgements M.B. and S.E.H. acknowledge the support of the French Agence Nationale de la Recherche (ANR), under grant ANR-20-CE45-0002 (project EMMI). S.E.H. is founded by a Roux-Cantarini fellowship from the Institut Pasteur (Paris, France). M.B. would like to acknowledge PRACE for awarding access to Piz Daint at CSCS, Switzerland.

Ethics declarations

The authors declare no conflict of interests.