Scholarly article on topic 'Blockade of Neuronal α7-nAChR by α-Conotoxin ImI Explained by Computational Scanning and Energy Calculations'

Blockade of Neuronal α7-nAChR by α-Conotoxin ImI Explained by Computational Scanning and Energy Calculations Academic research paper on "Biological sciences"

Share paper
Academic journal
PLoS Comput Biol
OECD Field of science

Academic research paper on topic "Blockade of Neuronal α7-nAChR by α-Conotoxin ImI Explained by Computational Scanning and Energy Calculations"

OPEN 3 ACCESS Freely available online


Blockade of Neuronal a7-nAChR by a-Conotoxin ImI Explained by Computational Scanning and Energy Calculations

Rilei Yu, David J. Craik, Quentin Kaas*

Division of Chemistry and Structural Biology, Institute for Molecular Bioscience, The University of Queensland, Brisbane, Queensland, Australia


a-Conotoxins potently inhibit isoforms of nicotinic acetylcholine receptors (nAChRs), which are essential for neuronal and neuromuscular transmission. They are also used as neurochemical tools to study nAChR physiology and are being evaluated as drug leads to treat various neuronal disorders. A number of experimental studies have been performed to investigate the structure-activity relationships of conotoxin/nAChR complexes. However, the structural determinants of their binding interactions are still ambiguous in the absence of experimental structures of conotoxin-receptor complexes. In this study, the binding modes of a-conotoxin ImI to the a7-nAChR, currently the best-studied system experimentally, were investigated using comparative modeling and molecular dynamics simulations. The structures of more than 30 single point mutants of either the conotoxin or the receptor were modeled and analyzed. The models were used to explain qualitatively the change of affinities measured experimentally, including some nAChR positions located outside the binding site. Mutational energies were calculated using different methods that combine a conformational refinement procedure (minimization with a distance dependent dielectric constant or explicit water, or molecular dynamics using five restraint strategies) and a binding energy function (MM-GB/SA or MM-PB/SA). The protocol using explicit water energy minimization and MM-GB/SA gave the best correlations with experimental binding affinities, with an R2 value of 0.74. The van der Waals and non-polar desolvation components were found to be the main driving force for binding of the conotoxin to the nAChR. The electrostatic component was responsible for the selectivity of the various ImI mutants. Overall, this study provides novel insights into the binding mechanism of a-conotoxins to nAChRs and the methodological developments reported here open avenues for computational scanning studies of a rapidly expanding range of wild-type and chemically modified a-conotoxins.

Citation: Yu R, Craik DJ, Kaas Q (2011) Blockade of Neuronal a7-nAChR by a-Conotoxin ImI Explained by Computational Scanning and Energy Calculations. PLoS Comput Biol 7(3): e1002011. doi:10.1371/journal.pcbi.1002011

Editor: James Briggs, University of Houston, United States of America

Received November 9, 2010; Accepted January 5, 2011; Published March 3, 2011

Copyright: © 2011 Yu et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

