Cocaine is a member of tropane alkaloids family, one of eight stereoisomers of 8-methyl-8-azabicyclo (3.2.1) octane-2-carboxylic methyl ester [1]. This compound has many physiological effects such as anesthetic [2], vasoconstriction [3], heart rate and blood pressure increaser [4], stimulant of sympathetic nervous system [5] and others, this alkaloid has been extensively used in medicine [6-9]. Health and treatments monitoring, diagnosis of overdoses, forensic investigation and other disciplines are interested in cocaine research [10-15].

In this work, the application of a molecular modeling method to study peptides selection and cocaine–peptide interaction was evaluated with the aim of introducing a streamlined technique to find new biomimetic receptors for future application in biotechnological fields.

The use of molecular modeling software implemented in a virtual screening methodology has offered a considerable assistance for the rational design of molecular systems with high selectivity and affinity, also experimentally proved [16-21]. In the past decade laboratorybased researches between simulated binding scores obtained from a database of candidates and experimental results were reported, by showing the ability of synthetic receptors to bind a target, such as small peptides [20,22], small organic compounds [23-26] and DNA [27,28]. In particular the docking program FRED [29] used in this work was already proved to calculate binding scores in line with other software and with experimental data [17,20,22,28].

The premise of this work was to obtain synthetic peptides by mimicking the amino acid residues geometry of known biological protein receptors active site of a given target compound like cocaine. These amino acid residues were then combined to form polypeptides in a semi-combinatorial way in order to decrease, by orders of magnitude, the number of combinations to be screened, and lowering significantly the molecular complexity in terms of sequence length and threedimensional structure.

The goal was to obtain a series of viable synthetic molecular traps in a short period of time using modest computational resources. These possible candidates obtained can be proposed especially for biotechnological purposes reducing the number of experimental laboratory trials and therefore lowering research costs.

The methodology is based on the assumption that small peptides can be designed by initially mimicking cocaine binding site of larger biological structures and then increasing the cocaine-peptide affinity by addition or mutation of amino acids in the sequence.

The virtual screening approach used here was based on the idea of making a rough selection of small peptides and filter them afterwards, thus allowing the analysis of large amounts of structures by reducing computational time. The peptides sequence length was restricted to 5-6 amino acids long in order to prevent the appearance of secondary structure motifs that increase the system complexity, extending the calculation time, thus decreasing the number of possible structures to be analyzed. Understanding the intermolecular interaction between peptide-cocaine binding affinity was very useful to confirm the functionality of these biomimetic receptors for successive biotechnological applications.

This approach is spreading because of the necessity to use cost- and time- effective ways to identify lead compounds [30-35]. However, such approach is not exempt of problems and results are not always reliable. Incorporating in virtual screening calculations, the dynamic nature of receptor structures is a well-known problem [36-39]. In many common molecular docking algorithms, target protein geometry is kept rigid in a single low-energy conformation, and only conformational and positional flexibility of a ligand is considered [40,41]. Small peptides, however, can have different conformational states with similar energies [42-45]. In many cases, binding site conformation of a receptor exhibits significant motion including rearrangements of side chains and backbone upon ligand binding. Using a single receptor conformation in docking experiments can lead to errors in identification of binding modes and errors in prediction of binding affinities [46]. This can significantly reduce the chances of finding new ligands.

To overcome this issue, here, the shape variation was considered for both peptides and cocaine target. The peptide flexibility was studied generating their conformers and determining the minimum number of them required to make a statistically significant representation of each peptide conformational space. Computing mobility for both receptor and ligand assured a more reliable computed binding score. Finally, the in-silico results were tested in a preliminary in vivo experiment, using two different computationally designed hexapeptides as selective SPE sorbent for cocaine.

Software used

All calculations and molecular modeling experiments were performed using a desktop PC with a 3.4 GHz Intel Core I7-2600 processor having 8 GBytes DDR3 RAM with 1333 MHz bus, running Microsoft Windows 7 Professional 64 Bits.

