Virtual Screening for Dipeptide Aggregation: Toward Predictive Tools for Peptide Self-Assembly.

Several short peptide sequences are known to self-assemble into supramolecular nanostructures with interesting properties. In this study, coarse-grained molecular dynamics is employed to rapidly screen all 400 dipeptide combinations and predict their ability to aggregate as a potential precursor to their self-assembly. The simulation protocol and scoring method proposed allows a rapid determination of whether a given peptide sequence is likely to aggregate (an indicator for the ability to self-assemble) under aqueous conditions. Systems that show strong aggregation tendencies in the initial screening are selected for longer simulations, which result in good agreement with the known self-assembly or aggregation of dipeptides reported in the literature. Our extended simulations of the diphenylalanine system show that the coarse-grain model is able to reproduce salient features of nanoscale systems and provide insight into the self-assembly process for this system.

Self-assembling peptide materials are currently widely studied for potential applications in biomedicine and nanotechnology.[1] They can form various nanostructures from 1D fibers[2] and tubes[3−5] to sheets[6] and 3D networks with tunable supramolecular functionalities such as electronic conductivity[7] and bioactivity.[8] These features, combined with their benign toxicology and ease of preparation, position peptide-based nanostructures at the forefront of next-generation, soft nanomaterials design. The high degree of flexibility in the types of nanostructures that can be produced is due to the vast chemical space that is available when considering the 20 gene-encoded amino acids and various functional groups that can be combined to create self-assembling peptide based units. As such, a method to rapidly screen these various combinations in the computer before investing time and resources in the laboratory is highly beneficial to the future development of novel materials.

Dipeptides, which consist of only two amino acids, have been chosen to apply our virtual screening procedure to, as they are the smallest of known peptide self-assemblers. They are also of interest in their own right and have shown interesting morphologies and functionalities.[9,10] The most famous and well-studied example is diphenylalanine (FF), as first described by Reches and Gazit in 2003[9] and recently reviewed by Yan et al.[10] FF is known to form tubular structures upon dilution from solvents into water or simply sonication in deionized water[11] but can be manipulated to follow different supramolecular order in other solvents.[12,13] Other dipeptide examples include FW[9] and IF,[14] which form tubular structures and fibrous hydrogels, respectively. In a systematic survey, Görbitz showed that many more dipeptides can form crystalline structures with microporous unit cells.[15,16] However, the formation of crystalline structures in saturated solutions does not necessarily imply self-assembly at low concentrations in solvent. Here we report a computational approach that allows us to determine, in an initial rapid screening, which peptides will exhibit aggregation in water. The aggregation of the peptides is a necessary precondition for their self-assembly, which can then be studied using longer, larger, and more detailed simulations.

For many amyloid peptide fragments, different simulation approaches have shown aggregation of individual oligopeptide fragments.[17−20] Wu and Shea have recently reviewed a number of coarse-grain methods for protein aggregation in the light of amyloid-related diseases and observed several mechanisms for fibrillization,[21] and progress was made toward understanding the general requirements for self-assembly.[22,23] In an effort to screen for peptide fibrillogenesis, lattice models have also been employed to determine contributions from sequence factors.[24] Of specific interest for this study, for the phenylalanine dipeptide and tripeptide, replica-exchange MD was performed, which showed the formation of open or ring-like peptide networks together with a significant β-sheet content.[25] However, these simulations are often limited to small systems sizes (e.g., 12 dipeptides were used in the replica-exchange study) or short simulation times due computational cost issues, although aggregation of self-assembling monomers often occurs on a microsecond time scale.[26] Villa et al. have showed that one can use implicit solvent and coarse-grained (CG) models to improve the accessible time frame without losing agreement with all-atom MD simulations.[27,28]

In the current study, we have employed the coarse grain MARTINI force field for biomolecular simulations[29,30] to model the aggregation ability for all possible combinations of two amino acids. This force field utilizes a 4:1 atom/CG-bead mapping to represent protein backbone and side chains (and a (2 to 3):1 atom/CG-bead mapping for ring systems). As an example, Figure 1 displays the diphenylalanine peptide structure and CG-bead representation of the same peptide. There is a significant loss of atomic detail in coarse-graining, and this removes the ability of the method to reproduce features such as explicit H-bonding interactions that fine-tune the relative positioning of one peptide to another, which is important, for example, in the formation of β-sheets. However, the polarity, shape, and nonbonded interaction potential of the amino acids are implicitly included through the parametrization of functional groups based on comparison with experimental results.[29,30] Therefore, although the detailed supramolecular structure is not available, coarse grain methods do provide insight into the driving force for the aggregation and consequently the self-assembly of peptides.

Figure 1

(a) Structure of the FF dipeptide and (b) MARTINI coarse-grain representation. Red beads represent the peptide backbone and white beads represent the side chain.

