Free Energy Simulations of Binding of HsTx1 Toxin ... - ACS Publications

Jan 7, 2014 - The scorpion toxin HsTx1 is an equally potent blocker of Kv1.3, which is ...... Toxin from the Caribbean Sea Anemone Stichodactyla helia...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Free Energy Simulations of Binding of HsTx1 Toxin to Kv1 Potassium Channels: the Basis of Kv1.3/Kv1.1 Selectivity M. Harunur Rashid and Serdar Kuyucak* School of Physics, University of Sydney, Sydney, New South Wales 2006, Australia ABSTRACT: The voltage-gated potassium channel Kv1.3 is an established target for treatment of autoimmune diseases. Hence, there are intense efforts to develop immunosuppressant drugs from Kv1.3 blockers. ShK toxin from sea anemone is the most advanced peptide in this regard, but its lack of selectivity for Kv1.3 over Kv1.1 is an ongoing concern. The scorpion toxin HsTx1 is an equally potent blocker of Kv1.3, which is also selective for Kv1.3. It is of interest to understand the molecular basis of this selectivity as the lessons learned may suggest new avenues for enhancing the selectivity of other Kv1.3 blockers. Here we construct accurate models of Kv1.x−HsTx1 complexes using docking and molecular dynamics simulations. For each complex, the binding free energy of HsTx1 is determined from the potential of mean force calculations. Good agreement is found between the computed and experimental binding free energies, which increases confidence in the complex models. Comparison of the binding modes of HsTx1 with Kv1.1 and Kv1.3 reveals that the lower affinity of HsTx1 for Kv1.1 is due to its inability to come close to the pore domain, which prevents the pore inserting lysine from making proper contacts with the tyrosine carbonyls in the selectivity filter.



INTRODUCTION Voltage-gated potassium (Kv) channels are involved in many cellular functions.1 Dysfunction or proliferation of Kv channels interferes with normal cell function, resulting in many diseases ranging from neurological and cardiovascular to more complex metabolic and autoimmune diseases.2 Such diseases may be treated by blocking or modulating the particular Kv channel associated with a disease.3,4 In the case of autoimmune diseases, such as multiple sclerosis, rheumatoid arthritis, and type 1 diabetes, overexpression of Kv1.3 channels has been implicated in their onset.5−7 Thus, potent and selective blockers of Kv1.3 channels could be used in treatment of autoimmune diseases.8−10 ShK toxin from sea anemone11 has a very high affinity for Kv1.3 (IC50 = 11 pM)12 and therefore has attracted the most attention in this regard. However, ShK has a similarly high affinity for Kv1.1 (IC50 = 16 pM),12 and to prevent side effects, it is essential to find analogues of ShK that are selective for Kv1.3 over Kv1.1. This program has been vigorously pursued over the past decade, resulting in some potential drug leads, e.g., ShK-Dap,12 ShK-186,13 and ShK-192.14 Unfortunately, use of nonprotein components in these analogues has been found to be a limiting factor in their further development as therapeutic agents, and it is desirable to find new Kv1.3selective analogues of ShK containing only natural amino acids. The scorpion toxin HsTx115 also has a very high affinity for Kv1.3 (IC50 = 11 pM),16 similar to ShK, but it also exhibits nearly thousand-fold selectivity for Kv1.3 over Kv1.1. Thus, HsTx1 avoids the selectivity problem that has been a major obstacle for ShK, and could be a candidate for treatment of autoimmune diseases. Unfortunately, blocking properties of HsTx1 for Kv1 channels were discovered much later than ShK. As a result, its physiological and pharmacological properties are © 2014 American Chemical Society

not as well described as ShK to consider it for drug development at present. Nevertheless, it is important to understand how nature has achieved such an exquisite Kv1.3/ Kv1.1 selectivity. Insights gathered from the binding modes of Kv1.x−HsTx1 complexes may be used in development of ShK or other peptide analogues which have enhanced selectivity for Kv1.3. Toxin binding to ion channels has been pursued in many computational studies in the past decade (see refs 17 and 18 for reviews). The initial focus was on determination of complex structures using a variety of methods from docking19,20 to Brownian dynamics21,22 and molecular dynamics23,24 (MD) simulations. Since then, docking combined with MD simulations for refinement has been established as an optimal method for complex structure determination,25 and used in many toxin binding studies.26−43 Validation of a predicted complex structure is an important issue, which has rarely been addressed in earlier computational work due to lack of data and difficulty of free energy calculations. Mutagenesis data are ideal for validation purposes, but such data are available in very few cases. Alternatively, one can use binding free energieswhich are widely availablebut their accurate computation for peptides has been a problem.44−46 This problem has been overcome in a recent computational study of binding of charybdotoxin to a KcsA potassium channel surrogate.31,32 The complex structure was determined from NMR,47 so there was no uncertainty with regard to the structure for this complex. The potential of mean force (PMF) of the toxin was Received: November 7, 2013 Revised: January 2, 2014 Published: January 7, 2014 707

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

Table 1. Pore Domain Residues That Differ among the Kv1.x Channelsa mKv1.1 rKv1.2 mKv1.3 a

E350 D352 D375

E351 E353 D376

A352 R354 P377

E353* D355* S378*

H355 Q357 G380

S357 P359 N382

S369 S371 T394

Y379* V381* H404*

V381 T383 V406

Residues that play a role in selectivity are indicated by an asterisk.

1, it consists of an α helix sitting on a β sheet with two strands. The structure of HsTx1 is stabilized by four disulfide bridges

constructed using umbrella sampling MD simulations, from which the binding free energy was determined within chemical accuracy.32 Following this development, binding free energies have been used for validation purposes in many studies of toxin binding to ion channels.33−41,43 We note, in particular, the study of binding of ShK to Kv1 channels,40 which forms the basis of this work. Here we present a computational study of the binding of HsTx1 to the Kv1.1, Kv1.2, and Kv1.3 channels that will help to understand the molecular basis of its selectivity for Kv1.3. Complex structures are determined from docking followed by refinement with MD simulations. Because there are no mutagenesis data for the Kv1.x−HsTx1 complexes, we validate the models using the binding free energy data. For each complex, we construct the PMF for dissociation of HsTx1 from umbrella sampling MD simulations. The standard binding free energies are determined from the integral of the PMFs and compared to the experimental values. It is hoped that the present work will rekindle experimental interest in HsTx1 as a high-affinity Kv1.3 blocker, and spur further activity for its development as a therapeutic agent for treatment of autoimmune diseases.



