Article pubs.acs.org/JCTC
Terms of Use CC-BY
Energy Landscapes and Global Optimization of Self-Assembling Cyclic Peptides Mark T. Oakley and Roy L. Johnston* School of Chemistry, University of Birmingham, Edgbaston, Birmingham, B15 2TT, United Kingdom ABSTRACT: Self-assembled cyclic peptide nanotubes have attracted much attention because of their antimicrobial properties. Here, we present calculations on the formation of cyclic peptide dimers using basin-hopping and discrete path sampling. We present an analysis of the basin-hopping move sets that most efficiently explore the conformations of cyclic peptides. Group rotation moves, in which sections of the ring are rotated as a rigid body, are the most effective for cyclic peptides containing up to 20 residues. For cyclic peptide dimers, we find that a combination of group rotation intramolecular moves and rigid body intermolecular moves performs well. Discrete path sampling calculations on the cyclic peptide dimers show significant differences in the dimerization of hexa- and octapeptides. nanotubes have been shown to form across22,23 or along24 the membrane, depending on the peptide sequence. Here, we take a different approach and consider only the stationary points on the energy landscape. The solvent is treated implicitly, which substantially reduces the computational resources required. Locating minima on the energy landscape then becomes a global optimization problem that is split into two parts: searching the conformations of individual cyclic peptides and assembly of cyclic peptide aggregates. Several methods are available to search the conformations of small flexible molecules. Those that have been used on cyclic peptides include distance geometry methods,25 Monte Carlo,26 LowMode molecular dynamics,27 iterative stochastic search,28 dihedral angle sampling29 and diffusion equation evolutionary programming simulated annealing.30 Optimization of clusters of flexible molecules is a more difficult problem. An evolutionary algorithm has been used to optimize the dimer of kanamycin A.31,32 Basin-hopping has been used to optimize the small aggregates of the amyloid-β peptide.33,34 Here, we use the basin-hopping algorithm35 to locate the minima on the cyclic peptide potential energy surfaces. To understand the conformations of molecules fully we need to know both the minima and the transition states linking them. Discrete path sampling36−38 is an efficient method for exploring energy landscapes and we have previously used it to study the conformations of peptides.39,40 Here, we present energy landscapes for the peptide dimers from discrete path sampling calculations.
1. INTRODUCTION Self-assembly is the process by which small molecules aggregate to form large, ordered structures. Some head-to-tail cyclic peptides are known to self-assemble to form nanotubes.1−5 Self-assembling cyclic peptide nanotubes have attracted attention for their possible biomedical uses, including activity against bacteria.6,7 Cyclic peptide nanotubes have been observed to form across lipid bilayers.3,8,9 The hydrogenbonding arrangement of cyclic peptide nanotubes is similar to that seen in amyloid fibrils, and it has been shown that cyclic peptides can disrupt amyloid formation.10 Cyclic peptide− polymer conjugates have potential uses as self-assembling materials.11,12 Self-assembling cyclic peptides have sequences comprising alternating D- and L-peptide residues because this arrangement places the side-chains equatorial to the peptide rings with the peptide groups aligned axially to hydrogen bond to other cyclic peptides. Cyclic octapeptides are the most widely studied because they have the highest tendency to form nanotubes.1−4 However, dimerization of cyclic hexapeptides has been observed,13 and nanotubes of larger cyclic peptides are known.9,14 Nanotubes with antiparallel hydrogen bonding arrangements are more stable than those with parallel hydrogen bonds.15 Several molecular dynamics (MD) investigations have been undertaken to model the properties of cyclic peptide nanotubes. The formation of dimers has been modeled with MD.16 The composition of cyclic peptides has a substantial effect on the stability and flexibility of their nanotubes.17 Longer nanotubes are more stable than shorter nanotubes.18,19 MD has been used to study the motion of water20 and other small molecules21 through cyclic peptide nanotubes. Cyclic peptide−polymer conjugates have also been studied with MD.12 All-atom MD is very computationally demanding, and modeling the formation of nanotubes is prohibitively expensive. For this reason, coarse-grained models have been used to study the aggregation of cyclic peptides in lipid bilayers, and © 2014 American Chemical Society
2. METHODS The energies of all structures were evaluated using the AMBER FF03 force field.41,42 The effect of solvent was included with the Generalized Born implicit solvation model.43 Calculations Received: January 3, 2014 Published: February 25, 2014 1810
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
were performed in water (εr = 80) or in a low-dielectric solvent (εr = 2) to model a membrane environment. Global optimization of single cyclic peptides and dimers was performed using the basin-hopping algorithm35 as implemented in GMIN.44 In this algorithm, structural deformations are followed by a gradient-driven quench to a local minimum before a Monte Carlo acceptance test. Thus, all attempted moves are between local minima on the potential energy surface. For efficient searching, the moves must be sufficiently large to escape the basin of attraction of the initial local minimum. However, moves that enter basins that are not physically accessible should be avoided. Inversion of the chirality of the Cα atoms is not possible in experiment and has a barrier of about 70 kcal mol−1 in the Amber force field.45 Therefore, any moves that lead to structures where the chirality of one or more Cα atoms are inverted are rejected. We consider three types of intramolecular move: atom displacements, MD moves, and group rotations. A series of basin-hopping optimizations were performed to benchmark each move class. For cyclic peptides comprising 10 or fewer peptide residues, the global minimum structures are found by all searches after fairly short basin-hopping runs. In these cases, the performance was measured in terms of the mean time to the first encounter of the global minimum from 100 independent searches. For larger cyclic peptides, the conformational spaces are too large for us to be confident that the lowest minima we find are the global minima. In these cases, basinhopping runs were run for 10 000 Monte Carlo steps and the performance was measured in terms of the lowest energies found at the end of each run. The initial structures for all of these searches were taken from a molecular dynamics simulation at 1600 K with snapshots taken every 10 ps. The parameters used for all move classes were optimized by performing sets of 10 searches with different values of the parameters for peptides containing up to 10 residues. These initial 10 searches are included in the set of 100 searches used to calculate the mean first encounter time. Unless noted, the same search parameters are used for all of the cyclic peptides studied here. The atom displacement moves (STEP keyword) involve displacement of all the Cartesian coordinates of all atoms by a up to 1.5 Å. In the MD moves (AMBERMDSTEPS keyword), short MD runs of 1000 steps are performed to escape the basins of attraction of local minima. These runs are performed at high temperature of 1600 K to allow the move to rapidly surmount large barriers, such as those between cis and trans isomers. In the group rotation moves (GROUPROTATION keyword), two atoms in the peptide backbone are selected and all atoms between them are rotated around the axis defined by those two atoms (Figure 1).46 Four types of group rotation are used (Table 1): cis−trans, peptide, dipeptide, and tripeptide moves. Cis−trans and peptide moves are attempted for all cyclic peptides. Dipeptide moves are attempted for cyclic peptides with six or more residues. The performance of tripeptide moves is evaluated for cyclic peptides comprising 12 or more residues. The magnitudes of all group rotation moves are in the range ±180°. The probability of selecting a group for group rotation is defined so that each basin-hopping step comprises on average 1/3 cis−trans moves, three peptide moves and, where appropriate, one of each of the larger moves. A combined set of intramolecular and intermolecular moves is needed to search the structures of dimers, and we consider two classes of move. MD moves use a short, high-temperature
Figure 1. Atoms involved in a dipeptide group rotation move. The axis is defined by two Cα atoms (yellow). All of the intermediate atoms (purple) are rotated around this axis. Some hydrogen atoms are omitted for clarity.
Table 1. Atoms Involved in the Each Type of Group Rotation Move move cis−trans peptide dipeptide tripeptide
ends Cα,i Cα,i Cα,i Cα,i
Ni+1 Cα,i+1 Cα,i+2 Cα,i+3
MD run to perform both the intra- and intermolecular moves. We also use a combination of group rotation intramolecular moves and rigid-body47 intermolecular moves. For both types of move, a hybrid minimization scheme is used (Figure 2). The second step is particularly important when using MD moves because the two cyclic peptides can drift apart during the MD run.
Figure 2. Basin-hopping scheme with hybrid minimization.
We evaluate the move sets on two types of cyclic peptides. Cyclic oligoglycines do not have any side chains and allow us to study the effect of the moves on the peptide backbones. We have not studied the dimerization of these peptides because they are not known to form nanotubes. Cyclic peptides with alternating D -Ala and L-Ala residues were selected as representative cyclic D-, L-peptides. Other peptide sequences form more stable nanotubes,17 but alanine was selected to reduce the effort of modeling the peptide side chains. The energy landscapes of some of the cyclic peptides and dimers were explored with discrete path sampling,36−38 as implemented in PATHSAMPLE.48 The databases of stationary points were seeded with minima from the basin-hopping calculations. Pathways between minima were found using OPTIM,49 with candidate transition states found by doubly nudged elastic band calculations50 and refined with hybrid eigenvector following.51−53 Pairs of minima were initially selected for connection by the missing connections algorithm54 and later using the UNTRAP method to remove artificial frustration.55 The resulting databases of stationary points are visualized as disconnectivity graphs.56,57 Stationary points that are permutation-inversion isomers are grouped together in 1811
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
second-lowest structure several times before locating the global minimum. The displacement moves require a similar amount of CPU time to MD moves when n ≤ 8. We have also tested the performance of basin-hopping on a series of cyclic peptides with alternating D-Ala and L-Ala residues. We have previously shown that the global minimum of cyclo-((D-Ala-L-Ala)2) includes one peptide in a cis conformation (Figure 4a).39 The global minima for all of the larger cyclo-((D-Ala-L-Ala)n) have all of their peptide bonds in trans conformations (Figure 4).
these disconnectivity graphs. All disconnectivity graphs show the free energy surface at 298 K. The free energies of all stationary points are calculated using the harmonic superposition approximation.58,59 We use the pattern of hydrogen bonds as a metric to measure the parallel or antiparallel alignment of these dimers. Hydrogen bonds in the dimers were counted using the definition from the Dictionary of Secondary Structure Prediction.60 Intermolecular hydrogen bonds were then assigned as parallel or antiparallel based on the angles between the two vectors joining the Cα atoms on either side of the peptide groups involved. The alignment is defined as npar−nantipar.
3. RESULTS 3.1. Intramolecular Moves. We have tested the performance of basin-hopping moves on a series of cyclic oligoglycines. For consistency with previous studies,39,40 all basin-hopping searches on single cyclic peptides were performed using water as a solvent. We have previously shown that the global minimum for cyclo-(Gly4) has all four peptide groups aligned axially, with no intramolecular hydrogen bonds (Figure 3a).39 The global minima for all of the larger cyclic oligoglycines exhibit several intramolecular hydrogen bonds (Figure 3).
Figure 4. Global minima of cyclo-((D-Ala-L-Ala)n) in water.
The performance of the displacement moves is poor for these molecules. Cartesian displacements that are sufficiently large to step between basins on the potential energy surface also lead to the inversion of at least one of the chiral centers in about 45% of the attempted moves. For comparison, MD moves never lead to inversion of chiral centers and group rotation moves invert in 5−10% of attempts. In most cases, the group rotation moves give the most efficient optimization. The MD moves locate the global minimum of cyclo-((D-Ala-LAla)4) (Figure 4c) very rapidly (Table 3). However, this seems to be a particularly easy global minimum to locate because the mean first encounter times with all move types are faster than the corresponding times for cyclo-((D-Ala-L-Ala)3). Basin-hopping optimization has also been performed on several larger cyclic peptides comprising up to 20 residues (table 4). For these molecules, we only consider the MD and group rotation moves because of the poor performance of the displacement moves in the previous calculations. The performance was assessed in terms of the mean energy of the lowest structure found in each basin-hopping run and the lowest energy found across all basin-hopping runs. For all of the cyclic peptides, the lowest minimum found using group rotation moves up to and including the dipeptide move is more stable than the lowest found with molecular dynamics moves. By both measures, group rotation with dipeptide moves out-performs MD moves for all cyclic oligoglycines. The inclusion of tripeptide moves makes only a small difference to the lowest energies for the smaller rings but gives a substantial improvement for cyclo-(Gly20). The CPU time taken for each Monte Carlo step is also substantially longer for the MD moves. The energy landscape of cyclo-((D-Ala-L-Ala)3) is funnelled, with small downhill barriers separating all local minima from the global minimum (Figure 5). The conformation with the lowest free energy in water is similar to that seen in cyclic peptide nanotubes. In nonpolar solvent, the lowest free energy structure exhibits some intramolecular hydrogen bonding. In cyclo-((D-Ala-L-Ala)3), conformers containing cis peptides are particularly unstable and are at least 8.8 kcal mol−1 above the global minimum in water and 10.3 kcal mol−1 higher in
Figure 3. Global minima of cyclo-(Glyn) in water.
For all four of the cyclic oligoglycines, the group rotation moves give the most efficient optimization in terms of the CPU time and the number of energy evaluations (Table 2). For n ≤ 8, the MD moves require a similar number of quenches to the group rotation moves but need more CPU time because of the overhead of performing the short MD run. When n = 10, molecular dynamics moves perform poorly and visit the Table 2. Mean Times to the First Encounter of the Global Minimum for 100 Basin-Hopping Runs on Cyclic Oligoglycines mean first encounter time molecule
move
cyclo-(Gly4)
displacement MD GR displacement MD GR displacement MD GR displacement MD GR
cyclo-(Gly6)
cyclo-(Gly8)
cyclo-(Gly10)
energy evaluations
minimizations
CPU time (s)
× × × × × × × × × × × ×
99 32 18 90 21 51 680 190 260 1000 1700 420
2.6 4.7 0.45 7.9 6.8 3.8 140 120 40 370 1900 120
1.8 3.7 3.2 2.6 2.7 1.3 3.0 2.9 8.8 5.2 3.0 1.8
104 104 103 104 104 104 105 105 104 105 106 105
1812
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
Table 3. Mean Times to the First Encounter of the Global Minimum in 100 Basin-Hopping Runs on Cyclic Oligo-D-,L-alanines mean first encounter time molecule cyclo-((D-Ala-L-Ala)2)
cyclo-((D-Ala-L-Ala)3)
cyclo-((D-Ala-L-Ala)4)
cyclo-((D-Ala-L-Ala)5)
move
energy evaluations
minimizations
CPU time (s)
displacement MD GR displacement MD GR displacement MD GR displacement MD GR
× × × × × × × × × × × ×
180 13 44 5200 1100 300 2600 39 130 12000 3900 530
11 3.7 2.7 1000 670 67 1100 48 63 8100 7700 430
4.5 1.6 1.1 1.9 1.5 1.2 1.2 6.2 6.9 6.4 6.6 3.3
4
10 104 104 106 106 105 106 104 104 106 106 105
Table 4. Lowest Energies Found in 20 × 10 000 Step Basin-Hopping Runs on Cyclic Peptides with 12 or More Residues molecule cyclo-(Gly12)
cyclo-(Gly14)
cyclo-(Gly16)
cyclo-(Gly18)
cyclo-(Gly20)
cyclo-((D-Ala-L-Ala)6)
cyclo-((D-Ala-L-Ala)7)
cyclo-((D-Ala-L-Ala)8)
cyclo-((D-Ala-L-Ala)9)
cyclo-((D-Ala-L-Ala)10)
move
mean lowest energy
lowest energy
MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide)
−10.076 −10.514 −10.436 −15.066 −15.770 −15.795 −17.997 −19.552 −19.007 −22.833 −21.548 −23.480 −25.779 −25.788 −27.465 −12.207 −13.243 −13.252 −17.690 −18.380 −18.291 −25.577 −26.638 −26.438 −29.764 −30.590 −30.540 −35.782 −35.578 −35.470
−10.370 −10.514 −10.514 −15.983 −16.540 −16.540 −19.140 −21.067 −21.174 −24.282 −26.437 −24.674 −27.779 −29.236 −29.471 −13.252 −13.252 −13.252 −18.637 −18.637 −18.637 −26.480 −26.814 −26.814 −31.668 −31.668 −31.668 −37.079 −38.493 −38.434
mean CPU time (s) 2.1 5.8 7.3 3.1 8.0 1.0 4.7 1.1 1.5 5.8 1.9 1.9 6.8 2.5 2.5 4.3 1.3 1.6 5.9 1.9 2.5 8.5 2.6 3.4 1.1 3.3 4.3 1.2 4.3 5.6
× × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×
104 103 103 104 103 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 105 104 103 105 104 104
the dimer of cyclo-((D-Ala-L-Ala)3) is a parallel nanotube with six intermolecular hydrogen bonds (Figure 7a). The antiparallel nanotube structure (Figure 7b) lies 9 kcal mol−1 above this. The parallel structure is also the global free energy minimum, with the antiparallel structure less stable by 2.8 kcal mol−1 . In both cases, the peptide groups are tilted away from the optimal alignment for nanotube growth, which is consistent with the dimer being the largest aggregate that has been observed in experiments on cyclic hexapeptides.2,13 The free energy landscape of the cyclo-((D-Ala-L-Ala)3) dimer shows two funnels that correspond to parallel and antiparallel structures
nonpolar solvent. This represents a substantial destabilization of the cis isomers, which are about 4 kcal mol−1 above the trans isomers in acyclic peptides and unstrained cyclic peptides.40 The free energy landscape of cyclo-((D-Ala-L-Ala)4) is also funnelled (Figure 6). The lowest free energy structures in both solvents have some intramolecular hydrogen bonds. Cis isomers are destabilized to a lesser extent than in cyclo-((DAla-L-Ala)3) but are still at least 6.3 kcal mol−1 above the lowest all-trans structures. 3.2. Self-assembly. Cyclic peptides with repeating D- and L-residues show a tendency for self-assembly, and we focus on these in this section. The global potential energy minimum for 1813
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
7a). Perfect antiparallel and parallel structures are 10.8 and 27.1 kcal mol−1 less stable, respectively. The disconnectivity graph of the cyclo-((D-Ala-L-Ala)4) dimer (Figure 8) comprises three main funnels that become
Figure 5. Disconnectivity graphs showing the free energy landscape of cyclo-((D-Ala-L-Ala)3) at 298 K in water (left) and nonpolar solvent (right).
Figure 6. Disconnectivity graphs showing the free energy landscape of cyclo-((D-Ala-L-Ala)4) at 298 K in water (left) and nonpolar solvent (right).
Figure 8. Low-lying minima of the cyclo-((D-Ala-L-Ala)4) dimer. Energies are in kcal mol−1 and are relative to the global potential or free energy minimum. Also shown is a disconnectivity graph of the free energy landscape at 298 K. The 610 minima and 832 transition states accessible by transition states lower than 13 kcal mol−1 from the global minimum are shown. Coloring is according to hydrogen bonding pattern from parallel (blue) to antiparallel (red).
mutually accessible about 12 kcal mol−1 above the global minimum. The funnel containing the global minimum contains several structures with both inter- and intramolecular hydrogen bonds. The next deepest funnel contains the antiparallel nanotube and is competitive with the global minimum, with the lowest structures 0.4 kcal mol−1 above the global free energy minimum. The third funnel contains parallel structures and has minima at least 1.4 kcal mol−1 higher than the global free energy minimum. Note that several distorted structures are lower in energy than the perfect nanotubes in both of these funnels. We have tested the performance of MD and group rotation moves by measuring the mean time to the first encounter of the global potential energy minimum over 50 basin-hopping runs (Table 5). For the cyclo-((D-Ala-L-Ala)3) dimer, both move sets require a similar number of quenches to locate the global minimum. However, for cyclo-((D-Ala-L-Ala)4), MD moves require substantially more quenches to locate the global minimum. For both cyclic peptide dimers, the computational overheads of the MD moves lead to each attempted move taking 2−3 times as much CPU time as the group rotation moves.
Figure 7. Low-lying minima of cyclo-((D-Ala-L-Ala)3) dimers in nonpolar solvent. Energies are in kcal mol−1 and are relative to the global potential or free energy minimum. Also shown is a disconnectivity graph of the free energy landscape at 298 K. The 542 minima and 810 transition states accessible by transition states lower than 13 kcal mol−1 from the global minimum are shown. Coloring is according to hydrogen bonding pattern from parallel (blue) to antiparallel (red).
(Figure 7). The barrier for the parallel to antiparallel conversion is 12.8 kcal mol−1. For cyclo-((D-Ala-L-Ala)4), the global potential energy minimum does not resemble a nanotube and has both intermolecular and intramolecular hydrogen bonds (Figure 1814
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
University’s research community. See http://www.birmingham. ac.uk/bear for more details.
Table 5. Mean Times to the First Encounter of the Global Minimum in 50 Basin-Hopping Runs on Cyclic Peptide Dimers
■
mean first encounter time molecule
move
energy evaluations
cyclo-((D-Ala-L-Ala)3)
MD GR MD GR
1.6 5.8 8.9 1.9
cyclo-((D-Ala-L-Ala)4)
minimizations
CPU time (s)
× × × ×
880 850 4100 2200
2700 1000 27000 6200
106 105 106 106
(1) Ghadiri, M. R.; Granja, J. R.; Milligan, R. A.; McRee, D. E.; Khazanovich, N. Nature 1993, 366, 324−327. (2) Hartgerink, J. D.; Granja, J. R.; Milligan, R. A.; Ghadiri, M. R. J. Am. Chem. Soc. 1996, 118, 43−50. (3) Kim, H. S.; Hartgerink, J. D.; Ghadiri, M. R. J. Am. Chem. Soc. 1998, 120, 4417−4424. (4) Zhu, J.; Cheng, J.; Liao, Z.; Lai, Z.; Liu, B. J. Comput. Aided Mol. Des. 2008, 22, 773−781. (5) Chapman, R.; Danial, M.; Koh, M. L.; Jolliffe, K. A.; Perrier, S. Chem. Soc. Rev. 2012, 41, 6023−6041. (6) Fernandez-Lopez, S.; Kim, H.-S.; Choi, E. C.; Delgado, M.; Granja, J. R.; Khasanov, A.; Kraehenbuehl, K.; Long, G.; Weinberger, D. A.; Wilcoxen, K. M.; Ghadiri, M. R. Nature 2001, 412, 452−455. (7) Dartois, V.; Sanchez-Quesada, J.; Cabezas, E.; Chi, E.; Dubbelde, C.; Dunn, C.; Granja, J.; Gritzen, C.; Weinberger, D.; Ghadiri, M. R.; Parr, T. R. Antimicrob. Agents Chemother. 2005, 49, 3302−3310. (8) Ghadiri, M. R.; Granja, J. R.; Buehler, L. K. Nature 1994, 369, 301−304. (9) Granja, J. R.; Ghadiri, M. R. J. Am. Chem. Soc. 1994, 116, 10785− 10786. (10) Richman, M.; Wilk, S.; Chemerovski, M.; Wärmländer, S. K. T. S.; Wahlström, A.; Gräslund, A.; Rahimipour, S. J. Am. Chem. Soc. 2013, 135, 3474−3484. (11) ten Cate, M. G. J.; Severin, N.; Börner, H. G. Macromolecules 2006, 39, 7831−7838. (12) Bertran, O.; Curcó, D.; Zanuy, D.; Alemán, C. Faraday Discuss. 2013, 166, 59. (13) Sun, X.; Lorenzi, G. P. Helv. Chim. Acta 1994, 77, 1520−1526. (14) Khazanovich, N.; Granja, J. R.; McRee, D. E.; Milligan, R. A.; Ghadiri, M. R. J. Am. Chem. Soc. 1994, 116, 6011−6012. (15) Kobayashi, K.; Granja, J. R.; Ghadiri, M. R. Angew. Chem., Int. Ed. Engl. 1995, 34, 95−98. (16) Khurana, E.; Nielsen, S. O.; Ensing, B.; Klein, M. L. J. Phys. Chem. B 2006, 110, 18965−18972. (17) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. Phys. Chem. Chem. Phys. 2012, 14, 15135. (18) Vijayaraj, R.; Sundar Raman, S.; Mahesh Kumar, R.; Subramanian, V. J. Phys. Chem. B 2010, 114, 16574−16583. (19) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. J. Phys. Chem. B 2012, 116, 9922−9933. (20) Comer, J.; Dehez, F.; Cai, W.; Chipot, C. J. Phys. Chem. C 2013, 117, 26797−26803. (21) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. Phys. Chem. Chem. Phys. 2013, 15, 1260. (22) Khurana, E.; DeVane, R. H.; Kohlmeyer, A.; Klein, M. L. Nano Lett. 2008, 8, 3626−3630. (23) Khalfa, A.; Treptow, W.; Maigret, B.; Tarek, M. Chem. Phys. 2009, 358, 161−170. (24) Khalfa, A.; Tarek, M. J. Phys. Chem. B 2010, 114, 2676−2684. (25) Bonnet, P.; Agrafiotis, D. K.; Zhu, F.; Martin, E. J. Chem. Inf. Model. 2009, 49, 2242−2259. (26) Büttner, F.; Norgren, A. S.; Zhang, S.; Prabpai, S.; Kongsaeree, P.; Arvidsson, P. I. Chem.Eur. J. 2005, 11, 6145−6158. (27) Labute, P. J. Chem. Inf. Model. 2010, 50, 792−800. (28) Rayan, A.; Senderowitz, H.; Goldblum, A. J. Mol. Graphics Modell. 2004, 22, 319−333. (29) Rezai, T.; Bock, J. E.; Zhou, M. V.; Kalyanaraman, C.; Lokey, R. S.; Jacobson, M. P. J. Am. Chem. Soc. 2006, 128, 14073−14080. (30) Goldtzvik, Y.; Goldstein, M.; Benny Gerber, R. Chem. Phys. 2013, 415, 168−172. (31) Dieterich, J. M.; Hartke, B. Mol. Phys. 2009, 108, 279−291. (32) Dieterich, J. M.; Gerstel, U.; Schröder, J.-M.; Hartke, B. J. Mol. Model. 2011, 17, 3195−3207. (33) Strodel, B.; Wales, D. J. J. Chem. Theory Comput. 2008, 4, 657− 672.
4. CONCLUSIONS We have tested the performance of several move classes in basin-hopping optimization of cyclic peptides and their dimers. Simple Cartesian displacement moves perform poorly because of their tendency to invert chiral centers. MD moves are better, but the computational cost of a single basin-hopping step is significantly higher. The best performance, in terms of CPU time and the effectiveness of each basin-hopping step, is obtained with group rotation moves. When searching the conformations of cyclic peptide dimers, a hybrid minimization strategy, where quenches are performed treating each molecule as a rigid body before an all-atom quench, is very effective. The combination of group rotation and rigid body moves is more efficient than MD moves, both in terms of the number of basin-hopping steps required to locate the global minimum and the computational cost of each step. Analysis of the energy landscape of the cyclo-((D-Ala-L-Ala)3) dimer shows that the parallel structure is the global minimum, but also that this structure is not well-aligned to seed the growth of a nanotube. For the cyclo-((D-Ala-L-Ala)4) dimer there are several competing low-energy structures, and some of these could proceed to form nanotubes. In the future, we will consider the formation of larger cyclic peptide nanotubes. From some preliminary calculations, we find that basin hopping searches using all atom models are practical for systems up to at least four cyclic peptides. For nanotubes that are too large to study with all-atom models, we will use a course-grained potential, such as the Martini model.61 We have used an implicit solvation model in this study. Cyclic peptide nanotubes are large enough to contain small molecules, and it is possible that explicit solvent models will lead to different hydrogen bonding patterns. Basin-hopping calculations including a small number of explicit solvent molecules could be a computationally feasible way of capturing this effect.
■
REFERENCES
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors thank Prof. David Wales for helpful discussions and Dr. Chris Whittleston for implementing some of the Monte Carlo move classes in GMIN. We acknowledge the Engineering and Physical Sciences Research Council, U.K. (EPSRC) for funding under Programme Grant EP/I001352/1. The computations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service, which provides a High Performance Computing service to the 1815
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
(34) Strodel, B.; Lee, J. W. L.; Whittleston, C. S.; Wales, D. J. J. Am. Chem. Soc. 2010, 132, 13300−13312. (35) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111− 5116. (36) Wales, D. J. Mol. Phys. 2002, 100, 3285−3305. (37) Wales, D. J. In Energy Landscapes; Wales, D. J., Ed.; Cambridge University Press: Cambridge, U.K., 2003; pp 397−409. (38) Wales, D. J. Mol. Phys. 2004, 102, 891−908. (39) Oakley, M. T.; Johnston, R. L. J. Chem. Theory Comput. 2013, 9, 650−657. (40) Oakley, M. T.; Oheix, E.; Peacock, A. F. A.; Johnston, R. L. J. Phys. Chem. B 2013, 117, 8122−8134. (41) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999−2012. (42) Lee, M. C.; Duan, Y. Proteins 2004, 55, 620−634. (43) Onufriev, A.; Bashford, D.; Case, D. A. Proteins 2004, 55, 383− 394. (44) Wales, D. J. GMIN: A program for finding global minima and calculating thermodynamic properties from basin-sampling. http:// www-wales.ch.cam.ac.uk/GMIN/(accessed Oct. 11, 2013). (45) Somani, S.; Wales, D. J. J. Chem. Phys. 2013, 139, 121909. (46) Mochizuki, K.; Whittleston, C. S.; Somani, S.; Kusumaatmaja, H.; Wales, D. J. Phys. Chem. Chem. Phys. 2014, 16, 2842−2853. (47) Kusumaatmaja, H.; Whittleston, C. S.; Wales, D. J. J. Chem. Theory Comput. 2012, 8, 5159−5165. (48) Wales, D. J. PATHSAMPLE: A program for refining and analysing kinetic transition networks. http://www-wales.ch.cam.ac.uk/ PATHSAMPLE/(accessed Oct. 11, 2013). (49) Wales, D. J. OPTIM: A program for characterizing stationary points and reaction pathways. http://www-wales.ch.cam.ac.uk/ OPTIM/ (accessed Oct. 11, 2013). (50) Trygubenko, S. A.; Wales, D. J. J. Chem. Phys. 2004, 120, 2082− 2094. (51) Henkelman, G.; Jónsson, H. J. Chem. Phys. 1999, 111, 7010− 7022. (52) Munro, L. J.; Wales, D. J. Phys. Rev. B 1999, 59, 3969−3980. (53) Kumeda, Y.; Munro, L. J.; Wales, D. J. Chem. Phys. Lett. 2001, 341, 185−194. (54) Carr, J. M.; Trygubenko, S. A.; Wales, D. J. J. Chem. Phys. 2005, 122, 234903. (55) Strodel, B.; Whittleston, C. S.; Wales, D. J. J. Am. Chem. Soc. 2007, 129, 16005−16014. (56) Becker, O. M.; Karplus, M. J. Chem. Phys. 1997, 106, 1495− 1517. (57) Wales, D. J.; Miller, M. A.; Walsh, T. R. Nature 1998, 394, 758− 760. (58) Krivov, S. V.; Karplus, M. J. Chem. Phys. 2002, 117, 10894− 10903. (59) Evans, D. A.; Wales, D. J. J. Chem. Phys. 2003, 118, 3891−3897. (60) Kabsch, W.; Sander, C. Biopolymers 1983, 22, 2577−2637. (61) Marrink, S. J.; Tieleman, D. P. Chem. Soc. Rev. 2013, 42, 6801− 6822.
1816
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Terms of Use CC-BY
Energy Landscapes and Global Optimization of Self-Assembling Cyclic Peptides Mark T. Oakley and Roy L. Johnston* School of Chemistry, University of Birmingham, Edgbaston, Birmingham, B15 2TT, United Kingdom ABSTRACT: Self-assembled cyclic peptide nanotubes have attracted much attention because of their antimicrobial properties. Here, we present calculations on the formation of cyclic peptide dimers using basin-hopping and discrete path sampling. We present an analysis of the basin-hopping move sets that most efficiently explore the conformations of cyclic peptides. Group rotation moves, in which sections of the ring are rotated as a rigid body, are the most effective for cyclic peptides containing up to 20 residues. For cyclic peptide dimers, we find that a combination of group rotation intramolecular moves and rigid body intermolecular moves performs well. Discrete path sampling calculations on the cyclic peptide dimers show significant differences in the dimerization of hexa- and octapeptides. nanotubes have been shown to form across22,23 or along24 the membrane, depending on the peptide sequence. Here, we take a different approach and consider only the stationary points on the energy landscape. The solvent is treated implicitly, which substantially reduces the computational resources required. Locating minima on the energy landscape then becomes a global optimization problem that is split into two parts: searching the conformations of individual cyclic peptides and assembly of cyclic peptide aggregates. Several methods are available to search the conformations of small flexible molecules. Those that have been used on cyclic peptides include distance geometry methods,25 Monte Carlo,26 LowMode molecular dynamics,27 iterative stochastic search,28 dihedral angle sampling29 and diffusion equation evolutionary programming simulated annealing.30 Optimization of clusters of flexible molecules is a more difficult problem. An evolutionary algorithm has been used to optimize the dimer of kanamycin A.31,32 Basin-hopping has been used to optimize the small aggregates of the amyloid-β peptide.33,34 Here, we use the basin-hopping algorithm35 to locate the minima on the cyclic peptide potential energy surfaces. To understand the conformations of molecules fully we need to know both the minima and the transition states linking them. Discrete path sampling36−38 is an efficient method for exploring energy landscapes and we have previously used it to study the conformations of peptides.39,40 Here, we present energy landscapes for the peptide dimers from discrete path sampling calculations.
1. INTRODUCTION Self-assembly is the process by which small molecules aggregate to form large, ordered structures. Some head-to-tail cyclic peptides are known to self-assemble to form nanotubes.1−5 Self-assembling cyclic peptide nanotubes have attracted attention for their possible biomedical uses, including activity against bacteria.6,7 Cyclic peptide nanotubes have been observed to form across lipid bilayers.3,8,9 The hydrogenbonding arrangement of cyclic peptide nanotubes is similar to that seen in amyloid fibrils, and it has been shown that cyclic peptides can disrupt amyloid formation.10 Cyclic peptide− polymer conjugates have potential uses as self-assembling materials.11,12 Self-assembling cyclic peptides have sequences comprising alternating D- and L-peptide residues because this arrangement places the side-chains equatorial to the peptide rings with the peptide groups aligned axially to hydrogen bond to other cyclic peptides. Cyclic octapeptides are the most widely studied because they have the highest tendency to form nanotubes.1−4 However, dimerization of cyclic hexapeptides has been observed,13 and nanotubes of larger cyclic peptides are known.9,14 Nanotubes with antiparallel hydrogen bonding arrangements are more stable than those with parallel hydrogen bonds.15 Several molecular dynamics (MD) investigations have been undertaken to model the properties of cyclic peptide nanotubes. The formation of dimers has been modeled with MD.16 The composition of cyclic peptides has a substantial effect on the stability and flexibility of their nanotubes.17 Longer nanotubes are more stable than shorter nanotubes.18,19 MD has been used to study the motion of water20 and other small molecules21 through cyclic peptide nanotubes. Cyclic peptide−polymer conjugates have also been studied with MD.12 All-atom MD is very computationally demanding, and modeling the formation of nanotubes is prohibitively expensive. For this reason, coarse-grained models have been used to study the aggregation of cyclic peptides in lipid bilayers, and © 2014 American Chemical Society
2. METHODS The energies of all structures were evaluated using the AMBER FF03 force field.41,42 The effect of solvent was included with the Generalized Born implicit solvation model.43 Calculations Received: January 3, 2014 Published: February 25, 2014 1810
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
were performed in water (εr = 80) or in a low-dielectric solvent (εr = 2) to model a membrane environment. Global optimization of single cyclic peptides and dimers was performed using the basin-hopping algorithm35 as implemented in GMIN.44 In this algorithm, structural deformations are followed by a gradient-driven quench to a local minimum before a Monte Carlo acceptance test. Thus, all attempted moves are between local minima on the potential energy surface. For efficient searching, the moves must be sufficiently large to escape the basin of attraction of the initial local minimum. However, moves that enter basins that are not physically accessible should be avoided. Inversion of the chirality of the Cα atoms is not possible in experiment and has a barrier of about 70 kcal mol−1 in the Amber force field.45 Therefore, any moves that lead to structures where the chirality of one or more Cα atoms are inverted are rejected. We consider three types of intramolecular move: atom displacements, MD moves, and group rotations. A series of basin-hopping optimizations were performed to benchmark each move class. For cyclic peptides comprising 10 or fewer peptide residues, the global minimum structures are found by all searches after fairly short basin-hopping runs. In these cases, the performance was measured in terms of the mean time to the first encounter of the global minimum from 100 independent searches. For larger cyclic peptides, the conformational spaces are too large for us to be confident that the lowest minima we find are the global minima. In these cases, basinhopping runs were run for 10 000 Monte Carlo steps and the performance was measured in terms of the lowest energies found at the end of each run. The initial structures for all of these searches were taken from a molecular dynamics simulation at 1600 K with snapshots taken every 10 ps. The parameters used for all move classes were optimized by performing sets of 10 searches with different values of the parameters for peptides containing up to 10 residues. These initial 10 searches are included in the set of 100 searches used to calculate the mean first encounter time. Unless noted, the same search parameters are used for all of the cyclic peptides studied here. The atom displacement moves (STEP keyword) involve displacement of all the Cartesian coordinates of all atoms by a up to 1.5 Å. In the MD moves (AMBERMDSTEPS keyword), short MD runs of 1000 steps are performed to escape the basins of attraction of local minima. These runs are performed at high temperature of 1600 K to allow the move to rapidly surmount large barriers, such as those between cis and trans isomers. In the group rotation moves (GROUPROTATION keyword), two atoms in the peptide backbone are selected and all atoms between them are rotated around the axis defined by those two atoms (Figure 1).46 Four types of group rotation are used (Table 1): cis−trans, peptide, dipeptide, and tripeptide moves. Cis−trans and peptide moves are attempted for all cyclic peptides. Dipeptide moves are attempted for cyclic peptides with six or more residues. The performance of tripeptide moves is evaluated for cyclic peptides comprising 12 or more residues. The magnitudes of all group rotation moves are in the range ±180°. The probability of selecting a group for group rotation is defined so that each basin-hopping step comprises on average 1/3 cis−trans moves, three peptide moves and, where appropriate, one of each of the larger moves. A combined set of intramolecular and intermolecular moves is needed to search the structures of dimers, and we consider two classes of move. MD moves use a short, high-temperature
Figure 1. Atoms involved in a dipeptide group rotation move. The axis is defined by two Cα atoms (yellow). All of the intermediate atoms (purple) are rotated around this axis. Some hydrogen atoms are omitted for clarity.
Table 1. Atoms Involved in the Each Type of Group Rotation Move move cis−trans peptide dipeptide tripeptide
ends Cα,i Cα,i Cα,i Cα,i
Ni+1 Cα,i+1 Cα,i+2 Cα,i+3
MD run to perform both the intra- and intermolecular moves. We also use a combination of group rotation intramolecular moves and rigid-body47 intermolecular moves. For both types of move, a hybrid minimization scheme is used (Figure 2). The second step is particularly important when using MD moves because the two cyclic peptides can drift apart during the MD run.
Figure 2. Basin-hopping scheme with hybrid minimization.
We evaluate the move sets on two types of cyclic peptides. Cyclic oligoglycines do not have any side chains and allow us to study the effect of the moves on the peptide backbones. We have not studied the dimerization of these peptides because they are not known to form nanotubes. Cyclic peptides with alternating D -Ala and L-Ala residues were selected as representative cyclic D-, L-peptides. Other peptide sequences form more stable nanotubes,17 but alanine was selected to reduce the effort of modeling the peptide side chains. The energy landscapes of some of the cyclic peptides and dimers were explored with discrete path sampling,36−38 as implemented in PATHSAMPLE.48 The databases of stationary points were seeded with minima from the basin-hopping calculations. Pathways between minima were found using OPTIM,49 with candidate transition states found by doubly nudged elastic band calculations50 and refined with hybrid eigenvector following.51−53 Pairs of minima were initially selected for connection by the missing connections algorithm54 and later using the UNTRAP method to remove artificial frustration.55 The resulting databases of stationary points are visualized as disconnectivity graphs.56,57 Stationary points that are permutation-inversion isomers are grouped together in 1811
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
second-lowest structure several times before locating the global minimum. The displacement moves require a similar amount of CPU time to MD moves when n ≤ 8. We have also tested the performance of basin-hopping on a series of cyclic peptides with alternating D-Ala and L-Ala residues. We have previously shown that the global minimum of cyclo-((D-Ala-L-Ala)2) includes one peptide in a cis conformation (Figure 4a).39 The global minima for all of the larger cyclo-((D-Ala-L-Ala)n) have all of their peptide bonds in trans conformations (Figure 4).
these disconnectivity graphs. All disconnectivity graphs show the free energy surface at 298 K. The free energies of all stationary points are calculated using the harmonic superposition approximation.58,59 We use the pattern of hydrogen bonds as a metric to measure the parallel or antiparallel alignment of these dimers. Hydrogen bonds in the dimers were counted using the definition from the Dictionary of Secondary Structure Prediction.60 Intermolecular hydrogen bonds were then assigned as parallel or antiparallel based on the angles between the two vectors joining the Cα atoms on either side of the peptide groups involved. The alignment is defined as npar−nantipar.
3. RESULTS 3.1. Intramolecular Moves. We have tested the performance of basin-hopping moves on a series of cyclic oligoglycines. For consistency with previous studies,39,40 all basin-hopping searches on single cyclic peptides were performed using water as a solvent. We have previously shown that the global minimum for cyclo-(Gly4) has all four peptide groups aligned axially, with no intramolecular hydrogen bonds (Figure 3a).39 The global minima for all of the larger cyclic oligoglycines exhibit several intramolecular hydrogen bonds (Figure 3).
Figure 4. Global minima of cyclo-((D-Ala-L-Ala)n) in water.
The performance of the displacement moves is poor for these molecules. Cartesian displacements that are sufficiently large to step between basins on the potential energy surface also lead to the inversion of at least one of the chiral centers in about 45% of the attempted moves. For comparison, MD moves never lead to inversion of chiral centers and group rotation moves invert in 5−10% of attempts. In most cases, the group rotation moves give the most efficient optimization. The MD moves locate the global minimum of cyclo-((D-Ala-LAla)4) (Figure 4c) very rapidly (Table 3). However, this seems to be a particularly easy global minimum to locate because the mean first encounter times with all move types are faster than the corresponding times for cyclo-((D-Ala-L-Ala)3). Basin-hopping optimization has also been performed on several larger cyclic peptides comprising up to 20 residues (table 4). For these molecules, we only consider the MD and group rotation moves because of the poor performance of the displacement moves in the previous calculations. The performance was assessed in terms of the mean energy of the lowest structure found in each basin-hopping run and the lowest energy found across all basin-hopping runs. For all of the cyclic peptides, the lowest minimum found using group rotation moves up to and including the dipeptide move is more stable than the lowest found with molecular dynamics moves. By both measures, group rotation with dipeptide moves out-performs MD moves for all cyclic oligoglycines. The inclusion of tripeptide moves makes only a small difference to the lowest energies for the smaller rings but gives a substantial improvement for cyclo-(Gly20). The CPU time taken for each Monte Carlo step is also substantially longer for the MD moves. The energy landscape of cyclo-((D-Ala-L-Ala)3) is funnelled, with small downhill barriers separating all local minima from the global minimum (Figure 5). The conformation with the lowest free energy in water is similar to that seen in cyclic peptide nanotubes. In nonpolar solvent, the lowest free energy structure exhibits some intramolecular hydrogen bonding. In cyclo-((D-Ala-L-Ala)3), conformers containing cis peptides are particularly unstable and are at least 8.8 kcal mol−1 above the global minimum in water and 10.3 kcal mol−1 higher in
Figure 3. Global minima of cyclo-(Glyn) in water.
For all four of the cyclic oligoglycines, the group rotation moves give the most efficient optimization in terms of the CPU time and the number of energy evaluations (Table 2). For n ≤ 8, the MD moves require a similar number of quenches to the group rotation moves but need more CPU time because of the overhead of performing the short MD run. When n = 10, molecular dynamics moves perform poorly and visit the Table 2. Mean Times to the First Encounter of the Global Minimum for 100 Basin-Hopping Runs on Cyclic Oligoglycines mean first encounter time molecule
move
cyclo-(Gly4)
displacement MD GR displacement MD GR displacement MD GR displacement MD GR
cyclo-(Gly6)
cyclo-(Gly8)
cyclo-(Gly10)
energy evaluations
minimizations
CPU time (s)
× × × × × × × × × × × ×
99 32 18 90 21 51 680 190 260 1000 1700 420
2.6 4.7 0.45 7.9 6.8 3.8 140 120 40 370 1900 120
1.8 3.7 3.2 2.6 2.7 1.3 3.0 2.9 8.8 5.2 3.0 1.8
104 104 103 104 104 104 105 105 104 105 106 105
1812
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
Table 3. Mean Times to the First Encounter of the Global Minimum in 100 Basin-Hopping Runs on Cyclic Oligo-D-,L-alanines mean first encounter time molecule cyclo-((D-Ala-L-Ala)2)
cyclo-((D-Ala-L-Ala)3)
cyclo-((D-Ala-L-Ala)4)
cyclo-((D-Ala-L-Ala)5)
move
energy evaluations
minimizations
CPU time (s)
displacement MD GR displacement MD GR displacement MD GR displacement MD GR
× × × × × × × × × × × ×
180 13 44 5200 1100 300 2600 39 130 12000 3900 530
11 3.7 2.7 1000 670 67 1100 48 63 8100 7700 430
4.5 1.6 1.1 1.9 1.5 1.2 1.2 6.2 6.9 6.4 6.6 3.3
4
10 104 104 106 106 105 106 104 104 106 106 105
Table 4. Lowest Energies Found in 20 × 10 000 Step Basin-Hopping Runs on Cyclic Peptides with 12 or More Residues molecule cyclo-(Gly12)
cyclo-(Gly14)
cyclo-(Gly16)
cyclo-(Gly18)
cyclo-(Gly20)
cyclo-((D-Ala-L-Ala)6)
cyclo-((D-Ala-L-Ala)7)
cyclo-((D-Ala-L-Ala)8)
cyclo-((D-Ala-L-Ala)9)
cyclo-((D-Ala-L-Ala)10)
move
mean lowest energy
lowest energy
MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide) MD GR (dipeptide) GR (tripeptide)
−10.076 −10.514 −10.436 −15.066 −15.770 −15.795 −17.997 −19.552 −19.007 −22.833 −21.548 −23.480 −25.779 −25.788 −27.465 −12.207 −13.243 −13.252 −17.690 −18.380 −18.291 −25.577 −26.638 −26.438 −29.764 −30.590 −30.540 −35.782 −35.578 −35.470
−10.370 −10.514 −10.514 −15.983 −16.540 −16.540 −19.140 −21.067 −21.174 −24.282 −26.437 −24.674 −27.779 −29.236 −29.471 −13.252 −13.252 −13.252 −18.637 −18.637 −18.637 −26.480 −26.814 −26.814 −31.668 −31.668 −31.668 −37.079 −38.493 −38.434
mean CPU time (s) 2.1 5.8 7.3 3.1 8.0 1.0 4.7 1.1 1.5 5.8 1.9 1.9 6.8 2.5 2.5 4.3 1.3 1.6 5.9 1.9 2.5 8.5 2.6 3.4 1.1 3.3 4.3 1.2 4.3 5.6
× × × × × × × × × × × × × × × × × × × × × × × × × × × × × ×
104 103 103 104 103 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 104 105 104 103 105 104 104
the dimer of cyclo-((D-Ala-L-Ala)3) is a parallel nanotube with six intermolecular hydrogen bonds (Figure 7a). The antiparallel nanotube structure (Figure 7b) lies 9 kcal mol−1 above this. The parallel structure is also the global free energy minimum, with the antiparallel structure less stable by 2.8 kcal mol−1 . In both cases, the peptide groups are tilted away from the optimal alignment for nanotube growth, which is consistent with the dimer being the largest aggregate that has been observed in experiments on cyclic hexapeptides.2,13 The free energy landscape of the cyclo-((D-Ala-L-Ala)3) dimer shows two funnels that correspond to parallel and antiparallel structures
nonpolar solvent. This represents a substantial destabilization of the cis isomers, which are about 4 kcal mol−1 above the trans isomers in acyclic peptides and unstrained cyclic peptides.40 The free energy landscape of cyclo-((D-Ala-L-Ala)4) is also funnelled (Figure 6). The lowest free energy structures in both solvents have some intramolecular hydrogen bonds. Cis isomers are destabilized to a lesser extent than in cyclo-((DAla-L-Ala)3) but are still at least 6.3 kcal mol−1 above the lowest all-trans structures. 3.2. Self-assembly. Cyclic peptides with repeating D- and L-residues show a tendency for self-assembly, and we focus on these in this section. The global potential energy minimum for 1813
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
7a). Perfect antiparallel and parallel structures are 10.8 and 27.1 kcal mol−1 less stable, respectively. The disconnectivity graph of the cyclo-((D-Ala-L-Ala)4) dimer (Figure 8) comprises three main funnels that become
Figure 5. Disconnectivity graphs showing the free energy landscape of cyclo-((D-Ala-L-Ala)3) at 298 K in water (left) and nonpolar solvent (right).
Figure 6. Disconnectivity graphs showing the free energy landscape of cyclo-((D-Ala-L-Ala)4) at 298 K in water (left) and nonpolar solvent (right).
Figure 8. Low-lying minima of the cyclo-((D-Ala-L-Ala)4) dimer. Energies are in kcal mol−1 and are relative to the global potential or free energy minimum. Also shown is a disconnectivity graph of the free energy landscape at 298 K. The 610 minima and 832 transition states accessible by transition states lower than 13 kcal mol−1 from the global minimum are shown. Coloring is according to hydrogen bonding pattern from parallel (blue) to antiparallel (red).
mutually accessible about 12 kcal mol−1 above the global minimum. The funnel containing the global minimum contains several structures with both inter- and intramolecular hydrogen bonds. The next deepest funnel contains the antiparallel nanotube and is competitive with the global minimum, with the lowest structures 0.4 kcal mol−1 above the global free energy minimum. The third funnel contains parallel structures and has minima at least 1.4 kcal mol−1 higher than the global free energy minimum. Note that several distorted structures are lower in energy than the perfect nanotubes in both of these funnels. We have tested the performance of MD and group rotation moves by measuring the mean time to the first encounter of the global potential energy minimum over 50 basin-hopping runs (Table 5). For the cyclo-((D-Ala-L-Ala)3) dimer, both move sets require a similar number of quenches to locate the global minimum. However, for cyclo-((D-Ala-L-Ala)4), MD moves require substantially more quenches to locate the global minimum. For both cyclic peptide dimers, the computational overheads of the MD moves lead to each attempted move taking 2−3 times as much CPU time as the group rotation moves.
Figure 7. Low-lying minima of cyclo-((D-Ala-L-Ala)3) dimers in nonpolar solvent. Energies are in kcal mol−1 and are relative to the global potential or free energy minimum. Also shown is a disconnectivity graph of the free energy landscape at 298 K. The 542 minima and 810 transition states accessible by transition states lower than 13 kcal mol−1 from the global minimum are shown. Coloring is according to hydrogen bonding pattern from parallel (blue) to antiparallel (red).
(Figure 7). The barrier for the parallel to antiparallel conversion is 12.8 kcal mol−1. For cyclo-((D-Ala-L-Ala)4), the global potential energy minimum does not resemble a nanotube and has both intermolecular and intramolecular hydrogen bonds (Figure 1814
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
University’s research community. See http://www.birmingham. ac.uk/bear for more details.
Table 5. Mean Times to the First Encounter of the Global Minimum in 50 Basin-Hopping Runs on Cyclic Peptide Dimers
■
mean first encounter time molecule
move
energy evaluations
cyclo-((D-Ala-L-Ala)3)
MD GR MD GR
1.6 5.8 8.9 1.9
cyclo-((D-Ala-L-Ala)4)
minimizations
CPU time (s)
× × × ×
880 850 4100 2200
2700 1000 27000 6200
106 105 106 106
(1) Ghadiri, M. R.; Granja, J. R.; Milligan, R. A.; McRee, D. E.; Khazanovich, N. Nature 1993, 366, 324−327. (2) Hartgerink, J. D.; Granja, J. R.; Milligan, R. A.; Ghadiri, M. R. J. Am. Chem. Soc. 1996, 118, 43−50. (3) Kim, H. S.; Hartgerink, J. D.; Ghadiri, M. R. J. Am. Chem. Soc. 1998, 120, 4417−4424. (4) Zhu, J.; Cheng, J.; Liao, Z.; Lai, Z.; Liu, B. J. Comput. Aided Mol. Des. 2008, 22, 773−781. (5) Chapman, R.; Danial, M.; Koh, M. L.; Jolliffe, K. A.; Perrier, S. Chem. Soc. Rev. 2012, 41, 6023−6041. (6) Fernandez-Lopez, S.; Kim, H.-S.; Choi, E. C.; Delgado, M.; Granja, J. R.; Khasanov, A.; Kraehenbuehl, K.; Long, G.; Weinberger, D. A.; Wilcoxen, K. M.; Ghadiri, M. R. Nature 2001, 412, 452−455. (7) Dartois, V.; Sanchez-Quesada, J.; Cabezas, E.; Chi, E.; Dubbelde, C.; Dunn, C.; Granja, J.; Gritzen, C.; Weinberger, D.; Ghadiri, M. R.; Parr, T. R. Antimicrob. Agents Chemother. 2005, 49, 3302−3310. (8) Ghadiri, M. R.; Granja, J. R.; Buehler, L. K. Nature 1994, 369, 301−304. (9) Granja, J. R.; Ghadiri, M. R. J. Am. Chem. Soc. 1994, 116, 10785− 10786. (10) Richman, M.; Wilk, S.; Chemerovski, M.; Wärmländer, S. K. T. S.; Wahlström, A.; Gräslund, A.; Rahimipour, S. J. Am. Chem. Soc. 2013, 135, 3474−3484. (11) ten Cate, M. G. J.; Severin, N.; Börner, H. G. Macromolecules 2006, 39, 7831−7838. (12) Bertran, O.; Curcó, D.; Zanuy, D.; Alemán, C. Faraday Discuss. 2013, 166, 59. (13) Sun, X.; Lorenzi, G. P. Helv. Chim. Acta 1994, 77, 1520−1526. (14) Khazanovich, N.; Granja, J. R.; McRee, D. E.; Milligan, R. A.; Ghadiri, M. R. J. Am. Chem. Soc. 1994, 116, 6011−6012. (15) Kobayashi, K.; Granja, J. R.; Ghadiri, M. R. Angew. Chem., Int. Ed. Engl. 1995, 34, 95−98. (16) Khurana, E.; Nielsen, S. O.; Ensing, B.; Klein, M. L. J. Phys. Chem. B 2006, 110, 18965−18972. (17) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. Phys. Chem. Chem. Phys. 2012, 14, 15135. (18) Vijayaraj, R.; Sundar Raman, S.; Mahesh Kumar, R.; Subramanian, V. J. Phys. Chem. B 2010, 114, 16574−16583. (19) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. J. Phys. Chem. B 2012, 116, 9922−9933. (20) Comer, J.; Dehez, F.; Cai, W.; Chipot, C. J. Phys. Chem. C 2013, 117, 26797−26803. (21) Vijayaraj, R.; Van Damme, S.; Bultinck, P.; Subramanian, V. Phys. Chem. Chem. Phys. 2013, 15, 1260. (22) Khurana, E.; DeVane, R. H.; Kohlmeyer, A.; Klein, M. L. Nano Lett. 2008, 8, 3626−3630. (23) Khalfa, A.; Treptow, W.; Maigret, B.; Tarek, M. Chem. Phys. 2009, 358, 161−170. (24) Khalfa, A.; Tarek, M. J. Phys. Chem. B 2010, 114, 2676−2684. (25) Bonnet, P.; Agrafiotis, D. K.; Zhu, F.; Martin, E. J. Chem. Inf. Model. 2009, 49, 2242−2259. (26) Büttner, F.; Norgren, A. S.; Zhang, S.; Prabpai, S.; Kongsaeree, P.; Arvidsson, P. I. Chem.Eur. J. 2005, 11, 6145−6158. (27) Labute, P. J. Chem. Inf. Model. 2010, 50, 792−800. (28) Rayan, A.; Senderowitz, H.; Goldblum, A. J. Mol. Graphics Modell. 2004, 22, 319−333. (29) Rezai, T.; Bock, J. E.; Zhou, M. V.; Kalyanaraman, C.; Lokey, R. S.; Jacobson, M. P. J. Am. Chem. Soc. 2006, 128, 14073−14080. (30) Goldtzvik, Y.; Goldstein, M.; Benny Gerber, R. Chem. Phys. 2013, 415, 168−172. (31) Dieterich, J. M.; Hartke, B. Mol. Phys. 2009, 108, 279−291. (32) Dieterich, J. M.; Gerstel, U.; Schröder, J.-M.; Hartke, B. J. Mol. Model. 2011, 17, 3195−3207. (33) Strodel, B.; Wales, D. J. J. Chem. Theory Comput. 2008, 4, 657− 672.
4. CONCLUSIONS We have tested the performance of several move classes in basin-hopping optimization of cyclic peptides and their dimers. Simple Cartesian displacement moves perform poorly because of their tendency to invert chiral centers. MD moves are better, but the computational cost of a single basin-hopping step is significantly higher. The best performance, in terms of CPU time and the effectiveness of each basin-hopping step, is obtained with group rotation moves. When searching the conformations of cyclic peptide dimers, a hybrid minimization strategy, where quenches are performed treating each molecule as a rigid body before an all-atom quench, is very effective. The combination of group rotation and rigid body moves is more efficient than MD moves, both in terms of the number of basin-hopping steps required to locate the global minimum and the computational cost of each step. Analysis of the energy landscape of the cyclo-((D-Ala-L-Ala)3) dimer shows that the parallel structure is the global minimum, but also that this structure is not well-aligned to seed the growth of a nanotube. For the cyclo-((D-Ala-L-Ala)4) dimer there are several competing low-energy structures, and some of these could proceed to form nanotubes. In the future, we will consider the formation of larger cyclic peptide nanotubes. From some preliminary calculations, we find that basin hopping searches using all atom models are practical for systems up to at least four cyclic peptides. For nanotubes that are too large to study with all-atom models, we will use a course-grained potential, such as the Martini model.61 We have used an implicit solvation model in this study. Cyclic peptide nanotubes are large enough to contain small molecules, and it is possible that explicit solvent models will lead to different hydrogen bonding patterns. Basin-hopping calculations including a small number of explicit solvent molecules could be a computationally feasible way of capturing this effect.
■
REFERENCES
AUTHOR INFORMATION
Corresponding Author
*E-mail: [email protected]. Notes
The authors declare no competing financial interest.
■
ACKNOWLEDGMENTS The authors thank Prof. David Wales for helpful discussions and Dr. Chris Whittleston for implementing some of the Monte Carlo move classes in GMIN. We acknowledge the Engineering and Physical Sciences Research Council, U.K. (EPSRC) for funding under Programme Grant EP/I001352/1. The computations described in this paper were performed using the University of Birmingham’s BlueBEAR HPC service, which provides a High Performance Computing service to the 1815
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
Journal of Chemical Theory and Computation
Article
(34) Strodel, B.; Lee, J. W. L.; Whittleston, C. S.; Wales, D. J. J. Am. Chem. Soc. 2010, 132, 13300−13312. (35) Wales, D. J.; Doye, J. P. K. J. Phys. Chem. A 1997, 101, 5111− 5116. (36) Wales, D. J. Mol. Phys. 2002, 100, 3285−3305. (37) Wales, D. J. In Energy Landscapes; Wales, D. J., Ed.; Cambridge University Press: Cambridge, U.K., 2003; pp 397−409. (38) Wales, D. J. Mol. Phys. 2004, 102, 891−908. (39) Oakley, M. T.; Johnston, R. L. J. Chem. Theory Comput. 2013, 9, 650−657. (40) Oakley, M. T.; Oheix, E.; Peacock, A. F. A.; Johnston, R. L. J. Phys. Chem. B 2013, 117, 8122−8134. (41) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; Caldwell, J.; Wang, J.; Kollman, P. J. Comput. Chem. 2003, 24, 1999−2012. (42) Lee, M. C.; Duan, Y. Proteins 2004, 55, 620−634. (43) Onufriev, A.; Bashford, D.; Case, D. A. Proteins 2004, 55, 383− 394. (44) Wales, D. J. GMIN: A program for finding global minima and calculating thermodynamic properties from basin-sampling. http:// www-wales.ch.cam.ac.uk/GMIN/(accessed Oct. 11, 2013). (45) Somani, S.; Wales, D. J. J. Chem. Phys. 2013, 139, 121909. (46) Mochizuki, K.; Whittleston, C. S.; Somani, S.; Kusumaatmaja, H.; Wales, D. J. Phys. Chem. Chem. Phys. 2014, 16, 2842−2853. (47) Kusumaatmaja, H.; Whittleston, C. S.; Wales, D. J. J. Chem. Theory Comput. 2012, 8, 5159−5165. (48) Wales, D. J. PATHSAMPLE: A program for refining and analysing kinetic transition networks. http://www-wales.ch.cam.ac.uk/ PATHSAMPLE/(accessed Oct. 11, 2013). (49) Wales, D. J. OPTIM: A program for characterizing stationary points and reaction pathways. http://www-wales.ch.cam.ac.uk/ OPTIM/ (accessed Oct. 11, 2013). (50) Trygubenko, S. A.; Wales, D. J. J. Chem. Phys. 2004, 120, 2082− 2094. (51) Henkelman, G.; Jónsson, H. J. Chem. Phys. 1999, 111, 7010− 7022. (52) Munro, L. J.; Wales, D. J. Phys. Rev. B 1999, 59, 3969−3980. (53) Kumeda, Y.; Munro, L. J.; Wales, D. J. Chem. Phys. Lett. 2001, 341, 185−194. (54) Carr, J. M.; Trygubenko, S. A.; Wales, D. J. J. Chem. Phys. 2005, 122, 234903. (55) Strodel, B.; Whittleston, C. S.; Wales, D. J. J. Am. Chem. Soc. 2007, 129, 16005−16014. (56) Becker, O. M.; Karplus, M. J. Chem. Phys. 1997, 106, 1495− 1517. (57) Wales, D. J.; Miller, M. A.; Walsh, T. R. Nature 1998, 394, 758− 760. (58) Krivov, S. V.; Karplus, M. J. Chem. Phys. 2002, 117, 10894− 10903. (59) Evans, D. A.; Wales, D. J. J. Chem. Phys. 2003, 118, 3891−3897. (60) Kabsch, W.; Sander, C. Biopolymers 1983, 22, 2577−2637. (61) Marrink, S. J.; Tieleman, D. P. Chem. Soc. Rev. 2013, 42, 6801− 6822.
1816
dx.doi.org/10.1021/ct500004k | J. Chem. Theory Comput. 2014, 10, 1810−1816