Utilizing the scripting tools in VMD,[31] we have analyzed the output configuration of every simulated system defining the aggregation propensity (AP) as the ratio of the solvent-accessible surface area (SASA) of the dipeptide molecules in the initial minimized box to the SASA of the final configuration of the simulation. A ratio greater than two is interpreted as a high degree of aggregation. From Figure 2, which shows the AP for all 400 dipeptide combinations, it is apparent that F and W residues promote aggregation, whereas small or charged residues decrease the score. In addition, we note that the graph is roughly symmetric. There are differences for some combinations, however, most notably SF and FS, where the AP of the former is 0.6 higher. (See the Supporting Information for details.) Most dipeptides have a score close to 1, indicating that they do not exhibit a propensity to aggregate, and thus their self-assembly in aqueous conditions can be excluded.

Figure 2

Two-dimensional grids indicating the AP score for dipeptides. Horizontal and vertical axes show the amino acid one-letter codes for the first (N-terminus) and second (C-terminus) amino acid, respectively.

The AP scores were compared to known experimental work on nine different dipeptides in aqueous solution, as summarized in Table 1. Although preparation methods vary throughout the literature, our data accurately predicts the self-assembly/aggregation, or lack there of, for those dipeptides for which experimental data is available. From this table, the data in Figure 2 and the final configurations of the MD simulations, we propose that an AP > 2 indicates a good candidate for further investigation into whether the system is likely to self-assemble. Clearly, this is an empirical, observation-based assignment and it will be very interesting to test more complex examples and those on the borderline (i.e., 1.5 ≤ AP ≤ 2.5) in an effort to break this rule. For example, HS is an interesting candidate with an AP score of 1.66 (see Supporting Information for the full list of AP scores) and does not contain any hydrophobic amino acids. Nonetheless, within the systems studied thus far the AP score provides a useful and apparently robust indicator of the desired property.

Table 1

Comparison of Self-Assembly Score with Experiments

dipeptide	reference	assessment method	observed structure	APa
FE	(32)	rheology	none	1.1
FF	(9−11,33)	TEM, AFM	tubes, vesicles	3.2
FK	(32)	rheology	none	1.2
FW	(9,34)	TEM	tubes, aggregates	3.5
IF	(14)	TEM, SEM	fibers	2.3
VF	(14)	TEM	none	1.8
WF	(9)b	TEM	aggregatesc	3.5
WW	(9)b	TEM	aggregatesc	3.2
WY	(9)b	TEM	aggregatesc	2.1

Aggregation propensity score is the ratio of the solvent-accessible surface area (SASA) of the dipeptide molecules in the initial minimized box to the SASA of the final configuration of the simulation.

Supporting Information.

Amorphous.

Comparison of the SASA for the initial and final configurations of the simulations provides a clear indication of whether aggregation occurs and as such whether self-assembly may occur. However, the AP score provides no details of the nature of the nanostructures. Globular, fibrous, and tubular structures have been observed among the results, but these features cannot be distinguished by this score. Work is continuing in our lab to identify a more suitable, generic descriptor for quantifying the supramolecular structure.

We have shown that the tendency for dipeptides to aggregate can be determined within a 400 ns simulation window on the relatively small (300 dipeptide) systems used in the virtual screening procedure. However, to test whether the CG model is able to converge (either qualitatively, or even semiquantitatively) toward the nanostructures observed experimentally, the simulation time and system size should be further increased to allow these structures to form. We have explored this possibility in the case of diphenylalanine by studying a larger system (1600 dipeptides in a box with double the dimensions of the original box, Figure 3a) to see whether the structures observed experimentally in an aqueous environment[9−11,33] are in agreement with the final structure from the simulation.

Figure 3

Snapshots at different time points in the extended FF dipeptide MD simulations. Red: backbone beads, white: side chain beads. Water beads are omitted. (a) 0 μs; randomly placed dipeptides. The periodic box (blue lines) is indicated in the Figure. (b) 0.2 μs; formation of sheet-like aggregates. (c) 0.5 μs; vesicle formation by sheet folding. (d) 1.5 μs; fused vesicles forming a hollow tube. The end of the tube is cut off to show water beads inside (blue).

In the final snapshot of the 1.5 μs simulation (Figure 3d), the FF monomers have generated a tubular nanostructure in agreement with that observed experimentally (Table 1). The 1D nanotube has an outer diameter of 60 ± 7 Å and an inner diameter of 29 ± 5 Å (Figure 3d). The size of the water core (inner diameter) of our simulated tube agrees reasonably well with the X-ray diffraction analysis by Kim et al.[11] and Görbitz[15] of crystallized FF nanotubes (van der Waals diameter of inner channel ∼10 Å). Moreover, we also note that the dihedral angle between the side chains rotates from 180° in the starting monomer structures (Figure 1) to an average value of 0° in the final tube like nanostructure (Figure 4). The rotation of the side chains away from the ideal 180° value for the dipeptide is also in agreement with the X-ray structures observed by Görbitz.[16] In the initial system (i.e., after minimization), the backbone dihedral angles are clustered around 120 and 180° (green columns in Figure 4). However, as the simulation continues, we observe that the side chains rotate during the aggregation process to adopt a dihedral angle between the side chains centered around 0°. After 1 ns of simulation, the majority of the dipeptides have already reached this equilibrium value (blue columns in Figure 4), although the distribution is not as sharp as is achieved in the final snapshot at 1.5 μs (red columns in Figure 4). Although a direct comparison of the distance and dihedral parameters is not possible because of the lack of atomic detail in the CG model, the similar size and observed rotation in both the experimental and simulated systems is compelling.

