Gentle and fast all-atom model refinement to cryo-EM densities via Bayes’ approach

Better detectors and automated data collection have generated a flood of high-resolution cryo-EM maps, which in turn has renewed interest in improving methods for determining structure models corresponding to these maps. However, automatically fitting atoms to densities becomes difficult as their resolution increases and the refinement potential has a vast number of local minima. In practice, the problem becomes even more complex when one also wants to achieve a balance between a good fit of atom positions to the map, while also establishing good stereochemistry or allowing protein secondary structure to change during fitting. Here, we present a solution to this challenge using Bayes’ approach by formulating the problem as identifying the structure most likely to have produced the observed density map. This allows us to derive a new type of smooth refinement potential - based on relative entropy - in combination with a novel adaptive force scaling algorithm to allow balancing of force-field and density-based potentials. In a low-noise scenario, as expected from modern cryo-EM data, the Bayesian refinement potential outperforms alternatives, and the adaptive force scaling appears to also aid existing refinement potentials. The method is available as a component in the GROMACS molecular simulation toolkit.

coordinates of individual atoms into the maps [5]. This enables understanding of e.g. binding site properties, interactions with lipids or other subunits, structural rearrangements between alternative conformations, and in particular it makes it possible to model structural dynamics on nanosecond to microsecond time scales via molecular dynamics (MD) simulation [6]. If the interaction descriptions (force fields) used in these simulations were perfect and one had access to infinite amounts of sampling, computational methods should be able to further improve the structure just by starting refinement from a rough initial density, but in practice both force fields and sampling have shortcomings. Nevertheless, it remains an attractive idea to combine the best of both worlds by using cryo-EM data as large-scale constraints while force fields are employed to fine-tune details -in particular details such as local stereogeometry or interactions on resolution scales that go beyond what the cryo-EM data can resolve. Cryo-EM data and stereochemical constraints have been combined favorably in the past to aid structure modeling into three-dimensional cryo-EM densities either by adding force field terms enforcing desired stereochemistry to established modeling tools [7][8][9][10][11] or by adding a heuristic density-based biasing potential to molecular dynamics simulations [12,13] or elastic network models [14].
In practice, it is not straightforward how to best combine experimental map data with simulations and achieve both satisfactory results and rapid convergence. Density-based biasing potentials can in principle achieve arbitrarily good fits to a map, but it comes at a cost of distorting the protein structure. To address the challenges of balancing desired stereochemical properties with cryo-EM data, refinement protocols have been expanded to include secondary structure restraints [12], multiple resolution ranges [11,15], as well as multiple force constants [11,16]; the latter two either consecutively in individual simulations [11] or via Hamiltonian replica exchange [15,16]. A common challenge of all these approaches is the increase in the ruggedness of the applied bias potential function as the resolution of the cryo-EM density maps increases, and how to correctly balance molecular mechanics and forces from the biasing potential. This leads to an apparent modeling paradox that further improving structural models for cryo-EM densities with molecular dynamics appears to be harder the more high-resolution data is available for these models.
To attack this challenge from a fundamental standpoint, Bayes approach has been used to derive probabilities for all-atom structural models given a cryo-EM density [17] and to weigh cryo-EM data influence against other sources of data [18]. These modeling approaches offer valuable insight into the data content in cryo-EM maps and provide promising new ways to model cryo-EM densities by treating them as generic experimental data. However, they do not reflect the underlying physics of data acquisition and density reconstruction from micrographs and have previously not yielded refinement potentials of a new quality to be applicable, e.g., in molecular dynamics simulations.
One way of circumventing the number of model assumptions that are necessary to reflect the reconstruction of three-dimensional cryo-EM densities is to employ Bayesian models that directly connect micrographs and all-atom ensembles [19]. These attempts have previously proven to be prohibitively costly as a way to derive driving forces for molecular dynamic simulation because projections of model atom coordinates onto millions of cryo-EM particle images (i.e., images of molecules) are required for a single force evaluation.
In this work, we show how it is possible to borrow the highly successful approach to density reconstruction and use Bayesian modeling of cryo-EM density maps from structures to derive a new biasing potential that is smooth, long-ranged, and provides fewer barriers to refinement than established potentials based on cross-correlation [11,13] or inner product (equivalent to using potentials proportional to inverted cryo-EM density [12]). This provides a number of advantages, including an ability to overcome density barriers and in particular avoid excessive biasing forces resulting from large local gradients in cryo-EM density maps.
It also avoids the need for constraints e.g. on secondary structure and rather allows the simulation to freely explore local conformational space, while the experimental data is used to bias sampling to experimentally favored regions of the global conformational space.
We further demonstrate how better balancing between the force field and cryo-EM density components can be achieved by adaptive force scaling derived from thermodynamic principles. This allows refinement with a single fixed parameter at low computational cost for a range of system sizes and initial model qualities. Additionally, to evaluate biasing potentials based on model to cryo-EM density comparison, we suggest a transformation of all-atom structure to model density that reflects cryo-EM specific characteristics while minimizing computational effort.
We investigate advantages and drawbacks of the newly derived potential in practical applications when compared to established inner-product and cross-correlation biasing potentials in a noise-free and experimental cryo-EM data. Finally, we show how the proposed refinement methods rectify a distorted initial model with cryo-EM data. A full open-source implementation is freely available as part of the GROMACS molecular dynamics simulation code [20].

