Choice of proteins

Tandem Src homology 3 (SH3) domains (Fig 2a) are widely found in signal transduction proteins and they share functions such as mediating protein-protein interactions and regulating ligand binding [25]. Kinetic and thermodynamic properties of native and all the possible circular permutations of SH3 single domain have been well characterized [26]. Two different circular permutant constructs of the sequence are known to fold to a circularly permuted native conformation (PDB accession codes are 1TUC and 1TUD) that is similar to the wild-tpe (WT) protein [26].

With a similar function to the SH3 domains, Src homology 2 (SH2) domains (Fig 2b) are also involved in the mediation of intra- and intermolecular interactions that are important in signal transduction [27]. The SH2 domains are well-known from crystallographic analysis to form metastable domain-swapped dimers [28, 29].

Fibronectin type III (fn3) domains (Fig 2c) are highly abundant in multidomain proteins, and often involved in cell adhesion. We have chosen to study the third fn3 domain of human tenascin (TNfn3), which has been used as a model system to study the mechanical properties of this family. Single-molecule AFM experiments revealed that a small fraction (∼ 4%) of domains in native tenascin (i.e. the full tenascin protein containing both TNfn3 and other fn3 domains) [8], with a similar signature to that observed for I27. Subsequently, misfolding events have been identified in a polyprotein consisting of repeats of TNfn3 only [9]. Interestingly, a structure has been determined for a domain-swapped dimer of TNfn3 involving a small change of the loop between the second and third strand [30].

PDZ domains (Fig 2d) are one of the most common modular protein-interaction domains [31], recognizing specific-sequence motifs that occur at the C-terminus of target proteins or internal motifs that mimic the C-terminus structurally [32]. Naturally occurring circularly permuted PDZ domains have been well studied [33–35], and domain-swapped dimers of PDZ domains have been characterized by NMR spectroscopy [36, 37].

Titin (Fig 2e) is a giant protein spanning the entire muscle sarcomere [38]. The majority of titin’s I-band region functions as a molecular spring which maintains the structural arrangement and extensibility of muscle filaments [39]. The misfolding and aggregation properties of selected tandem Ig-like domains from the I-band of human Titin (I27, I28 and I32) have been extensively studied by FRET experiments [7, 24]. In the earlier work on tandem repeats of I27 domains, around 2% misfolding events were reported in repeated stretch-release cycles in AFM experiments [8]. A slightly larger fraction (∼ 6%) of misfolded species was identified in single-molecule FRET experiments and rationalized in terms of domain swapped intermediates, captured by coarse-grained simulations [7, 11].

In contrast, with the above misfolding-prone systems, there are certain polyprotein chains have been shown be resistant to misfolding, according to pulling experiments. For instance little evidence for misfolding was identified in a polyprotein of GB1 [14] (Fig 2g), with more than 99.8% of the chains (GB1)₈ folding correctly in repetitive stretching–relaxation cycles [14].

Lastly, we consider polyUbiquitin (Fig 2f), for which there is conflicting experimental evidence on misfolding. Initial force microscopy studies showed only the formation of native folds [15], with no misfolding. Later work suggested the formation of collapsed intermediates [22], however the signature change in molecular extension of these was different from that expected for fully domain-swapped misfolds. A separate study using a lock-in AFM [21] found Ubiquitin to conform closely to expectations for a two-state folder, without evidence of misfolding. For this protein, there is a strong imperative to avoid misfolding, since Ubiquitin is initially expressed as a tandem polyUbiquitin chain in which adjacent domains have 100% sequence identity, yet this molecule is critical for maintaining cellular homeostasis [40].

Coarse grained simulation model

A coarse grained structure-based (Go-like) model similar to the earlier work is employed for the study here [7, 41]. Each residue is represented by one bead, native interactions are attractive and the relative contact energies are set according to the Miyazawa–Jernigan matrix. The model is based on that described by Karanicolas and Brooks [41], but with native-like interactions allowed to occur between domains as well as within the same domain, as described below [7]. All the simulations are run under a modified version of GROMACS [42]. For the seven species we studied in this work, the native structures of single domains that were used to construct the models for SH3, SH2, PDZ, TNfn3, Titin I27, GB1 and Ubiquitin correspond to PDB entries 1SHG [43], 1TZE [44], 2VWR, 1TEN [45], 1TIT [46], 1GB1 [47] and 1UBQ [48] respectively. For the single domains of SH3(1SHG), TNfn3(1TEN) and GB1(1GB1), additional linker sequences of Asp-Glu-Thr-Gly, Gly-Leu and Arg-Ser, respectively, are added between the two domains to mimic the constructs used in the corresponding experiments [9, 14, 26]. Construction of the Titin I27 model was described in our previous work [7].

In order to allow for domain-swapped misfolding, the native contact potentials within a single domain are also allowed to occur between corresponding residues in different domains, with equal strength. Specifically, considering each single repeat of the dimeric tandem that has L amino acids, given any pair of residues (with indices i and j) that are the native interactions within a single domain, the interaction energy for the intradomain interaction (E_i,j(r)) is the same as the interdomain interaction between the residue (i or j) and the corresponding residue (j + L or i + L) in the adjacent domain, i.e. E_i,j(r) = E_i+L,j(r) = E_i,j+L(r) = E_i+L,j+L(r).

