The ribosome lowers the entropic penalty of protein folding

Protein expression and purification

DNA constructs of FLN5 were previously described^7,11. Coding sequences for titin I27 and HRAS were introduced into the pLDC-17 vector using standard procedures. Further mutations were introduced using site-directed mutagenesis; for ¹⁹F labelling, amber stop codons were introduced⁶ in position 32 in HRAS, and residue 14 with an additional K87H point mutation in I27. FLN5 variants were expressed as His-tagged proteins and isotopically labelled in Escherichia coli BL21 DE3-Gold cells as previously described^6,7; an identical protocol was used to produce purified samples of I27 and HRAS. RNC constructs comprised an arrest-enhanced variant of the SecM stalling sequence, FSTPVWIWWWPRIRGPP, as previously described⁶. Purification of isolated FLN5 A₃A₃ was performed by affinity chromatography followed by size-exclusion chromatography in the presence of 6 M urea prior to buffer exchange into Tico buffer (10 mM Hepes, 30 mM NH₄Cl, 12 mM MgCl₂, 1 mM EDTA). The full protein sequence of the FLN5 A₃A₃ is deposited together with its chemical shift assignment on the BMRB (entry 51023). For the RDC, pulse-field gradient NMR (PFG-NMR) and PRE-NMR experiments, the additional mutation C747V (referred to FLN5 A₃A₃ V747) was introduced to yield a cysteine-less construct for site-specific spin labelling. The protein concentration was determined using the BCA assay according to the manufacturer’s instructions. RNCs were expressed, isotopically labelled uniformly with ¹⁵N, or site-specifically with ¹⁹F, and purified as previously described^6,7. For samples for intermolecular PRE-NMR experiments involving ribosome labelling, we generated modified E. coli BL21 strains with cysteine mutations in uL23 and uL24 using CRISPR as previously described². RNC samples were prepared in Tico buffer for experiments. Western blot analyses were undertaken with an anti-hexahistidine horseradish peroxidase-linked antibody (Invitrogen, 1:5,000 dilution).

Fluorescent and PEG-maleimide labelling of 70S and RNC samples

Ribosomes and RNCs were first reduced using 2 mM TCEP overnight at 277 K, then buffer exchanged into labelling buffer. For fluorescein-5-maleimide and PEG-maleimide, labelling was performed in Tico at pH 7.5. ABD-MTS labelling was performed in labelling buffer (50 mM HEPES, 12 mM MgCl₂, 20 mM NH₄Cl, 1 mM EDTA, pH 8.0). Samples were labelled using a 10x molar excess of ABD-MTS, or fluorescein-5-maleimide. Cysteine mass-tagging by PEGylation was performed as previously described with 10,000-fold molar excess of PEG over sample⁷. ABD-MTS and PEGylation reactions were analysed using 12% Bis-Tris SDS–PAGE gels⁶¹. The fluorescein-labelled reactions were run on a 20% Tricine SDS–PAGE gel, modified from ref. ⁶².

NMR spectroscopy

All NMR experiments were recorded with Topspin 3.5pl2. NMR experiments of FLN5 A₃A₃ were performed in Tico buffer at pH 7.5 and 283 K. Chemical shifts were previously assigned⁷ and obtained from data recorded on a Bruker Avance III operating at 700 and 800 MHz equipped with TCI cryoprobes. All samples contained 10% (v/v) D₂O and 0.001% (w/v) DSS as a reference. Data were processed analysed using NMRPipe⁶³ (v11.7), CCPN⁶⁴ (v2.4) and MATLAB (R2017a, Mathworks).

Amide ¹H and ¹⁵N chemical shifts were obtained from two-dimensional ¹H–¹⁵N SOFAST-HMQC experiments⁶⁵ using an acquisition time of 50 ms in the direct dimension. The inter-scan delay was 50 ms. Cα chemical shifts were obtained from 3D BEST-HNCA experiments recorded at 800 MHz with acquisition times of ~50 ms and inter-scan delays of 150 ms. C’ chemical shifts were obtained from BEST HNCO experiments recorded at 700 MHz using acquisition times of ~50 ms and inter-scan delays of 200 ms. RNC samples were doped with 20 mM NiDO2A (Ni(ii) 1,4,7,10-tetraazacyclododecane-1,7-bis(acetic acid)) to enhance sensitivity⁶⁶. Cosine-squared window functions were used in processing the spectra.

For PRE-NMR experiments, we used a cysteine-less construct with the C747V mutation and introduced six and eight labelling sites in the isolated and ribosome-bound protein, respectively. Samples were reduced overnight at 277 K in Tico supplemented with 2 mM TCEP. TCEP was then removed by buffer exchange into labelling buffer (50 mM HEPES, 12 mM MgCl₂, 20 mM NH₄Cl, 1 mM EDTA, pH 8.0) and subsequently labelled overnight at 277 K with 10× molar excess of MTSL. Following labelling, excess MTSL was removed by buffer exchanging the sample back into Tico buffer for NMR. The same labelling protocol was used for isolated protein and RNC samples. To measure the PREs, we recorded the signal intensities with MTSL in the paramagnetic and diamagnetic state. Direct measurements of relaxation rates proved not feasible for RNC samples due sensitivity limitations. 2D ¹H–¹⁵N SOFAST-HMQC experiments⁶⁵ were recorded at 800 MHz and 283 K using ~100 μM of protein or ~10 μM of RNC. Experiments were recorded with an acquisition time of 100 ms and 35 ms, in the direct and indirect dimension, respectively. The inter-scan delay was 450 ms to allow for complete relaxation. To acquire the diamagnetic data, the sample was reduced with 2.5 mM (RNC) or 100× molar excess (isolated) sodium ascorbate. Following complete reduction, the same HMQC experiment was recorded. To extract the PREs, spectral peaks were first fitted to a Lorentzian shape in both the direct and indirect dimension using NMRPipe⁶³. Errors were obtained from the spectral noise (RMSE). From the fitted peaks, intensity ratios of I_para/I_dia were calculated and converted to PRE rates for Bayesian ensemble reweighting by numerically solving equation S34 (see Supplementary Notes 3–4) for Γ₂. Sample integrity was monitored using interleaved ¹H,¹⁵N SORDID diffusion measurements as previously described⁷.