A. Deriving refinement forces via Bayes approach
Our canonical algorithm to refine all-atom models into a cryo-EM density map ρ with molecular dynamics is based on initially roughly aligning density map and target structure, generating a model density from coordinates x, and then determining a dimensionless similarity score S between the generated model density and the target cryo-EM density (Fig. 1).
The similarity is used to derive fitting forces, which are then combined with the force field potential U ff based on a heuristic balance between the density-derived forces and force field determined by a force constant k. The combined driving forces are determined by the total potential energy, (1) We find that applied density forces imply a similarity score between the model structure and target density, and vice-versa. Assuming that a single configuration of atoms gives rise to the observed cryo-EM density, Bayes' approach quantifies the probability density that the given model describes the cryo-EM data as [17] p( x|ρ) ∝ p( x)p(ρ| x) . ( Boltzmann inversion at temperature T connects the left-hand sides of Eqs. (1) and (2) [21] where c is an arbitrary potential energy offset In this formalism, the force field provides the prior p( x) that would have determined the model without any additional observations, while the similarity measure provides the conditional probability that a particular given structure yields a target density p(ρ| x), scaled with force-constant k.  Note that this does not assume any particular form of the similarity, but it provides a general relationship that also relates established similarity measures like cross-correlation or inner-product to their implicit assumptions about the likelihood function above, and in turn enables the construction of new similarity scores that drive refinement procedures depending on the assumptions about the underlying measurement process.
B. Bayes approaches yield negative relative entropy as similarity score To derive a new refinement potential from the likelihood of measuring a density given the structure, we assume that cryo-EM densities are linked to atom-electron scattering probabilities, where electron-atom interaction leads to a phase shift in electrons. With this assumption, two steps are necessary to calculate the density likelihood p(ρ, x) from coordinates. First, an electron-scattering probability density ρ s is created from a given structure.
Second, this density is compared to a given measured density. Cryo-EM micrographs are typically normalized to unity variance around the particle region [22]. As a consequence, cryo-EM densities are scaled by a free parameter, and they may thus be rescaled asρ v = rρ v .
Including this re-scaling in our model, the likelihood to observe a density, given coordinates x is p(ρ| x) = r=0 ∞p(r)p(ρ|ρ s ( x), r). Two further assumptions enable the derivation of new similarity scores.
First, we assume that measured scattering probabilities per voxel are independent of other voxel values. This does not exclude spatial correlation between density data but states that the scattering process in one voxel does not influence the electron interaction in other voxels.
With this result, it suffices to define a probability distribution at each voxel ρ s v . For the per-voxel scattering probability, we present two different sets of assumptions, leading each to a refinement potential in their own right. Many more assumptions may be laid out; here we choose to present those that go beyond previous modeling, yet are well-treatable in model complexity and integratable into a molecular dynamics framework where forces have to be calculated numerically stable and fast. Therefore, we choose not to integrate the full image formation from the microscope detector to three-dimensional density but rather start the modeling process with a three-dimensional density and some additional assumptions on what each voxel represents. By presenting two different approaches to what voxel values in a density present, it should become even clearer how to adapt our modeling to more complex models.
Using the assumption that an incident electron will scatter at one voxel and that the probability at which voxel to scatter is proportional to the density values, the scaling factor r for the density is fixed to r = v ρ. With a scaling factor s that describes the expected number of interactions in the simulated density, this turns into a Dirichlet distribution to describe the scattering process as From this follows In this picture the reported density is treated as a probability density, requiring the removal of negative values and normalization to unity. The resulting potential of this modeling approach is proportional to the Kullback-Leibler divergence between simulated and experimental density with a free scaling parameter. This potential in turn can be seen as an inner-product based potential where density is replaced by its logarithm.
In an alternative picture, reported cryo-EM densities at each voxel represent interaction counts. They are assumed to be proportional to the number of interactions of N incident electrons that each interact with probability ρ s v . This assumes that vitreous ice is not visible and does not contribute to the scattering, which is commonly achieved by shifting the offset of cryo-EM densities so that water density is represented with voxel values that fluctuate around zero. Only accounting for positive density, we describe this scattering interaction process by a Poisson distribution with parameter λ = N ρ s v . While it is theoretically possible to expand the model to include noise fluctuations and negative densities, we omit this for the sake of reducing model complexity.
With these assumptions (the detailed algebraic transformations are laid out in the Supporting Information S1 Appendix), we obtain a similarity score between simulated model density and cryo-EM density proportional to the negative relative entropy, or Kullback- To distinguish both relative-entropy based potentials, we call the potential derived from the first set of assumptions "swapped relative-entropy" and keep referring to the latter as relative entropy.
In a full Bayesian picture the unknown scaling factor r would be treated as a model parameter whose marginal distribution is integrated out. However, the prior distribution for this parameter is unknown and hard to estimate correctly due to the complex density reconstruction algorithm. Even if we were to determine this scaling correctly, a "correct" force scale in such a Bayesian setting would most likely result in very large forces for conformations with significant deviants at the start of simulations, resulting in numerical instabilities in a molecular dynamics setting. The unknown prior distribution p(r) results in an unknown scaling in the force constant, which can be seen from so-called re-gauging with an additional Here, we choose a scaling that ensures that the cryo-EM density has the properties of a probability distribution and choose our arbitrary re-gauging such that vρ v = 1. The force constant balancing cryo-EM data vs. force field / stereo-geometry is thus a free parameter and is chosen adaptively with a protocol described below.
The newly derived relative-entropy-based similarity score has a domain of [−∞, 0] with perfect agreement at zero. Due to the log ρ s v term, it differs prominently from established similarity scores like cross-correlation [11] and inner-product (formulated as a force following the gradient of a smoothed inverted density which is equivalent in this approach [12]; see Supporting Information S1 Appendix). In contrast, the relative-entropy based score receives the largest contribution from voxels where cryo-EM data has no corresponding model density data.
This leads to a different behavior from established similarity scores with local minima for locally good agreement with cryo-EM data while the relative-entropy based potential will only have minima where there is good global agreement between structure and density. As a consequence, the relative-entropy based density potential is expected to perform better in situations where other potentials cannot escape local minima, at the cost of higher sensitivity to noise in the data, and especially additional density data that is not accounted for in the atomic model.

