Insights into the Folding of Disulfide-Rich μ-Conotoxins.

The study of protein conformations using molecular dynamics (MD) simulations has been in place for decades. A major contribution to the structural stability and native conformation of a protein is made by the primary sequence and disulfide bonds formed during the folding process. Here, we investigated μ-conotoxins GIIIA, KIIIA, PIIIA, SIIIA, and SmIIIA as model peptides possessing three disulfide bonds. Their NMR structures were used for MD simulations in a novel approach studying the conformations between the folded and the unfolded states by systematically breaking the distinct disulfide bonds and monitoring the conformational stability of the peptides. As an outcome, the use of a combination of the existing knowledge and results from the simulations to classify the studied peptides within the extreme models of disulfide folding pathways, namely the bovine pancreatic trypsin inhibitor pathway and the hirudin pathway, is demonstrated. Recommendations for the design and synthesis of cysteine-rich peptides with a reduced number of disulfide bonds conclude the study.

Introduction

Conotoxins are neuropeptides from the venom of marine cone snails, which interact with a wide range of biological targets (e.g., ion channels, transmembrane receptors, and transporters) and hence are of pharmaceutical interest and of great potential as molecular probes to study the specific subtypes of ion channels and receptors.[1,2] Conotoxins consist of approximately 10–50 amino acid residues and are classified according to their cysteine patterns.[3,4] The typical CC–C–C–CC pattern defines the framework for conotoxins of the M-superfamily comprised by ψ-, μ-, and κM-conotoxins.[3,5,6] The family of the 27 currently known μ-conotoxins[4,7] selectively binds to the ion channel pore of the voltage-gated sodium channels, thus blocking the influx of sodium ions into the cell. Therefore, these peptides also gained interest as useful tools for research studies in electrophysiology.[100−9] μ-Conotoxins are cysteine-rich peptides consisting of 6 cysteines which can give rise to 15 conformational isomers of different disulfide connectivities in various combinations of disulfide bonds.[9] However, the dominant isomer among them bears the disulfide linkage of C1–C4, C2–C5, and C3–C8, often referred to as the “native fold”.[9] The molecular principles underlying the folding bias contributing to one particular dominant isomer is as yet unclear. This subject is of great importance because an accumulation of cysteines may also occur in distinct regions of larger peptides and proteins such as that observed in defensins,[10,11] resistins,[12] Kunitz serine protease inhibitors,[13] and various growth factors.[12,14,15] In such cases, any information on their preferred disulfide connectivity would be relevant especially in the absence of actual chemical or structural data.

Conotoxins represent a promising tool for studying the impact of disulfide bonds on the folding process owing to their small to medium size which makes them intermediates between peptides and proteins and also their high disulfide bond content.[3,5,16] We studied the μ-conotoxins GIIIA, KIIIA, PIIIA, SIIIA, and SmIIIA (Figure ) using unbiased all-atom molecular dynamics (MD) simulations performed on their NMR structures. The simulation output was correlated with the in vitro data from the oxidation reaction for all the μ-conotoxins that had been described earlier.[100−9] The individual μ-conotoxins could be grouped based on the reaction product yield and the side product formation,[17] which allowed for a functional comparison with the computational study. For the MD studies, a rather unconventional approach was pursued by analyzing the process of refolding; that is, the behavior and stability of an individual peptide was observed upon successive opening of disulfide bonds in the folded peptide NMR structures. The analyses conducted on the resulting MD trajectories gave rise to inferences on characteristic factors contributing to the conformational stability of the folded conopeptides. On the other hand, we used the distances sampled by the sulfur atoms between the reduced cysteine residues to observe and determine if the peptide with a disulfide bond removed in silico tries to refold or not. This enabled drawing a distinction between the influence of the disulfide bond and the properties exerted by the rest of the sequence to the maintenance of close-to-native backbone conformations in peptides with computationally reduced disulfide bonds. This consequently simulates the state at which the peptide is at its final stage of folding and in its near-native conformation.

Figure 1

Structure and surface representations of the investigated μ-conotoxins. (A) NMR structures of the five μ-conotoxins (GIIIA, KIIIA, PIIIA, SIIIA, and SmIIIA) used in this study represented as a cartoon. The secondary structure elements, α-helix (purple), 3–10 helix (blue), turn (cyan), and coil (white) were generated by STRIDE[19] in visual molecular dynamics (VMD). The cysteine residues forming the disulfide bonds (yellow) were labeled. (B) Molecular surface was generated by SURF[20] in VMD, indicating the hydrophobic (white), basic (blue), acidic (red), and hydrophilic (green) regions. All structures were taken from the ConoServer database.[7,21]

