Homology modeling and molecular dynamics simulation of odonthubuthus doriae (Od1) scorpion toxin in comparison to the BmK M1.

All of the α-subgroups share similarity in their sequence and structure but different in the toxicity to various voltage-gated sodium channels (VGSCs). We modeled the first 3D structural model of the Od1 based on BmK M1 using homology modeling. The reliability of model for more investigation and compare to BmK M1 has been examined and confirmed. Then the model structure is further refined by energy minimization and molecular dynamics methods. The purpose of this modeling and simulation is comparison toxicity of two mentioned toxins by investigation structural feature of functional regions including core domain, 5-turn and C-terminal which make NC domain. In the one hand, it is intriguing that Od1 in comparison to BmK M1 shows same solvent accessible surface area (SASA) in 5-turn region but a little more exposed and feasibility (more SASA) in C-terminal region and key functional residues of C-terminal such as positive residues Arg58, lys62 and Arg (His)64. These data suggested that Od1 has similarity with BmK M1 but has more toxicity to sodium channel. In the other hand 5-turn proximity of C-terminal to 5-turn in BmK M1with cis peptide bond is less than Od1 without cis peptide bond which is a confirmation with experimental data about BmK M1.A better understanding of the 3-D structure of Od1and comparison to BmK M1 will be helpful for more investigation of functional characters action of natural toxins with a specialized role for VGSCs.

Background

α-like toxins are scorpion toxin induce a prolongation of the action potential of excitable cells such as nerves and muscle by slowing down the inactivation of the voltage gated sodium channel (VGSC). Most α-like toxins are single chain composed of 60-70 amino acid cross-linked by 4-disulfide bridges [1, 2]. In the field of Iranian scorpion toxicology, the Iranian scorpion Odonthubuthus doriae has received mark attention with reference to isolation and pharmacologically characterization of its venom toxins. A number of peptide toxins from its venom have been reported. They are Od1, an α-like toxin [3], OdK1 [4], and OdK2 [5] with an effect on voltage gated potassium channels. The pharmacological effect of Od1 showed that the inactivation process of the insect channel para/tipe, was severely hampered by 200nm Od1 (EC50=80) ±14nm), whereas Nav1.2 /β1, was not affected on concentration up to 5µm. Nav1.5/β1 was influenced only at micromolar concentrations up to 5µm [3]. The potent modulation of the VGSCc Nav1.7 was seen by Od1 at low nanomolecular concentrations. The inactivation of Nav1.7 is dramatically impaired with an EC50 value of 4.5nm. The pharmacological evidences show that Od1 is a high affinity ligand with a unique selectivity pattern on Nav1.7 which found mainly in peripheral sensory and sympathetic neurons and plays a critical role in short-term and inflammatory pain (Wood et al., 2007). This selective high affinity modulation of Nav1.7 currents is seen with well-known VGSC modulators, local anesthetic and class I antiarrhythmic drugs such as lidocain in a frequency dependent matter with an EC50 value of 450µm [6]. The isolated analgesic peptide from scorpion venom such as BmKITAP, an excitatory insectselective toxin [7]; BmkdITAP3, a depressant insect-selective [8], BmKAGAP, an antitumor analgesic peptide showing inhibitory effect on visceral and somatic pain [9] BMkAngp1 [8]; BmkAngM1[10]; BMKAs [11]; BMKIT2 [12] could target VGSCs by virtue of their primary structure homology and similar scaffold among long-chain sodium channel toxins [13]. Od1 exhibits different and specific sensitive for insect or mammalian VGSCs unlike other pharmacologically studied scorpion toxins such as Lqh-3 and Lqh-2 [14, 15], in addition to the tarantula peptides protx-1 and protx-11 [16]. In the present study, the three-dimensional (3D) structure of Od1 was identified based on the available template structural homologues of BmK M1[17], and the model validation with standard parameters [18, 19]. This study may useful in further functional characters action of natural toxins with a specialized role for VGSCs involved in pain pathway and inflammatory pain.

Methodology

Homology modeling of OD1:

