PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
Mol Simul. Author manuscript; available in PMC 2010 November 23.
Published in final edited form as:
Mol Simul. 2009 September; 35(10-11): 986–997.
doi:  10.1080/08927020902902742
PMCID: PMC2990536
NIHMSID: NIHMS247424

Relative free energy of binding between antimicrobial peptides and SDS or DPC micelles

Abstract

We present relative binding free energy calculations for six antimicrobial peptide–micelle systems, three peptides interacting with two types of micelles. The peptides are the scorpion derived antimicrobial peptide (AMP), IsCT and two of its analogues. The micelles are dodecylphosphatidylcholine (DPC) and sodium dodecylsulphate (SDS) micelles. The interfacial electrostatic properties of DPC and SDS micelles are assumed to be similar to those of zwitterionic mammalian and anionic bacterial membrane interfaces, respectively. We test the hypothesis that the binding strength between peptides and the anionic micelle SDS can provide information on peptide antimicrobial activity, since it is widely accepted that AMPs function by binding to and disrupting the predominantly anionic lipid bilayer of the bacterial cytoplasmic membrane. We also test the hypothesis that the binding strength between peptides and the zwitterionic micelle DPC can provide information on peptide haemolytic activities, since it is accepted that they also bind to and disrupt the zwitterionic membrane of mammalian cells. Equilibrium structures of the peptides, micelles and peptide–micelle complexes are obtained from more than 300 ns of molecular dynamics simulations. A thermodynamic cycle is introduced to compute the binding free energy from electrostatic, non-electrostatic and entropic contributions. We find relative binding free energy strengths between peptides and SDS to correlate with the experimentally measured rankings for peptide antimicrobial activities, and relative free energy binding strengths between peptides and DPC to correlate with the observed rankings for peptide haemolytic toxicities. These findings point to the importance of peptide–membrane binding strength for antimicrobial activity and haemolytic activity.

Keywords: antimicrobial peptides, IsCT, binding free energy calculations, molecular dynamics simulations

1. Introduction

The number of reports of drug-resistant strains of bacteria that cause infectious diseases such as pneumonia, tuberculosis, syphilis and gonorrhoea is rapidly rising [15]. Antimicrobial peptides (AMPs) are amongst new, attractive candidates for combating bacterial infections. AMPs are recognised as important components of the immune defence system in many species including insects, plants and mammalian organisms. They are structurally diverse and many of them are highly selective against bacteria, including current antibiotic resistant strains. The bactericidal effects of these peptides are usually fast, lysing bacterial cells within few minutes [6,7].

Numerous experimental and computational studies have been providing insight into the molecular mechanism of action of AMPs. It is widely accepted that many of the naturally occurring AMPs, which are predominantly cationic, kill pathogens through non-specifically targeting and disrupting the anionic bacterial lipid membrane[822].

Despite their promise, an important constraint for using known AMPs in clinical applications is their haemolytic capacity against human erythrocytes. This toxic profile of AMPs limits their systemic administration for therapeutic purposes. This may be puzzling since eukaryotic cytoplasmic membranes are rich in a mixture of zwitterionic phospholipids [23], in contrast to bacterial cells which have a negative charge density on their outer surface due to a large fraction of anionic lipids [24,25]. Contrary then to the expectation that cationic AMPs would penetrate and disrupt only bacterial membranes, they also disrupt mammalian cell membranes. These results point to the complexity of the mechanism of antimicrobial activity and toxicity as well as to the inherent experimental limitations on obtaining a mechanistic picture of exactly how these peptides interact and disrupt the membranes of both bacterial and mammalian membranes. Indeed, more than 30 years after the discovery of first AMPs there is no universally agreed upon model for antimicrobial activity, let alone toxicity.

Here, we investigate the viability of calculating relative free energy of binding between AMPs and model membrane mimics to explain variable AMP antimicrobial activity and haemolytic toxicity. We use the scorpion derived AMP IsCT [26] and two of its analogues that exhibit variable activity and toxicity profiles as a test case. These are attractive model systems because there is unambiguous structural and functional information: the structure of the peptides has been determined with NMR, their activity has been tested against many bacterial species and their hemolycity has been quantified against human erythrocytes.

IsCT (ILGKI WEGIK SLF) is a non-selective linear peptide, isolated from the venom of the African scorpion Opisthacanthus madagascariensis. The molecular weight of IsCT is 1501.9 Da and its C-terminal is found to be naturally amidated. It assumes an α-helical structure with an amphipathic conformation when immersed in tri-fluoroethanol (TFE) [26]. The short length of its amino acid sequence provides an attractive target for mutagenesis studies of AMPs. The wild-type IsCT is toxic to both bacterial and eukaryotic cells [26]. Lee et al. [27] engineered several analogues of IsCT with point mutations and studied their antibacterial and haemolytic activities. The antibacterial activity measurements were obtained using the broth microdilution method, while haemolysis percentage measurements using absorbance at 405 nm were used to evaluate the haemolytic toxicity of these analogues [27]. The percent haemolysis was measured for IsCT and the two mutants, IsCT2 (ILGKI WKGIK SLF) and IsCT3 (ILGKI WKPIK KLF) at varying concentrations. At 100 μM concentration, the haemolytic toxicities of the three peptides were found to be approximately 78, 52 and 20%, respectively. The antimicrobial activity of IsCT and its analogues was also examined against both Gram-positive and Gram-negative bacteria by measuring the minimal inhibitory concentration (MIC) of the peptides in different microbial cultures including Escherichia coli, Pseudomonas aeruginosa, Salmonella typhimurium, Staphylococcus epidermidis and Staphylococcus aureus [24]. The antimicrobial activities of the three peptides are shown in Table 1. The activities are partially dependent on the organism. Overall, IsCT2 has the highest antimicrobial activity, closely followed by IsCT3. The only significant difference between IsCT2 and IsCT3 is their antimicrobial activity against S. aureus, where the MICs are measured to be 1 and 2 μM, respectively. The first analogue, IsCT2, was synthesised to increase the net positive charge by replacing Glu7 by Lys7, whereas the second analogue, IsCT3, was synthesised to provide the hinge region in the central position of the peptide as well as increasing the net positive charge of the peptide by replacing Glu7, Gly8 and Ser11 by Lys7, Pro8 and Lys11, respectively.