C. The potential energy landscapes based on relative entropy are smooth
The proposed relative entropy density-to-density similarity measure has favorable properties in one-dimensional model refinement of one and two particles to a reference density (Fig. 2). Both newly derived potentials have the same analytical form for our onedimensional model case.
In contrast to cross-correlation and inner product similarity measures that have a steep and sudden onset for refinement forces in one dimension, the relative-entropy similarity score has a harmonic shape with long-ranged interactions that allow for efficient minimization.
Using relative-entropy, the particle to be refined is attracted by a harmonic spring-like potential to the best-fitting position; far away from the minimum forces are large, but their magnitude decreases monotonically as the minimum is approached. Inner-product and crosscorrelation based fitting potentials, however, exert almost no force on the particle outside the Gaussian spread width, while exerting a suddenly increasing force when moving closer to the Gaussian center, and are only insignificant very close to the minimum.
For the refinement of two particles, this advantage is only maintained for one of the newly derived potentials, where the relative-entropy-based potential energy landscape is less rugged and has fewer pronounced features and minima than the corresponding landscapes for the inner-product and cross-correlation based potentials (Fig. 2). Only a single diagonal barrier is found in the relative-entropy-based potential landscape, corresponding to a "swapping" of particle positions, which alleviates the search for a global minimum. The inner-productbased free energy landscape has its minimum at a configuration where both particles are at the same position at the highest density. This issue can be alleviated in practical applications through a force-field prior that would enforce a minimum distance between the atoms (e.g. through van der Waals interactions). The swapped relative entropy potential on the other hand exhibits behavior similar to the inner-product, with similar minima but an overall smoother energy landscape.
To model the influence of a force field, the two particles were connected with a harmonic bond with increasing influence. The balance between density-based forces and bond strongly determines the shape of the resulting energy landscape, but here too relative entropy provides Figure 2. Similarity score determines ruggedness of the effective refinement potential energy landscape, also when balancing it with structural bias. From top to bottom: a One-dimensional refinement of a single particle (black circle) towards a Gaussian-shaped density (gray) with innerproduct (purple), cross-correlation (ochre), relative-entropy swapped (dark blue) and relativeentropy (green) as similarity scores. b Expanded model with two particles (black circles, x 1 smaller and x 2 larger) with two amplitude peaks in a one-dimensional density and target distribution (gray), and the resulting two-dimensional effective potential energy landscapes for inner-product, crosscorrelation, swapped relative-entropy and relative-entropy similarity measures. c Combination of the similarity measure and force field contribution to the potential energy landscape, exemplified by a harmonic bond that keeps particles at half the distance between the Gaussian centers. For all relative weights of the contributions of the refinement potential and bond potential energy landscape (ratio 1:2 upper panel, 2:1 middle panel, as illustrated by the scale on the left), the relative entropy similarity score produces smooth landscapes with minima at the positions that are expected from the model input.
a smoother landscape less sensitive to the specific relative weight of refinement and bond potentials.
D. Adaptive force scaling reduces work exerted during refinement and allows for comparison of density-based potentials For now, we cannot derive the force constant for density-guided simulations automatically from the cryo-EM density alone and it is thus set heuristically. Established protocols where the force constant has to be determined manually require an iterative trial-and-error approach. We address this by introducing an adaptive force-scaling as depicted in Fig. 3a to automatically balance force-field and density-based forces during the refinement.
Cryo-EM refinement simulations are non-equilibrium simulations with the aim to drive a system from an initial model state to a final state that is as similar to the cryo-EM density as possible while avoiding structural distortions that result e.g. from unphysical paths. To avoid or at least reduce the latter during refinement, a heuristic protocol has been devised that aims to minimize work exerted on the system while still requiring as little time as possible for the refinement. To minimize the exerted work formulated as the adaptive-force scaling starts from a low force constant k. This is then increased if similarity decreases and conversely decreased if the similarity is increasing (Fig. 3b). Any feedback protocol of this type is guaranteed to exert less work to reach the same similarity score than keeping the force constant fixed at the final value of the adaptive scaling protocol, given that the score is monotonically increasing.
A one-dimensional Brownian diffusion model system (Fig. 3c) is used to test the performance of the concrete scaling protocol as described in the methods section of this paper. In this model, the similarity score simply increases with increasing particle coordinate value.
Biasing the system towards increasing coordinate values with adaptive force scaling in contrast to a constant force allows for the particle to reach a given coordinate value in the same average first-passage time at a much lower average work input. Without any coupling of the free energy landscape to the adaptive force scaling protocol other than through the particle  Table 6 in Supporting Information S1 Appendix). For convenience, we approximate scattering amplitudes by unity for all heavy atoms and zero for hydrogens. The magnitude of thermal fluctuation of atoms at cryogenic temperatures determines the spread width σ.
In practice, these limitations to the model resolution are superseded by the finite performance of the measurement instrument and the reconstruction process where structural heterogeneity, detector pixel size, microscope lenses, and particle alignment limit the resolution. We do not account for structural heterogeneity, because it is an ensemble effect.
A connection between the approach presented here and an ensemble model may be made though by employing a probability distribution p( x) instead of x in Eq.
(2) and leveraging ensemble simulations [23]. Other resolution-limiting effects are approximated by additional convolution of the generated maps with a Gaussian kernel. Rather than aiming to reproduce the same blur as in the experimental map, we strive to preserve as much information as possible from the physical model.
A balance between information loss due to under-sampling on the grid on the one hand and information loss due to coarse blurring is found where the Gaussian width at half maximum height equals the resolution. The maximum representable resolution on a grid corresponds to twice the Nyqvist frequency δ (corresponding to the pixel and voxel size) so that the Gaussian width σ is approximated in refinement simulations from the highest local resolution or, where that data is not applicable, from the voxel-size, For computational efficiency Gaussian spreading is truncated at 4σ for all simulations in this publication, accounting for more than 99.8% of the density (Fig. 1 in Supporting Information S1 Appendix). The small limitation on the maximal distance between the initial model structure and the cryo-EM density through this cutoff has proven to be irrelevant for all practical purposes, as density-based forces will "pull" structures into densities as soon as there is minimal overlap between model density and cryo-EM density, which can easily be achieved with an initial alignment. Interestingly, this approach results in a smaller Gaussian spreading width than previously applied ones that aim to reproduce a density map with the same overall resolution as the experimental cryo-EM density. As a result, it maintains as much structural information as possible in the model density while still reducing the computational costs.

