In this work, we illustrate our approach for the protein calmodulin in complex with the smooth muscle Myosin Light Chain Kinase (smMLCK) peptide. Throughout, we will use a previously solved crystal structure of this complex [PDB ID: 1cdl (Meador et al., 1992)] as a reference.

In previous studies, specific nitroxide spin-labeled target peptides that bind to calmodulin were used; in this manner it was possible to map out the orientation of the peptide with respect to the protein (Zhang et al., 1995; Yuan et al., 2004). In this work, we collected PRE data for a total of ten spin-labelled protein sites (Figure 1). Nine of the sites on the protein were chosen to be solvent-exposed and within secondary structure elements by manual inspection of predicted secondary structure and solvent exposure. The remaining site, C149, was a single-residue extension of the C-terminus. To better simulate the process for a system without a previously determined structure, the protein’s known structure was not used in choosing the spin-labelled sites. Indeed, we learned later that several of the selected sites provided little useful information because they are either distant from the rest of the protein or they provided information that is mostly redundant with that obtained from other labeling sites.

To simulate the limited availability of isotopically labelled peptide, either due to cost or difficulty of production, only four of the ten spin-label data sets (chosen randomly) were collected with isotopically labelled peptide. The remaining six data sets were collected without labelled peptide, which results in the peptide being present but unlabeled and undetectable in the ¹H, ¹⁵N HSQC NMR experiments.

The ten single-cysteine point calmodulin (CaM) mutants (S17C, T34C, N42C, N53C, R86C, T110C, T117C, E127C, Q143C, 149C) were made by standard site-directed mutagenesis methods for attachment of the thiol-specific nitroxide spin label (1-oxyl-2,2,5,5-tetramethyl-δ-3-pyrroline-3-methyl) methanethiosulfonate (MTSL, Toronto Research Chemicals). Correctness of the mutations was confirmed by DNA sequencing. Calmodulin contains no Cys residues, so highly site-specific labeling can be obtained in this manner. MTSL is a relatively compact, yet highly reactive molecule compared to other commercially available nitroxide spin-labels that have been used to modify Cys residues; its shorter more rigid structure would be an advantage for the PRE studies [for discussion see for example (Guo et al., 2008; Fawzi et al., 2011)]. ¹³C and ¹⁵N-labeled CaM was expressed in M9 minimal medium with 99.9% ¹⁵NH₄Cl and ¹³C₆-glucose (0.5 gr/L and 3 gr/L, respectively; Cambridge Isotopes Laboratories) as isotope sources. Proteins were expressed and purified as described previously (Liu and Vogel, 2012; Ishida et al., 2016).

We followed a standard protocol for attaching the nitroxide spin-label to each single-cysteine CaM mutants with the spin-labelling reagent MTSL (Battiste and Wagner, 2000). To prepare the CaM/smMLCK complex sample, a 1.2-fold excess of either unlabeled or labelled peptide was mixed with each CaM mutant protein. All preparations were divided into two NMR samples. One sample was reduced to inactivate the spin-label by adding a 3-fold excess of ascorbic acid.

A construct with a 6xHis-KSI (D38A) fusion-protein tag was generated for smMLCK peptide expression in Escherichia coli (Jaroniec et al., 2005; Ishida and Vogel, 2010). The ketosteroid isomerase (KSI) coding sequence generates an insoluble protein, and this directs the protein-peptide fusion directly into inclusion bodies, where they are protected from proteolytic cleavage (Hwang et al., 2014). A linkage sequence “GGGGSSDP” with the Asp-Pro acid cleavage site was designed between the KSI protein and the sequence of the smMLCK peptide. The entire 6xHis-KSI-GGGGSSDP-smMLCKp gene sequence was inserted between the NdeI and XhoI sites of the pET15b(+) plasmid (Novagen), which was subsequently transferred into BL21(DE3) E. coli cells for protein expression. The cells were grown in either LB media (for unlabeled peptide) or minimal M9 media (containing ¹³C₆-glucose and ¹⁵NH₄Cl isotope to produce isotope-labeled peptide) and they were induced at OD₆₀₀ = 0.6 with 1 mM IPTG for 4 h at 37°C. A cell lysate was prepared as previously described. The insoluble fusion protein was separated after one hour of centrifugation (18,000 rpm) and then resuspended in 6 M guanidine hydrochloride. Impurities were removed before the insoluble proteins can be extracted with metal chelate chromatography on a nickel affinity column. After extensive dialysis with double distilled H₂O, the precipitated insoluble protein was collected and the Asp-Pro bond was cleaved in 10% formic acid at 80°C for 90 min (Hwang et al., 2014). The protein-peptide mixture was flash frozen with liquid nitrogen and lyophilized. Insoluble proteins and other impurities were removed after the lyophilized mixture was resuspended in a 20 mM Tris-HCl buffer (pH = 8.0). Finally, the unlabeled and isotope-labeled smMLCK peptides were purified with reverse-phase HPLC (COSMOSIL 5C₁₈-AR-300, Nacalai United States). All purified peptides were lyophilized and stored at −20°C for further use. The final peptide sequence after cleaving is PARRKWQKTGHAVRAIGRLSS. The N-terminal proline is not observable in the NMR experiments and was not included in modeling with MELD.