PFG-NMR experiments were used to measure the diffusion coefficients and the R_h of FLN5 variants. 1D ¹H,¹⁵N-XTSE diffusion measurements were recorded at 700 (FLN5, FLN5 Δ6) and 800 MHz (FLN5 A₃A₃). Eight to sixteen gradient strengths ranging linearly from 5% to 95% of the maximum gradient strength of 0.556 T m^−1 were used. By measuring the signal intensity at each gradient strength, diffusion coefficients could be obtained by fitting the data to the Stejskal–Tanner equation⁶⁷, which were converted to R_h using the Stokes–Einstein equation.

RDCs for isolated FLN5 A₃A₃ were measured in Tico buffer at 283 K and pH 7.5 in a PEG/octanol mixture⁶⁸. RDCs are reported as the splitting of the isotropic splitting subtracted from the aligned splitting, corrected for the negative gyromagnetic ratio of ¹⁵N. RDCs were measured by preparing a solution containing 4.6% (w/w) pentaethylene glycol monooctyl ether (C₈E₅), 1-octonal (molar ratio 1-octanol:C₈E₅ = 0.94) and 110 μM of protein. Alignment was confirmed by measuring the D₂O deuterium splitting at 283 K (17.6 Hz). All RDC NMR experiments were acquired on a Bruker Avance III HD 800 MHz spectrometer equipped with a TCI cryoprobe. A set of four different RDCs (¹D_NH, ¹D_CαCO, ¹D_CαHα and ²D_HNCO) was measured per sample (isotropic and anisotropic) using the 3D BEST HNCO (JCOH and JCC) or BEST-HNCOCA (JCAHA) experiments^69,70,71. The one-bond ¹H−¹⁵N coupling was determined by recording two ¹⁵N-HSQC sub-spectra, in-phase (IP) and anti-phase (AP). For the measurement of the ¹H-¹³CO coupling constants a BEST HNCO-JCOH experiment was used with an introduced DIPSAP filter. Such J-mismatch compensated DIPSAP spin-state filter offers an attractive approach for accurate measurement of small spin–spin coupling constants⁷². For that, three separate experiments were recorded with different filter lengths (2τ = 1/J) for each anisotropic and isotropic media, where the sub-spectra associated to the separated spin states (two in phase and one anti-phase) are combined using a linear relation k (IP) + (k − 1) (IP) ± (AP) with k = 0.73, the theoretical optimized scaling factor. The spectra were recorded with 144 × 104 × 1,536 complex points in the ¹³C(t₁)/¹⁵N (t₂)/¹H (t₃) dimensions, respectively, and with the spectral widths set to 15,244 Hz (¹H), 2,070 Hz (¹⁵N) and 1,510 Hz (¹³C) for the HNCO-JCOH. For the HNCO-JCC and HNCOCA-JCAHA 256 × 200 × 1,536 complex points were acquired in the ¹³C(t₁)/¹⁵N (t₂)/¹H (t₃) dimensions, with spectral widths of 15,244 Hz (¹H), 1,900 Hz (¹⁵N) and 1,214 Hz/5050 Hz (¹³C). The recycle delay was set to 200 ms, the acquisition time to 100 ms with 16 scans per increment, and the data was acquired in the non-uniform sampling format (2246 points for HNCO-JCOH and 7680 for the HNCO-JCC/HNCOCA-JCAHA experiments were sampled using the schedule generator from the web portal nus@HMS (http://gwagner.med.harvard.edu/intranet/hmsIST/). The time domain data was converted into the NMRPipe⁶³ format and reconstructed using the sparse multidimensional iterative lineshape-enhanced method (SMILE)⁷³. Coupling constants were obtained from line splitting in the ¹³C or ¹⁵N dimension obtained with CCPN analysis software⁶⁴.

¹⁹F NMR experiments were recorded on a 500 MHz Bruker Avance III spectrometer equipped with a TCI cryoprobe at 298 K (unless otherwise indicated) using a 350 ms acquisition time and 1.5–3 s recycle delay as previously described⁶. We used an amber-suppression strategy to incorporate the unnatural amino acid tfmF, as previously described⁶. Multiple experiments were recorded in succession to monitor sample integrity over time also as previously described⁶. Data were processed using NMRPipe⁶³. Spectra were baseline corrected, peaks were fit to Lorentzian functions and errors of the linewidths and integrals (that is, populations) were estimated using bootstrapping (200 iterations, calculating the standard error of the mean), or from the spectral noise for states whose resonance was not detectable, in MATLAB⁶. ¹⁹F-translational diffusion experiments were performed as previously described⁶.

Thermodynamic parameters of folding (ΔH, ΔS and ΔC_p) were obtained from a nonlinear fit to a modified Gibbs–Helmholtz equation, assuming ΔC_p remains constant across the experimental temperature range:

K_eq is the equilibrium constant, T is the temperature in Kelvin and T₀ is the standard temperature (298 K). We also fitted the data to the linear van’t Hoff equation (assuming ΔC_p = 0).

The Scipy package with optimize.curve_fit function was used to perform the fits⁷⁴ and errors were estimated as one s.d. from the diagonal elements of the parameter covariance matrix. All parameters (ΔH, ΔS, ΔC_p) generally showed strong correlations with each other (r ≥ 0.8), and thus, their uncertainties correlate also. These parameter correlations are expected⁷⁵. The magnitudes of ΔH and −TΔS are also expected to correlate because we study the temperature dependence of folding in a range where ΔG of folding is close to 0.

Folding free energies were calculated from the experimental populations using ΔG = −RTln(K). The folding free energy of the FLN5+67 wild-type RNC, where no unfolded state is observable, was estimated on the basis of two destabilizing mutants FLN5(V664A/F665A) and FLN5(V707A). The stability of these mutants was measured using ¹⁹F NMR on and off the ribosome. The FLN5+67 wild-type folding energy (ΔG_N-U) was then calculated as the average from the V664A/F665A and V707A mutants using ΔG_WT,+67 = ΔG_mut,+67 − ΔΔG_mut-WT,iso, where ΔΔG_mut-WT,iso is the experimentally measured destabilization in isolation. Given that at FLN5+34, both mutants show a weaker destabilization than in isolation, we reasoned that this estimate of the FLN5+67 wild-type folding free energy is its lower bound (most negative).

¹⁹F transverse relaxation rate (R₂) measurements were recorded using a Hahn-echo sequence and acquired as pseudo-2D experiments with relaxation delays of 0.1 to 200 ms. Data were processed using NMRPipe and analysed using MATLAB. Data were fit to lineshapes and R₂ was obtained by fitting the integrals to single exponential functions. We also orthogonally determined R₂ from linewidth measurements of spectra acquired by 1D ¹⁹F pulse-acquire experiments, which showed excellent correlations. The lineshape-derived R₂ values also showed good correlation with previously determined rotational correlation times²⁵ (τ_C). We additionally determined the S²Ï„_C of FLN5 in 60% glycerol at 278 K by measurements of triple quantum build-up and single quantum relaxation as previously described^76,77. Thus, our R₂ values can be used to determine rotational correlation times (τ_C,exp). The obtained τ_C,exp was used to estimate the bound population as \(\frac{{\tau }_{{\rm{C}},\exp }-{\tau }_{{\rm{C,iso}}}}{{\tau }_{{\rm{C,bound}}}-\,{\tau }_{{\rm{C,iso}}}}\), where τ_C,iso is the rotational correlation time of the isolated protein²⁵ (7.7 ns at 298 K) and τ_C,bound is the expected rotational correlation time of the bound state. τ_C,bound is taken as the rotational correlation time of the ribosome itself (~3,000 ns at 298 K) for a fully rigid bound state. From the bound populations (p_B), the resulting change in the folding free energies of the intermediates was calculated as ΔΔG_I-U,RNC-iso = RT(ln(1 − p_B)). We report the estimate for a fully rigid bound state in the main text (\({S}_{{\rm{bound}}}^{2}=1.0\)) but note that even one order of magnitude more flexibility in the bound state (\({S}_{{\rm{bound}}}^{2}=0.1\)) only accounts for up to 1.1 ± 0.6 and 0.4 ± 0.1 kcal mol^−1 of stabilization for I1 and I2 on the ribosome at FLN5+47, respectively. These estimates still cannot account for the >4 kcal mol^−1 of intermediate stabilization observed on the ribosome⁶.

All NMR experiments of RNCs are recorded and continuously interleaved with a series of 1D ¹H/19 F spectra and ¹H,¹⁵N/¹⁹F diffusion measurements^6,7,61,78. These provide the most sensitive means to assess changes in the sample, and when alterations in signal intensities or linewidths (that is, transverse relaxation rates), chemical shifts or translational diffusion measurements of the nascent chain are observed, data acquisition is halted. Only data corresponding to intact RNCs are summed together and subjected to a final round of analysis. Where signal-to-noise remains low, datasets from multiple samples are compared to ensure identical spectra, before summation together into a single NMR spectrum. Biochemical assays provide an orthogonal means to assess nascent chain attachment to the ribosome. Identical samples incubated in parallel with NMR samples are analysed by SDS–PAGE (under low pH conditions⁶¹) and detected with nascent-chain-specific antibodies. Ribosome-bound species migrate with an addition ~17-kDa band-shift relative to released nascent chains due to the presence of the tRNA covalently linked to the nascent chain. Combined with time-resolved NMR measurements, these analyses confirm that the reported NMR resonances originate exclusively from intact RNCs.

Mass spectrometry

FLN5 A₃A₃ was buffer exchanged into 100 mM (NH₄)₂CO₃ at pH 6.8 (using formic acid for pH adjustment). Analyses were run on the Agilent 6510 QTOF LC–MS system at the UCL Chemistry Mass Spectrometry Facility. Samples contained ~10–20 μM of protein and 10 μl were injected onto a liquid chromatography column (PLRP-S, 1,000 Å, 8 μm, 150 mm × 2.1 mm, maintained at 60 °C). The liquid chromatography was run using water with 0.1% formic acid as mobile phase A and acetonitrile with 0.1% formic acid as phase B with a gradient elution and a flow rate of 0.3 ml min^−1. ESI mass spectra were continuously acquired. The data were processed to zero charge mass spectra with the MassHunter software, utilizing the maximum entropy deconvolution algorithm.

Small-angle X-ray scattering

We measured SAXS of an isolated FLN5 A₃A₃ C747V sample in Tico buffer supplemented with 1% (w/v) glycerol. Data collection was performed at the DIAMOND B21 beamline (UK)⁷⁹ with a beam wavelength of 0.9408 Å, flux of 4 × 10¹²â€‰photons s^−1 and an EigerX 4 M (Dectris) detector distanced at 3.712 m from the sample. A capillary with a 1.5 mm diameter kept at 283 K was used for data acquisition. We acquired SAXS data at multiple protein concentrations (5.5, 2.75, 1.38, 0.69, 0.34 and 0.17 mg ml^−1) to assess whether the sample exhibited signs of aggregation or interparticle interference. At 5 mg ml^−1, we observed weak signs of interparticle interference in the low q region of the scattering profile, which is also reflected in the R_g obtained by Guinier analysis (using the autorg tool from ATSAS⁸⁰; Supplementary Table 1). Data were recorded as a series of frames, non-defective frames were averaged, and buffer subtracted with PRIMUS⁸⁰. Size-exclusion chromatography–SAXS (SEC–SAXS) experiments were additionally performed in Tico buffer with 1% (w/v) glycerol using a KW402.5 (Shodex) column to confirm the monodispersity of the sample. We chose the 2.75 mg ml^−1 dataset as the final dataset to compare with our molecular dynamics simulations. This dataset exhibited the highest signal to noise ratio and did not show signs of interparticle interference, and accordingly, the R_g obtained from the 2.75 mg ml^−1 dataset is consistent with the value obtained from lower concentrations and the main SEC–SAXS peak (Supplementary Table 1).

Circular dichroism spectroscopy

The circular dichroism (CD) spectrum of isolated FLN5 A₃A₃ V747 was acquired in 10 mM Na₂HPO₄ pH 7.5 at 283 K. A Chirascan-plus CD spectrometer (Applied Photophysics), a protein concentration of 44 μM and a cuvette with a 0.5 cm pathlength were used.

HRAS refolding experiments

HRAS refolding experiments were performed with the HRAS G-domain (residues 1–166). The protein was unfolded overnight at 298 K in Tico buffer with 2 mM β-mercaptoethanol, 8 M urea and protein concentration of 15 μM. The protein was then refolded by rapidly diluting into Tico buffer (supplemented with 2 mM β-mercaptoethanol and 50 μM GDP) to reach final urea and protein concentrations of 0.94 M and 1.76 μM, respectively, and allowed to incubate at 298 K for 24 h. For NMR analyses of refolded samples, we prepared 18 μM of refolded protein with the same urea concentrations and dilutions.

We assayed the functional/activity state of HRAS using GDP/GTP nucleotide exchange (‘activity’) assay⁸¹ with fluorescently labelled GTP that exhibits higher fluorescence when bound to HRAS than free in solution (BODIPY FL GTP, ThermoFisher). 0.4 μM of HRAS, 0.01 μM of BODIPY GTP and 1 μM of SOS_cat (the catalytic domain of Son of sevenless) were incubated at room temperature and the maximum (plateau) fluorescence recorded and normalized by the signal of the buffer (signal/noise ratio). SOS_cat was produced as previously described⁸². Fluorescence measurements were performed using the CLARIOstar microplate reader (BMG Labtech) with excitation and emission wavelengths set to 488 and 514 nm, respectively.

The proteolytic stability of HRAS was assayed with thermolysin at a concentration of 0.05 mg ml^−1 incubated with HRAS samples over the course of 5 h in vitro and 9 h in rabbit reticulocyte lysate (RRL, TNT coupled reticulocyte lysate, Promega). Reactions were quenched with 23 mM EDTA. Timepoints were analysed by western blot analysis using a pan-RAS polyclonal antibody (ThermoFisher, 1:1,000 dilution), utilizing an anti-rabbit IgG horseradish peroxidase-linked secondary antibody (Cell Signaling Technology, 1:1,000 dilution). Densitometry analyses were performed with ImageJ⁸³. For the RRL experiments, refolding reactions were performed in RRL for 24 h at 298 K and a final HRAS concentration of 1.6 μM followed by pulse proteolysis and we quantified the relative band intensities (refolded/control) for each time point to account for increased background on the western blot during the proteolysis reaction.

Molecular dynamics simulations

We used the FLN5 A₃A₃ C747V sequence for all simulations. A reliability and reproducibility checklist is provided in Supplementary Table 8. GROMACS (version 2021)⁸⁴ was used for all all-atom molecular dynamics simulations in explicit solvent. We employed the Charmm36m force field in combination with the CHARMM TIP3P water model (C36m) and the CHARMM TIP3P water model with an increased water hydrogen LJ well-depth (denoted here as C36m+W)⁸⁵. We also used the a99sb-disp force field together the a99sb-disp TIP4P-D water model⁸⁶. Default protonation states were used in all cases. Starting from a random extended conformation, for all force field combinations the system was solvated in a dodecahedron box with 151,135 water molecules and 12 mM MgCl₂ (resulting in 455,116 atoms and an initial box volume of 4,688 nm³). Systems were then energy minimized using the steepest-decent algorithm. For the following dynamics simulations, we used the LINCS algorithm⁸⁷ to constrain all bonds connected to hydrogen and a timestep of 2 fs using the leap-frog algorithm for integration. Nonbonded interactions were calculated with a cut-off at 1.2 nm (including a switching function at 1.0 nm for van der Waals interactions) and the particle mesh Ewald (PME) method⁸⁸ was used for long-range electrostatic calculations. We then equilibrated the systems in two phases. First, we performed a 500 ps equilibration simulation in the NVT ensemble with position restraints on all protein heavy atoms. The temperature was kept at 283 K using the velocity rescaling algorithm⁸⁹ and a time constant of 0.1 ps. Next, we further equilibrated the systems for 500 ps in the NPT ensemble at 283 K and a pressure of 1 bar with a compressibility of 4.5 × 10^−5 bar^−1 using the Berendsen barostat⁹⁰. Following equilibration, we relaxed our initial structure for 100 ns at 283 K without any position restraints using the Parrinello–Rahman algorithm⁹¹ and then picked five structures from this simulation for production simulations. We ran a total of 5× 2 μs (with different initial coordinates and velocities) yielding a total of 10 μs of sampling per force field. For the C36m+W combination we ran an additional 5× 2 μs starting from 5 new starting structures yielding 20 μs in total.

We also generated a prior ensemble with a physics-based coarse-grained (C-alpha) model. We generated the C-alpha model template from the FLN5 crystal structure using SMOG 2.3⁹², where all bonded terms have a global energy minimum at the values taken in the crystal structure⁹³. Nonbonded van der Waals interactions were modelled using a 10–12 Lennard-Jones potential with σ and λ parameters described in the M1 parameter determined by Tesei et al.⁹⁴ (equation (S3)). We used the arithmetic mean of two residues to determine σ and λ. Electrostatic interactions were modelled using the Debye–Hückel theory with parameters described previously⁷. Interactions between Cα beads separated by less than four residues were excluded. We ran initial simulations at a range of reduced temperatures to determine the effect on the average compactness and ran final simulations at a reduced temperature of 1.247 (150 K in GROMACS) as we did not observe a significant increase in average R_g beyond this temperature. Simulations were run for a total of 3 × 10⁹ steps with GROMACS (v2018.3).

After simulations, the coarse-grained ensemble was backmapped to an all-atom structure using PULCHRA (v3.06)⁹⁵.

RNC simulations were parameterized using the C36m+W force field/water model combination^85,96. We modelled the ribosome using the structure PDB 4YBB⁹⁷ as a template, which we previously refined against a cryo-EM map containing an FLN5 RNC⁹⁸. As in our previous work, we only retained ribosome atoms around the nascent chain exit tunnel and accessible surface outside the vestibule⁷. The FLN6 linker and SecM sequence were initially modelled using a cryo-EM map of a FLN5+47 RNC (Mitropoulou et al., manuscript in preparation). The rest of the nascent chain (MHHHHHAS N-terminal tag and FLN5) was then built using PyMol version 2.3 (The PyMol Molecular Graphics System, Schrödinger) and we generated a random initial starting structure with a short simulation using a structure-based force field, SMOG2.3⁹², without native contacts. The FLN5+31 A₃A₃ RNC (containing the C747V mutation) complex was then centred in a dodecahedral box, solvated using 1,030,527 water molecules and neutralized with 706 Mg²⁺ ions, resulting in a final system size of 3,163,127 atoms. The initial box volume was 32,117 nm³. The large box size was necessary to accommodate the highly expanded unfolded state. We then used the same cut-offs and simulation methods as for the isolated protein. We initially also ran a 500 ps equilibration simulation in the NVT ensemble using position restraints on all heavy atoms using a force constant of 1,000 kJ mol^−1 nm² in along the x, y and z axes. We used a temperature of 283 K, which was held constant using the velocity rescaling algorithm⁸⁹ and a time constant of 0.1 ps. Then, we ran a 500 ps equilibration simulation in the NPT ensemble at 283 K using the same position restraints. The pressure was kept at 1 bar with a compressibility of 4.5 × 10^−5 bar^−1 using the Berendsen barostat⁹⁰. The position restraints for all nascent chain atoms (except the terminal residue at the PTC in the ribosome) were then removed, while the ribosome atoms kept being position restrained. In this setup, we ran a 1 ns equilibration simulation at 283 K and 1 bar, using the Parrinello–Rahman algorithm⁹¹. All production simulations were performed using position restraints for the ribosome atoms and C-terminal nascent chain residue at the PTC. Using the equilibrated configuration, we then ran two simulations of ~100 ns to picked ten starting structures for production simulations. Then, ten production simulations of 1.5 μs each (15 μs) were initiated from these different starting structures using random initial velocities. Before the production simulation, each structure was re-equilibrated at 283 K and 1 bar with a 500 ps NVT and 500 ps NPT simulation.

Lastly, to compare our C36m+W simulations with a model that only considers steric exclusion as a nonbonded interaction, we also ran simulations of a simple all-atom model, based on a structure-based model template⁹². We used the FLN5 crystal structure to define the energy minima of all bond and dihedral angles and removed all native contacts. Simulations of isolated and ribosome-bound FLN5 A₃A₃ were run for 1 × 10⁹ steps and 100,000 frames were harvested for analysis. This ensemble was used to compare the expansion of the ensemble, ribosome interactions and conformational entropy with the C36m+W simulations.

Calculation of PREs

The transverse PRE rates of backbone amide groups, Γ₂, were back-calculated from the ensembles using the Solomon–Bloembergen equation^99,100

where μ₀ is the permeability of space, γ_H is the gyromagnetic ratio of the proton, g_e is the electron g-factor, γ_B is the Bohr magneton, S is the proton nuclear spin and J(ω₀) is the generalized spectral density function. For flexible spin labels attached via rotatable bonds the spectral density can be expressed as in equation (S5)¹⁰¹.

where \(\langle {r}^{-6}\rangle \) is the average of the electron–hydrogen distance (r) distribution, S² is the generalized order parameter for the electron–hydrogen interaction vector, τ_C is the correlation time defined in terms of the rotational correlation time of the protein (τ_r) and the electron spin relaxation time (τ_s):

τ_t is the total correlation time defined as:

τ_i is the internal correlation time of the spin label. Since for nitroxide labels electron spin relaxation occurs on a much slower timescale than rotational tumbling^101,102, τ_C can be approximated to τ_r such that expression for τ_t simplifies to

Given that τ_C is not known a priori, we iteratively scanned τ_C values in the range of 1 to 15 ns to find a value for which optimal agreement with the experimental data is achieved (as judged by the reduced χ²)^94,103. The spin label correlation time¹⁰⁴, τ_i was set to 500 ps, in agreement with molecular dynamics simulations¹⁰⁵ and electron spin resonance measurement¹⁰⁶.

The generalized order parameter S² for the electron–hydrogen interaction vector can be decomposed into its radial and angular components¹⁰⁷:

where the individual components are defined as

and \({Y}_{2}^{m}\) are the second order spherical harmonics and Ω^mol are the Euler angles in the frame. A weighted ensemble average of S² can be calculated by taking a weighted ensemble average of the individual radial and angular components.

A previously published rotamer library containing 216 MTSL rotamers¹⁰⁸ was used to explicitly model the flexibility of the spin label, similar to other existing methods^109,110. The rotamer library was aligned to all employed labelling sites for each conformer using the backbone atoms of the labelling site and Cys-MTSL moiety. Clashing rotamers were discarded, where a steric clash between the rotamer and the protein was defined using a 2.5 Å cut-off distance. Only backbone and Cβ atoms were considered for the protein, assuming sidechains can rearrange to accommodate the MTSL rotamer¹¹¹. For MTSL, only the sidechain was included (heavy atoms beyond the Cβ atom). Protein frames for which at least one labelling position cannot sterically allow any MTSL rotamers were discarded. The rotamer library was used to calculate a weighted ensemble-averaged Γ₂ over the rotamer ensemble for each protein conformer in the protein ensemble using equations (S3−S11). The protein ensemble average can then be calculated by averaging Γ₂ over the ensemble.

PRE intensity ratios were then calculated from the ensemble-averaged PRE rate, \(\langle {\varGamma }_{2}\rangle \), using

where R₂ is the linewidth in the proton dimension (residue-specific), R_2,MQ is the linewidth in the nitrogen dimension (multiple-quantum term) and Δ is the delay time in the HMQC experiment (5.43 ms). See Supplementary Note 3 for additional details.

For RNCs, we considered that that ribosome tethering may increase the correlation time of the electron–amide interaction vector due to restricted molecular tumbling near the exit tunnel. We therefore calculated an order parameter, \({S}_{{\rm{NC}}}^{2}\), which quantifies the motion of the electron-interaction vector over the entire nascent chain conformer ensemble (\({S}_{{\rm{NC}}}^{2}\) is distinct from the order parameter S² that quantifies the motion of the MTSL rotamer library attached to a labelling site for a specific protein conformer; equation (S9)). S² is given by

where \({S}_{{\rm{NC}},{\rm{angular}}}^{2}\) and \({S}_{{\rm{NC}},{\rm{radial}}}^{2}\) are given by

and \({Y}_{2}^{{\rm{m}}}\) are the second order spherical harmonics and Ω^mol are the Euler angles in the frame. We approximated the position of the free electron with the Cα atom of the labelling site in this case. A \({S}_{{\rm{NC}}}^{2}\) value of 0 indicates that the vector tumbles completely independent of the ribosome and that the correlation time of the electron–amide vector is the same as for the isolated protein, τ_C,iso. A \({S}_{{\rm{NC}}}^{2}\) value of 1 means that the vector tumbles with the same rotational correlation time as the ribosome (τ_r,70S = 3.3 μs per cP, as determined by fluorescence depolarization¹¹², and τ_r,70S = 4.3 μs at 283 K in H₂O). The effective correlation time, τ_C,eff, of each amide-electron vector is given by

We used a value of 3 ns for τ_C,iso, which was the optimal value determined for isolated FLN5 A₃A₃. Generally, τ_C (equation (S6)) is approximated as τ_C ≈ τ_r because the electron spin relaxation time, τ_s, occurs on a much slower timescale. In fact, measurements of the spin relaxation time of nitroxides have been measured to be on a timescale from hundreds of nanoseconds to several microseconds^113,114,115. The calculated values of τ_C,eff are predominantly below 100 ns except for labelling sites C744, uL23 G90C and uL24 N53C, where values of up to ~250 ns are observed (Supplementary Tables 5 and 6). Thus, we still expect τ_C to be dominated by τ_r and make use of the τ_C ≈ τ approximation.

Finally, reference PRE profiles for a fully extended peptide were calculated from a linear polyalanine chain using a τ_C of 5 ns and R_2,H/R_2,MQ of 100 Hz.

Bayesian inference reweighting

We performed ensemble refinement by reweighting the molecular dynamics-derived ensembles against the experimentally deduced Γ₂ rates using the Bayesian Inference of Ensembles (BioEn) software and method described in the corresponding paper^116,117. These calculations were performed using in-house scripts of the software with the modification to incorporate upper and lower bound restraints in addition to regular restraints with gaussian errors. To this end, these inequality restraints were treated as normal gaussian restraints but subjected to a conditional statement. Lower bound restraints (Γ₂ > 64.5 s^−1 for isolated FLN5 A₃A₃; Γ₂ > 96.0 s^−1 for the RNCs) were applied only if the back-calculated Γ₂ was below the lower bound value. Similarly, upper bound restraints (Γ₂ < 2.2 s^−1 for isolated FLN5 A₃A₃; Γ₂ < 3.7 s^−1 for the RNCs) were applied only if the back-calculated average was above the upper bound. This effectively allows the back-calculated value to vary freely above the lower bound and below the upper bound but imposes a penalty if the inequality condition is not met. The errors of the lower and upper bound values were taken as the combined relative error of that datapoint (that is, the intensity ratio).

As described by Köfinger et al.¹¹⁷, the reweighting optimization problem can be efficiently solved by minimizing the negative log-posterior function (L).

θ expresses the confidence in the initial ensemble, N is the ensemble size, M is the number of experimental restraints, w_α is the vector of weights for the conformers in the ensemble, \({y}_{i}^{\alpha }\) is the back-calculated experimental value i, Y_i is the experimental restraint i, σ_i is the uncertainty of experimental restraint i, and S_KL is the Kullback–Leibler divergence defined as

\({w}_{\alpha }^{0}\) is the vector of initial weights (which were uniform). We used the log-weights method to minimize the negative log-posterior¹¹⁷ and performed reweighting calculations for a range of θ values, as the optimal value of θ cannot be known a priori. Therefore, we conduct L-curve analysis^117,118 by plotting S_KL (entropy) on the x axis and the goodness of fit, quantified by the reduced χ² value, on the y axis. The reduced χ² was calculated against the experimental intensity ratios (I_para/I_dia). This is an effective method to prevent overfitting and introducing a minimal amount of bias into the prior ensemble^117,119. After reweighting, we also calculated the effective fraction of frames contributing to the ensemble average¹¹⁹ as an indication of the extent of fitting.

For RNCs, we used the same approach with an additional modification. Since the PRE depends on τ_C,eff and \({S}_{{\rm{NC}}}^{2}\) which are a function of the weights of individual structures in the ensemble, this consequently leads to changes in τ_C,eff and \({S}_{{\rm{NC}}}^{2}\) when reweighting is performed. Therefore, the conformer-specific PRE values that were used for reweighting are not the same anymore after reweighting. To account for this, we performed 20 iterative rounds of reweighting where each additional round receives input weights and τ_C,eff from the previous round. We found that this leads to convergence of the weights and conformer-specific PREs.

We found that for the ribosomal labelling sites, uL23 G90C and uL24 N53C, the reweighting results are sensitive to the specific ribosome structure used to fit the MTSL rotamer library to, since small variations in the local structure of the labelling site can lead to different rotamer distributions. We tested two different rotamer distributions for the ribosomal labelling sites (Extended Data Fig. 5a), finding that one of them (referred to as R2) gives better agreement with the intermolecular PRE data after reweighting and fits better into the expected density or MTSL rotamers when rotamers are fitted to ten high-resolution ribosome structures (Extended Data Fig. 5a). The R2 rotamer distribution is more representative of the expected variation from structural changes in the labelling sites and was therefore used for our final reweighting calculations.

Calculation of RDCs, R
_h and chemical shifts from molecular dynamics simulations

To back-calculate the R_h from static structures we used an approximate relationship between R_g and R_h values¹²⁰, the latter being calculated from the programme HYDROPRO¹²¹. Thus, we calculated the R_g from Cα atoms using MDAnalysis¹²² and then converted it to R_h using

N is the number of amino acids, α₁ takes a value of 0.216 Å^−1, α₂ takes a value of 4.06 Å, and α₃ has a value of 0.821. The estimated value of R_h (relative to the HYDROPRO calculation) has an average relative uncertainty¹²⁰ of 3%. HYDROPRO itself has a relative uncertainty of ±4% with respect to experimental values¹²¹. Therefore, we treat the back-calculated ensemble-average R_h with a total relative uncertainty of ±5%. The ensemble average was calculated as previously described by Ahmed et al. for back-calculation of PFG-NMR derived values¹²³ of R_h

Chemical shifts were calculated using the SHIFTX2 software¹²⁴ and RDCs were calculated using the global alignment prediction method implemented in PALES¹²⁵. We then scaled the magnitude (that is, the extent of alignment) of the calculated RDCs by a global factor to optimize the Q-factor for each ensemble.

Calculation of SAXS profiles from molecular dynamics simulations

We used Pepsi-SAXS¹²⁶ to compute the theoretical scattering profiles of each conformer in the molecular dynamics ensembles. We treated the contrast of the hydration layer (δ_p) and the effective atomic radius (r₀) as global parameters and used values of 3.34 e^− nm^−3 and 1.025 × r_m (r_m = average atomic radius of the protein) in line with previous work that showed these parameters to well suited for flexible proteins¹²⁷. The constant background and scale factor were also fitted globally using least-squares regression^103,127. The goodness of fit was assessed using the reduced χ² metric, where n is the number of datapoints, q is the scattering angle, \({I}_{q}^{{\rm{calc}}}\) and \({I}_{{q}}^{\text{exp}}\) are the calculated and experimental scattering intensities, respectively, and σ_q is the experimental error:

Structural analysis

The Python package MDAnalysis¹²² and MDTraj¹²⁸ were used for general analysis of the ensembles involving atomic coordinates. For native contact analysis, we calculated the fraction of native contacts (relative to the native FLN5 crystal structure) as¹²⁹

where r_i,j and \({r}_{i,j}^{0}\) are the distances between atoms i and j in frame X and the template structure, respectively, β modulates the smoothness of the switching function (default value 5 Å^−1 used) and λ is a factor allowing for fluctuations of the contact distance (default value 1.8 used).

Asphericity was calculated using MDAnalysis¹²² as defined by Dima and Thirumalai¹³⁰:

\(\bar{\lambda }\) represents the mean eigenvalue obtained from the inertia tensor \(\bar{\lambda }=\,\frac{{\rm{trT}}}{3}\).

For the intrachain contact analysis, we defined contacts between Cα–Cα distances of less than 10 Å. The contact features qualitatively were unchanged when using lower cut-off values or when calculating contacts between all heavy atoms. Secondary structure populations were calculated using DSSP¹³¹ implemented in MDTraj. The SASA was calculated using GROMACS⁸⁴. Clustering was also performed in GROMACS using the GROMOS algorithm¹³² and Cα RMSD cut-offs in the range of 1.2–1.8 nm.

Error analysis from ensembles

Errors from the molecular dynamics ensembles were estimated using a block analysis of the full concatenated ensembles (composed of multiple statistically independent simulations). We performed block analysis for the concatenated ensembles to verify that the estimate of the standard error of the mean (s.e.m.) plateaus/fluctuates at block sizes larger than blocks corresponding to the individual trajectories. The final block size was chosen either in the plateau region of block analysis plots or corresponding to the blocks of the statistically independent simulation (10 independent simulations were run for the isolated and RNC systems with the C36m+W force field, and thus 10 blocks were chosen for block analysis and error estimation). The error after reweighting with PRE-NMR data was calculated the same way using a weighted standard error, where blocks are weighted according to the weights obtained from reweighting with PRE-NMR data. Exemplar block analysis plots are shown in Supplementary Fig. 8.

Energetic analyses from structural ensembles

The conformational entropy was calculated as defined by Baxa et al.¹³³. Proline, glycine and alanine entropies were calculated from the backbone probability distribution P_i(Φ,Ψ). Residues with a maximum of two sidechain torsion angles, X_n, the entropy was calculated from the probability distribution P_i(Φ,Ψ,X₁,X₂), while residues with more sidechain torsion angles was calculated from the sum of entropies obtained using the P_i(Φ,Ψ,X₁), and P_i(X_n), after subtraction of the entropy obtained from P_i(X₁). Entropies were calculated from probability distributions using \(S={-k}_{{\rm{B}}}{\sum }_{i=1}^{n}{P}_{i}{\rm{ln}}({P}_{i})\). We used a block analysis from the pooled ensembles (i.e., all individual trajectories concatenated together) to check that the entropy difference between on and off the ribosome is robust with respect to sampling by calculating entropy changes with increasing amounts of total sampling (from the 15 and 20μs of concatenated sampling for the RNC and isolated protein, respectively). The errors were then also estimated from the same sampling/block sizes up to 7.5 μs of molecular dynamics sampling. This is because the estimate of entropy differences trend increases up to total sampling times of 7.5 μs (Extended Data Fig. 6g).

The energetic contributions due to changes in solvation were estimated based on empirical relationships between changes in the polar and apolar accessible surface area^75,134 (ΔASA_polar and ΔASA_apolar). The apolar and polar surface area of the protein were defined based on the atomic partial charges in the C36m force field⁸⁵. Atoms with an absolute charge of less than or equal to 0.3 were defined as apolar. The change in heat capacity of hydration is related to these quantities by

where α and β are 0.34 ± 0.11 and −0.12 ± 0.12 cal mol^−1 K^−1 Å^−2, respectively. We obtained these values as an average and standard deviation of parameters previously reported in the literature as summarized in ref. ¹³⁵ to account for the uncertainty of the parameters in addition to the uncertainty coming from conformational sampling in our simulations. The enthalpy change due to solvation is then obtained from⁷⁵

T is the temperature and γ and δ constants taking on values of −8.44 and 31.4 cal mol^−1 Å^−2, respectively. While we are not aware of alternative parameter sets for the solvation enthalpy (equation (S26)) in the literature, we treated these parameters with a relative uncertainty of 50% to show that even with such high levels of uncertainty our conclusions are not affected. Finally, the solvation entropy and change in free energy are then calculated using

where T_apolar and T_polar are the temperatures at which ΔS_solv,apolar and ΔS_solv,polar are 0 (385 K and 335 K, respectively). Our previous work indicated ribosome solvation changes during coTF is not a major factor in coTF thermodynamics (see Supplementary Note 9), we estimated the above quantities using surface areas calculated excluding the ribosome. We regard these absolute quantities as an estimated upper bound for ΔG_solv because it is likely that folding intermediates and the native state also interact with the ribosome⁶, thus effectively cancelling out any reduction in SASA of the unfolded state due to ribosome interactions. However, in the following section we describe an alternative, more direct approach for the solvation entropy that does not rely on this assumption.

Calculation of solvation entropy changes using the 2PT method

The water and solvation entropy changes were also assessed more directly from molecular dynamics simulations using the two-phase thermodynamic (2PT) method¹³⁶ implemented in the DoSPT code (https://dospt.org/index.php/DoSPT)¹³⁷. For these calculations, we chose five snapshots from our isolated FLN5 A₃A₃ V747 simulations detailed above (that is, with different initial protein conformations and solvent configuration) and use these to initiate short molecular dynamics simulations for entropy calculations. We first re-equilibrated the boxes for 10 ns at the target temperature in the NPT ensemble at 1 bar using the velocity rescaling algorithm⁸⁹ and the Parrinello–Rahman algorithm⁹¹ as detailed above and the velocity Verlet integration algorithm (md-vv in GROMACS⁸⁴). Production simulations were then run in the NVT ensemble at 283 K and 298 K (to assess the effect of temperature on the water entropy calculations) for 20 ps using the md-vv integrator and saving coordinates and velocities for analysis every 4 fs. Control simulations of pure TIP3P (CHARMM TIP3P) water in a cubic box with a box vector length of 5 nm, resulting in 4,055 water molecules. Five independent simulations were performed by first energy minimizing the system using the steepest-decent algorithm. Then, using a 2 fs timestep and thermostat/barostat settings as for the protein and the md-vv integrator we equilibrated the water box first in the NVT ensemble for 1 ns, followed by 1 ns in the NPT ensemble using the Berendson barostat⁹⁰. The water box was then further equilibrated in the NVT ensemble for 1 ns prior to the production simulation in the NVT ensemble for 20 ps, saving coordinates and velocities every 4 fs. These production simulations were also performed at 283 K and 298 K and then used to calculate the molar entropies of pure water at these temperatures with DoSPT.

For water entropy calculations in the protein system, we first analysed the radial distribution function water surrounding the protein molecule using our 15 μs and 20 μs molecular dynamics ensembles of the isolated protein and RNC and the GROMACS rdf functionality⁸⁴ to identify the region of the first two hydration shells that show significantly reduced water dynamics. Using this analysis, we chose a distance cut-off of 3.5 Å between the protein and water centre of mass to define the hydration layer around the protein. With this criterion we then calculated the probability distribution and average number of water molecules in the hydration layer to assess the difference in solvation on and off the ribosome. Water molecules that remain within a defined distance range from the protein during the entire 20 ps production simulation were then selected to calculate the average molar entropy per molecule of water in different environments with DoSPT. The accessible volume for this subsystem was estimated by using the average volume occupied per water molecule in a pure water box under identical conditions multiplied by the number of molecules. To obtain the change in solvation entropy (difference between the RNC and isolated system, ΔS_solv,RNC-iso), we used

where N_diff is the average difference in the number of water molecules in the hydration layer (RNC-iso) and ΔS_solv,water is the entropy difference between water molecules in the hydration layer (0–3.5 Å from the protein) and water molecules in bulk solution (defined here as 36–46 Å from the protein).

Reporting summary

Further information on research design is available in the Nature Portfolio Reporting Summary linked to this article.