Accurate prediction of terahertz spectra of molecular crystals of fentanyl and its analogs

Fentanyl is a potent synthetic opioid pain reliever with a high bioavailability that can be used as prescription anesthetic. Rapid identification via non-contact methods of both known and emerging opioid substances in the fentanyl family help identify the substances and enable rapid medical attention. We apply PBEh-3c method to identify vibrational normal modes from 0.01 to 3 THz in solid fentanyl and its selected analogs. The molecular structure of each fentanyl analog and unique arrangement of H-bonds and dispersion interactions significantly change crystal packing and is subsequently reflected in the THz spectrum. Further, the study of THz spectra of a series of stereoisomers shows that small changes in molecular structure results in distinct crystal packing and significantly alters THz spectra as well. We discuss spectral features of synthetic opioids with higher potency than conventional fentanyl such as ohmefentanyl and sufentanil and discover the pattern of THz spectra of fentanyl analogs.

To facilitate THz identification of fentanyl, it is advantageous to assemble a database of fentanyl THz spectra. This can be done, for instance, by quantum mechanical (QM) density-functional theory (DFT) calculations in solid state. Such efforts will enable chemometrics techniques for understanding and validation of reflectance spectroscopy experiment, and the fast prediction techniques by machine learning (ML) methods. In this work, we will establish a library of THz spectra of fentanyl and its analogs via QM calculation.
The detection of fentanyl by THz spectra is possible due to unique phonon properties (molecular motion). Each phonon mode in the given wavelength range corresponds to a hindered molecular rotation (libration), coupled to intramolecular vibrations. Differences in the molecular structure of the fentanyl analogs leads to distinctly different crystalline packing and hence THz spectra. In general, long-distance intermolecular interaction determines the strength and the position of peaks on THz spectra in the solid state. Specific intermolecular interactions such as hydrogen bonding and π-π stacking tend to have significant influence of THz frequency shift 45 .
In order to accurately describe long-distance dispersion interactions for organic molecular solids, we have decided to apply Grimme's PBEh-3c method 46 , based on corrected Perdew-Burke-Ernzerhoff (PBE) 47 hybrid generalized-gradient-approximation functional. Conformational analyses of fentanyl and its analogs ohmefentanyl 48-50 , 3-methylfentanyl 49 , BIYTAF 48 , sufentanil 48 , and alfentanil 48,51 have been studied previously in different groups by molecular dynamics and molecular docking to understand the docking site of fentanyl analogs to μ-opioid receptor and by DFT calculations to understand stable conformations in gas and aqueous phase. To the best of our knowledge, this is the first time THz spectra of fentanyl and its analogs in solid state had been predicted by DFT calculation. Other factors which may affect measured THz spectra (such as mixture of different fentanyl analogs and hydrates) can be aggressed in our future work by ML methods.

Computational methods
All the crystal structures of the various fentanyl molecules were obtained from the Cambridge Structural Database 52 . The structures were preprocessed using cif2cell 53 to convert crystallographic information file to the proper input file format. Some of the fentanyl analogs were missing hydrogen atoms, which were added geometrically using Molden program 54 . Marvin 55 was used for drawing, displaying, and characterizing two-dimensional chemical structures and Mercury 56 was used for visualizing three-dimensional geometry. Post-processing of the output files to produce the graphical representation of the absorption spectra and vibrational normal modes was done using CRYSPLOT 57 . Microsoft Excel 2019 was used for presenting absorption spectra 58 . CRYSTAL17 59 , was utilized for solid-state DFT calculations with periodic boundary conditions. Sufficiently large supercell (between 53 and 134 atoms) was constructed to ensure long range vibrations were accounted for. The multiplicity of the unit cell was determined by the ratio of the volume of the simulation supercell to the volume of a primitive cell. Corrections for basis-set superposition error and long-range dispersion interactions were included within PBEh-3c method 46 . This method is based upon the modified global hybrid functional PBEh, where percentage of Hartree-Fock exchange (42%) is adjusted to get correct average bond length. Modified double-ζ Gaussian basis set dubbed def2-mSVP 46,60 is defined within the method as well. The benchmark study, validating the accuracy of the chosen computational method is described in detail in Supporting Information.
Full geometry optimization, including both atom coordinates and unit cell parameters (a, b, c, α, β, and γ), was applied to fully relax the geometry within the constraints of the space group symmetry. Threshold of convergence for crystal structure, based on the energy change between consecutive geometry optimization steps was set to 10 −8 Bohr. Vibrational normal modes 61 and IR intensities were calculated using the Berry phase method 62 . For the calculated vibrational frequencies, the scaling factor of 0.95 was applied for all of the spectral predictions 46 . The damping factor for the calculation of the Lorentzian linewidth at half maximum of absorption peak was set to 2.0 cm −1 , in order to match the spectral shapes recorded at room temperature. Threshold on Self Consistent Field energy convergence was 10 −10 Hartree. Truncation criteria for bielectronic integrals were as follows: overlap threshold for Coulomb integrals was 10 −7 , penetration threshold for Coulomb integrals was 10 −7 , overlap threshold for HF exchange integrals was 10 −7 , pseudo-overlap (HF exchange series) were 10 −7 and 10 −14 . These criteria are based on well established guidelines which were described in Ref 63 .