Figure 1. NMR structure of HsTx1 oriented with the pore inserting lysine (K23) pointing downward. Four disulfide bridges between C3− C24, C9−C29, C13−C31, and C19−C34 are indicated explicitly as well as the hydrogen bonds between the residues A1−C29 and C3− R27.

METHODS

between C3−C24, C9−C29, C13−C31, and C19−C34 as well as hydrogen bonds between the N-terminal and the β sheet. The first three disulfide bridges are similar to other scorpion toxins with three disulfide bridges and are responsible for the common α/β scaffold motif.53 The fourth disulfide bridge has been shown to keep the β sheet structure in a twisted configuration.54 Its removal changes the β sheet to a nontwisted configuration and has a substantial effect on its affinity for Kv1.3, reducing it by 50-fold.54 Thus, the β sheet structure of HsTx1 is critical for its high-affinity binding to Kv1.3. This can be readily surmised from Figure 1, where HsTx1 is presented with the pore inserting lysine (K23) pointing downward. The binding interface with the channel is seen to be formed mainly by the β sheet. The initial configurations of Kv1.x−HsTx1 complexes are generated using the protein−protein docking program HADDOCK.55,56 The pore inserting lysine is used as a restraint, and an ensemble of 10 NMR structures of HsTx1 is employed in the docking calculations. A unique consensus model is found for each complex, which is then embedded in a lipid bilayer and solvated. The resulting systems are relaxed in MD simulations following the protocols established in Kv1− ShK simulations.40 The equilibrated systems are run for a further 20 ns, and the trajectory data collected are used in checking the stability of the complexes and analysis of the binding mode. In all three cases, the complex structures have been found to be very stable. As a test of the docking program, one can compare the predicted distance between the centers of mass (COM) of the protein and the ligand with the average distance obtained from the MD simulations. In all three complexes, the HADDOCK prediction for the channel-COM to toxin-COM distance is found to be within 1 Å of the MD average, which is about 27 Å.

Modeling of Kv1 Channels. The crystal structure of the Kv1.2 channel is available (PDB ID 2R9R),48 and those of Kv1.1 and Kv1.3 have been constructed from the Kv1.2 structure as detailed in our previous work on binding of ShK to Kv1 channels.40 Because the homology among these Kv1 channels is over 90% and there is one-to-one correspondence between the pore domain residues (see Table 1), construction of the Kv1.1 and Kv1.3 models is straightforward. The only point to be noted is that the side chain of H404 in Kv1.3 which replaces V381 in Kv1.2forms an H-bond with the side chain of D402 in the neighboring monomer. Correct modeling of the H404−D402 cross-linking in the pore mouth was shown to be essential40 for reproducing the alanine scanning mutagenesis data for the Kv1.3−ShK complex.49 The simulation system for each channel is constructed using the VMD software50 as previously described.40 The channel protein is embedded in a bilayer consisting of 125 POPC molecules and solvated with a 0.1 M KCl solution. Each system is relaxed and equilibrated following the protocols developed in previous MD simulations of potassium channels.51 Stability and integrity of the channel models have been well established during the binding studies of ShK to Kv1 channels.40 Modeling of Kv1−HsTx1 Complexes. The structure of HsTx1 was determined from solution NMR (PDB ID 1QUZ).52 It is a 34-residue peptide with the sequence ASCRTPKDCA DPCRKETGCP YGKCMNRKCK CNRC. The natural form of HsTx1 has an amidated C-terminal, which has 5-fold more affinity for Kv1.3 compared to the synthetic carboxylated form.15 Amidated HsTx1 contains four Arg, five Lys, two Asp, and one Glu residues with a total charge of +7e, which is the same as in ShK. However, the structure of HsTx1 is quite different from that of ShK. As shown in Figure 708

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

MD Simulations and Umbrella Sampling. The NAMD software57 with the CHARMM36 force field is used to carry out the MD simulations.58 An NpT ensemble is used with periodic boundary conditions. The temperature and pressure are maintained at 300 K and 1 atm, respectively, using Langevin coupling with damping coefficients of 5 ps−1. The particle-mesh Ewald algorithm is employed to calculate the long distance electrostatic interactions. Lennard-Jones interactions are switched off within a distance of 10−12 Å. The list of nonbonded interactions are truncated at 13.5 Å. A time step of 2 fs is used in MD simulations, and the trajectory data is written at 1 ps intervals. The reaction coordinate (Kv1.x-COM to HsTx1-COM distance along the z axis) is recorded at every time step in umbrella sampling simulations. Umbrella sampling simulations are performed using the protocols from our earlier work on toxin binding to potassium channels.31,32,40 Umbrella windows are generated at 0.5 Å intervals along the channel axis starting from the binding site and moving the toxin in both directions. A harmonic force with k = 30 kcal/mol/Å2 is applied to the COM of HsTx1 backbone atoms, which is pulled at a speed of 5 Å/ns for 0.1 ns in steered MD simulations. After each steering simulation, HsTx1 is equilibrated for 0.3 ns with the same force constant to relax the toxin. PMFs are extended up to a distance of 15 Å from the binding pocket, which is found to be sufficient to obtain a flat region of PMF signaling the bulk conditions. We require 5% overlap between the neighboring window samples,32 and include an extra window in between if this condition is not satisfied. Each window is simulated using a force constant of 30 kcal/mol/Å2, and the reaction coordinates of HsTx1 are saved for PMF calculations. The data from all the windows are combined and unbiased using the weighted histogram analysis method59 to construct the PMF of HsTx1. Simulations are continued until we have established the convergence of the PMFs from block data analysis. The binding constant of HsTx1 is determined from the 3-D integral of the PMF, W(r), as Keq =

∫site e−W (r)/k T d3r B

(1)

where it is assumed that W = 0 in the bulk. Because computation of a 3-D PMF is prohibitively expensive, a 1-D approximation is invoked, and the binding constant is determined by integrating the PMF, W(z), along the z axis60 Keq = πR2

∫z

z2

e−W (z)/ kBT dz

1

Figure 2. (A) Distribution of the transverse (x and y) COM coordinates of HsTx1 in the binding pocket of Kv1.3 relative to the average COM position. Both distributions are fitted using a Gaussian with σ = 0.48 Å corresponding to an R value of 0.68 Å. (B) Deviation of the average transverse COM position of HsTx1 from the initial one in the binding pocket as the toxin is pulled out. (C) Changes in the effective radius R as the toxin is pulled out from the binding site to the bulk.