Kinetic folding simulation of dimeric tandem

To investigate the folding kinetics of the dimeric tandem, a total of 1024 independent simulations are performed on each system for a duration of 12 microseconds each. Different misfolding propensities are observed at the end of the simulations. With the exception of Ubiquitin and GB1, the vast majority of the simulations reached stable native states with separately folded domains. A small fraction of simulations form stable domain-swapped misfolded states. All the simulations are started from a fully extended structure, and run using Langevin dynamics with a friction of 0.1 ps⁻¹ and a time step of 10 fs.

Folding reaction coordinates

We note that all the generated domain-swapped misfolding structures, containing the central and terminal domains, can be monitored by a reaction coordinate based on circular permutated native-like contact sets. Each circularly permuted misfold can be characterized according to the loop position K in sequence where the native domain would be cut to form the circular permutant (K = 0 corresponds to the native fold). If a native contact C_native = (i,j) exists between residues i and j in the native fold, the corresponding native-like contacts for the central (C_in(K)) and terminal domains (C_out(K)) of the domain swapped conformation are generated as where Θ(x) is the Heaviside step function and L is the length of each single domain (plus interdomain linker). S_in,K is the set of native-like contacts C_in of the central domain, and S_out,K is the set of all the native-like contacts C_out of the terminal domain. S_in,K and S_out,K can be used to define a contact-based reaction coordinate to analyze the kinetics of the dimeric tandem misfolding. The corresponding fraction of contacts for the central domain could be calculated by: (1) where N is the total number of domain swapped contacts, S_K = S_in,K ∪ S_out,K (equal to the total number of native contacts), r_ij(χ) is the distance between residue i and j in the protein configuration χ. is the corresponding distance in the native structure for native-like contacts, β = 50 nm⁻¹ and λ = 1.2 is used to account for fluctuations about the native contact distance.

Equilibrium properties and free energy surfaces

The equilibrium properties of a single domain of each system are obtained from umbrella sampling along the native contacts Q as the reaction coordinate. The obtained melting temperature of each system is listed in Table A in S1 Text. A temperature at which the folding barrier ΔG_f of approximately ∼ 2.5 k_BT is chosen for the 2-domain tandem simulations for reasons described below. The stability ΔG_s is calculated as (2) where k_B and T are the Boltzmann constant and temperature respectively. Q_‡ is the position of the barrier top in F(Q), separating the folded and unfolded states and F(Q) represents the free energy profile on Q. Barrier heights ΔG_f were simply defined as ΔG_f = G(Q_‡) − G(Q_u), where Q_u is the position of the unfolded state free energy minimum on Q.

Relative contact order

We calculated the relative contact order [49], RCO_K of different circular permutants K via (3) where L is the length of the single domain, and N is the total number of the native like contacts (the same for different K). S_in,K is the contacts set of the circular permutant corresponding to the “central domain” of the misfolded state. Note that the contact order calculation here is using residue-based native contacts (the same ones defined as attractive in the Gō model), instead of all atom native contacts.

Ising-like theoretical model

An Ising-like model was built based on the native contact map, in which each residue is considered either folded or unfolded and so any individual configuration can be specified as a binary sequence, in a similar spirit to earlier work [50–52]. Interactions between residues separated by more than two residues in the sequence are considered. To simplify the analysis, we also consider that native structure grows only in a single stretch of contiguous native residues (native segment), which means the configurations such as …UFFFUUUUU… or …UUUUUFFFU… are allowed, however, …UFFFUUUFFFU… is not allowed (“single sequence approximation”)[50]. Each residue which becomes native incurs an entropy penalty ΔS, while all possible native contacts involving residues within the native segment are considered to be formed, each with a favourable energy of contact formation ϵ.

The partition function for such a model can be enumerated as: where k_B and T are the Boltzmann constant and temperature. G(χ) is the free energy determined by the number of native contacts n(χ) in the configuration χ, and the number of native residues, N_f(χ). The distribution of the microstates (χ) can be efficiently generated by the Metropolis-Hastings method with Monte Carlo simulation. In each iteration, the state of one randomly chosen residue (among the residues at the two ends of the native fragment and their two neighbouring residues) is perturbed by a flip, from native to unfolded or from unfolded to native, taking the system from a microstate χ₁ with energy E₁ to a microstate χ₂ with energy E₂. The new microstate is subject to an accept/reject step with acceptance probability (4)

To mimic the folding stability difference between native and circular permutant folds, a penalty energy term E_p has been added whenever the native fragment crosses the midpoint of the sequence from either side (the function θ(χ) above is 1 if this is true, otherwise zero). That situation corresponds to formation of a domain-swapped structure, in which there is additional strain energy from linking the termini, represented by E_p. We only use the Ising model here to investigate formation of the first domain (either native or circular permutant), by rejecting any proposed Monte Carlo step that would make the native segment longer than the length of single domain, L.

First passage simulations of misfolding in multidomain proteins

