Modelling dynamics in protein crystal structures by ensemble refinement

1Since the dawn of structural biology there have been experimental observations of dynamic motion in proteins and other biomolecules. Multiple biophysical methods have firmly established that such atomic ‘wigglings and jigglings’ play an inherent role in both protein structure and function; and, in conjunction with high-resolution structures insights into dynamics aid the understanding of biomolecular functions in catalysis, ligand or drug binding and macromolecular interactions. Presently X-ray diffraction and NMR spectroscopy are the primary source of data for high-resolution protein structures. Whereas microscopy methods may provide information regarding long-range conformational changes, NMR characterizes fluctuations at atomic detail. However, due to the challenging nature of such experiments the number of dynamics studies is relatively sparse in contrast with the wealth of structural information available in the Protein Data Bank (PDB). The majority of entries in the PDB derived from X-ray diffraction data are presented as static, single, structures, although there is often extensive disorder resulting from protein dynamics and crystal-lattice distortions. Extracting atomic fluctuations from these diffraction data would dramatically increase the scope for dynamics studies of biomolecules and potentially reveal atomic details of structure–function–dynamic mechanisms that have previously been obscured.

2The diffraction data of proteins are affected by multiple sources of disorder, notably arising from atomic vibrations, concerted motions of protein domains and inter-molecular lattice distortions. Structural models of proteins should account for both anisotropic and anharmonic distributions around the mean atomic positions to reproduce the observed Bragg intensities accurately. However, explicit modelling of such distributions in macromolecules using current methods requires extensive parameterization inappropriate for the diffraction quality of a typical protein crystal. Multi-conformer structures represent both anisotropic and anharmonic disorder, but despite numerous attempts at automating the inclusion of minor conformations, 95% of all protein residues in the Protein Data Bank (PDB) derived from diffraction data are modelled with a single conformation. As opposed to multiple discrete models, a MD simulation with time-averaged restraints results in a population of structures in which the individual models are interrelated by a Boltzmann-weighted energy function. This method introduced by and implemented in macromolecular crystallography by, showed a reduction in R-value. However, cross-validation introduced subsequently revealed chronically over-fitted models.

3Here, we present an ensemble-refinement method that restricts the number of structures modelled and thereby prevents over-fitting of the data. We model large-scale motions, attributable to, for example, lattice distortions, by an underlying global disorder model. This approach allows MD simulations to sample local atomic fluctuations only, without the need for sampling large-scale global disorder. We show that the method yields reproducible ensembles with improved fit to the X-ray data, as validated by cross validation, Rfree, and stereochemical analyses. Analyses of the ensembles show that detailed features are observed indicating atomic fluctuations that may be relevant for the biological function of the macromolecules.

4We performed MD simulations, in which the model was restrained by a time-averaged X-ray, maximum-likelihood target function (see ‘Materials and methods’). The X-ray restraint optimized 〈Fcalc(hkl)〉 against Fobs(hkl), where 〈Fcalc(hkl)〉 are computed as rolling averages from the structures in the MD trajectory, with the length of the averaging window determined by the relaxation time τx. This approach contrasts with the traditional crystallographic refinement approach, where Fcalc(hkl) are computed from a single structure and optimized against Fobs(hkl).

5Prior to the simulations we approximated the large-scale disorder by an overall TLS model derived from the atomic B-factors of the refined single structure. Using one TLS group per protein molecule or domain, we iteratively fitted TLS parameters to the atomic B-factors of the protein atoms excluding atoms with large deviations in B-factor from the TLS-derived B-factor (the parameter pTLS described the percentage of atoms included in TLS-fitting; see ‘Materials and methods’). The resulting TLS model was applied to all atoms throughout the simulation. Effectively, this TLS model of the protein core excludes the effects of hyper-flexible surface loops and, hence, describes the global disorder that may be attributed to inter-molecular lattice distortions and overall intra-molecular breathing or domain shifts.