F. Refinement against noise-free data
To separate additional noise effects in experimental data and possible limitations in the above model, we first assess refinement with ideal data where a small straight helix model system [14] has been refined against a synthetically generated target density of the same helix in a kinked configuration. As illustrated in Fig. 4a, adaptive force scaling and relativeentropy as similarity score efficiently fit the helix into the synthetic cryo-EM density [24].
The combination of adaptive force scaling and different similarity scores achieved a consistent global fit when the helix was aligned to the density, with some fluctuations of the results (Fig. 3c) due to the stochastic nature of molecular dynamics simulations. Simulations starting from both the aligned and unaligned starting positions get stuck in local minima, which can result in bad fits. The average total RMSD of all replicates was lowest for relative-entropy starting from the aligned position and further improved when the helix was initially unaligned (Tab. 1 and Tab. 2 in Supporting Information S1 Appendix). The relative-entropy based potential shows markedly better results for the unaligned refinement and achieved a fit with less 1Å RMSD in 6 out of 7 replicates while inner-product, cross- correlation, and relative-entropy swapped based potentials in some instances completely fail to align the helix (Fig. 5 in Supporting Information S1 Appendix). For a single helix this is a slightly artificial case, but in a large structure undergoing significant transitions, it will be common for some secondary structure elements to not overlap with the target density in an initial phase of the refinement.