In order to characterize the potential misfolding properties of each type of domain, we have used a Gō-type energy function based on the native structure. Such models have successfully captured many aspects of protein folding, including ϕ-values [53, 54], dimerization mechanism [55, 56], domain-swapping [57–60], and the response of proteins to a pulling force [61, 62]. More specifically, a Gō type model was used in conjunction with single-molecule and ensemble FRET data to characterize the misfolded states and misfolding mechanism of engineered tandem repeats of Titin I27 [7, 12]. We have therefore adopted the same model. Although it is based on native-contacts, it can describe the type of misfolding we consider here, which is also based on native-like structure. Note that this model effectively assumes 100% sequence identity between adjacent domains, the scenario that would most likely lead to domain-swap formation. It is nonetheless a relevant limit for this study, as there are examples in our data set of adjacent domains having identical sequences which do misfold (e.g. titin I27) and those which do not (e.g. protein G).

For each of the folds shown in Fig 2, we ran a large number of simulations, starting from a fully extended, unfolded chain, for sufficiently long (12 μs each) such that the vast majority of them reached either the correctly folded tandem dimer, or a domain-swapped misfolded state similar to that shown in Fig 1e for titin. In fact, for each protein, a number of different misfolded topologies are possible, illustrated for the Src SH3 domain in Fig 3. Each of these domains, shown in conventional three-dimensional cartoon representation in the right column of Fig 3 and in a simplified two-dimensional topology map in the left column, consists of two native-like folded (or misfolded) domains. For convenience, we call the domain formed from the central portion of the sequence the “central domain” and that from the terminal portions the “terminal domain”. We have chosen to characterize each topology in terms of the position, K, in sequence after which the central domain begins. Thus, the native fold has K = 0, and all the misfolded states have K > 0. Typically, because of the nature of domain swapping, K must fall within a loop. Of course, there is a range of residues within the loop in question that could be identified as K and we have merely chosen a single K close to the centre of the loop. This position, and the central domain, are indicated for the Src SH3 misfolded structures in Fig 3. We note that each of these central domains can also be considered as a circular permutant of the native fold, in which the ends of the protein have been joined and the chain has been cut at position K.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Folded and misfolded topologies of Src SH3.

(a) Schematic of Src SH3 fold, in which the three-dimensional β sheet structure (shown in (b)) is unrolled into two dimensions, for each domain (N-terminal and C-terminal in blue and red respectively). On the N-terminal domain are indicated the sequence positions K ∈ {0, 18, 37, 46} characterizing the possible circularly permuted “central domains”, with K = 0 corresponding to the native fold and K > 0 indicating the approximate starting residue for the “central domain” misfold. (c), (e), (g): two dimensional representations of the observed misfolded topologies of Src SH3. In each case, the residue K characterizing the misfold is indicated by the bullet point and the central domain is enclosed by a broken rectangle. (d), (f), (h): three-dimensional representations of the misfolds shown in (c),(e),(g) respectively.

https://doi.org/10.1371/journal.pcbi.1004933.g003

With this nomenclature in hand, we can more easily describe the outcome of the folding simulations for the seven domain types considered in terms of the fraction of the final frames that belonged to the native fold, versus each of the possible misfolded states. These final populations are shown in Table 1. We see that for five of the domains (SH3, SH2, PDZ, TNfn3, Titin I27), misfolded structures are observed, with total populations ranging from 5–10%. For the remaining two domains, Ubiquitin (UBQ) and protein G (GB1), no misfolded population is observed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Summary of misfolding statistics and central domain properties.

K labels the type of fold/misfold (see text; K = 0 is native); RCO is relative contact order [49]. ΔG_f and ΔG_s are the folding barrier and stabilities of a single folded/misfolded domain. Population is frequency of each state at the end of the 1024 trajectories. Maximum standard error on populations is 1.6% for a sample size of 1024. Numbers in brackets are rank correlations with folded/misfolded populations.

https://doi.org/10.1371/journal.pcbi.1004933.t001

Consistency with existing experimental data

The ability to capture domain-swapped misfolds with simple coarse-grained simulations potentially allows us to investigate the origin of the misfolding, and its relation, if any, to the topology of the domain in question. However, we also need to benchmark the accuracy of the results against experiment as far as possible, in order to show that they are relevant. There are two main sources of information to validate our results. The first is the overall degree of domain-swapped misfolding for those proteins where it has been characterized, for example by single molecule AFM or FRET experiments. Qualitatively we do observe good agreement, where data is available: in experiment, domains which have been shown to misfold are TNfn3 (AFM) and Titin I27 (AFM, FRET), which are both found to misfold here, while there is no detectable misfolded population for protein G (AFM), again consistent with our results. We also do not observe any misfolding for Ubiquitin, consistent with the lack of experimental evidence for fully domain-swapped species for this protein [15–23].

Quantitatively, the fractional misfolded population is also consistent with the available experimental data. For instance, the frequency of misfolded domains in native tenascin is ∼ 4% as shown by previous AFM experiments [8], the misfolded population of I27 dimers is ∼5% in single-molecule FRET experiments [7] while the misfolded population of GB1 domains in polyproteins (GB1₈) is extrememly low (< 0.2%) [14]. Even though the observed misfolding population of the misfolded tandem dimer is low, it is potentially a problem considering that many of the multidomain proteins in nature have large number of tandem repeats, such as Titin which contains twenty-two I27 repeats [63]. Recent FRET experiments on I27 tandem repeats have shown that the fraction of misfolded proteins increases with the number of repeats. For the 3- and 8-domain polyproteins, the fraction of misfolded domains increases by a factor of 1.3 and 1.8, respectively, relative to a tandem dimer [12].