The pathways for disulfide folding have been classified into two extreme models despite exhibiting a high degree of diversity.[18] One model is represented by the bovine pancreatic trypsin inhibitor (BPTI)-like folding, where there is a predominance of native intermediates at various steps down the folding funnel, and the other extreme, the hirudin pathway, is defined by highly heterogeneous non-native intermediates. Interestingly, conotoxins are so far placed in between these two extreme models in a hybrid BPTI–hirudin model.[18] In the present study, we thus try to improve the clarity regarding the classification of the peptides used herein considering the existing models. Finally, an attempt to relate the simulation inferences from these peptides with open disulfide bonds to their propensity (or) favorability to retain close-to-native conformations is presented via a qualitative grouping.

Results and Discussion

Studies of conotoxin folding include, in general, experimental approaches such as regioselective oxidation strategies, spontaneous oxidative self-folding, and optimization of the folding methods tested, recombinant expression, and so far available information from the biosynthesis of conotoxins in the venom duct of cone snails. However, it is still unresolved how exactly the cone snails produce properly folded peptide and protein toxins. In order to approach the mechanisms behind conotoxin folding, several studies focused on using computational strategies to provide valuable insights into the structural features important for the folding process, but often a combination of both experimental and theoretical work is missing or only one conotoxin was in the focus of most of these studies.

Oxidative Self-Folding of μ-Conotoxins

In order to provide the in-house experimental data for comparison with the computational results, equal amounts of each of the selected five μ-conotoxins (Figure , Table ) were used for the oxidation reaction according to a protocol earlier described yielding an undirected and spontaneous formation of disulfide bonds.[100−9] Under the conditions applied, oxidation of μ-SmIIIA and μ-PIIIA resulted in several peaks of fully oxidized product(s) as confirmed by mass spectrometry; that is, the disulfide connectivity of the individual fractions is different from the native (major) fold.[9] In contrast, oxidation of all other peptides, that is, μ-GIIIA, μ-KIIIA, and μ-SIIIA, resulted in one major product confirming earlier reports (Figure ).[17,22,23]

Table 1

Comparison of Sequence Characteristics of the μ-Conotoxins Investigated in This Studya

Residues are highlighted according to their character: basic (blue), acidic (red), polar uncharged (green), and cysteine (yellow) (Z: pyroglutamate, O: 4-hydroxyproline). *All peptides were used as amides. Native μ-SmIIIA occurs as the C-terminal acid, however, is usually used as an amide.[100,8] We used μ-SmIIIA as an amide for reasons of comparison because both structures were found to be identical.[30]#In general, the IC50 values determined for the toxins ion channel blocking activity at the skeletal muscle NaV1.4 expressed in Xenopus oocytes are given, with the exception of SmIIIA where only KD is available.

Figure 2

High-performance liquid chromatography (HPLC) elution profiles of μ-conotoxins studied. Linear, reduced precursors of μ-conotoxins (A) and folded crude mixtures after 1 h (B), with the main product marked with an asterisk. Respective electrospray ionization (ESI) mass spectra of the oxidized peptides are shown in (C).

The formation of one major product for μ-GIIIA and μ-SIIIA can be explained by a rapid collapse into the favored native fold with the disulfide connectivity C1–C4/C2–C5/C3–C6 as introduced earlier.[17] For μ-KIIIA, such a rapid collapse also results in one main product; however, the connectivity C1–C5/C2–C4/C3–C6 is preferred here, and the pattern C1–C4/C2–C5/C3–C6 is only present as a minor fraction. Small differences in the elution profiles of buffer-oxidized μ-KIIIA (the crude product), that is, product formation, might result from the differences in batch size and composition of the oxidation buffer compared to the results of Khoo et al.[25] In the case of μ-PIIIA and μ-SmIIIA, a different folding mechanism indicative of a slower rearrangement results in the formation of several isomers. This can be attributed to more diverse noncovalent interactions and electrostatic forces compared to μ-GIIIA and μ-SIIIA. Here, it was suggested that the native isomer accumulates via reshuffling of disulfide bonds during the folding process and is dependent on the thermodynamic stability of the isomer formed.[17,31] Although μ-SIIIA and μ-SmIIIA have a high sequence similarity (Table ), the higher number of basic residues in μ-SmIIIA (six Arg) might cause the formation of multiple isomers, whereas in case of μ-SIIIA (two Arg and one Lys), only one isomer is preferred.[17,32] μ-PIIIA forms multiple isomers possibly for the same reason.[9] In contrast, the structure of μ-GIIIA tolerates a high number of basic residues and forms only one major product compared to that of μ-PIIIA and μ-SmIIIA.[17] Several aspects influencing the folding, such as the number of hydroxyproline residues, the loop size between the two linked cysteines, or amidation of the C-terminus, are discussed controversially in the literature without a clear preference indicating the uniqueness of each sequence and the respective biological activity.[17,28,33,34]