Table 1
The antimicrobial activities of IsCT, IsCT2 and IsCT3 are measured by the MICs of the peptides in different microbial cultures.

Our goal here is to investigate if relative binding affinities of peptides to dodecylphosphatidylcholine (DPC) and sodium dodecylsulphate (SDS) micelles can reproduce the rankings of the observed differential haemolytic toxicities toward eukaryotic cells and antimicrobial activities against bacterial cells. We use the molecular mechanics Poisson–Boltzmann surface area (MM-PBSA) protocol [28] to compute the relative binding free energies of IsCT and two of its analogues to SDS and DPC micelles. These are considered simple, computationally inexpensive models of bacterial and mammalian membranes, respectively [15,19].

We report that the computed relative binding free binding energies of IsCT analogues to SDS and DPC micelles are indeed correlated with the rankings of the available experimental measurements on antimicrobial activity and haemolytic toxicity, respectively. Theoretical and implementation details of the calculation methods are discussed in the following two subsections and in the subsection on relative binding free energy calculations.

2. Methods

2.1 Initial structures

The starting coordinates of the DPC micelle–water complex were extracted from simulations carried out by Wymore et al. [29]. This structure was obtained after extensive minimisation and dynamics of about 1 ns in a cubic simulation cell. The force field parameters for DPC molecules were obtained from the work of Wong and Kamath [30]. Similarly, the parameters for the force field and the initial coordinates of the SDS micelle–water complex were taken from molecular dynamic (MD) simulations carried out by MacKerell [31].

The experimental structures of IsCT and its analogues were acquired using NMR spectroscopy at a peptide molar concentration of 1.0 mM dissolved in 150 mM SDS micellar solution [27]. The initial coordinates of IsCT and its analogues were obtained from the RCSB Protein Data Bank (pdb ids: 1T51, 1T52 and 1T55). The ionisable groups were protonated or deprotonated assuming the pH for the unbound peptide and micelle–peptide solutions to be equal to 7. All glumatic acid groups were deprotonated, while lysine and NH2 terminal groups were protonated. The C-termini were blocked by neutralising its negative charge through amidation. Thus, the net charges of IsCT, IsCT2 and IsCT3 are equal to +2e, +4e and +5e, respectively.

2.2 Molecular dynamic simulations

MD simulations for each one of the three peptides were carried out in water, SDS and DPC environments. These simulations were conducted using the CHARMM program [32] version c30b2 and the all hydrogen force field CHARMM27 [33,34]. Water was modelled using the TIP3P force potential [35]. In the simulations of peptide–micelle complexes, we placed the micelle consisting of 60 DPC molecules and 4377 waters in a cubic simulation box of cell size 56.153 Å3, and the micelle consisting of 60 SDS molecules and 4375 water molecules in a cubic box of cell size 54.153 Å3. The cell dimensions were setup to obtain the equilibrium bulk water number density of 0.033 mol Å−3 several angstroms away from the water–micelle interface.

The simulations of individual peptides were performed using a box of size 53.93 Å3 containing 5314 water molecules. To attain charge neutrality in each simulation box and account for the salinity of the physiological medium (~150 mM), appropriate numbers of sodium and chloride ions were randomly distributed in the aqueous phase.

In earlier studies, we have shown that simulations converge to their final equilibrium state with respect to the distance between the centres of mass of the peptide and the micelle irrespective of the initial peptide position [14,15], and recently, we demonstrated that simulations are in excellent agreement with FTIR experiments [19]. Without loss of generality, here each peptide in the peptide–SDS simulations was placed in the simulation box with its centre of mass coinciding with that of the SDS micelle. In the peptide–DPC simulations, each peptide was placed in the aqueous phase parallel to the DPC micelle interface. Ultimately, we see the peptides position equilibrate on the surface of the micelles. Since, both SDS and DPC micelles are initially spherically symmetric, the initial orientation of the peptide with respect to the radial direction is irrelevant. The fast relaxation times of SDS and DPC molecules render computationally feasible the repositioning of each peptide to its final conformation.

To remove initial close contacts between a peptide and an SDS or DPC micelle as well as to prevent the penetration of water during equilibration, all systems were subjected to several stages of minimisation using the steepest descent method. Initially, each peptide and the bulk water were kept under weak harmonic constraints with spring constants of 10 and 5 kcal/mol1 Å2, respectively. The constraints were gradually decreased every 20,000 steepest decent minimisation cycles. Consequently, each system was minimised for an additional 20,000 minimisation steps without any constraints. After the minimisation stage, each system consisting of more than 16,000 atoms was gradually heated to 303.15 K. Thereafter, each MD simulation was carried out under NPT conditions (1 atm and 303.15 K) using the extended system pressure [36] and the Hoover temperature control [37], respectively. For the extended system pressure control algorithm, all the components of the system mass were set to 500 amu. The electrostatic interactions were approximated using the particle mesh Ewald (PME) summation method [38,39] without truncation and a real space Gaussian width of 0.25 Å−1, a β-spline of order 4 and an FFT grid of approximately one point per angstrom. The SHAKE algorithm was applied to constraint the fastest degrees of freedom involving bonds with hydrogen atoms, thus eliminating high-frequency vibrations and allowing 2 fs integration time steps. The non-bonded van der Waals interactions were smoothly switched off over a distance of 3.0 Å, between 9 and 12 Å. The simulations were carried out using periodic boundary conditions and the leap-frog integrator. The simulation times were approximately 20 ns for the individual peptide, 40 ns for the peptide–SDS complexes and 50 ns for the peptide–DPC complexes.

2.3 Relative binding free energy calculations

The binding free energy to form a complex is the reversible work required for bringing the two solvated binding molecules from an infinite distance to the bound state. There is no direct means for calculating the free energy of a molecular system’s state. Approximate binding free energy calculation methods, such as thermodynamic integration and perturbation methods, perform well for small systems, but are of limited use for larger systems [40].

An attractive and computationally plausible approach to estimate the binding free energy for larger systems, like the ones studied here (Figure 1), is using the end point component analysis and continuum solvent models [28]. This method takes advantage of the fact that the free energy is a state function; as such, free energy differences between two states are independent of the path the system follows between these two states. A convenient path, or thermodynamic cycle, is as follows: the system of solvated molecules is first transferred to the gas phase, the binding can occur in the gas phase and then the complex is solvated back into the solvent.