The second type of evidence comes from experimental structures of domain-swapped dimers. For several of the proteins, bimolecular domain-swapped structures have been determined experimentally. While no such structures have yet been determined for single-chain tandem dimers, we can compare the misfolded states with the available experimental data. For each experimental example, we are able to find a corresponding misfolded species in our simulation with very similar structure (related by joining the terminis of the two chains in the experimental structures). The domain swapped dimers solved obtained from experiments (Fig 4a, 4c, 4e and 4g) are strikingly similar to the domain swapping dimeric tandem from simulations, which are the domain swapped SH3 domains when K(sequence position after which the central domain begins) = 37 (Fig 4b), SH2 with K = 72 (Fig 4d), TNfn3 with K = 28 (Fig 4f) and PDZ with K = 23 (Fig 4h). Most of these states have relatively high population among all the possible misfolds as observed from the simulations (“Population” in Table 1). While the coverage of possible domain swaps is by no means exhaustive, the observed correspondence gives us confidence that the misfolded states in the simulations are physically plausible.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Comparison of domain-swapped misfolds with experimental structures.

Selected misfolded dimeric tandems obtained from the simulations (right column) are compared with corresponding experimental structures (solved by crystallography or NMR) of domain-swapped dimers involving two separate protein chains (left column). The proteins are, from top to bottom (a),(b): SH3, (c),(d): SH2, (e),(f): TNfn3 and (g),(h): PDZ domains The PDB accession codes are 1I07, 1FYR, 2RBL and 2OSG respectively.

https://doi.org/10.1371/journal.pcbi.1004933.g004

Circular permutants as models of misfolding intermediates

Having shown that the misfolding propensities we obtain are qualitatively consistent with experimental evidence (and in the case of Titin I27, in semi-quantitative agreement with single-molecule FRET), we set out to establish some general principles relating the properties of each domain to its propensity to misfold in this way. We can start to formulate a hypothesis based on the alternative folding and misfolding pathways illustrated in Fig 1. Native folding has as an intermediate a state in which either the N- or the C-terminal domain is folded. In contrast, on the misfolding pathway, the first step is formation of the central domain, followed by that of the terminal domain. This parallel pathway scheme suggests that a descriptor of the overall misfolding propensity may be obtained from the rate of formation of a single correctly folded domain, relative to that of the central domain (neglecting back reactions, because this are rarely seen in our simulations). We can study the central domain formation in isolation, since these structures are just circular permutants of the native fold, i.e. the two proteins have the same sequence as the native, but with the position of the protein termini moved to a different point in the sequence, as is also found in nature [35]. These structures can be thought of as originating from the native by cutting a specific loop connecting secondary structure elements (the free energy cost of splitting such an element being too high), and splicing together the N- and C- termini. In the context of the tandem dimers, the position at which the loop is cut is the same K that defines the start of the central domain in sequence.

We investigate the role of the central domain by characterizing the free energy landscape of the single domain of each system, as well as all of its possible circular permutants, using umbrella sampling along the reaction coordinate Q_K. Q_K is exactly analogous to the conventional fraction of native contacts coordinate Q [64], but defined using the corresponding (frame-shifted) contacts in the circular permutant pseudo-native structure. The index K indicates the position along the sequence of the WT where the cut is made in order to convert to the circular permutant.

The free energy surfaces F(Q_K) of two representative systems, SH3 and Ubiquitin, are shown in Fig 5, with the data for the remaining proteins given in the Fig A in S1 Text. The free energy barrier height for folding ΔG_f and the stability ΔG_s are listed in the Table 1. The free energy plots indicate that the single domains of Ubiquitin and GB1 are stable only for the native sequence order, and not for any of the circular permutants. Based on the type of misfolding mechanism sketched in Fig 1, one would expect that unstable circular permutants would result in an unstable central domain, and consequently no stable domain-swappping misfolding would occur in the dimer folding simulations, as we indeed observe. This is also consistent with previous studies of polyproteins of GB1 and Ubiquitin using using AFM experiments, which reveal high-fidelity folding and refolding [14, 65, 66]. We note that only under very strongly stabilizing conditions is any misfolding observed for ubiquitin dimers: running simulations at a lower temperature (260 K), we observe a very small (1.3%) population of misfolded states from 1024 trial folding simulations. At a higher temperature of 295 K, once again no misfolding is observed.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Free energy profile of WT and its circular permutant domains.

The structures of SH3 and Ubiquitin are shown in (a) and (c), with the “cut” positions K in the WT to form circular permutant labeled with crosses. (b) is the free energy surfaces F(Q) of WT SH3 as well as its circular permutants at 300K. (d) is the F(Q) of WT Ubiquitin and its circular permutants at 290K. The labels K indicate the residue index of the cut position. The free energy curves of the circular permutant cases are shifted vertically for visual clarity, and coloured using the colours corresponding to the crosses in (a) and (c). The free energy plots of the other systems: GB1, SH2, TNfn3 and PDZ are shown in Fig A in S1 Text.

https://doi.org/10.1371/journal.pcbi.1004933.g005