(2)

Here the integration limits are chosen such that z1 is in the binding site and z2 is in the bulk where W vanishes. The factor πR2 measures the average cross-sectional area of the binding pocket as explored by the COM of HsTx1. The effective radius R is determined from the transverse fluctuations of the COM of HsTx1 in the binding pocket. Because there are different methods of choosing R in the literature, we briefly justify our choice here. Inspection of the transverse (i.e., x and y) coordinates of the toxin-COM in the binding pocket of Kv1.3 (Figure 2A) shows that they have a Gaussian distribution given by ρ(x , y) = ρ0 e−(x

2

2

+ y )/2σ

ρ(x , y) = ρ0 e−W (x , y)/ kBT

(4)

one can see that the transverse part of the PMF can be represented by a harmonic potential W (x , y ) =

1 k(x 2 + y 2 ) 2

(5)

2

with k = kBT/σ . Assuming that the total PMF can be written as W(r) = W(x, y) + W(z), the x and y integrals in eq 1 separate and can be evaluated as

2

(3)

with equal widths for the x and y distributions, i.e., σx = σy = σ. Comparing with the Boltzmann distribution

2

2

∫ e−k(x +y )/2k T dx dy = 2πkBT /k = 2πσ 2 B

709

(6)

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

Figure 3. Snapshots of the Kv1.x−HsTx1 complexes showing the strongly interacting residues (top panels, Kv1.1; middle, Kv1.2; bottom, Kv1.3). To give the full picture, two views of the cross sections of the complex depicting the monomers C and A (left panel) and B and D (right panel) are presented. Only the residues involved in the binding are indicated explicitly.

Using the definition of R, R2 = σx2 + σy2 = 2σ2, in eq 6, and substituting this result in eq 1, one obtains eq 2 for the binding constant. For the validity of this approximation, the toxin’s average transverse position and fluctuations around it should not deviate much from those in the binding pocket while the toxin is pulled out. This is shown in Figure 2B and C, which are obtained from the umbrella sampling simulation of Kv1.3− HsTx1. The average axial position of HsTx1 remains well within 1 Å of its position in the binding pocket in all umbrella windows (Figure 2B). The transverse fluctuations of the COM of HsTx1 have an approximate Gaussian distribution in all windows. The R values determined from the widths of the Gaussian distributions are seen to rise from 0.68 Å in the

binding site to slightly above 1 Å as the toxin is pulled out (Figure 2C). From the 2-D diffusion equation, R is expected to grow like R = (4Dt)1/2 in the bulk, where D is the diffusion coefficient of the toxin. Due to the large inertia of the toxin, D is quite small; hence, in ns-MD simulations, R changes little from its value in the binding site. We also note that R remains nearly constant in the vicinity of the binding pocket, where the major contribution to the integral in eq 2 comes from. In fact, performing the integral for the binding constant using a variable R as in Figure 2C gives essentially the same result as using a fixed R in eq 2. The values of R obtained from restraint-free MD simulations of the complexes are 0.71, 0.72, and 0.68 Å, respectively, for Kv1.1−HsTx1, Kv1.2−HsTx1, and Kv1.3− HsTx1. 710

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

The standard binding free energy of HsTx1 is determined from the binding constant using G b = −kBT ln(KeqC0)

Table 2. List of the Strongly Interacting Residues in the Kv1.x−HsTx1 Complexesa

(7)

where C0 is the standard concentration of 1 M (i.e., 1/C0 = 1661 Å3).



RESULTS AND DISCUSSION Characterization of the Kv1.x−HsTx1 Complexes. Structures of the Kv1.x−HsTx1 complexes are shown in Figure 3. The only common motif between the three binding modes is the K23 side chain making contact with the filter carbonyls. In all other aspects, the Kv1.3−HsTx1 binding mode is quite different from the other two. For example, HsTx1 has no contacts with the residues in the filter periphery in Kv1.1 and Kv1.2. It has only contacts with the turret residues E353 in Kv1.1 and D355 in Kv1.2. This situation is reversed in Kv1.3. The residue corresponding to E353/D355 is S378 in Kv1.3 (see Table 1), which is involved in a hydrogen bond. All the other contacts in Kv1.3 are with the residues in the filter periphery. We note that such a clear distinction between the binding modes is not observed in the Kv1.x−ShK complexes (see Figure 3 in ref 40), where both the turret and filter periphery residues contribute to the binding modes. A more careful comparison of the Kv1.1−HsTx1 and Kv1.1− ShK complexes shows that HsTx1 is not as close to the pore domain of Kv1.1 as ShK. In particular, K23 is unable to fully insert into the filter and make contacts with three Tyr carbonyls. Instead, K23 makes contacts with one Tyr carbonyl and two Gly carbonyls on one side of the pore (Figure 3, top). Interaction of the pore inserting Lys with the Tyr carbonyls in the filter has been a uniform feature of the Kv1.x−ShK complexes. Its loss in the Kv1.1−HsTx1 complex, together with lack of contacts in the filter periphery, is expected to have a detrimental effect on the binding free energy of HsTx1. A more quantitative characterization of the binding modes is provided by the average atom−atom distances for strongly interacting residues (Table 2). Here we also list the HADDOCK results to see how close the initial docking pose is to the final MD complex. Overall, HADDOCK does a good job in predicting the initial structure, especially in Kv1.3− HsTx1. A few problems are observed in the interactions of the Arg side chains of HsTx1 with the turret residues, which fail to make contact. These may be improved by including further snapshots of HsTx1 from MD simulations in the ensemble docking. A further problem occurs in Kv1.1−HsTx1, where HADDOCK predicts full contact of the K23 side chain with three Tyr residues in the filter, which is not feasible according to the MD simulation results. Again, it is instructive to compare the results in Table 2 with similar ones for ShK (Table 2 in ref 40). In Kv1.1, the main difference is that ShK interacts with two Y379 residues near the filter, which are missing in the binding mode of HsTx1. The other charge contacts are similar, but as pointed out above, the K23 side chain is not properly inserted into the filter, which lessens its impact. Three factors appear to be contributing to the poor docking of HsTx1 to Kv1.1: (i) The binding interface of HsTx1 is formed by a long β sheet which limits formation of contacts with the channel residues. (ii) The Y379 side chains in Kv1.1 are unable to make cross-links (like H404 in Kv1.3) and project out, causing steric barriers. (iii) Formation of three strong charge contacts between the Arg residues of HsTx1 and

