Cardiotoxins: Functional Role of Local Conformational Changes

Article Cite This: J. Chem. Inf. Model. 2017, 57, 2799-2810

pubs.acs.org/jcim

Cardiotoxins: Functional Role of Local Conformational Changes Anastasia G. Konshina,† Nikolay A. Krylov,†,‡ and Roman G. Efremov*,† †

Shemyakin−Ovchinnikov Institute of Bioorganic Chemistry, Russian Academy of Sciences, 16/10 Miklukho-Maklaya str., 117997 GSP, Moscow V-437, Russia ‡ Joint Supercomputer Center, Russian Academy of Sciences, Leninsky prospect, 32a, Moscow 119991, Russia S Supporting Information *

ABSTRACT: Cardiotoxins (CTs) from snake venoms are a family of homologous highly basic proteins that have extended hydrophobic patterns on their molecular surfaces. CTs are folded into three β-structured loops stabilized by four disulfide bridges. Being well-structured in aqueous solution, most of these proteins are membrane-active, although the exact molecular mechanisms of CT-induced cell damage are still poorly understood. To elucidate the structure−function relationships in CTs, a detailed knowledge of their spatial organization and local conformational dynamics is required. Protein domain motions can be either derived from a set of experimental structures or generated via molecular dynamics (MD). At the same time, traditional clustering algorithms in the Cartesian coordinate space often fail to properly take into account the local large-scale dihedral angle transitions that occur in MD simulations. This is because such perturbations are usually offset by changes in the neighboring dihedrals, thus preserving the overall protein fold. States with a “locally perturbed” backbone were found in experimental 3D models of some globular proteins and have been shown to be functionally meaningful. In this work, the possibility of large-scale dihedral angle transitions in the course of long-term MD in explicit water was explored for three CTs with different membrane activities: CT 1, 2 (Naja oxiana) and CT A3 (Naja atra). Analysis of the MD-derived distributions of backbone torsion angles revealed several important common and specific features in the structural/dynamic behavior of these proteins. First, large-amplitude transitions were detected in some residues located in the functionally important loop I region. The K5/L6 pair of residues was found to induce a perturbation of the hydrophobic patterns on the molecular surface of CTsreversible breaking of a large nonpolar zone (“bottom”) into two smaller ones and their subsequent association. Second, the characteristic sizes of these patterns perfectly coincided with the dimensions of the nonpolar zones on the surfaces of model two-component (zwitterionic/anionic) membranes. Taken together with experimental data on the CT-induced leakage of fluorescent dye from such membranes, these results allowed us to formulate a two-stage mechanism of CT−membrane binding. The principal finding of this study is that even local conformational dynamics of CTs can seriously affect their functional activity via a tuning of the membrane binding site − specific “hot spots” (like the K5/L6 pair) in the protein structure.

■

into a continuous hydrophobic “bottom”was defined.13−15 In addition, based on MD simulations of CT 1 (Naja naja) with an explicit model membrane, two other erythrocyte membraneactive sites containing positively charged residues were proposed.16 At the same time, despite the large body of accumulated data, the exact molecular mechanisms of CTinduced cell damage remain unclear. This is caused by difficulties in studying protein−membrane interactions using experimental techniques. In this situation, computer simulations represent a prospective alternative, given that atomistic modeling is capable of providing a structural/dynamic picture of CTs binding to membranes. MD simulation is the most common and powerful tool for this.17 It should be noted that only a moderate number of works devoted to MD simulations of CTs are known. Most of them have focused on the modeling of the toxins’ interaction

INTRODUCTION Cardiotoxins (CTs) are 60-amino acid length excretory proteins with a rigid β-structured core. CTs belong to the family of highly basic membrane active protein toxins from snake venom.1,2 All CTs have similar three-dimensional (3D) structures adopting a “three-fingered” fold stabilized by four disulfide bonds.3,4 The “fingers” of CT loops I−III are formed by strands of antiparallel β-structures. At present, CTs are known to be cytotoxic to various types of cells.5−12 Unlike neurotoxins and other representatives of the three-fingered proteins specifically interacting with membrane receptors, the cell membrane per se is the most likely target for CTs. The interaction of the latter with lipid membranes (and membrane mimetics like micelles and liposomes) as well as CT-induced membrane damage followed by its lysis (reported by dye leakage measurements) were clearly demonstrated in a number of nuclear magnetic resonance (NMR) and fluorescent spectroscopy studies.13 The principal membrane-binding motif of CTsmostly apolar tips of loops I−III assembled © 2017 American Chemical Society

Received: June 28, 2017 Published: October 20, 2017 2799

DOI: 10.1021/acs.jcim.7b00395 J. Chem. Inf. Model. 2017, 57, 2799−2810

Article

Journal of Chemical Information and Modeling with implicit18 or explicit model membranes16,19 without detailed analysis of the observed conformational states of proteins. Short-length (nanoseconds) MD simulations of CTs in water were performed studying toxin behavior in different environments (water, detergent micelles),20 the existence of tightly bound waters in the loop II region21 and the unfolding stability of CTs at various temperatures.22 Deciphering the structure−function relationships for CTs is especially important and interesting since these proteins have a very similar and stable spatial structure, although their activity can seriously vary. In contrast to many other membrane-targeting polypeptides (e.g., linear antimicrobials23) CTs do not undergo significant conformational rearrangements upon binding to membranes. 1 Therefore, it can be assumed that the functionality of the toxins is determined by certain fine-tuned and (probably) local aspects of their dynamics in solvent. The conformational transitions in proteins associated with their functional activity are especially interesting and important. Processes of biomolecular recognitionfirst of all ligand− receptor interactions and enzyme catalysisassume the existence of functional motions of a biomolecule and ensemble of its conformations.24−26 Single or just a few experimental structures (X-ray or NMR models) are usually not sufficient to describe the functionally important dynamic properties of a protein molecule. In most cases, more detailed information is required, including rather small conformational changes, which are limited by one or two residues. The functional meaning of such a local conformational dynamics and the methods to detect it are among issues that can hardly be captured in an experiment. At the same time, short-lived and/or highly flexible transition states can be one of the key elements required to understand the molecular mechanisms of protein function. Recent advances in the development of conformational sampling algorithms and related software have made it possible to study protein motions on micro- and millisecond time scales where many critical biochemical processes occur.17,25,27,28 To probe the protein conformational heterogeneity observed in such simulations, one needs to extract from MD trajectories a rational set of clusters with defined structural characteristics. A straightforward rationalization of MD data is often difficult to obtain because of the high dimensionality of the protein conformational space. The greater the amplitude of protein domain movements, the easier it is to allocate its different conformational states. It is quite straightforward to determine the mutual rearrangement of large protein domains and hinge residues responsible for the movements with standard techniques such as root-mean-squared coordinate deviations (RMSD) clustering algorithms29 and principal component analysis (PCA).30,31 The task becomes less trivial in the case of a large number of small-amplitude changes as well as highly flexible unfolded proteins. Among the variety of commonly used pairwise distance metric clustering algorithms there is no “best” one and the clustering results depend to a large degree on MD data, atom sets for the pairwise comparison, distance cutoff, and other input parameters.29,32 In view of a great number of low-amplitude movements observed in atomistic MD simulations, many clusters may be heterogeneous, containing models with different types of conformational behavior inside the same group.33 But lowering the RMSD cutoff value strongly increases the number of clusters and complicates the analysis. Clustering in the backbone dihedral angles space instead of Cartesian coordinates was successfully applied to achieve a