6Ensemble refinement was tested using 20 diffraction datasets from the PDB and started from either the PDB or PDB_REDO structures (‘Materials and methods’). Upper resolution limits of the datasets ranged from 1.1 to 3.1 Å resolution and structures had 50 to 1,004 amino-acid residues in the asymmetric unit (Table 1). The simulations were run at an effective temperature of 300 K for the protein atoms, using a temperature bath (Tbath) slightly below 300 K to allow for heating due to the non-conservative nature of the time-averaged X-ray restraint modulated by its weight wx-ray (‘Materials and methods’). Explicitly modelled solvent atoms were added and/or removed intermittently during the simulation dependent on the corresponding electron-density and difference maps (‘Materials and methods’). Bulk solvent effects were accounted for by an averaged Flat Bulk-Solvent Model (‘Materials and methods’). The parameters pTLS, τx and the Tbath and wx-ray pair were optimized in a grid search resulting in a shallow optimum scored by Rfree (Figure 1A). After a period of equilibration, the trajectory of structures was acquired over an extensive period of time (40 times τx).

7Figure 1B shows the R-values as they developed over the simulation time for a structure with PDB code 1UOY, for which the largest improvement in Rfree was observed among the datasets tested (possibly due to the high degree of anisotropic and anharmonic side-chain motion for this case). The R-values started at a high value and remained high (∼35%) for the individual structures, which is in agreement with the observation that the derived global TLS B-factor model is not optimal for fitting a single structure to the data. Averaging the structure factors over the relaxation time τx of 1 ps (corresponding to the rolling average structure factors used in the X-ray restraint) dropped the Rwork and Rfree to ∼11% and ∼15% respectively. The Rwork and Rfree of the collected ensemble of structures (corresponding to unweighted averaged structure factors) monotonically decreased to 10.3 and 13.7% respectively. over the acquisition period of 40 ps. The improvement in R-values from the ensemble model with respect to the single-structure model spanned the entire resolution range of the data (Figure 1C). Acquisition over 40 times τx yielded a highly redundant set of structures. We reduced the number of structures by calculating the minimum number of structures, that is 167 in the case of 1UOY, required to reproduce the R-value of the trajectory (Figure 1D and ‘Materials and methods’).

8Analysis of all 20 datasets showed that ensemble refinement improved the Rfree by between 0.3 and 4.9 percentage points compared to single structures re-refined using the same program package, that is Phenix, with a mean improvement of 1.8% in Rfree values (Table 1, Figure 2A). The effect of the starting structure on ensemble refinement was assessed by using alternative refinement programs, phenix.refine, Refmac, and Buster, to generate varying input models. No significant differences were observed due to the different starting models (Tables 3 and 4). The improvement in Rfree, number of structures in the final ensemble and the averaging time τx tended to increase with resolution (Figure 2A–C). The optimum values for the parameters pTLS and Tbath are not correlated with resolution (Figure 2D,E). Concomitant with the reduction in R-values, the ensemble models reduced electron-density differences, decreasing rms fluctuations in difference maps by 0 to 41% with an average of 12% improvement (Table 2). The difference electron-density maps for the single-structure and ensemble models indicated improvements throughout the asymmetric unit cell, as exemplified in Figure 1E.

9We used the high-quality experimental phases available to high resolution for 1YTT of mannose-binding protein to validate the ensemble-refinement method. The overall correlation coefficient between the electron-density map from the ensemble model (obtained without experimental phases) and the experimentally phased electron-density map was 0.903, compared with 0.873 and 0.895 for the published and re-refined single structures. These seemingly small improvements in overall quality indicators allow for significant local improvements. Real-space correlation coefficients highlighted marked local improvements for flexible residues in particular (Figure 3A) with 11 residues improving by more than 0.1 in correlation coefficient. This observation was consistent with local improvements in electron-density differences in regions of flexible or disordered side chains (Figure 3B). Moreover, the ensemble model contained structural details previously identified in a multiple-model approach by, as shown for the anisotropic distribution for the side chain of Phe121 (Figure 3C) and diffuse water shells around hydrophobic residues (Figure 3D). Figure 3E shows that the most occupied water sites in the ensemble correlated with low atomic B-factors for waters in the single-structure model.

