In silico Evaluation of Cucurbit[6]uril as a Potential Detector for Cocaine and Its Adulterants Lidocaine, Caffeine, and Procaine

Illicit drugs and their trafficking require worldwide efforts in investigation, detection, and control. Colorimetric tests are often applied to identify drugs. Cocaine has some well-known adulterants that can provide a false positive response. Cucurbit[6]uril (CB[6]) has been suggested as a potential detector for cocaine and other illicit drugs. This work uses in silico methods to evaluate the use of CB[6] to detect cocaine and these interfering substances. More specifically, this work analyzes different possibilities of CB[6] complexation with cocaine, lidocaine, caffeine, and procaine and compares the results achieved for cocaine and its adulterants. Different methodologies were employed: quantum chemistry was investigated through DFT B3LYP/TZVP (density functional theory-Becke, three-parameter, Lee-Yang-Parr with triple zeta valence plus polarization basis set) and the semi-empirical methods Austin model 1 (AM1), parametric methods 3, 6, and 7 (PM3, PM6, PM7), and Recife model 1 (RM1). We used these methodologies intending to compare the reasonability and reproducibility of the results in the gas phase condition. Solvent influence was studied by molecular dynamics (MD) simulations. Results showed that CB[6] does not bind to these substances, as judged from the positive values of binding free energy obtained with all methods. DFT and MD were the most reliable methods whereas semiempirical ones were not reproductible in describing these systems. Results also showed that interactions are not specific, so CB[6] does not provide a good response for cocaine detection.


Introduction
Drug trafficking and consumption are worldwide concerns, and integrated efforts are needed to combat them. Among illicit drugs, cocaine has a potential supplydriven market in expansion, while coca bush cultivation has increased concomitantly with the global drug manufacture. 1 Correct drug identification is essential to ensure human rights: it can prevent the arrest of innocents and avoids the miscarriage of justice. 2,3 It is also important in harm reduction projects because there are advantages in testing drugs before they are consumed. Many harm reduction agencies are exploring techniques for illicit drug testing to identify and, when possible, to quantify their constituents, thereby allowing their users to make informed decisions. 3,4 Drug characterization comprises two steps: presumptive (in situ) and confirmatory (in lab) tests. Presumptive tests are less time-consuming and less expensive. However, most of these tests consist of colorimetric tests and can give doubtable results. 5 The Scott test is the most widely employed for cocaine identification in law enforcement actions. Nevertheless, it can return false positives because cocaine is not sold in its pure form. On the other hand, confirmatory tests are more trustful, despite being more expensive as they require specialized equipment, trained examiners, and complex sample preparation, besides they

In silico Evaluation of Cucurbit[6]uril as a Potential Detector for Cocaine and Its
Adulterants Lidocaine, Caffeine, and Procaine can be destructive. 2,6 Therefore, new materials and methods are constantly being searched in an attempt to obtain trustful identification at a low cost. [7][8][9][10][11] Cucurbiturils (CB) can form host-guest complexes with various drugs consisting of small organic or inorganic molecules. 12 Hydrophobic effects facilitate drug encapsulation within the cucurbituril cavity, and further hydrogen bonding or ion-dipole interactions with the cucurbituril portal stabilize the resulting complex. 13 The defining structural features of the CB[n] family are their highly symmetric structure bearing negatively charged carbonyl rims and a hydrophobic cavity. 14,15 Cucurbituril (cucurbit [6]uril, or CB [6]) is a hexameric macrocyclic compound originating from self-assembly during an acid-catalyzed condensation reaction between glycoluril and formaldehyde. Curcubituril is commonly abbreviated as CB [6] or even CB6, Q [6], Q6, or Cuc6 in some cases; '6' represents the number of glycoluril units in the macrocycle. 16,17 The CB[n] cavity can recognize amines and azaheterocyclic compounds via hydrophobic effects, size effects, and ion-dipole interactions at the CB[n] portals. 18,19 The rigid structure and the ability to form stable complexes with molecules and ions make CB [6] an attractive building block to construct supramolecular architectures. 20 CB [6] can also facilitate cy c l o a d d i t i o n . CBs have been used in drug identification studies for opioids, marijuana and cocaine, 3,4-methylenedioxymethamphetamine (MDMA) and amphetamine-type stimulants. 14,[21][22][23] Previous in silico methods based on the semi-empirical methods PM7 (parametric method 7) 22 and RM1 (Recife model 1) 23 have been applied to study cucurbituril complexes. GAFF (generalized Amber force field) have also been used to examine how acyclic cucurbituril (aCBs) and opioids interact. 21 The DFT-B3LYP (density functional theory-Becke method, three-parameter, Lee-Yang-Parr) hybrid functional and the 6-31G(d) basis set have been applied to study reactive blue 19 dye adsorption onto CB [6] and CB [8]. 24 Microsecond time scale molecular dynamics simulations have aided the investigation into the binding enthalpies of CB [7] with different guests in aqueous solution. 25 Adsorption of hydrogen sulfide by CB [7] has been studied by ab initio van der Waals density functional (vdW-DF) calculations. 26 This work aimed to employ in silico methods to evaluate the CB [6] ability to form a complex with cocaine and some of its adulterants: lidocaine, caffeine, and procaine. Quantum chemistry and molecular dynamics approaches were used. Our intention was to verify the differential energy involved with CB [6] and the molecules of interest as an indicator of the interaction between these compounds.