All structures were designed and cleaned up with Hyperchem 8.0.5. The interactions between ligand and protein atoms were determined with RCSB PDB Ligand Explorer, from Protein Data Bank (PDB). In order to find sequence homology amongst proteins ClustalX was used.

OpenEye Scientific Software package under academic license were used. The energy minimization process was carried out using SZYBKI 1.5.1 in its default parameterization [47]. OMEGA 2.4.3 was used to generate conformers to both cocaine and peptides. The program was used with MMFF as the force field [48]. The docking software FRED 2.2.5, was used to conduct the virtual screening applying default parameters [29]. VIDA 4.1.1 was used for visualization, postcalculations analysis and representation [49].

The entire process was scripted, automated and executed using AutoIT V3, a freeware BASIC-like scripting language.

Peptide scaffold

Peptides were designed in zwitter ionic mode, using only the 20 natural amino acids.

Amino acid residues participating in cocaine interactions were identified using the RCSB PDB Ligand Explorer. The cut off distance (H bonds, electrostatic, hydrophobic, etc.) was set to 4Å as a maximum from any cocaine atom. A multiple sequence alignment was done with ClustalX in order to find amino acids homologies in cocaine-protein binding sites. The most important conserved amino acid residues interacting with cocaine in protein binding sites were selected for building up the first pentapeptides library.

Peptide libraries: screening process, docking and structural analysis

The screening process was divided in 4 steps. In each step a peptide library was generated by considering the peptide conformers flexibility, or introducing mutations in selected peptide sequences:

1. A first library of 768 pentapeptides, without conformers, was designed using as pivots the most conserved amino acid residues, found in the 4 proteins binding cocaine.

2. A second library of 1050 pentapeptides was generated by using 10 conformers for each structure of 105pentapeptides from the first library, selected on the basis of the binding score: the 55 pentapeptides having the better scores vs cocaine; 25 pentapeptides with scores close to average; and 25 pentapeptides with the worst scores.

3. A third library of 3000 hexapeptides was then generated, without conformers, using as backbones the 25 higher ranked pentapeptides from the second library, inserting in each position of the sequence the 20 natural amino acids one by one.

4. A fourth library of 1450 hexapeptides was designed and docked by using 10 conformers for each structure of 145 hexapeptides from the third library, selected on the basis of the binding score: the 115 hexapeptides having the better scores vs cocaine; 15 hexapeptides with scores close to average; and 15 hexapeptides with the worst scores.

The first and third libraries were roughly screened without considering peptide flexibility. A calculation refinement was considered in second and fourth libraries, including conformers for each structure and therefore considering the shape variation and flexibility of peptides. In all cases every library was obtained from results obtained in previous docking run.

Boxes defining the active site were generated for each peptide. The box size was comprised from 4500ų to 7200 ų ranged in the 90% of all cases. The peptide was inside the box considering the whole structure as possible binding site for cocaine. The time required for each peptide structure, from the initial design to final docking, was about 5 minutes.

In-vivo experimental testing

Peptide SPE sorbents were from EspiKemsrl (Sesto Fiorentino, FI, Italia), QHWWDW-resin (Nova Syn TGA), with a peptide substitution level of 0.17 mmol g^-1 and ESSIDH-resin (Nova Syn TGA), with a peptide substitution level of 0.20 mmol g^-1. The cartridges (volume 1 ml) were packed with 30 mg of modified peptide resin dissolved in 5 mL of an ethanol/water solution (80:20, v/v). The cartridges were conditioned and equilibrated by washing with ethanol. The extraction procedure was performed in four steps: -Conditioning of the stationary phase with Tris-HCl (pH=7.5); -Sample loading (1 mL); -Washing with 1 mL of ultrapure water; -Elution with 1 mL of formic acid 5 mM in methanol, this fraction was analyzed by LC-MS/MS as reported in details in another work along with all chemicals used [50].

Peptide scaffold