10Next, we analysed the stereochemistry of the computed ensemble models. The robustness of the observed atomic distributions was tested by repeating ensemble refinements 10 times using different random starting velocities. Figure 4A shows that the observed distributions are highly reproducible. With data extending to 1.5-Å resolution correlations above 0.99 were observed between residue distributions from separate runs. At lower resolutions, the majority of residues showed correlations above 0.95, with occasionally correlations dropping below 0.8 in very flexible regions (see Figure 4B,C). Clearly, the ensembles contain multiple values for each geometrical term that form a distribution, instead of a single stereochemical value obtained from a single structure; Figure 5A–D presents examples of side-chain distributions (by χ1 and χ2 angles) observed in the ensembles along with standard deviations computed from the 10 repeats. Averaged over all 20 cases, the rms deviations from idealized bond lengths, bond angles and dihedral angles for the re-refined single-structures were 0.010 Å, 1.26° and 15.2° respectively. (Figure 6—source data 1). These deviations decreased for the ensemble models by 0.002 Å, 0.26° and 6.6° respectively, when considering the centroids of the observed stereochemical distributions. Taking all fluctuations around the centroids (i.e. complete distributions) into account, these values increased by 0.002 Å, 0.33° and 4.0° respectively compared to the statistics from single structures. This indicates high stereochemical quality for the ensemble model, but that the ensemble of structures contained fluctuations exhibiting larger deviations from ideality. Figure 6A,B shows that high-energy conformations, as indicated by for example non-favourable Ramachandran φ,ψ-angle combinations, occurred transiently and were concentrated in regions of structural flexibility. Counting the most frequent Ramachandran classification for each φ,ψ-angle showed that the ensembles have a similar percentage for ‘allowed’ and an increased number of ‘outliers’ compared to the single structures (Figure 6C and Figure 6—source data 2). These analyses illustrate that in ensemble refinement conformations were sampled, rather than optimized to a single configuration as in single-structure refinement. Similar to Brünger, we observe that lower Rfree values correlate with better quality of the Ramachandran statistics (Figure 6D).

11The presence of non-crystallographic symmetry (NCS) allowed for crystallographically-independent observations of atomic fluctuations in multiple copies of a protein molecule (Figure 7A–C). In some cases, the applied global TLS models differed significantly between NCS-related copies (Figure 7B). Nevertheless, we observed atomic fluctuations similar both in magnitude and location for related copies in areas not affected by crystal packing (Figure 7B,C; additional cases of NCS are presented in Figure 7—figure supplement 1–4). Apparently, variations in overall disorder arising from packing differences of NCS copies (as indicated by different B-factor distributions) were well accounted for by the applied global TLS models. Similarly, a global increase in disorder present in a dataset collected at ambient temperature vs an isomorphous dataset collected under cryo-conditions was fully accounted for by an increase in global TLS (Figure 8A). These data indicate that the derived atomic fluctuations are molecular traits and that the global TLS model accounts for overall disorder, which includes for example lattice or packing effects.

12Inspection of the obtained ensembles showed that most proteins, as expected, are characterized by well-ordered residues in the protein core and flexible residue side chains and loops on the outside (an example is given in Figure 9A). However, three cases exhibited marked flexibility of residue side chains on the inside of the molecule. 1BV1, major birch pollen allergen, has a large forked solvent channel with multiple disordered side chains and a diffuse water network (Figure 9B). The cavity is consistent with its putative role as a general plant steroid carrier. Presumably, the flexible internal residues play a role in binding the diverse ligands. More surprising are the disordered cores in 1X6P of PAK pilin (Figure 9C) and 3K0N of the enzyme proline isomerase (Figure 9D); in both these cases, the datasets were recorded at ambient temperatures. Multiple (16) aliphatic and aromatic side chains are highly flexible, forming a molten core in the pilin molecule. These flexible residues, which are extremely well conserved in the type IV pilin family, create the central interface between the characteristic long α-helix and β-sheet of this protein fold. We hypothesize that this monomeric pilin structure represents an intermediate molten state, which becomes stabilized upon protomer filament formation. The third case with flexible residues on the interior is 3K0N of proline isomerase. As with the pilin structure, several (11) aromatic and aliphatic residues showed large side chain fluctuations, yielding a molten core of the protein structure. However, the same protein under cryogenic conditions (3K0M) showed mostly well-ordered side chains in the ensemble (Figure 9D, right-hand side), indicating that at cryogenic conditions the molten core has been annealed to its ground state configuration. As discussed in more detail in the next paragraph, the observed flexibility of the core residues at ambient temperature is likely of functional relevance for the enzyme. Thus, the computed ensemble models highlighted a hitherto unnoticed phenomenon of molten cores in folded proteins, which are likely relevant for the biological function of these molecules.