All NMR experiments were carried out on a 600 MHz Bruker AVANCE spectrometer with a field strength of 14.1 T. Backbone resonance assignments for the protein and the bound peptide were confirmed with the following 3D experiments: HNCO, HNCA, HNCOCA, HNCACB, and CBCA(CO)NH, as described previously (Liu and Vogel, 2012). All data were processed using NMRPipe (Delaglio et al., 1995) and analyzed with the program NMRView (Johnson and Blevins, 1994). All chemical shifts from these experiments were used to obtain backbone torsion angles from the program TALOS+ (Cornilescu et al., 1998; Shen et al., 2009). Secondary structure elements as identified through the assigned chemical shifts were as expected based on the known structure.

Two ¹H, ¹⁵N HSQC spectra were obtained for each spin-label construct. Each system contained each ¹⁵N-labeled protein and either unlabeled or ¹⁵N-labeled peptide, depending on the spin-label site (S17C, N53C, T127C, and 149C had isotopically labeled peptide). One HSQC was collected with active spin-label, whereas the other HSQC was collected with reduced, inactivated spin-label. The distances between the spin-label and the affected nuclei were calculated using the two-time point method (Iwahara et al., 2007).

Finally, to supplement the PRE experiments, we obtained RDC measurements for the amide groups in the complex with a sample where both the protein and peptide are isotopically labelled. Residual dipolar couplings (RDC) were measured for the CaM/smMLCK complex sample in a partially aligned media, which contains 2 mM bis-Tris (pH = 7.0), 300 mM KCl and 16 mg/ml Pf1 bacteriophage (Asla Biotech Ltd.). The IPAP-HSQC experiment was used for the RDC measurements (Ottiger et al., 1998). In these experiments the effects of dipole-dipole cross-correlated relaxation can impact the accuracy of ¹J_NH splitting measured from the spectra introducing a small residual bias in the RDCs. While these systematic errors can be eliminated by using a selectively-decoupled sequence (Yao et al., 2009), the errors are small relative to the magnitude of the measured RDCs and are expected to have a minimal effect on structure determination (Yao et al., 2009). Our work uses only a single RDC alignment. Notably, a mutant of Calmodulin is capable of selective binding to lanthanides, which provides a strategy for the measurement of multiple RDCs (Bertini et al., 2009). A quantitative assessment of protein mobility/heterogeneity by RDC would require the use of multiple alignments (Barbieri et al., 2002; Tolman, 2002; Bouvignies et al., 2006; Higman et al., 2011; Guerry et al., 2013; Andrałojć et al., 2015). However, as discussed further below, it is not currently possible to conduct such an analysis with the MELD approach, as MELD compares individual structures, rather than ensembles of structures, to the experimental data.

Here, we employ MELD, a physics-based Bayesian approach for structural determination to infer the ensemble of structures most consistent with the known physics of protein structure and experimental data (MacCallum et al., 2015; Perez et al., 2015). MELD uses a Bayesian framework to combine a physics-based prior distribution with a data likelihood function to make statistically consistent inferences about conformations that explain the experimental data.