Cocaine search in PDB web site returned 52 structures. These structures were then filtered by ligands and redundancies, giving a final set of 3 antibody-cocaine complexes (murine-human chimeric Fab of GNC92H2, anti-cocaine monoclonal antibody M82G2 and antibody 7A1) [51-53] and 1 acetylcholine-binding protein (AChBP) -cocaine complex [54].

These 4 complexes were used for conformational search of protein binding pockets, identifying similarities between structures to determine conserved residues. The type of interactions present in the complexes was also identified. The sequence alignment carried out with ClustalX showed that zones binding to cocaine were not conserved. All active sites found involved 2 or 3 chains with different geometrical arrangements. Figure 1 shows an example of highly conserved chains and low homology amongst sequences. In both cases the residues participating in the interaction did not match. Articles referencing PDB structures confirmed this result [51-54]. Most of interactions were hydrophobic, stacking and hydrogen bonds, although the later were few.

Figure 1: ClustalX sequence alignment using GONNET series between light chains (top) and heavy chains (bottom) from anti-cocaine antibody M82G2 (1Q72) and cocaine catalytic antibody 7A1 FAB' (2AJV). The stars (*) denote a position having a single, fully conserved residue in the chains.

According to chemical function and interatomic distances the cocaine molecule was subdivided in 5 regions. As reported in figure 2 the zones were the benzene ring, the ester directly attached to benzene, the ester fragment attached to aliphatic ring, the amine group and the cyclo-heptane ring.

Figure 2: According to chemical function and interatomic distances the cocaine molecule was subdivided in 5 regions: aromatic R1 (benzene ring), ester R2 (the ester directly attached to benzene), ester R3 (ester fragment attached to aliphatic ring), amine R4 and aliphatic R5 (the cycloheptane ring).

All interacting amino acid residues throughout the 4 analyzed cocaine complexes were referred to these regions. The interaction of cocaine with amino acid residues in the 4 complexes was reported using number of occurrences and the percentage ratio respect to all contacts found in the complex region. The results were summarized in table 1. Tyrosine with 40% and tryptophan with 20% of total contacts were the amino acids more represented, highlighting the aromatic residue group significance in the cocaine binding interactions.

Table 1: The interaction of cocaine molecule with amino acid residues in the 4 biological complexes. In bold-italic the aminoacids selected for generating the first 768 pentapeptides library. Res.=amino acid residue; #=number of contacts; %=percent of contacts in cocaine region.

According to the data of table 1, the most relevant amino acid residues for each region were chosen. The minimum sequence length of initial peptide was set to 5 amino acids long. This chain length allowed peptides to include in their sequences at least one of the most representative amino acids for each cocaine region, reproducing a significant part of interactions found in biological association complexes.

In order to reduce the number of structures to be tested, the first 768 pentapeptides library was generated using only the first 3 amino acids for first position, comprised over 50% of all interactions for this region by representing aromatic, polar and negatively charged amino acids. All 4 amino acids were selected in second and fourth positions. The first 3 amino acids plus serine in third position; serine was chosen because it was the only polar amino acid among the ones with lower frequencies in that region. The first 4 amino acids for the fifth position were considered, discarding leucine for his aliphatic properties.

Using the combination of these amino acids in their respective positions the pentapeptides library resulted in 768 elements decreasing in four orders of magnitude the number of combinations to be tested in a combinatorial approach.

Peptide libraries: screening process and docking

The docking results for the 768 pentapeptides are shown in figure 3. The lines intersection marked a strong change of the variation in the docking score, used as cut-off for the selection of the 55 pentapeptides forming the more stable complexes with cocaine, which corresponded to the 7.2% of all structures in the library. The score values were calculated using the consensus of multiple standard scoring functions, where lower values represented higher peptide-cocaine affinity. The structures of the higher ranked 55 pentapeptides had as most representative amino acids, in first position, tyrosine with 40% of occurrences, in second position, arginine with 56% of occurrences in third, fourth and fifth position and tryptophan with 53, 36 and 36% respectively. These results highlighted the presence of aromatic and charged amino acids motifs in the better scored structures.