Funding: This work was supported by a grant from the Australian Research Council (ARC DP1093115) ( RY is a recipient of a scholarship from the China Scholarship Council (CSC). DJC is a National Health and Medical Research Council Professorial Fellow (Grant ID 569539). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Competing Interests: The authors have declared that no competing interests exist.

* E-mail:


Nicotinic acetylcholine receptors (nAChRs) are a large family of ligand-gated ion channels that mediate rapid synaptic transmission in the central and peripheral nervous system [1,2]. nAChRs are implicated in disorders such as Alzheimer's diseases, schizophrenia, depression, hyperactivity disorders and tobacco addiction [3— 6]. All nAChRs are comprised of five homologous subunits, which are divided into a large N-terminal extracellular ligand-binding domain (LBD), a transmembrane domain, and an intracellular domain [7] (Figure 1). The nAChR subtypes include hetero- or homo-pentamers of a1-10, c, P1-4, S and/or e subunits. These subtypes differ in their pharmacological and kinetic properties, as well as their localization [8,9]. For example, the a7-nAChR is widely expressed in the brain, whereas the a3P2-nAChR is mostly expressed in the cerebellum and spinal cord [10].

Conotoxins are disulfide-rich toxins produced in the venom gland of marine cone snails [11,12]. Each of the >500 species in the Conus genus produces hundreds of different conotoxins [1315], which together form a large pool of many thousands of

bioactive peptides. Conotoxins target a diverse range of membrane receptors and ion channels to rapidly and efficiently immobilize prey [13]. The a-conotoxin family specifically and potently inhibits nAChR subtypes and, consequently, these conotoxins are useful tools in neurophysiological studies. The ability to specifically target nAChRs has also attracted interest for the development of drugs, and several conotoxins or derivatives are currently in clinical trials for the treatment of pain [16,17]. The majority of known a-conotoxins display a similar topology, as shown in Figure 1. This topology includes four cysteines arranged in a common sequence pattern -CCXmCXnC-, where X is any non-cysteine residue, and n and m are the numbers of inter-cysteine residues. Disulfide bonds connect cysteines I-III and II-IV [18,19].

ImI is one of the shortest a-conotoxins, with a loop spacing topology of m =4, n = 3 [20] and, initially, was reported to specifically interact with a7- and a9-nAChRs [21]. Later, the a3P2-nAChR was also found to be blocked by ImI [22]. ImI has been extensively studied: its structure has been determined using NMR [23-25], and its interaction with the a7-nAChR has been

Author Summary

Conotoxins are peptide toxins extracted from the venom of carnivorous marine cone snails. Members of the a-conotoxin subfamily potently block nicotinic acetylcholine receptors (nAChRs), which are involved in signal transmission between two neurons or between neurons and muscle fibers. nAChRs are important pharmacological targets due to their involvement in the transmission of pain stimuli and also in numerous neurone diseases and disorders. Their potency and specificity have led to the development of a-conotoxins as drug leads, and also to their use in the investigation of the role of nAChRs in various physiological processes. The most studied con-otoxin/nAChR system, ImI/a7, was modeled in this study, and several computational methods were tested for their ability to explain the perturbations observed experimentally after introducing single point mutations into either ImI or the a7 receptor. The aim of this study was to establish a theoretical basis to rapidly identify new a-conotoxin mutants that might have improved specificity and affinity for a given receptor subtype. Furthermore, hundreds of thousands of conotoxins are predicted to exist, and computational methods are needed to help streamline the discovery of their molecular targets.

probed by several mutational studies [26-31]. In the absence of a crystallographic structure of any nAChR, several early structural models of the binding of Iml to the LBD of a7-nAChR were generated [22,32], but they are now superseded because better templates, additional experimental data and improved modeling methods are available [33-35].

In this study, an improved model of the interaction of a7-nAChR with wild-type ImI has been developed and the structural and energetic impact of more than 30 mutations of Iml and of selected positions of the receptor were investigated. We describe for the first time a model able to explain the majority of mutation studies. Optimal methods to predict relative mutational energies were investigated, and an approach that used energy minimization produced excellent correlations with experimental values, producing R2 values of 0.74. Finally, an energy decomposition of the mutational energies was done and showed that different terms of the energy function played distinct roles. Although we focus here on conotoxin ImI, experimental mutational studies have been carried out on a range of other conotoxins, in a first step toward their development as drugs [30,31,36]. In silico mutational studies such as those described here could dramatically accelerate the development of conotoxin-based drugs and also help identify wildtype toxins with interesting pharmacological activity among the thousands of conotoxins that are predicted to exist.


In the first step of this study a model of the complex between ImI and the human a7-nAChR LBD was generated. The crystallographic structures of ImI bound to an Aplysia californica acetylcholine binding protein (AChBP), which is distantly related to nAChRs, (PDB ID: 2c9t) [34] and of bungarotoxin bound to the LBD of an isolated subunit a1 (PDB ID: 2qc1) [35] were used to build the initial model, using comparative modeling. Another crystallographic structure of ImI bound to AChBP (PDB ID: 2byp) [33] has been determined but was not used as a template because the coordinates of some amino acid atoms are missing. The secondary structure elements and the location of ImI binding site

1 5 10


Figure 1. Nicotinic acetylcholine receptor structure and a-conotoxin ImI. (A) Nicotinic acetylcholine receptors (nAChRs) are ligand gated-ion channels. Their structure is composed of a ligand-binding domain (red), a transmembrane domain (blue), and an intracellular domain (white). nAChRs are permeable to Na+ and K+ and, for some isoforms, Ca2+. The opening of the channel is triggered by acetylcholine or nicotine. One of the acetylcholine binding sites is indicated as a blue star. (B) a-conotoxin ImI comprises 12 residues and is C-terminally amidated (indicated by * in the sequence). The structure features a short a-helix and two disulfide bonds that link cysteines I-III and II-IV.


in our model are displayed in Figure 2 on the sequence and on the lowest energy model of a7-nAChR. This model was then subjected to 10 ns of molecular dynamics simulation. Thirty-four single point mutations of ImI/a7-nAChR that have been experimentally described in previous studies [26-30] were then generated in a series of models extracted over the 10 ns simulation. Finally, 14 different strategies were compared to evaluate the mutational energies of single point mutants.

Conformational variability of a7-nAChR in apo state and bound to ImI

Two series of 10 ns molecular dynamics simulations of the a7-nAChR, either in the apo state or bound to ImI, are

Figure 2. Sequence and structure of a7-nAChR ligand-binding domain. Sequence alignment of Homo sapiens a7-nAChR ligand-binding domain (LBD) (UniProtKB/SwissProt P36544), Mus musculus a1-nAChR LBD (PDB ID: 1pq1), and Aplysia californica AChBP (PDB ID: 1tg9), which is structurally analogous to nAChRs. Below the alignment, the secondary structure elements and acetylcholine binding sites are shown on the lowest energy three-dimensional model of the a7-nAChR nAChR LBD obtained by comparative modeling. Residues in the sequence alignment are numbered according to the a7-nAChR sequence. The conserved positions between the three sequences are on a dark green background, whereas the positions presenting amino acids shared by only two sequences are on a light green background. The secondary structure elements are the a-helix hi and the b-strands pi-10. The LBD is a pentamer of five subunits. The acetylcholine binding sites, indicated by star symbols, are located at the interface between the subunits. These binding sites mainly comprise the C-loop from one subunit, which is designated as the principal subunit, and the beta strands 01, 02, 03, 05' and 06 from another subunit, which is designated as the complementary subunit. The secondary structures of one subunit are highlighted in the side view, and the arrangement of the subunits and of the binding sites is shown on the top view. In the alignment, the residues of AChBP in contact with ImI in the crystal structure 2c9t are underlined in blue for positions in the principal subunit and in white for positions in the complementary subunit. doi:10.1371/journal.pcbi.1002011.g002

summarized in Figure 3. The a carbon root-mean-square deviation (RMSD) to the initial conformation became stable after 2000 ps for both simulations, indicating that they had reached equilibrium (Figure 3A,B). Indeed, the largest fluctuation, which is displayed by the third subunit, is < 1 A over the last 8000 ps of the simulation. The a-carbon root-mean-square

fluctuations (RMSF) indicate that the b-strand regions are conformationally stable, but that the C-loop and Cys-loop regions are flexible (Figure 3C,D). The dynamic property of the C-loop is particularly interesting, as the change of conformation of this loop is thought to be vital for the physiological role of nAChRs [33,37-40]. It has been shown that the interaction

of agonists with nAChRs causes the C-loop to adopt a closed conformation and this change of conformation has been hypothesized to trigger the opening of the channel [41]. According to this hypothesis, competitive antagonists stabilize the C-loop in an open conformation, potentially preventing the channel from opening. Interestingly, in our study, the C-loop in the apo model fluctuates significantly (Figure 3E), whereas the C-loop of the a7-nAChR in complex with ImI is stabilized in an open conformation (Figure 3F). It can therefore be

concluded that ImI stabilizes the C-loop in an open conformation, which, according to previous studies, should inhibit channel activity.

Molecular dynamics simulation significantly refined the conformation of the a7-nAChR/ImI model. Indeed, after 10 ns molecular dynamics, the conformation of the C-loop of the a7-nAChR/ImI model is stable and different from that of the two templates. As shown in Figure 4, the C-loop of the a7-nAChR/ ImI is more closed than the C-loop of AChBP in complex with ImI

Figure 3. Analysis of the stability of a7-nAChR over 10 ns molecular dynamic simulations in the apo (A,C,E) and ImI-bound (B,D,F) states. b strand a carbon root-mean-square deviations (RMSD) of each of the subunits over the molecular dynamics simulations to the starting frame for the apo (A) and ImI-bond models (B). a carbon root-mean-square fluctuation (RMSF) of each subunit of the apo (C) and ImI-bond (D) models. Fluctuation of the distance between the sulfur atom of a7-C190 side chain and the a carbon of a7-Y32 in the apo (E) and ImI-bond (F) models. This distance characterizes the closure of the C-loop. The RMSD is calculated using Ca atoms in b strands. The RMSD and distances were averaged using a 16 ps window.


Figure 4. Comparison of the binding site of AChBP/ImI complex (PDB ID: 2c9t), a7-nAChR/ImI complex (model, this work) and a1-nAChR/bungarotoxin complex (PDB ID: 2qc1). In the a1-nAChR/bungarotoxin structure, only one subunit was crystallized, and the bungarotoxin is not shown. The model displaying a7-nAChR was obtained by a combination of comparative modeling and molecular dynamics, and the displayed conformation corresponds to energetically minimized frames after 10 ns of simulations. The C-loop, the principal subunit, and the complement subunit are indicated. In the three first panels and from left to right, the conformation of the C-loop increasingly reduces the volume of the binding site. The fourth panel, on the right, shows a superimposition of the AChBP and nAChR subunits, highlighting the different C-loop conformations between the model and the two experimental templates. doi:10.1371/journal.pcbi.1002011.g004

but more opened than that of a1-nAChR subunit in complex with a-bungarotoxin, which is a classical antagonist of nAChR. The positions of the b-sheets are conserved between the template AChBP crystal structure and the a7-nAChR/ImI model. The h1 a-helices occupy slightly different positions, with the a7-nAChR a-helices being closer from the center of the pore than the AChBP ones (not shown).

Comparison to previous modeling and pairwise interaction studies

Our model of a7/ImI significantly differs from those [22,32] that were developed before the publication of the crystallographic structures of AChBP/ImI [33,34]. In the previous studies, models were built by homology with crystallographic structures of AChBP with the C-loop in a closed conformation, but several recent studies suggest that this C-loop conformation is incompatible with the nAChR inactive state [38,41]. Moreover, the previous studies tentatively tried to justify the binding mode of ImI using weak mutational energy couplings revealed by mutant cycle analyses, which were interpreted as pairwise interactions [28]. It proved to be impossible to reproduce all the pairwise interactions identified by this method [28]. Recently, Gleitsman et al. [42] measured similar weak mutational energy couplings occurring between residues located 60 A from each other, one being in the C-loop of an nAChR and the other in the middle of the trans-membrane domain. That study demonstrated that weak couplings are not evidence of direct interaction. On the contrary, a strong coupling was observed between a7-Y195 and ImI-R7 [28], and in our model, the side chains of these two residues are tightly packed together, as is apparent in Figure 5.

Recently Armishaw et al. [30] docked ImI into a structural model of the a7-nAChR derived by comparative modeling, using one of the AChBP/ImI crystallographic structures. Their strategy involved the mutation of a7-Y93 to Ala before performing the docking procedure, and finally the "back'' mutation of position 93 into Tyr. Presumably, the docking strategy did not succeed to place the conotoxin without this mutation step. Indeed, docking

molecules onto a structure derived by comparative modeling is a challenging task because the low accuracy of the receptor conformation either causes steric hindrance or does not allow

Figure 5. Analysis of the binding mode of ImI to the a7-nAChR.

The structure of the binding pocket occupied by ImI after molecular dynamics simulation is displayed and positions discussed in the text are highlighted. The a7 principal subunit is in orange, the a7 complementary subunit is in pale yellow, and ImI is in violet. Nitrogens are in blue, oxygens are in red and sulfurs are in yellow. doi:10.1371/journal.pcbi.1002011.g005

side chains to be tightly packed around the docked ligand [43]. The model presented by Armishaw et al. [31] is very similar to the final conformation of our molecular dynamics, despite the use of different strategies. Their model was not compared to previous experimental mutation studies, but we here provide qualitative and qualitative explanations to those mutation studies.

Structural explanation of mutational studies

The binding of Iml to the a7-nAChR has been investigated experimentally and the impact of mutations of a7 and/or Iml on the affinity (Kd) or inhibition activity (IC50) are known [26-29]. Here we investigate structural explanations for the influence of single point mutations on a7/ImI affinity through an analysis of models of the mutated complex. Mutations involving unnatural residues have not been considered here because their parameters are less refined than those for standard amino acids. The aim of our study is to compare different methods to predict the impact of single point mutations on binding affinities between conotoxin ImI and a7-nAChR; the use of unnatural residues would complicate the interpretation of those comparisons as the deviation between computed and experimental mutational energies could arise from inaccuracy in the parameters as well as from the methodological differences. The a7/ImI model will be referred to as the "wildtype model'', whereas the models of the complexes presenting mutations are referred to as ''mutated models''. Three positions of ImI, i.e., D5, P6, and R7, have been found experimentally to be important for the interaction [26]. Four receptor positions, a7-N111, a7-Q117, a7-P120 and a7-153, have some influence on the affinity of the complex but are not directly in contact with ImI in our model [28].

ImI-D5. The mutation of ImI-D5 to Asn experimentally decreases the affinity of the complex [27-29]. Residues ImI-D5, ImI-R7, a7-D197 and a7-P196 are proximate in the wild-type model, as shown in Figure 6. ImI-D5 is probably involved in charge and hydrogen bond interactions with ImI-R7, which is in turn possibly involved in both a charge and hydrogen bond interactions with the side chain of a7-D197 and the backbone of a7-P196. In the mutated model ImI-D5N, displayed in Figure 6, the side chain of ImI-R7 does not contact a7-D197 and a7-P196 as it does in the wild-type model. Presumably, ImI-D5 plays a significant role to stabilize the conformation of ImI-R7, which allows ImI-R7 to interact with both a7-D197 and a7-P196. Thus, the disruption of the interaction between ImI-D5 and ImI-R7, which is not at the interface, indirectly causes a decrease in affinity by weakening interface interactions between ImI-R7, a7-P196 and a7-D197.

ImI-P6. According to our model, ImI-P6 is tightly packed in a hydrophobic pocket formed by a7-W55, a7-Y93, a7-L119 and a7-W149 (Figure 5). Mutations of a7-Y93, a7-W149 and to a lesser extent a7-L119 to a partly hydrophobic Thr residue decreased ImI affinity, and accordingly, the mutated models displayed reduced packing around ImI-P6 (not shown). Consistent with this analysis, an increase in affinity was achieved by introducing a Pro derivative with increased hydrophobicity [30]. Mutations of ImI-P6 to Gly, Ala or Val reduced the affinity of ImI for a7 [26,28], which is also consistent with the models. However, caution should be exercised in the interpretation of those mutations, as the mutation of P6 induces dramatic conformational changes of ImI backbone [44]. The method we used to model the mutated models cannot take into account such dramatic conformational changes, and therefore the study of ImI-P6 mutants was not carried out.

ImI-R7. In the wild-type model, the aliphatic part of the ImI-R7 side chain is in contact with a7-Y93, a7-Y195 and ImI-P196,

whereas the positively charged guanidinium moiety is proximate to the negatively charged headgroups of a7-D197 and ImI-D5 (Figure 6). Mutation of ImI-R7 to Gln breaks the charge interactions with a7-D197 and ImI-D5, which is consistent with the decrease in affinity observed experimentally [28]. Moreover, as discussed previously, the mutation of ImI-D5 to Asn also influences the conformation of R7. The mutation of a7-D197 to Asn decreases the affinity for ImI [28], corresponding to a loss of charge interaction (not shown).

The involvement of a7-Y195 in van der Waals interactions with R7 is supported by the observation that a7-Y195T decreases the binding affinity, whereas a7-Y195F does not [28]. An interaction between ImI-R7 and a7-Y195 has also been deduced by double mutant cycle analysis [28].

Mutations of a7-R186 to Ala, Glu, Gln and Val increased the affinity [28]. In the wild-type model, a charge repulsion interaction occurs between a7-R186 and ImI-R7 (Figure 6), and this unfavorable interaction is removed by mutating position 7 into a non-charged or negatively charged residue.

a7-N111. a7-N111 is the last position of the P5'-strand and is not in direct contact with ImI. The model of the mutant a7-N111S (Figure 6) features a longer P6 strand and a change in conformation of the P5'-P6 loop compared to the wild-type. Several side chains located in the binding site, including R79, Q117 and H115, have a slightly different orientation in the mutated models. Although our models show that position 111 has an influence on the binding of ImI, it is difficult to provide simple qualitative explanations of the increase in affinity measured experimentally [28].

a7-Q117. Mutation of a7-Q117 to Ala or Ser increases affinity for ImI [28]. The a7-Q117A model shows that the conotoxin is closer to the backbone of the P5' and P6-strands compared to the wild-type (Figure 6) and the buried surface area is increased by 80 A2. This mutation decreases the size of side chain at position 117 and therefore allows ImI-W10 to have better packing at the interface. The increased affinity resulting from the mutation a7-Q117S, which also decreases the size of 117 position side chain, is explained similarly.

a7-P120. Mutation of a7-P120 to Ala decreases the affinity [28]. As shown in Figure 6, mutation of a7-P120 to Ala caused a local conformational change of the backbone, which resulted in the rearrangement of the neighboring side chain at position 119. In the wild-type model, the side chain of a7-Ile-119 closely stacks with ImI-P6, but in the mutated model the side chain of a7-Ile-119 has fewer contacts with ImI. Thus, a7-P120 indirectly contributes to the affinity by influencing the conformation of a side chain at the interface.

a7-G153. Mutation of position a7-153 from a Gly to a Ser, a bulkier residue, results in a drastic decrease in affinity [28]. In the mutated model the C-loop adopts a more open conformation than in the wild-type model. It is likely that a steric exclusion between (7-S153 and (7-P194 forces the C-loop to change its position relative to the binding site, decreasing the number of interactions at the interface, and therefore accounting for the drop in affinity measured experimentally.

Comparison of methods to compute mutational energies

Mutational energies of single point mutants were computed using two energy functions: molecular mechanics generalized Born surface area (MM-GB/SA) and molecular mechanics Poisson-Boltzmann surface area (MM-PB/SA) energy functions. The mutated models were first refined using either the minimization based approach (MBA) or the molecular dynamics simulation based approach (MDBA). For the MBA, mutations were

Figure 6. Superimposition of wild-type and mutated models. The models of the mutants shown were refined using molecular dynamics and the conformations shown in this figure are the last frames of the molecular dynamics trajectories. The arrows highlight local conformational changes. doi:10.1371/journal.pcbi.1002011.g006

introduced in 15 frames extracted from the wild-type 10 ns molecular dynamics simulation and the mutated models were minimized using either explicit water (EWM) or a distance-dependent dielectric constant (DDDCM). For the MDBA, mutations were created on the model of the last frame of the wild-type 10 ns molecular dynamics simulation, and the mutated models were subjected to 500 ps molecular dynamics. Because the complex showed only small conformational variations during the last 8 ns of the simulation, the last frame can be chosen as representative of the wild-type structure. The energies were

averaged on the minimized mutated models for MBA, and on 50 frames extracted from the last 100 ps of the 500 ps molecular dynamics for MDBA.

Energy predictions using MBA. The energies of 16 ImI single point mutants were computed (Table 1) and are compared with experimental values in Figure 7. Using the simple DDDCM, MM-GB/SA gave better predictions of mutational energies than MM-PB/SA, as shown by the correlation coefficient R2 of 0.71 and 0.58 between experimentally derived energies and energies computed with generalized Born and Poisson-Boltzmann,

Table 1. Calculated and experimental mutational energies (kcal/mol) of Iml mutants.


mutant EXPC

S4A 0.2 ±0.4 0.8 ± 0.5 0.1 ±0.3 0.3 ±0.3 3.2 ±1.2 9.9 ± 1.4 2.5±1.1 10.4: ±1.1 0.8±0.9 4.9 ±1.0 -0.3 ±0.4 —4.4±0.3 4.3 ±0.3 7.6±0.3 0.3 ±0.2

D5A 3.8±0.9 10.5 ±1.5 3.7±0.6 4.7±0.7 3.8 ±1.0 22.5: ±1.2 2.0±0.8 10.1: ±1.4 4.7±1.0 4.5 ±1.0 4.3 ±0.4 8.3 ±0.4 2.2 ±0.3 4.7±0.2 2.1 ±0.1

D5K 29.6±3.7 50.9 ±5.3 12.3 ±1.6 1 7.0 ±1.8 8.3 ±1.4 20.4: ±1.6 —0.4±1.2 8.1 ± 1.2 — 1.9± 1.0 2.2 ±1.6 37.4±0.6 22.9 ±0.3 20.6±0.2 12.6±0.3 4.0±0.2

D5N 6.0±0.9 8.1 ± 1.4 3.5 ±0.8 2.2 ±0.9 10.9±1.4 23.8: ±1.4 7.6 ±1.2 10.6: ±1.2 6.6 ±1.2 1.6 ± 1.4 11.8±0.4 10.7 ±0.4 5.7±0.3 4.0±0.3 2.8±0.2

R7A 11.9±1.3 18.4 ±1.3 14.0±0.8 16.8±0.9 8.5 ±1.7 13.2: ±1.5 10.1 ±1.3 19.2: ±1.2 14.5±1.0 15.6±1.0 42.5 ±0.3 16.4 ±0.4 15.9±0.3 18.4±0.2 3.0±0.3

R7E 26.8 ±1.6 66.0 ±3.6 15.5 ±1.0 33.0±2.6 0.2 ±1.3 21.1 : ±1.8 14.4 ± 1.3 43.0: ±1.3 10.0±1.0 20.3 ±1.1 36.3 ±0.4 79.2 ±0.3 18.1 ±0.3 36.4±0.4 3.9±0.2

R7K 11.0±1.0 14.0 ±1.3 11.2 ±0.8 13.4±1.4 13.3±1.3 16.9: ±1.3 9.9 ±1.2 15.7: ±1.1 9.8 ±1.0 8.4 ±1.2 23.1 ±0.3 26.8 ±0.4 7.9±0.4 5.5 ±0.4 2.9±0.2

R7L 9.4±0.8 17.6 ±1.6 10.6 ±1.2 16.5 ±1.2 9.0 ±1.4 17.9: ±1.6 5.2 ±1.0 15.8: ±1.4 3.0±1.2 6.0 ±1.3 13.8±0.3 36.2 ±0.3 9.7±0.3 12.8±0.4 3.3 ±0.6

R7Q 13.5±1.2 24.5 ±1.8 8.1 ±1.2 13.1 ±1.3 2.0± 1.5 15.1 : ±1.3 1.2±1.1 9.7 ± 1.4 6.0 ±1.0 3.1 ±1.1 3.3±0.4 5.7±0.4 9.1 ±0.3 12.1 ±0.3 3.5 ±0.2

A9S -0.5 ±0.3 2.2 ± 0.7 -4.5 ±0.8 -2.7 ±0.9 3.6±1.3 6.0 ± 1.1 -3.2 ±1.0 -0.7 ±1.0 -1.6±1.1 -5.1 ±1.0 — 6.8±0.3 — 7.0±0.3 —2.3±0.3 —2.8±0.2 -0.1 ±0.2

W10A 3.5±0.7 7.3 ± 0.8 6.5 ±0.6 3.7±0.7 7.1 ±1.4 13.9: ±1.5 2.8±1.2 8.0 ± 1.2 3.2 ±1.0 -1.1 ±1.0 3.0±0.3 —6.9±0.4 6.3 ±0.3 3.2 ±0.4 1.1 ±0.1

W10F 1.3 ±0.9 2.3 ± 0.4 3.1 ±0.5 1.9±0.6 — 0.5 ±1.3 7.3 ± 1.4 —2.6±1.2 2.6± 1.2 -3.0±1.3 — 5.5 ±1.3 1.7±0.3 -9.3 ±0.5 5.3 ±0.4 1.2 ±0.3 0.8±0.5

W10T 4.3±0.7 2.1 ± 0.6 4.8±0.5 1.7±0.3 4.7 ±1.3 11.0: ±1.2 5.8±1.1 4.9 ± 1.1 1.9±1.2 — 3.9± 1.0 13.9±0.4 4.3 ±0.4 6.8±0.3 1.6±0.4 2.0±0.2

W10Y 1.1 ±0.4 1.5± 0.5 —0.7±0.7 -0.9 ±0.9 -0.7 ±1.2 12.7: ±1.3 —2.8±1.1 2.2 ± 1.0 -1.1 ±1.0 0.7±0.8 —14.8±0.3 —23.0±0.4 -0.3 ±0.2 —3.2±0.2 0.1 ±0.3

R11A 1.0±0.6 3.8± 0.7 2.2 ±0.8 3.9±0.7 4.0 ±1.5 6.0 ± 1.2 2.9±1.3 10.2: ±1.1 —2.8± 1.0 -4.5 ±0.9 -16.2 ±0.4 —2.8±0.4 1.1 ±0.3 2.6±0.3 -0.2 ±0.2

R11Q 1.1 ±0.6 4.2 ± 0.9 3.1 ±0.9 3.7±0.7 4.8± 1.1 21.3: ±1.2 4.4 ±1.0 7.9 ± 0.9 -4.3 ±1.1 -4.9 ±0.8 —5.5±0.5 —8.6±0.4 2.7±0.3 2.0±0.4 0.3 ±0.2

R2 0.71 0.58 0.74 0.64 0.16 0.40 0.30 0.34 0.44 0.45 0.70 0.59 0.71 0.51

Minimization based approach (MBA), which uses either a distance-dependent dielectric constant minimization (DDDCM) or an explicit water minimization (EWM). The standard deviations of the mutational energies were computed using 15 frames.

bMolecular dynamics simulation based approach (MDBA). Five different strategies, noted (i) to (v), were used to refine the conformation of the structural model, (i) all receptor atoms >6 A from the conotoxin were restrained to their position; (ii) all receptor atoms >4.5 A from the conotoxin were restrained to their position; (iii) all the atoms located >6 A from the mutated residue were restrained to their position; (iv) all the atoms from the receptor were restrained to their position, and all the atoms from the conotoxin mutants were free to move; (v) all residues were restrained to their position. The standard deviations of the mutational energies were computed using 50 frames. Experimental mutational energies (EXP). The values were derived from experimental Kd values [26-28] using the equation AAGbinding = EXP = —RT ln[Kd (mutant)/K(J (wild-type)] at 298 K. dMutational energy computed using molecular mechanics generalized Born (GB) surface area (MM-GB/SA) at 298 K. eMutational energy computed using molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) at 298 K. doi:10.1371 /journal.pcbi.1002011 ,t001

Figure 7. Correlation between experimentally derived and calculated mutational energies of the ImI mutants in the minimization based approach MBA. Mutational energies were computed using either molecular mechanics generalized Born (GB) surface area (MM-GB/SA) or molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) energy functions at 298 K. The mutated models were refined using MBA with either distance dependent dielectric constant minimization (DDDCM) or explicit water minimization (EWM). Experimental mutational energies (AAG Exp) were derived using the corresponding Kd values of ImI wild-type/(7-nAChR and ImI mutants/(7-nAChR [26-28]. doi:10.1371/journal.pcbi.1002011.g007

respectively. The more computationally intensive EWM led to a subtle improvement for both energy functions, and MM-GB/SA still performed better than MM/PB-SA (R2 of 0.74 and 0.64, respectively).

It is a priori surprising that the generalized Born method gave better predictions because generalized Born parameters were optimized to reproduce the more computationally demanding results from finite-difference Poisson-Boltzmann [45]. A possible explanation for this phenomenon is that the approximations introduced in the generalized Born method are able to partly ameliorate the inaccuracy of theoretical models, whereas the Poisson-Boltzmann method cannot. Indeed, it has been shown that MM-GB/SA provides better ranking of models generated by docking than MM-PB/SA [46-48]. Few previous studies of protein/peptide complexes used simultaneously MM-PB/SA and MM-GB/SA [49] and to our knowledge none made an extensive comparison of the ability of the two energy functions to rank mutational energies.

Regarding the refinement method, the small superiority of EWM over DDDCM mainly arises from the prediction of two mutational energies, R7E and D5K, which both result in a reversal of charges. Therefore, models minimized using explicit water representation seem to be slightly more accurate that the ones obtained from implicit solvation methods, like DDDCM. Nevertheless, it should be noted that energy computations using DDDCM were, on average, four times faster than EWM.

Energy predictions using MDBA. For the MDBA, short molecular dynamics simulations were carried out to refine the models of the mutants, and the mutational energies were evaluated using MM-GB/SA or MM-PB/SA. The resulting energies are shown in Table 1 and are compared to experimentally derived energies in Figure 8.

To achieve the simulations within a practical computational time, only short simulations could be carried out for each mutant. Because of this short simulation time, the molecular dynamics trajectories only partially sample the accessible conformational space. As a result, the conformations of the positions that are not influenced by the mutations are sampled differently in the simulation of the wild-type and mutant complexes, and those differences create artificial noise in the computation of the mutational energies. To overcome this problem, some atoms were restrained to their initial location by a quadratic force, and five strategies were employed to select the atoms that should be restrained. In strategy (i), all of the receptor atoms located further than 6 A from the conotoxin were restrained to their position. This led to a poor correlation with experimental values, with R2 = 0.16 for MM-GB/SA and R2 = 0.40 for MM-PB/SA (Figure 8). To reduce the possible detrimental influence of positions located out of the binding site, additional restraints were added to the system in strategy (ii) by lowering the distance cut-off to 4.5 A. This change improved the prediction made with MM-GB/SA (R2 = 0.30) but worsened the one computed using MM-

Figure 8. Correlation between experimentally derived and calculated mutational energies of ImI mutants in the molecular dynamics based approach MDBA. Mutational energies were computed by using either molecular mechanics generalized Born (GB) surface area (MM-GB/SA) or molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) approaches at 298 K. In the MDBA five alternative position restraint strategies were employed: (i) all receptor atoms >6 A from the conotoxin were restrained to their position; (ii) all receptor atoms >4.5 A from the conotoxin were restrained to their position; (iii) all the atoms located >6 A from the mutated residue were restrained to their position; (iv) all the atoms from the receptor were restrained to their position, and all the atoms from the conotoxin mutants were free to move; and (v) all residues were restrained to their position. Experimental mutational energies (AAG Exp) were derived using the corresponding Kd values of ImI wild-type/a7-nAChR and ImI mutants/a7-nAChR [26-28]. doi:10.1371/journal.pcbi.1002011.g008

PB/SA (R2 = 0.34). The change in cut-off allowed the decrease in the number of atoms without restraints from about 1000 to 800. In strategy (iii), the mutated residue was taken as the center of the distance cut-off. With a 6 A distance cut-off, only about 400 atoms were without restraints. The correlation with experimental values improved for both MM-GB/SA and MM-PB/SA compared to strategy (ii), with R2 reaching 0.44 and 0.45, respectively. In strategy (iv), all the atoms from the receptor were restrained to their position, but all the atoms from the conotoxin mutants (about 170 atoms) were free to move. Surprisingly, the addition of restraints on the whole receptor significantly improved the correlation of both MM-GB/SA and MM-PB/SA, as the R2 reached 0.70 and 0.59, respectively. In the last strategy (v), all residues were restrained, which allowed only subtle changes to the atom positions. The agreement with experimental binding energies was similar to strategy (iv) as R2 = 0.71 for generalized Born and R2 = 0.51 for Poisson-Boltzmann. By only allowing the atoms to have local moves, the conformational sampling of strategies (iv) and (v) have similar effect to energy minimization, and the correlation with experimental energies are closer to the MBA than strategies (i)-(iii). Nevertheless, despite being considerably more time consuming, the molecular dynamics refinement approach used in the MDBA did not produce a better prediction than the minimization approach used in the MBA. We therefore conclude that the MBA is the best method to use.

Mutational energy of receptor mutants. To further validate the accuracy of the MBA, 48 additional mutational energies using DDDCM and EWM were computed by mutating positions of the receptor (Table 2). The correlation coefficient between the experimental and predicted mutational energies does not change significantly by including the predictions made for the

receptor mutants (Figures S2). The stability of the correlation coefficient upon addition of new data demonstrate that our models and methods could be used to predict the relative binding affinities of other single point mutations of the ImI/a7-nAChR system.

Binding energy and mutational energy decompositions

To better understand the energetic components stabilizing the ImI/a7-nAChR complex, the free energies of the system and the mutational free energies of the mutated complexes, computed using EWM, were decomposed into entropic, electrostatic, van der Waals, and hydrophobic contributions. The solute entropic contribution to the binding energy has been neglected in our previous calculations, but it was estimated using normal-mode analysis in this section. As shown in Table 3, the van der Waals interactions and the hydrophobic effects stabilize the complex, whereas the electrostatic contribution is destabilizing. The observation that the van der Waals interactions and hydrophobic effect are predominant over the electrostatic interactions correlates with a statistical analysis of interface features carried out over the 10 ns of the molecular dynamics simulation. During the simulation, the average buried surface area of the wild-type complex was of 1150 A2, which is twice as large as the average 500 A2 of peptide/protein interfaces [50], and can be associated with the important van der Waals and hydrophobic effect energies. ImI and a7-nAChR form, on average over the simulation, three hydrogen bonds and one salt-bridge, and this small number of electrostatic interactions is consistent with average values for a-helical peptides [50].

The decomposition of the mutational free energies are displayed in Table 4 and the correlation between different contributions and experimentally derived mutational energies are shown in Figure 9.

Table 2. Calculated and experimental mutational energies (kcal/mol) of the receptor mutants.


R2 0.71 0.58 0.64 0.74

aDistance dependent dielectric constant minimization (DDDCM). The standard deviations of the mutational energies were computed using 15 frames. bExplicit water minimization (EWM). The standard deviations of the mutational energies were computed using 15 frames.

Experimental mutational energies (EXP) that were derived from experimental Kd values [26-28] using the equation AAGbinding = EXP = —RT ln[Kd (mutant)/ Kd (wild-type)] at 298 K.

dMutational energy computed using molecular mechanics generalized Born (GB) surface area (MM-GB/SA) at 298 K.

eMutational energy computed using molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) at 298 K. doi:10.1371/journal.pcbi.1002011.t002

The electrostatic contribution has by far the best agreement with experimental energies (R2 = 0.62 with MM-GB/SA) and is therefore the major contributor to the specificity between mutants. The van der Waals interactions also participates, but to a lesser extent (R2 = 0.35). On the contrary, the hydrophobic effect does not allow differentiation of the mutants (R2 = 0.01). With a correlation coefficient of 0.23, the solute entropy term has only a small influence on specificity. Furthermore, the correlation coefficient between the experimentally derived and predicted mutational energies does not change substantially after including the solute entropic component (Figure 7 and 9). This justifies the proposal that solute entropic contributions, which are computationally demanding, can be neglected when predicting the ranking of single point mutant binding affinities based on the computation of binding energies.

Concluding remarks

In this study, an extensive computational analysis of a complex between an a-conotoxin, Iml, and an nAChR, a7-nAChR, was carried out. In the absence of the crystal structure of a complete nAChR extracellular domain, modeling the interaction of an inhibitor of an nAChR is difficult. Nevertheless, we successfully studied the binding mode between a-conotoxin Iml and a7-nAChR using a combination of comparative modeling and molecular dynamics simulation. Using this model, we have explained the effect of mutations described in previous experimental studies.

The structures of 16 mutated ImI/a7-nAChR complexes were refined using MBA or MDBA, and the binding energies were predicted using MM-PB/SA and MM-GB/SA. To our knowl-

edge, this study constitutes the first attempt to use these energy functions to study the binding of a range of a-conotoxin variants to an nAChR. The approach using a simple minimization to refine the model (MBA) led to the best agreement between predicted mutational energies and experimental values.

Another important conclusion of our study is that affinity between ImI analogues and a7-nAChR was mainly governed by van der Waals and non-polar desolvation energies, whereas the electrostatic interactions were mainly important for the specificity. Interestingly, the entropy had little influence on the mutational energy of single point mutants. Because a-conotoxins share the same tightly packed structural fold, our observations on the energy decomposition are likely to help in the rational optimization of a-conotoxins pharmacological properties in general.

In order to perform extensive computational scanning of a-conotoxins, a fast and accurate approach is necessary. In this respect, we have identified that the best method to achieve this goal is to refine the mutated models by minimization using explicit solvation and to compute mutational energies using MM-GB/SA.

Materials and Methods

Comparative modeling of a7-nAChR

In the absence of the crystal structure of a7-nAChR, the human a7-nAChR LBD was modeled using a comparative approach, following a strategy described previously [51-53]. The crystal structure of an isolated Mus musculus muscle type extracellular domain of the a1-nAChR subunit in complex with the inhibitor a-bungarotoxin was solved at 1.94 A resolution (PDB ID: 2qc1) [35]. The Mus musculus a1-nAChR subunit shares 38% sequence identity with the Homo sapiens a7 subunit and superimposes with, on average, 2.9 A rmsd with the AChBP subunits. An electron microscopy structure of a complete muscle type nAChR of Torpedo marmota (PDB ID: 2bg9) revealed a similar arrangement of subunits as the one presented by AChBP. As the electron microscopy structure is of low resolution (4 A), the AChBP structures (PDB ID: 2c9t) were employed as structural templates in our comparative modeling strategy to orient the five a7 subunits in the pentamer. The orientations of the side chains were modeled according to the a1 template (PDB ID: 2qc1) due to its overall higher sequence identity to the a7 subunit than AChBP. A sequence alignment between the two structural templates and a7-nAChR LBD is displayed in Figure 2. The secondary structure elements and ligand-binding sites observed on the experimental structures and predicted for a7-nAChR are also shown in Figure 2. The modeling of the nAChR C-loop required special attention as its change in conformation allows the binding site to accommodate ligands of different sizes [33]. In our model, the structure of AChBP in complex with ImI (PDB ID: 2c9t) was used to derive restraints in the C-loop region because AChBP has locally higher sequence identity than a1, and because the C-loop conformation in the AChBP structure allows ImI to fit in the binding site. Conversely, the Cys-loop, the P1-P2 loop, the A-loop, and the B-loop were modeled using information from the a1 template because it displays higher sequence identity than AChBP. Multiple sequence alignment between Aplysia californica AChBP and the LBD of a1, a2, a3, a4, P1, P2, a6, a7, a9 and a10 was generated using MUSCLE with default parameters [54]. The multiple alignment between AChBP, a7 and a1 was then manually adjusted based on structural superimpositions of the crystal structures of AChBP (PDB ID: 2c9t) and a1 (PDB ID: 2qc1). The comparative modeling program Modeller [55] (version 9v7) was then employed to generate 100 three-dimensional structural models of the a7-nAChR complex. The Cys-loop region was

Table 3. Decomposition of the binding free energy (kcal/mol) of Iml and mutants.

Mutant ent

Generalized Born GB

Poisson-Boltzmann PB


29.2 ±1.4 26.6 ±1.6 25.4±1.2 26.0±2.2

24.1 ±0.8

25.2 ±1.7 27.0±1.1 25.5 ±1.5 26.7±1.3 25.5 ±1.3

25.2 ±1.2 26.1 ±1.7 27.1 ±1.7 27.8±1.6

27.3 ±1.4 27.1 ±1.6 24.7 ±1.2

-97.0±0.9 -94.0±0.8 -95.9±0.7 -95.7±0.8 -86.3±0.6 -89.8±0.8 -93.3±0.8 -90.9±0.7 -91.6±0.7 -97.9±0.5 -91.3±0.6 -95.6±0.7 -93.6±0.8 -96.1 ±0.9 -95.0±0.8 -95.3±0.9 -97.2±0.8



-432.1 ±9.2



-49.2 ±12.7

-283.1 ±10.3


-136.3 ±10.7




-267.2 ±8.8





294.8 ±7.5 463.7 ±7.4 462.7 ±8.6

674.2 ±9.5

161.3 ±14.0 85.3 ±11.0 318.6±8.6

162.5 ±10.7 166.7 ±9.2 305.7±7.7

289.6 ±8.0 290.6 ±7.8 296.2 ±7.7 301.0±7.0

193.9 ±5.2 197.9 ±5.6 299.4±7.1

-12.4 ±0.1 -12.5 ±0.1 -12.6 ±0.1 -13.2 ±0.1 -11.2±0.1 -12.3 ±0.1 -12.4 ±0.1 -11.8±0.1 -12.1 ±0.1 -12.8 ±0.1 -11.8±0.1 -12.3 ±0.1 -12.1 ±0.1 -12.7 ±0.1 -12.1 ±0.1 -12.3 ±0.1 -12.5 ±0.1

28.0±2.3 28.8 ±1.7 30.6 ±1.6 39.8 ±1.6 30.0±2.5

36.2 ±2.9 35.5 ±2.4 31.9±2.3

30.3 ±2.3 24.7±2.3 28.2±2.1 29.6±2.1

29.1 ±2.1 26.6±2.2 27.8±1.7 29.3 ±1.7

28.2 ±2.0

282.4 ±6.9

451.2 ±6.8 450.1 ±7.7 661.0 ±8.0 150.1 ±15.5 73.0 ±12.5

306.1 ±9.2

150.7 ± 11.6 154.6± 10.1 292.9 ±7.7 277.8±7.9

278.3 ±7.8

284.2 ±7.7

288.4 ±6.9

181.8 ±5.3 185.5±5.8

286.9 ±7.2

296.7±7.0 466.5 ±7.7 462.9±8.7 680.2 ±10.2 166.2 ±13.4

104.2 ± 11.2 322.4±8.5 170.1 ±10.8

173.3 ±9.4 308.9 ±7.6 288.5 ±7.9 290.9 ±7.9 294.8 ±7.8 302.3±7.0 197.3 ±4.9 200.1 ±5.7 301.0±7.2

-12.1 ±0.1 -12.3±0.1 -12.2±0.1 -12.6±0.1 -11.2±0.2 -11.8±0.1 -12.0±0.1 -11.5±0.2 -11.8±0.1 -12.2±0.1 -11.4±0.1 -11.9±0.1 -11.7±0.1 -12.1 ±0.1 -11.8±0.1 -12.0±0.1 -12.1 ±0.1

30.0±2.5 31.5 ± 1.7 30.8± 1.7 45.8±2.4 34.9±2.5 55.1 ±2.5 39.3±2.4 39.5±2.5 37.0±2.4 28.0±2.6 27.0±2.3 30.0±2.2 27.6±2.3 27.9±2.3 31.3±2.1 31.5±2.0 29.8±2.1

284.6±6.6 454.7±7.1 450.8±7.7 667.6±8.8 155.0±14.9 92.5 ±12.7 310.3±9.2 158.5 ± 11.9 161.6±10.3 296.7 ±7.89 277.1 ±7.9 279.0±8.1

283.1 ±7.7

290.2 ±7.0 185.5 ±5.0 188.1 ±5.9 288.9 ±7.3

-52.2±1.6 -51.2± 1.5 -52.5±1.5 -43.1 ±1.5 -43.3 ±1.6 -40.7±2.2 -43.2 ±1.4 -45.4± 1.7 -46.6± 1.5 -60.5 ±1.6 -49.8± 1.6 -52.3±1.5 -49.5 ±1.4 -54.3 ±1.6 -52.0± 1.3 -51.2± 1.3 -56.7 ±1.6

-50.0±2.1 -48.2 ±1.6 -51.9±1.6 -36.5±±1.i -38.5±2.0 -21.3±3.8 -39.0±1.6 -37.5±2.2 -39.6±1.8 -56.7±2.1 -50.6±1.8 -51.5 ± 1.6 -50.6±1.7 -52.6±1.8 -48.3 ±1.8 -48.6 ±1.6 -54.7±1.8

decomposition of the binding free energy (AG) into entropy (ent), van der Waals (vdw), Coulombic electrostatic (Coul), polar desolvation (epol) and non-polar desolvation (SA) components. Two energy functions were used: molecular mechanics generalized Born (GB) surface area (MM-GB/SA) and molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA). The temperature was set at 298K. The standard deviations of the binding free energies were computed using 15 frames. bAG(GB) is the free energy computed using MM-GB/SA. CAG(PB) is the free energy computed using MM-PB/SA. dElectrostatic interaction (ele = epol + Coul).

eDesolvation contribution to the the binding free energy (des = epol + SA). doi:10.1371/journal.pcbi. 1002011 .t003

Table 4. Decomposition of the mutational energies (kcal/mol) of Iml mutants.

Generalized Born GB Poisson-Boltzmann PB

Mutant vdwa enta elea SAa elea SAa AAG (GB)b AAG (PB)c

S4A 0.1±0.6 4.561.6 -0.260.2 0.160.1 0.260.2 0.060.2 4.661.4 4.861.5

D5A 3.260.6 1.961.4 0.660.3 -0.160.1 1.760.4 -0.260.1 5.661.8 6.562.1

D5K 1.560.8 1.361.0 11.661.6 -0.760.2 16.062.2 -0.560.2 13.661.6 18.363.3

D5N 1.260.6 0.761.8 2.461.0 -0.160.2 1.060.4 -0.160.1 4.261.3 2.961.8

R7A 10.960.7 -0.661.7 1.860.9 1.360.2 5.161.1 0.960.3 13.461.4 16.362.5

R7E 7.4 60.9 0.561.7 8.061.8 0.2 60.1 25.361.6 0.3 60.1 16.0 62.0 33.5 6 2.9

R7K 3.960.7 2.361.8 7.360.9 0.160.1 9.560.9 0.160.1 13.661.7 15.862.3

R7L 6.360.5 0.862.0 3.760.9 0.760.2 9.760.9 0.660.2 11.461.8 17.262.3

R7Q 5.660.7 2.061.7 2.161.1 0.460.1 7.261.0 0.460.1 10.162.3 15.161.8

A9S -0.760.6 0.861.7 -3.561.1 -0.360.1 -1.860.4 -0.160.1 -3.761.9 -1.961.5

W10A 5.860.5 0.561.8 0.060.5 0.760.2 -2.860.3 0.760.2 7.061.7 4.261.1

W10F 1.560.5 1.462.2 1.460.3 0.260.1 0.260.3 0.260.1 4.561.3 3.261.4

W10T 3.560.6 2.462.3 0.960.4 0.460.1 -2.260.5 0.460.1 7.362.4 4.262.4

W10Y 1.060.6 3.161.9 -1.660.9 -0.260.1 -2.060.4 0.060.1 2.462.0 2.261.9

R11A 2.260.5 2.661.2 -0.460.9 0.460.2 1.461.0 0.360.2 4.761.3 6.561.7

R11Q 1.860.6 2.461.8 1.160.8 0.260.1 1.760.8 0.160.1 5.561.9 6.161.9

decomposition of the mutational energy (AAG) into van der Waals (vdw), entropic (ent), electrostatic (ele) and non-polar desolvation (SA) components. Two energy functions were used: molecular mechanics generalized Born (GB) surface area (MM-GB/SA) and molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA). The temperature was set at 298K. The standard deviations of the mutational energies and energy components were computed using 15 frames. bAAG(GB) is the mutational free energy computed using MM-GB/SA. cAAG(PB) is the mutational free energy computed using MM-PB/SA. doi:10.1371/journal.pcbi.1002011.t004