13NMR spectroscopy has previously revealed specific dynamics for active-site residues of proline isomerase. The solvent-exposed residues Arg55 and Met61 in the active site showed disorder in 3K0M, where data were collected at 100 K. For 3K0N collected at 288 K, a number of additional residues with multiple conformations were observed (Figure 8A,B). These included Ser99, Phe113, which are part of the substrate-binding pocket together with Arg55 and Met61 (Figure 10A), and Leu98, which neighbours the flexible residue Ser99 but points into the hydrophobic core. Ensemble refinement of the 288 K data revealed a large number of residue side chains in the core to be flexible. This flexibility in the core appears to be linked to the dynamics of active-site residues through the intervening β-sheet. In particular, the main-chain H-bonding network (C=O·HN) of neighbouring β-strands within the sheet was flexible, as indicated by anisotropy in the C=O bonds of residues 55-62-113-98 (with the largest anisotropy observed for 55 and 62; see Figure 10B). Analysis of the side-chain conformations for the active-site residues Arg55, Met61, Ser99 and Phe113 showed a 10% minor conformation for the four active-site residues (Figure 10C), which is in good agreement with NMR relaxation data (see Figure 2 in ). Mutation of Ser99 to Thr (>14 Å away from the catalytic Arg55) affects the side-chain distributions and lowers the activity ∼300-fold, similar to an Arg55Lys mutation of the catalytic residue. Thus, the ensemble refinement results support the notion put forward by Eisenmesser et al. and Fraser et al. that side chain dynamics play a critical role in the enzymatic function of proline isomerase and, moreover, expand upon this theme to reveal mechanistic insights arising from the underlying detailed dynamics.

14Ligand binding to HIV protease is known to have marked effects on the enzyme structure. We compared HIV protease in its apo form, 2PC0, and bound to ligand JE-2147, 1KZK. As for proline isomerase, HIV protease exhibited flexible, moldable, substrate-binding pockets in the apo state. Enthalpic and entropic binding of the ligand with high affinity (KD = 41 pM) reduced the flexibility in the substrate binding pockets by protein–ligand H-bond interactions and van der Waals stacking (Figure 11A,B). Upon ligand binding, Asp90 became ordered through H-bonding with the ligand in P2, whereas its dimeric partner lacked a H-bonding partner in P2′ and remained flexible as in the unbound state. The canonical aspartic protease catalytic residues, Asp25 of both monomers, became ordered upon ligand binding. Concomitantly, we observed significant changes in dynamics of specific core residues (Figure 11C). Some residues, most notably Thr26, ‘froze’ (Thr26 is part of the conserved Asp25-Thr26-Gly-27 sequence). In contrast, the side chains of other residues, most notably Cys95 and Leu97 opposite of Thr26, became highly disordered in the bound state, whereas they were relatively ordered in the unbound state. This observation supports NMR data that showed that conformational variability increases upon inhibitor binding for Leu97 amongst others. These data suggest that the entropy lost by the catalytic aspartates upon ligand binding is compensated with an increase in disorder of specific core residue side chains. This type of dynamic modulation was also observed for Ca2+ binding to calmodulin, where this effect was dubbed entropy compensation. Similar to the molten core dynamics for proline isomerase, the structure ensembles generated by the ensemble refinement method revealed specific core dynamics for HIV protease, in particular a conformational exchange that is likely functionally relevant.