In contrast to the situation for GB1 and Ubiquitin, all of the circular permutants of the SH3 domain in Fig 5 are in fact stable, although less so than the native fold. The destabilization of circular permutants relative to native is in accord with the experimental results for the Src SH3 domain [26] (rank correlation coefficient stabilities is 0.80). The other domains considered also have stable circular permutant structures. This is consistent with the fact that all of these domains do in fact form some fraction of domain-swapped misfolded states. The simplest view of the misfolding mechanism would be as a kinetic competition between the correctly folded intermediates versus the domain-swapped intermediates with a central domain folded (i.e. a “kinetic partitioning” mechanism [67]). In this case one might naively expect that the propensity to misfold would be correlated with the relative folding rates of an isolated native domain and an isolated circular permutant structure. However, the folding barriers ΔG_f projected onto Q (for native) or Q_K (for circular permutants) show little correlation to the relative frequency of the corresponding folded or misfolded state, when considering all proteins (Table 1). Since this barrier height may not reflect variations in the folding rate if some of the coordinates are poor (yielding a low barrier) or if there are large differences in kinetic prefactors, we have also directly computed the folding rate for the circular permutants of those proteins which misfold, and confirm that the rates of formation of the native fold and circular permutants are similar. We indeed obtain a strong correlation between the folding rate of the isolated circular permutant and the folding barrier ΔG_f (Table B in S1 Text), which implies Q is a sufficiently good reaction coordinate here. We have also considered the relative contact order (RCO) as a proxy for the folding rate, since it has been found to correlate with folding rates for two-state folding proteins [49, 68]. However, the RCO calculated based on the native or circularly permuted folds did not correlate with either the barrier height for single domain or circular permutant folding, or with the extent of misfolding in dimeric tandem proteins (Table 1). Since the folding rates do not explain misfolding propensities by themselves, another possibility is that the reverse reactions have to be considered. However, once they had formed, in most cases we did not observe unfolding of the first native domain, or of the intermediate with central domain folded, indicating that back reactions should not be needed, at least to explain the simulation data. This lack of refolding is a consequence of the significant stability of the native folds, which controls the relative folding and unfolding rates (and indirectly, those of the circular permutants). Under these conditions, given that folding rates are much higher, once a native fold (or circular permutant misfold) has formed, it is much more likely that a second domain will fold, rather than the first domain unfolding. Our choice of stabilizing conditions was motivated by the fact that misfolding is observed in experiment under conditions where the folded state is much more stable, and the stabilities (ΔG_s) of the folded single domains in our simulations are generally comparable to those in experiment (experiment vs simulation, UBQ: 6.1 vs 4.2, GB1: 5.3 vs 3.1, PDZ: 7.5 vs 4.5, SH3: 4.1 vs 9.2, Titin: 7.5 vs 8.1, Tnfn3: 5.3 vs 8.1 kcal/mol) [26, 69–72].

On the other hand, we did note that there was a significant, and unexpected, correlation between the population of the final folded or misfolded states and the stability ΔG_s of the corresponding intermediate. Spearman rank correlation coefficients between the folded stability ΔG_s of the intermediate structure and the frequency of folded/misfolded states were 0.63, 0.94, 0.74, 0.81, 0.86 for the SH3, PDZ, TNfn3, SH2 and Titin I27 domains respectively. We note that there is also a reasonable correspondence between the relative stabilities of circular permutants in simulation and experiment, where data are available [12, 26]. How can the correlation with stabilities rather than folding rates of the isolated domains be understood? The resolution lies in the difference between the folding to either type of intermediate represented in Fig 1, and folding of the single domain “models” for these species, namely that the intermediates fold in the context of the full sequence. This is important because a large fraction of native (or native-like) contacts are shared between the native fold and the various misfolded domains. As such, the native and misfolded states can be considered as belonging to the same folding funnel, with differentiation between the two occurring at a late stage of folding. This scenario is illustrated schematically in Fig 6, in which folding to either a state with one native domain folded (on left), or one possible domain-swapped misfolded intermediate (on right) are considered. The states of the proteins are represented by very coarse-grained contact maps (e.g. representing contacts between pairs of β-strands [73], rather than between residues). As can be seen, dividing this funnel into the separate funnels by considering only native contacts for the native or circularly permuted fold would be misleading (green and red funnels respectively), since the two funnels share several configurations, and many of their states can be converted to one in the other funnel by flipping a single coarse-grained “residue” between folded and unfolded states.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Folding/misfolding funnel.

Illustration of relation between folding funnels for native and domain-swapped domains. (A) Example native contact map, highly coarse-grained for simplicity. (B) Map of all possible native-like contacts for a two-domain protein, showing native contacts in black and domain-swapped contacts in blue. (C) In the context of the two-domain sequence, the folding funnels for a single native domain (green broken line) and domain-swapped domain (red broken line) are interconnected, forming part of a single global funnel (black line). States are considered part of the native funnel if all contacts formed belong to the native state, and to the domain-swap funnel if all contacts formed belong to the domain-swapped structure. Note that only a subset of possible states are shown, for clarity (e.g. other domain-swapped species are possible). Only states with a single native-like stretch of residues are considered, whose length does not exceed that of a single folded domain. Arrows connect states differing by a single coarse grained residue flipping between native and non-native.

https://doi.org/10.1371/journal.pcbi.1004933.g006