modeled using the a1 subunit template only, whereas the C-loop and B-loop were modeled based on AChBP. The model selected according to the DOPE score [56] was analyzed using MolProbity [57] and 94% residues were in the favorable region of the Ramachandran plot, which is acceptable for a comparative model

ImI/a7-nAChR complex

A model of the structure of the complex ImI/a7-nAChR was obtained by comparative modeling. An X-ray diffraction structure of the complex between ImI and AChBP (PDB ID: 2c9t) was used to provide restraints between ImI and the nAChR and also structural restraints to ImI conformation. The structure of a7-nAChR was modeled using the same sequence alignment/ structure described previously to model the apo state. The use of a comparative modeling approach is justified by the fact AChBP and a7-nAChR are likely to have very similar binding modes because they share a high level of sequence identity in their binding sites (52% identity according to alignment in Figure 2).

Molecular dynamics simulation

Molecular dynamics simulations (MD) were performed using Gromacs 3.3.1 package [55] and the 53a6 forcefield. The ImI/a7-nAChR model was solvated in a cubic box with an edge length of 11.4 nm solved by adding 40,773 SPC water molecules. 74 Na+ and 27 Cl— ions were added to simulate a physiological NaCl concentration of 0.1 M and to neutralize the system. The system was minimized using 1000 steps of steepest descent algorithm. The temperature was progressively raised from 0 K to 300 K over 100 ps of constant pressure and temperature (NPT) MD simulation with all the protein atoms restrained to their initial position. Ten nanosecond NPT MD was then performed on the