MELD uses Bayes’ theorem: p ( x | D ) ∝ p ( x ) p ( D | x ) , (1) where x represents the atomic coordinates and D represents the data. The physics-based prior, p ( x ) , specifies which structures are more likely a priori and determines the distribution of structures in the absence of data. In the present study the physics-based prior is given by the Amber ff14SB force field (Maier et al., 2015) with a grid-based torsion potential (Perez et al., 2015) and the OBC generalized-Born implicit solvent model (Onufriev et al., 2004). The likelihood function, p ( D | x ) , captures the compatibility between the data and some structure x . In MELD, the likelihood function takes the form of a unique restraint function (MacCallum et al., 2015), explained in more detail below. The goal of Bayesian inference is to compute the posterior distribution, p ( x | D ) , which is the most statistically consistent inference given the prior, likelihood, and data.

As discussed in the Results section below, the term ensemble is highly overloaded in structural biology and care is required in interpretation. MELD belongs to the class of methods where a single structure, rather than entire ensemble of structures, is considered in the likelihood function, such that each member of the ensemble individually agrees with the experimental data. Any conformational heterogeneity (e.g., flexible loops) may represent true intrinsic heterogeneity, but may also simply reflect a lack of data. As such, MELD produces a form of uncertainty ensemble in the terminology of Bonomi et al., 2017.

The input to our approach is: 1) the protein sequence, 2) TALOS+ secondary structure predictions derived from backbone chemical shifts (Cornilescu et al., 1998; Shen et al., 2009; MacCallum et al., 2015), 3) distance restraints derived from PRE measurements, and 4) orientational restraints derived from RDC measurements. We have recently demonstrated the success of a similar approach for PRE measurements in solid-state NMR (Perez et al., 2019).

PRE data is often both noisy and sparse (Kim et al., 2014), which makes structural inference challenging. Despite collecting data for ten spin-label positions, we can derive only a few distance restraints that are short in spatial distance (in this case, we define short as <12 Å) (Figure 2). Of these short spatial distances, only a small number correspond to residues that are distant in sequence, which would provide the most information about folding (Gottstein et al., 2012). Furthermore, as stated previously, to simulate the limited availability of isotopically labelled peptide, only four of ten datasets (S17C, N53C, T127C, and 149C) had labelled peptide, and there are no short distance PREs between the peptide and the protein. This leaves the peptide’s correct placement to be dictated by longer, less informative spatial distance restraints and the physical model, which makes accurate inference more challenging.

Our first step was to develop a consistent method to convert ensemble-averaged PRE measurements into distance restraints. PRE data were turned into approximate distances using the Solomon-Bloembergen equations following the standard approach (Battiste and Wagner, 2000; Iwahara et al., 2007).

For nitroxides, Curie-spin relaxation is negligible and the transverse relaxation enhancement, Γ₂, is dominated by direct dipole-dipole interactions (Clore and Iwahara, 2009). In this case, Γ 2 is related to the distance between the paramagnetic center and the observed nucleus, r: Γ 2 = K r 6 [ 4 τ c + 3 τ c 1 +   ω H 2 τ c 2 ] where ω_H is the Larmor frequency of the proton, and τ_c is the correlation time for the electron-nuclear interaction defined as τ c =     ( τ r − 1   +   τ s − 1 ) − 1 where τ_r is the rotational correlation time, and τ_s is the electron relaxation time. Previous experiments (Lee et al., 2002) have determined that τ c ≈ 9.5   n s for Calmodulin in complex with smMLCK peptide at 25°C. For nitroxides, the electron relaxation time is long (Iwahara et al., 2007) ( τ s > 10 − 7 s ), so the rotational correlation time dominates and τ c ≈ τ r ≈ 9.5   n s . K is given by: K = 1 15   ( μ 0 4 π ) 2 γ 1 2 g 2 μ B 2 S ( S + 1 ) where µ₀ is the permeability of vacuum, γ₁ is the nuclear gyromagnetic ratio, g is the electronic g factor, µ_B is the Bohr magneton, and S is the electron spin quantum number. Γ₂ can be estimated using a two-point time measurement (Iwahara et al., 2007): Γ 2 = R 2 , p a r a − R 2 , d i a = 1 Δ T ln I d i a ( T b ) I p a r a ( T a ) I d i a ( T a ) I p a r a ( T b ) where Δ T is a time delay chosen to minimize the error in Γ₂, and I_dia and I_para are the peak intensities for the diamagnetic and paramagnetic samples, respectively (Iwahara et al., 2007). In this work, we use Δ T = 20   m s .