We can see this explicitly by plotting some representative folding transition paths from the Src SH3 dimer simulations. In Fig 7 top row, we show a folding event for a simulation which forms a native fold (at the N-terminus), and in the bottom row, for a simulation which forms a circularly permuted central domain with K = 18. Each event is projected onto two different reaction coordinates,Q_K, for K = 0 (standard native Q) and K = 18 (the Q when the circular permutant for K = 18 is considered as “native”). As is evident, a large fraction of the transition path looks very similar in Fig 7b and 7d, with contacts that could be considered equally as native-like or central domain-like being formed initially in the lower left part of each plot. Around Q₀ ≈ Q₁₈ ≈ 0.5, the first trajectory moves toward the native structure, where it terminates (Fig 7b). The second trajectory also deviates initially more toward the folded structure, but then switches back near the end to form the central domain structure instead (Fig 7d). A similar branching of folding pathways has also been proposed in a recent computational study of domain swapped dimer formation [74].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Transition paths for the formation of the first (folded or misfolded) domain in tandem SH3 dimers.

Folding (a) and misfolding (c) kinetics are projected along the reaction coordinate Q_in, where two different kinds of Q_in are chosen depends on K, which is for the native fold and the circular permutated misfold when K = 0 and 18, respectively. Transition-path segments are defined as being between Q_K = 0.1 and Q_K = 0.9. In the right panels, the same trajectories are projected onto the Q_{K = 0} and Q_{K = 18} (panel (b) for trajectories in (a) and (d) for those in (c).

https://doi.org/10.1371/journal.pcbi.1004933.g007

The common funnel picture helps to explain why the stability of the isolated native or circular permutant domains may be correlated with their frequency of formation in the context of the full length sequence in which either could potentially be formed. Initially, nucleation of folding could occur by formation of native contacts anywhere in the sequence. Indeed, they are most likely to form near the centre of the chain. However, as more native/native-like structure is accumulated, the nascent, partially folded protein will be biased to form the contacts leading to the lower free energy structure, and so the folding nucleus is likely to move towards one of the termini of the protein. We note that while a previous study suggested that the stability of the individual domains might be affected by conjugation to another folded domain [75], this is unlikely to be relevant because in our case the misfolding is controlled by formation of the first domain, while the second domain is still unfolded.

Further insight into how the above free energy bias influences the outcome of the folding kinetics can be obtained by considering the progressive formation of folded structure. In order to characterize the location of nascent folded structures, we define a new order parameter representing the average position of native contacts along the sequence, (5) where (i, j) is the native or native-like contact formed by the residues i and j in the configuration χ, and S(χ) is the set of all such contacts which are formed in χ. We can locate the position of nascent structure in the sequence by plotting the distributions of for χ drawn from the equilibrium distribution at selected values of the global coordinate Q, defined as the fraction of native contacts in the native dimer structure (i.e. Q = 0.5 corresponds to a single folded or misfolded domain; both native and native-like contacts are counted, and divided by the total number of contacts in the native state). Fig 8 shows that early in folding, at low Q values (shaded histograms in Fig 8), the distribution of is broad, and centered in the middle of the sequence. This implies that folding could potentially begin at many positions along the sequence, with no initial preference for folded or circularly permuted structure. However, as folding proceeds closer to formation of a complete domain, develops two maxima, one in the N-terminal and one in the C-terminal part of the chain, corresponding to native domain formation. The nascent native-like structure thus naturally migrates towards the termini to avoid the free energy penalty of forming a circularly permuted misfolded intermediate.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Distribution of the “folding nucleus” location from the tandem dimer simulations (Table 1).

The of the (a) (SH3)₂ (b) (SH2)₂ (c) (TNfn3)₂ and (d) (PDZ)₂ are extracted at two different Q on the folding pathway (see individual figure legends for Q values). Note that the Q ∼ 0.5 corresponds to the structure with the first domain fully formed. The spread of contacts in sequence, within a given conformation, also becomes narrower with increasing Q (Fig B of S1 Text).

https://doi.org/10.1371/journal.pcbi.1004933.g008

Ising-like theoretical model

The results from the previous sections show that the misfolding propensity is highly correlated with the the stability of the isolated native domain and its circular permutants. To further explain how this might occur, we investigate the dependence of the misfolding propensity on the stability of the central domain in the context of full sequence (dimeric tandem). We have constructed an even simpler simulation model for formation of the first intermediate (native, or circularly permuted), by using a simulation of a Wako-Saito-Muñoz-Eaton Ising-like model [51, 76]. In the version we consider here, each residue is considered either to be folded or unfolded, so that each configuration can be described as a binary string. Furthermore, we impose the single-sequence approximation, namely that all native-like structure forms in a single segment of contiguous residues. We also restrict the number of folded residues to be at most one half of the dimer sequence length, so that only a single folded or misfolded domain can form, the aspect we are most interested in. To model the stability difference between native and circularly permuted domains, we introduce an additional energy penalty E_p for any folded segment which crosses the midpoint of the dimer sequence. Such a folded segment must be forming a circular permutant misfold and as such will incur some additional “strain” energy from joining the termini of the original fold.

We show results from a typical Monte Carlo trajectory for this model in Fig 9. We have used two parameters to characterize the results, the fraction of native or native-like contacts, Q_res, and (Eq 5). Q_res equals to the number of residues which are in the native-like state divided by the total number of residues of one domain (L). The projection of a trajectory for the model onto in Fig 9a shows that the most stable states occur for in the center of either the first of the second natively folded domain. Nonetheless, there are other stable states at intermediate , which correspond to the circular permutant intermediates. These have a lower stability, because a value E_p > 0 was used in this instance. The effect of the stability penalty for the circular permutants is illustrated by the two-dimensional free energy surfaces in Fig 9b–9d. In all cases, there are minima at low Q_res corresponding to unfolded structures and at high Q_res for folded (native or circular permutants). If the penalty E_p = 0 (Fig 9b), in addition to the stable native folds at and , there are a variety of other free energy minima at high Q_res corresponding to circular permutants, which have essentially identical free energy to the native fold. However, as E_p is increased, the relative population of these misfolded states decreases (Fig 9c and 9d), as expected.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Ising-like model.

a) Monte Carlo simulation trajectory segment when E_p = 4.0 kcal/mol. The free energy profile of vs Q_res changes when the difference of the stability between the WT and the circular permutant become larger and larger, in which cases the E_p are b) 0.0, c) 4.0 and d) 6.0 respectively. All the free energy plots are at the temperature T = 525 K.