Methodology
In silico methods include several methodological approaches. A combination of theoretical and experimental work is usually used to find agreement among them. Despite the importance of experimental evaluation, computational methods can provide valuable and detailed information that is not possible to obtain directly from empirical observation. 27 However, the theoretical tool must provide an accurate result to reproduce experimental findings. That is why we decided to carry out these analyses by combining different methods: semiempirical, DFT and MD (molecular dynamics) simulations.

Quantum chemistry study
Crystallographic structures of the molecules studied herein were obtained from The Cambridge Structural Database (CCDC). [28][29][30][31] Figure 1 presents the studied system. The CCDC structures were previously optimized with UFF (universal force field) in the Avogadro v.1.1.1 32,33 software. This force field can optimize the geometry of all elements and can properly handle inorganic materials and organometallic materials. 34 After optimization with UFF, the molecules were submitted to semi-empirical optimization with the MOPAC software by applying five different methods: Austin model 1 (AM1), parametric methods 3, 6, and 7 (PM3, PM6, PM7), and Recife model 1 (RM1). [35][36][37][38][39][40] All the molecules were reoptimized with the hybrid DFT-B3LYP 41 with a single zeta Ahlrichs basis set, and then single point refinement was carried out with a triple zeta Ahlrichs basis set (TZVP). 42 All calculations were performed by ORCA v4.0.1.2 software. 43 To characterize the stationary state, we have used the crystallographic structures for all molecules as a starting point for further optimization. No imaginary frequencies were found, indicating the achievement of a minimum energy structure for each molecule. [44][45][46][47][48] The thermodynamic functions and the final energies of the complexes were obtained. 24 To evaluate the tendency of the molecules to form a complex with CB [6], the hypotheses listed in Table 1 were investigated from a theoretical viewpoint. The table also defines the respective abbreviations. In all cases, the energy tendency was evaluated by the binding energy (∆E bind ) equation 1: 49,50 ∆E bind = E(complex) -(E(CB [6]) + E(Mol)) where Mol = cocaine, lidocaine, caffeine, or procaine.