Figure 3: The docking results of the 768 pentapeptides vs cocaine. The lines intersecting the score -13.0 were used as cutoff for the selection of the top ranked 55 peptides corresponding to the 7.2% of the library.

From the first 768 pentapeptides screening, four pentapeptides ranked in position 2^nd, 38^th, 360^th and 760^th were selected to study peptides flexibility, generating 10, 20, 50, 100 and 200 conformers for each structure. This conformational study was carried out to determine the minimum number of conformers required to make a statistically reliable computed score. The pentapeptides were chosen arbitrarily but within peptides group showing an overall trend in terms of score values, similar sequence patterns and amino acid composition. In figure 4, the results were expressed in percentage.

Figure 4: Conformers analysis using4 pentapeptides representing different positions in the 768pentapeptides library. Average binding scores of 1, 10, 20, 50, 100 and 200 pentapeptide conformers vs cocaine were reported in percentage. CFM=conformer; R=docking ranked position.

The standard deviations ranged from 10 to 18% depending on the amount of conformers. According to these results, 10 conformers were chosen for further docking studies because the overall tendency remained unchanged when using more conformers. The divergences observed between the score average using 10 conformers and the others are less than 10% of variance in the worst case. Because a significant divergence was observed comparing the scores obtained with or without conformers even if the trend remained unchanged, to reduce the false positives and/or false negatives, 105 pentapeptides were selected from the 768 pentapeptides library. Ten conformers were generated for each of the 105 pentapeptides, resulting from the fact that the conformational space flexibility of the peptide was not considered in first library. The selected peptides included 55 best ranked structures (7.2% of the library), 25 structures with scores close to global average and 25 structures having the worst scores. Considering the conformers average score, 78% of 55 structures having the lower conformers average scores, were comprised within the 55 best structures selected in the first screening step without conformers, 13% corresponded to the ones having the intermediate scores and just 9% had the worst scoring. Typical standard deviations were within 17% margin in 95% of the cases. Aromatic residues remained the most conserved residues. Considering the first 25 highest ranked structures reported in table 2, in the first position, tyrosine residue had 52% of occurrences, tryptophan appeared 64%, 76%, 32% and 40% of the cases in the 2^nd, 3^rd, 4^th and 5^th position respectively. The general tendency continued with charged and aromatic amino acids the most frequent in all highest ranked sequences.

Table 2: The 25 pentapeptides with the lower conformers average score selected to generate the 3000 hexapeptides library. The pentapeptides score and rank position in the first 768 library run using one conformer was also reported. AV=average; SD%=relative standard deviation.

To increase the possibilities of attaining more stable complex with cocaine, the 25 pentapeptides with the higher ranking average scores (Table 2) were selected to generate a 3000 hexapeptide library by inserting alternatively in each sequence position the 20 natural amino acids. Then, for a conformer screening refinement 145 hexapeptides were selected from this 3000 hexapeptide library: the 115 best ranked peptides, 15 peptides with scores close to average and the 15 worst ranked peptides. For each hexapeptide, 10 conformers were generated and then docked with cocaine. In figure 5 the conformers average scores of the 145 hexapeptides selected from the 3000 hexapeptides library.

Figure 5: Binding average score trend of the 145 hexapeptides vs. cocaine calculated using 10 conformers for each hexapeptide.

Results showed that trend scores remained sigmoidal like in the first 768 pentapeptides library behavior (Figure 3), having 9% with a score lower than -10 and 6% higher than zero. In 95% of cases, typical standard deviations were within the 15% margin. A decrease of almost 4 units was observed comparing the lowest conformers average score between hexapeptide (-13.7, shown in figure 5) and pentapeptide (-9.4, shown in table 2), revealing that, potentially, hexapeptides formed stronger complex with cocaine than pentapeptides.

The 20 best ranked hexapeptides, with conformers, were 100% from the best ranked hexapeptides without conformers. When considering the 50 best ranked peptides the matching percentage decreased but only to 93%.