HsTx1

Kv1.1

dock

MD aver.

R4-N2 R14-N2 K23-N1 K23-N1 R33-N2 HsTx1

E353-O2(C) E353-O1(B) Y375-O(C) G376-O(BC) E353-O2(A) Kv1.2

8.5 5.1 3.0 6.5 2.7 dock

3.0 ± 0.5 2.9 ± 0.3 2.7 ± 0.3 3.0 ± 0.2 2.7 ± 0.4 MD aver.

R14-N2 K23-N1 R27-N1 HsTx1

D355-O2(C) Y377-O(ACD) D355-O2(D) Kv1.3

5.3 2.7 5.4 dock

2.8 ± 0.3 2.7 ± 0.3 3.1 ± 0.4 MD aver.

T5-Cγ2 K7-N1 Y21-OH K23-N1 N26-Nδ2 R33-N2

M403-Cε(C) S378-O(B) D402-O2(B) Y400-O(ABC) D402-O(D) E373-O2(A)

5.1 6.0 2.7 2.7 3.0 2.7

3.7 2.7 2.7 2.7 2.8 2.7

± ± ± ± ± ±

0.3 0.5 0.2 0.3 0.4 0.3

a

The average atom−atom distances obtained from HADDOCK and MD simulations are given in the third and fourth columns (in units of Å). The N−O distances are shown for the charge interactions and the closest C−C distance for the hydrophobic ones. Bare C, N, and O refer to the backbone atoms, and the subscripted ones refer to the side chain atoms. The monomer identity is given in parentheses. For the filter Y and G residues, the average of the N−O distances is given.

the turret E353 residues prevents the toxin from getting closer to the pore domain. The Kv1.2 channel is quite homologous to Kv1.1; hence, one expects the binding modes of HsTx1 to be quite similar. Nevertheless, two mutations, indicated by stars in Table 1, cause some differences in the binding modes. The shorter side chains of D355 residues mean that it is not possible for two Arg side chains in HsTx1 to form a bridge across two monomers like in Kv1.1−HsTx1 (Figure 3, top). This results in loss of a charge contact but also allows HsTx1 to come closer to the pore domain. The second point is also helped by the less bulky V381 side chains, so that the K23 side chain is able to insert into the filter properly and makes contact with three Tyr carbonyls. Contrasting the binding modes of HsTx1 and ShK to Kv1.2, they are similar apart from HsTx1 having no contacts in the filter periphery. Again, we expect this to result in a lower affinity of HsTx1 to Kv1.2 compared to ShK. Though not as much as in Kv1.1 because K23 is properly inserted into the filter, and a second Arg−Asp charge interaction should reduce the effect of loss of a contact. The situation in the Kv1.3−HsTx1 complex is more complicated. Relative to the Kv1.3−ShK binding mode, two contacts are lost from the filter periphery. However, two of the remaining contacts with the D402 residues are very tight, which compensates for the loss of contacts to some degree. In particular, Y21 makes contact with the side chain of D402(B), which is cross-linked to the H404(A) side chain and hence provides a stable point of binding. In both complexes, there is only one notable charge interaction with the turret residues (E373−R33 in HsTx1 and S378−H19 in ShK). The former is clearly a much stronger interaction, which again helps to compensate for the lost contacts in HsTx1. The smaller number of higher quality contacts in HsTx1 is expected to yield a high affinity binding to Kv1.3 similar to ShK. These qualitative observations on the relative affinities of ShK and 711

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

where the bulk conditions prevailshow that they are similar to the bulk values indicated by the dotted line in Figure 4. We have also visually checked, by aligning the HsTx1 structures from the last umbrella windows with the NMR structure, that there are no deformations in the toxin structure. We have checked the overlaps of the reaction coordinate between the neighboring windows in all PMF calculations and inserted extra windows in between whenever the overlap dropped below 5%. Umbrella sampling simulations are continued until we are assured of the convergence of the PMF. Convergence studies of the PMFs for dissociation of HsTx1 from the Kv1.1, Kv1.2, and Kv1.3 channels are presented in Figure 5. Convergence occurs relatively fast in the case of Kv1.3 but much slower in Kv1.1 and Kv1.2. Comparing the results in Figure 5 with the corresponding ones for ShK (Figure 5 in ref 40), we see that the convergence patterns are similar in Kv1.3 but convergence of the HsTx1 PMFs is delayed in Kv1.1 and Kv1.2 relative to the ShK PMFs. A similar observation can be made for the binding affinities they are similar in Kv1.3, but the affinities of HsTx1 for Kv1.1 and Kv1.2 are reduced compared to those of ShK. These examples indicate that there is a correlation between the binding affinity and convergence of a PMF, at least for the channel−toxin complexes considered here. The PMFs of toxins with high affinity converge quite fast, while those with low affinity take a longer time to converge. The larger range of equilibration times in Figure 5 (from 1 to 4.5 ns) again highlights the importance of performing a proper convergence study before choosing the production data for the final PMF. For example, use of a uniform 1 ns data for equilibration would corrupt the PMFs in Kv1.1 and Kv1.2, resulting in larger well depths (and affinities) than those indicated by the final PMFs in Figure 5. The final PMFs for the Kv1.x−HsTx1 complexes are shown in Figure 6. We will comment on the behavior of the PMFs when we discuss the evolution of the contact distances below. The relative and standard binding free energies determined from the PMFs are listed in Table 3 and compared to the experimental values.16 The calculated standard binding free energies of HsTx1 agree with the experimental values within 1 kcal/mol. Thus, the chemical accuracy of the calculations is maintained in binding studies of HsTx1. Because HsTx1 has a different structure and binding modes compared to ShK, these results increase the confidence in the computational methods used to obtain the complex models and calculate the binding free energies. The binding free energies are also consistent with the deductions made by comparing the binding modes of HsTx1 with those of ShK. That is, HsTx1 has almost the same binding free energy for Kv1.3 as ShK but its affinity is reduced substantially for Kv1.1 and to a lesser degree for Kv1.2. Conversely, the observations made for the binding modes of HsTx1 above provide a mechanistic understanding for the calculated/measured selectivity of HsTx1 for Kv1.3 over Kv1.1. Analysis of the evolution of the contact distances in umbrella windows provides further insights for the binding modes; e.g., one can get a more realistic assessment for the strength of the interactions between pairs. Also, specific features of the PMFs can be explained by consulting the behavior of the pair distances. The results presented in Figure 7 are obtained from the average of the pair distances in each umbrella window used in the PMF calculations. We first consider the Kv1.1−HsTx1 complex, whose PMF is quite different from the PMFs we have seen so far for the potassium channel toxins. Instead of a sharp