Conformational Analysis Using Molecular Dynamics

The solution NMR structures used herein as initial conformations for starting MD simulations were derived from the ConoServer database (Figure , Table ).[4,7] Comparison, alignment, and structural differences of the μ-conotoxin structures were discussed previously by Yao et al. and Tietze et al.[9,28] All μ-conotoxins, except μ-KIIIA, possess conserved structures that align significantly better in the C-terminal part compared to the N-terminal region.[28] A high similarity was seen for the backbone conformations between the loop 2 and loop 3 regions of μ-KIIIA, μ-SIIIA, and μ-SmIIIA, which have a higher selectivity for blocking NaV1.2 over NaV1.4 channels (Table ). On the other hand, μ-GIIIA and μ-PIIIA, which prefer NaV1.4 over NaV1.2, superimpose well in the second loop between C2 and C5.[9] Besides contributing to the structural rigidity, the disulfide bridges cause cysteine residues to form a hydrophobic core, enveloped by other charged and hydrophilic residues (Figure ). This hydrophobic effect plays a key role in the stability of the native fold.[35]

With three disulfide bridges present, there are six ways to open them one by one. The order of disulfide bond formation during the synthesis (and also in vivo) is not known, but as mentioned before, the process appears to be guided by thermodynamic aspects. Regarding the simulation strategy, it was decided to open the longest disulfide bridge first (bridging the longest sequence in between, see Table S1) as it is expected to instantly introduce the highest level of flexibility into the peptide backbone. The intention for this opening strategy was to increase the conformational entropy of the reduced version. In the case of μ-KIIIA (non-native connectivity), the shortest disulfide bridge was opened first. This would serve as a means to gauge the effect of loop size in retaining a stable structure. The loops in between the three disulfide bonds differ in size (Table S1) and thus needed to be considered in the analysis process.

It was observed during the simulations at room temperature that all peptides retained their initial conformation as demonstrated by the root-mean-square deviation (RMSD) of Cα atoms compared to the chosen starting NMR structure of each peptide (Table ) with μ-KIIIA possessing the lowest RMSD of 1.1 Å and μ-SmIIIA the highest of 2.7 Å (Figure S1, Table S2). The RMSD of the Cα atoms, the root mean square of fluctuation (RMSF) of all atoms of each residue, and the radius of gyration (Rg) of the whole protein were computed from each simulation for all of the five peptides (Figure S1, Table S2).

The native structure of μ-GIIIA displayed a hydrophobic core formed by the cysteines and a salt bridge between R1 and D12. As the first disulfide bond C2–C5 was removed, the distance between the now unbound cysteines increased moderately. However, this did not affect the overall three-dimensional conformation of the backbone proved by only a 0.3 Å increase in the backbone RMSD, but it could be observed that its RMSD progression is inconsistent in the simulation especially between 30 and 50 ns. The 300 ns simulation showed that the peptide had a stable backbone indicated by a relatively unwavering RMSD curve between 100 and 300 ns of simulation time as shown in Figure .

Figure 3

RMSD plots of disulfide bond opened versions of the five μ-conotoxins. The comparison of backbone stability between the peptides with the C2–C5 disulfide bond removed (black) and both the C2–C5 and C3–C6 bridges removed (red) between 100 and 300 ns of simulation time: (A) μ-GIIIA, (B) μ-KIIIA, (C) μ-PIIIA, (D) μ-SIIIA, and (E) μ-SmIIIA. Above the plots is a representation of two cases of disulfide connectivity discussed. Here, red represents the version with the single C1–C4 disulfide bond and black represents the C1–C4/C3–C6 disulfide connectivity.