MD simulations
Molecule parametrization CB [6] and the studied drugs were parametrized according to a protocol described in more detail elsewhere. 51 Briefly, the molecules were subjected to structural optimization followed by single-point energy calculations in Gaussian 09; the HF/6-31(d) method and the Merz-Kollman scheme were employed. 52 Point charges were derived through the restricted electrostatic potential (RESP) fit available in the antechamber program from Amber 16 suite. 53 The molecules were then fully parametrized with the aid of the GAFF. 53,54 MD simulations and calculation of binding enthalpies For each CB[6]-drug system, different manuallygenerated conformations were subjected to MD simulations to determine their time stabilities and their respective binding enthalpies (∆H bind ). Each studied system was embedded into a cubic box and solvated in methanol to match the experimental conditions. 23,55-60 The box size was adjusted so that the minimum distance between the box walls and the solute surface was 10.5 Å and the number of methanol molecules was exactly 630 in all cases. The methanol explicit solvation model available in Amber 16 was used to conduct the MD simulations. 53,61 No counterions were added because the solute molecules were all neutral. Then, the systems were subjected to a sequence of energy minimization (EMs), equilibration, and production runs according to a modified version of a protocol reported by Fenley et al., 25 which allowed the ∆H bind values to be calculated through a multibox strategy. In addition, some other adaptations were made during the equilibration stages. Differently from the tutorial, isothermal-isochoric (NVT) equilibration was accomplished for 500 ps, while the solute heavy atoms were kept restrained with a spring constant (k) of 10 kcal Å -2 mol -1 . Heating was carried out by linearly increasing the temperature from 10 to 298.15 K; this final value was kept constant in all subsequent simulations. A subsequent 500-ps isothermal-isobaric (NPT) equilibration was performed by applying the same previous restraints on the solute heavy atoms, which were then decreased  Position Abbreviation Cocaine and CB [6] Cocaine in the center COC_P01 Amino group oriented toward the center COC_P02 Aromatic ring oriented toward the center COC_P03 Ester carbonyl group ester oriented toward the center COC_P04 Lidocaine and CB [6] Lidocaine in the center LID_P01 Amino group oriented toward the center LID_P02 Aromatic ring oriented toward the center LID_P03 Caffeine and CB [6] Caffeine in the center CAF_P01 Five-member ring oriented toward the center CAF_P02 Aromatic ring oriented toward the center CAF_P03 Procaine and CB [6] Procaine in the center PROC_P01 Amino group oriented toward the center PROC_P02 Aromatic ring oriented toward the center PROC_P03 CB [6]: cucurbit [6]uril.
down to k = 2 kcal Å -2 mol -1 in a stepwise fashion during three consecutive 500-ps NPT equilibrations (k = 8, 6, and 2 kcal Å -2 mol -1 , respectively). Next, in accordance with the tutorial, 40-ns unrestrained NPT equilibrations were conducted for all the systems to determine the respective average box sizes. Finally, 1 μs NVT production runs were accomplished for the stable systems after their box sizes were modified to match the respective average volumes calculated from the precedent 40-ns MD simulations. Three replicate MD simulations were run for each CB [6]-drug conformation by following the same sequence of steps previously described, and different random velocities were assigned during the heating process. All the EMs and heatings were carried out with Amber 16 pemd.MPI, while NPT equilibrations and production runs were conducted with Amber 16 pmemd.cuda. 53 Treatment of nonbonded interactions and temperature and pressure control, among others, were identical to those employed by Fenley et al. 25 For the free drugs and CB [6], the average blocking approach was applied to estimate the standard error of the mean (SEMs) of their average potential energies (E ptot ) required for ∆H bind calculations, as only one productive run was carried out for these systems. The gmx analyze program available in GROMACS v.5.1.3 was employed for the previous analyses. 62 The E ptot mean values for the CB[6]-drug poses were obtained by averaging the results for the replicate MD simulations and using only those frames in which the drug was encapsulated. SEMs were estimated using the following formula: SEM = SD/√3, where SD is the standard deviation of the mean values in the three replicate MD simulations. Finally, the errors were propagated in order to calculate the SEMs of the reported ∆H bind values.

Free energy calculations
The attach-pull-release (APR) protocol reported by Henriksen et al. 63 was used to calculate the binding free energies (∆G bind ) of the studied CB[6]-drug complexes that remained stable during, at least, one of their respective 1 μs MD simulation replicates. Last frames generated during such stable trajectories were employed as the starting structures to perform the APR simulations, which were automatically carried out by means of an updated set of python scripts that enable the use of various non-aqueous solvents, including methanol, during the system setup. Except for the solvent, the MD simulations were conducted by the default inputs of the APR protocol in most cases. The solvent input was set to methanol, and 900 molecules of the latter solvent were added to each solvation box. The production runs of each window were simulated for up to 40 ns, the first 10 ns being discarded from subsequent free energy calculations. Hydrogen mass repartitioning (HMR) was employed to increase the simulation time step to 4 fs. Radial conformational restraints, i.e., jacks, were also set at the entrance of the CB [6] cavity to allow smooth exit of the encapsulated drug molecule. 63 Three pairs of oxygen atoms (O,O) and six pairs of nitrogen atoms (N,N), all lying in radially opposite positions at the entrance ring where the drug was placed, were chosen to define the jacks. In all cases, the jack distance and the force constant were set to 10 Å and the default value (25 kcal Å -2 mol -1 ), respectively. The complexes with a six-member ring inserted in CB [6] were the exception: the jack distance was set to 11.5 Å due to the wider molecular shape of the enclosed moiety. The free energies and their uncertainties were determined by the thermodynamic integration (TI) approach and blocking analysis, respectively. The equilibration and production runs were performed with Amber 16 pmemd.MPI and pmemd.cuda, respectively. 53

Trajectory analysis and visualization
Root mean square deviations (RMSDs) for the drug heavy atoms were calculated throughout the concatenated heating, NPT equilibration and production runs for every trajectory, by fitting all frames to the CB [6] heavy atoms and taking the minimized starting structure as reference. As a complementary way to assess complex stability, the solvent accessible surface area (SASA) of CB [6] was monitored for every frame during all MD simulations of CB[6]-drug poses. The bondi set of van der Waals atomic radii and the MOLSURF program available in Amber 16 were employed in such calculations, which, together with those of RMSD values, are integrated into the cpptraj module. 53 All figures depicting the molecules were created with PyMOL v2.1.0. 64 Table 2 summarizes all the results for the calculations performed with DFT methods for those conformations that achieved convergence. No relevant information was obtained from semi-empirical methods because they did not reach convergence for most conformations. Semiempirical results are shown in the Supplementary Information (SI) section. We observe that the DFT calculation returned positive values for the Gibbs free energy, evidencing that encapsulation was not spontaneous even when exothermic values were obtained for enthalpic contribution.

