#### Abstract

It is not clear so far how humans can recognize odor. One of the theories regarding structure-odor relationship is vibrational theory, which claims that odors can be recognized by their modes of vibration. In this sense, this paper brings a novel comparison made between musky and nonmusky molecules, as to check the existence of correlation between their modes on the infrared spectra and odor. For this purpose, sixteen musky odorants were chosen, as well as seven other molecules that are structurally similar to them, but with no musk odor. All of them were submitted to solid theoretical methodology (using molecular mechanics/molecular dynamics and Neglect of Diatomic Differential Overlap Austin Model 1 methods to optimize geometries) as to achieve density functional theory spectra information, with both Gradient Corrected Functional Perdew-Wang generalized-gradient approximation (GGA/PW91) and hybrid Becke, three-parameter, Lee-Yang-Parr (B3LYP) functional. For a proper analysis over spectral data, a mathematical method was designed, generating weighted averages for theoretical frequencies and computing deviations from these averages. It was then devised that musky odorants satisfied demands of the vibrational theory, while nonmusk compounds belonging either to nitro group or to acyclic group failed to fulfill the same criteria.

#### 1. Introduction

##### 1.1. Theories regarding Odor Prediction and Recognition

Musks are responsible for bringing sensuality and warmth to a fragrance [1], as well as adding their excellent fixative properties [2, 3]. Presently, natural products extraction involves high costs and generates small amounts of raw material, so synthetic odorants are greatly demanded. To properly design a new odorant, previous information on how odors are read is necessary. There are major recent advances on neurobiology, biophysics, biochemical fields [4–7], though it is still unclear how humans recognize odor. At this scenario, quantitative structure-activity relationship models (QSAR models) are valuable for the proposal of new potential musky smelling molecules [8–12]. Such approach allows for the reduction in the number of costly and unnecessary chemical syntheses. Theories regarding odor recognition at the biological level have also been discussed in the literature. One of them, known as the odotope theory [13–15], states that noncovalent bonds (i.e., weak repulsive and attractive interactions) between the molecule and the olfactory receptor (OR) proteins (responsible for decoding the molecule and starting transduction) would elicit a unique response in the cerebral cortex.

A second theory regarding structure-odor relationship is the vibrational theory, which relates the molecule’s vibrational modes to its odor. It was firstly proposed by Dyson in 1938, followed by Wright, in 1977, and refined by several authors [16]. It claims that ORs would be able to recognize odorant vibrational modes. In 1996 [17], Turin suggested that such reading would occur, at the protein-odorant level, through a mechanism close to inelastic electron tunneling spectroscopy (IETS), an inelastic nonoptical vibrational spectroscopy [18]. After building a specific algorithm to mimic vibrational modes at protein level, the author could state that* musks*—for instance—will present four bands near 700, 1000, 1500 (or 1750 cm^{−1} for nitromusks), and 2200 cm^{−1}, and that each band will have 400 cm^{−1} width [19]. There is much to be understood concerning odor recognition, supported on vibrational theory’s ideas, as indicated by the variety of papers dealing with the subject [20–24].

##### 1.2. Paper’s Aims

In order to explore vibrational theory, musks and non-musks (see Figure 1) have been equally modeled so to come to their theoretical IR spectra. Firstly molecules have been computationally analyzed by molecular mechanics and molecular dynamics (MM/MD), in order to ascertain the conformational space available to them. After, semiempirical calculations employing Neglect of Diatomic Differential Overlap Austin Model 1 (NDDO AM1), Hamiltonian took place for geometry corrections. Then final calculations for energy optimization and IR spectra collection were set with density functional theory (DFT), using both Gradient Corrected Functional Perdew-Wang generalized-gradient approximation (GGA/PW91) and hybrid Becke, three-parameter, Lee-Yang-Parr (B3LYP) functionals, with the double numerical basis set with polarization (DNP) for GGA/PW91 and 6-311+G(d,p) for B3LYP. Methodology adopted is up-to-date with the literature and was mainly based on Lerner and coworkers paper [26]. For a complete discussion on computed harmonic vibrational frequencies and their scaling factors and accuracy, see data obtained by Scott and Radom [27].