Importantly, the residues significant for bioactivity experienced none to a very minimal increase in fluctuations (Figure , Table S4). Upon subsequent removal of the second disulfide bridge between C3–C6, an obvious stretch in the overall shape of the peptide was observed (Figure S3). The peptide had an excessive flexibility and adopted close to completely unfolded conformations during different intervals in the simulation (Figure S3). Interestingly, the cysteines forming the C3–C6 bond moved much closer to each other than the C2–C5 cysteine residues in the two disulfide bond opened structure. This is shown by the decrease in the RMSD of the C2–C5 and C3–C6 removed peptide between 60 and 70 ns (Figure S3). In longer time scales (300 ns), it was observed that the RMSD values dropped close to the ones with just one disulfide bond (C2–C5) removed, indicating that μ-GIIIA tends to fold back to retain its preferred native conformation (Figure ). The Rg of μ-GIIIA followed an almost identical pattern of progression to the RMSD, peaking between 40 and 65 ns before falling back toward its initial values, representing an unfolding-(re)folding event. The fluctuations of the residues K8, K11, R13, K16, and R19, which were reported to be responsible for ion channel binding,[3] did not show a significant change compared to the structure with all disulfide bonds intact (Table S3). With two disulfide bridges removed, the C-terminus including R19 and K16 showed a significant increase in fluctuations. It has been reported that K16 has a low priority for biological activity. Moreover, the exchange of this residue increased the binding affinity compared to the native toxin.[36] Our findings suggest that if the C2–C5 disulfide bond alone was removed, the structure may adopt a conformation still representing a structure close to the native fold.

μ-KIIIA has the lowest RMSD of 1.1 Å on average with respect to the chosen starting structure among the five peptides in their natively folded form (Table S2). In 100 ns of simulation time, the removal of the C2–C4 disulfide bond, thereby altering the C4 residue, results in the loss of its helix and increases the RMSD by 0.7 Å. However, at the 300 ns time scale, conformations sampled by μ-KIIIA showed the reappearance of its native helix (Figure ). The remaining C1–C5 and C3–C6 disulfide bonds were sufficient to retain the backbone stability and conformation, respectively, of μ-KIIIA. The progression of RMSD and Rg for μ-KIIIA followed the same scheme as observed for μ-GIIIA, that is, with the structure possessing two reduced disulfides revealing the largest variation in the conformational flexibility. It was observed from the 300 ns simulations that unlike μ-GIIIA, the two disulfide-deficient versions of μ-KIIIA did not regain the backbone stability of its one disulfide bond removed counterpart as shown in Figure .

Figure 5

Grouping of the peptides based on the favorability of two disulfide bond stability. The C2–C5 disulfide-deficient conformations (cartoon representations colored to distinguish secondary structural elements) of the studied μ-conotoxins (A) GIIIA, (B) KIIIA, (C) SIIIA, (D) SmIIIA, and (E) PIIIA superimposed on their energy-minimized native structure (black cartoon representation). From the 300 ns trajectory, five conformations (one every 60 ns) have been used. The average RMSD of these conformations in comparison to the reference native structure is shown in Å. Besides displaying the regions of similarity and dissimilarity between the native and the C2–C5-deficient versions, the figure also provides a grouping for the five peptides in terms of favorability of the disulfide-deficient version retaining structural characteristics of the native peptide (based on RMSD, RMSF, and Rg).

Although the average RMSD of native μ-PIIIA was a decent 2.3 Å compared to the selected starting structure over the course of the simulation, the backbone was constantly subjected to changes as can be seen from the RMSD plot (Figure S1). The removal of the C2–C5 disulfide bond resulted in a lesser fluctuating RMSD progression, although it came at the expense of a 1.5 Å increase in RMSD in the first 100 ns of simulation time. Conformations sampled by this peptide between 100 and 300 ns of simulation time showed the reappearance of its native helices (α-helix between O8 and S14 and 3_10 helix between L3 and C5). The structure of μ-PIIIA with both the C2–C5 and C3–C6 disulfide bonds removed showed the largest extent of structural variation among the five conopeptides with an average RMSD of 4.9 Å (Figure S5). The peptide did not tend to refold within 300 ns simulation time. With one disulfide bond removed, none of its functionally significant residues showed a pronounced increase in fluctuations, the highest of which was a 1.3 Å increase for R14 compared to the natively folded peptides. Unlike μ-GIIIA, the two disulfide-deficient versions of μ-PIIIA did not regain the backbone stability of its one disulfide bond removed counterpart. Meanwhile, the single disulfide-deficient version of μ-PIIIA adopted a very stable conformation, with the RMSD curve almost flat lining between 200 and 300 ns of the 300 ns simulation (Figure S3).