HsTx1, deduced from a comparison of the binding modes, are in conformity with the measured binding constants of ShK49 and HsTx1 presented in Table 3. Because there are no mutagenesis data available for the Kv1.x−HsTx1 complexes, we rely on the binding free energies for validation of the complex models, which we consider next. Table 3. Comparison of Binding Free Energies for the Kv1.x−HsTx1 Complexesa complex

ΔGwell

Gb(PMF)

Gb(exp.)

Kv1.1−HsTx1 Kv1.2−HsTx1 Kv1.3−HsTx1

−13.5 ± 0.6 −12.3 ± 0.6 −17.2 ± 0.6

−10.1 ± 0.6 −8.9 ± 0.6 −14.0 ± 0.6

−11.1 ± 0.1 > −9.5 ± 0.1 −14.9 ± 0.2

a

The relative binding free energies obtained from the well depth in the PMFs are shown in the second column. The standard binding free energies determined from the PMFs using eqs 2 and 7 (third column) are compared to the experimental values in the last column.16 Errors in the binding free energies are estimated from the block data analysis of the PMF data. All energies are in kcal/mol.

Umbrella Sampling Simulations and PMFs of HsTx1. For each Kv1.x−HsTx1 complex, umbrella sampling MD simulations are performed as described in the Methods. An important issue in biased simulations is whether the bias forces applied during the pulling cause any permanent deformation of the toxin.31,32 Because HsTx1 is stabilized by four disulfide bonds, it is expected to have a relatively robust structure that will not be affected by the umbrella forces. This is demonstrated in Figure 4, where the average RMSD values of

Figure 4. Average RMSD of the HsTx1 backbone atoms with respect to the NMR structure obtained at each window of the umbrella sampling simulations. The results for Kv1.1, Kv1.2, and Kv1.3 are shown with red, green, and blue, respectively. The dotted line indicates the average RMSD obtained from the MD simulations of HsTx1 in the bulk.

the backbone atoms of HsTx1 relative to the NMR structure are plotted as a function of the channel−toxin distance. It is seen that, unlike in ShK, binding does cause a slight suppression of the RMSD fluctuations in HsTx1 in all three channels. We attribute this to the fact that the binding interface is formed by a β sheet, which is more flexible in the bulk and is stabilized upon binding. The important issue for the PMF calculations is whether binding causes any permanent deformation in the HsTx1 structure. The RMSDs obtained from the last window of the umbrella sampling simulations 712

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

Figure 6. Comparison of the PMFs for the unbinding of HsTx1 from the Kv1.1, Kv1.2, and Kv1.3 channels. The z distance is measured from the COM of HsTx1 to the COM of Kv1.x.

bulk region (Figure 6). From Figure 7, we see that K23 is the first residue to detach at z = 30 Å, indicating its relative weakness. It is followed by R14 at z = 31.5 Å. The remaining two residues (R4 and R33), on the other hand, do not completely break off until z > 36 Å. The side chains of R4 and R33 swing back and forth in several windows, which presumably contributes to the slow convergence of the PMF and its gradual rise. The Kv1.2−HsTx1 PMF exhibits the common behavior of the potassium channel toxins noted above (Figure 6). The pore inserting Lys (K23) keeps contact up to z = 31 Å, which is the typical range found in other toxins as well. Its breaking causes the kink at this point. From Figure 3 or Table 2, one could not distinguish between the strengths of the R14 and R27 charge interactions with D355. However, inspection of the distances in Figure 7 makes it clear that R27 detaches quickly, whereas the R14 contact persists up to z = 33 Å. Thus, R14 clearly provides the strongest interaction in the Kv1.2−HsTx1 binding mode. The R14−D355 contact is also responsible for the gradual rise of the PMF in the range z = 31−33 Å. Overall, Kv1.2−HsTx1 has fewer and lower quality contacts consistent with its lower affinity compared to the other complexes. The behavior of the Kv1.3−HsTx1 PMF is in between those of Kv1.1 and Kv1.2 (Figure 6). All the contacts are broken after z = 32 Å, except for R33 which detaches only after z = 36 Å (Figure 7). Thus, the PMF does rise sharply up to z = 32 Å, but after that, it keeps rising more gradually instead of exhibiting a shoulder as in Kv1.2. A distinguishing feature of the Kv1.3− HsTx1 binding mode is the even distribution of the contacts among the four monomers. This leads to a more stable complex structure as attested by the preservation of the contact distances during the first 3 Å pulling of the toxin (Figure 7). From the persistence lengths, we can determine the relative strengths of the various contacts. The R33 contact persists for 9 Å and is clearly the strongest interaction. It is followed by K7, which keeps contact for 5 Å. The other three residues (Y21, K23, and N26) follow similar trajectories, detaching after 3 Å, and hence should make similar contributions to the binding free energy. Comparison of these results with the persistence lengths obtained from the Kv1.3−ShK simulations (Figure 9 in ref 40) shows that indeed the quality of the contacts in Kv1.3−HsTx1 is much higher. Thus, fewer but higher quality contacts in Kv1.3−HsTx1 compared to Kv1.3−ShK are able to yield a

Figure 5. Convergence study of the Kv1.x−HsTx1 PMFs from block data analysis. To reduce fluctuations in PMFs, we use a large sampling size (2 ns), which is slided in 0.5 ns steps over the range of the data. The PMFs drop monotonically during equilibration and then fluctuate around a baseline after the system is equilibrated. The production data used in the final PMFs are indicated in parentheses in the figures.

rise followed by a shoulder region associated with the breaking of all the contacts, it exhibits a gradual rise all the way to the 713

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

similar binding free energy (−14.0 kcal/mol for HsTx1 to be compared with −14.2 kcal/mol for ShK40).



CONCLUSIONS



AUTHOR INFORMATION