GGA/PW91 and B3LYP theoretical spectra were compared to the experimental data (when available in free data banks as NIST and SDBS) and then were evaluated concerning the expression of the four characteristic musk bands (centered at 700, 1000, 1500 or 1750, and 2200 cm^{−1}). Now, these four bands were properly indicated using a specific algorithm built to reproduce IETS at biological level; however, the vibrational modes primarily responsible for such bands were not nominated. To clarify this question, a comparison was carried between IR theoretical data and IET data published previously, since both spectra show the same vibrational modes at the same frequencies, and it is much simpler to estimate theoretical IR spectra. Aside from the equivalence on estimating vibrational modes, there are no other comparisons that could fit in, once mode intensities are differently calculated in one technique and another, generating bands that can be either more prominent or even absent when comparing spectra. Besides, using IR spectra to predict odor or olfactory receptor affinity has failed repeatedly [28]. Finally, vibrational modes responsible for “musk bands” in the spectra are not to be taken as solely responsible for musk odor. It is well known that there are several molecular descriptors involved at major biological phenomena, such as carrying the odorant from the environment to an OR and docking odorant-OR, and they are equally relevant when considering the vibrationally assisted odor recognition mechanism.

The sixteen musks** 1–16** are already in use in perfumes, so to have their activity is certified. The seven molecules** 17–23** are structurally similar to the previous molecules but do not possess any musky odor; besides, nonmusks** 17** and** 18** are* woody*. Absence of odor is not to be taken as a different odor character but (i) as a molecule’s inability to bind its OR, due to particularities of its structure, or (ii) as the lack of specific vibrational modes in their minimum quantity, precluding proper energy decrease on electron scattering and its subsequent tunneling [13, 14, 16]. So, regarding adopted criteria relying on the choice of such molecules, it addresses the main hypothesis here being tested: only musk odorants are supposed to express vibrational modes under mentioned ranges, and no other molecule is supposed to present the same vibrational pattern. All compounds have their IUPAC name and Chemical Abstracts Service (CAS) registry number listed in Table 1.

#### 2. Methods and Methodology

##### 2.1. Computational Methods

Calculation programs used were those available in the Materials Studio Package Program [29] and Gaussian 09 [30] for B3LYP calculations. Molecules had their conformational space studied undervacuum environment through MM/DM calculations, using the CFF99/COMPASS [31] force field implemented in the Discover Simulation program [32, 33]. Main parameters chosen for MM routine were Conjugate Gradient method [34] with the Polak-Ribiere algorithm [35], running until it achieved convergence at 1.0 × 10^{−4 }kcal·mol^{−1}·Ǻ^{−1}. MD trajectories were carried out at 650 K, during 1.0 ns, with a time step of 1.0 fs (therefore, collecting 1.000 frames), in the canonical ensemble. After previous scanning in which temperature and trajectory time were changed, parameters chosen showed to be effective on expressing molecules flexibility, allowing the collection of a representative subset of initial conformations.

Selected conformers from MM/DM procedure were then minimized with Vamp program [36], using NDDO AM1 [37–39], as a previous optimization step towards DFT calculation. Spin multiplicity was automatically determined, and a convergence of 0.1 kcal·mol^{−1}·Ǻ^{−1} for the heat of formation value was set. Self-consistent field (SCF) cycles had a convergence tolerance of 5.0 × 10^{−7 }eV·atom^{−1}.

Final geometries were minimized using the DMol^{3} program [40–42] for DFT calculations. Basis set chosen was DNP [43], with GGA/PW91 [44] as exchange functional. Total energy convergence of 1.0 × 10^{−5} Ha was required, with a 3.7 Å global cutoff. Final cutoff of the SCF cycles was 5.0 × 10^{−7 }eV·atom^{−1}, and the orbital cutoff quality was kept equal to 0.1 eV·atom^{−1}. The same starting structure used in GGA/PW91 calculation was employed for hybrid functional B3LYP [45] calculation with 6-311+G(d,p) basis set, now using Gaussian 09. Molden’s processing program [46] was applied, along with Xming emulator [47], to achieve IR continuum spectra.

##### 2.2. Development of a Frequency-Weighted Average and Other Useful Tools

A mathematical formula had to be developed in this paper, in order to consider each and every mode belonging to the same 400 cm-1 range, as well as their intensities. Therefore, the use of weighted averages over theoretical data was chosen; its numerical result was called* central frequency* (cf). With this central frequency, it became possible to consider all theoretical peaks occurring in the vicinity of the diagnostic frequencies, instead of just affirming the existence of the expected band only by naming the most intense theoretical peak. Central frequencies were estimated as follows:

This mathematical approach (which simply generates a frequency-weighted average) was applied to theoretical frequencies found at ±200 cm^{−1} around 700, 1000, and 1500 (or 1750) cm^{−1} and resulted in central frequencies for compounds** 1–23**, which were computed in Table 4 (for GGA/PW91 data) and in Table 5 (for B3LYP data). The mentioned range has been chosen so as to enable each band to complete its 400 cm^{−1} width, as established in the vibrational theory.

#### 3. Results and Discussion

Data will be presented and discussed in three successive sections. In the first section the collection of theoretical IR spectra will be presented. In the second section they will be compared to experimental IR spectra to discuss their precision and accuracy. The third section will deal with the rationalization between musk spectra and musk odor activity. In order to establish such relation, best theoretical IR spectra will be chosen from this paper’s second section.

##### 3.1. Part 1: Computation of IR Spectra

Concerning* ab initio* methods, convolved GGA/PW91 spectra are shown in Figure 2, while B3LYP spectra are gathered in Figure 3. For GGA/PW91 spectra, carbon-oxygen single bond stretches occur on frequencies between 1018 cm^{−1} (for nitro analog** 21**) and 1220 cm^{−1} (for (4*S*,7*R*)-galaxolide** 7**); carbon-nitrogen single bond stretches occur under the frequency of 1200 cm^{−1} for nitro compounds; carbon-sulfur single bond stretches are observed at 1054 cm^{−1} for tonkene** 10**. Carbon-carbon double bonds stretches happened at 1687 cm^{−1} for* cis*-globanone** 4** and at 1706 cm^{−1} for* trans*-globanone** 5**, while all other compounds expressing aromatic C=C bonds had their stretching modes between 1543 cm^{−1} (as in musk ketone** 11**) and 1618 cm^{−1} (for (4*S*,7*R*)-galaxolide** 7**). Nitro groups have presented their symmetric stretching modes near 1310–1340 cm^{−1} and their nonsymmetric stretching modes at 1520–1560 cm^{−1}. For aldehydes, ketones, ethers, and esters, the carbon-oxygen double bond stretching represents an important mark for confirming spectra: here, the lowest value for such mode is found for tetralin** 9**, at 1680 cm^{−1}, while the highest value, of 1777 cm^{−1}, belongs to acyclic analog** 20**. Carbon-hydrogen bending modes start at very low frequencies and end at 1493 cm^{−1} (for helvetolide** 14**); simultaneously, the band of carbon-hydrogen stretching modes starts at 2898 cm^{−1} (for (4*S*,7*R*)-galaxolide** 7**) and finishes at 3198 cm^{−1} (for nitro analog** 22**).

For B3LYP spectra, the lowest frequency found for the carbon-oxygen single bond stretching was 1006 cm^{−1} (for romandolide** 15**), while its highest frequency was 1403 cm^{−1} (for tonkene** 10**); the frequency of the carbon-oxygen double bond stretching was found between 1738 cm^{−1} (for musk** 8**) and 1810 cm^{−1} (for musk** 14**). The carbon-carbon double bond stretching occurred on the frequency of 1707 and 1718 cm^{−1} for compounds** 4** and** 5** and between 1699 and 1712 cm^{−1} for compound** 16**; finally, the carbon-carbon double bond stretching for aromatic compounds occurred between 1531 (for non-musk** 21**) and 1650 cm^{−1} (for musk** 7**). The stretching for nitrogen-oxygen in nitro groups was found at 1250–1417 cm^{−1} for symmetric modes and at 1531–1618 cm^{−1} for asymmetric modes. Finally, carbon-hydrogen bending modes began at very low frequencies and went up to 1540 cm^{−1}, at most (for musk** 13**), while carbon-hydrogen stretching occurred between 2953 cm^{−1} (for compound** 6**) and 3229 cm^{−1} (for nitromusk** 13**).

##### 3.2. Part 2: Comparing GGA/PW91 and B3LYP Spectra to Experimental Data