Figure 1
IsCT–DPC complex.

In particular, considering the association of a peptide, P, to a surfactant micelle, M, to form a bound peptide–micelle complex (Figure 1), the binding free energy of the complexation process in aqueous solution, Paq + Maq [right arrow over left arrow] PMaq, is

ΔGbind=GPMGPGM.
(1)

A direct computation of each one of these free energy terms is not computationally feasible. Figure 2 shows the convenient thermodynamic cycle that is considered to decompose the free energy, G, of a molecule into gas phase molecular mechanical potential energy (UMM), solvation free energy (ΔGsolv) and solute entropy (SMM) contributions as follows:

G=UMM+ΔGsolvTSMM.
(2)
Figure 2
An example of a thermodynamic cycle that can be used to compute binding free energy of a peptide to a surfactant micelle.

The angle bracket, left angle bracket···right angle bracket, indicates an average over a set of snapshots along computed molecular dynamics trajectories. The molecular potential energy, UMM, is a sum of intermolecular and intramolecular interaction terms, calculated from a classical molecular force field (CHARMM [33]) and given by

UMM=Ubond+Uangle+Udihedral+Uvdw+Uelec,
(3)

where Ubond, Uangle, Udihedral, Uvdw and Uelec are the bond, angle, dihedral angle, van der Waals and electrostatic terms of the potential energy, respectively.

The binding free energy of a peptide–micelle complex can then be computed using the thermodynamic cycle (Figure 2) as follows:

ΔGbind=ΔΔGsolv+ΔGgas,
(4)

where the solvation free energy difference, ΔΔGsolv, is given by

ΔΔGsolv=ΔGsolvPMΔGsolvPΔGsolvM.
(5)

The gas phase binding free energy, ΔGgas, is computed using the following formula,

ΔGgas=(UMMPMUMMPUMMM)T(SMMPMSMMPSMMM).
(6)

The entropic part of the free energy can be decomposed into translational, rotational, conformational and vibrational contributions [41,42]. Translational and rotational entropies are computed by utilising polyatomic ideal gas formulas derived using statistical mechanics principles. Specifically, the vibrational entropy is calculated using normal modes frequencies [43], while the loss of conformational entropy of each peptide side chain upon complexation is estimated using the empirical scale of Pickett and Sternberg [44]. The electrostatic contributions to the solvation free energies of each complex component and the complex itself are calculated using the Poisson–Boltzmann method [45]. The non-polar contributions to the free energy are approximated using a linear solvent-accessible surface area relationship [28,46]. The latter implies that the solvation free energy can be estimated using electrostatic and non-polar contributions as follows:

ΔGsolv=ΔGsolvelec+γΔSASA,
(7)

where ΔSASA is the difference in the solvent accessible area between the complex and its components, while γ is an empirical solvation parameter which is set to 7.2 (cal mol−1 Å−2) [28,46]. The non-polar contribution to the binding free energy accounts for the loss of solvent entropy due to the formation of a cavity and due to the presence of the solute as well as other effects such as the van der Waal interactions between the solvent and solute [46].

Ultimately, of interest is the difference between the binding free energy of a mutant peptide to a micelle and that of the wild-type IsCT to the same micelle, ΔΔGA·M,

ΔΔGA·M=ΔGbindA·MΔGbindISCT·M
(8)

where ΔGbindISCT·M and ΔGbindA·M are the binding free energy of the wild-type IsCT and one of its analogues to one of the micelles, respectively.

Although it is difficult to accurately compute absolute binding free energies, the rankings of the relative binding free energy, given by Equation (7), are expected to be accurate since it is presumed that systematic errors may cancel out [41,47,48]. For example, systematic errors include errors that arise due to the rotation or translation of a given molecule with respect to the grid coordinates when the PB equation is numerically solved [49], and these are expected to have approximately the same value for all systems studied.

In the following sections we describe in detail the relative binding free energy computation protocol.

2.4 Solvation free energy calculations

The relative binding affinities were calculated from two separate simulations of each peptide and the corresponding complex. This will help to account for the flexibility of the peptides and other conformational changes upon binding. Snapshots were extracted every 20 ps intervals from the last 10 ns of the trajectories of peptide–DPC complexes and the last 5 ns of the trajectories of peptide–SDS complexes. For the simulations of the peptides in water, snapshots were extracted every 10 ps from the last 3 – 4 ns of the trajectories of individual peptides. The solvent accessible surface areas were calculated using a probe radius of 1.4 Å. No cut-off distance was applied in evaluating the non-bonded and electrostatic parts of the molecular mechanics energy.

The electrostatic contributions to the solvation and desolvation free energies were obtained by solving the linearised Poisson–Boltzmann (PB) equation on a 1713 grid and using the PBEQ module [50] as implemented in the CHARMM program. In the PB model, the system consisting of a macromolecule (e.g. peptide or peptide–micelle complex), water molecules and salt ions can be represented by a dielectric cavity with explicit fixed charges on the macromolecule, a high dielectric constant and continuous mobile charges, respectively. The atomic coordinates and vdW radii of the peptides and micelles were used to define the dielectric cavity region as the region inaccessible to contact by a sphere of radius 1.4 Å rolling over the molecular surface [41]. The ion exclusion radius was set to 1.5 Å and the ionic strength was set to 0.15 M in the solvated state and to 0.0 in the reference gas phase state. The electrostatic potential was computed using a hexahedral mesh with a grid space of 0.5 Å and the Debye–Hückel approximation to estimate the boundary conditions [51].

There have been multiple values reported in the literature of the dielectric constant for proteins and lipid bilayers [41,5255]. Since the charge distribution of the peptides, micelles and peptide–micelle complexes are accounted for explicitly by averaging over many configurations, all were assigned a low dielectric constant of 1.0 [56,57], whereas the solvent was assigned a dielectric constant of 80. It is worth mentioning here that a dielectric constant of 1.0 for the peptides and peptide–micelle complexes does not account for the induced polarisability of these molecules and molecular complexes.

2.5 Entropic analysis

