|Home | About | Journals | Submit | Contact Us | Français|
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.
The number of reports of drug-resistant strains of bacteria that cause infectious diseases such as pneumonia, tuberculosis, syphilis and gonorrhoea is rapidly rising [1–5]. 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[8–22].
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 , 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  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) . 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 . Lee et al.  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 . 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 . 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.
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  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.
The starting coordinates of the DPC micelle–water complex were extracted from simulations carried out by Wymore et al. . 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 . 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 .
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 . 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.
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  version c30b2 and the all hydrogen force field CHARMM27 [33,34]. Water was modelled using the TIP3P force potential . 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 . 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  and the Hoover temperature control , 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.
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 .
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 . 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.
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 PMaq, is
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:
The 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 ) and given by
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:
where the solvation free energy difference, ΔΔGsolv, is given by
The gas phase binding free energy, ΔGgas, is computed using the following formula,
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 , while the loss of conformational entropy of each peptide side chain upon complexation is estimated using the empirical scale of Pickett and Sternberg . The electrostatic contributions to the solvation free energies of each complex component and the complex itself are calculated using the Poisson–Boltzmann method . 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:
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 .
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,
where and 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 , 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.
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  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 . 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 .
There have been multiple values reported in the literature of the dielectric constant for proteins and lipid bilayers [41,52–55]. 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.
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  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 , 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.
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.
Although the length of our simulations is significantly longer than many of the reported simulations of peptide–micelle systems [63–66], 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 [67–69] 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 3. The ψ2 and 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, 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.
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.
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 [70–74]. 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 .
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 . 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.
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.
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].
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.
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.