more accurate discrimination between protein conformations.33 Also, PCA in dihedral angle space is considered advantageous due to the efficient separation of internal and global (translational and rotational) motions, and as a result, the better exploration of the free energy landscape of a protein, especially in the case of intrinsically disordered proteins. For the latter, the fitting procedure is nontrivial or even impossible due to the significant structure rearrangements of the molecules.34,35 Of particular interest are cases when significant perturbations of dihedral angles are not accompanied by large conformational changes expressed in terms of RMSD values. Some combinations of the changed backbone dihedrals can also result in very similar structures in the Cartesian space. In rigid and well-structured proteins, like CTs, thermal fluctuations of not only adjacent backbone dihedrals but also of more distant torsions along the loop chain are interdependent and correlated. This so-called “latent dynamics” is supposed to prevent large steric clashes and distortion of the secondary structure elements but it can be unimportant in a functional sense.36 At the same time, analysis of significant changes of backbone dihedrals was applied to a data set of 459 high-resolution pairs of protein structures.37 As a result, large-amplitude angle transitions were detected in many proteins. It was found that some local protein motions revealing subtle RMSDs, may have high-amplitude dihedral transitions and these point structural changes are involved in intermolecular interactions. The peptide plane flip represents an extreme case of compensated large-amplitude changes of ψ and φ angles of two adjacent residues so that the intervening peptide plane flips upside down with little effect on the orientation of the flanking regions including side chains.38 Transitions of a β-sheet structure via peptide plane flips of the mature amyloid protein to an α-sheet conformationa proposed amyloidogenic prefibrillar intermediatewere detected in MD simulations.39,40 This is a rare published example of such angle transitions. Also, the peptide plane flips were detected earlier in loop regions of human lysozyme upon its MD studies in water.41 Taking all of the above into account, it is reasonable to assume that the backbone dihedral angle analysis of MD data is well-suited for uncovering potentially important local motions with a high degree of dihedral angle perturbations in CTs. Are there conformational states of CTs associated with significant dihedral angle transitions that can affect their membrane-binding activity? Is it possible to determine some flexibility pattern common to all CT representatives that are in general rigid and highly stable molecules? In the present work, these questions were addressed by means of computational simulations. For this purpose, three homologous CTs: CT 1, CT 2 (Naja oxiana), and CT A3 (Naja atra) with different biological activities were simulated in water using the modern all-atom force field AMBERff99SB-ILDN.42 Two independent long-term MD simulations were performed for each of the toxins and the general reproducibility of its MD behavior was proved. Based on the backbone dihedral (φ/ψ) angle distribution analysis of MD configurations, amino acid residues of CTs were mapped according to their structural/dynamic parameters. Stable conformational states that are common for all three CTs but differ from their starting structures were discovered. One of them is characterized by large-amplitude angle changes in the K5/L6 pair of residues, thus leading to serious reorganization of 2800


Article

Journal of Chemical Information and Modeling

Gaussian curves similarly to Watanabe et al.53 Two essentially different variants of such distributions were found: two or more well-separated (P ∼ 0) data sets or single continuous set of data (P > 0). In the latter case, multiple forms of the curve shape were detected. Between the two extreme cases: a simple Gausslike P function and a curve with some well-resolved peaks, various intermediate variants with more or less strongly overlapping peaks were also elucidated. The scheme illustrating classification of residues based on the analysis of their dihedrals distributions is shown in Figure S1. For many angles, a broad band with one prominent peak (Figure 1A: φ-7, Figure S2A: φ-

the functionally important hydrophobic pattern (the so-called “bottom”) on the molecular surface of CTs−its reversible breaking into two smaller ones and their subsequent association. The comparison of the latter’s sizes in different conformational states and in a set of model two-component phosphatidylserine (PS) containing lipid membranes revealed a correlation with fluorescent measurements of CT-induced leakage of PS-containing vesicles.43 As a result, a two-stage mechanism of CT−membrane binding was hypothesized. It assumes that interactions of the toxins with cell membranes are regulated by complementarity of surface hydrophobic/hydrophilic organization of both partners.

■

METHODS MD Simulations. MD simulations in water solution were carried out for CT A3 from Naja atra, CT 1, and CT 2 from Naja oxiana. The initial models were taken from the Protein Data Bank (PDB): X-ray structure 1H0J44 and NMR-models 1RL545 and 1CB9,46 respectively. The overall quality score (from 0 to 10) of the spatial models of CTs was evaluated by means of the structure validation web service: www.prosess.ca. For each of CTs, two independent trajectories with the length of 1 and 6 μs were obtained. The all-atom Amber99SBILDN force field42 and the TIP3P water model47 were used. Brief discussion related to selection of the force field is given in the Supporting Information (SI). MD simulations were performed using the GROMACS 4.5.2 suite.48 The systems were set up and relaxed with the same computational protocol. Protein molecules were solvated with a minimum box buffer of 12 Å. 3D periodic boundary conditions were imposed. To keep the system electrically neutral, Cl− counterions were added. For better equilibration of the protein−water system, a stepwise energy minimization was carried out by the steepest descent method, fixing, at first, coordinates of all protein atoms, then backbone atoms, and, finally, Cα atoms only. Further, energy of the entire system was minimized without restrains. Then, the system was gradually heated from 5 K up to 300 K during 200 ps. Resulting configuration of the system was taken as a starting one for MD production runs at 300 K in an NPT (constant pressure and temperature) ensemble with a 2 fs time step. A double 10/12 Å cutoff and particle mesh Ewald (PME) algorithm49 were used to treat van der Waals and electrostatic interactions, respectively. During MD production run, temperature and pressure were maintained using V-rescale thermostat50 and Berendsen barostat,51 respectively. MD Data Analysis. MD snapshots for subsequent analysis were extracted from the production runs at intervals of 100 ps. Standard MD data processing including RMSD, secondary structure, and backbone dihedral angle calculations were made using utilities supplied with the GROMACS package. Structural variations of the whole protein molecule as well as its individual domains (tips of loop I residues 4−12; loop II 27−35; loop III 40−49, and loop IV 15−19) were characterized through the changes in backbone RMSDs. 3D visualization of molecular structures was made using PyMOL52 tools. Analysis of φ and ψ Dihedral Transitions. To detect and characterize local conformational changes in the toxin structures, φ and ψ values for all residues in each MD frame were calculated. The resulting data on dihedrals distribution were converted into normalized histograms with a set of bins of 1° width each. Hereinafter, we use the term “angle” for the backbone φ and ψ dihedrals. The probability distribution function (P) of each angle was approximated as a sum of

Figure 1. Examples of different types of backbone dihedral angle transitions detected in the course of 6-μs MD simulations of CTs in water. (A) Normalized histograms of dihedral angle distributions. The type of backbone dihedral (ϕ or ψ) and the corresponding residue number are given above the plots. “b1”-Types of residues have two or more separate bands for at least one of their backbone dihedrals. (B) Time dependence of φ-6 and ψ-4 values for the loop I’ residues being of b1-type for all three CTs and for CT 1 only, respectively (see text for details). Data are given for every CT molecule in the order: CT 1 (top), CT 2 (middle), and CT A3 (bottom). The local functionally important structural changes of the b1-residue pair K5/L6 are directly coupled to large-scale transitions of backbone dihedrals (ψ-5, φ-6) marked as “state 1”, “state 2”, and “state 0” (state close to the starting one) on the upper panel. The thin plot corresponds to the angle ψ-4. See the text for details.

26) was detected. The residues that have single Gaussian-like distributions for both of their backbone dihedrals were considered to be rigid and named as “core-type”. The residue type was estimated taking into account one of its backbone dihedrals that has a higher dispersion. For diffuse indistinguishable peaks, the results of Gaussian fitting procedure strongly depend on the completeness degree of the available protein conformational space. Precise analysis of these multiple states was excluded from the present work. The corresponding residues were assigned to the flexible ones and named as “bending-type”. It should be noted that the focus was made on the large-scale angle transitions that are easily 2801


Article