In contrast, almost all residues of μ-SIIIA displayed marginally higher RMSF values for the one and two bond-removed structures in comparison to the native structure containing all the three disulfide bonds. In the one disulfide-deficient version of μ-SIIIA, although the marginal increase in residue fluctuations was observed, the functionally significant restudies W12, R14, and H16 revealed only a minimal increase in the overall residue fluctuations. Another key residue R18, however, showed a larger mobility between the native and one disulfide-deficient versions. More importantly, the one disulfide-deficient version of μ-SIIIA retained its α-helical motif between K11 and H16, which is a significant aspect in terms of targeting sodium channels.[37]

We focused further on μ-GIIIA that forms a single oxidation product and μ-SmIIIA that forms multiple oxidation products during the synthesis to illustrate the phenomenon described above (Figure S3). The overall RMSD between the native connectivity and the structure with one disulfide bond removed was observed to be low in μ-GIIIA (2.1 Å) and high in μ-SmIIIA (2.9 Å) among the five peptides. Figure S3 shows average conformations for all the three 100 ns simulations of μ-GIIIA and μ-SmIIIA compared with their corresponding RMSD plots. An interesting observation from both peptides with an opened disulfide bond was the formation of new secondary structure elements that were not present in the native state. The structure of C2–C5-deficient μ-GIIIA achieved a reasonable equilibration between 200 and 300 ns of the simulation, and the inspection of the trajectory revealed that the conformations sampled by this peptide had a 3_10 helix between K16 and Q18 (Figure ). In μ-SmIIIA with one opened disulfide bond, an α-helix was formed between residues R13 and H18 (Figure ).

The stability of this helix through the entire course of the simulation can be accounted for by a combination of hydrogen bond formation and the presence of the bonded C15 in the center of the helix. In comparison with the native fold, the structure with the C2–C5 bond opened appeared well ordered. The distribution of the hydrogen bonds around the helix can be seen in Figure .

Figure 4

Comparison of μ-SmIIIA native fold and the structure with one disulfide bond opened. (A) Structure of the native peptide (completely oxidized, three disulfide bonds) with 100 conformations of its basic residues (red) superimposed with three residues marked as important for binding activity. The distribution of hydrogen bonds shows a sparse black area which indicates that the region surrounding it is relatively flexible. (B) Structure of the peptide containing one opened disulfide bond (C2–C5) showing (blue) its well-formed α-helix and the dense well-ordered hydrogen bonding illustrated as black cylinders. Higher rigidity inducing an improved structural stability of the peptide in (B) in comparison to the fully oxidized peptide in (A) is apparent from the reference RMSF plot (C). Despite the rigidity of the backbone, the orientations of the basic residues differ from the native structure.

The μ-SmIIIA structure with two disulfide bonds reduced formed a less stable 3_10 helix between residues R16 and H18, increasing the flexibility of the conformation. However, the functionality of a conotoxin is dependent on a stable backbone structure coupled with the favorable orientation of basic side chain residues for binding to their target and not solely on the flexibility of a distinct part of the peptide.[3] In μ-SmIIIA, though the 3_10 helix formation reduced the Rg of the peptide, the overall peptide structure drifted significantly from that of the native state, and the orientations adopted by the side chains of its basic residues varied largely when compared to either the native or the structure with one disulfide bond opened. Even with a single disulfide bond removed, the functionally significant R7 showed a moderate increase in fluctuations (Figure , Table S4).

The 300 ns simulations for the disulfide-deficient species of both μ-SIIIA and μ-SmIIIA revealed similar trends based on the observed RMSD progression (Figure ). In the C2–C5 disulfide removed version, between the two, μ-SmIIIA had a marginally more stable backbone than μ-SIIIA. However, with the C3–C6 disulfide bond also removed, both peptides had equally unstable backbones that did not show signs of refolding. Herewith, it is demonstrated how the individual disulfide bridges and the residues that occur within the loops affect the conformational stability on the basis of observed fluctuations in the backbone and the side chain residues.