The loss of translational, rotational, conformational and vibrational entropies upon complexation for a peptide–micelle system is a particularly challenging problem due to the flexibility of the peptides and the fluid nature of the micelles, both of which lead to relatively large fluctuations of the conformations of the bound peptide around its equilibrium state. Therefore, an extensive sampling of all degrees of freedom is required to accurately compute the association entropy. In this work and for each snapshot, the bound peptide–micelle complex and the micelle itself are considered rigid molecules in the gas phase. This allows utilising polyatomic ideal gas formulas to compute the translational and rotational entropies. The conformational entropy loss upon complexation for the peptides was approximated by time average of the loss of side-chain conformational entropy which was estimated using the empirical scale of Pickett and Sternberg [44] for each snapshot in the trajectories. A more detailed analysis for computing the change in the translational and rotational entropy upon complexation based on the principal RMS fluctuations in the distance of centre of masses and Euler angles used to describe the orientation of the peptide with respect to the micelles [58], is not expected to change the rankings of the final relative binding energy values.

The change in the vibrational entropy upon complexation is usually calculated in the gas phase or with a distance dependent dielectric constant because it is generally difficult to account for the loss of entropy for the solvent from simulations [57,59,60]. The VIBRAN module of the CHARMM program was used to compute the force constant matrix to determine the normal mode frequencies. Table 2 shows that there are significant differences between the vibrational entropies obtained using a distance dependant and constant dielectric constant after minimisation of structures extracted from the last snapshot as well as the vibrational entropy computed using the time average of all snapshots using a constant dielectric and without performing minimisation. The difference between the vibrational entropy obtained using normal mode analysis with a constant dielectric and minimised peptide structures, and vibrational entropy computed using the time average of all snapshots without performing minimisation is relatively small. Although the values of the vibrational entropy are likely to be underestimated, we believe that adopting the time average values of vibrational entropies in the gas phase for our relative binding free energy calculations is a reasonable approximation. We evaluated the vibrational entropies of the peptide–micelle complexes by averaging the values of all snapshots collected every 100 ps from the peptide–micelle simulation trajectories, whereas the vibrational entropies of the peptides where evaluated by averaging the values of all snapshots collected every 10 ps from the peptide simulation trajectories.

Table 2
Vibrational Entropy values of IsCT and its analogues obtained using normal mode analysis.

2.6 Stability of simulated systems

In the previous work, we have demonstrated that peptide–micelle systems reach equilibrium when experimental information for the peptide structures is used to start the simulations [14,61,62]. In the simulations of the SDS complexes, the three peptides diffuse from the core of the SDS micelle to water–micelle interface in 4–6 ns (Figure 3(a)). The positions of the centre of mass of the three peptides with respect to the centre of mass of the micelles have reached equilibrium after about 30 ns. Figure 4(a)–(c) indicates that the root mean squared deviations (RMSDs) from the initial minimised structures for the Cα atoms of the three peptides in the presence of water, SDS or DPC environments are typical with fast initial rise. The RMSDs of the three peptides in the SDS system relatively stabilised after approximately 35 ns. Figure 3(b) indicates that the position of the centre of mass of three peptides in the DPC complexes have reached equilibrium after about 35 ns. IsCT takes more than 30 ns to reach its equilibrium position with respect to the centre mass of the DPC micelle possibly due to the weaker electrostatic interactions of IsCT with the DPC micelle. Furthermore, IsCT in DPC complex is very flexible. It undergoes fast conformational changes especially in the first 30 ns of the simulation as indicated by the results in Figure 4(c). Equilibrium is assumed to have been reached after about 40 ns for the simulations of the three peptides in DPC complexes based on the stable distance between the centre of mass of the peptides and the micelles, number of backbone dihedral angle clusters and the RMSDs of the peptides. Similarly, we also assume that the system is at equilibrium after about 17 ns for the simulations of the three peptides in water (see Figure 4(a)). Table 3 shows the electrostatic and van der Waals interaction energies between peptides and the micelles which are calculated using the last 10 ns of the simulations of DPC complexes and the last 5 ns of the SDS complexes.

Figure 3
(a) Time evolution of the positions of the centre of mass of the three peptides with respect to the centre of mass of the SDS micelle. (b) Time evolution of the positions of the centre of mass of the three peptides with respect to the centre of mass of ...
Figure 4
Time evolution of the RMSD from the initial minimised structures for the Cα atoms of the three peptides in the presence of (a) water, (b) SDS and (c) DPC environments, respectively.
Table 3
Electrostatic and van der Waal interaction energies of IsCT and its analogues to SDS and DPC micelles.

Although the length of our simulations is significantly longer than many of the reported simulations of peptide–micelle systems [6366], it is difficult to determine the simulation time required to sample the conformational space of these peptides and hence obtain reliable results for binding free energy calculations. To establish the reliability of the relative binding free energy results and to further investigate if the peptides have reached equilibrium, we utilised the ART-2 clustering algorithm [6769] as implemented in CHARMM. ART-2 is a statistical technique of conformational clustering. Generated conformations are grouped together defining a measure of conformational similarity based on dihedral angle distances. The generated clusters can be thought of as ‘conformational states’. The smaller the number, the more stable the peptide conformationally. ART-2 helps determine the simulation production periods by estimating the number of available conformations at different time intervals based on the number of clusters of the dihedral angles of the backbone of the three peptides found in a given interval. Time series of each peptide dihedral angles were extracted from the trajectories starting at different initial times, tini, and ending at the final time of each trajectory. tini was sampled at equal time intervals of 1 ns for peptides in SDS and DPC environments and 200 ps for peptides in water. The results for this clustering analysis for the three peptides in SDS and using a cluster radius of 35° are shown in Figure 5. Both IsCT and IsCT2 form one cluster beyond tini = 35 ns. However, IsCT3 forms two clusters which are close to each other except for two dihedral angles, ψ2 and [var phi]3. The ψ2 and [var phi]3 values are equal to (−101.79, −178.57), (53.74, −107.93) for the first and second cluster, respectively. The percentage of the number of members in the first cluster to the total number of structures in the production period is approximately 86%. Similar behaviour was observed for IsCT3 in DPC with two clusters centred at three different dihedral angles; ψ2, [var phi]3 and ψ3, which are equal to (−69.91, 122.91, 175.150) and (148.60, −163.80, −149.81) for the first and second cluster, respectively. The percentage of the number of members in the first cluster to the total number of structures in the production period is approximately 96%. This means that the average values of the various contributions to the relative binding free energy of IsCT3 to DPC micelle are dominated by one conformation. The stabilisation of the number of the clusters suggests that equilibrium has been reached, where many states with high energy have the same free energy as few states with low energy. Therefore, it is reasonable to start the production periods at 17, 35 and 40 ns for the unbound peptides, the SDS and DPC complexes, respectively.