https://doi.org/10.1371/journal.pcbi.1004933.g009

Predicting circular permutant stability using alchemical free energy method

Knowing that misfolding correlates with the relative stability of the native single domain and its circular permutants is useful because it suggests a means to predict the likelihood of misfolding, provided one has an estimate of the circular permutant stability. While one could attempt to determine this experimentally, by synthesizing the circular permutants, or computationally, as we have done, it would be very helpful to have a quick method to estimate this stability a priori [77]. Here, we have developed such a method, based on an alchemical transformation from the native to the domain-swapped misfolded state: the overall conversion between the native and circular permutant can be expressed as the sum of the free energy changes of two steps as shown in Fig 10: firstly, joining the N- and C- termini of WT (Fig 10a) to form a cyclic intermediate state (I) (Fig 10b), in which procedure the free energy change is (ΔG_J). Note that this step is the same for all circular permutant folds. The second step is to cut different loops on the cyclic configuration I to form a circular permutant CP (Fig 10c) with the change of free energy (ΔG_C):

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Alchemical transformation from native to circular permutant.

(a) Native structure; (b) cyclized structure; (c) circular permutant after cutting another loop.

https://doi.org/10.1371/journal.pcbi.1004933.g010

Assuming the native and circular permutant unfolded states to have the same free energy, the overall change in stability between the circular permutant (CP) and the wild-type (WT) is:

For the first step (Fig 10a and 10b), in order to join the two termini, they must be sufficiently close. In general, bringing them closer together will require peeling off a small part of the native structure, starting from the termini. If we imagine that all of the structure between residues i at the N-terminus and j at the C-terminus remains native, then the change in free energy for linking the termini for cyclization, ΔG_J, can be split into energetic and entropic components:

Assuming the states of the residues p ∈ {0 .. i − 1} and q ∈ {j + 1 .. L}, which are on the N- and C- termini respectively, change from the native state to non-native state (joint loop), the total energy increase will be ∑_{x ∈ σ(p)} ϵ_px + ∑_{x ∈ σ(p)} ϵ_qx, which is the summation of all the native contact energy, in the Go model over the sets of residues σ(p), σ(q) involving residues p and q respectively. x represents the residues that form the native contacts with either p or q. We approximate the entropy gained per residue as δs, where δs = ∑_native(i,j) ϵ_ij/(TN), where N is the number of residues and T is the folding temperature. The gained entropy is set to 0 if residue p or q does not have any contact with other part of the protein except for the neighboring residues, and the number of such residues is denoted by κ. The average length contribution (r₀) of peeling off each residue from the native structure is set to 3.5 Å here. The topological requirement of joining the two termini by peeling off residues 1 to i − 1, and j + 1 to L from the native state is that the linear distance between the residues i and j (d(R_i,R_j)) on the native structure is shorter than the effective length contributed by the joint parts: (6) Note that if N- and C- termini point in opposite directions, such as the TNfn3, Titin I27, UBQ and GB1 domains (Fig 2), around six residues (three on each side) of the two termini will form the turn of the joint loop which does not contribute to the the effective length. Therefore, an offset number M = 6 is used in this case. This is justified because turns in proteins are usually defined by four residues (or 3 residue-residue bonds) [78]. For SH3, SH2 and PDZ domains (Fig 2), whose N- and C- termini align to the same direction, M is set to 0. With the above condition (Eq 6), the minimum overall change of ΔG_J by adjusting i and j could be given by: where

Analogously, for the second step (Fig 10b and 10c), assume the loop is cut at the position between residue position K and K + 1, the states of the residues on each side of the cutting point p′ ∈ {i′,.., K} and q′ ∈ {K + 1,.. j′}, will change from the native state to the non-native state. The gained entropy per residue is δs. ∑_{x ∈ σ(p′)} ϵ_p′x and ∑_{x ∈ σ(q′)} ϵ_q′x are the summation of all the energy of the native contacts which are broken due to the cutting. Therefore, by comparing different combinations of i′ and j′, the minimum change of stabilities ΔG_C in this step is: where

The ΔΔG_tot calculated using this alchemical free energy method is very well correlated with the stability of the circular permutant ΔG_s obtained by umbrella sampling simulations (Table 1) as shown in Fig 11. It is consistent with the experimental results that GB1 and Ubiquitin have the most unstable circular permutant folds in general. The main contribution of ΔΔG_tot is from ΔG_J, since the enthalpy penalty is large when many native contacts are broken by joining the terminis such as in the case of UBQ and GB1. The free energy cost of cutting the loop (ΔG_c) is relatively small and is similar for all the circular permutants during the transformation (Fig 10).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. ΔΔG_tot from alchemical model vs ΔG_s (Table 1).