15The development of new small molecule therapeutics is often supported by the use of macromolecular structure, typically X-ray crystallography of complexes between target proteins and drug candidates. These complexes are typically interpreted as static structures, and the impact of dynamics, if considered at all, is probed using computational methods. Our new ensemble refinement approach makes it possible to study the role of dynamics in drug–target complexes in the context of the experimental data. Therefore, we analysed two structures of Abl kinase in complex with Imatinib (also known as Gleevec), that is 1IEP, and PD173955, that is 1M52. These compounds bind the Abl kinase with high affinity, 37 nM and 100 nM respectively. The ensembles provide insights into the flexibility of the protein residues and the ligand moieties in the complex. Figure 12A shows the variation in H-bonding observed in Abl kinase–Imatinib. Variable H-bonding interactions were observed for the hydrophilic N-methylpiperazine moiety with the backbone carbonyl atoms of Ile360 and His361. In contrast, the ensemble displayed a well-ordered H-bond between the anilino-NH and Thr315 ‘gatekeeper’ side-chain. Moreover, the ordered water network between Glu286, Lys271 and the pyrimidine moiety of Imatinib was reproduced in the ensemble model (Figure 12B). We observed that the Abl kinase adopts two different states in these crystal structures. In the Imatinib complex the activation loop, residues 381–402, is highly disordered (Figure 12C), which was confirmed by comparison to previously published NMR data. In general, the ensemble models indicate details of tight and highly ordered drug–target interactions on one side vs disordered interactions elsewhere, which are indicative of less tight interactions, that may suggest which sites to modify in a drug-optimization cycle.

16We have shown that far more structural information can be reliably extracted from protein diffraction data than is achieved to date by traditional single-structure modelling methods. Our ensemble refinement method samples distributions that reflect structural details of protein dynamics. The resulting ensemble models provide a more comprehensive description of the molecules and allow interpretation of the molecular function in terms of both the three-dimensional arrangements of the protein residues and their flexibilities. Moreover, ensemble models minimize the risk of structural over-interpretation associated with the seemingly rigid single-structure models. We found comparative analyses of protein molecules in different states to be very useful for identifying detailed changes in structural dynamics that may be mechanistically relevant for the molecular function.

17Partitioning large-scale disorder into a global model separates intermolecular variations of protein packing in the crystal from the detailed intra-molecular atomic fluctuations. Effectively, the X-ray gradient dictates the MD sampling to yield featureless, (mFobs−DFmodel)exp[iφmodel], electron-density difference maps, while the global disorder model is accounted for by taking BTLS into account when computing the atomic densities. In this way, the ensemble of structures is generated to model the anisotropic and anharmonic electron-density distributions precisely, while being restrained by the bonded and non-bonded energy terms used in the MD simulation. The separation of global disorder and local atomic fluctuations contrasts the original approach by, where the MD sampling had to account for both the large scale global disorder and local fluctuations leading to very long relaxation times τx of 16 ps. In the current work much shorter relaxation times of 0.25–2 ps can be used, thereby limiting potential over-fitting markedly. The method is applicable to data with a wide range of upper resolution limits. We see marked improvements in Rfree for datasets ranging from 1.5 to 2.6-Å upper-resolution limit. A detailed interpretation of the ensembles is allowed, supported by the very high local correlations between independent ensemble refinements. However, at lower resolution limits and for highly disordered loops the local correlation between independent runs drops and detailed interpretation is not feasible. Thus, even though the number of independent parameters in an ensemble model is not clearly defined (and therefore the parameter-to-observation ratio is unclear), the gain in Rfree and the very high local correlations between independent runs indicate a high reliability of the ensemble models. However, the method is not a panacea for highly disordered protein regions. In the absence of ordered conformations for a certain region of the protein (as implicitly defined by the diffraction data) the ensemble refinement will sample diverse conformations in order to prevent the build up of negative peaks in the electron-density difference map. In other words, if the data ‘says’ that a region is disordered, ensemble refinement will generate diverse conformations for that region. Furthermore, dataset pathologies caused, for example, by radiation damage may have confounded effects that obscure the dynamics inherent to the protein molecule. Thus, perhaps somewhat counter-intuitively, this modelling method that accounts for inherent protein dynamics does not help to resolve structural details of disordered regions, but is particularly suited to resolve dynamical fluctuations in ordered parts of the protein structure.