The distances derived from PRE data correspond to ensemble averages with an r ⁻⁶ weighting, but in MELD (and most other structure determination software), restraints are applied to single structures rather than ensembles. To account for conformational heterogeneity of both the protein and the flexible spin-label, the PRE-derived distances are turned into flat bottomed harmonic restraints that allow for a range of distances without penalty. This approach is a tradeoff that ensures that individual structures are not erroneously over-restrained but this can allow discrepancies between the measured and modelled ensemble averages. Our aim is to produce an approximately correct ensemble starting from an extended chain. If desired, the resulting ensemble can be further refined using a variety of ensemble approaches (Boomsma et al., 2014; Hummer and Köfinger, 2015; Bonomi et al., 2016; Gaalswyk et al., 2018).

We divided the data into short, medium, and long distances with corresponding upper and lower bounds (Table 1). Short and long distances are difficult to quantify with precision because the peak is either completely broadened for residues close to the spin-label or barely changes intensity for those that are far away. These distances are turned into broad restraints that either start from zero or extend to infinity for short and long, respectively. Medium distances correspond to peak intensity changes that can be quantified more precisely and are turned into restraints centered around the predicted value. All distances include a buffer of ±5 Å of the measured distance to account for the flexibility of the spin-label and noise in the experimental data (Perez et al., 2019).

Due to noise in the experimental data, partially overlapping peaks, and instantaneous fluctuations in both the protein structure and the position of the spin label, we observed that restraints are sometimes violated even with a ±5 Å buffer. To mitigate this issue, we used MELD’s unique ability to require that only a certain fraction of the restraints must be satisfied by each structure. We set this active fraction to 0.9. Essentially, as long as 90 percent of the restraints are satisfied, the resulting restraint energy will be zero. We treat the remaining restraints as being derived from spurious data, so they are entirely ignored. Every timestep, MELD decides which restraints are active based on the current structure. Further details can be found in the SI and in MacCallum et al. (2015).

In our approach, the various hyperparameters (boundaries between short/medium/long, size of buffer, active fraction) are fixed. One potential improvement would be to place a hyperprior on these values and infer them using an extended Bayesian approach like Inferential Structure Determination (Rieping et al., 2005). This would allow the data and physical model to determine the most likely values of these hyperparameters, rather than requiring their specification a priori. As MELD does not currently support inference of hyperparameters, we chose the simpler approach of setting a wide buffer and lower active fraction, which potentially sacrifices a small amount of information.

The spin-label was modeled using virtual sites (Banci et al., 1996) following the approach of Islam and Roux (Islam et al., 2013; Islam and Roux, 2015). These virtual sites represent the spin-label as a non-interacting dummy particle to simplify the simulation without losing relevant information for structural refinement. These simplified dummy nitroxide spin-labels are parameterized to match the spin-labels’ 3D spatial distribution and dynamics in all-atom simulations. The virtual sites are non-interacting, allowing us to account for all ten spin labels in a single simulation without the risk of interactions between them.

Secondary structure restraints were derived from TALOS+ (Cornilescu et al., 1998; Shen et al., 2009) and used to restraint MELD simulations as previously described (MacCallum et al., 2015). Our approach works by first breaking the protein into overlapping 5-residue fragments. If 4/5 of the residues in the fragment are predicted to be helical or extended, then the fragment is restrained using a combination of torsion and distance restraints (MacCallum et al., 2015). All secondary structure restraints are then combined into a collection with an active fraction of 0.95, which allows 5 percent of fragments to differ from their predicted secondary structure.

The traditional approach to incorporate RDCs into simulations is based on solving for the optimal alignment tensor, which requires solving a system of equations every time step using singular value decomposition (SVD) or related methods, which can be computationally intensive (Losonczi et al., 1999). We found this to be particularly problematic in the GPU-accelerated framework of MELD, where this traditional approach led to a 300 percent increase in run time (data not shown), primarily due to the extreme speed of the rest of the force/energy calculations and the challenge of efficiently parallelizing SVDs for small systems of equations on a GPU. To mitigate this issue, we instead followed the approach in Habeck, Nilges, and Riepling (Habeck et al., 2008), which we implemented using an OpenMM CustomCentroidBondForce (Eastman et al., 2017). In our implementation, the alignment tensor elements are encoded in two non-interacting dummy particles coupled to the rest of the system through an additional energy term. This approach has two benefits. First, it is dramatically faster than the standard approach on GPUs, with negligible cost compared to the calculation of the non-bonded forces. Second, this approach accounts for uncertainty and produces a joint distribution of alignment tensors and structures, providing a Bayesian posterior estimate of the conformational ensemble that better reflects uncertainty. A full explanation of our implementation can be found in the SI. To account for uncertainty in the experimental data and to avoid erroneously over-restraining individual structures to the ensemble average data, we use a flat-bottomed restraint where the energy is zero if the computed RDC is within 1.5 Hz of the measured value. Another approach that avoids the need to solve for the alignment tensor is given in Camilloni and Vendruscolo (2015).