https://doi.org/10.1371/journal.pcbi.1004933.g011

Effect of linker length

From the alchemical method we can see that the difference of stability ΔΔG largely depends on the native contacts that are broken in the procedure when joining the N- and C- termini. However, ΔΔG could also be lowered by extending the linkers between domains, for instance, by adding extra residues at the two termini. If the loop formed by the linker is long enough, fewer native contacts will need to be disrupted, so that the circular permutant folds would be more stable. Therefore we have investigated the stability of circular permutant folds as a function of the length of C-terminal extension, by adding Gly-Ser repeats (forming no native contacts). This extra peptide corresponds to the linker between the tandem domains. The stability of circular permutants, obtained from simulations using umbrella sampling, as a function of linker length (ll) is shown in Fig 12 (raw potentials of mean force on Q_K in Figs C and D in S1 Text). As expected, longer linkers between the tandem repeats give more stable circular permutants. The relative change from ll = 0 to ll = 20 is roughly the same for all circular permutants of a given protein, as expected since the change in all cases is the same loop extension. Note, however, that the effect is much larger for ubiquitin than for TNfn3. To investigate the consequences of the change of central domain stability for the misfolding propensity of tandem repeats, we carried out first passage simulations of a tandem dimer of Ubiquitin with linker lengths of 5,8 and 10 residues respectively. The setup of the dimer simulations was the same previously. For each linker length, 1024 independent simulations were run from fully extended structures. No domain-swapped misfolding was found for ll = 5 and ll = 10, however, we indeed obtained three domain-swapped misfolding events for ll = 8. Two of the misfolds belong to the K = 61 (Fig 12a) type and the other one is K = 36 type domain-swapping. As one can see from Fig 12a, the circular permutants K = 36 and K = 61 are the ones which are most stable with ll = 8. However, they are still somewhat unstable, explaining the small fraction of misfolded states obtained. In this case it is clear that the length of the linker between the termini of ubiquitin is one way in which domain-swapped misfolding is avoided in this protein: since ubiquitin is synthesized initially as an N-C linked polyubiquitin chain [40], it is essential to avoid such misfolding, given the importance of this protein to cellular homeostasis. It should be noted though that the influence of the linker depends very much on the protein, as might be expected, from the much smaller effect on the stability of circular permutants of TNfn3 than those of ubiquitin in Fig 12. In experiments on titin I27, the misfolded population was, within error, the same with and without the addition of a four residue RSEL linker [7].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Stability of different CP of Ubiquitin (a) and TNfn3 (b) with different linker lengths.

The linker sequence composition is (GS)₂-S, (GS)₅, (GS)₇-S and (GS)₁₀ giving ll of 5, 10, 15 and 20 respectively.

https://doi.org/10.1371/journal.pcbi.1004933.g012

Lastly, we comment briefly on the effect of linker composition. Although we treat the linkers as structureless chains, not forming native contacts, there may be some effect of linker flexibility, arising from the backbone dihedral potential in our model. To test for this effect, we have carried out an additional 1024 independent simulations with the dimeric tandem repeat of the SH3 domain using a different four residue linker composition, GGGG, rather than the original DETG, as used in the original circular permutant studies by Serrano et al [26]. With the new linker GGGG, the observed misfolded populations are 94.3%, 1.3%, 2.3% and 2.1% for K = 0, 18, 37 and 46 respectively. The differences are not statistically significant compared to the results with the original linker.

Conclusions

We have investigated the factors which favour formation of domain-swapped misfolded states in multidomain proteins, by building on knowledge of the folding/misfolding mechanism. Counter to our original expectations, the misfolding yield does not depend primarily on the relative folding rates of the native single-domain protein and its circular permutants, representing intermediates for correct folding and misfolding respectively. Although the folding rates of wild-type and circular permutants may often be quite similar, the fraction of misfolded protein is much smaller than this comparison would suggest. Instead, it appears that misfolding is correlated with the stability of the native single-domain protein relative to its circular permutants. This can be understood because the rate of formation of the first intermediate (native-like or misfolded) occurs in the background of the full-length sequence. In this context, while folding may be initiated at any point in the chain, the nascent structure will tend to migrate towards the N- or C-terminus because of the free energy bias towards the native fold; circular permutants invariably pay a cost in stability for joining the protein termini. Thus the folding rate of isolated circular permutants relative to wild-type protein may not be a good proxy for these rates in the context of the full length sequence, whilst the domain stability is a better guide as to the free energy bias towards a particular structure. This suggests that the rates of formation of these domains inferred from single-molecule experiments [12] should be interpreted as the rates in the context of the full length sequence. In our analysis, we have neglected the effect of back-reactions. Since these occurred rarely in the simulations, they were not needed to explain the results. We have also quantified the effect of linker length on domain swapping, finding that sufficiently long linkers can permit misfolded species to form in cases where they did not for the native spacing. Finally, we have developed a simple model for predicting the stability of misfolded intermediates (circular permutants of native), which should prove useful for determining whether a given protein may be susceptible to this type of misfolding.