Sequence Variability Analysis. Amino acid sequences of snake venom CTs were taken from the UniProt database (http://www.uniprot.org). A nonredundant data set of 83 mature sequences of CTs from different species was used to build multiple sequence alignment with the help of the ClustalW software.58 Variability of sequence positions was calculated using the PVS server (http://imed.med.ucm.es/ PVS/). A consensus sequence was then generated based on Shannon entropy (H) function calculations,59 with the cutoff for highly conserved positions (H ≤ 1.0).

delineated by the appearance of an individual peak separated from the starting angle population. So that, the overlap region between the two sets of data is fully absent (P = 0) or minimal (P ∼ 0) (Figures 1A, S1, and S2). In the latter case, the distance between their mean values (Pmax) exceeds more than 3 times the halfwidth of the fitted Gaussian (σ) (for at least one peak) (Figure 1A: ψ-4 CT 1). All residues with a φ/ψ angle revealing such large-scale transitions during MD simulations were designated as “b1-type” residues. As a result, we obtained three principal types of dihedrals distribution: (i) a single band well fitted with one Gaussian; (ii) superposition of several differently populated overlapping states; (iii) separated or partly overlapping bands, which are easily decomposed into Gaussians. Corresponding residues are attributed to core-, bending-, and b1-types, respectively. The last two types are defined as flexible residues. In general, the dihedrals distributions for each protein were rather close in both 1- and 6-μs MD runsthe large majority of their residues belonged to the same type. However, the main attention was given to the longer MD simulations that provided more representative sampling of the dihedral angles space. In these simulations, significant angle changes were observed more often and occurred in a number of residues after the time points, where shorter MD simulations had finished. More details are presented in the SI. Analysis and visualIzation of Hydrophobic Clusters on Molecular Surfaces of CTs and Two-Component Model Lipid Bilayers. Identification of hydrophobic clusters formed on molecular surfaces of CTs and lipid membranes was accomplished using the molecular hydrophobicity potential (MHP) approach.54 For each MD frame (extracted with a 20 ps step), MHP values were calculated in the solvent-accessible surface points. Time-averaged surface MHP distributions were visualized in terms of the two-dimensional projection maps. In the case of CTs, this was done using the protein surface topography (PST) method.55 For lipid bilayers, MHP values were projected on a plane perpendicular to the membrane normal, as described elsewhere (e.g., ref 56). The following one- and two-component hydrated lipid bilayers composed of dioleoylphosphatidylcholine (DOPC) and dioleoylphosphatidylserine (DOPS) were investigated: DOPC, DOPC/ DOPS15%, DOPC/DOPS20%, DOPC/DOPS30%, DOPC/ DOPS50%, DOPC/DOPS70%, DOPC/DOPS85%, DOPS (percent of DOPS molecules is indicated, see ref 57 for more details). Corresponding MD data were provided by Darya V. Pyrkova. MHP analysis was carried out for the final 40 ns equilibrated parts of MD trajectories. For each system, statistical error in evaluation of the hydrophobic cluster sizes distribution was calculated based on four consecutive 10 ns fragments of the equilibrated part of MD-trajectory. The obtained values over four data sets were averaged and standard error was defined. Delineation of hydrophobic/hydrophilic clusters on molecular surfaces was made as described earlier.57 MHP calculations, fitting of dihedral angles distribution curves with Gaussians as well as delineation of long-living tightly bound water molecules inside the protein loops were performed with the original in-house software. MHP calculations and visualization of hydrophobic/hydrophilic properties onto the solvent-accessible molecular surfaces of CTs were partially accomplished using the PLATINUM program (freely available on the web site http://model.nmr. ru/platinum).

■

RESULTS AND DISCUSSION Mapping of Amino Acid Residues of CTs Based on Backbone Dihedral Angle Analysis of MD Data. It is known that φ and ψ angles unambiguously describe a polypeptide fold. Also, backbone conformational changes have to be reflected in distributions of the corresponding angles. The localized changes in a simulated protein structure are quite difficult to trap (RMSD is small). Also, these events are rare in MD simulations. In order to collect a representative set of protein conformers for their exhaustive conformational analysis, three homologous CTs from snake venom were subjected to microsecond-scale MD simulations in water. Similarity of the found conformational changes and their similar localization in the highly homologous molecules increase the reliability of the conclusions made. To simplify the task of conformational analysis of long MD trajectories by reducing the number of variables and to focus exclusively on strong local changes in CTs’ backbone, such an analysis was carried out in the dihedral angles space instead of the Cartesian one. Conformational mobility of CTs expressed in terms of RMSD values, along with the examples of strong dihedrals transitions found in MD is illustrated in Figure S3. To describe conformational plasticity of every residue, a particular shape of its φ/ψ probability function was inspected, giving the main attention to residues with the mostly dispersed angles distributions. The focus was made only on these two dihedrals as their influence on global protein motions is the most considerable one. Large-amplitude transition of the angle ω corresponds to trans−cis isomerization of a peptide bond that is a very rare eventit was not observed in our simulations. Also, conformations of amino acid side chains were omitted from consideration. Examples of MD statistics related to side chain dihedrals distribution can be found elsewhere (see, e.g., Watanabe et al.53). First, all residues were divided into a small group of the structurally rigid ones (core-residues), where both angles remain in a single state (well fitted to one Gaussian) and a large group (bending residues, b-type) of conformationally flexible residues, where at least one (φ- or ψ-) distribution is fitted with more than one Gaussian (see Methods and Figure S1). Indeed, the curve fitting procedure clearly showed that most of the angles reveal multiple states and flexible regions of a molecule can be easily identified due to enrichment of these regions with the b-type residues. At the same time, the main attention was given to the cases with fully separated populations of the angles values, regarding them as different conformational states of such a residue (b1-type). This is because strong local perturbations of backbone dihedrals in conformational dynamics of CTs can potentially affect the membrane-binding potency (and, therefore, biological activity) of the toxins. 2802


Article


Figure 2. Structural/dynamic types of amino acid residues of CTs defined based on the results of 6-μs MD simulations. Core- and b1-type residues (see text for details) are given with black bold and blue color, respectively. The b1-residues of p-f lip and ψ-only types are boxed with continuous and dashed lines, respectively. Residues in β-structure are given above the sequences with “b” symbols. Unstructured tips of the loops I−IV are also indicated on the top. Toxins’ names are marked on the right. A consensus sequence (“cons”) given below was generated by analyzing all the aligned CT sequences extracted from the UniProt database. In this row, the presence of a one-letter residue code for a sequence position indicates conservation of around 80% (Shannon’s entropy value 1.0) residues within CT sequences pool, which have only homologous replacements, are marked “p” and “f” for polar and hydrophobic ones, respectively. Positions of variable residues that undergo several mutations are marked with “-” symbols.

MD-based mapping of amino acid residues of the three CTs made according to the proposed classification is shown in Figure 2. Generally, the members of CT family reveal significant homology of their amino acid sequences. As seen in Figure 2, more than half of the toxin residues are highly conserved. However, in each CT, the core-residues (Figure 2, black font, in bold type) with a single conformational state are randomly distributed in sequences and compose no more than 30% of all residues. The totality of MD data suggests that the majority of CT residues are attributed to the bending-type, whereas the highly conservative sequence fragment (residues 33−44) is enriched with the core-residues. This region includes β-strand 35−39 of the loop II and the sequence motif IDV (residues 39−41)earlier, it was considered as a characteristic sequence signature of the CT family.2 As expected, highly flexible b1residues occupy unstructured regions of the molecules (Figure 2). But attribution of a residue to the b1-type based on the sequence variability analysis only is hardly ever possible. As shown in Figure 2, the b1-residues were detected both in variable and in highly conservative pools. On the contrary, in the course of MD simulations, b1-residues are easily identified and their new conformational states usually represent timeisolated MD-populations (Figures 1B and S2) segregated from the starting ones. Due to this, their structural features can be characterized in detail. According to statistical analysis of large changes of backbone dihedrals in the nonredundant data set of experimentally determined pairs of protein structures, the most frequently observed are the so-called “p-flip” and “ψ-only” transitions.37 The former one involves changes in ψi and φi+1 angles of residues i and i + 1, respectively, if |Δψi| > 120° and |Δφi+1| > 120°. Similar conditions (|Δψi| + |Δφi+1| > 200°) together with the following second inequality|Δψi + Δφi+1 | < 50°were used to search for the peptide plane flips among the data set of pairs of X-ray models.38 All b1-residues in the loops II and IV of CTs (except CT 2) belong to the p-flips type and satisfy the both inequalities given above. It should be noted that the simultaneous agreement of the two inequalities indicates largescale rotation of the peptide plane between the two almost unaffected side chains of the corresponding pair of b1-residues (Figure 3B and C: CT 1, CT A3). Besides the p-flip changes, the ψ-only transitions in b1residues (|Δψi| > 120°) were also detected in MD simulations of CT 2 (K18, N19, P30), CT 1, and A3 (K44) (Figures 2 and 3C: CT 2). Owing to strong perturbation of CO, N−H vectors associated with the p-flip and ψ-only transitions, local alterations in backbone H-bond patterns were observed (Figure

Figure 3. Large-amplitude backbone angles transitions (b1-type of residues) observed in MD simulations of CTs. Conformations of loop fragments (only the backbone traces are shown) including one or two pairs of b1-residues (all heavy atoms are drawn) represent the starting states and MD-states colored in light-green and light-blue, respectively. To compare MD-conformations of the b1-residues (carbon atoms are cyan-colored) with their starting structures, the latter ones are displayed in green. Conformational states (in terms of dihedral angles values) that are common for all three CTs were detected in two residue pairs K5/L6 (A) and E(A)16/G17 (B) only. b1-Residues (28−32) that have different MD-conformations for every CT are located in the loop II (C). b1-Residues responsible for the indicated MD-states are marked in large font. Nonconserved residues are named with one-letter code in the order: CT 1/CT 2/CT A3. Hydrogen bonds in starting- and MD-structures are shown with black and blue dashed lines, respectively. The residues involved in hydrogen bonding are named. Tightly bound waters that appeared in the state 2 of loop I (A) of all three CTs are shown in ball representation.

3). Also, in case of large perturbation of only one backbone angle in a pair of adjacent residues, the more pronounced in this pair are the backbone distortion and side chains reorientation. Moderate-size ( 120°) of only φ-L6 (state 1) or ψ-K5 2803


Article


mainly formed by apolar tips of the loops I−III represents a principal membrane-binding motif of CTs (Figure 5A).

(state 2) angle (Figure 3A). So, taking into account the centralpoints values (Pmax) of the angles populations, no one or only the first inequality (|Δψi| + |Δφi+1| > 200°)38 was satisfied upon the transitions “state 0 → state 1” and “state 0 → state 2”, respectively. In the whole MD sets, the value of |Δψi + Δφi+1| varies and thus shows more or less degree of compensation of the side chains rotation. More noticeable side chain reorientation of L6 in the state 2 compared to the state 1 is seen in Figure 3A. Summarizing, five high mobility sites of backbone dihedrals (residues: 5−6, 28−32, 44, and 16−19 in the loops I−IV, respectively) were identified in CT molecules (Figures 2 and 4). From them, only two pairs of residues (conservative K5/L6

Figure 5. MD-derived conformations of CTs with uniform (A) and broken (B) “hydrophobic bottom” connecting the tips of the loops I− III. Hydrophobic and hydrophilic solvent-accessible surfaces of CT molecules in the conformational states 0 (A) and 2 (B) are colored green and blue according to the molecular hydrophobicity potential (MHP) values. Hydrophobic patterns formed by the tips of the loops I−III in different states are shown with dashed ellipses. Location of the residue L6, which is responsible for disintegration of the hydrophobic bottom, is indicated.

The pair of b1-residues K5/L6 is located exactly at the interface of the two hydrophobic molecular surface fragments belonging to the loops I and II (Figures 4 and 5A). In the course of MD simulation of CTs, these residues were shown to determine three conformational states of the loop I (residues: 4−12) (Figure 3A). Main characteristics of these states are given in Table 1. All three states appeared to be stable and long-lived (>1 μs) at least in one of CTs (Figure 1B, Table 1). Depending on the toxin, these states were populated quite differently. Any of the states 1 and 2 has definite H-bonding patterns that are common for all CTs. Thus, in the state 1, two H-bonds were observed: 5K:OHN:7 V and 6L:NHO:36R, but they were absent in the state 0 (Table 1, Figure 3A). Interestingly, the latter H-bond connects the loops I and II. In state 2, complete breaking (CT 1, 2) or pronounced destabilization (CT A3) of the conservative H-bond 5K:NHO:10A/F/F was observed (Figure 3A). This leads to shortening of the loop I β-sheet (residues: 2−4 and 11−13) from 3 to 2 residue pairs in the more or less significant part of the state 2 (Table 1). The same pair of terminal β-structured residues (4, 11) of CT 1 (from Indian cobra) was found to lose their H-bond during MD simulations of the toxin/membrane complex.16 Simultaneously, tightly bound water molecules forming 2−3 hydrogen bonds with amide protons (K5−V7) and backbone oxygens of A/F/F10 appeared inside the loop (Figure 3A). Compared to other toxins, the lifetime of state 2 in CT 1 was the greatest one (∼1 μs). Note that only in CT 1 the preceding residue (N4) was attributed to the b1-type with two rather distant populations of the ψ-angles (Figure 1A). In CTs, the pair of spatially adjacent residues 4 and 11 borders the unstructured hydrophobic extremity of the loop I. New stable (nearly microsecond) backbone conformations of these residues were observed in MD simulation of CT 1 only after transitions of K5/L6 dihedrals into the state 2 (Figure 1B). Thus, it seems reasonable to assume that the MD-observed long-living b1-residue conformational states promote dihedral mobility of the neighboring residues and make conformational changes in this region more pronounced (Table 1). Nevertheless, irreversible alterations in the loop structure of CT 1 do

Figure 4. Structural/dynamic types of CT residues based on the mobility of their backbone dihedrals in the course of long-term MD simulations. Spatial models of CT 1 (PDB code: 1RL5), CT 2 (1CB9), and A3 (1H0J) superimposed over backbone atoms of the βstructured residues are shown in ribbon representation. Core- and b1type residues (see text for details) are designated with their side chain atoms in line and stick representation, respectively. The color scheme corresponds to three shades (dark (CT A1), medium CT 2), light (CT A3)) of gray and blue for ribbon structure and side chain atoms of b1residues, respectively. Flexible residues attributed to the b1-type in all three CTs and located in the same sequence positions are named with one-letter codes (nonconserved residues are given in the order: CT 1/ CT 2/CT A3).

in the loop I and E(A)16/G17 in the loop IV) were attributed to the b1-type with the same transition states (corresponding Pmax values are close) and were observed in MD simulations of all three CTs (Figures 3A, 3B, S4). However, loop IV is located on the opposite side from the membrane-active domain of CTs, and assumptions about its role in membrane binding seem unlikely. On the contrary, the remaining pair of conservative b1-residues (K5/L6) probably has a functional meaning. These residues lie at the interface of the structured domains of the loops I and II. Hence, local conformational rearrangements associated with changes of these dihedrals can modify the functionally active pattern on molecular surface of CTs. No correlations between the conformational changes provoked by b1-residues pairs of the loops I and IV were found. In the next section, the site of structural changes formed by the key pair of b1-residues K5/L6 and its possible role in modulation of membrane-binding activity of CTs are analyzed in detail. Conformational Plasticity of the “Hydrophobic Bottom” of CTs. Experiments on model membranes14,15 have demonstrated that the extended hydrophobic bottom 2804


Article


Table 1. Parameters of MD-Observed Conformational States of Loop I in CTs Promoted by Large-Amplitude Dihedral Angle Transitions of Residues K5/L6 toxin

Pmax ψ-K5a

Pmax φ-L6a

CT A3 CT 2 CT 1

179 ± 8 (179)b 178 ± 8 (−170) 178 ± 10 (−167)

−56 ± 10 (−42) −61 ± 9 (−50) −64 ± 11 (−43)

CT A3 CT 2 CT 1

123 ± 13 116 ± 12 119 ± 11

60 ± 7 60 ± 7 60 ± 7

CT A3 CT 2 CT 1

−39 ± 14 −32 ± 14 −26 ± 15

−140 ± 16 −139 ± 16 −140 ± 11

occupancy (%)c

β-sheet (2−4,11−13) of loop I (%)d

86 87 4 11 1.5 36 1 9 49

State 0 99 99 97 State 1 96 97 99 State 2 80 36 15

lifetimee

H-bonds (%)d 5K:NH−O:10A/F/F; 5K:O− HN:7 V; 6L:NH−O:36R

>1 μs >1 μs 0.1 μs 1 μs

99; 86; 60 92; 87; 89 98; 91; 72

>0.01 μs >0.5 μs >1 μs

58; 0; 0 0; 0; 0 1; 0; 0

Mean ± SD are shown. bCorresponding dihedral angles values in the experimental 3D structures of CTs used as starting models in MD simulations. cOccupancy (%) of a given conformational state in the corresponding MD-trajectory. dNumber of representatives of a given state, which are characterized by the indicated structural parameter. eMaximal values are shown.

a

Figure 6. Structural rearrangements of CT 1, CT 2, and A3 observed in the course of MD simulations. Results of the PST calculations.55 MDaveraged molecular hydrophobicity potential (MHP) spherical projection maps for the MD-conformational states 0, 1, 2. MHP values on the protein solvent-accessible surface are color-coded according to the scale on the right. Projections of centers of mass of side chain atoms are given with the corresponding residue names and numbers. Only the main hydrophobic residues forming the membrane-binding motif (“bottom”) of CTs are marked. The site of hydrophobic linkage combining the loop I and loop II−III regions into a single hydrophobic bottom is marked with a white arrow.

The second flexible domain in the hydrophobic bottom of CTs includes tips of the loops II and III (Figure 4). In the most variable loop II (residues: 27−35), distinct structural transitions were observed for each of CTs. Almost all residue positions in this loop are assigned with nonhomologous mutations (Figure 2) including b1-residues site in the center of the loop (Figure 4). In all MD simulations, loop III of CTs possessed the lowest fluctuations near the starting structure (data not shown). Unstructured region of the loop III (residues: 40−48) contains mainly conserved residues but only three of them (47−49) are donating to the hydrophobic bottom (Figure 5). Dihedral transitions in the remote residue K44 of CT 1 and A3 (Figure S2) were found to promote a reversible long-living conformational change with a high level of local deformation of the loop (data not shown). However, alternative conformations of the b1-residues of both loops induce local changes in this region without affecting integrity of the hydrophobic bottom. The conservative (among the CT family) omega-shaped structure of the loop II with H-bonded water molecule in its cavity also remains intact. The stabilizing role of such waters found in a

not occur, and the loop I returns back nearly to its starting conformation. At last, the most intriguing finding is related to rearrangements of molecular surfaces of CTs in the state 2. The local change of the backbone geometry together with reorientation of the side chain of L6 in the protein segment 4−7 grows in the row of states: 0 < 1 < 2 and leads to interruption of the continuous hydrophobic bottom connecting the tips of the loops I−III (Figures 5, 6). In Figure 6, spherical projections of the entire toxins’ molecular surfaces with mapped hydrophobic/hydrophilic properties that are calculated and averaged over MD trajectories are presented. These projections were generated using the “protein surface topography” (PST) method.55 PST permits the mapping physicochemical properties (e.g., hydrophobicity) on a molecular surface and visualize them by projection onto a sphere. It is seen that in all CTs, a hydrophilic “groove” splitting two apolar patterns of the loops I and II starts to appear in the state 1 and becomes very prominent in the state 2. 2805


Article


their surface polarity patterns match those appearing on the membrane surface well. To characterize the latter, MHP analysis for a set of one- and two-component lipid bilayers was carried out (see Methods). The main goal was to evaluate the characteristic dimensions of dynamic nonpolar domains (clusters) on the membrane surface and compare them with the amphiphilic surface patterns observed for CTs (see above). Analysis of MD-data obtained for zwitterionic (DOPC) membranes with varying amounts (from 0 to 100%) of anionic (DOPS) lipids leads to the following conclusions (Figure 7A):

number of NMR and X-ray models of CTs was previously discussed in Konshina et al.1 Summarizing, the membrane-binding site of CTs is formed by three spatially distant segmentsunstructured tips of the loops I−III, each of which includes residues that undergo largeamplitude dihedrals transitions in the course of MD simulations. The local long-living structural changes of b1residues, which are inherent in all CTs, were delineated only in their loops I. It is important that such events were stably reproduced in all MD runs. One of the most interesting results obtained in this work is the detected reorganization of the surface hydrophobic pattern in the membrane-binding site of CTsthis rearrangement is directly coupled to large transitions of the backbone dihedrals within the residue pair K5/L6. The possible functional impact of this phenomenon is discussed below. Complementarity of Hydrophobic Patterns on Molecular Surfaces of CTs and Model Lipid Bilayers. Being frequently observed in cytolytic peptides,60 surface amphiphilicity is considered as one of the major players determining their membrane activity. That is why, delineation of the factors regulating amphiphilicity “portrait” of a protein is required for elucidation of its intimate structure/activity relationships. All experimentally determined spatial structures of CT 1, 2, and A3 manifest on their molecular surfaces an extended continuous hydrophobic patch (or bottom) formed by the tips of loops I− III and framed with a belt of cationic residues. At the same time, such a pattern is not staticin the course of MD simulations it undergoes significant rearrangements determined by conformational dynamics of CTs. The most prominent changes are related to breaking of the hydrophobicity zone into two nonpolar fragments formed by the loop I and the loops II− III without a hydrophobic “bridge” between them. Such transformations are caused by the high-amplitude backbone dihedrals transitions in the pair of residues K5, L6 and subsequent reorientation of the side chain of L6 (Figures 5 and 6). On average, the aforementioned two apolar patches observed in MD states of CTs with the “broken bottom” (state 2, see Table 1) have areas of their two-dimensional projections c.a. 100 and 200 Å2, respectively (Table 2). In MD states with the continuous apolar bottom formed by all three loops (state 0, Table 2), the corresponding average area is about 270−350 Å2.

Figure 7. Size distributions of hydrophobic clusters on molecular surfaces of model lipid bilayers. (A) Results of MD simulations of the fully hydrated one- and two-component model DOPC/DOPS membranes.57 Bilayer compositions are shown on the right. The cluster size corresponds to the area of a hydrophobic pattern obtained via projection of the membrane surface properties onto a plane. (B) Dependence of CT-induced efficiency of leakage of fluorescent dye on the PS content (%) in POPC liposomes (adapted from ref 43). Vertical bars indicate standard error.

Table 2. Areasa of the Major Hydrophobic Pattern (Bottom) on Molecular Surfaces of CTs

(i) Hydrophobic patterns with an area up to 100 Å2 appeared in all water−lipid systems with high probability. (ii) A number of large (with area >200 Å2) hydrophobic clusters was delineated only in lipid bilayers containing at least 20% of DOPS. Such nanometer-scale lateral heterogeneities in PC/PS mixtures have been described earlier.57 Also, it has been shown that the hydrophobic clusters can efficiently promote membrane binding of external agents, like membrane active peptides.65 In this case, protein−membrane interaction is driven by the self-accommodation of amphiphilic peptides on the heterogeneous hydrophobic/hydrophilic membrane surface. So, under certain conditions, the hydrophobic cluster sizes of the two interacting partnersthe CT and lipid bilayer perfectly fit each other. Based on the observed correlations, we propose the following two-stage mechanism of CT−membrane binding. In the beginning, the CT molecule embeds through the tip of its loop I. The modest hydrophobicity of this loop

two patterns

CT A3 CT 2 CT 1

single pattern: loops I−III hydrophobic bottom, Å2

loop I, Å2

loops II− III, Å2

311 ± 22 346 ± 23 266 ± 24

95 ± 23 101 ± 22 98 ± 9

213 ± 32 220 ± 40 161 ± 25

a

MD-averaged value (±SD) of the area calculated for planar projection of the pattern.

Probably, such a dynamic amphiphilicity of CTs surface is important for the binding of these proteins to lipid bilayers surfaces of the latter also have a fluctuating “mosaic” nature.56,61,62 Because intermolecular interactions often require a complementarity of surface hydrophobic−hydrophilic patterns of the partners,63,64 it was hypothesized that efficient membrane binding of CTs is accomplished only in cases when 2806


Article


sequence were attributed to distinct structural/dynamic types (core or bending) in different CTs (Figure 2). Based on MDderived conformational mobility data for the selected molecular regions, noticeable differences in the flexibility of loop I (decreasing in the row: CT 1 ≥ CT 2 > CT A3) were delineated (Figure S5). Of note, the conformational transitions associated with the changes of the backbone dihedrals in residues K5/L6 are less pronounced in MD simulations of CT A3. So, the lifetimes of the MD-state 2 in this toxin (∼10 ns) were much shorter than those observed for CT 1 (1 μs; Figure 1B, Table 1). It is remarkable that such a tendency of the loop I conformational mobility to decrease correlates with increasing membrane-damaging activity of the given CTs with respect to liposomes.43 It also agrees with our earlier modeling results concerning the relative efficiency of the membrane binding of CT 1 and CT 2.67,68 However, when postulating that the strength of membrane−CT interaction correlates with the degree of loop I conformational mobility, one should keep in mind the following potential shortcomings of the computational approach: (1) Minor conformational differences in the starting protein models may have major effects on their backbone dynamics in the course of MD simulations. This is a well-known problem of the reproducibility of MD-based structural parameters compared with the experimentally determined ones, in particular, with NMR-derived order parameters.69 Obviously, the larger the system, the more computational resources are required for its relaxation and the less adequate the conformational sampling possible. It has been mentioned that the high-resolution X-ray structures taken as the starting point usually appeared less flexible in MD calculations than NMRmodels.70−72 Indeed, compared with both starting NMRmodelsCT 1 and CT 2the X-ray structure of CT A3 demonstrated a more stable behavior, especially in the loop I fragment (Figure S5). Interestingly, the amino acid sequences of the loop I regions in CT 1 and A3 are practically identical (Figure 2). Instead, their mobility scales were found to be significantly different. Therefore, it can be explained by the individual structural properties of the toxins and also depends on the quality of experimentally determined starting structures. Note, the X-ray model of CT A3 has the higher index of overall quality among the studied CTs (3.5, 3.5, and 5.5 for CT 1, 2, and A3, respectively). (2) Comparing the membrane-active sites of CTs and estimating their influence on the strength of membrane binding, one should keep in mind the major effect of the loop II ending. Being the most variable in the sequence, it usually contains positively charged residues with extended side chains (K31 in CT A3), whose length facilitates the adjustment of the hydrophobic membrane binding site despite the charged locus. However, CT 1 has D29 here and shows the lowest lytic activity among the CTs under consideration (Figure 7B).43 The spatial structure of this toxin was solved in both aqueous (PDB code: 1RL5) and membranous (1ZAD, in DPC micelles) solutions. As seen from the comparison of the structure of loop II in these two spatial models, the adaptation of the corresponding structural fragment to the micelle environment was observed.74 Some variations of the backbone dihedrals, including large-amplitude angle

and the apolar clusters frequently observed on the surface of all the considered lipid bilayers have a similar size. At the next stage, the toxin binds to the membrane with a bigger hydrophobic fragment (loops II−III). The suggested hypothesis was validated by comparing MD results with the experimental data on the lytic efficiency of given CTs against POPC liposomes containing different amounts (20, 35, 50, and 70%) of anionic lipid (PS) (Figure 7B).43 The membranedamaging activity of CTs was measured by the release of the liposome-entrapped calcein fluorescent dye. For all three toxins, the lytic effect was found to be the lowest in neutral POPC liposomes and highest in the mixed liposomes containing c.a. 20% of PS. We explain this well-defined maximum by appearance of a large number of lipid hydrophobic clusters with a size similar to that of the apolar bottom of the loops II− III. Lipids in highly charged (>20% PS) model bilayers were also found to form big (>200 Å2) hydrophobic clusters on the membrane surface (Figure 7A). At higher PS contents, the lytic effect of CT 2 and A3 decreased but still remained at a substantial level. Obviously, the role of electrostatic interactions between the positively charged CT molecules and the anionic liposomes has to be considerable. Energetically favorable contacts of protein groups with lipid heads may have an ambivalent effect on the CT membrane activity: (i) initial attraction and retention of toxin molecules by the lipid surface; (ii) subsequent prevention of deep insertion of CT into the membrane. Due to the formation of a number of ionic contacts, especially in case of highly charged liposomes, the toxin molecule “sticks” to the bilayer surface and may adopt an orientation unfavorable to the effective recognition of the hydrophobic zones and further embedding. To summarize, we suggest that the initial step of the CT− membrane association includes long-range electrostatic attraction of the basic toxins to the negatively charged membrane surface followed by highly complementary hydrophobic interactions between apolar clusters of a similar size promoting the interface-inserted mode of toxin binding. This initial stage of protein−membrane recognition can be strengthened by the specific short-range “trapping” of anionic lipid heads by particular toxin’s sites as discussed by Konshina et al.66 In the present study, we do not consider CT binding to membranes with atomistic detailsthis work is in progress now. At the same time, the proposed two-stage model and the importance of mutual adaptation of amphiphilicity patterns of the interacting partnersthe toxin and the lipid bilayerperfectly agree with similar mechanisms determining membrane binding of linear cationic peptides.62,65 MD-Derived Conformational Differences among CTs have no Unambiguous Structure−Activity Interpretation: Limitations of the Computational Approach. One of the intriguing challenges is that structurally rigid CT molecules exhibit only subtle conformational differences in the known experimental spatial models, although they differ greatly in their biological activities.5,10,43 However, during the discussion we deliberately did not emphasize the differences revealed in the dynamic behavior of the toxins. In particular, each of the three CTs possesses a unique ensemble of conformational states associated with the discovered potentially important dihedrals transition in the K5/L6 pair of residues (Figure 1, 3A, Table 1). Also, new states of some b1-residues found in one CT were not observed in MD simulations of the others (Figure 2). Finally, in a number of cases, residues with similar positions in the 2807


Article


here, a detailed inspection of protein dihedrals is a powerful and useful complement to standard MD analysis.

transitions in the key charged residue, were found in MD simulations (Figure 3C). But, strong conformational changes that we can see in the corresponding NMRmodel of CT 1 (1ZAD) are likely to require a membranous environment and have not been revealed in long-term MD simulations in water. To summarize, all the differences in the structural/dynamic behavior of CTs detected in MD simulations are interesting and may have a biological meaning. On the other hand, some pitfalls such as the insufficient conformational sampling and the influence of the starting structure remain uncertain and require careful case-by-case consideration. Taking into account all these potential shortcomings, in this study we concentrated on the search and interpretation of only the most prominent common trends in the conformational dynamics of CT molecules. Molecular simulations are useful in exploring new conformational states that deviate more or less from the initial structure. However, the experimental validation of conformations is often difficult. Standard NMR experiments only provide an averaged molecular structure near its ground state and do not detect the so-called “minor impurities” of low-populated alternative structures. Modern complex methodologies, in particular relaxation dispersion NMR spectroscopy,73 make it possible to conduct detailed structural studies of such poorly presented conformational states of a protein. Although the experimental work is laborious and requires 15N−13C labeled toxin samples, it provides an opportunity for the experimental confirmation of MD-elucidated toxin conformations in the future.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications Web site. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jcim.7b00395. Selection of the force fields: Brief discussion concerning the selected force field advantages. Figure S1: backbone dihedral analysis scheme underlying the proposed residue classification. Figure S2: Examples of different types of backbone dihedral angle transitions detected in the course of 6-μs MD simulations of CTs in water. Figure S3: Time dependence of backbone RMSDs and some backbone dihedrals, which undergo large-amplitude transitions during simulations. Reproducibility of structural/dynamic characteristics of CTs in MD simulations: Details of common and different features of short (1 μs) and long (6 μs) MD simulations of each cardiotoxin. MD-derived conformational mobility of CTs in terms of their backbone dihedrals distributions (Figure S4) and root-mean-square fluctuations of backbone coordinates (Figure S5) (PDF)

■

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected].