whole system without restraints with Berendsen temperature bath coupling set at 300 K and an isotropic molecule based scaling setup at 1 atm [58]. The electrostatic interaction between non-covalent atoms was computed with particle-mesh Ewald method [59] with a distance cutoff of 10 A. The LINCS algorithm [60] was used to constrain all bonds and the time step of the simulation was set to 2 fs. The simulation of the apo state a7-nAChR was prepared using the same procedure. Ten nanosecond MD simulations were performed twice for ImI/a7-nAChR and three times for the a7-nAChR in apo state systems. Stability and conformational variabilities of those five simulations are provided in Figure 3 and in supplementary material figure S1.

Models of single point mutants

In the MBA, 15 frames were extracted every 500 ps in the interval between 3-10 ns of the 10 ns molecular dynamics simulation trajectory of the wild-type model. Each frame was minimized using AMBER10 [61] (with the AMBER ff03 forcefield) by 2000 steps of steepest descent algorithm followed by 2000 steps of conjugate gradient algorithm with the backbone of the complex restrained. Two thousand five hundred steps of steepest descent minimization and 2500 steps of conjugate gradient minimization were then performed without restraints. The side chains in the ligand were mutated using Modeller and all the residues (including residues of the ligand) were minimized using AMBER10 [61]. In the DDDCM approach, e = 4r was used, whereas in the EWM approach, the protein was solvated in a water box with a minimum of 8 A between the solute and the side of the box.