Molecular dynamics
We defined each starting conformations of the analyzed systems according to Table 1. In this section we will assess the stability of the encapsulated drugs by using MD simulations and thermodynamic calculations based on the generated ensembles. Figure 2 shows the outcomes of the MD simulations conducted for the different poses of cocaine in complex with CB [6]. Note that poses COC_P01 and COC_P03 converged to equivalent conformations during their respective MD simulation replicates that remained stable up to 1 μs. DFT calculations evidenced the same behavior (Table 2 and Figure S1, SI section). In fact, for the CB[6]cocaine complexes, the only binding mode that remained stable in at least one long MD simulation replicate was the mode bearing the phenyl group inside the host cavity. In the three remaining poses, cocaine either dissociated from CB [6] or converged to the same previous stable structures.
MD simulations carried out for different initial poses of the CB[6]-lidocaine complex identified a stable binding mode up to the completion of the simulation time in some replicates, while others showed the same binding mode dissociating from the host. Figure 3 depicts structural representations of the initial (t = 0) and final frames of the MD simulations in each case, together with the SASA of the CB [6] molecule and RMSD time profiles calculated with respect to the initial position of the drug heavy atoms. The insets in certain graphs demonstrated that RMSD values varied widely after restraints applied during equilibration were removed. Out of the three poses analyzed in this work, two of them (LID_P01 and LID_P02) converged to the same final structures in certain replicates, in which the diethylamine moiety of lidocaine was inserted in CB [6]. DFT calculations revealed similar final convergence to LID_P01. Conversely, the lidocaine pose in which the aryl group was placed into CB [6] readily dissociated during a shorter 42-ns unrestrained MD simulation (LID_P03). Figure 4 displays the results for caffeine. In this case, we identified a stable conformation of the CB[6]-caffeine complex in at least one MD simulation replicate, in which the drug's six-membered ring was inserted in CB [6], whereas the other two poses were unstable during their respective MD simulations (Figure 4).
On the other hand, procaine presented two stable conformations according to Figure 5, each of them bearing a different group lying within the CB [6] cavity. Note that poses PROC_P01 and PROC_P02 converged to equivalent conformations during all their respective MD simulation replicates. The diethylamine moiety inserted in CB [6] was obtained from two different starting poses during independent MD simulations.
After analyzing stability, we selected the stable conformations of each host-guest system to conduct thermodynamic calculations. We subjected five complexes to the APR protocol to calculate their respective ∆G bind values (see Figures S3 and S4, SI section). In addition, we accomplished complete thermodynamic characterization of such complexes by estimating their ∆H bind values through the independent multibox approach described in Methodology section, which in turn allowed us to calculate the entropic contributions (T∆S bind ). Table 3 summarizes the results of the previous calculations for all the analyzed systems.
As noted in Table 3, all the ∆G bind values were positive, indicating that encapsulation of the studied drugs was not spontaneous. The same behavior was evident from DFT   and semiempirical results. Furthermore, both MD and DFT calculations predicted that all the binding processes were endothermic, except for the slightly exothermic insertion of the diethylamine moiety of procaine in CB [6]. Therefore, we did not expect any interactions with a net stabilizing effect to arise during host-guest binding. On the other hand, note that, except for the CB[6]-caffeine complex in MD calculation, the formation of all the other  complexes was entropically unfavorable according to both DFT and MD calculations. Finally, having unfavorable complexes that displayed time stability during at least one of their respective long MD simulation replicates was not contradictory. In fact, the analyzed poses were manually inserted in CB [6], and the eventual exit process had to overcome a large energetic barrier associated with host ring opening. Hence, the systems remained trapped in unstable states during 1 μs MD simulations, which seems to depend on the initial velocities assigned to the atoms, as for almost all systems, except for procaine, we were able to sample dissociation events in that time scale. We expect that, in principle, longer simulation times will lead to the eventual dissociation of all complexes and no encapsulation will be observed afterwards. These simulations, however, are not required because the estimated free energies already predict that the analyzed drugs are not encapsulated by CB [6].