CONCLUSIONS The mapping of the conformational mobility of CTs in terms of backbone dihedrals φ/ψ transitions for every residue allows to delineate a number of new conformational states alternative to the experimentally determined ones. These states were well reproduced in MD simulations of all three CTs. Large-scale dihedral transitions and associated local conformational changes were observed within the well-known functionally active region, which has been well-established based on the analysis of experimental models of the toxins. It has never been shown before that the local transformations of only a pair of residues can play a crucial role in membrane binding. Thus, reversible large-scale transitions of backbone dihedrals in residues K5 and L6 result in corresponding breaking/ association of the large hydrophobic bottom on CTs surface. This, in turn, regulates the toxins interactions with polar/ nonpolar zones of the membranes. Hence, compatible with experiment, our results meaningfully complement the existing picture. The identification of the large-amplitude backbone angle transitions in the course of long-term MD simulations may be useful for the detection of peculiar hot spots in protein conformational changes which are functionally relevant but extremely local. On one side, these structural features may be missing in usual clustering procedures based on the RMSD differences of numerous MD-models. This is mainly due to structurally rigid proteins possessing different activities. On the other side, the prediction of functionally active residues based on sequence alignment and the search for conserved residues in potentially important regions are nonobvious and quite difficult, especially in the case of highly homologous proteins, where many residues are conserved. Hence, according to these criteria, CTs represent a class of rather complex molecules to address the aforementioned issues. In this situation, as shown