The term “ensemble” is highly overloaded in structural biology, with a variety of meanings in different contexts (Bonomi et al., 2017; Andralojc and Ravera, 2018; Gaalswyk et al., 2018). Care must be taken to ensure correct interpretation.

MELD samples from a well-defined conformational ensemble (Gaalswyk et al., 2018), specifically the Bayesian posterior distribution given by Eq. 1. Interpretation of this ensemble is straightforward: it is the statistically consistent posterior distribution inferred from the prior, likelihood, and data. How should one select or report structures from this ensemble? The approach we take here is simply to report all structures, as this fully captures the heterogeneity of the distribution. If there is a limit to the number of structures reported, one simple, correct approach is to select a subset of structures at random. Alternatively, one could cluster the structures and report the cluster medoids and populations along with some measure of the variance of structures within the cluster. A variety of approaches are supported by the PDB-Dev archival system which is being developed for structural models obtained using integrative modeling (Vallat et al., 2018). However, we note that since MELD samples structures with the correct posterior probabilities, it is incorrect to further select structures based on other criteria, such as selecting the lowest energy structures.

A second consideration in ensemble interpretation is the nature of the likelihood function. The experimental measurements are averages over a thermodynamic ensemble. The most correct modeling approach is to use a likelihood function that considers an entire ensemble of models, ensuring that the predicted average quantities match their corresponding experimental measurements (Bonomi et al., 2017; Andralojc and Ravera, 2018; Gaalswyk et al., 2018). This is an ill-defined inverse problem (Bonomi et al., 2017; Andralojc and Ravera, 2018), so regularization is required, typically in the form of physical modelling and entropy maximization (Bonomi et al., 2017; Andralojc and Ravera, 2018; Gaalswyk et al., 2018). While conceptually appealing, ensemble likelihood methods are complex with high computational requirements. Alternatively, most methods in structural biology, including the MELD approach described here, use single-structure likelihoods (Boomsma et al., 2014). These methods are overly restrictive, as they require each member of the ensemble to be consistent with the data to within some tolerance. In the current approach, we use relatively wide tolerances, but this still does not guarantee that that the computed ensemble accurately models the true distribution.

The primary issue is that for a given set of experimental measurements, there are many possible ensembles that could produce it. The ensemble that MELD generates ensures that each structure is in reasonable agreement with the data and allows for a reasonable degree of flexibility. However, the MELD average might not precisely match the experimental measurement due to the use of wide tolerances. Furthermore, the true ensemble could be “broader” than the one generated by MELD—the true ensemble could have many structures that are individually in poor agreement with the data, while still having the same ensemble average, see Figure 1 of Gaalswyk et al. (2018) for a simple illustration. In the terminology of Bonomi et al. (2017), MELD produces an uncertainty ensemble, where heterogeneity in the calculated ensemble could represent true heterogeneity in the system or could simply reflect a lack of data for some part of the protein. The single-structure approach is likely reasonable when the true ensemble has only a modest amount of heterogeneity, e.g., small fluctuations around an average structure, but could be expected to break down for highly heterogenous systems, e.g., systems containing intrinsically disordered regions.

Although we do not pursue it here, a promising approach would be to use a method like MELD to compute an initial approximate ensemble that could then be used as a starting point for ensemble approaches (Boomsma et al., 2014; Hummer and Köfinger, 2015; Bonomi et al., 2016; Gaalswyk et al., 2018).