Figure 5
The number of clusters of each peptide dihedral angles found in the interval tini < t < tfinal.

In Figure 6, the structures of the peptides are shown at the initial minimised structures and the final structures. IsCT3 and IsCT2, remain helical, whereas IsCT loses helicity in the presence of DPC environment. In SDS complexes, IsCT and IsCT2 are helical, while, IsCT3 has a bend in the middle followed by an α-helix.

Figure 6
A comparison between the initial minimised structures of the three peptides (top) and the final structures (bottom) as extracted from the last snapshot of each simulation in water, SDS and DPC environments, respectively.

2.7 Micelles as membrane mimics

There has been a long standing interest of determining the validity of anionic and zwitterionic micelles as mimics of bacterial and mammalian membranes. Although the correlation between the binding free energy of AMPs to model micelles and the binding free energy of AMPs to phospholipid membranes is not rigorously understood, SDS and DPC micelles have been widely used in NMR spectroscopic studies of AMPs to resemble zwitterionic mammalian lipid and anionic bacterial membrane interfaces, respectively [7074]. The main advantages of using micelles over lipid bilayers are the faster relaxation time scales of DPC and SDS molecules and the smaller domains of simulation. Moreover, both micelles have many of the essential physicochemical properties needed to modulate peptide–membrane interactions [14,15,61,62,75]. Ultimately, the specificity of AMPs towards bacterial cells is believed to arise from the differences in the electrostatic nature of the surface of bacterial (anionic) and mammalian (zwitterionic) interfaces. These differences in the modes of interactions of AMPs with anionic and zwitterionic interfaces can be captured in their differential modes of interaction with SDS and DPC micelles. A particular shortcoming of using micelles to study AMPs is their spherical geometry, which is fundamentally different than planar geometry of lipid bilayers. The interfacial curvature can potentially induce peptide structures which may never be observed near lipid bilayers. However, experimental structures of peptides formed upon association with detergent micelles are sometimes found to be comparable to the experimental structures of the same peptides bound to phospholipid bilayers, which more accurately resemble the biological membrane [76].

3. Results and discussion

To investigate the relationship between relative binding affinities of IsCT analogues to the surfactant micelles and their biological function (activity and hemolycity), the various contributions to the binding free energy were evaluated using the sampling windows discussed earlier; Table 4 shows the values of these contributions for IsCT2 and IsCT3. The haemolysis percentage measurements of IsCT, IsCT2 and IsCT3 were found to be 78, 52 and 20%, respectively [27]. The relative binding energy of IsCT2 to the DPC micelle is found to be equal to (8 kcal mol−1), while the relative binding affinity of IsCT3 is found to be equal to (65 kcal mol−1). Both values correlate with the ranking of their haemolytic toxicities relative to the haemolytic toxicity of the wild-type IsCT. The positive values of the relative binding free energies indicate that both peptides are indeed correlated with the relative haemolytic toxicity of these peptides.

Table 4
Components of the relative binding free energy of IsCT2·DPC and IsCT3·DPC.

The evaluation of the relative binding free energies of IsCT analogues to SDS micelle was initially more problematic. On the one hand, the relative binding free energy of IsCT2 was (−13 kcal mol−1), a stronger binding which correlates with the higher antimicrobial activity of IsCT2 relative to that of IsCT. On the other hand, for the IsCT3 analogue, as shown in Table 6, the computed binding free energy difference of 64 kcal mol−1 does not correlate with its measured antimicrobial activities. What may explain this nonetheless is the fact that one of the SDS molecules dissociated from the peptide–micelle complex after 33.6 ns. The aggregation number of the SDS micelle hence changed from 60 to 59. Therefore, an immediate comparison between simulated peptide–micelle complexes is not appropriate.

Table 6
Components of the relative binding free energy of IsCT2·SDS and IsCT3·SDS.

To derive a reasonable estimate of the binding free energy of IsCT3 to SDS60 micelle, we use the thermodynamic cycle shown in Figure 7. The binding free energy of an SDS molecule to an IsCT3·SDS59 complex, ΔG1, was approximated using the binding free energy of an SDS monomer to SDS59 micelle, which is assumed to be comparable to the binding free energy of an SDS monomer to IsCT3·SDS59 complex. It was evaluated using the average effective binding free energies per SDS molecule in an SDS60 micelle (Table 5). The values for the bound SDS molecule are computed using a 20 ns MD simulation of an SDS60 micelle with the same simulation parameters and conditions of a peptide–SDS complex. The same values were also computed for the dissociated SDS monomer extracted from the IsCT3·SDS60 complex trajectory. ΔG2 is computed using a thermodynamic cycle similar to the thermodynamic cycle used in Figure 2. It is worth mentioning here that the various energetic and entropic values for the free SDS monomer are not needed to compute the binding free energy of IsCT3 to SDS60 through the path in which SDS60 loses a monomer and then binds to IsCT3 to form IsCT3·SDS59 complex and a free SDS monomer, then followed by the SDS monomer binding to IsCT3·SDS59 to form the IsCT3·SDS60 complex. The reason for this is that these contributions will cancel out when adding up ΔG1 and ΔG2. The net relative binding affinity of IsCT3 to SDS micelle and the individual contributions as estimated using the aforementioned assumption are summarised in Table 6. The relative binding free energy value is consistent with the fact that IsCT3 is more antimicrobial active than IsCT. Finally, we report the values of the standard binding free energies of IsCT and its analogues as well as the various contributions to them as shown in Table 7. These results indicate that the change in the molecular mechanical energies and the non-polar contributions to the solvation energy favour the binding, while the changes in the electrostatic contributions to the solvation energy are not favourable to the binding. The calculations also suggest that the net association entropy values have a small contribution to relative binding affinities. In particular, the loss of translational, rotational and side-chain conformational entropy was found to oppose binding to SDS micelle, while the vibrational entropy change was found to favour binding in the three peptide complexes. The binding free energy differences are largely determined by the balance between the solvation free energy differences and interaction energy differences. The binding affinity values reported in this study are relatively large, indicating a strong binding affinity. The relative binding affinities are also large in both SDS and DPC micelles. Such large differences in the relative binding free energy are not surprising since comparable differences have been observed in flexible peptide–protein complexes [54,77].