A closer means of observing the folding behavior of the peptides in this study is to track the movement of the opened cysteine residues in the simulation. By tracking the distances spanned by the Sγ atoms of the cysteines, the tendency of the peptide to fold back to its original conformation or to explore completely new conformations can be identified. This behavior reflects on the underlying folding model that the peptide prefers to adopt. On the basis of this idea, we were able to find a clear correlation between the observations from the synthesis (Figure ) and simulation (Figure ). The peptides that preferred forming distinct main products (μ-GIIIA and μ-KIIIA) in the synthesis (Figure ) also exhibited their preference to fold back to their original conformations as seen from the C2–C5 Sγ distance profiles in Figure . This is an indicator to the preference of the BPTI-like folding pathway. On the other hand, μ-PIIIA and μ-SmIIIA that formed a mixture of products in the synthesis (Figure ) preferred moving away from their original conformations and exploring new conformations as indicated by the C2–C5 Sγ distance profiles in Figure . This reflects on the Hirudin-like folding pathway preferred by these peptides. The behavior of μ-SIIIA could be placed between these two extreme cases. The videos from the MD simulations where the distances between the opened CYS residues were tracked are provided for the cases representing the opposite extremes, namely, μ-GIIIA and μ-PIIIA are provided as a part of the Supporting Information.

Figure 6

Relation between the unbound cysteine distances and the underlying folding pathway. Top: the distances between the Sγ atoms of the cysteine residues from the peptides with a single disulfide opened plotted over the 300 ns simulation for (A) μ-GIIIA, (B) μ-KIIIA, (C) μ-PIIIA, (D) μ-SIIIA, and (E) μ-SmIIIA. Bottom: a schematic representation of the placement of the five μ-conotoxins within the established BPTI-like and hirudin-like models.[18] The NMR structures of hirudin (the gray cartoon) and BPTI (the blue cartoon) with their three disulfide bridges shown as yellow sticks. This classification is based on the observations from this study.

Conclusions

The folding of smaller disulfide-rich peptides and oligopeptides is a less well-understood folding event because of a much higher degree of flexibility and often a lower extent of structure-forming elements. MD simulations are increasingly used to assist the experimental work for understanding and predicting the folding process. Additionally, they have been routinely used in structure–activity relationship studies, drug discovery, and design pipelines.[38−42] Previous studies using the five μ-conotoxins investigated herein gave insights into the folding and binding modes adopted by these peptides.[9,25,28,32,36,43] Simulating the complete oxidative folding pathway following the formation of non-native disulfide intermediates until the native disulfide bonds are formed as reported by previous studies[44,45] using coarse-grained models is not within the scope of this study. Our work, however, aimed at determining how a particular disulfide bond contributes to the stability of the peptide. Consequently, this approach reviews the validity of the logic that the removal of a disulfide bridge, that is, herein C2–C5, represents a reduction of the backbone stability when considering RMSD and RMSF values.

With the C2–C5 disulfide bond removed, only μ-SmIIIA revealed a noticeable increase in the average fluctuations of its functionally significant residues. In addition, the fact that secondary structural elements such as α-helices were formed in some peptides containing only two disulfide bonds suggests that in distinct cases (e.g., μ-GIIIA), a greater structural rigidity of the backbone may be achieved if one disulfide bridge is removed. This helix-induced stability while strengthening the backbone might reduce the extent of overall fluctuations of the basic residues responsible for binding activity. The biological activity and selectivity of disulfide-deficient mutants might differ from the native conformation as shown for μ-GIIIA recently. Fifteen different disulfide isomers are possible for a peptide containing six cysteines, and still three different isomers (ribbon, bead, and globular) might occur in case of four cysteines. It has, however, not been mentioned in the report by Han et al.[46] which isomer of the μ-GIIIA analogs has been tested because the structural characterization of the respective products was not performed. Apart from the reports on μ-GIIIA regarding the disulfide-deficient variants, another study by Khoo et al.[47] provides an insight into the removal of disulfide bridge C1–C9 in μ-KIIIA, resulting in only a minimal change in the biological activity against NaV1.2 and NaV1.4. In contrast, there are no experimental data for the disulfide-deficient species of μ-conotoxins PIIIA, SIIIA, and SmIIIA available so far.