G. Refinement against experimental cryo-EM data
Experimental cryo-EM densities of aldolase and a GroEl subunit were used to test the performance of adaptive-force scaling in combination with different refinement potentials on experimental data with increasing amounts of density that is not accounted for in our model description and noise that cannot be fully accounted for by the model assumptions.
By using adaptive-force scaling refinement of a previously published X-Ray structure of rabbit-muscle-aldolase [25] against a recently published independently determined cryo-EM structure [26], we consistently achieve accurate refinement throughout all potentials with good stereochemistry (Tables 3 & 5 in Supporting Information S1 Appendix). Figure 5a shows the final models of refinement with a global deviation of less than 1Å heavy-atom root mean square deviation (RMSD) from the deposited model using inner-product and crosscorrelation measures and slightly above 1Å for the relative-entropy based density potential ( Table 3 in Supporting Information S1 Appendix).
The close agreement with the cryo-EM data is reflected in the FSC of the models refined against the density (Fig. 5b) being nearly indistinguishable from the deposited model. The relative-entropy-based potential emphasizes agreement with global features at the cost of local resolution (Fig. 5 in Supporting Information S1 Appendix), while still providing good agreement to the cryo-EM density.
The unweighted FSC average [27] serves as an established similarity score that is not related to the biasing potentials which were used to refine the system (Fig. 5b). All underlying potentials appear to lead to refinement simulations that converge in less than 2 ns, as shown in Fig. 5c. The less rugged and long-range potential properties of relative-entropy based density forces are reflected in a rapid rigid-body like initial fit to global structural features, while the other potentials show gradual improvements in fit.
The refinement of a GroEl subunit in two different conformations as determined by cryo-EM [28] stretches the limits of the model assumptions of our refinement potential by refining it against a more noisy model with imperfect map-to-model correspondence. Similar to aldolase refinement, adaptive force-scaling allows for rapid and reliable refinement into the model density, as shown in Fig. 6. However, the relative-entropy based potential propensity to taking all density into account leads to deviations from the published model in regions with density that has no correspondence in the all-atom model ( Fig. 9 in Supporting Information S1 Appendix).