18Ensemble refinement of 20 protein datasets highlighted global dynamics features of protein molecules. Surprisingly, in some cases the ensembles indicated the existence of folded protein structures that display molten cores. Most likely, such molten cores may indicate intermediates of protein molecules that function in larger complexes (such as PAK pilin), or alternatively these molten cores support dynamical fluctuations that are needed for ligand binding and enzyme functioning (as for birch pollen allergen and proline isomerase respectively). Furthermore, the ensembles show details of specific order–disorder transitions, or conformational exchanges, between active site and core residues (as for HIV protease in the unbound and bound state) that suggest a mechanism of entropy compensation to support the enzymatic activity. The difference in dynamics observed between the ensembles of proline isomerase at cryo and ambient temperatures indicates that flash freezing of a crystal anneals local conformational fluctuations and thereby removes protein dynamics that may be functionally relevant.

19In conclusion, this new method of modelling X-ray diffraction data reveals a wealth of detailed information about the dynamics of biomolecules that complements the high-resolution structural information already available from the crystallographic experiment. In depth understanding of structure–dynamics in biomolecules will enhance our insights into the molecular mechanisms that underlie biological processes.

20The method of ensemble refinement was implemented in the Phenix software. Adaptations and new procedures developed for ensemble refinement are given in section ‘Ensemble refinement methods’. Simulations were performed as described in section ‘Ensemble refinement protocol’. Details of the single-structure re-refinements used for comparison with the ensemble models are given in section ‘Single structure re-refinements’. Validation of the global disorder TLS model, the dependency of ensemble refinement on the starting structure and additional ensemble refinement calculations are given in ‘Additional ensemble refinement calculations’.

21Ensemble refinement was performed using phenix.ensemble_refinement, as will be made available in the next release of Phenix.

22The overall model structure factors are calculated as (1), defined by, incorporating overall anisotropic scaling and bulk solvent contributions.where, k is the overall scale factor, h is the column vector with Miller indices, A is the orthogonalisation matrix, Bcart is the anisotropic scale matrix, Fcalc is the structure factors calculated from atomic model, ksol and Bsol are the parameters for the flat bulk solvent model and Fmask are the structure factors calculated from bulk solvent mask.

23In order to restrain the instantaneous structures produced during the MD simulation with time and spatially averaged X-ray data, time-averaged restraints are used. This produces time-averaged (or rolling-average) structure factors such that (1) becomes (2).

24This is a time-dependent memory function, that is a ‘rolling’ average, where the size of the averaging window is controlled by the τx parameter (typically 1 ps). This averaging function is updated with the current individual structure every 10 time-steps (∆t) during the simulation and is implemented as in (3).

25Due to the stochastic behaviour of solvent molecules and the number of partially disordered or low occupancy sites, explicitly modelled solvent atoms are repositioned every 250 time-steps. Electron-density and difference density maps are generated using Fmodelt, excluding reflections in the free R set. Water oxygen atoms with an electron-density peak >1.0 σ in the 2mFobs − DFmodel map or a peak >3.0 σ in the mFobs − DFmodel map are preserved, otherwise the atom is removed. New water atoms are added for positions which have a 2mFobs − DFmodel peak >1.0 σ and a mFobs − DFmodel peak >3.0 σ, and are between 1.8–3.0 Å in distance to an existing atom. For high-resolution cases these criteria are adjusted to include mFobs − DFmodel map peaks >2.5 σ. Newly positioned atoms are assigned a random, Boltzmann-weighted, velocity. Explicitly modelled solvent atoms contribute to the atomic model (Fcalc).

26Bulk solvent is modelled using a solvent mask. The mask structure factors (Fmask) are averaged in the same manner as the atomic model (Fcalc) (4).

27The ksol and Bsol bulk solvent parameters and Bcart scaling parameters used for the duration of the simulation are calculated from the starting structure as described previously, they are re-optimized for the final ensemble.