In the MDBA, a water cap with a radius of 16 A from the center of the binding pocket was added. MD was only performed on the mutants of the last frame obtained from MBA above.

Figure 9. Correlation between calculated mutational energy components of ImI mutants and experimentally derived mutational energies. The explicit water minimization approach was employed to compute the Gibbs free energy (AAG) using either molecular mechanics generalized Born (GB) surface area (MM-GB/SA) or molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) approaches at 298 K. The energies were decomposed into van der Waals (vdw), electrostatic (ele), surface area (SA, only shown for GB) and entropic components (not shown). Experimental mutational energies (AAG Exp) were derived using the corresponding Kd values of ImI wild-type/a7-nAChR and ImI mutants/a7-nAChR [26-28]. doi:10.1371/journal.pcbi.1002011.g009

Before performing MD, the system was minimized using 2000 steps of steepest descent minimization followed by 2000 steps of conjugate gradient minimization. The water box was equilibrated by increasing the temperature from 0 to 300 K while maintaining the solute under constraints, and then further maintaining the simulation at 300 K for 40 ps. In the production phase, the restraints in the binding site were removed and 500 ps MD was performed with a 2 fs time step. The non-bonded cutoff was set to 12 A and SHAKE was applied for all the bonds involving hydrogen atoms. For strategy (i), water molecules and residues within 6 A of the ligand were flexible; for strategy (ii), water molecules and residues within 4.5 A of the ligand were flexible; for strategy (iii), water molecules and residues within 6 A of the mutated residues were flexible; for strategy (iv), all the atoms of the ligand and water molecules were flexible; and for strategy (v), only the water molecules were flexible. For every strategy, 50 frames