Figure 7
A thermodynamic cycle used to recalculate the relative binding free energy of IsCT analogue, IsCT3.
Table 5
Components of the standard free energy of SDS60 are shown in the first column. The second column indicates the average of each component per SDS monomer, while the third column shows these values for a free SDS monomer.
Table 7
Components of the standard binding free energy of IsCT·SDS, IsCT2·SDS and IsCT3·SDS.

4. Conclusions

The objectives of this study were to develop a reliable method for computing the relative binding free energy of AMPs to surfactant micelles and to determine whether computed relative binding affinities correlate with the relative antimicrobial activities as well as the haemolytic toxicities of AMPs. Although the calculations relied on a number of assumptions, we showed that the relative binding free energies can be estimated using the MM-PBSA approach in a way that they correlate with biological function. It is important to point out that we did not fit any model parameters to obtain the experimental rankings of antimicrobial activity and cytotoxicity of the IsCT analogues. We have applied the standard MM-PBSA framework that is computationally tractable, although still resource demanding, since more than 330 ns of molecular dynamics simulation was carried out. The correlation of rankings of the relative binding affinity values with the experimental antimicrobial activity and cytotoxicity data is remarkable. Certainly, the study of only three peptides with two simplified membrane mimics does not provide conclusive correlation evidence, with adequate statistical significance. Because structural information of peptides is largely lacking, the presented protocol cannot be currently employed with a large AMP dataset, which would be ideal. The results are encouraging, nonetheless, and point to the importance of the initial binding events between a peptide and a membrane on the tendency of the peptide to disrupt the membrane’s integrity.

Acknowledgments

This work was supported by a grant from the NIH (GM 070989). The authors gratefully acknowledge the computational support from the Minnesota Supercomputing Institute. This work also was partially supported by the National Computational Science Alliance under TG-MCA04N033. A. Sayyed-Ahmad also appreciates the partial funding from the Consortium for Bioinformatics and Computational Biology, University of Minnesota, Minneapolis.

References