To determine how the experimental data should be incorporated, we performed several simulations varying in their set up (Trial1–Trial4). We explored various ways of combining the restraints into collections (Table 2). In MELD, at every timestep, the restraints in a collection are sorted by energy, and the active fraction with the lowest energy are “active” and contribute their forces and energy to the system, while the remainder are “inactive” and ignored. The division of restraints into collections matters because it determines how MELD decides which restraints are active and which are ignored. For example, Trial1 combines all of the restraints into a single collection. In this case with an active fraction of 0.9, MELD can freely ignore any 10 percent of the restraints, which could be, for example, ignoring one of the ten spin labels entirely. Trial2 separates the restraints by spin-label and into short, medium, and long-distance ranges, resulting in 30 collections. Now MELD can only ignore 10 percent from each spin label/distance combination, while the remaining 90% will be active. Trial3 extends Trial2 by adding the RDC restraints. Trial1–Trail3 start from an extended conformation generated by the tleap tool from the AmberTools suite (Case et al., 2021). Trial4 follows the same protocol as Trial3 but starts from the reference crystal structure as a control.

For each trial, we ran a 2.5 µs replica exchange simulation using 48 replicas. The temperature and the force constant for each restraint collection varied across replicas (see SI for details). The last 0.5 µs of the lowest replica was used for analysis.

We compare the trials using kernel density estimation plots (KDE; see Supporting Information for details) of the backbone root mean square deviation (RMSD) to the reference structure [PDB: 1cdl (Meador et al., 1992)], excluding the flexible tails at the N and C terminals of the protein which are not present in the reference (Figure 3).

Trial1 is the most straightforward approach and combines all of the data into a single collection. Many of the structures have relatively low RMSD to the reference (<4Å), which is promising considering the rather limited experimental data, but there are also structures with much higher RMSDs of ∼5 Å and ∼7.5 Å RMSD. There are various explanations for these high RMSD conformations, but perhaps the simplest is that this way of grouping all restraints into a single collection allows MELD to ignore short spatial distance restraints that would otherwise eliminate these conformations. As previously stated, the utility of a restraint depends on its spatial distance. Shorter distances provide highly constraining information. However, this highly constraining nature means that these restraints are more difficult to form, leading MELD ignore them in favour of more easily satisfied restraints.

To test this hypothesis, in Trial2, we further subdivided the restraints by separating the short, medium, and long restraints from each dataset into separate collections, resulting in 30 total collections.

The resulting RMSD distribution is centered at a modest RMSD of ∼4 Å, which is slightly worse than the mode from Trial1. However, this method of combining restraints into collections has wholly eliminated the high RMSD conformations.

Although the RMSDs obtained are only modest (∼4 Å), these results were obtained with a very sparse dataset with only one spin-label per 17 amino acids. This equates to 6.4 total restraints per residue, and only 0.8 short-distance restraints per residue. For context, NOE-based structures from fully protonated samples typically have >15 NOE restraints per residue, all with short distances.

Based on visual examination, several of our spin-label sites appear to give restraints that are largely uninformative, either because they are far from the remainder of the protein (e.g., R86C) or because they are mostly redundant with other spin-label positions (e.g., T110C), see Figures 1, 2. Our results could be improved with a more judicious choice of the 10 label sites, but it is unclear how to do this without pre-existing knowledge of the structure. The results could likely be improved further by adding additional spin-label sites using the calculated structural ensemble to optimize probe location, although we do not pursue this here. Such an iterative strategy could be a viable approach to improve model accuracy but comes at an additional experimental cost. More rigid spin-labels (Fawzi et al., 2011) could also improve results, as MTSL still displays significant conformational heterogeneity that results in less precise distance restraints.

Despite collecting data for 10 different spin-label sites, few yielded informative short spatial distance, high sequence distance restraints (Figure 2), limiting the models’ achievable accuracy to relatively modest RMSDs of around 4 Å. Rather than collect additional PRE data, we instead chose to explore the utility of combining PRE information with residual dipolar couplings (RDCs) measured for the amide groups.

Residual dipolar couplings provide information about the orientation of amide NH bonds complementary to the distance information from PRE experiments. In Trial3, we combined PRE information (using the same strategy as Trial2) with RDC data. The inclusion of RDC data led to a substantial improvement in the RMSD (Figure 3). The RMSD ranges from approximately 1.6–4.0 Å with an average RMSD of 2.8 Å, including both Calmodulin and the smMLCK peptide (Figure 4). This improvement of RMSD upon inclusion of RDC data is consistent with previous studies showing that RDC data provides valuable information on the relative orientation of the two lobes of calmodulin (Mal et al., 2002; Gifford et al., 2011).