For an overview over GGA/PW91, B3LYP, and experimental [47] IR data, Figure 4 and Table 2 were built, in which vibrational modes for the most important functional groups are organized in accordance with theoretical method adopted. Each band expresses the range in which labeled vibrational mode occurs, considering the maximum and the minimum found for all compounds. Thus, it becomes possible to state that B3LYP and GGA/PW91 results for C–H stretching modes start at 90 and 140 cm^{−1}, respectively, values that are higher than the ones provided by experimental data; simultaneously, the same band ends at +50 cm^{−1} for GGA/PW91 and at +90 cm^{−1} for B3LYP, in comparison to the experimental data. C–H bending modes finish at +30 cm^{−1} for GGA/PW91 and at +80 cm^{−1} for B3LYP when compared to experimental values. Aromatic C=C stretches occurred in an expected range for GGA/PW91, with a +20 cm^{−1} shift for alkene compounds; B3LYP showed a continuous band for alkene and aromatic C=C stretches, which started at +60 cm^{−1} and ended at +30 cm^{−1} over experimental data. C–O stretching modes are found within an expected range for both GGA/PW91 and B3LYP methods; for tonkene** 10**, GGA/PW91 predicted C–O mode at 1381 cm^{−1}, while the same mode is found at 1403 cm^{−1} in B3LYP data. For C=O stretching modes, GGA/PW91 showed a +20 cm^{−1} increment over lower and higher values found for experimental data, and B3LYP shifted over experimental range in +60 cm^{−1}. Finally, nitro N=O stretching modes have had both symmetric and asymmetric modes accurately and precisely estimated for GGA/PW91; for B3LYP, symmetric nitro stretching showed a broader band than that of GGA/PW91 and of experimental data, starting at a +30 cm^{−1} shift and finishing at a −10 cm^{−1} shift, while their asymmetric modes were correctly predicted.

Therefore, while GGA/PW91 final results were all fitted into experimental reports, considering the allowed ±30 cm^{−1} shifts, B3LYP results overestimated C–H stretching and bending modes, alkene C=C stretching modes, and aldehyde and ketone C=O stretching modes.

To quantify collected results, experimental IR spectrum given for each particular molecule was compared to its theoretical spectra, as shown in Table 3. In the literature—public and personal IR database—only spectra for compounds** 1**–**5**,** 8**,** 10**, and** 12** were found. The vibrational modes chosen as internal patterns of comparison were as follows, when occurring: (i) the two most intense peaks regarding the C–O stretching mode; (ii) the most intense peak for the C=O stretching mode; (iii) the most intense peaks for the symmetric and asymmetric C–H stretching modes; (iv) the highest wavenumber found for the C–H bending modes; (v) the highest intensity for the C=C stretching mode; and (vi) the most intense peaks for the symmetric and asymmetric N=O stretching modes. They are altogether able to express the* accuracy* relying on theoretical data, once all functional groups have had their modes listed within chosen criteria, and the whole spectrum has been covered.

Extra information comes on the so-called DIFF columns in Table 3. Consider where the difference established between experimental and theoretical values was plotted (for both GGA/PW91 and B3LYP results), as shown in (2). In such equation, for a certain mode , the difference (DIFF) is obtained through the subtraction of its theoretical value (TheoWavenumber) from its experimental value (ExpWavenumber). The average of the collected differences on Table 3 was estimated as the values of absolute difference were divided by the number of occurring modes, disregarding all C–H stretch modes, as follows:

This choice was made due to the fact that all C–H stretches were overestimated in more than 100 cm^{−1} for both theoretical methods, as discussed before. This last mode, therefore, had its deviation separately estimated as follows:
in which the absolute is collected, is summed for all C–H stretching modes within the same theoretical method, and is divided by the number of occurring modes.

So finally GGA/PW91 IR results are found to be more accurate than B3LYP results, with calculated deviations that are equal to 13 cm^{−1} and 43 cm^{−1}, respectively. Nevertheless, they correspond to a 3.25% (DFT) and a 10.75% (B3LYP) shift over the 4000–0 cm^{−1} that represents the spectrum, which ensures that both methods are powerfully predictive. Simultaneously, both methods display great shift over C–H stretching modes, corresponding to an increase of 124 cm^{−1} and 145 cm^{−1}, respectively. Considering that one of this paper’s aims is to compare theoretical IR data with postulates of the vibrational theory on musky odor, the last mentioned modes are irrelevant, as they are not supposed to be responsible for muskiness.