With respect to drug design and synthesis, the simplification to two disulfide bonds would be a clear benefit for disulfide-rich peptides and proteins. A similar study by Yu et al. on α-conotoxin cVc1.1 complements our idea of the reduction of the number of disulfide bonds.[49] In this respect, we can conclude from our MD simulations that two disulfide bridges could be sufficient to maintain a stable backbone for the majority of the μ-conotoxins studied. However, it is important that the deficient structure is sufficiently supported by at least one pair of cross-linked disulfide bridges that span to almost either ends of the sequence. From the results obtained, a rank order of the five peptides can be provided: μ-GIIIA and μ-KIIIA fall in the highly favorable category, μ-SIIIA falls in the moderately favorable group, and μ-PIIIA and μ-SmIIIA fall into the least favorable group (Figure ). We also conclude that the C3–C6 disulfide bridge plays the greatest role in retaining the backbone stability for the current strategy of disulfide bond removal employed.

As expected, the removal of two disulfide bridges led to an increased backbone flexibility, the formation of a series of intermediate conformations, and a less stable peptide. The increased fluctuations of basic side chain residues responsible for the interactions with the sodium ion channels are in this case unlikely to stay in favorable orientations for binding. Furthermore, the formation of different disulfide isomers for μ-KIIIA, μ-PIIIA, and μ-SmIIIA in the experimental self-folding approach indicates a difference in their folding behavior,[9,17,25,29] which cannot be explained unambiguously by simulations and disulfide removal. Finally, the observations of the structural stability of the backbone observed and the extent to which the peptide tries to fold back gives us clues on which disulfide folding model a particular peptide tends to likely prefer. It must be said that this distinction is still not completely in black and white but can be ranked or ordered relative to each other between the extreme models (BPTI and hirudin). Our attempt to classify the peptides between the BPTI and hirudin folding models is illustrated in Figure . The correlation of observations from the experiment and simulation as discussed earlier in the results adds further validity to this classification proposed in Figures and 6. Thus, we have been able to demonstrate the usefulness of molecular simulations in applications beyond the observation of the structural behavior of a peptide in solution to being used as a tool for the generalized assignment of peptides to established folding models.

Materials and Methods

MD Simulations

Disulfide bonds were systematically removed to yield a partially folded conformation as the starting structure for the simulation. From an atomic perspective, this translates to changing the bonded cysteines to nonbonded cysteines by protonating the sulfur atoms. This step was done by pdb2gmx program within the GROMACS 5.1.4 package.[50,51] This approach induces the least changes in the coordinate file to create a disulfide-deficient species as opposed to the usually performed mutation studies where the Cys residues are replaced by the Ala or Ser residues. It further eliminates the errors arising from the manual manipulation of the coordinate file. Four out of five μ-conotoxins in this study possess the native disulfide connectivity (C1–C4/C2–C5/C3–C6), whereas μ-KIIIA adopts a C1–C5/C2–C4/C3–C6 connection as the stable conformation.[25] The following structures were used herein: PDB ID 1TCG (μ-GIIIA),[22] PDB ID 2LXG (μ-KIIIA),[25] S00159 (μ-PIIIA),[9] BMRB 20025 (μ-SIIIA),[28] and PDB ID 1Q2J (μ-SmIIIA).[29] The NMR ensembles of μ-GIIIA, μ-KIIIA, μ-SIIIA, and μ-SmIIIA had 20 structures and μ-PIIIA had 15 structures in the respective PDB file. The first model was selected as the best representative structure except for SmIIIA, whereas model 13 was chosen as mentioned in the PDB file as the best model for this peptide.