RDCs provide information about how the amides are oriented, which, when combined with secondary structure restraints (derived from the measured backbone chemical shifts) and distance restraints (derived from PREs), serves to dramatically limit the possible structures that simultaneously agree with the experimental data and the physical model.

To assess the potential quality of sidechain packing, we examined the single best structure obtained during our simulations, which has an RMSD of 1.6 Å (Supplementary Figure S1). For this “best” structure, the RMSD for sidechain heavy atoms is 2.1 Å for all sidechains and 1.4 Å for core sidechains. This is notable as there are no restraints on the side chains themselves, only between the spin labels and the backbone amide protons. This packing phenomenon with MELD has been noted previously and can be attributed to the accuracy of the physical model (Perez et al., 2019). However, we note that the sidechain and backbone RMSDs are generally correlated, and this structure has a lower backbone RMSD than average, so the average side chain RMSDs will be higher than these figures.

As noted previously, the experimental data contained no short distance PREs to the peptide, so placement of the peptide is dependent on a combination of medium and long restraints with the physical model. Furthermore, only 4 of 10 experiments contained labelled peptide, with the peptide undetected in the remaining experiments. Nevertheless, the combination of available data and the physical model was still able to routinely position the peptide correctly (the peptide is included in the RMSD calculations shown in Figure 3). The structure of calmodulin depends on the peptide and its binding (Barbato et al., 1992; Ikura et al., 1992), so correct placement of the peptide is critical.

Calmodulin consists of two lobes connected by a flexible linker that becomes structured upon peptide binding (Barbato et al., 1992). Examination of each lobe individually shows that our modeled ensembles are tightly clustered (Figure 5), indicating that most of the heterogeneity in our calculated ensemble arises from the relative motion of the two lobes. If we consider the RMSD of each lobe to the reference individually, the results are consistently lower than for the whole protein (Figure 6). The RMSD of the C-lobe to the reference structure is ∼2.1 Å (Figure 6B), which is consistent with typical RMSDs for small globular proteins seen in MD simulations. The results for the N-lobe are similar. The resulting heterogeneity in the relative orientation of the two domains should be interpreted with caution, due to the use of a single-structure likelihood, as discussed above.

A previous study (Carlon et al., 2019) examining the joint X-ray/NMR refinement of Calmodulin in complex with the Death-Associated Protein kinase (DAPk) peptide revealed poor agreement between X-ray and NMR data due to large interprotein contacts in the crystal that stabilize a conformation that is in poor agreement with the solution NMR data. The crystal structure of the full-length DAPk protein in complex with Calmodulin lacks these contacts and is in much better agreement with the NMR data. These results highlight the need for caution when comparing structures determined by X-ray crystallography and NMR, particularly in cases where flexibility can be expected.

A study of 109 pairs of NMR and crystal structures (Sikic et al., 2010) showed that typical C⍺ RMSDs range from ∼0.5 to 4 Å with a mean of ∼2.0 Å when using the DALI (Holm and Sander, 1993) alignment. The typical variability of models within a given NMR ensemble was similar (Sikic et al., 2010). Our results for the individual lobes of calmodulin give similar average RMSDs, indicating that our approach is producing results comparable to typical NMR structures using NOEs and fully protonated samples despite the substantial sparsity in our data. Our results for the full-length complex produce a slightly higher average RMSD, which reflects heterogeneity in the exact relative placement of the two lobes.

To further test our predictions’ quality, we also ran calculations using the same protocol as Trial3 but starting from the reference crystal structure rather than from an extended chain (Figure 6), which sets a bound on the possible accuracy that could be obtained. The resulting RMSD distributions are similar to our predictions. This indicates that given: 1) the available experimental data, 2) potential limitations of the physical model used, 3) the use of a single-structure rather than ensemble likelihood, and 4) the challenges of comparing with a static crystal structure, the results obtained using MELD are essentially as good as they could be.

Each calculation was over 48 replicas for 2.5 µs, which required approximately 6 days on 48 GTX 1080Ti GPUs. However, examination of the RMSD over time (Supplementary Figure S2) shows that the simulations appear to be converged after ∼500 ns. In hindsight, the simulation length could have been reduced to 1 µs without a loss in quality, which would reduce simulation time to 2.5 days. While computationally expensive, our approach is readily feasible with access to advanced research computing or cloud computing resources.