ORCID

Roman G. Efremov: 0000-0002-5474-4721 Author Contributions

A.G.K. performed data analysis and wrote the paper. N.A.K. developed original analysis algorithms. R.G.E. supervised the research, performed data analysis, and wrote the paper. Funding

This work was supported by the Russian Science Foundation (grant 14-50-00131). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Access to computational facilities of the Supercomputer Center “Polytechnical” of St. Petersburg Polytechnic University and the Joint Supercomputer Center of RAS (Moscow) is greatly appreciated. We thank Ms. Anastasia Efremova for editing the manuscript.

■

ABBREVIATIONS CTs, cardiotoxins; DOPC, dioleoylphosphatidylcholine; DOPS, dioleoylphosphatidylserine; MD, molecular dynamics; MHP, molecular hydrophobicity potential; RMSD, root mean-squared deviation; PCA, principal component analysis; PDB, Protein Data Bank; POPC, palmitoyloleoylphosphatidylcholine; PC, phosphatidylcholine; PS, phosphatidylserine; 3D, three-dimensional

■

REFERENCES

(1) Konshina, A. G.; Dubovskii, P. V.; Efremov, R. G. Structure and Dynamics of Cardiotoxins. Curr. Protein Pept. Sci. 2012, 13, 570−584.

2808