Here we have presented a computational study of HsTx1 binding to Kv1 channels and provided a molecular-level description of its binding modes and a mechanistic basis for its Kv1.3/Kv1.1 selectivity. The binding free energies calculated from the PMFs agree with the experimental values within chemical accuracy, which provide strong support for the accuracy of the Kv1.x−HsTx1 complexes constructed. HsTx1 has already an excellent selectivity margin for Kv1.3 over Kv1.1. If necessary, this selectivity can be further improved through mutations of selected residues inferred from the binding modes. Thus, HsTx1 appears to provide an excellent basis for development of therapeutics for autoimmune diseases and warrants further study of its physiological and pharmacological properties to see if this goal is achievable. A uniform feature of the Kv1−toxin complexes studied here and elsewhere32,39−41,43 is that the pore inserting Lys forms Hbonds with the carbonyls of the Tyr residues in the filter. The Kv1.1−HsTx1 complex is unique in this regard because the pore inserting Lys (K23) is unable to fully insert into the filter and form H-bonds with the Tyr carbonyls. Three factors appear to have come together to yield this result; the β-sheet interface of HsTx1, the lack of cross-linking between Y379 and D377 in Kv1.1, which results in the bulky Tyr side chains sticking out of the pore; and the strong coupling between the three Arg side chains with the Asp side chains in the turret. Together, these steric factors and charge interactions prevent HsTx1 from coming close enough to the pore domain in Kv1.1 to enable a proper insertion of the K23 side chain into the filter. The residues corresponding to Tyr in Kv1.1 are Val in Kv1.2 and His in Kv1.3. As noted before, the H404 side chains in Kv1.3 form cross-links with D402 and the V404 side chains in Kv1.2 are not as bulky as those of Tyr. In addition, in neither the Kv1.2 nor the Kv1.3 complex three Arg side chains couple to the turret residues, holding it back as it happens in the Kv1.1 complex. Thus, despite the presence of the β-sheet, HsTx1 can still come sufficiently close to the pore domain in Kv1.2 and Kv1.3 to enable a proper insertion of the K23 side chain. These observations have some implications on the search for a Kv1.3-selective ShK analogue. In our previous study of ShK binding to Kv1 channels,40 we suggested that reducing the charge contacts between ShK analogues and Kv1.1 through the mutations K18A and R29A should improve its selectivity for Kv1.3. The results from HsTx1 binding suggest that preventing a proper insertion of the pore inserting Lys in Kv1.1 could also help with the selectivity problem. This may be achieved by including bulky residues on the binding interface of ShK, which will prevent the ShK analogue from coming near the pore domain of Kv1.1. Figure 7. Average distance between the interacting pairs of HsTx1 and Kv1.x plotted as a function of the window position. All the important pairs identified in Table 2 are considered. The trajectory data are taken from the umbrella sampling simulations that have been used in the construction of the final PMFs in Figure 6. Error bars are not shown to avoid cluttering, but they are quite small at contact distances (see Table 2).

Corresponding Author

*E-mail: [email protected]. Phone: ++61 2-90365306. Notes

The authors declare no competing financial interest. 714

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B



Article

(17) Gordon, D.; Chen, R.; Chung, S. H. Computational Methods of Studying the Binding of Toxins from Venomous Animals to Biological Ions Channels: Theory and Applications. Physiol. Rev. 2013, 93, 767− 802. (18) Rashid, M. D.; Mahdavi, S.; Kuyucak, S. Computational Studies of Marine Toxins Targeting Ion Channels. Mar. Drugs 2013, 11, 848− 869. (19) Shiau, Y. S.; Lin, T. B.; Liou, H. H.; Huang, P. T.; Lou, K. L.; Shiau, Y. Y. Molecular Simulation Reveals Structural Determinants of the Hanatoxin Binding in Kv2.1 Channels. J. Mol. Model. 2002, 8, 253−257. (20) Andreotti, N.; di Luccio, E.; Sampieri, F.; De Waard, M.; Sabatier, J. M. Molecular Modeling and Docking Simulations of Scorpion Toxins and Related Analogs on Human SKCa2 and SKCa3 Channels. Peptides 2005, 26, 1095−1108. (21) Fu, W.; Cui, M.; Briggs, J. M.; Huang, X.; Xiong, B.; Zhang, Y.; Luo, X.; Shen, J.; Ji, R.; Jiang, H.; Chen, K. Brownian Dynamics Simulations of the Recognition of the Scorpion Toxin Maurotoxin with the Voltage-Gated Potassium Ion Channels. Biophys. J. 2002, 83, 2370−2385. (22) Yu, K.; Fu, W.; Liu, H.; Luo, X.; Chen, K. X.; Ding, J. P.; Shen, J.; Jiang, H. Computational Simulations of Interactions of Scorpion Toxins with the Voltage-Gated Potassium Ion Channel. Biophys. J. 2004, 86, 3542−3555. (23) Eriksson, M. A.; Roux, B. Modeling the Structure of Agitoxin in Complex with the Shaker K+ Channel: A Computational Approach Based on Experimental Distance Restraints Extracted from Thermodynamic Mutant Cycles. Biophys. J. 2002, 83, 2595−2609. (24) Huang, X.; Dong, F.; Zhou, H. X. Electrostatic Recognition and Induced Fit in the κ-PVIIA Toxin Binding to Shaker Potassium Channel. J. Am. Chem. Soc. 2005, 127, 6836−6849. (25) Alonso, H.; Bliznyuk, A. A.; Gready, J. E. Combining Docking and Molecular Dynamic Simulations in Drug Design. Med. Res. Rev. 2006, 26, 531−568. (26) Zarrabi, M.; Naderi-Manesh, H. The investigation of Interactions of κ-Hefutoxin1 with the Voltage-Gated Potassium Channels: A Computational Simulation. Proteins 2007, 71, 1441− 1449. (27) Yi, H.; Qiu, S.; Cao, Z.; Wu, Y.; Li, W. Molecular Basis of Inhibitory Peptide Maurotoxin Recognizing Kv1.2 Channel Explored by ZDOCK and Molecular Dynamic Simulations. Proteins 2008, 70, 844−854. (28) Yi, M. G.; Tjong, H. T.; Zhou, H. X. Spontaneous Conformational Change and Toxin Binding in α7 Acetylcholine Receptor: Insight into Channel Activation and Inhibition. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 8280−8285. (29) Pietra, F. Docking and MD Simulations of the Interaction of the Tarantula Peptide Psalmotoxin-1 with ASIC1a Channels Using a Homology Model. J. Chem. Inf. Model. 2009, 49, 972−977. (30) Ander, M.; Luzhkov, V. B; Aqvist, J. Ligand Binding to the Voltage-Gated Kv1.5 Potassium Channel in the Open State - Docking and Computer Simulations of a Homology Model. Biophys. J. 2008, 94, 820−831. (31) Chen, P. C.; Kuyucak, S. Mechanism and Energetics of Charybdotoxin Unbinding from a Potassium Channel from Molecular Dynamics Simulations. Biophys. J. 2009, 96, 2577−2588. (32) Chen, P. C.; Kuyucak, S. Accurate Determination of the Binding Free Energy for KcsA-Charybdotoxin Complex from the Potential of Mean Force Calculations with Restraints. Biophys. J. 2011, 100, 2466− 2474. (33) Khabiri, M.; Nikouee, A.; Cwiklik, L.; Grissmer, S.; Ettrich, R.; Charybdotoxin, R. Unbinding from the mKv1.3 Potassium Channel: A Combined Computational and Experimental Study. J. Phys. Chem. B 2011, 115, 11490−11500. (34) Chen, R.; Robinson, A.; Gordon, D.; Chung, S. H. Modeling the Binding of three Toxins to the Voltage-Gated Potassium Channel (Kv1.3). Biophys. J. 2011, 101, 2652−2660.

