**2. Computational Methods**

There are only a few crystal structures for complexes of membrane proteins. Thus the first step in a computational study of toxin binding to ion channels is the construction of complex structures. Here accuracy of the model structure is of utmost importance because without an accurate complex model, free energy calculations in the next step have no chance of succeeding. Accordingly, we first discuss the computational methods used in structure prediction followed by the free energy methods.

### *2.1. Complex Structure Prediction from Docking and MD Simulations*

In order to find the structure of a channel-toxin complex, one first needs the individual structures of the channel and the toxin. Those of toxins can be determined using NMR in a straightforward manner, and many toxin structures are available from the protein data bank. Structures of channel proteins are determined from X-ray crystallography, and because it is much harder to crystallize membrane proteins, not many channel structures are available. Hence, one has to rely on homology modeling in most cases. The situation is relatively better in potassium channels where several crystal structures exist [12], including the mammalian voltage-gated potassium channel Kv1.2 [13]. Thus one can construct homology models of other Kv1 channels relatively easily, although it is still very important to validate such model channels using available functional data. In sodium channels, the crystal structure of a bacterial channel was determined recently [14], which has opened the way for homology modelling of the mammalian Nav1 channels. Due to substantial differences between the bacterial and mammalian sodium channels, proper validation of the homology models is even more important in this case.

Two main methods for prediction of protein-ligand complexes are docking and MD simulations. Docking programs allow fast screening of many ligands for a given target [23,24], but their accuracy is limited [25]. Conversely, MD simulations provide accurate representation of the protein-ligand interactions but they are too slow to predict the complex structure from scratch. Combination of the two methods, where the initial binding poses predicted by a docking method are refined in subsequent MD simulations, offers the most practical approach for finding complex structures. Initial applications of this approach to small ligands (< 50 atoms) using common docking programs such as AUTODOCK [26] and ZDOCK [27] produced promising results [18,28–30]. Its feasibility for larger and more flexible peptide ligands was first shown for the KcsA potassium channel-charybdotoxin complex [31,32]. The structure of this complex was determined from NMR experiments [33], so it could be used for testing the accuracy of the docking plus MD approach and the effectiveness of the docking programs [32]. Using a more sophisticated docking program such as HADDOCK [34,35], which allowed flexibility and ensemble docking, was found to give superior results compared with rigid docking programs [28,32].

Docking programs rank the complex poses according to their energy score. Inspection of the top 10 or 20 poses is usually sufficient to make a decision on an initial pose. A systematic study of potassium channel-toxin complexes using HADDOCK has shown that a consensus complex is obtained in most cases [32]. We assume this to be the case in the following, but if more than one pose is found from docking, each pose needs to be refined with MD. After an initial configuration for the complex structure is chosen from docking, one needs to prepare the simulation system for refinement with MD. There are well-established protocols for this purpose [36,37], which can be used for the complex model as well. Typically, the complex model is embedded in a lipid bilayer and solvated with a salt solution using the VMD program [38]. The resulting system is gradually relaxed in MD simulations until it reaches equilibrium. Care needs to be exercised to ensure that all the disulfide and hydrogen bonds in the complex structure are preserved during relaxation. Once the system is well equilibrated in MD simulations, a production run is performed for analysis of the binding mode. Snapshots of the complex system can be used to get a visual picture of the binding mode. A more quantitative description can be obtained by calculating the average pair distances for the strongly interacting residues. Charge interactions, where the N–O distance between the charged residues is less than 3 Å, and hydrophobic interactions involving the benzyl groups provide the strongest couplings (> 2 kcal/mol). Hydrogen bonds and charge interactions at larger distances are of intermediate strength (1–2 kcal/mol). Where available, comparison of the alanine scanning mutagenesis data with the predictions of the complex model provides the best validation for the proposed model. If no mutation data are available, one has to rely on binding free energies, which is discussed in the next section.

There are several MD programs that one can use for refinement of the complex. In the examples discussed here, the NAMD code [39] is used with the CHARMM force field [16] including the CMAP correction for the dihedral terms [40]. The NpT ensemble is the most suitable one for MD simulations of biomolecules and has been adapted in most MD studies. The temperature and pressure can be maintained at the standard values of 300 K and 1 atm using temperature and pressure coupling. Employment of the periodic boundary conditions avoids artefacts arising from using small boundary boxes and also facilitates computation of the long-range electrostatic interactions. Neglecting the long-range electrostatic interactions using cut-off distances causes errors and is not recommended. The current practice is to include them using the particle-mesh Ewald algorithm. The short-range Lennard-Jones interactions can be switched off within a distance of 10–13.5 Å without causing errors. A typical time step used in the MD simulations is 2 fs. Using longer time steps results

in accumulation of systematic errors while shorter ones require more computing time and are not used unless extreme accuracy is desired. For details of the basic formalism of MD simulations, we refer to the monographs [41–43]. A recent review of MD simulations of membrane proteins can be found in [44].