Article

Journal of Chemical Information and Modeling (2) Kumar, T. K. S.; Jayaraman, G.; Lee, C. S.; Arunkumar, A. I.; Sivaraman, T.; Samuel, D.; Yu, C. Snake Venom Cardiotoxins Structure, Dynamics, Function and Folding. J. Biomol. Struct. Dyn. 1997, 15, 431−463. (3) Kini, R. M.; Doley, R. Structure, Function and Evolution of Three-Finger Toxins: Mini Proteins with Multiple Targets. Toxicon 2010, 56, 855−867. (4) Galat, A.; Gross, G.; Drevet, P.; Sato, A.; Menez, A. Conserved Structural Determinants in Three-Fingered Protein Domains. FEBS J. 2008, 275, 3207−3225. (5) Dufton, M. J.; Hider, R. C. Structure and Pharmacology of Elapid Cytotoxins. Pharmacol. Ther. 1988, 36, 1−40. (6) Harvey, A. L.; Marshall, R. J.; Karlsson, E. Effects of Purified Cardiotoxins from the Thailand Cobra (Naja Naja Siamensis) on Isolated Skeletal and Cardiac Muscle Preparations. Toxicon 1982, 20, 379−396. (7) Hodges, S. J.; Agbaji, A. S.; Harvey, A. L.; Hider, R. C. Cobra Cardiotoxins. Purification, Effects on Skeletal Muscle and Structure/ Activity Relationships [published errtum appears in Eur J Biochem 1988 Feb 1;171(3):727]. Eur. J. Biochem. 1987, 165, 373−383. (8) Ownby, C. L.; Fletcher, J. E.; Colberg, T. R. Cardiotoxin-1 from Cobra (Naja Naja Atra) Venom Causes Necrosis of Skeletal Muscle in vivo. Toxicon 1993, 31, 697−709. (9) Stevens-Truss, R.; Hinman, C. L. Activities of Cobra Venom Cytotoxins toward Heart and Leukemic T-Cells Depend on Localized Amino Acid Differences. Toxicon 1997, 35, 659−669. (10) Feofanov, A. V.; Sharonov, G. V.; Dubinnyi, M. A.; Astapova, M. V.; Kudelina, I. A.; Dubovskii, P. V.; Rodionov, D. I.; Utkin, Y. N.; Arseniev, A. S. Comparative Study of Structure and Activity of Cytotoxins from Venom of the Cobras Naja Oxiana, Naja Kaouthia, and Naja Haje. Biochemistry (Moscow) 2004, 69, 1148−1157. (11) Wang, C. H.; Monette, R.; Lee, S. C.; Morley, P.; Wu, W. G. Cobra Cardiotoxin-Induced Cell Death in Fetal Rat Cardiomyocytes and Cortical Neurons: Different Pathway but Similar Cell Surface Target. Toxicon 2005, 46, 430−440. (12) Chen, L. W.; Kao, P. H.; Fu, Y. S.; Lin, S. R.; Chang, L. S. Membrane-Damaging Activity of Taiwan Cobra Cardiotoxin 3 is Responsible for its Bactericidal Activity. Toxicon 2011, 58, 46−53. (13) Dubovskii, P. V.; Konshina, A. G.; Efremov, R. G. Cobra Cardiotoxins: Membrane Interactions and Pharmacological Potential. Curr. Med. Chem. 2013, 21, 270−287. (14) Dubovskii, P. V.; Dementieva, D. V.; Bocharov, E. V.; Utkin, Y. N.; Arseniev, A. S. Membrane Binding Motif of the P-Type Cardiotoxin. J. Mol. Biol. 2001, 305, 137−149. (15) Dauplais, M.; Neumann, J. M.; Pinkasfeld, S.; Menez, A.; Roumestand, C. An NMR Study of the Interaction of Cardiotoxin Gamma from Naja Nigricollis with Perdeuterated Dodecylphosphocholine Micelles. Eur. J. Biochem. 1995, 230, 213−220. (16) Gorai, B.; Sivaraman, T. Delineating Residues for Haemolytic Activities of Snake Venom Cardiotoxin 1 from Naja Naja as Probed by Molecular Dynamics Simulations and in vitro Validations. Int. J. Biol. Macromol. 2017, 95, 1022−1036. (17) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H. F.; Shaw, D. E. Biomolecular Simulation: A Computational Microscope for Molecular Biology. Annu. Rev. Biophys. 2012, 41, 429−452. (18) Lazaridis, T. Implicit Solvent Simulations of Peptide Interactions with Anionic Lipid Membranes. Proteins: Struct., Funct., Genet. 2005, 58, 518−527. (19) Su, Z. Y.; Wang, Y. T. Coarse-Grained Molecular Dynamics Simulations of Cobra Cytotoxin A3 Interactions with a Lipid Bilayer: Penetration of Loops into Membranes. J. Phys. Chem. B 2011, 115, 796−802. (20) Efremov, R. G.; Volynsky, P. E.; Polyansky, A. A.; Nolde, D. E.; Arseniev, A. S. Protein-Membrane Interactions: Lessons from in silico Studies. In Structure and Biophysics−New Technologies for Current Challenges in Biology Beyond; Puglisi, J. D., Ed.; Springer Netherlands: Dordrecht, 2007; Vol. 231, Chapter 3, pp 19−39. (21) Sue, S. C.; Jarrell, H. C.; Brisson, J. R.; Wu, W. G. Dynamic Characterization of the Water Binding Loop in the P-Type