were extracted every 2 ps in the interval comprised between 400500 ps of the 500 ps MD simulation.

Additional MD simulations of some mutants

To provide qualitative explanations to the effect of the mutations of ImI-D5N, ImI-R7Q, a7-Q117A, a7-R186V, a7-N111S, a7-S113A, a7-P120A and a7-G153S, additional MD were performed for at least 500 ps. In those MD, a similar water cap was used, as described previously, and residues within 6 A of the mutated side chain were flexible.

Computation of mutational energy

The values of the binding free energy (AG binding) for each mutant were calculated based on the following equation:

AGbinding Gcomplex Gligand Greceptor

The free energy can be decomposed into three components:

G= <Gsolute> + <Gepol> + < GsA >

where G solute is the solute Gibbs free energy, G epol represents the polar contribution to the solvation energy and G SA represents non-polar contribution to the solvation energy. Polar contribution to the solvation energy is determined by solving the Poisson-Boltzmann Equation using the PB module implemented in AMBER10 [61], or the GB approach implemented in AMBER10 [62]. The non-polar contribution to the solvation energy is calculated using:

Gsa = y * SASA + a

where solvent-accessible surface area (SASA) was determined using the Molsurf [63] algorithm with a probe radius of 1.4 A. The surface tension y and constant parameter a in equation (3) were taken to their default values 0.0072 kcal/mol 1 A2 and, 0 kcal/ mol_1 A2 respectively. The effect of residue mutation on the binding energy was computed using:

AAGbinding = DGbinding(mut) - AGbinding (WT)

where AAG binding was defined as mutational energy that is the binding energy difference between the wild-type ligand (AG binding (WT)) and its mutants (AG binding (mut)).

The entropy contribution was estimated using normal-mode analysis, which employed the atomic fluctuation matrix produced from a normal mode calculation [64]. Alternatively, for some calculations, we made the approximation that the wild-type and mutated complexes have similar entropies. Using this approximation:

< GSolute > (mut) - < GSolute > (WT) = < Emm > (mut) - < Emm > (WT)

where E MM is the molecular mechanical energies of the proteins as given by the molecular mechanics potential. This equation was used to compute the difference of internal Gibbs free energy for the complex, the ligand and the receptor. The binding free energy of the mutants was calculated by solving equations (4) and (5) using the script, which is part of the AMBER10 distribution. The Poisson-Boltzmann equations were solved using

internal dielectric and external dielectric constants set to 2.0 and 80.0, respectively, a probe radius of 1.4 A, a grid spacing set to 0.5 A and ionic strength set to 0.15 M/L. For the GB algorithm, the salt concentration was set to 0.15 M/L.

Supporting Information

Figure S1 Rmsd, Rmsf and distance plots of the apo-state model and ImI/a7-nAChR complex over the 10 ns molecular dynamics simulations. Three a7-nAChR apo-state and two a7-nAChR/ImI simulations were performed in total. In the first row, b strand a carbon root-mean-square deviations (RMSD) of each of the subunits over the molecular dynamics simulations to the starting frame. In the second row, a carbon root-mean-square fluctuation (RMSF) of each subunit over the 10 ns molecular dynamics simulation ensemble. In the third row, fluctuation of the distance between the sulfur atom of a7-C190 side chain and the a carbon of a7-Y32. This distance characterizes the closure of the C-loop. (TIFF)

Figure S2 Correlation between the experimentally derived mutational energies and calculated mutational energies of ImI and receptor mutants. Mutational energies were computed using


1. Karlin A (2002) Emerging structure of the nicotinic acetylcholine receptors. Nat RevNeurosci 3: 102-114. doi:10.1038/nrn731.

2. Connolly CN, Wafford KA (2004) The Cys-loop superfamily of ligand-gated ion channels: the impact of receptor structure on function. Biochem Soc Trans 32: 529-534. doi:10.1042/BST0320529.

3. Arneric SP, Holladay M, Williams M (2007) Neuronal nicotinic receptors: a perspective on two decades of drug discovery research. Biochem Pharmacol 74: 1092-1101. doi:10.1016/j.bcp.2007.06.033.

4. Levin ED, Rezvani AH (2007) Nicotinic interactions with antipsychotic drugs, models of schizophrenia and impacts on cognitive function. Biochem Pharmacol 74: 1182-1191. doi:10.1016/j.bcp.2007.07.019.

5. Romanelli MN, Gratteri P, Guandalini L, Martini E, Bonaccini C, et al. (2007) Central nicotinic receptors: structure, function, ligands, and therapeutic potential. ChemMedChem 2: 746-767. doi:10.1002/cmdc.200600207.

6. Taly A, Corringer P, Guedin D, Lestage P, Changeux J (2009) Nicotinic receptors: allosteric transitions and therapeutic targets in the nervous system. Nat Rev Drug Discov 8: 733-750. doi:10.1038/nrd2927.

7. Unwin N (2005) Refined structure of the nicotinic acetylcholine receptor at 4A resolution. J Mol Biol 346: 967-989. doi:10.1016/j.jmb.2004.12.031.

8. Gotti C, Moretti M, Gaimarri A, Zanardi A, Clementi F, et al. (2007) Heterogeneity and complexity of native brain nicotinic receptors. Biochem Pharmacol 74: 1102-1111. doi:10.1016/j.bcp.2007.05.023.

9. Grady SR, Moretti M, Zoli M, Marks MJ, Zanardi A, et al. (2009) Rodent habenulo-interpeduncular pathway expresses a large variety of uncommon nAChR subtypes, but only the a3b4* and a3P3P4* subtypes mediate acetylcholine release. J Neurosci 29: 2272-2282. doi:10.1523/jneurosci.5121-08.2009.

10. Gotti C, Zoli M, Clementi F (2006) Brain nicotinic acetylcholine receptors: native subtypes and their relevance. Trends Pharmacol Sci 27: 482-491. doi:10.1016/

11. Olivera BM, Gray WR, Zeikus R, McIntoshJM, VargaJ, et al. (1985) Peptide neurotoxins from fish-hunting cone snails. Science 230: 1338-1343. doi:10.1126/science.4071055.

12. Olivera BM, Rivier J, Clark C, Ramilo CA, Corpuz GP, et al. (1990) Diversity of Conus neuropeptides. Science 249: 257-263. doi:10.1126/science.2165278.

13. Terlau H, Olivera BM (2004) Conus venoms: a rich source of novel ion channel-targeted peptides. Physiol Rev 84: 41-68. doi:10.1152/physrev. 00020.2003.

14. Buczek O, Bulaj G, Olivera B (2005) Conotoxins and the posttranslational modification of secreted gene products. Cell Mol Life Sci 62: 3067-3079. doi:10.1007/s00018-005-5283-0.

15. Kaas Q, Westermann J, Craik DJ (2010) Conopeptide characterization and classifications: an analysis using ConoServer. Toxicon 55: 1491-1509. doi:10.1016/j.toxicon.2010.03.002.

16. Craik DJ, Adams DJ (2007) Chemical modification of conotoxins to improve stability and activity. ACS Chem Biol 2: 457-468. doi:10.1021/cb700091j.

17. Halai R, Craik DJ (2009) Conotoxins: natural product drug leads. Nat Prod Rep 26: 526-536. doi:10.1039/b819311h.

18. McIntosh JM, Santos AD, Olivera BM (1999) Conus peptides targeted to specific nicotinic acetylcholine receptor subtypes. Annu Rev Biochem 68: 59-88. doi:10.1146/annurev.biochem.68.1.59.

either molecular mechanics generalized Born (GB) surface area (MM-GB/SA) or molecular mechanics Poisson-Boltzmann (PB) surface area (MM-PB/SA) energy functions at 298 K. The mutated models were refined using MBA with either distance dependent dielectric constant minimization (DDDCM) or explicit water minimization (EWM). (TIFF)


We gratefully acknowledge access to the facilities of the Australian Research Council (ARC) Special Research Centre for Functional and Applied Genomics. We thank Professor Adrian. E. Roitberg and Dr. Jason M. Swails for their latest version of script and excellent technical assistance.

Author Contributions

Conceived and designed the experiments: DJC QK. Performed the experiments: RY QK. Analyzed the data: RY DJC QK. Contributed reagents/materials/analysis tools: DJC QK. Wrote the paper: RY DJC QK.

19. Dutton JL, Craik DJ (2001) a-Conotoxins: nicotinic acetylcholine receptor antagonists as pharmacological tools and potential drug leads. Curr Med Chem 8: 327-344.

20. McIntoshJM, Yoshikami D, Mahe E, Nielsen DB, Rivier JE, et al. (1994) A nicotinic acetylcholine receptor ligand of unique specificity, a-conotoxin ImI. J Biol Chem 269: 16733-16739.

21. Johnson DS, Martinez J, Elgoyhen AB, Heinemann SF, McIntoshJM (1995) a-Conotoxin ImI exhibits subtype-specific nicotinic acetylcholine receptor blockade: preferential inhibition of homomeric a 7 and a 9 receptors. Mol Pharmacol 48: 194-199. doi: 10.1021/bi048918g.

22. Ellison M, Gao F, Wang H, Sine SM, McIntoshJM, et al. (2004) a-conotoxins ImI and ImII target distinct regions of the human a7 nicotinic acetylcholine receptor and distinguish human nicotinic receptor subtypes. Biochemistry 43: 16019-16026. doi:10.1021/bi048918g.

23. Lamthanh H, Jegou-Matheron C, Servent D, Menez A, Lancelin JM (1999) Minimal conformation of the a-conotoxin ImI for the a7 neuronal nicotinic acetylcholine receptor recognition: correlated CD, NMR and binding studies. FEBS Lett 454: 293-298. doi:10.1016/S0014-5793(99)00831-5.

24. Maslennikov IV, Shenkarev ZO, Zhmak MN, Ivanov VT, Methfessel C, et al. (1999) NMR spatial structure of a-conotoxin ImI reveals a common scaffold in snail and snake toxins recognizing neuronal nicotinic acetylcholine receptors. FEBS Lett 444: 275-280. doi:10.1016/S0014-5793(99)00069-1.

25. Gehrmann J, Daly NL, Alewood PF, Craik DJ (1999) Solution structure of a-conotoxin ImI by 1H nuclear magnetic resonance. J Med Chem 42: 2364-2372. doi:10.1021/jm990114p.