ACKNOWLEDGMENTS Calculations were performed using the HPC facilities at the National Computational Infrastructure (Canberra). Discussions on HsTx1 toxin with Ray Norton are acknowledged with thanks.



REFERENCES

(1) Hille, B. Ionic Channels of Excitable Membranes, 3rd ed.; Sinauer Associates: Sunderland, MA, 2001. (2) Ashcroft, F. M. Ion Channels and Disease: Channelopathies; Academic Press: San Diego, CA, 2000. (3) Wulff, H.; Zhorov, B. S. K+ Channel Modulators for the Treatment of Neurological Disorders and Autoimmune Diseases. Chem. Rev. 2008, 108, 1744−1773. (4) Wulff, H.; Castle, N. A.; Pardo, L. A. Voltage-Gated Potassium Channels as Therapeutic Targets. Nat. Rev. Drug Discovery 2009, 8, 982−1001. (5) Wulff, H.; Calabresi, P. A.; Allie, R.; Yun, S.; Pennington, M.; Beeton, C.; Chandy, K. G. The Voltage-Gated Kv1.3 K+ Channel in Effector Memory T Cells as New Target for MS. J. Clin. Invest. 2003, 111, 1703−1713. (6) Beeton, C.; Wulff, H.; Standifer, N. E.; Azam, P.; Mullen, K. M.; Pennington, M. W.; Kolski-Andreaco, A.; Wei, E.; Grino, A.; Counts, D. R.; et al. Kv1.3 Channels are a Therapeutic Target for T CellMediated Autoimmune Diseases. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 17414−17419. (7) Cahalan, M. D.; Chandy, K. G. The Functional Network of ion Channels in T Lymphocytes. Immunol. Rev. 2009, 231, 59−87. (8) Schmitz, A.; Sankaranarayanan, A.; Azam, P.; Schmidt-Lassen, K.; Homerick, D.; Hansel, W.; Wulff, H. Design of PAP-1, a Selective Small Molecule Kv1.3 Blocker, for the Suppression of Effector Memory T Cells in Autoimmune Diseases. Mol. Pharmacol. 2005, 68, 1254−1270. (9) Norton, R. S.; Pennington, M. W.; Wulff, H. Potassium Channel Blockade by the Sea Anemone Toxin ShK for the Treatment of Multiple Sclerosis and Other Autoimmune Diseases. Curr. Med. Chem. 2004, 11, 3041−3052. (10) Chi, V.; Pennington, M. W.; Norton, R. S.; Tarcha, E. J.; Londono, L. M.; Sims-Fahey, B.; Upadhyay, S. K.; Lakey, J. T.; Iadonato, S.; Wulff, H.; et al. Development of a Sea Anemone Toxin as an Immunomodulator for Therapy of Autoimmune Diseases. Toxicon 2011, 59, 529−546. (11) Castaneda, O.; Sotolongo, V.; Amor, A. M.; Stocklin, R.; Anderson, A. J.; Harvey, A. L.; Engstrom, A.; Wernstedt, C.; Karlsson, E. Characterization of a Potassium Channel Toxin from the Caribbean Sea Anemone Stichodactyla helianthus. Toxicon 1995, 33, 603−613. (12) Kalman, K.; Pennington, M. W.; Lanigan, M. D.; Nguyen, A.; Rauer, H.; Mahnir, V.; Paschetto, K.; Kem, W. R.; Grissmer, S.; Gutman, G. A. ShK-Dap22, a Potent Kv1.3-Specific Immunosuppressive Polypeptide. J. Biol. Chem. 1998, 273, 32697−32707. (13) Beeton, C.; Pennington, M. W.; Wulff, H.; Singh, S.; Nugent, D.; Crossley, G.; Khaytin, I.; Calabresi, P. A.; Chen, C. Y.; Gutman, G. A.; Chandy, K. G. Targeting Effector Memory T Cells with a Selective Peptide Inhibitor of Kv1.3 Channels for Therapy of Autoimmune Diseases. Mol. Pharmacol. 2005, 67, 1369−1381. (14) Pennington, M. W.; Beeton, C.; Galea, C. A.; Smith, B. J.; Chi, V.; Monaghan, K. P.; Garcia, A.; Rangaraju, S.; Giuffrida, A.; Plank, D.; et al. Engineering a Stable and Selective Peptide Blocker of the Kv1.3 Channel in T Lymphocytes. Mol. Pharmacol. 2009, 75, 762−773. (15) Lebrun, B.; Romi-Lebrun, R.; Martin-Eauclaire, M. F.; Yasuda, A.; Ishiguro, M.; Oyama, Y.; Pongs, O.; Nakajima, T. A Four DisulfideBridged Toxin with High Affinity towards Voltage-Gated K+ Channels Isolated from Heterometrus Spinnifer (Scorpionidea) Venom. Biochem. J. 1997, 328, 321−327. (16) Regaya, I.; Beeton, C.; Ferrat, G.; Andreotti, N.; Darbon, H.; De Waard, M.; Sabatier, J. M. Evidence for Domain Specific Recognition of SK and Kv Channels by MTX and HsTx1 Scorpion Toxins. J. Biol. Chem. 2004, 279, 55690−55696. 715

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716