GROMACS 5.1.4[50,51] was used for all the MD simulations in this study. An individual peptide was placed in the center of a cubic box that evolved to a final volume of 2.5 × 2.5 × 2.5 nm3. TIP3P[52] water model was used as the solvent to fill the box. Appropriate amounts of Cl– ions were added to balance the positive charge of the μ-conotoxins. Simulations were run using the AMBER99SB-ILDN[53] force field, which was chosen based on its better agreement with the NMR data and an accurate modeling of helical proteins in comparative studies.[54−56] In the process of preparing the peptide for the production MD simulation, energy minimization simulations were with 10 000 steps of the steepest descents minimization protocol and convergence reached when the maximum force on any atom is no greater than 100 kJ/mol/nm. A thermal equilibration in the NVT ensemble at 300 K using the velocity-rescaling modified Berendsen thermostat[57] and a constant pressure equilibration in the NPT ensemble using the Parrinello–Rahman barostat[58,59] at 1 atm were carried out for 20 ns each, prior to production MD. During both the temperature and the pressure ensemble simulations, positional restraints on the peptides were applied using the LINCS[60] algorithm. Each peptide was subjected to three production runs, and on each run, preprocessing and equilibration were performed independently. First, the conformation with all the three disulfide bonds was considered as the control simulation. In the second simulation, the link between C2–C5 (C2–C4 for KIIIA) was removed leaving the two other disulfide bonds intact. In the third simulation, both the C2–C5 (C2–C4 for KIIIA) and the C3–C6 disulfide bonds were removed, leaving the peptide constrained originally by three disulfide bonds now supported by a single disulfide bonds. The production MD was done for 100 ns for all the peptides, and for instances with the opening of a single or two disulfide bonds, a total of 300 ns of simulations were carried out by running an extended 200 ns simulation from the final checkpoint of the 100 ns trajectory. The extended 200 ns simulations were conducted to allow sufficient sampling for the observation of events from the refolding process such as the possibility of the unbound cysteines coming close to each other. All simulations were conducted with a 2 fs time step and data written to the logs and trajectory at every 5 ps. Periodic boundary conditions were applied to the system. Long-range electrostatics were accounted by the particle mesh Ewald method.[61,62] For every 100 ns of simulation, 20 000 frames were written to the trajectory. The effect of periodic boundary conditions was adjusted by the suppression of center of mass movement from the trajectory prior to analysis. Visualizations of conformations for the analysis and creation of images were performed using VMD.[63] The RMSD, RMSF, and Rg plots were created using the program Grace (version 5.1.25), whereas the distances between the unbound cysteines were plotted using the tools within VMD. The RMSD curves were plotted for every 10 ps (10 000 frames), whereas the distance curves were plotted for every 5 ps (20 000 frames).

Chemical Synthesis and Purification of μ-Conotoxins

Peptides were produced by an automated solid-phase peptide synthesis using a standard Fmoc (N-(9-fluorenyl) methoxycarbonyl)-protocol and an EPS 221 peptide synthesizer (Intavis Bioanalytical Instruments AG, Cologne, Germany) as described earlier and purified by preparative reversed-phase (RP) HPLC (Shimadzu LC-8A system, Duisburg, Germany). The gradient used was 0–50% eluent B in 120 min with 0.1% trifluoroacetic acid (TFA) in water (eluent A) and 0.1% TFA in acetonitrile/water (9:1) (eluent B) on a C18 column (Knauer Eurosphere 100, Berlin, Germany) with the dimensions 50 mm × 300 mm (5 mm particle size, 100 Å pore size). Reduced and oxidized peptides were analyzed on a LC–ESI micrOTOF-Q III mass spectrometer (Bruker Daltonics GmbH, Bremen, Germany) coupled with Dionex Ultimate 3000 (Thermo Scientific, Dreieich, Germany) equipped with a EC100/2 Nucleoshell RP18 Gravity 2.7 μm column (Macherey-Nagel, Düren, Germany). Analysis of the MS data was performed using Bruker Compass Data Analysis 4.1. The LC conditions used were as follows: eluent A was water with 0.1% acetic acid, whereas eluent B was acetonitrile containing 0.1% acetic acid. A gradient of 0–60% of eluent B in 12 min was used, and detection was at 220 nm.

Oxidation of Reduced μ-Conotoxin Precursors

Oxidative folding of the linear μ-conotoxins GIIIA, KIIIA, PIIIA, SIIIA, and SmIIIA in a buffer system containing redox agents was performed as described earlier. Each μ-conotoxin (1 mg) was subjected to oxidation, and fractions of the reaction mixture were monitored over time by RP HPLC using a Shimadzu LC-10AT system (Duisburg, Germany) equipped with a C18 column (Vydac 218TP54, Worms, Germany, 4.6 mm × 25 mm, 5 mm particle size, 300 Å pore size) and the gradient 0–60% eluent B in 60 min with 0.1% TFA in water (eluent A) and 0.1% TFA in acetonitrile (eluent B). Reaction control was performed over a time period of 24 h, and oxidation was stopped by adding 1% TFA in water. Monitoring revealed that the oxidation reactions were completed within the first 60 min of the reaction time. The fractions were collected for each peptide and subjected to LC–ESI mass spectrometry for the confirmation of the molar mass corresponding to the oxidized products.