1. Yoshikawa TT. Antimicrobial resistance and aging: beginning of the end of the antibiotic era? J Am Geriatr Soc. 2002;50:226–229. [PubMed]
2. Hiramatsu K, Hanaki H, Ino T, Yabuta K, Oguri T, Tenover FC. Methicillin-resistant Staphylococcus aureus clinical strain with reduced vancomycin susceptibility. J Antimicrob Chemother. 1997;40:135–136. [PubMed]
3. Kam KM, Luey KY, Fung SM, Yiu PP, Harden TJ, Cheung MM. Emergence of multiple-antibiotic-resistant Streptococcus pneumoniae in Hong Kong. Antimicrob Agents Chemother. 1995;39:2667–2670. [PMC free article] [PubMed]
4. Pallares R, Linares J, Vadillo M, Cabellos C, Manresa F, Viladrich PF, Martin R, Gudiol F. Resistance to penicillin and cephalosporin and mortality from severe pneumococcal pneumonia in Barcelona, Spain. N Engl J Med. 1995;333:474–480. [PubMed]
5. Grubb W. Genetics of MRSA. Rev Med Microbiol. 1998;9:153–162.
6. Orsolya T. Antimicrobial peptides: new candidates in the fight against bacterial infections. Peptide Sci. 2005;80:717–735. [PubMed]
7. Zasloff M. Antimicrobial peptides of multicellular organisms. Nature. 2002;415:389–395. [PubMed]
8. Langham AA, Khandelia H, Schuster B, Waring AJ, Lehrer RI, Kaznessis YN. Correlation between simulated physicochemical properties and hemolycity of protegrin-like antimicrobial peptides: predicting experimental toxicity. Peptides. 2008;29:1085–1093. [PMC free article] [PubMed]
9. Bolintineanu DS, Langham AA, Davis HT, Kaznessis YN. Molecular dynamics simulations of three protegrin-type antimicrobial peptides: interplay between charges at the termini, Î2-sheet structure and amphiphilic interactions. Mol Simulation. 2007;33:809–819. [PMC free article] [PubMed]
10. Hong RW, Shchepetov M, Weiser JN, Axelsen PH. Transcriptional profile of the Escherichia coli response to the antimicrobial insect peptide cecropin A. Antimicrob Agents Chemother. 2003;47:1–6. [PMC free article] [PubMed]
11. Gidalevitz D, Ishitsuka Y, Muresan AS, Konovalov O, Waring AJ, Lehrer RI, Lee KYC. Interaction of antimicrobial peptide protegrin with biomembranes. PNAS. 2003;100:6302–6307. [PubMed]
12. Yechiel S. Mode of action of membrane active antimicrobial peptides. Peptide Sci. 2002;66:236–248. [PubMed]
13. Langham AA, Kaznessis YN. Effects of mutations on the C-terminus of protegrin-1: a molecular dynamics simulation study. Mol Simulation. 2006;32:193–201.
14. Langham AA, Khandelia H, Kaznessis YN. How can a Beta-sheet peptide be both a potent antimicrobial and harmfully toxic: molecular dynamics simulations of protegrin-1 in micelles. J Peptide Sci. 2006;84:219–231. [PubMed]
15. Khandelia H, Langham AA, Kaznessis YN. Driving engineering of novel antimicrobial peptides from simulations of peptide-micelle interactions. Biochim Biophys Acta (BBA) Biomembranes. 2006;1758:1224–1234. [PMC free article] [PubMed]
16. Peschel A. How do bacteria resist human antimicrobial peptides? Trends Microbiol. 2002;10:179–186. [PubMed]
17. Hancock REW, Chapple DS. Peptide antibiotics. Antimicrob Agents Chemother. 1999;43:1317–1323. [PMC free article] [PubMed]
18. Qu XD, Harwig SS, Oren AM, Shafer WM, Lehrer RI. Susceptibility of Neisseria gonorrhoeae to protegrins. Infect Immun. 1996;64:1240–1245. [PMC free article] [PubMed]
19. Langham A, Waring A, Kaznessis Y. Comparison of interactions between beta-hairpin decapeptides and SDS/DPC micelles from experimental and simulation data. BMC Biochem. 2007;8:11. [PMC free article] [PubMed]
20. Langham AA, Ahmad AS, Kaznessis YN. On the nature of antimicrobial activity: a model for protegrin-1 pores. J Am Chem Soc. 2008;130:4338–4346. [PMC free article] [PubMed]
21. Sayyed-Ahmad A, Kaznessis YN. Determining the orientation of protegrin-1 in DLPC bilayers using an implicit solvent-membrane model in PLoS one. Public Libr Sci. 2009;4:e4799. [PMC free article] [PubMed]
22. Bolintineanu DS, Sayyed-Ahmad A, Davis HT, Kaznessis YN. Poisson-Nernst-Planck models of nonequilibrium ion electrodiffusion through a protegrin transmembrane pore. PLoS Comput Biol. 2009;5:e1000277. [PMC free article] [PubMed]
23. Yorek MA. Biological distribution. In: Cevc G, editor. Phospholipids Handbook. Marcel Dekker, Inc; New York, NY: 1993.
24. Nelson D, Lehninger A, Cox M. Principle of Biochemistry. Worth Publishing; New York, NY: 2000.
25. Yeaman MR, Yount NY. Mechanisms of antimicrobial peptide action and resistance. Pharmacol Rev. 2003;55:27–55. [PubMed]
26. Dai L, Yasuda A, Naoki H, Corzo G, Andriantsiferana M, Nakajima T. IsCT, a novel cytotoxic linear peptide from Scorpion Opisthacanthus madagascariensis. Biochem Biophys Res Commun. 2001;286:820–825. [PubMed]
27. Lee K, Shin SY, Kim K, Lim SS, Hahm KS. Antibiotic activity and structural analysis of the scorpion-derived antimicrobial peptide IsCT and its analogs. Biochem Biophys Res Commun. 2004;323:712–719. [PubMed]
28. Chong LT, Duan Y, Wang L, Massova I, Kollman PA. Molecular dynamics and free-energy calculations applied to affinity maturation in antibody 48G7. PNAS. 1999;96:14330–14335. [PubMed]
29. Wymore T, Gao XF, Wong TC. Molecular dynamics simulation of the structure and dynamics of a dodecylphosphocholine micelle in aqueous solution. J Mol Struct. 1999;195:195–210.
30. Wong T, Kamath S. Molecular dynamics simulations of Adrenocorticotropin (1–24) peptide in a solvated dodecylphosphocholine (DPC) micelle and in a Dimyistoylphosphatidylcholine (DMPC) bilayer. J Biomol Struct Dyn. 2002;20:39–57. [PubMed]
31. MacKerell A. Molecular Dynamics Simulation Analysis of a sodium dodecyl sulfate micelle in aqueous solution: decreased fluidity of the micelle hydrocarbon interior. J Phys Chem. 1995;99:1846–1855.
32. Brooks BR, Bruccoleri RE, Olafson BD, States DJ, Swaminathan S, Karplus M. A program for macromolecular energy, minimization, and dynamics calculations. J Comput Chem. 1983;4:187–217.
33. MacKerell AD, Bashford D, Bellott M, Dunbrack RL, Evanseck JD, Field MJ, Fischer S, Gao J, Guo H, Ha S, et al. Allatom empirical potential for molecular modeling and dynamics studies of proteins. J Phys Chem B. 1998;102:3586–3616. [PubMed]
34. Feller SE, MacKerell JAD. An improved empirical potential energy function for molecular simulations of phospholipids. J Phys Chem B. 2000;104:7510–7515.
35. Jorgensen WL, Chandrasekhar J, Medura JD, Impey RW, Klein ML. Comparison of simple potential function for simulating liquid water. J Chem Phys. 1983;79:926–935.
36. Andersen H. Molecular dynamics simulations at constant pressure and/or temperature. J Chem Phys. 1980;72:2384–2393.
37. Hoover W. Canonical dynamics: equilibrium phase-space distributions. Phys Rev A. 1985;31:1695–1697. [PubMed]
38. Darden T, York D, Pedersen L. Particle mesh Ewald: an N log(N) methods for Ewald sums in large systems. J Chem Phys. 1993;98:100089–110092.
39. Essmann U, Perera L, Berkowitz M, Darden T. A smooth particle mesh Ewald method. J Chem Phys. 1995;103:8577–8593.
40. Straatsma T, McCammon J. Multiconfiguration thermodynamic integration. J Chem Phys. 1991;95:1175–1188.
41. Noskov SY, Lim C. Free energy decomposition of protein-protein interactions. Biophys J. 2001;81:737–750. [PubMed]
42. Tidor B, Karplus M. The contribution of vibrational entropy to molecular association: the dimerization of insulin. J Mol Biol. 1994;238:405–414. [PubMed]
43. McQuarrie DA. Statistical Mechanics. Harper Collins; New York, NY: 1976.
44. Pickett SD, Sternberg MJE. Empirical scale of side-chain conformational entropy in protein folding. J Mol Biol. 1993;231:825–839. [PubMed]
45. Honig B, Nicholls A. Classical electrostatics in biology and chemistry. Science. 1995;268:1144–1149. [PubMed]
46. Chandler D. Interfaces and the driving force of hydrophobic assembly. Nature. 2005;437:640–647. [PubMed]
47. Pathiaseril A, Woods RJ. Relative energies of binding for antibody-carbohydrate-antigen complexes computed from free-energy simulations. J Am Chem Soc. 2000;122:331–338. [PMC free article] [PubMed]
48. Volkhard Helms RCW. Free energies of hydration from thermodynamic integration: comparison of molecular mechanics force fields and evaluation of calculation accuracy. J Comput Chem. 1997;18:449–462.
49. Feig M, Onufriev A, Lee MS, Im W, Case DA, Brooks CL., III Performance comparison of generalized Born and Poisson methods in the calculation of electrostatic solvation energies for protein structures. J Comput Chem. 2004;25:265–284. [PubMed]
50. Im W, Beglov D, Roux B. Continuum solvation model: computation of electrostatic forces from numerical solutions to the Poisson-Boltzmann equation. Comput Phys Commun. 1998;111:59–75.
51. Sayyed-Ahmad A, Tuncay K, Ortoleva PJ. Efficient solution technique for solving the Poisson-Boltzmann equation. J Comput Chem. 2004;25:1068–1074. [PubMed]
52. Mottamal M, Zhang J, Lazaridis T. Energetics of the native and non-native states of the glycophorin transmembrane helix dimer. Proteins. 2006;62:996–1009. [PubMed]
53. Laitinen T, Kankare JA, Peräkylä M. Free energy simulations and MM-PBSA analyses on the affinity and specificity of steroid binding to antiestradiol antibody. Proteins. 2004;55:34–43. [PubMed]
54. Ganoth A, Friedman R, Nachliel E, Gutman M. A molecular dynamics study and free energy analysis of complexes between the Mlc1p protein and two IQ motif peptides. Biophys J. 2006;91:2436–2450. [PubMed]
55. Mark LC, Olson A. Molecular docking of superantigens with class II major histocompatibility complex proteins. J Mol Recognit. 1997;10:277–289. [PubMed]
56. Warshel A, Sharma PK, Kato M, Parson WW. Modeling electrostatic effects in proteins. Biochim Biophys Acta (BBA) Proteins Proteomics. 2006;1764:1647–1676. [PubMed]
57. Holger Gohlke DAC. Converging free energy estimates: MM-PB(GB)SA studies on the protein-protein complex Ras-Raf. J Comput Chem. 2004;25:238–250. [PubMed]
58. Carlsson J, Aqvist J. Absolute and relative entropies from computer simulation with applications to ligand binding. J Phys Chem B. 2005;109:6448–6456. [PubMed]
59. Lee MS, Olson MA. Calculation of absolute protein-ligand binding affinity using path and endpoint approaches. Biophys J. 2006;90:864–877. [PubMed]
60. Zoete V, Meuwly M, Karplus M. Study of the insulin dimerization: binding free energy calculations and per-residue free energy decomposition. Proteins: Struct Function Bioinf. 2005;61:79–93. [PubMed]
61. Khandelia H, Kaznessis YN. Molecular dynamics simulations of helical antimicrobial peptides in SDS micelles: what do point mutations achieve? Peptides. 2005;26:2037–2049. [PubMed]
62. Khandelia H, Kaznessis YN. Molecular dynamics simulations of the helical antimicrobial peptide ovispirin-1 in a zwitterionic dodecylphosphocholine micelle: insights into host-cell toxicity. J Phys Chem B. 2005;109:12990–12996. [PubMed]
63. Lague P, Roux B, Pastor RW. Molecular dynamics simulations of the influenza hemagglutinin fusion peptide in micelles and bilayers: conformational analysis of peptide and lipids. J Mol Biol. 2005;354:1129–1141. [PubMed]
64. Xinfeng Gao TCW. Molecular dynamics simulation of adrenocorticotropin (1–10) peptide in a solvated dodecylphosphocholine micelle. Biopolymers. 2001;58:643–659. [PubMed]
65. Wymore T, Wong TC. Molecular dynamics study of substance P peptides partitioned in a sodium dodecylsulfate micelle. Biophys J. 1999;76:1213–1227. [PubMed]
66. Bond PJ, Sansom MSP. Membrane protein dynamics versus environment: simulations of OmpA in a micelle and in a bilayer. J Mol Biol. 2003;329:1035–1053. [PubMed]
67. Carpenter GA, Grossberg S. Self organization of stable category recognition codes for analog input patterns. Appl Opt. 1987;26:4919–4930. [PubMed]
68. Karpen ME, Tobias DJ, Brooks CL., III Statistical clustering techniques for the analysis of long molecular dynamics trajectories: analysis of 2.2-ns trajectories of YPGDV. Biochemistry. 1993;32:412–420. [PubMed]
69. Pao Y. Adaptive Pattern Recognition and Neural Networks. Addison-Wesley; New York, NY: 1989.
70. Schibli DJ, Hwang PM, Vogel HJ. Structure of the antimicrobial peptide tritrpticin bound to micelles: a distinct membrane-bound peptide fold. Biochemistry. 1999;38:16749–16755. [PubMed]
71. Jing W, Hunter HN, Hagel J, Vogel HJ. The structure of the antimicrobial peptide Ac-RRWWRF-NH2 bound to micelles and its interactions with phospholipid bilayers. J Pept Res. 2003;61:219–229. [PubMed]
72. Prates MV, Sforca ML, Regis WCB, Leite JRSA, Silva LP, Pertinhez TA, Araujo ALT, Azevedo RB, Spisni A, Bloch C., Jr The NMR-derived solution structure of a new cationic antimicrobial peptide from the skin secretion of the Anuran Hyla punctata. J Biol Chem. 2004;279:13018–13026. [PubMed]
73. Lindberg M, Biverstahl H, Graslund A, Maler L. Structure and positioning comparison of two variants of penetratin in two different membrane mimicking systems by NMR. Eur J Biochem. 2003;270:3055–3063. [PubMed]
74. Shankaramma SC, Athanassiou Z, Zerbe O, Moehle K, Mouton C, Bernardini F, Vrijbloed JW, Obrecht D, Robinson JA. Macrocyclic hairpin mimetics of the cationic antimicrobial peptide protegrin I: a new family of broad-spectrum antibiotics. Chem Bio Chem. 2002;3:1126–1133. [PubMed]
75. MacKerrell A. Molecular dynamics simulation analysis of a sodium dodecyl sulfate micelle in aqueous solution: decreased fluidity of the micelle hydrocarbon interior. J Phys Chem. 1995;99:327–341.
76. Rozek A, Friedrich CL, Hancock REW. Structure of the bovine antimicrobial peptide indolicidin bound to dodecylphosphocholine and sodium dodecyl sulfate micelles. Biochemistry. 2000;39:15765–15774. [PubMed]
77. Archontis G, Simonson T, Karplus M. Binding free energies and free energy components from molecular dynamics and Poisson-Boltzmann calculations. Application to amino acid recognition by aspartyl-tRNA synthetase. J Mol Biol. 2001;306:307–327. [PubMed]