H. Model rectification by combining force-field and cryo-EM data
To assess performance in larger structural transitions, we repeated the aldolase refinement when starting from initial model structures that have been distorted by heating with partially unfolded secondary structure elements (Fig. 7, as described in Methods). Figure 7b shows the final relative-entropy based model of the refinement procedure that achieved 1.13Å heavy-atom RMSD from the manually built model. Structural details at map resolution match in secondary structure elements. In contrast to refinement of the undistorted X-ray structure, the relative-entropy based potential gains less from the long-rangedness of the potential and the rapid alignment of large-scale features, because structural rearrangements were needed on all length scales. The adaptive force scaling protocol alleviates differences between density-based potentials in refinement speed and allows for refinement with good structural agreement in less than 3 ns.
The adaptive force scaling protocol allows the modeling to be more steered by cryo-EM data and reach model structures that would not have been accessible by modeling using stereochemical information from the force-field alone.

III. DISCUSSION
While defining a purely empirical similarity measure can sometimes suffice to fit structures to cryo-EM densities, connecting the similarity measure to the underlying measurement process of the target density enables derivation of natural similarity measures. From very few assumptions, this results in the density-based potential derived via Bayes' approach that coincides with the relative-entropy between a model density generated from model/simulation atom coordinates and the target cryo-EM density.
The newly defined potentials have favorable features, with the Poisson statistics-based relative entropy potential most prominently exhibiting long range and low ruggedness. These avoid local minima that do not correspond to desired configurations during refinement and allow rapid alignment of large-scale features, and they perform superior to established refinement in the zero noise setting with synthetic density maps. The noise content in current cryo-EM densities is likely still too high to be handled with the current minimalistic model assumptions, but as the quality of cryo-EM and other low-to-medium resolution techniques continues to rapidly improve, we believe there will be even more advantages to models that do not depend on smoothing. In addition, the adaptive force-scaling provides a surprisingly simple way to tackle the inaccessibility of the balance between force-field and density-based forces within our model assumptions. It allows for parameter-free refinement that is one to two orders of magnitude faster than currently established protocols. Conceptually it is orthogonal to, and easily combined with, multi-resolution protocols [15].
Another illustration of the usefulness of the Bayesian framework for handling force-field vs. fitting forces is how it enables us to deduce a close-to-optimal model density spread for refinement, and even more so that this value is not identical to the common practice of setting it equal to the experimental resolution. While many of these factors could still be tuned manually, removing them as free parameters means fewer arbitrary settings that avoid over-or underfitting, which will be even more important when trying to combine e.g. multiple sources of experimental data. For trial structure refinement against recent cryo-EM data, we show that we achieve excellent fits independent of initial model quality.
One limitation of the current formulation is that it does not explicitly take more information from the cryo-EM density reconstruction process into account. A first step to broaden the approach presented here is to account for the local resolution information from the cryo-EM map, which may be seen as a measure of the noisiness of the density data. In practice, the local resolution will still influence the fit since low local resolution will correspond to smoother regions of the map, and lower-magnitude gradients will lead to lower-magnitude fitting forces in those regions. However, a formally more correct way to address the problem is likely to treat the target cryo-EM density as a statistical distribution with a variance that is spatially resolved -this is something we intend to pursue in the future to see whether it can further improve the issue of the relative-entropy potential with more noisy data.
The algorithms proposed in this work are freely available, integrated, and maintained as part of GROMACS [29]. Overall, three independent building blocks are provided to aid the modeling of cryo-EM data that each may be individually implemented in current modeling tools: A new refinement potential, a new criterion for how to calculate the model density, both based on reasoning via Bayes' approach, and adaptive force scaling to gently and automatically bias stereochemistry and cryo-em data influence. The implementation also provides tools to monitor the refinement process. Although it can still be difficult for any automated method to compete with manual model building by an experienced structural biologist, we believe these methods provide new ways to extract as much structural information as possible from cryo-EM densities at minimal human and computational cost, which is particularly attractive e.g. for fully automated model building.

