We present a coarse-grained simulation model that is capable of simulating the minute-timescale dynamics of protein translocation and membrane integration via the Sec translocon, while retaining sufficient chemical and structural detail to capture many of the sequence-specific interactions that drive these processes. The model includes accurate geometric representations of the ribosome and Sec translocon, obtained directly from experimental structures, and interactions parameterized from nearly 200 μs of residue-based coarse-grained molecular dynamics simulations. A protocol for mapping amino-acid sequences to coarse-grained beads enables the direct simulation of trajectories for the co-translational insertion of arbitrary polypeptide sequences into the Sec translocon. The model reproduces experimentally observed features of membrane protein integration, including the efficiency with which polypeptide domains integrate into the membrane, the variation in integration efficiency upon single amino-acid mutations, and the orientation of transmembrane domains. The central advantage of the model is that it connects sequence-level protein features to biological observables and timescales, enabling direct simulation for the mechanistic analysis of co-translational integration and for the engineering of membrane proteins with enhanced membrane integration efficiency.
Most integral membrane proteins (IMPs) are co-translationally inserted into the membrane during biosynthesis via the Sec translocon, a multiprotein complex [1–4]. In this process, a ribosome docks to the cytosolic opening of the Sec translocon and feeds a nascent polypeptide chain (NC) into the translocon channel. Secretory proteins, or the soluble domains of IMPs, translocate across the lipid membrane by passing through the translocon channel [1,2]. Alternatively, the transmembrane domains (TMDs) of IMPs integrate directly into the lipid membrane via the translocon lateral gate (LG). Integration is facilitated by a conformational change in the channel that separates the two LG helices to create an opening between the channel interior and the hydrophobic core of the membrane [5–7]. The likelihood of integration or translocation of polypeptide segments depends on residue-specific chemical features of the nascent polypeptide chain, such as its hydrophobicity and charge [8–12], but is also governed by the dynamics of protein synthesis on the minute timescale [13,14]. To reach a stable folded structure, IMPs must integrate into the membrane with the correct topology (i. e. , orientation of each TMD with respect to the membrane), which depends sensitively on the properties of both the NC and the translocon itself [3,15]. Even single mutations to an IMP amino-acid sequence can disrupt integration and induce disease phenotypes [16] or decrease protein expression [17–19]; similarly, mutations to the translocon channel can inhibit IMP folding [8,20–23]. The important role for IMPs in cellular functions, such as signal transduction, the transport of nutrients, and cell adhesion, motivates the understanding of the effect of NC and translocon properties on the efficiency of co-translational integration. However, a detailed understanding of this process presents challenges for both theory and experiment due to the long range of timescales (from nanoseconds to minutes) that are involved. Experimental studies have elucidated many aspects of the structure and function of the Sec translocon, although their ability to directly probe the non-equilibrium co-translational integration process is limited. Structural characterization has revealed many of the components of the translocon complex in both eukaryotes [24–28] and prokaryotes [6,7, 29–32], while biophysical assays have investigated the functional effects of NC hydrophobicity [8,9], charges flanking TMDs [10–12], soluble loop length [13,14], and the forces exerted on a NC during translation [4,33,34]. Despite these findings, mechanistic details of the co-translational integration process remain in question [4] because most experiments are limited to probing final protein distributions—such as the fraction of protein in a specific topology [14] or the fraction of protein integrated in the membrane [35]—and do not typically resolve NC dynamics. Atomistic-scale molecular dynamics simulations can be used to probe detailed aspects of co-translational integration, with recent simulations providing insight into the energetics of TMD integration [36,37], the dynamics of water inside the translocon [38], the effect of NC properties on LG opening [39], the dynamics of a NC during the initial stages of translation [40–42], and the dynamics of IMP integration in simplified system representations [43,44]. However, the separation of timescales relevant to co-translational integration poses a significant challenge to conventional simulation methods: notably, ribosomal translation requires seconds to minutes to complete the biosynthesis of typical polypeptides [45–48], while conformational fluctuations of the NC occur on the nanosecond timescale. Currently available simulation approaches either fail to reach the biological timescales of ribosomal translation [38,40,41] or lack sufficient detail to describe detailed features of the NC-translocon interactions and NC conformational dynamics [43,44,49]. The model presented here overcomes these limitations, allowing direct comparison with a broad range of available experiments. In previous work, a highly coarse-grained (CG) model of Sec-facilitated IMP integration was developed in which all system coordinates are projected onto a two-dimensional plane passing through the translocon LG [43]. This 2D-CG model includes an explicit representation of NC translation, translocon LG conformational gating, and a sufficiently simple system description to enable minute-timescale unbiased trajectories. Previous work has demonstrated that the 2D-CG model correctly predicts the distribution of topologies obtained by TMDs as a function of C-terminal soluble loop length [43], the probability of membrane integration as a function of TMD hydrophobicity [43], the effect of charge mutations on the topology of the dual-topology protein EmrE [50], and the effect of sequence modifications on the integration efficiency of the multispanning protein TatC [19]. The 2D-CG model was also used to demonstrate a link between IMP integration efficiency and expression levels for TatC [19], enabling the computational prediction of amino-acid sequence modifications that improve IMP expression. These successes illustrate the potential for using CG methods to capture the essential physics of the co-translational protein translocation and membrane integration processes. However, several shortcomings of the 2D-CG model have been identified. In particular, the ribosome and translocon are modeled without detailed structural features, sequence-specific ribosome and translocon chemical features are not mapped directly to the CG representation, and interactions between the NC and the translocon are independent of NC sequence. These shortcomings limit the ability of the 2D-CG model to investigate phenomena arising from sequence-specific structural and chemical features, such as variations among homologs of the Sec translocon [6,7] or interactions between the NC and translocon [51,52]. In the current work, we describe a refined CG model that enables simulation of the long time- and length-scales that are relevant to co-translational protein integration, while preserving sequence-specific properties of the NC and translocon and capturing the structure of the ribosome-translocon complex. The new 3D-CG model extends the 2D-CG model by providing a realistic three-dimensional representation of the ribosome-translocon complex mapped directly from high-resolution structural data [6,25]. Additionally, the model is parameterized via a bottom-up approach to reproduce sequence-specific NC-translocon interactions, and it includes a protocol for directly mapping any input amino-acid sequence to a simulation representation, enabling simulation of any polypeptide using only the amino-acid sequence as input. The improved 3D-CG model is validated by reproducing experimental measurements of TMD integration efficiency [51] and signal peptide topogenesis [14]. The model further reproduces the “biological hydrophobicity” scale derived by von Heijne and co-workers [51], capturing the effects of single-residue mutations on stop-transfer efficiency. The strong agreement between simulation and experiment indicates that the 3D-CG model produces simulation predictions that can be confirmed by direct experimental analogues. The new model provides a framework for performing mutagenesis studies of the NC and ribosome-translocon complex to obtain a detailed mechanistic understanding of interactions that impact TMD integration and topogenesis, potentially enabling the prediction of IMP sequence modifications with enhanced membrane integration efficiency and stability.
Fig 1A presents the components of the 3D-CG model compared to an image of the ribosome-translocon complex obtained from a cryo-EM structure [25]. The SecYEG translocon (grey/green), ribosome (brown), and the NC (cyan/red) are represented with explicit CG beads, while the implicit membrane is drawn as a shaded region. As in the 2D-CG model [43], each CG bead has a diameter of σ = 0. 8 nm, the Kuhn length of a polypeptide chain [43,44], and represents three amino-acid residues; σ sets the length scale for the 3D-CG model. The coordinate system is defined such that the origin is placed at the geometric center of the translocon channel Cα atoms, the implicit membrane spans the x-y plane with its midplane located at z = 0σ, and the axis of the translocon is aligned with the z-axis (Fig 1C). The geometry of the Sec translocon is obtained by mapping all amino-acid residues of the translocon onto CG beads in a ratio of three amino acids to one CG bead, where the CG bead is positioned at the center of mass of the Cα atoms for each consecutive triplet of amino-acid residues in the translocon primary sequence. Triplets of amino acids with a net positive charge are assigned a +1 charge, and triplets of amino acids with a net negative charge are assigned a -1 charge. To determine the net charge of a triplet of amino acids the charges of the amino acids are summed, with arginine and lysine counted as +1, and aspartate and glutamate counted as -1 (see S2 Appendix for further discussion). The translocon is modeled in two distinct conformations, with the LG either closed or open (Fig 1B). CG bead coordinates for both conformations are obtained from residue-based coarse-grained simulations of the Methanocaldococcus jannaschii SecYEG translocon (PDB ID: 1RHZ) [6] (see S2 Appendix). The 3D-CG model of the translocon is oriented such that the y-axis of the simulation coordinate system passes between the helices of the LG when the translocon is in the open conformation (Fig 1C). The geometry of the ribosome is obtained by mapping the ribosome-translocon complex from a recent high-resolution cryo-EM structure (PDB ID: 3J7Q) onto CG beads [25]. Amino-acid residues are mapped onto CG beads in a 3: 1 ratio following the same procedure used for the translocon. Each RNA nucleotide in the ribosome is mapped onto two CG beads; one bead represents the sugar-phosphate backbone, while the other bead represents the nucleobase. This mapping is used to capture the excluded volume and the rigidity of the RNA scaffold and is consistent with previous work on coarse-grained DNA/RNA simulations [53–55]. Each CG bead representing a RNA sugar-phosphate backbone in the ribosome is assigned a -1 charge and each CG bead representing a nucleobase is neutral. Only the portion of the ribosome near the translocon channel is explicitly represented as CG beads in the final simulation system (Fig 1A; additional details are in S2 Appendix). Ribosome CG bead positions are identical for both translocon conformations. To characterize whether the ith NC bead, with position xi = (xi, yi, zi), is located in the implicit membrane region, we define the characteristic function S mem (x i) = [ 1 - S (x i, y i) ] S (z i), (1) which assumes a value of 1 in the membrane and 0 elsewhere. S (x, y) and S (z) are smooth switching functions, S (x, y) = 1 4 1 + tanh x 2 + y 2 + 1. 5 σ 0. 25 σ 1 - tanh x 2 + y 2 - 1. 5 σ 0. 25 σ, (2) and S (z) = 1 4 1 + tanh z + 2 σ 0. 25 σ 1 - tanh z - 2 σ 0. 25 σ, (3) where x 2 + y 2 is the radial distance from the coordinate system origin in the x-y plane. S (x, y) is approximately 1 for the range - 1. 5σ < x 2 + y 2 < 1. 5σ and 0 elsewhere, while S (z) is approximately 1 for the range -2σ < z < 2σ and 0 elsewhere (Fig 1C). Eqs 1–3 are used in Eq 8 to define the solvation of a NC bead. The potential energy function for the 3D-CG model is expressed U (xn, xc; q, g) =Ubond (xn) +Uexcl (xn) +Uelec (xn, xc; q) +Usolv (xn; g) +Uchan (xn, xc; g) +Uribo (xn), (4) where xn indicates the set of NC bead positions, xc indicates the set of channel and ribosome bead positions, q is the set of all bead charges, and g is the set of all NC bead transfer free energies. All interactions in the 3D-CG model are defined using an energy scale given by ϵ = kBT, where the temperature, T, is fixed at 310 K to represent physiological conditions. Bonded interactions between consecutive NC beads are described using the finite extension nonlinear elastic (FENE) potential, U bond (x n) = - 1 2 K 0 R 0 2 ∑ b ∈ Bonds ln 1 - r b 2 R 0 2, (5) where the sum runs over all bonds in the NC, rb is the distance between the NC beads that share bond b, K0 = 5. 833 ϵ/σ2, and R0 = 2σ. Short-ranged excluded volume interactions between pairs of NC beads are modeled using a purely repulsive Lennard-Jones (LJ) potential [56], U excl (x n) = ∑ i, j ∈ NC 4 ϵ i j σ i j r i j 12 - σ i j r i j 6 + ϵ i j, r i j < 2 1 / 6 σ i j 0, r i j ≥ 2 1 / 6 σ i j, (6) where the sum runs over all pairs of NC beads, rij is the distance between NC beads i and j, and ϵij = ϵ, and σij = σ. Electrostatic interactions are described using the Debye-Hückel potential, U elec (x n, x c; q) = ∑ i, j ∈ All l B q i q j ϵ r i j exp - r i j κ, (7) where the sum runs over all pairs of charged beads, lB is the Bjerrum length, qi is the charge of CG bead i in the NC, translocon, or ribosome, and κ is the Debye length. Assuming that electrostatic interactions are screened by physiological salt concentrations [57,58], the electrostatic length scales are approximated by κ = lB = σ. NC bead interactions with the implicit solvent are described using a position-dependent potential, U solv (x n; g) = ∑ i ∈ NC g i S mem (x i), (8) where xi is the position of NC bead i, and gi is the transfer free energy for partitioning NC bead i from water to the membrane. Residue-specific interactions between NC beads and translocon beads are given by U chan (x n, x c; g) = ∑ i ∈ NC [ 1 - S mem (x i) ] U chan aq (x i, x c; g i) + [ S mem (x i) ] U chan mem (x i, x c; g i). (9) Eq 9 smoothly interpolates between NC bead-translocon interactions for which NC bead i is positioned in aqueous solution inside the channel (U chan aq (x i, x c; g i) ) or positioned in the membrane near the channel exterior (U chan mem (x i, x c; g i) ). The exact functional forms of U chan aq (x i, x c; g i) and U chan mem (x i, x c; g i) are described in the section Parameterization of NC-translocon interactions. Interactions between NC beads and ribosome beads are included in the Uchan (xn, xc; g) potential energy term (Eq 9). Contrary to interactions between NC beads and translocon beads, interactions between NC beads and ribosome beads are not bead-type specific; they are described by a repulsive soft-core LJ potential (Eq 17), with ϵij = ϵ and σj = 1. 2σ. To prevent the NC from moving into the part of the ribosome that is not explicitly included in the simulations (see 3D-CG Model Geometry), a repulsive sphere is centered at (-10σ, -0. 5σ, 1. 0σ) (Fig 1C). Repulsive interactions with this sphere are described using U ribo (x n) = ∑ i ∈ NC 4 ϵ σ r i r - 2 σ 12 - σ r i r - 2 σ 6 + ϵ, r i r - 2 σ < 2 1 / 6 σ 0, r i r - 2 σ ≥ 2 1 / 6 σ, (10) where rir is the distance of the NC bead i from the center of the sphere. The time evolution of the NC beads is modeled using overdamped Langevin dynamics with a first-order Euler integrator [59], x n (t + Δ t) = x n (t) - β D ∇ x n U (x n (t), x c (t); q, g) Δ t + 2 D Δ t R (t), (11) where xn (t) are the positions of the NC beads at time t, U (xn (t), xc (t); q, g) is the 3D-CG model potential energy function (Eq 4), β = 1/kBT, D is an isotropic diffusion coefficient, and R (t) is a random number vector drawn from a Gaussian distribution with zero mean and unit variance. The timestep, Δt = 300 ns, permits stable integration of the equations of motion with a diffusion coefficient of D = 253. 0 nm2/s (see S3 Appendix for discussion and Table S2 in S3 Appendix for robustness with respect to timestep). Ribosome CG bead coordinates are fixed throughout the simulations. Translocon CG beads undergo stochastic transitions between fixed configurations associated with the open versus closed lateral gate. NC-dependent conformational gating of the translocon is attempted at every simulation timestep. The probability that the translocon transitions from a closed (x c closed) to open (x c open) conformation, popen (xn; q, g), is p open (x n; q, g) = 1 τ LG exp - β Δ G open (x n; q, g) 1 + exp - β Δ G open (x n; q, g) Δ t, (12) and the probability that the translocon transitions from an open to closed conformation, pclose (xn; q, g), is p close (x n; q, g) = 1 τ LG 1 1 + exp - β Δ G open (x n; q, g) Δ t. (13) The timescale for attempting translocon conformational changes, τLG = 500 ns, is obtained from prior molecular dynamics simulations [39,43]. The total free energy change for switching the translocon from the closed to open conformation, ΔGopen (xn; q, g), is given by Δ G open (x n; q, g) = Δ G empty + U (x n, x c open; q, g) - U (x n, x c closed; q, g), (14) where ΔGempty = 3ϵ is the free energy penalty for opening a closed channel in the absence of a substrate [60], U (x n, x c open; q, g) is the 3D-CG model potential energy function (Eq 4) with the channel in the open configuration, and U (x n, x c closed; q, g) is the 3D-CG model potential energy function (Eq 4) with the channel in the closed configuration. Previous simulations have found the translocon to exhibit both closed and open lateral-gate conformations [39], and the timescale needed to perform this conformational switch is relatively small (500 ns) in comparison to the other timescales modeled in the 3D-CG model [40]. Therefore, as in the 2D-CG model [43], the lateral-gate conformational changes in the 3D-CG model are described in terms of instantaneous switches between the closed and open conformations. If an attempted conformational change is accepted, all bead positions in the translocon are immediately switched to the positions corresponding to the new channel conformation. The equations of motion described by Eqs 11–14 rigorously obey detailed balance. Translation of the NC is modeled by adding CG beads to the C-terminus of the NC during a simulation trajectory. At the initiation of the trajectory, the C-terminal NC bead is fixed at the translation insertion point (Fig 1C). For each simulation timestep in which translation is performed, the C-terminal bead is moved in the +z direction by a distance equal to σΔt/ttrans, where ttrans is the timescale for translating a single CG bead. ttrans is set to 0. 6 seconds to reproduce a translation rate of 5 residues/second [45–48] unless otherwise specified. The C-terminal NC bead is otherwise held fixed, although all interactions between the C-terminal NC bead and other NC beads are included in Eq 4. The translation of the C-terminal bead is completed after a period of ttrans and its dynamics are described using Eq 11 for the remainder of the simulation trajectory. The next CG bead in the NC sequence is then positioned at the translation insertion point and the process is repeated until all NC beads have been translated. For the combined dynamics of the ribosome-translocon-NC system, a series of five steps is iterated at each trajectory timestep: (i) forces acting on each NC bead are calculated, (ii) NC bead positions are time-evolved using Eq 11, (iii) conformational gating of the translocon is attempted (Eqs 12 and 13), (iv) ribosomal translation is performed if all NC beads have not yet been translated, and (v) the simulation is terminated if user-defined conditions are met. Specific protocols for initializing and terminating simulation trajectories are provided for each workflow described in the Results. While the system geometry, 3D-CG model dynamics, and most terms in the 3D-CG model potential energy function (Eq 4) are fully described in the Methods, the functional forms of the NC-translocon interaction potentials, U chan aq (x i, x c; g i) and U chan mem (x i, x c; g i) in Eq 9, have yet to be specified. Here, we describe the protocol for obtaining these potentials, which determine sequence-specific NC bead-translocon interactions. First, we define a protocol for assigning an effective water-membrane transfer free energy, gi, and charge, qi, to a NC bead, based on available experimental data. Second, potentials of mean force (PMFs) for translocating model tripeptide substrates across the translocon channel are calculated using the MARTINI residue-based coarse-grained force field. Finally, sequence-specific NC bead-translocon interactions in the 3D-CG model are parameterized by reproducing the MARTINI PMFs using the 3D-CG potential energy function. The interactions between a general NC bead and the rest of the system is defined by four parameters: gi, qi, λo (gi), and λ i c (g i). These parameters are determined as described in detail in section 3D-CG Model Parameterization. Specifically, the NC bead transfer free energy, gi, is equal to the sum of the transfer free energies of the three amino-acid residues associated with the bead according to the Wimley-White hydrophobicity scale (Table S1 in S3 Appendix). For each residue that does not form secondary structure, gi is increased by 1. 78ϵ, the cost for partitioning a peptide bond that lacks hydrogen bonds. The CG bead charge, qi, is equal to the sum of the charges of the three associated amino-acid residues. The N- and C-terminal CG beads are assigned an additional +1 and -1 charge, respectively, and have 6ϵ added to their transfer free energies to account for the additional charge [66]. The scaling parameters for NC-channel interactions, λo (gi) and λc (gi), are determined from gi using the piecewise-linear interpolation scheme shown in Fig 2C. Fig 3 demonstrates the mapping procedure for an example amino-acid sequence. To start a 3D-CG simulation, both an input amino acid-sequence and a secondary structure assignment for this sequence must be provided. For the membrane integration simulations, the secondary structure of the experimental sequence is reported in the UniProt database and is assigned in the model directly from the available information [72]. For simulations of TMD topology, the secondary structure is not available through the UniProt database and is instead assigned using the PSIPRED secondary structure prediction server [73].
TMDs typically contain a large number of hydrophobic residues to improve stability within the lipid membrane [74]. von Heijne and co-workers measured the probability with which a designed segment (H-segment) of the leader peptidase (Lep) protein integrates into the membrane, demonstrating that the translocon is more likely to integrate hydrophobic NC segments [51]. It was found that increasing the hydrophobicity of a poly-alanine H-segment, through mutation of alanine residues to leucine residues, monotonically increased the probability of H-segment membrane integration. Previous simulations using model sequences and the 2D-CG simulation model revealed that this trend is caused by local equilibration of the H-segment across the translocon lateral gate [43]. Reproducing the same assay using the 3D-CG model, with full structural detail and an direct mapping of the NC amino acid sequence, provides a first means to quantitatively validate model predictions. To simulate the H-segment membrane integration assay with the 3D-CG model, the Lep protein sequence is mapped to CG beads following the procedure described in Mapping amino-acid sequence properties to CG beads. Three helical secondary structure elements, including the H-segment are identified via the UniProt database (ID: P00803). Eight 19-residue H-segments are studied. Each H-segment contains between 0 to 7 leucine residues and the remaining H-segment residues are alanine [51]. All trajectories are initialized from configurations in which the two N-terminal TMDs are already translated. To reduce computational cost, simulations are initiated with the second TMD pre-inserted in the lipid membrane (Fig 4A). The simulated sequences are limited to 90 CG beads in length, corresponding to a continuous stretch of amino acids starting from the second TMD (see S2 Appendix for all simulated sequences). Simulations are terminated when all CG beads of the H-segment either diffuse at least 2σ away from the translocon and span the membrane (integration, Fig 4A, S1 Movie) or when all CG beads have translocated to the lumenal side of the membrane (translocation, Fig 4A, S2 Movie). The probability of membrane integration is defined as the fraction of simulation trajectories that terminate by H-segment integration. Fig 4B shows the comparison of the experimental versus the simulated probability of H-segment membrane integration as a function of the number of leucine residues in the H-segment. The results of the experimental assay [51] are plotted in black squares and the shaded region indicates outcomes within 1 kcal/mol of the experimental measurement as determined by a best fit of the apparent free energy of integration via a sigmoidal curve [51]. The calculated results from the 3D-CG model simulations are plotted in red circles. In agreement with the experiments, the 3D-CG model shows that H-segment integration increases with the number of leucines. Although slightly shifted to the right of the experimental curve, the simulation results recover the same sigmoidal dependence of integration on leucine content and are within 1 kcal/mol accuracy of the experiment [51]. These results indicate that the 3D-CG model correctly predicts trends in NC membrane integration using only information about the protein sequence as input. Fig 4C and 4D investigate the issue of mapping from trios of amino-acid residues to a single CG bead. There are three possible CG representations (frameshifts) of the NC sequence that arise from the 3: 1 mapping of amino-acid residues to CG beads as shown in Fig 4C. Since there is no basis for choosing any one frameshift over the other two, each of the possible frameshifts is simulated, and the calculated membrane integration probabilities shown in Fig 4B is the averaged value over all three frameshifts. For each frameshift and for each of the eight H-segment sequences, 100 trajectories are calculated (ranging from 20–3000 s in time) leading to 2,400 total simulations which required a total of 15,520 CPU hours on 2. 6–2. 7 GHz Intel Xeon processors. All CG bead sequences used in the simulations are provided in S2 Dataset. Fig 4D shows the membrane integration probability for the H-segment sequences for each individual frameshift. Results based on individual frameshifts are comparable, with a notable discrepancy for the 7 leucine H-segment in Frame 1 where the particular grouping of amino acids into triplets resulted in an H-segment for which the integration probability is relatively low. This sensitivity to the choice of triplets is addressed by simply averaging the results over all three frameshifts, which is done for the results in Fig 4B. As shown in Fig 4B, experiments and the 3D-CG model simulations both show that increasing the hydrophobicity of a H-segment by mutating alanine residues to leucine residues increases the probability of H-segment membrane integration. von Heijne and co-workers have extended this analysis by determining the effect of all twenty amino acids on the probability of H-segment membrane integration in the context of the Lep construct [51]. Assuming that there is an effective two-state equilibrium between the integration and translocation outcomes, the probability of integration can be converted into an apparent free energy of integration, ΔGapp, defined by [51] Δ G app = - k T ln p (integration) / p (secretion). (18) By mutating the central residue of the H-segment in the same Lep construct used in the section Probability of membrane integration for NC segments of varying hydrophobicity, von Heijne and coworkers measured Δ G app aa, or the single-residue apparent free energy of integration, for all twenty naturally occurring amino-acid residues, thus deriving a “biological hydrophobicity scale” in analogy to other hydrophobicity scales [63]. Calculating the probability of membrane integration of the same set of H-segments with the 3D-CG model provides a means to validate the predicted effect of single amino-acid residue mutations. The simulation procedure for calculating the biological hydrophobicity scale is the same as illustrated in Fig 4A). To determine Δ G app aa for all 20 amino acids, 22 experimentally studied constructs of the mutated Lep sequence are mapped to a CG representation. Results are averaged over all three frameshifts for each of the 22 constructs, requiring a total of 66 CG bead sequences. All CG bead sequences modeled are provided in S2 Dataset. The probability of H-segment membrane integration is calculated from an ensemble of 200 trajectories (ranging from 20–2000 s in time) per sequence, leading to a total of 13,200 simulations which required a total of 77,003 CPU hours on 2. 6–2. 7 GHz Intel Xeon processors. The probability of H-segment membrane integration is converted to a Δ G app aa following the procedure of von Heijne and coworkers described below [51]. The Δ G app aa for alanine and leucine are determined first from a linear fit of ΔGapp for H-segments with 3 to 7 Leucine residues from the simulated membrane integration probability curves (Fig 4B) using Δ G app = n Leu Δ G app Leu - Δ G app Ala + 19 Δ G app Ala. (19) Δ G app aa for alanine and leucine are found to be 0. 13 kcal/mol and -0. 43 kcal/mol respectively. Experimentally determined values for alanine and leucine are 0. 1 kcal/mol and -0. 6 kcal/mol respectively. The difference in Δ G app aa between simulation and experiment for leucine gives rise to the slight rightward shift of the simulated membrane integration probability curve compared to the experiment in Fig 4B. To obtain Δ G app aa for the remaining amino acids, we employ [51] Δ G app aa = Δ G app x [ aa ] x - Δ G app x [ ref ] x + Δ G app ref. (20) ΔGx[aa]x is the apparent free energy of integration for an H-segment construct with the probed amino acid (aa) at the midpoint of the H-segment ΔGx[ref]x is the apparent free energy of integration for the same H-segment where the probed amino acid is replaced by a reference amino acid with a known apparent free energy of integration, Δ G app ref. The reference amino acids employed match those used in Ref. [51] and are specified in S2 Dataset. The H-segment constructs were chosen to have a leucine content such that the probability of membrane insertion for the sequence is nearly 50% to yield maximum sensitivity in the experimental assay [51]. For cysteine and methionine, we added two additional leucines to the simulated H-segment constructs compared to the experimental constructs to yield additional sensitivity in the computation. Fig 5 compares the values of Δ G app aa determined experimentally to the values of Δ G app aa calculated using the 3D-CG model. Each point represents a single amino acid. Points are colored by grouping amino-acid residues as charged (black), polar (red), aromatic (blue), or non-polar (green). The solid line is a linear fit to the data, while the dashed line illustrates a perfect correlation as a guide to the eye. Each Δ G app aa value is calculated from the average of three frameshifts (defined as in Fig 4). The average standard deviation between the frameshift results is 0. 2 kcal/mol, the error bars indicate the standard error of the mean. Individual frameshift values are reported in Table S3 in S3 Appendix. The experimental and 3D-CG simulation scales are highly correlated (r = 0. 88), confirming that the 3D-CG model reproduces trends in Δ G app aa with high fidelity. The data points largely lie above the dashed line, indicating that the 3D-CG simulations slightly overestimate the experimentally observed degree of integration. These results thus indicate that the 3D-CG is capable of reproducing the effect of single-residue mutations in good agreement with available biophysical measurements, although the quantitative agreement with experiments may still be improved via further model refinements. In addition to determining whether NC segments integrate into the membrane as TM domains, the translocon regulates the orientation with which TM segments integrate (Fig 6A) [14,33,75]. In particular, Spiess and co-workers found that an engineered TM signal anchor (H1Δ22) integrates in either the NER/Ccyt (i. e. Type 1) or the Ncyt/CER (i. e. Type 2) topology; it was also found that decreasing the rate of ribosomal translation by adding cycloheximide increases the preference for the Type 2 topology [14]. Furthermore, increasing the length of the soluble loop flanking the C-terminus of the TM segment also increases the probability that the TM segment obtains the Type 2 topology until the probability eventually plateaus for a sufficiently long loop length. Previous work using the 2D-CG model qualitatively captured both these trends and revealed that the mechanistic basis for the kinetic effect is flipping of the NC from the Type 1 topology to the Type 2 topology as a function of time [43]. However, due to the lack of residue-specific interactions in the 2D-CG model, this work employed model sequences. Additionally, due to the simplified geometric representation of the 2D-CG model, it predicted that p (Type 2) plateaus at shorter C-terminal lengths than is observed in the experiments. While the 2D-CG model can provide mechanistic insights [43], quantitative agreement with the experiments is poor compared to the 3D-CG model when directly mapping the amino-acid sequence (Fig S2 of S3 Appendix and corresponding discussion). Here, we test the 3D-CG model for predicting TMD topogenesis and the effect of translation kinetics on topology. The simulation approach for modeling TMD topogenesis is summarized in Fig 6A (see S3 and S4 Movies for example trajectories). The H1Δ22 sequence is mapped to CG beads, and the results are averaged over all three frameshifts. Nine different lengths of the C-terminal soluble loop are mapped directly from the experimental constructs used in [14]. The first 99 residues of the sequence are assumed to be part of helical domains based on secondary structure predictions from the PSIPRED server [73,76]. Simulations are initialized from configurations in which four CG beads are translated and have not yet entered the translocon. Simulations are terminated when CG beads of the TMD have all integrated into the lipid bilayer in either an Type 1 or Type 2 topology and diffuse 10σ away from the translocon. The final TMD topology is determined from the position of the C-terminal CG bead relative to the membrane upon simulation termination (Fig 6B). All simulations are performed with either the default translation rate of 5 residues/second (fast translation) or with a reduced translation rate of 1. 25 residues/second (slow translation) to model the effect of adding cyclohexamide in the experimental assay. 200 trajectories (ranging from 25–1200 s in time) are simulated for each of the three frameshifts and for each of the nine loop lengths and at both translation rates, leading to a total of 10,800 trajectories which required a total of 149,009 CPU hours on 2. 6–2. 7GHz Intel Xeon processors. Fig 6B compares the simulated and experimental results for the probability with which the TMD obtains the Type 2 topology as a function of the length of the C-terminal soluble loop. The results of the experimental assay are plotted on the right for reference. Results for the normal translation rate are in solid black lines, while results for the reduced translation rate are in dashed red lines. The simulation results correctly reproduce the trends observed in the experiments, including the increased probability of the Type 2 topology for longer C-terminal loop lengths and the eventual plateau in the probability of the Type 2 topology at long C-terminal loop lengths. Furthermore, like the experimental results, the CG model predicts a significant shift to greater Type 2 integration upon reducing the rate of ribosomal translation. We present a refined CG model for co-translational membrane protein integration via the Sec translocon that captures the detailed three-dimensional geometry of the ribosome-translocon complex from high-resolution structural data [6,25] and that describes residue-specific interactions between the NC and translocon based on detailed MD simulations. The bottom-up parameterization approach utilized here employs extensive residue-based coarse-grained simulations to inform the model parameters without the need for additional experimental inputs. In this work, the 3D-CG model is applied to calculate the membrane integration efficiency and topology of TMDs, where the only required input is the amino-acid sequence and NC secondary structure. The 3D-CG model captures the experimentally observed [51] sigmoidal dependence of the probability of TMD integration on substrate hydrophobicity. We extend this analysis to study the effect of all twenty amino-acids on the membrane integration probability yielding values of residue-specific TMD membrane integration probabilities in good agreement with the experimentally observed “biological hydrophobicity” scale [51]. These results demonstrate that the 3D-CG model successfully combines factors that are known from previous work to affect TMD integration at the translocon, such as interactions of the nascent chain and the translocon channel interior [37,38,40], the non-equilibrium nature of peptide elongation [37,43], and the sequence context of the TMD [77]. This suggests that the 3D-CG model is well suited for future applications to investigate phenomena such as the experimentally observed position dependence of the biological hydrophobicity scale [35] and the dependence of the observed hydrophobicity values on the amino-acid residues flanking the TMD [77]. The specific interactions between the NC and the translocon, determined as part of this study, already suggest a mechanism by which flanking residues can affect TMD integration; the high barrier for the translocation of charged residues limits translocation, resulting in more integration. Finally, the 3D-CG model accurately describes the experimentally observed effect of translation rate and C-terminal loop length on TMD topogenesis [14]. The 3D representation of the model ensures the correct ribosome-translocon geometry and volume scaling behavior necessary to capture the C-terminal length dependence of TMD topology, an effect not captured in a previous 2D model [40]. The main advantage of the 3D-CG model presented here, compared to previous work, is that it requires few assumptions. NC properties are directly mapped from the underlying amino acid sequence, the ribosome/translocon geometry is mapped from available structural data, and there is no projection onto a two-dimensional subspace. Provided only with an amino acid sequence and a secondary structure assignment, the 3D-CG model obtains striking agreement with experiment, validating the ability of the 3D-CG model to predict key aspects of Sec-fascilitated protein translocation and membrane integration. We additionally emphasize that the 3D-CG model provides a refinable framework for simulating IMP co-translational membrane integration via the Sec translocon. Currently, the bottom-up parameterization strategy uses MARTINI PMFs for four distinct tripeptide substrates as input information. The 3D-CG model parameterization could be refined, either by calculating the PMF of other substrates using the MARTINI force field, by considering the role of changes in substrate protonation state in the channel interior, or by calculating PMFs using an atomistic force field. Furthermore, improved methods for parameterization and uncertainty quantification can be employed to determine parameter sets consistent with the available data [78]. All of these refinements can be made within the current 3D-CG model framework, and they will enable incorporation of additional information and improved quantitative prediction. This framework can also be naturally extended to include additional complexity, such as NC secondary and tertiary structure, other proteins that are part of the Sec translocon complex, and a heterogeneous translation rate. Future studies aimed at the prediction of multispanning IMP topology will guide model development. The 3D-CG model presented here broadens the capability of computer simulation approaches for future studies of the TMD membrane insertion process. In particular, by providing residue-specific NC-translocon interactions the current model enables direct comparison to biophysical measurements of forces between the NC and the translocon due to hydrophobic and electrostatic forces [34,52]. Furthermore, the realistic representation of the structure and interactions enables future mutational studies and comparison of species-specific features of the ribosome-translocon complex to obtain a detailed understanding of key residues that impact TMD integration and topogenesis. The encouraging agreement between 3D-CG model simulation outcome and experiments for single-spanning TMDs displays the capabilities of the 3D-CG framework. It enables the calculation of minute-timescale trajectories in three dimensions, facilitating computational studies that are not possible using existing models with less detail, or existing models that are unable to reach the biologically relevant timescales. The 3D-CG model, with initial model parameters obtained here using a bottom-up strategy, provides a systematically improvable framework for the simulation of co-translational membrane protein integration via the Sec translocon.
Ubiquitous across all kingdoms of life, the Sec translocon is an essential piece of molecular machinery for protein biosynthesis. The translocon is a transmembrane channel that enables the secretion of newly synthesized proteins across the lipid membrane, as well as the integration of protein domains into the membrane interior. Understanding the function and regulation of the translocon is necessary for developing a refined view of early stage protein folding and targeting in the cell. Although computational methods are well suited to elucidating the interactions of the translocon with newly synthesized proteins, conventional simulation techniques are unable to reach the exceedingly long timescales (i. e. , minutes) that are relevant for protein biosynthesis. In this work, we present a novel coarse-grained approach that realistically models the ribosome/translocon/nascent-protein complex, while also allowing for the efficient simulation of biological timescales. The coarse-grained model is parameterized on the basis of extensive molecular dynamics simulations and enables the simulation of any nascent protein with only amino-acid sequence information as input. To validate the model, we perform over 26,000 simulations of protein biosynthesis, enabling direct comparison and demonstrating good agreement with a range of experimental studies describing this minute-timescale process.
chemical compounds
aliphatic amino acids
organic compounds
simulation and modeling
membrane proteins
leucine
amino acids
cellular structures and organelles
thermodynamics
research and analysis methods
lipids
proteins
chemistry
biophysics
cell membranes
ribosomes
free energy
physics
biochemistry
biochemical simulations
cell biology
organic chemistry
integral membrane proteins
biology and life sciences
physical sciences
computational biology
biophysical simulations
Structurally detailed coarse-grained model for Sec-facilitated co-translational protein translocation and membrane integration