Figure 4

Side-chain–backbone–backbone–side-chain dihedral angle distribution (binwidth of 5°) at the start of the simulation (0 μs, green), after 0.1 μs (blue), and at the end of the 1.5 μs simulation (red).

In addition to the final tubular structure, the CG simulation also provides new information with respect to the mechanism by which the nanotubes are formed. From the initial random distribution of dipeptide monomers (Figure 3a), we observe an initial formation of sheet-like aggregates (Figure 3b) after only 0.2 μs. These sheet-like structures grow and accumulate to a point where they are then able fold up to form hollow vesicles after ∼0.5 μs (Figure 3c). Extending the simulation time to up to a total of 1.5 μs results in the continuous accumulation of these vesicles, which fuse together into 1D nanotubes (Figure 3d). Therefore, whereas the FF simulation shows that after 400 ns the final nanostructure that results from the self-assembly process may not be reached, the initial aggregation stages of self-assembly are already in process and as such the virtual screening protocol outlined above is effectively able to screen for the possibility of the system to self-assemble. Indeed, for dipeptides IF, WF, and WW, the final structures of their virtual screening simulations spanned one of the box dimensions to form fibrous assemblies, whereas FW also showed a tendency to form tubular structures, in agreement with the literature in Table 1. Results of these simulations can be found in the Supporting Information (Figure S1). The AP for all experimentally studied dipeptides after extended, 4 μs simulations (Figure S2; see the Supporting Information) is displayed in Table S1 (see Supporting Information) and does not show significant changes with respect to the values in Table 1, supporting the validity of the 400 ns value of the AP.

In conclusion, the MARTINI coarse-grain force field has been applied to predict the aggregation ability of all 400 dipeptides of the 20 gene-encoded amino acids. Good agreement was found comparing the predicted aggregation propensity and supramolecular structure to experimental results from literature for various dipeptides including diphenylalanine. From the FF simulation, it is apparent that equilibrium has not been reached after the 400 ns used in calculating the AP score; however, this is not the goal of the initial screening protocol, which is intended to identify systems that can potentially self-assemble and as such warrant further investigation. In all MD simulations, care has to be taken when considering simulation time and system size because some effects may be time- or concentration-dependent. The AP score after 400 ns provides a convenient, affordable, and apparently robust screening model for the ability of dipeptides to aggregate and hence their potential for self-assembly. For more accurate structural information and insights into the mechanism of formation of those systems identified as having a potential for self-assembly, longer and larger simulations can be performed to allow the system to converge.

The self-assembly of biomolecules often occurs on a time scale that exceeds the microsecond region, and as such atomistic methods are too computationally expensive to model multiple large systems. Moreover, as was observed in the FF simulation leading to the nanotube formation, there is a critical size requirement for the system to be modeled before self-assembly can be observed. That is, in the case of FF, we see initial sheet formation; the sheets are only able to form vesicles (the precursor stage to tube formation) once they has reached a sufficient size to fold on themselves. Atomistic simulations will not typically be of sufficient size to reach safely and surpass these critical size regimes because of the cost associated with such large systems at this level of theory. Thus, although detailed information may be lost, CG models are a very useful tool in the discovery of new biomaterials because they can easily be employed to screen a large amount of different molecules for their self-assembly properties. The results of the initial screening and subsequent production simulations for systems of interest can then form the basis for more detailed atomistic studies.

Computational Methods

A simulation for each dipeptide (in their zwitterionic form) was set up using the GROMACS molecular dynamics package.[35] A cubic box with 300 dipeptides, placed randomly with a minimum distance of 3 Å between them, was solvated in standard MARTINI CG water (four water molecules per bead) to a final concentration of ∼0.4 M. A Berendsen thermostat and barostat[36] were used to keep the temperature at 303 K and pressure at 1 bar, respectively. Bond lengths in aromatic side chains and the backbone–side-chain bonds in I, V, and Y were constrained using the LINCS algorithm.[37] All boxes were energy minimized using the steepest descent integrator and then equilibrated for 4 × 106 time steps of 25 fs.[38−41] The extended FF simulation was performed in duplicate, and reported tube sizes were averaged over 10 different positions in either simulation. For all experimentally studied dipeptides, a similar simulation over 40 × 106 steps was performed to study the effect of a longer simulation on the AP.

The total screening simulation time equates to 100 ns, but because of the smoothness of the CG potentials, this roughly equates to an effective 400 ns of atomistic simulation time.[29,42] Throughout this Letter, all times refer to the simulation time speeded up by this factor of 4. The size of our systems, the time scales required, and the number of systems that are assessed in a virtual screen excludes the possibility of performing these MD simulations using all-atom models on modest computing resources. The complete, 400 ns simulation for a single dipeptide system described above can be carried out in ca. 4 h on a four-core 2.6 GHz Opteron compute node. More details of the simulation can be found in the Supporting Information.