A. Calculating density-based forces
For ease of implementation and computational efficiency the derivative of Eq. (4) is decomposed into a similarity measure derivative and a simulated density model derivative, summed over all density voxels v Though the convolution in eq. (11) might be evaluated with possible performance benefits in Fourier space, we have implemented the more straightforward real-space approach.
The forward model ρ s is calculated using fast Gaussian spreading as used in [30]; the integral over the three-dimensional Gaussian function over a voxel is approximated by its function value at the voxel center v at little information loss (Fig. 2 in Supporting Information S1 Appendix). Amplitudes of the Gaussian functions [31] have been approximated with unity for all atoms except hydrogen. The explicit terms that follow for S(ρ, ρ s ) and ∇ r ρ s v (r) are stated in the Supporting Information S1 Appendix.

B. Multiple time-stepping for density-based forces
For computational efficiency, density-based forces are applied only every N fit steps. The applied force is scaled by N fit to approximate the same effective Hamiltonian as when applying the forces every step while maintaining time-reversibility and energy conservation [32,33].
The maximal time step should not exceed the fastest oscillation period of any atom within the map potential divided by π. This oscillation period depends on the choice of reference density, the similarity measure, and the force constant and has thus been estimated heuristically.

C. Adaptive force scaling
Adaptive force constant scaling decreases the force-constant when similarity increases by a factor 1 + α, with α > 0, and reversely increases it by a factor 1 + 2α when similarity decreases. The larger increase than decrease factor enforces an increase in similarity over time.
To avoid spurious fast changes in force-constant, similarity decrease and increase are