26. Quiram PA, Sine SM (1998) Structural elements in a-conotoxin ImI essential for binding to neuronal a7 receptors. J Biol Chem 273: 11007-11011. doi:10.1074/ jbc.273.18.11007.

27. Servent D, Thanh HL, Antil S, Bertrand D, Corringer PJ, et al. (1998) Functional determinants by which snake and cone snail toxins block the a7 neuronal nicotinic acetylcholine receptors. J Physiol Paris 92: 107-111. doi:10.1016/S0928-4257(98)80146-0.

28. Quiram PA, JonesJJ, Sine SM (1999) Pairwise interactions between neuronal a7 acetylcholine receptors and a-conotoxin ImI. J Biol Chem 274: 19517-19524. doi:10.1074/jbc.274.28.19517.

29. Rogers JP, Luginbuhl P, Pemberton K, Harty P, Wemmer DE, et al. (2000) Structure-activity relationships in a peptidic a7 nicotinic acetylcholine receptor antagonist. J Mol Biol 304: 911-926. doi:10.1006/jmbi.2000.4247.

30. Armishaw C, Jensen AA, Balle T, Clark RJ, Harpsoe K, et al. (2009) Rational design of a-conotoxin analogues targeting a7 nicotinic acetylcholine receptors: improved antagonistic activity by incorporation of proline derivatives. J Biol Chem 284: 9498-9512. doi:10.1074/jbc.M806136200.

31. Armishaw CJ, Singh N, Medina-Franco JL, Clark RJ, Scott KCM, et al. (2010) A synthetic combinatorial strategy for developing a-conotoxin analogs as potent a7 nicotinic acetylcholine receptor antagonists. J Biol Chem 285: 1809-1821. doi:10.1074/jbc.M109.071183.

32. Dutertre S, Nicke A, Tyndall JDA, Lewis RJ (2004) Determination of a-conotoxin binding modes on neuronal nicotinic acetylcholine receptors. J Mol Recognit 17: 339-347. doi:10.1002/jmr.683.

33. Hansen SB, Sulzenbacher G, Huxford T, Marchot P, Taylor P, et al. (2005) Structures of Aplysia AChBP complexes with nicotinic agonists and antagonists

reveal distinctive binding interfaces and conformations. EMBO J 24: 3635-3646. doi:10.1038/sj.emboj.7600828.

34. Ulens C, Hogg RC, Celie PH, Bertrand D, Tsetlin V, et al. (2006) Structural determinants of selective a-conotoxin binding to a nicotinic acetylcholine receptor homolog AChBP. Proc Natl Acad Sci U S A 103: 3615-3620. doi:10.1073/pnas.0507889103.

35. Dellisanti CD, Yao Y, Stroud JC, Wang Z, Chen L (2007) Crystal structure of the extracellular domain of nAChR a1 bound to a-bungarotoxin at 1.94 A resolution. NatNeurosci 10: 953-962. doi:10.1038/nn1942.

36. Halai R, Clark RJ, Nevin ST, Jensen JE, Adams DJ, et al. (2009) Scanning mutagenesis of a-conotoxin Vc1.1 reveals residues crucial for activity at the a9a10 nicotinic acetylcholine receptor. J Biol Chem 284: 20275-20284. doi:10.1074/jbc.M109.015339.

37. Henchman RH, Wang H, Sine SM, Taylor P, McCammon JA (2003) Asymmetric structural motions of the homomeric a7 nicotinic receptor ligand binding domain revealed by molecular dynamics simulation. Biophys J 85: 3007-3018. doi:10.1016/S0006-3495(03)74720-1.

38. Cheng X, Wang H, Grant B, Sine SM, McCammon JA (2006) Targeted molecular dynamics study of C-loop closure and channel gating in nicotinic receptors. PLoS Comput Biol 2: e134. doi:10.1371/journal.pcbi.0020134.

39. Liu X, Xu Y, Li H, Wang X, Jiang H, et al. (2008) Mechanics of channel gating of the nicotinic acetylcholine receptor. PLoS Comput Biol 4: e19. doi:10.1371/ journal.pcbi.0040019.

40. Yi M, Tjong H, Zhou H (2008) Spontaneous conformational change and toxin binding in a7 acetylcholine receptor: insight into channel activation and inhibition. Proc Natl Acad Sci U S A 105: 8280-8285. doi:10.1073/ pnas.0710530105.

41. Yakel JL (2010) Gating of nicotinic ACh receptors: latest insights into ligand binding and function. J Physiol 588: 597-602. doi:10.1113/jphysiol. 2009.182691.

42. Gleitsman KR, Shanata JA, Frazier SJ, Lester HA, Dougherty DA (2009) Longrange coupling in an allosteric receptor revealed by mutant cycle qnalysis. Biophys J 96: 3168-3178. doi:10.1016/j.bpj.2008.12.3949.

43. Pons C, Grosdidier S, Solernou A, Perez-Cano L, Fernandez-Recio J (2010) Present and future challenges and limitations in protein-protein docking. Proteins 78: 95-108. doi:10.1002/prot.22564.

44. Kang TS, Radic Z, Talley TT, Jois SDS, Taylor P, et al. (2007) Protein folding determinants: structural features determining alternative disulfide pairing in a-and K/l-conotoxins. Biochem 46: 3338-3355. doi:10.1021/bi061969o.

45. Tsui V, Case DA (2000) Theory and applications of the generalized born solvation model in macromolecular simulations. Biopolymers 56: 275-291. doi:10.1002/1097-0282(2000)56:4<275::AID-BIP10024>3.0.C0; 2-E.

46. Rastelli G, Del Rio A, Degliesposti G, Sgobba M (2010) Fast and accurate predictions ofbinding free energies using MM-PBSA and MM-GBSA. J Comput Chem 31: 797-810. doi:10.1002/jcc.21372.

47. Zhang X, Li X, Wang R (2009) Interpretation of the binding affinities of PTP1B inhibitors with the MM-GB/SA method and the X-score scoring function. J Chem Inf Model 49: 1033-1048. doi:10.1021/ci8004429.

48. Hou T, Wang J, Li Y, Wang W (2010) Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses

generated from docking. J Comput Chem. (ahead of print);doi:10.1002/ jcc.21666.

49. Liu H, Yao X (2010) Molecular basis of the interaction for an essential subunit PA—PB1 in influenza virus RNA polymerase: insights from molecular dynamics simulation and free energy calculation. Mol Pharm 7: 75-85. doi:10.1021/ mp900131p.

50. London N, Movshovitz-Attias D, Schueler-Furman O (2010) The structural basis of peptide-protein binding strategies. Structure 18: 188-199. doi:10.1016/ j.str.2009.11.012.

51. Le Novere N, Grutter T, ChangeuxJ (2002) Models of the extracellular domain of the nicotinic receptors and of agonist- and Ca2+-binding sites. Proc Natl Acad Sci U S A 99: 3210-3215. doi:10.1073/pnas.042699699.

52. Costa V, Nistri A, Cavalli A, Carloni P (2003) A structural model of agonist binding to the a3^4 neuronal nicotinic receptor. Br J Pharmacol 140: 921-931. doi:10.1038/sj.bjp.0705498.

53. Perez EG, Cassels BK, Zapata-Torres G (2009) Molecular modeling of the a9a10 nicotinic acetylcholine receptor subtype. Bioorg Med Chem Lett 19: 251-254. doi:10.1016/j.bmcl.2008.10.094.

54. Edgar R (2004) MUSCLE: a multiple sequence alignment method with reduced time and space complexity. BMC Bioinformatics 5: 113. doi:10.1186/1471-2105-5-113.

55. Sali A, Blundell TL (1993) Comparative protein modelling by satisfaction of spatial restraints. J Mol Biol 234: 779-815. doi:10.1006/jmbi.1993.1626.

56. Morris AL, MacArthur MW, Hutchinson EG, Thornton JM (1992) Stereo-chemical quality of protein structure coordinates. Proteins 12: 345-364. doi:10.1002/prot.340120407.

57. Chen VB, Arendall WB, Headd JJ, Keedy DA, Immormino RM, et al. (2010) MolProbity: all-atom structure validation for macromolecular crystallography. Acta Crystallogr D Biol Crystallogr 66: 12-21. doi: 10.1107/ S0907444909042073.

58. Chowdhuri S, Tan M, Ichiye T (2006) Dynamical properties of the soft sticky dipole-quadrupole-octupole water model: a molecular dynamics study. J Chem Phys 125: 144513. doi:10.1063/1.2357117.

59. Sagui C, Darden TA (1999) Molecular dynamics simulations of biomolecules: long-range electrostatic effects. Annu Rev Biophys Biomol Struct 28: 155-179. doi:10.1146/annurev.biophys.28.1.155.

60. Hess B, Bekker H, Berendsen HJC, Fraaije JGEM (1997) LINCS: A linear constraint solver for molecular simulations. J Comput Chem 18: 1463-1472. doi:10.1002/(SICI)1096-987X(199709)18:12<1463::AID-JCC4>3.0.C0; 2-H.

61. Case DA, Darden TA, Cheatham III TE, Simmerling CL, Wang J, et al. (2008) AMBER 10. San Francisco: University of California.

62. Onufriev A, Bashford D, Case DA (2004) Exploring protein native states and large-scale conformational changes with a modified generalized born model. Proteins 55: 383-394. doi:10.1002/prot.20033.

63. Connolly ML (1983) Solvent-accessible surfaces of proteins and nucleic acids. Science 221: 709-713. doi:10.1126/science.6879170.

64. Wang J, Morin P, Wang W, Kollman PA (2001) Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA. J Am Chem Soc 123: 5221-5230. doi:10.1021/ja003834q.

Copyright of PLoS Computational Biology is the property of Public Library of Science and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.