Results and discussion
Geometry optimization, hydrogen bond, and charge distribution. Predicted unit cell parameters of fentanyl analogs are reported in Table S4. Their deviations from experimental values (collected in Table S3  Hirshfeld population analysis 64 describes molecular charge density in terms of atomic contribution and allows one to examine the nature of the charge transfer between various species in each crystal structure. The magnitude of Hirshfeld atomic charges is less dependent by basis set and is typically smaller than that of Mulliken charges 65 . Of all the fentanyl and its analogs, atomic charge of O atom in N-phenyl-propanamide is − 0.555 ± 0.016, N atom in N-phenyl-propanamide is − 0.363 ± 0.047, and N atom in piperidine ring is − 0.189 ± 0.027 (neutral) or 0.031 ± 0.059 (if protonated, H atom is 0.248 ± 0.067). Hence, the charge distribution in fentanyl and its analogs is consistent.
THz spectra of fentanyl and its analogs. Predicted THz absorption spectra of fentanyl and selected analogs are shown in Fig. 2. These correspond to polycrystalline samples (anisotropy is averaged out). Even one functional group difference results in a change of the crystal packing and thus the shape of absorption spectra. As one can see in Fig. 2, an additional methyl group for 3-methylfentanyl, a hydroxyl group for BIYTAF, and a carboxyl group for TACDOS results in different THz spectra. There are various medium peaks between 35 and www.nature.com/scientificreports/ 85 cm −1 of 3-methylfentanyl (Fig. 2a). Fentanyl, 3-methylfentanyl, and BIYTAF all have strongest peak above 85 cm −1 . Normal mode visualization indicates that two peaks for fentanyl at 45 and 58 cm −1 arise from the intermolecular wiggling of T-shaped packing of the phenyl rings, which is not found in the other 3 analogs (Fig. 2a,b). The addition of methoxymethyl group on piperidine ring results in the appearance of several weak peaks in THz spectra of R-30490 ( Fig. 2c) but no significant peaks like 45, 58, and 68 cm −1 in fentanyl and 68, 77, and 86 cm −1 in TACDOS (see Fig. 2a,b). Strong peaks arise from N…H H-bonds to piperidine ring in fentanyl and TACDOS. Comparison between TACDOS and TACDUY (with solvents water and methanol) is shown in Fig. S2. Medium to strong peaks between 25 and 75 cm −1 are from H-bonds between solvents water and methanol. For alfentanil and thiofentanil (Fig. 2d), the strong peaks at 56 and 62 cm −1 arise from Cl…H stretch between counterion and protonated piperidine ring. R4 group in both R-30490 and sufentanil (see Table 1) does not form strong H-bonds, and dispersion interactions produce small but sharp peaks below 80 cm −1 . Figure 3a shows the vibrational normal modes of fentanyl at 45, 58, and 67 cm −1 and of 3-methylfentanyl at 48, 57, and 62 cm −1 are visualized in Fig. 3b. The difference in one methyl group leads to a distinct crystal packing and normal modes of fentanyl and 3-methylfentanyl. Similarly, sufentanil where the phenyl ring is substituted by thiophene ring results in weak peaks near 57, 67, and 78 cm −1 , while alfentanil with the substitution to tetrazole ring leads to strong peaks at 56, 85, and 88 cm −1 .
THz spectra of fentanyl analogs with chiral centers. The mode of molecular packing influences THz spectra significantly. Even for fentanyl analogs of the same atomic composition, stereoisomers induce changes in the packing and location of librational centers. Figure 4 shows the example of different crystal packing formed by the 3 ohmefentanyl isomers from R2, R3, and R4 groups (see Table 1) along with b-or a-axis. Herringbone packing is visible when viewing cis-(2R,3R,4S) and cis-(2S,3R,4R) from 3 different directions in Fig. 4a,b,d. Orange arrow in Fig. 4b shows T-shaped interaction. However, parallel packing is seen when viewing trans-(2S,3R,4R) along a-and b-axes while herringbone packing is seen along c-axis (Fig. 4c). The orientation of OH in R2 and CH 3 in R3 groups influences the packing significantly.
THz absorption spectra of stereoisomers for selected fentanyl analogs are shown in Fig. 5. In ohmefentanyl (Fig. 5a) with stereoisomeric variations cis-(2R,3R,4S) and cis-(2S,3R,4R), similar THz spectra shows both the crystal packing include T-shaped interaction (Fig. 4a,d). There is merely a slight wavelength shift and it represents the packing dominates THz spectra. Intermolecular H-bond between ohmefentanyl molecules is absent even in the presence of OH group in R2 (see Table 1). Conversely, ohmefentanyl trans-(2S,3R,4S) exhibits a different spectrum with the peak at 79 cm −1 because of π-π stacking interaction from phenyl ring in R1, which is absent in cis- (2R,3R,4S) and cis-(2S,3R,4R).   Fig. 5b. Both of them have several small peaks with one strong peak at 53 or 78 cm −1 . Because of the solvent molecules and/or counterion in FOPFIZ and CEWDEN10, small and medium peaks are present to reflect the formation of multiple H-bonds. The strongest peak of FOPFIZ at 53 cm −1 shows the T-shaped phenyl ring interactions (see Fig. 6a). The strongest peak of CEWDEN10 at 78 cm −1 shows interactions between the dangling ethyl groups instead of phenyl rings (CH 3 …π interaction, Fig. 6b).
JINPAX (2R,5S) and VEYCIL (2S,5S) in Fig. 5c have similar situation that combine several small peaks with one or two strong peaks at 95 or 98 cm −1 (shown in Fig. 6c,d). THz spectra therefore distinguishes stereoisomers in these examples. JINPAX contains hydrogen sulfate that can form H-bond with nitrogen atom in piperidine ring. CH…π and parallel/T-shaped phenyl ring interactions determines the appearance of medium and strong peaks on both JINPAX and VEYCIL spectra.
All the other THz spectra of fentanyl analogs FBPIPA, BEGZUJ, IFIGIN, ZIZPON, and XALTAF can be found in Fig. S3. Strongest peaks above 75 cm −1 are based upon parallel π-π stacking, and T-shaped CH…π, CH 3 …π interactions. Note that there is no significant H-bond interaction in these analogs. Medium and small peaks below 75 cm −1 arise from long-distance, displaced stacking π-π and T-shaped interactions or from solvent molecules and counterions.

Conclusions
Spectra of fentanyl and its analogs in THz range have specific characteristics. First, the generic picture of THz spectra of fentanyl and its analogs is that between 75-100 cm −1 , there are 2-4 strong peaks from significant π-π interactions and/or H-bonds. Second, THz range between 20-75 cm −1 consists of small to medium peaks because of non-ideal π-π interactions, H-bond, and dispersion interaction. Third, the existence of π-π interactions or H-bonds that belong to solvent molecules or counterions leads to medium or strong peaks between 20-75 cm −1 .
Compared with compounds that are in white powder form (see Fig. S1), THz spectra of fentanyl and its analogs have the preference of peaks at certain THz range. The three characteristics may help to apply ML to predict unknown street drugs in the future and enable the detector design for synthetic opioids.  License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.