Discussion
The Scientific Working Group for the Analysis of Seized Drugs (SWGDRUG), 4,65 provides the most widely accepted and followed recommendation for forensic illicit drug evaluation. According to them, there are categories of techniques for forensic characterizations ( Table 4). The main recommendation is that, when a validated technique belonging to Category A is used, at least one other technique from the other categories should also be employed. Moreover, when a technique belonging to Category A cannot be used, at least three other different techniques must be applied, being two of them of Category B (and not correlated). Color tests are in Category C and they cannot be used uniquely to detect drugs, since they can provide doubtable results (false positives or false negatives). Electrochemical techniques, otherwise, are no longer considered in this recommendation.
Cucurbiturils have been used as a material able to encapsulate many substances, among them drugs of abuse. The potential of fluorescent acyclic cucurbituril (aCB) to bind opioids has been investigated. 21 Marijuana and cocaine detection with piezoelectric devices chemically modified with CB [6] has been studied. 22 CB[6]-modified electrodes have been shown to detect MDMA by voltammetry. 23 CB [7] has been applied as a sensor to detect amphetaminetype stimulants. 14 Some of these studies use cucurbituril compounds as an alternative to colorimetric tests. However, there are no studies that make a differential approach to compare drugs with their common adulterants. Therefore,  this work intended to study cocaine in a comparative fashion with lidocaine, caffeine, and procaine. Results for theoretical methods showed that the semiempirical approach did not produce reliable and reproducible results for this specific system. That is why we have decided to improve the information with DFT and MD simulations. Despite there are papers published 23,38 with semiempirical approach to experimental observation, we have demonstrated that for this specific system semi-empirical is not accurate enough to support laboratory observation.
The results concerning final conformations in DFT and MD calculations provided similar results for cocaine. Calculations converged to the same structures, showing that encapsulation preferentially occurred with the aromatic ring inserted in the cavity. Semi-empirical results were not effective since most of the positions did not reach into convergence (see SI section). For lidocaine, DFT and MD calculations also led to similar results. In this case, the amino group inserted in the CB [6] cavity was the most probable conformation.
DFT and MD calculations gave opposite results for caffeine. Convergence to CAF_P01 and CAF_P02 in the gas phase had similar energies. CAF_P03 was the only conformation possible in the presence of the solvent. For procaine we have similar results for both calculations. For DFT, the final conformation when the molecule was in the center (PROC_P01) resembled the conformation obtained for both PROC_P01 and PROC_P02 by MD.
In terms of electronic energy, the energy was similar for the cocaine final DFT conformations. The results for enthalpies showed that cocaine had similar behavior for both DFT and MD calculations. Lidocaine had positive values for DFT and MD calculations. However, the DFT enthalpic contribution was higher when compared with cocaine. For MD, the opposite trend emerged. For caffeine, the energetic results were not comparable because DFT and MD calculations provided different convergences. For DFT and MD calculations, procaine enthalpy values had opposite trends regarding dimension, but they gave essentially an endothermic behavior (slightly negative for PROC_P01 and PROC_P02 during MD).
As for free energy, both methods, as well as all the converged structures returned positive values, indicating an unfavorable binding process. DFT showed that PROC_P03 had the lowest energy values. All the molecules in the center of the CB [6] cavity displayed the highest energies. MD calculations provided the lowest value for LID_P01 and LID_P02, which gave the same final conformation when lidocaine left the center of the cavity. Except for caffeine results obtained by MD, none of the entropic contributions were favorable.

Conclusions
Investigation is only effective when scientific methods are applied to the real case. No doubt must remain about the nature of the suspicious substance. In this paper, the main objective was to investigate how CB [6] is able to differentiate cocaine and its common adulterants. It is important since there are many studies in literature that consider these substances to detect drugs as an alternative to colorimetric tests. We used in silico methods to study the energetic interaction between CB [6] and cocaine or its adulterants lidocaine, caffeine, and procaine. There was partial agreement between the MD and DFT calculations for cocaine and lidocaine. Nevertheless, the same final conformations were achieved, in which the process was found to be endothermic, not spontaneous, and entropically unfavorable. For caffeine and procaine, the opposite behavior emerged, showing that solvent influence made the difference in these cases. MD calculations provided slightly exothermic value for procaine and a favorable entropic contribution for caffeine. Calculations revealed the aromatic ring inside the CB [6] being the most probable interaction. No calculation showed that cucurbituril was specific for cocaine. Therefore, despite published works on cocaine detection with CB [6], this material should not be considered as a suitable detector for this drug. Overall, these results suggest the inability of CB [6] to encapsulate the drugs studied herein.

Supplementary Information
Supplementary data are available free of charge at https://jbcs.sbq.org.br as PDF file.