The protein sequence of Od1 was previously indicated but there is no available 3-D structure of Od1. We utilized MODELLER [20] for modeling Od1. Sequence similarity of known structures to Od1 has been searched using profile build of MODELLER. The result of sequence alignment showed homologues sequence with identity range between 60- 75%, and selection of best model has been done based on parameters including lowest value of the DOPE assessment score together with the highest GA341 assessment score. BmK M1 (PDB ID: 1DJT) [17] – an α- like toxin and in same family with Od1- with sequence similarity and e-value about 75% and 0 respectively, was selected as the template. Final 3-D structure of Od1 model has been exhibited in the (Figure 1b). Also Sequence alignment of template and target was calculated with CULSTALW (Figure 1a) [21]. To assess reliability of model, we executed online VERIFY 3D (http://nihserver.mbi.ucla.edu/SAVS/) and PROCHECK V 3.4.4 software [18].

Figure 1

(A) alignment of two toxins using CLUSTALW. A line above the alignment is used to mark strongly conserved positions. Three marks (“*”, “:” and “.”) are used: “*” mark shows fully conserved residue, “:” mark shows ‘strong’ groups is fully conserved and “.” mark shows ‘weaker’ groups is fully conserved; (B) Final 3-D structure of the of Od1 model. 5 turn (8 to 12) has been shown with red turn and C-ter segment (residues 58 to 64) with blue coil which these two regions make NC domain (both blue coil and red turn). Plots for checking reliabity of homology model; (C) This figure illustrates reliability of homology model provided that residues in homology model (black line) like template (gray line) show score greater than 0.2,; (D) Ramachandran plot of homology model of Od1 shows good stereo-chemical quality of the model. Most favored regions [A, B, L], Additional allowed regions [a, b, l, p], generously allowed regions [~a, ~b, ~l, ~P], another regions are disallowed regions. Glycin residues have been shown as triangles.

Simulation Protocol:

The molecular dynamics simulation was carried out for both toxins, homology model of Od1 and BmK M1.The system was neutralized by 1Cl- and 2 Cl- respectively, we utilized cutoff of 12 Å for short-range non-bonded interactions. Energy minimization was carried out in 10000 step in 20ps.All simulation was performed at NPT at 310k and a pressure of 1 atm and set 2fs for time step. The MD simulation was performed for 10ns using NAMD Version 2.8 [22] and Charmm 27 force field [23] and the model visualized by visual molecular dynamic (VMD) [24].

Result

Homology model:

Analysis of the model Verify 3-D score [19] was depicted for the template and target proteins, which has been used for assessment of the reliability of the structure of the model. The (Figure 2a) showed that 100% of residues in homology model like template show score greater than 0.2, this imply that this model is reliable model for more investigation (The percentage of residues with a score > 0.2 should be more than 80% for a reliable model) [19]. Furthermore we plot Ramachandran plot of the homology model by using PROCHECK [18] as a further test of model, we found that 90.3% of the residues placed in most favored regions and 9.3% of amino acids placed in additional allowed regions and there is no residues (except Gly) in disallowed residues. In general, a score close to 100% implies good stereo-chemical quality of the model (Figure 2b).

Figure 2

(a) Root mean square deviation (RMSD). RMSD plot of two toxins which black line is Od1 and gray line is BmKM1; (b) Distance between 5-turn and C-ter and black line is Od1 and gray line is BmKM1.

General structural feature of Od1 modeling:

Since all toxins have similar secondary structure but different in toxicity and selectivity function, assignment of secondary structure of Od1 was evaluated by STRIDE [25] as implement in VMD and contains 3-beta sheets and alpha helix in βαββ motif that residues 2-9, 34-40, 43- 52 form 3 β-sheets anti parallel respectively and alpha helices composed of residues 20-29 with coils and turns that connect these regions (Figure 1b). Also Superimpose of Od1 and BmK M1 in Figure 3c show similarity structure of them.

Figure 3

Solvent accessible surface area (SASA). (A) SASA of Cterminal (residue 57 to 64); (B) SASA of 5turn (residue 8 to 12) are presented and in both of the plot black line is Od1 and gray line is BmKM1; (C) Superimposition illustration of OD1 (the red ribbon) and BMKM1 (the blue ribbon) along with SASA of residue 64 using VMD is exhibited. The red surface is Arg64 in Od1 but the blue surface is His 64 in BmKM1.

Molecular dynamic simulation:

RMSD The root-mean-square deviation (RMSD) is the measurement of the average distance between the atoms which is calculated for protein structure to check its stability and we measured RMSD for all atoms of the two proteins backbone during 10 ns simulations (Figure 3a). Square deviation (SD) of protein reached 0.2 and 0.1 Å respectively, which shows stability of two proteins during simulation and are equilibrium for more investigation.

Functional site:

The two structural regions of the alpha-like scorpion toxins consisting of the NC domain which composed reverse turn (residue 8-12) in the combination with the C-tail (57-65) and core domain (residues: 40-42) have notable attention in the binding preference for the distinct VGSCs (Figure 1b) [26]. The five residue of the reverse turn in the scorpion toxin consist of some conserved and variable residues in 3 subgroups α-toxins. In the our study, alignment of these two α-like toxins shows that three residues in position 8, 9 and 10 are different and more acidic in Od1 including Asp (Figure 1a). Since it has been shown that substitutions at the five-residue turn and stabilization of the 9-10 bond in the cis conformation show no significant effect on mammalian and insect sodium channel [26-28]. Consequently we can suggest that these two acidic residues make little difference between these two α-like toxins. Granier et al found that antigenic segments of scorpion toxins are located at exposed regions of the molecular surface [29] so measurement of the solvent accessible surface area (SASA) of functional region is an important comparison structure-function analysis between Od1 and BmK M1. For this SASA (Figure 3a & Figure 3b) of these two functional regions has been evaluated and it show that SASA of 5 turn is similar in both toxins but SASA of C-terminal (58 to 64) of Od1 has a little more exposure rather than SASA of same region in BmK M1. Also It has been shown that positive residues in position 58, 62 and 64 in BmK M1strongly affect affinity of binding to sodium channels and these residues are more exposed and accessible in BmK M1 rather than buried as in the classical a-toxins [28-30]. Od1 similar BmK M1contains same residues in this position except position 64 which it has been occupied by Arg and His in Od1 and BmK M1 respectively (Figure 1a). So we measured solvent accessible surface area (SASA) of these residues after simulation, the data in the table 1 show that Arg is more exposed than His Table 1 (see supplementary material) (Figure 3c). Since the key functional residues such as 58 and 62 are same but SASA of them exhibit difference (Table 1). In the one hand, this imply that more exposing residues can interact stronger and better to sodium channels which may be more favorable than that of in BmK M1 and increased toxicity of Od1, in the other hand this feature demonstrate that Od1 as an α-like toxin has exposed Arg in position 58 .

NC domain at the experimental studies:

It has been shown that α-toxins have a unique tertiary structure in C-terminal (58-64) and 5-turn (8-12) which caused by disulfide bond between Cys12 and Cys 63, and constitute NC domain (Figure 1b). This region is functionally site for binding to sodium channel [31]. Existence cis peptide bond between Pro9 and His10 in BmK M1mediates the spatial relationship between the turn 8-12 and the C-terminal stretch 58-64 and cause a proximity of C-terminal and 5-turn which forms a local tertiary structure unique for this functional site [28, 31]. Distance between C-terminal and 5-turn is an important indicator for structure-function studies. We calculated distance between two regions during simulation using VMD DisRg plug-in [32], and plot of distance (Figure 3b) of two studied α- like toxins during 10 ns of simulation show that latter distance in Od1 is more than BmK M1maybe due to absence of cis peptide bond in the Od1.

Discussion

As showed and verified previously, two opposite molecular surface, face A and face B are conservative among BmK M1, M4, M8 and Aah2, even along with well known difference in their primary sequence. These two faces in OD1 are same with described toxins. This similarity may highlight the conserved general conformation of these two faces among alpha-toxins and alpha-like toxins groups. The hydrophobic residue of face A consist Tyr 35, Trp47, Trp38 and Tyr5 together with face B residues Tyr21, 42 and14 are same in the alpha-like and alphatoxins. SASA deformation of residues Tyr42 of OD1 showed significant difference with BmK M1.These residues between OD1 and BmK M1 are identical and their SASA are same only in Tyr42 located in functional core domain of α-like toxin (core domain: residues 40-43) (Table 1). Since Trp38 and Tyr42 (located in the β2-sheet and loop between β2-sheet and β3-sheet) are involved in the pharmacological effect of the toxin, so this position maybe affect the difference binding and affinity of α- like toxins OD1 and BmK M1[31, 33]. Arg58 has a key role for binding to VGSCs. The accessibility and feasibility of this position has been shown to be important for pharmacological characteristics [28, 29]. However this position is identical in both toxins. Their SASA of Arg is 87.5 and 37.66 in OD1 and BmK M1 respectively. This difference may responsible for two pharmacological characteristics profile of alpha-like toxins Table 1 (see supplementary material).

Conclusion

The α-like toxin, Od1 has two specialized domains in its structure that are responsible for binding to S3-S4 of IV domain of sodium channel 1.5 and 1.7 called NC domain and core domain. In this paper the first, we modeled Od1 using homology modeling based on another α-like toxin BmK M1 and then, MD simulation was carried out for two toxins to refine model and comparison structural feature and toxicity of them. The presented model will be potentiate and facilitate structural and functional investigations to use these two toxins in targeted drug delivery.