E. Comparing refined structures to cryo-EM densities
Fourier shell correlation curves and un-weighted Fourier shell correlation averages [27] were calculated at 4 ps intervals from structures during the trajectories by generating densities from the model structures using a Gaussian σ of 0.45Å, corresponding to a resolution as defined in EMAN2 [24] to 2Å.
F. Map and model preparation before refinement

Helix
Noise-free helix density maps at 2Å simulated resolution on a 1Å voxel grid were generated from an atomic model using "molmap" as provided by chimera [34]. Two frames taken from an equilibration simulation of the helix model were used as starting models for subsequent density-guided refinement.

Aldolase
A simulation box of the exact dimensions as aldolase density map EMD-21023 was used. The corresponding aldolase model (PDB id 6V20) was treated as a ground truth for RMSD calculations. A previously determined X-ray model (PDB id 6ALD) was used as starting structure. The system was subjected to energy minimization before fitting. No symmetry constraints were used in simulations.

GroEL
The density-guided simulations of GroEL were performed between conformational states within an individual oligomer [28]. State 1 (PDB id 5W0S(1)) was used as starting model and fit toward conformational state 3 (EMD-8750, additional map 3). A section of the density map corresponding to a single oligomer was used as the target. The target map was created using ChimeraX [35] and by including any volume data within a 5Å range of the previously determined model corresponding to state 3 (PDB id 5W0S (3)). The target density was centered in a map and the map size was set to 100³ using relion image handler [36].

G. Generation of a distorted model
To generate a distorted starting model, the aldolase protein X-ray was heated to 433 1 3 K over a period of 5 ns. During heating, the pressure was controlled with the Berendsen barostat, favoring simulation stability over thermodynamic considerations. To disentangle effects from decreasing the temperature and fitting, the distorted structure was subjected to 5 ns of equilibration at 300 K before starting the density-guided simulations with the same protocol as described above.

H. Molecular dynamics simulation
All simulations were carried out with GROMACS version 2021.3 [29] and the CHARMM27 force-field [37,38] in an NPT and NVT ensembles with neutralized all-atom systems in 150 mM NaCl solution. The temperature was regulated with the velocity-rescaling thermostat at a coupling frequency of 0.2 ps to ensure rapid dissipation of excess energy from density-based potential, when structures are very dissimilar from the cryo-EM density, i.e., far from equilibrium. The pressure was controlled with the Parinello-Rahman barostat for aldolase simulations, helix simulations were carried out at constant volume. Aldolase and GroEL were aligned roughly to the density by placing their center of geometry in the center of the cryo-EM density box. Forces from density-guided simulations were applied every N fit = 10 steps according to the protocol described above [33]. All simulations to refine a structure against a density were carried out with adaptive force-scaling. We used a coupling constant of τ = 4 ps, balancing time to result with time for structures to relax. For aldolase simulations, the Gaussian spread width was determined by using a lower bound on the highest estimated local resolution of 1.83Å. Spread width for GroEl simulations was set to 1.0455Å, based on the 1.23Å voxel size in the map used for refinement. Periodic boundary conditions are treated as described in the Supporting Information S1 Appendix.
All simulation setup parameters and workflows have been made available.

V. DATA AVAILABILITY
Simulation starting structures, generated densities, setup parameters, and complete workflow setups via Makefiles and Python scripts to generate Fig. 2 as well as data for Fig. 3 and per-residue RMSD are available via Zenodo (https://doi.org/10.5281/zenodo.4556616).
The code to perform density-guided molecular dynamics simulations is maintained within GROMACS and publicly available in release 2021 and later, as well as in the repository at https://http://gitlab.com/gromacs/gromacs. Fourier shell correlation analysis of trajectories has been implemented on top of the GROMACS codebase following conventions in EMAN2 [24] and is available at https://gitlab.com/gromacs/gromacs/-/commits/fscavg. Python scripts to generate Fig. 2, as well as data for Fig. 3 and per-residue RMSD are available via Zenodo (https://doi.org/10.5281/zenodo.4556616).