28The overall scale factor, k, is constrained during the simulation. For the maximum-likelihood target function, as shown for acentric reflections, (5) during normalisation (6), the sum of the rolling-average structure factor array (2) is scaled to the sum of the structure factor array from the starting model (Fref) as shown in (7).where, Eobs and Emodel are the normalised structure factors, σA is the Sigma-A weighting factor, I0 is a modified Bessel function of order 0 and ε is the expected intensity factor.

29The simulations are performed such that the non-solvent atoms are at a target temperature (Ttarget) of 300 K, where the simulation is coupled to a velocity-scaled temperature-bath. The temperature bath is set to a value less than 300 K, typically 295–298 K. Because the X-ray restraints are computed from a time-dependent memory function, the X-ray energy term is non-conservative and thus heating occurs. During the equilibration phase the X-ray weight (wx-ray) is modulated by the temperature of the protein atoms (Tprotein) every 10 time-steps (∆t), such that the non-solvent atoms sample consistently at the target temperature (8).

30Thus, the thermostat offset controls the X-ray weight in a system independent manner whilst maintaining the target temperature. In the acquisition phase the X-ray weight is fixed to the averaged value used in the equilibration phase.

31The partitioning of inter-molecular disorder is performed before the start of the simulation using ADPs from the traditionally refined starting structure. TLS groups are assigned per molecule or domain as appropriate to model global packing disorder. For each group, TLS parameters are fitted to the ADPs of the starting structure for all non-solvent, non-hydrogen atoms. The agreement of the isotropic equivalents for the fitted TLS ADPs (Btls) and the reference ADPs (Bref) is scored as (9) for all non-solvent, non-hydrogen atoms.

32A percentile of atoms with the poorest fitting ADPs (pTLS) are excluded from the next round of TLS parameter fitting and repeated until the fitted TLS parameters converge. The converged TLS parameters are then applied to all atoms within that group for the duration of the simulation. Solvent atoms are assigned to the TLS group of the closest non-water atom, this assignment is updated every 250 time-steps. This TLS model produced lower R-values than using the ADP values from the re-refined single structure or using the one overall isotropic B-factor for all atoms in the model (Table 5).

33Structure factors for the final ensemble are calculated from the population of collected structure as in (2) where Fcalc and Fmask are defined as (10) and (11).

34The acquisition phase is split into several time blocks, in each of which 250 structures are typically stored. The R-values of all possible contiguous time blocks are calculated and the periods with the lowest Rwork are selected. This selection reduces the Rwork by 0–1.0% (mean improvement in 0.3%). For the 1YTT dataset with high quality experimental phases, the block selection for lowest Rwork corresponds well with the overall map correlation coefficient computed between the experimentally phased map and the map derived from the ensemble model (Figure 13). Next, to reduce the redundancy in the number of structures in the final ensemble (during the simulation thousands of structures are collected), we calculate the smallest number of structures that reproduce the Rfree within 0.1%. This is performed by iteratively parsing the stored structures with increasing time spacing (see Figure 1D). The overall and bulk-solvent scale factors are optimised for the final ensemble. The ensembles of structures are stored using the standard PDB format for multiple models, with B-factors listed as computed from the TLS model and overall B-factor scaling contributions.

35All atoms comprising the ensemble are assigned a probability (Pi) based on the positional likelihood of atom i in a given model relative to the complete ensemble of models. Fcalc electron-density maps are calculated for each model in the ensemble and 〈Fcalc〉 electron-density map is calculated for the complete ensemble as (10). Pi is calculated as (12).

36Calculating from an electron-density function allows for non-Gaussian distributions unlike RMSF, which is calculated from mean atomic position. These probabilities aid the visual inspection of the ensemble models and allow the observer to control the level of detail displayed (Figure 14).

37The starting structures were taken from the PDB server or from the PDB_REDO server if the Rfree was <0.25% than the equivalent PDB structure. We removed alternative positions and set corresponding occupancies to one. Overall anisotropic scale factors and solvent scale and B-factor (ksol and Bsol) were calculated based on these traditional single-structures (i.e. using the refined B-factor models). Next, the atomic B-factors were substituted by B-factors derived the global TLS disorder model (‘Materials and methods’—TLS approximation of the global disorder).