The Journal of Physical Chemistry B

Article

(35) Chen, R.; Chung, S. H. Engineering a Potent and Specific Blocker of Voltage-Gated Potassium Channel Kv1.3, a Target for Autoimmune Diseases. Biochemistry 2012, 51, 7775−7782. (36) Chen, R.; Chung, S. H. Structural Basis of the Selective Block of Kv1.2 by Maurotoxin from Computer Simulations. PLoS One 2012, 7, e47253. (37) Chen, R.; Chung, S. H. Binding Modes of μ-Conotoxin to the Bacterial Sodium Channel (NaVAb). Biophys. J. 2012, 102, 483−488. (38) Chen, R.; Chung, S. H. Conserved Functional Surface of Antimammalian Scorpion β-Toxins. J. Phys. Chem. B 2012, 116, 4796− 4800. (39) Chen, P. C.; Kuyucak, S. Developing a Comparative Docking Protocol for the Prediction of Peptide Selectivity Profiles: Investigation of Potassium Channel Toxins. Toxins 2012, 4, 110−138. (40) Rashid, M. H.; Kuyucak, S. Affinity and Selectivity of ShK Toxin for the Kv1 Potassium Channels from Free Energy Simulations. J. Phys. Chem. B 2012, 116, 4812−4822. (41) Pennington, M. W.; Rashid, M. H.; Tajhya, R. B.; Beeton, C.; Kuyucak, S.; Norton, R. S. A C-Terminally Amidated Analogue of ShK is a Potent and Selective Blocker of the Voltage-Gated Potassium Channel Kv1.3. FEBS Lett. 2012, 586, 3996−4001. (42) Rashid, M. H.; Heinzelmann, G.; Huq, R.; Tajhya, R. B.; Chang, S. C.; Chhabra, S.; Pennington, M. W.; Beeton, C.; Norton, R. S.; Kuyucak, S. A Potent and Selective Peptide Blocker of the Kv1.3 Channel: Prediction from Free-Energy Simulations and Experimental Confirmation. PLoS One 2013, 8, e78712. (43) Mahdavi, S.; Kuyucak, S. Why Drosophila Shaker K+ Channel is not a Good Model for Ligand Binding to Voltage-Gated Kv1 Channels. Biochemistry 2013, 52, 1631−1640. (44) Gilson, M. K.; Zhou, H. X. Calculation of Protein-Ligand Binding Energies. Annu. Rev. Biophys. Biomol. Struct. 2007, 36, 21−42. (45) Deng, Y.; Roux, B. Computations of Standard Binding Free Energies with Molecular Dynamics Simulations. J. Phys. Chem. B 2009, 113, 2234−2246. (46) Steinbrecher, T.; Labahn, A. Towards Accurate Free Energy Calculations in Ligand-Protein Binding Studies. Curr. Med. Chem. 2010, 17, 767−785. (47) Yu, L. P.; Sun, C. H.; Song, D. Y.; Shen, J. W.; Xu, N.; Gunasekera, A.; Hajduk, P. J.; Olejniczak, E. T. Nuclear Magnetic Resonance Structural Studies of a Potassium Channel-Charybdotoxin Complex. Biochemistry 2005, 44, 15834−15841. (48) Long, S. B.; Tao, X.; Campbell, E. B.; MacKinnon, R. Atomic Structure of a Voltage-Dependent K+ Channel in a Lipid MembraneLike Environment. Nature 2007, 450, 376−382. (49) Rauer, H.; Pennington, M.; Cahalan, M.; Chandy, K. G. Structural Conservation of the Pores of Calcium-Activated and Voltage-Gated Potassium Channels Determined by a Sea Anemone Toxin. J. Biol. Chem. 1999, 274, 21885−21892. (50) Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33−38. (51) Bastug, T.; Kuyucak, S. Importance of the Peptide Backbone Description in Modeling the Selectivity Filter in Potassium Channels. Biophys. J. 2009, 96, 4006−4012. (52) Savarin, P.; Romi-Lebrun, R.; Zinn-Justin, S.; Lebrun, B.; Nakajima, T.; Gilquin, B.; Menez, A. Structural and Functional Consequences of the Presence of a Fourth Disulfide Bridge in the Scorpion Short Toxins: Solution Structure ot the Potassium Channel Inhibitor HsTx1. Protein Sci. 1999, 8, 2672−2685. (53) Mouhat, S.; Jouirou, B.; Mosbah, A.; De Waard, M.; Sabatier, J. M. Animal Toxins Acting on Voltage-Gated Potassium Channels. Biochem. J. 2004, 378, 717−726. (54) Carrega, L.; Mosbah, A.; Ferrat, G.; Beeton, C.; Andreotti, N.; Mansuelle, P.; Darbon, H.; De Waard, M.; Sabatier, J. M. The Impact of the Fourth Disulfide Bridge in Scorpion Toxins of the α-KTx6 Subfamily. Proteins 2005, 61, 1010−1023. (55) Dominguez, C.; Boelens, R.; Bonvin, A. M. J. J. HADDOCK: A Protein-Protein Docking Approach Based on Biochemical or Biophysical Information. J. Am. Chem. Soc. 2003, 125, 1731−1737.

(56) de Vries, S.; J.; van Dijk, A. D. J.; Krzeminski, M.; van Dijk, M.; Thureau, A.; Hsu, V.; Wassenaar, T.; Bonvin, A. M. J. J. HADDOCK versus HADDOCK: New Features and Performance of HADDOCK 2.0 on the CAPRI Targets. Proteins 2007, 69, 726−733. (57) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. NAMD2: Greater Scalability for Parallel Molecular Dynamics. J. Comput. Chem. 2005, 26, 1781−1802. (58) Klauda, J. B.; Venable, R. M; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2012, 114, 7830−7843. (59) Kumar, S.; Bouzida, D.; Swensen, R. H.; Kollman, P. A.; Rosenberg, J. M. The Weighted Histogram Analysis Method for FreeEnergy Calculations on Biomolecules. J. Comput. Chem. 1992, 13, 1011−1021. (60) Roux, B.; Allen, T.; Berneche, S.; Im, W. Theoretical and Computational Models of Biological Ion Channels. Q. Rev. Biophys. 2004, 37, 15−103.

716

dx.doi.org/10.1021/jp410950h | J. Phys. Chem. B 2014, 118, 707−716