Tyr in the first position and Trp in the second, third and fourth position were always amongst the highest ranked structures. In these structures the W-W sequence starting from positions 2 or 3 were the much conserved sequence (76% of the cases). When aliphatic or uncharged amino acids were in sequence the ranking decreased dramatically.

Structural analysis and in-vivo experimental testing

The peptides were built taking into account the amino acid residues of protein binding pockets, but reducing the size to the minimum, in order to control the possible shape during the computational simulation.

The (Table 3) shows the structural analysis reporting the occurrence percentage in each structure position of the 73 higher ranked hexapeptides. The hexapeptides structure analysis was reported considering the higher ranked structures having scores less than -6 and corresponding to the 50% of the 145 hexapeptides analyzed in the last screening step.

Table 3: Structural analysis of the 73 higher ranked hexapeptides having scores less than -6 and corresponding to the 17 of the 145 hexapeptides considered. P=position in sequence; =occurrences percentage in hexapeptides.

The correlation between amino acids in best ranked and those in binding pockets were not very strong, the percentages were quite different and ranks of amino acids changed often. Tyr in first and last position was the most represented with 58.9% and 31.5% of occurrences. The Trp in second and third position was found with 52.1 and 83.3% of occurrences respectively, confirming the strong presence of aromatic residues. In fifth position a remarkably high presence of Asp was observed. The distribution of aliphatic and uncharged residues was rare in higher ranked sequences.

According to chemical function, binding score and interatomic interactions with the 5 cocaine regions, two hexapeptides were selected for a preliminary in vivo experiment. The peptide QHWWDW having a 10 conformers average score of -9.6 in binding cocaine, and the was selected because, in simulations, the complex formed with cocaine was, with or without conformers, always amongst the 20 best ranked structures.

On the other hand the ESSIDH hexapeptide, having a 10 conformers average score of -6.1 in binding cocaine, derived from pentapeptides which were ranked amongst the worst structures for both single and 10 conformers simulations.

Moreover considering the primary structure, the hexapeptides QHWWDW and ESSIDH had, respectively, the higher and lesser preserved amino acids residues participating in the cocaine interaction, identified in the 4 biological receptors.

Considering the two peptides response vs. cocaine, the default parameters used by FRED were quite in agreement with the SPE experimental results as reported in table 4.

Table 4: Experimental and simulated data comparison; the binding scores were reported also in percentage considering the results of the 145 hexapeptides library.

In term of the shape, the three dimensional hexapeptide structure in docking cocaine diverged widely from one to another, but despite the differences, cocaine was always in the center of the hexapeptide molecular surface, especially versus those having better binding scores. Figure 6 shows the simulated electrostatic molecular surfaces of the two hexapeptides used in the preliminary in-vivo experiments. Both, hexapeptides were able to bury cocaine in their electrostatic molecular surfaces, but the area was much larger for QHWWDW than ESSIDH. Reflecting the difference in binding score. Anyway, it should be noted that the in-vivo interaction can be due to a synergic cooperation of the amino acid residues that have a certain amount of freedom to move around the carbon backbone (much larger than in proteins), increasing the probability to interact with cocaine, especially when cocaine is found in the center of the complex.

Figure 6: Electrostatic molecular surfaces of the two experimentally tested hexapeptides in complex with cocaine.

Using a modest computing framework, it was possible to obtain a suitable set of peptides binding cocaine. The higher ranked peptides maintained some specific sequences in concordance with results obtained in biological proteins. According to our findings, sequences rich in Trp and Tyr gave high affinity for cocaine, but negative impact was observed in the presence of aliphatic or uncharged residues. Preliminary in vivo experimental results were in agreement with simulated data. Despite the results, higher ranked peptides were not supposed to be the best in binding cocaine, but the methodology proposed can be used in support to experimental tests, rationalizing and reducing by orders of magnitude the choice of molecular traps.

The authors thank the financial support of FP7-PEOPLE-IRSES project NANOSENS 230815.