38At t = 0 〈Fcalc〉 and 〈Fsol〉 are set to Fcalc and Fsol. Boltzmann-weighted velocities are assigned to the atoms, corresponding to T = 300 K. The bath temperature Tbath used for velocity scaling is coupled to the X-ray weight (wx-ray) calculation, resulting in a temperature of 300 K for all non-solvent atoms. The simulation time-step used is 0.5 fs and the force-field parameterisation is as described. Simulations are started in parallel with varying values of pTLS (e.g. 0.2, 0.6, 0.8, 0.9, 1.0), τx (e.g. 0.25, 0.5 and 1.0 ps) and Tbath (e.g. 295 and 299 K). Water positions are picked according to electron-density criteria and updated every 250 steps. Every 10 time-steps rolling average structure factors, 〈Fcalc〉 and 〈Fsol〉, are updated for use in the time-averaged X-ray restraints. σA Values are updated if the Rfree of the rolling average model improves by >0.25%. The simulations have an equilibration phase (20τx) in which the temperature, X-ray weight and averaged structure factors stabilize. This is followed by an acquisition phase (40τx) where the values for wx-ray and σA are fixed and the structures for the final ensemble model are collected.

39CPU time for a dataset at 2.0-Å resolution with 199 residues in the asymmetric unit is 25 hr for each simulation using a 1.9 GHz processor.

40The single structure re-refinements used the same starting structure as the ensemble refinements (alternative conformations were not removed) and were re-refined using phenix.refine (version 1.7.1) and Buster (version 2.10.0). Standard parameters were used with the exception of optimizing the target weights and increasing the number of macro-cycles to 8 in phenix.refine. Explicit water refinement was performed, anisotropic ADPs were used if present in starting structure and TLS parameters were defined as used in the starting structure. PDB_REDO models were used as deposited.

41Burling et al. had previously shown that the choice of ADPs for ensemble refinement can affect the resultant structures. Three alternative ADP models were tested for seven datasets, as shown in Table 5. ADP model 1, ‘Global isotropic B-factor’, uses one overall isotropic B-factor applied to all atoms in the simulation. Multiple trials were performed to establish the optimum single value. For comparison the Wilson B-factor of the data is listed. ADP model 2, ‘Refined ADPs’, uses the ADPs from the refined single-structures. Best results were obtained by multiplying the refined ADPs by given scale factor. ADP model 3, ‘Basal TLS ADPs’, uses the basal TLS model with one TLS group per chain (including all non-hydrogen, non-solvent atoms) obtained as described in ‘Materials and methods’—TLS approximation of the global disorder, where pTLS is the percentage of atoms included the iterative fitting procedure. The basal TLS model returns the lowest Rfree values in all test cases.

42To test the effect of the starting structure three datasets (1UOY, 3CA7 and 1BV1) were re-refined with Buster, phenix.refine and Refmac as given by the PDB_REDO server. Each of these re-refined structures was used as the input structure for ensemble refinement, using the same run-time parameters. Each ensemble refinement was repeated three times using a different random number to generate the initial atomic velocities. The results are shown in Table 3. The mean Rfree (averaged over the random number seed repeats) of the resulting ensembles from the three different input structures are within 0.5%. The Fmodel cross correlation of ensemble pairs (best representative from each program, selected by Rfree) was calculated and is shown in Table 4. All ensemble pairs exhibit a cross correlation of greater than 0.99.

43Because occupancy and B-factor are strongly coupled in a traditional refinement, the occupancies of bound ligands and ions are typically set to unity, while the corresponding B-factors are refined in a single-structure refinement. In ensemble refinement, the B-factors are not refined, but are derived from the global TLS model and the atomic fluctuations.

44All simulations were initially performed with full occupancy for bound ligands and ions. In several cases, this resulted in excessive sampling of the ligand or ion, as seen when inspecting the ensemble and reported by the kinetic energies during the simulation, which were far in excess of neighbouring protein atoms. These observations indicate that the corresponding occupancy of the bound ligand or ion is less than one. In these cases the occupancies were lowered and the simulations repeated until the kinetic energy of the ligand or ion were equivalent to the proximal protein components.