Cardiotoxin: Implication for the Role of the Bound Water Molecule. Biochemistry 2001, 40, 12782−12794. (22) Gorai, B.; Prabhavadhni, A.; Sivaraman, T. Unfolding Stabilities of Two Structurally Similar Proteins as Probed by TemperatureInduced and Force-Induced Molecular Dynamics Simulations. J. Biomol. Struct. Dyn. 2015, 33, 2037−2047. (23) Dubovskii, P. V.; Vassilevski, A. A.; Kozlov, S. A.; Feofanov, A. V.; Grishin, E. V.; Efremov, R. G. Latarcins: Versatile Spider Venom Peptides. Cell. Mol. Life Sci. 2015, 72, 4501−4522. (24) Fenwick, R. B.; Esteban-Martin, S.; Salvatella, X. Understanding Biomolecular Motion, Recognition, and Allostery by Use of Conformational Ensembles. Eur. Biophys. J. 2011, 40, 1339−1355. (25) Spyrakis, F.; BidonChanal, A.; Barril, X.; Luque, F. J. Protein Flexibility and Ligand Recognition: Challenges for Molecular Modeling. Curr. Top. Med. Chem. 2011, 11, 192−210. (26) Boehr, D. D.; Nussinov, R.; Wright, P. E. The Role of Dynamic Conformational Ensembles in Biomolecular Recognition. Nat. Chem. Biol. 2009, 5, 789−796. (27) Chong, L. T.; Saglam, A. S.; Zuckerman, D. M. Path-Sampling Strategies for Simulating Rare Events in Biomolecular Systems. Curr. Opin. Struct. Biol. 2017, 43, 88−94. (28) Maximova, T.; Moffatt, R.; Ma, B. Y.; Nussinov, R.; Shehu, A. Principles and Overview of Sampling Methods for Modeling Macromolecular Structure and Dynamics. PLoS Comput. Biol. 2016, 12, e1004619. (29) Shao, J. Y.; Tanner, S. W.; Thompson, N.; Cheatham, T. E. Clustering Molecular Dynamics Trajectories: 1. Characterizing the Performance of Different Clustering Algorithms. J. Chem. Theory Comput. 2007, 3, 2312−2334. (30) David, C. C.; Jacobs, D. J. Principal Component Analysis: A Method for Determining the Essential Dynamics of Proteins. Methods Mol. Biol. 2014, 1084, 193−226. (31) Amadei, A.; Linssen, A. B. M.; Berendsen, H. J. C. Essential Dynamics of Proteins. Proteins: Struct., Funct., Genet. 1993, 17, 412− 425. (32) Keller, B.; Daura, X.; van Gunsteren, W. F. Comparing Geometric and Kinetic Cluster Algorithms for Molecular Simulation Data. J. Chem. Phys. 2010, 132, 074110. (33) Frickenhaus, S.; Kannan, S.; Zacharias, M. Efficient Evaluation of Sampling Quality of Molecular Dynamics Simulations by Clustering of Dihedral Torsion Angles and Sammon Mapping. J. Comput. Chem. 2009, 30, 479−492. (34) Altis, A.; Nguyen, P. H.; Hegger, R.; Stock, G. Dihedral Angle Principal Component Analysis of Molecular Dynamics Simulations. J. Chem. Phys. 2007, 126, 244111. (35) Sittel, F.; Jain, A.; Stock, G. Principal Component Analysis of Molecular Dynamics: on the Use of Cartesian vs. Internal Coordinates. J. Chem. Phys. 2014, 141, 014111. (36) Omori, S.; Fuchigami, S.; Ikeguchi, M.; Kidera, A. Latent Dynamics of a Protein Molecule Observed in Dihedral Angle Space. J. Chem. Phys. 2010, 132, 115103. (37) Nishima, W.; Qi, G. Y.; Hayward, S.; Kitao, A. DTA: Dihedral Transition Analysis for Characterization of the Effects of Large MainChain Dihedral Changes in Proteins. Bioinformatics 2009, 25, 628− 635. (38) Hayward, S. Peptide-Plane Flipping in Proteins. Protein Sci. 2001, 10, 2219−2227. (39) Yang, M. F.; Lei, M.; Yordanov, B.; Huo, S. H. Peptide Plane Can Flip in Two Opposite Directions: Implication in Amyloid Formation of Transthyretin. J. Phys. Chem. B 2006, 110, 5829−5833. (40) Armen, R. S.; DeMarco, M. L.; Alonso, D. O. V.; Daggett, V. Pauling and Corey’s Alpha-Pleated Sheet Structure May Define the Prefibrillar Amyloidogenic Intermediate in Amyloid Disease. Proc. Natl. Acad. Sci. U. S. A. 2004, 101, 11622−11627. (41) Kitao, A.; Hayward, S.; Go, N. Energy Landscape of a Native Protein: Jumping-Among-Minima Model. Proteins: Struct., Funct., Genet. 1998, 33, 496−517. (42) Lindorff-Larsen, K.; Piana, S.; Palmo, K.; Maragakis, P.; Klepeis, J. L.; Dror, R. O.; Shaw, D. E. Improved Side-Chain Torsion Potentials 2809