##### 3.3. Part 3: IR Spectra and Musk Activity

Herein, both GGA/PW91 and B3LYP data will be taken into account so as to compare theoretical results with the bands described in the vibrational theory (700, 1000, 1500 or 1750, and 2200 cm^{−1}), which would assign, as supposed by the previous theory, a musky odor to these molecules.

In Tables 5 and 6 the information was gathered so as to dispose musks** 1** to** 23** in lines and their central frequencies (at 700, 1000 and 1500 or 1750 cm^{−1}) in columns. Besides, the compounds were gathered in groups, in accordance with their structural similarity (see first column on the left): MCM for macrocyclic musks, PCM for polycyclic musks, NM for nitromusks, and ACM for acyclic musks. Nonmusks are called MCa for macrocyclic analogs, Aa for acyclic analogs, and Na for nitro analogs. As no vibrational mode has occurred at 2200 ± 200 cm^{−1}, there is no such column.

For each musk group, the central frequencies averages (*ave*, (5)) and their standard deviations (*std*, (6)) were estimated: these values are plotted at the right side of the central frequency list. The next column on the right shows the differences between each central frequencies and the average of central frequencies (*diff*, (7)). Consider
where represents the central frequencies found for each group (MCM, PCM, etc.), having herein been defined lower and upper limits for , so that for MCM; for PCM; for NM; and for ACM. Consider

The values found for the and of the MCM group were the same as those found for the MCa group; those belonging to NM were repeated for Na; and data calculated for ACM were taken as parameters for Aa.

Finally, each value was allowed to be lower or two times its std_{group} value, as stated in

If this happens, it is possible to say that theoretical data fit vibrational theory’s demands over study band; if not, they do not agree. When repeating this process for all three bands under the spectra, values for the Total Deviations (TD) are collected. The zero input (0) means that all three bands obeyed (8). Logically, TD inputs between 1 and 3 indicate how many of the bands presented a value that was higher than the doubled value.

##### 3.4. Absence of the 2000–2400 cm^{−1} Band

None of the previous tables bring a mode among the 2000–2400 cm^{−1} range. The 2200 cm^{−1} band was described by Turin as the carbonyl stretching mode, which only occurred at this range when working with NDDO/AM1 Hamiltonian. Through* ab initio* calculation, the mentioned mode was observed at the range of 1680–1810 cm^{−1}, and it was only taken into account for cfs estimates when dealing with the 1550–1950 cm^{−1} band for nitro compounds. Regarding structural features, compounds** 6**,** 7**,** 12**, and** 13** do not possess the carbonyl group and are real musks, while compounds** 17**,** 18**,** 19**,** 20**, and** 23** contain the carbonyl group and are odorless. Therefore, carbonyl stretches cannot be taken as criteria for IR data-musky behavior profile.

##### 3.5. GGA/PW91 cf’s Plots and Musk Activity

Table 5 (with GGA/PW91 results) presents only two molecules, nonmusks** 20** and** 23**, which failed to express a complete match between theoretical data and the signed bands of the vibrational theory. For nonmusk** 20**, the cf value estimated near 1000 cm^{−1} was higher than its ave_{group} value: this compound had both O– stretching modes occurring at high intensities and frequencies below 1200 cm^{−1} (which are precisely the 1090 and 1144 cm^{−1} wavenumbers). On the other hand, musks** 14–16** and nonmusk** 19** have presented a stretch mode (or at least one of them) outside the 1000 ± 200 cm^{−1} range. Concerning nonmusk** 23**, it can be said that it showed a lower cf value, which remained around the 700 cm^{−1} band, and a higher cf value, near the 1000 cm^{−1} band. This can be explained by the fact that, (i) referring to the 700 cm^{−1} band,** 23** had an average for intensity modes that was lower than the average values for the other Na and NM compounds, and, therefore, lower* weight values* (see (4)) and lower cf; and (ii) referring to the 1000 cm^{−1} band, the O– stretch mode at 1073 cm^{−1} had an intensity value that was 8.7 times higher than the second most intense mode, generating the observed overestimated cf value.

That being said and according to mathematical approach adopted in this paper, no relationship was found between IR properties and the presence of musky odor through GGA/PW91. It does not mean that GGA/PW91 calculation leads to affirm that musks and nonmusks have indistinguishable spectra. It means that, considering the weighted-averages values resource in order to simplify the adoption of all theoretical peaks under analysis, resulting central frequencies are similar among musks and nonmusks. Such observation may be due the fact that GGA/PW91 spectra had narrower bands for each vibrational mode rather than B3LYP spectra. As GGA/PW91 proved to be a more accurate method, it requires a more accurate statistical tools, so to properly deal with all theoretical modes, as to affirm (or deny) a possible difference on musks and nonmusks IR spectra.

##### 3.6. B3LYP cf’s Plots and Musk Activity

On the other hand, Table 6 (with B3LYP data) shows that all musks from** 1** to** 16** have had theoretical central frequencies that fitted the frequencies demanded at 700, 1000, and 1500 (or 1750) cm^{−1} in the vibrational theory. Simultaneously, Aa and Na showed that at least one of their cfs was too shifted from their averages. Molecule** 23**, for instance, is the only nonmusk possessing an aldehyde functional group, and its C=O stretching mode is found at the 1550–1950 cm^{−1} range under study. Nevertheless, nitro musk** 11** also holds an aldehyde functional group, and though its cf at 1750 cm^{−1} is indeed higher than the other cfs for the NM, it is still 70 cm^{−1} lower than** 23**’s cf. Still concerning the analysis of Table 6, we can say that the acyclic analogs** 19** and** 20** have had their central frequencies shifted to lower values at the 1500 cm^{−1} diagnostic band, the only one in which there was no agreement between theoretical results and expectations deriving from the vibrational theory. At the 1500 ± 200 cm^{−1} range, hydrocarbon scissoring modes and alkene/aromatic C=C stretching modes are found, and** 19** and** 20** lack the gem-methyl at the cyclohexane that are observed in** 14** and** 15**; moreover, although in** 16** this gem-methyl is also absent, it possesses two alkene groups, for which its stretching modes are also found at this range. Finally, macrocyclic analogs showed the same spectra data as macrocyclic musks: as compounds** 18** and** 19** differ from MCM in that they lack two methylene groups, it is possible to suggest that such loss leads the molecule’s C–H vibrational mode sum to be too low to be detected by its OR, as suggested elsewhere [23]. Besides, both** 18** and** 19** present a* woody* smell; by being smaller, they may be able to bind to a variety of different ORs, which do not reassemble the behavior of the musk, detected by possibly only one OR or just a few more [48].

#### 4. Conclusions

This paper worked on developing a method which could state the following: (i) whether musky molecules had their IR theoretical spectra fitted into the four bands assigned by vibrational theory in 2002 Turin’s paper and (ii) whether nonmusky molecules did not have their IR theoretical spectra fitted into same criteria. It was first concluded that only three bands were useful for discriminating musk from nonmusk odorants. The fourth IR band refers to a carbonyl stretching mode, and data here presented showed that this mode cannot be taken as criteria when building a muskiness profile.

Present work was carried out with two different functionals applying DFT calculation level, and both were able to accurately predict IR spectra. A mathematical strategy has been adopted in order to summarize IR theoretical data and properly compare musk to nonmusk molecules. Using this approach, GGA/PW91 computed spectra showed no discriminating patterns between musks and nonmusks. B3LYP data showed theoretical frequencies and the three diagnostic bands to be matching, while nitro and macrocyclic nonmusks failed in at least one of such comparisons.

It was therefore possible to conclude that no specific mode alone is responsible for the discrimination of the musky odor. Instead, when a great amount of frequencies, expressing weak or medium intensities, were found along the 700–1700 cm^{−1} range, musk odor occurs. At this range, observed modes were CH_{2} wag and scissor bending modes, C=C stretches, and symmetric nitro stretching modes.

This trial scanning over IR vibrational modes and their relation to musk odor perception brings new informational towards the understanding and application of vibrational theory’s ideas and opens the possibility to a more extensive study, in which a broader number of molecules may be modeled so as to come to a quantitative proposal for a musky-like novel odorant, regarding their IR vibrational modes.

#### Conflict of Interests

The authors declare that there is no conflict of interests regarding the publication of this paper.

#### Acknowledgment

The authors thank Professor Inês Sabioni Resck (from Instituto de Química, Brasília, Brazil), for initial ideas and reliance.