Article

Journal of Chemical Information and Modeling for the Amber ff99SB Protein Force Field. Proteins: Struct., Funct., Genet. 2010, 78, 1950−1958. (43) Konshina, A. G.; Boldyrev, I. A.; Omel'kov, A. V.; Utkin, Y. N.; Efremov, R. G. Anionic Lipids: Determinants of Binding Cytotoxins from Snake Venom on the Surface of Cell Membranes. Acta Naturae 2010, 2, 88−95. (44) Forouhar, F.; Huang, W. N.; Liu, J. H.; Chien, K. Y.; Wu, W. G.; Hsiao, C. D. Structural Basis of Membrane-Induced Cardiotoxin A3 Oligomerization. J. Biol. Chem. 2003, 278, 21980−21988. (45) Dubovskii, P. V.; Lesovoy, D. M.; Dubinnyi, M. A.; Konshina, A. G.; Utkin, Y. N.; Efremov, R. G.; Arseniev, A. S. Interaction of ThreeFinger Toxins with Phospholipid Membranes: Comparison of S- and P-Type Cytotoxins. Biochem. J. 2005, 387, 807−815. (46) Dementieva, D. V.; Bocharov, E. V.; Arseniev, A. S. Two Forms of Cytotoxin II (Cardiotoxin) from Naja Naja Oxiana in Aqueous Solution - Spatial Structures with Tightly Bound Water Molecules. Eur. J. Biochem. 1999, 263, 152−162. (47) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (48) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845−854. (49) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: an Nlog(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (50) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling Through Velocity Rescaling. J. Chem. Phys. 2007, 126, 014101. (51) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684−3690. (52) The PyMOL Molecular Graphics System, Version 1.6; Schrödinger, LLC, 2013. (53) Watanabe, H.; Elstner, M.; Steinbrecher, T. Rotamer Decomposition and Protein Dynamics: Efficiently Analyzing Dihedral Populations from Molecular Dynamics. J. Comput. Chem. 2013, 34, 198−205. (54) Efremov, R. G.; Alix, A. J. Environmental Characteristics of Residues in Proteins: Three-Dimensional Molecular Hydrophobicity Potential Approach. J. Biomol. Struct. Dyn. 1993, 11, 483−507. (55) Koromyslova, A. D.; Chugunov, A. O.; Efremov, R. G. Deciphering Fine Molecular Details of Proteins’ Structure and Function with a Protein Surface Topography (PST) Method. J. Chem. Inf. Model. 2014, 54, 1189−1199. (56) Pyrkova, D. V.; Tarasova, N. K.; Pyrkov, T. V.; Krylov, N. A.; Efremov, R. G. Atomic-Scale Lateral Heterogeneity and Dynamics of Two-Component Lipid Bilayers Composed of Saturated and Unsaturated Phosphatidylcholines. Soft Matter 2011, 7, 2569−2579. (57) Pyrkova, D. V.; Tarasova, N. K.; Krylov, N. A.; Nolde, D. E.; Pentkovsky, V. M.; Efremov, R. G. Dynamic Clustering of Lipids in Hydrated Two-Component Membranes: Results of Computer Modeling and Putative Biological Impact. J. Biomol. Struct. Dyn. 2013, 31, 87−95. (58) Larkin, M.; Blackshields, G.; Brown, N.; Chenna, R.; McGettigan, G.; McWilliam, H.; Valentin, F.; Wallace, I.; Wilm, A.; Lopez, R.; Thompson, J.; Gibson, T.; Higgins, D. ClustalW and ClustalX Version 2. Bioinformatics 2007, 23, 2947−2948. (59) Litwin, S.; Jores, R. Shannon Information as a Measure of Amino Acid Diversity. In Theoretical and experimental insights into immunology; Perelson, A.; Weisbuch, G., Eds.; Springer: Berlin, Heidelberg, 1992; pp 279−287. (60) Kini, R. M.; Evans, H. J. A Common Cytolytic Region in Myotoxins, Hemolysins, Cardiotoxins and Antibacterial Peptides. Int. J. Pept. Protein Res. 1989, 34, 277−286. (61) Polyansky, A. A.; Volynsky, P. E.; Arseniev, A. S.; Efremov, R. G. Adaptation of a Membrane-active Peptide to Heterogeneous Environ-

ment. I. Structural Plasticity of the Peptide. J. Phys. Chem. B 2009, 113, 1107−1119. (62) Polyansky, A. A.; Volynsky, P. E.; Arseniev, A. S.; Efremov, R. G. Adaptation of a Membrane-Active Peptide to Heterogeneous Environment. II. The Role of Mosaic Nature of the Membrane Surface. J. Phys. Chem. B 2009, 113, 1120−1126. (63) Efremov, R. G.; Chugunov, A. O.; Pyrkov, T. V.; Priestle, J. P.; Arseniev, A. S.; Jacoby, E. Molecular Lipophilicity in Protein Modeling and Drug Design. Curr. Med. Chem. 2007, 14, 393−415. (64) Pyrkov, T. V.; Kosinsky, Y. A.; Arseniev, A. S.; Priestle, J. P.; Jacoby, E.; Efremov, R. G. Complementarity of Hydrophobic Properties in ATP-Protein Binding: A New Criterion to Rank Docking Solutions. Proteins: Struct., Funct., Genet. 2007, 66, 388−398. (65) Polyansky, A. A.; Ramaswamy, R.; Volynsky, P. E.; Sbalzarini, I. F.; Marrink, S. J.; Efremov, R. G. Antimicrobial Peptides Induce Growth of Phosphatidylglycerol Domains in a Model Bacterial Membrane. J. Phys. Chem. Lett. 2010, 1, 3108−3111. (66) Konshina, A. G.; Boldyrev, I. A.; Utkin, Y. N.; Omel’kov, A. V.; Efremov, R. G. Snake Cytotoxins Bind to Membranes via Interactions with Phosphatidylserine Head Groups of Lipids. PLoS One 2011, 6, e19064. (67) Efremov, R. G.; Volynsky, P. E.; Nolde, D. E.; Dubovskii, P. V.; Arseniev, A. S. Interaction of Cardiotoxins with Membranes: A Molecular Modeling Study. Biophys. J. 2002, 83, 144−153. (68) Lomize, A. L.; Pogozheva, I. D.; Lomize, M. A.; Mosberg, H. I. The Role of Hydrophobic Interactions in Positioning of Peripheral Proteins in Membranes. BMC Struct. Biol. 2007, 7, 44. (69) Zeiske, T.; Stafford, K. A.; Friesner, R. A.; Palmer, A. G. Starting-Structure Dependence of Nanosecond Timescale Intersubstate Transitions and Reproducibility of MD-Derived Order Parameters. Proteins: Struct., Funct., Genet. 2013, 81, 499−509. (70) Gilquin, B.; Bourgoin, M.; Menez, R.; Le Du, M. H.; Servent, D.; Zinn-Justin, S.; Menez, A. Motions and Structural Variability within Toxins: Implication for their Use as Scaffolds for Protein Engineering. Protein Sci. 2003, 12, 266−277. (71) Fan, H.; Mark, A. E. Relative Stability of Protein Structures Determined by X-Ray Crystallography or NMR Spectroscopy: A Molecular Dynamics Simulation Study. Proteins: Struct., Funct., Genet. 2003, 53, 111−120. (72) Cox, K.; Bond, P. J.; Grottesi, A.; Baaden, M.; Sansom, M. S. P. Outer Membrane Proteins: Comparing X-Ray and NMR Structures by MD Simulations in Lipid Bilayers. Eur. Biophys. J. 2008, 37, 131−141. (73) Korzhnev, D. M.; Religa, T. L.; Banachewicz, W.; Fersht, A. R.; Kay, L. E. A Transient and Low-Populated Protein-Folding Intermediate at Atomic Resolution. Science 2010, 329, 1312−1316. (74) Dubovskii, P. V., personal communication.

2810


Cardiotoxins: Functional Role of Local Conformational Changes

Recommend Documents