Melatonin/nanoclay hybrids for skin delivery

Melatonin is a hormone used for treating several disorders. However, its oral administration is problematic due to the variable absorption and extensive first-pass metabolism. The topical application of this drug does not present these disadvantages and can be used for therapeutic and skin protection purposes. The adsorption and release of melatonin in the montmorillonite clay was investigated experimentally and theoretically. The drug-clay interaction products were prepared and characterized, showing a modified and improved diffusion and release of melatonin. The crystal structure and spectroscopic properties of melatonin and melatonin-water-mineral materials were studied by molecular modeling, finding that the adsorption of this drug is energetically favourable, and the results were consistent with the experimental data. Ab initio molecular dynamics simulations with the metadynamics methodology showed the release of melatonin from the confined interlayer nanospace of montmorillonite to the liquid bulk water with a 20 kcal/mol free energy bar-


Introduction
Melatonin, N-acetyl-5-methoxytryptamine, is an indolamine derived from L-tryptophan. It is produced mainly by the pineal gland in brain [1], in a five-stage synthesis starting from tryptophan and with serotonin as an intermediate.
At present, melatonin is the focus of intense research because of their many physiological functions and pharmacological effects [2]. In addition, melatonin can have direct or indirect antioxidant action [3][4][5][6], neutralizing several toxic radicals or stimulating the synthesis of other antioxidants and preserving the enzymatic functionality [7]. Melatonin also regulates the electron transport chain in mitochondria, protecting mitochondrial DNA and proteins against oxidative damage [8,9], and also protecting the cells against ionizing radiation [10,11]. Several studies have demonstrated the benefits of its topical use for therapeutic purposes and skin protection [12][13][14]. In addition to providing skin protection, melatonin can be also a useful product for wound healing treatments [15][16][17][18][19]. Melatonin increase angiogenesis and stimulates interleukin-1, TNF-α cytokines and TGF-β1 [20,21]. It has been recently included in nanofiber wound dressings [22]. However, melatonin shows low water solubility and light degradation [23][24][25].
Nanoclays are good candidates to be used as vehicles of melatonin because of their widely known capabilities for the retention of organic molecules, and protecting the loaded actives from degradation and improving their release profiles once administered [26][27][28][29]. Nanoclays are well known ingredients of skin formulations [30] once they have accomplished safety and stability requirements [31]. The application of clays in skin drug delivery has also been studied [32,33]. Clays are ideal excipients for the development of intelligent skin drug delivery systems to be used in wound healing [32,34]. Zinc oxide nanoparticles confined inside halloysite nanotubules revealed their ability to interact with a wide portion of UV-vis radiation, which leads to enhanced sun screening performance [35]. In vitro biocompatibility (cytotoxicity and proliferation), antimicrobial efficacy and cell motility gap closure (assay of wound closure) of montmorillonite/chitosan nanocomposites loaded with silver sulfadiazine were studied [36]. A nanocomposite based on nanotubular clay mineral halloysite and chitosan oligosaccharides has also been proposed to enhance healing in the treatment of chronic wounds [37]. Polymer films loaded with a carvacrol clay hybrid for the modified delivery of carvacrol in infected skin ulcer treatment was designed by [38]. Electrospun scaffolds based on biopolymers and loaded with montmorillonite or halloysite have shown enhanced fibroblast cells attachment and proliferation with negligible proinflammatory activity [39]. Polysaccharide-based scaffolds (chitosan, chitosan/chondroitin sulfate and chitosan/hyaluronic acid loaded with norfloxacin (as a free drug or in montmorillonite hybrids) have been also developed [40,41]. The drug intercalated into montmorillonite interlayer and the hybrid was even more effective as antibacterial that the pristine drug. Montmorillonite, laminar mineral, can be used to modify the release of drugs [29,42], and to protect them from the degradation [43]. Moreover, it have been recently demonstrated that clays have a proliferative effect on fibroblasts [44], which can lead to a potential synergistic effect between melatonin and montmorillonite for wound heling and skin protection.
In order to advance the current technology of melatonin release systems for skin delivery treatments that need release control, it is also important to characterize the physical adsorption on the confined nanospaces of clay minerals. Important physical and chemical properties are present here, which study can shed light on the characteristics of clay minerals confined nanospaces and the gradually drug release from the melatonin-clay system. Therefore, through the long path that melatonin has to follow before going out of the clay mineral confined space, it will have different geometrical and energetic interactions with the surfaces, interlayer cations and water, which are necessary to unveil to understand the adsorption-diffusion-desorption phenomena occurring in these confined systems. This knowledge will reveal the structure and mechanisms, which will be useful for designing new delivery systems. To do so, we performed an experimental and computational investigation on the adsorption and release of melatonin by clay minerals, which has allowed rationalizing the excellent performance of clay minerals as modified and improved release media for melatonin.
With these premises, in this study the design of melatonin-nanoclays hybrids has been accomplished, including sorption and release experiments, solid state characterization of the systems and molecular modeling. Preparation of MEL-clay minerals interaction products Clay powder, VHS, was dispersed in 100 mL of acetonitrile solution of MEL under magnetic stirring at room temperature for 24 h, in order to ensure a complete interaction between the components, in such a way that the drug-clay ratios were 1:1 (MEL1-VHS1) and 1:5 (MEL1-VHS5) w/w. After 24 h, the solvent was evaporated at 40°C by means of a rotary evaporator (Buchi® R II) and the solid residue (MEL-VHS) was stored in a desiccator at room temperature with freshly activated silica-gel grafted with a moisture indicator of Co salt.

X-ray diffraction
Powder X-ray diffraction was performed with an X-Pert Pro® diffractometer (Marvel Panalytical, Madrid, Spain) and the CuKα radiation. The powdered samples were scanned in the range of 4º-50º of the 2θ angle, steps were of 0.008 of 2θ units and the counting time was of 10.16 sec/step. The information of diffraction was analyzed using X-Pert Data Viewer® software.

Thermal analysis
Differential Scanning Calorimetric analysis (DSC) and thermogravimetric analysis (TGA) were performed with a Mettler Toledo mod. TGA/DSC1 calorimeter, equipped with sensor and FRS5 microbalance (precision 0.1 µg) and FP89 software package. The heating rate was 5°C/min in the 30-950°C temperature range for TGA; and for DSC, the heating rate was 2°C/min in the 20-270°C range for additional specific DSC runs. Air was used in TGA and nitrogen was used as purge gas in DSC under flow of 15 mL/min. 2.1.2.5. Fourier transform infrared (FTIR) spectroscopy FTIR spectra were recorded in the range 4000-600 cm −1 with a 0.5 cm −1 resolution and a well-plate sampler (JASCO 6200 spectrophotometer with Spectra Manager II software).
2.1.2.6. In vitro release studies: Franz diffusion cell Experiments were performed by means of Franz diffusion cells by using 15 mg of MEL1-VHS5 interaction product, containing 2.5 mg of pure MEL drug. The in vitro release studies were applied in the MEL1-VHS5 sample, because the drug-clay ratio 1:5 w/w is the one that presents the greatest interaction between MEL and mineral. In addition, 2.5 mg of pure MEL was also tested as a reference. The powder of MEL and MEL-VHS were subjected to in vitro release studies with a Franz diffusion cell with 8-mm orifice diameter (diffusion area 0.5 cm 2 ), thermostated at 32±0.5°C and with isotonic saline solution at pH 7.4 (NaCl 9 g/L, as a dissolution medium simulating the plasma). Membrane filters of 0.45 µm (Millipore®) were positioned between the donor and the receptor compartments of the Franz cell. The donor compartment was filled with 400 µL of isotonic saline solution, in order to reduce dehydration of the membrane and the sample. The receptor compartment was filled with 5.8 mL of isotonic saline solution and continuously stirred with a magnetic bar. The isotonic saline solution was previously degassed under vacuum and ultrasound conditions. In addition, the Franz cells were connected to a peristaltic pump that constantly filled the cell with a flow of 1.54 mL/min with isotonic saline solution, ensuring compliance with sink conditions. At established time intervals until 10 hours, samples were withdrawn from dissolution medium, by filtering through 0.45 µm Milli-pore® membranes. The amount of drug dissolved in the samples was measured by HPLC, by using the equipment and methodology explained in the next section. The drug assays were repeated in 8 replicas to ensure that the results are reproducible.
2.1.2.7. HPLC Analysis Quantitative determination of MEL was performed by using a 1260 Infinity II Agilent HPLC system (Santa Clara, CA, USA) equipped with quaternary pump, autosampler, column oven and UV-VIS diode-array spectrophotometer. The stationary phase was a Kro-masil® C18 column, 5 µm, 250 × 4.6 mm (Teknokroma, Barcelona, Spain) and the mobile phase was a mixture of H2O and CH3CN (35:65 v/v). The flow rate was set at 0.8 mL/min with a 10 µL injection. We used a Perkin-Elmer Lambda-25 spectrophotometer detector at a 225 nm wavelength, with 5 min runtime for each analysis. Data were recorded and analyzed by using software LC Open LAB HPLC 1260 (Agilent). The analytical method was linear in the concentration range 0.1-15 mg/L of MEL in isotonic saline solution, resulting in correlation coefficients of 1.

Computational section 2.2.1. Atomistic Simulation Models
The crystal structure of MEL was taken from experimental crystallographic X-ray diffraction data analysis (MELATN02, ref. 1210960 of CSD database) [47]. One molecule of MEL was extracted from this crystal for studying the molecular structure with force fields (FF) and density functional theory (DFT) methods. In the adsorption calculations, this molecule was studied in a periodical box with the same dimension of the VHS crystal. Taking into account the chemical composition of the mineral used in experiments, the unit cell of our VHS model was Na(Al3.17Mg0.83)(Si7.83Al0.17)O20(OH)4 [42]. For the adsorption calculations a 3x2x1 supercell of VHS model was created with formula Na6(Al19Mg5)(Si47Al1)O120(OH)24. This supercell was designed by considering the molecular dimension of MEL for avoiding close contacts between the MEL molecules in neighbouring cells of the periodical boundary conditions model. The octahedral Mg substitutions were placed maintaining a maximal disperse distribution along the octahedral sheets, as it was considered in previous works [48,49]. The adsorption complex was created placing one, two or three MEL molecules in the interlayer space of 3x2x1 supercell montmorillonite and, besides, twelve molecules of water per supercell were included around MEL in the interlayer space.
Ab initio molecular dynamics (AIMD) simulations were performed with a 3x2x1 supercell of montmorillonite without tetrahedral substitutions with formula Na6(Al18Mg6)(Si48)O120(OH)24. We created a nanoparticle model of montmorillonite cleavaging (100) surfaces based, similar to models used on previous works [50]. This model was placed in a simulation orthorhombic supercell with dimensions a = 36.0 Å, b = 18.13 Å, c = 15.0 Å, vacuum space was filled with 206 water molecules along the x and y axes and the interlayer space along z-axis and placing the MEL in the interlayer space ( Figure  S1). In this way the resulting model is a MEL···VHS supercell in contact with a bulk water reservoir. All the dangling bonds as consequence of the cleavage were saturated with hydrogen atoms or OH/OH2 groups as shown in Figure S1. During all the AIMD simulations, the row of octahedral ions with the lowest x value was kept frozen, thus effectively 'pinning' the model to a certain (x,z) point but otherwise allowing it to 'tilt' around the axis defined by that point.

Computational details
Atomistic calculations based on empirical interatomic potentials were performed for the optimization of the molecular and crystal structure of MEL and the study of its infrared spectra. Several FF have been compared, such as Compass (CF) [51], Universal (UF) 5 of 23 [52], and INTERFACE [53,54]. The last one has provided good results in previous studies about the interactions of organic molecules with phyllosilicate surfaces [55,56]. To carry out these calculations, we used the Discover and Forcite programs within the Materials Studio package [57].
Ewald summation and atom methods were used to compute the Coulomb and Van der Waals interactions, respectively with a cut-off of 18.5 Å, which was validated in previous works [42,[58][59][60].
Simultaneously, we used the SPC model for water molecules (atomic charges values of +0.41 and -0.82 for H and O atoms, respectively) and also the TIP3P model (charges values +0.417 and -0.834 for H and O atoms, respectively) [61,62] in the calculations performed with INTERFACE FF.
Besides, quantum mechanical methods based on DFT were also used to study the isolated and crystalline structure of MEL with the DMol 3 [63] and CASTEP [64] codes, including the dispersion corrections of Tkatchenko-Scheffler [65] and Grimme [66]. DMol 3 is based on the linear combination of numerical atomic orbitals and we used Double-ζ polarized (DZP) basis sets, whereas CASTEP is based on plane-waves and we used a mesh cut-off energy value of 500 Ry. In both cases, we used the generalized gradient approximation (GGA) with the revised Perdew-Burke-Ernzerhof exchange-correlation functional (RPBE) [67], using only the Γ point of the Brillouin zone. Semi-core pseudopotentials (DSPP) were used in DMol 3 and ultrasoft pseudopotentials were used in CASTEP. The convergence threshold for the self-consistent field (SCF) was 10 −6 Ha. The normal vibration modes were studied within the harmonic approximation by applying steps of finite displacements on every atom. The frequencies of these vibrational modes were obtained by diagonalizing the dynamical matrix related with the geometrical changes generated. This procedure allows to localize the minima structure in the potential energy surface and to obtain the normal vibration mode frequencies to compare with experimental data. This procedure was validated in molecular aggregates [68] and crystalline systems recently [69,70].
The AIMD simulations were carried out with the module QUICKSTEP [71] as implemented in CP2K [72]. The electronic structure of the system was treated within the Restricted Kohn-Sham formalism, employing RPBE exchange-correlation functional [67] including Grimme dispersion corrections (DFT-D2) [66]. We used Goedecker-Teter-Hutter (GTH) pseudopotentials together with molecularly optimized short-range DZP basis sets (DZVP-MOLOPT-SR), and a 600 Ry density cut-off for the expansion of the auxiliary plane waves basis set. The AIMD was propagated employing the second-generation Car-Parrinello MD scheme [73], setting the parameters for the Langevin-like equation γL=0.002 fs -1 and γD= 4.5×10 -8 fs -1 with a temperature of 300 K and a 0.4 fs timestep. The system was equilibrated for more than 10 ps, during which we adjusted on-the-fly the value of γD. In order to observe the diffusion of the MEL from the interlayer space to the bulk water reservoir out of the mineral in a computationally affordable time, we exploit Metadynamics [74], which was used satisfactorily in these systems previously [75,76]. In this work, we used the Multiple-Walker implementation for a more efficient sampling of the free energy surface (FES) [77]. The employed collective variables (CVs) were the X and Z coordinates of the center of mass of the MEL, which effectively are equivalent to its coordinates with respect to the montmorillonite slab since the latter is pinned to a fixed (X, Z) point. The height of the Gaussians comprising the history-dependent potential was 0.0015 au and their widths were varied between 0.14 and 0.06 Å depending on the spacing of the walkers, in order to reach an efficient reconstruction of the FES. We started the simulation with 13 walkers, adding more walkers as they became increasingly spaced in the CVs space. Eventually, we used 49 walkers in total.

Preparation and solid characterization of the melatonin-clay hybrid complexes
After testing the stability of melatonin (MEL) in different solvents (see Supplementary Materials), we determined, via HPLC measurements, that MEL immediately degrades in ethanol, methanol, isopropanol and acetone. However, MEL remains stable in acetonitrile even after one month, with a solubility of approximately 41 mg/mL. Thus, the interaction products were prepared by solving MEL in acetonitrile and subsequently dispersing in montmorillonite (VHS), following the methodology described above, in drugclay ratios 1:1 (called MEL1-VHS1) and 1:5 w/w (called MEL1-VHS5).
These prepared hybrid complexes were characterized by X-ray diffraction (XRD), comparing with the initial MEL solid (Figure 1a). MEL showed the most intense reflection at around 16.5° 2θ and other reflections at 19.0°, 24.5°, 25.5° and 26.5° 2θ according with those previously reported [78]. These peaks have also been found in the spectrum calculated with the Reflex Powder Diffraction module [57] of the experimental crystal polymorph of melatonin ( Figure S2) [47]. After the interaction of MEL with VHS (Figure 1b), the peak of MEL was observed at 16.5° (2θ units). This peak was observed at the oriented aggregate of XRD with a special high intensity in the MEL1-VHS1 interaction product. This higher intensity peak could be because of a preferred orientation of MEL and a higher proportion of MEL than the MEL1-VHS5 interaction product and also due to the complete interaction of the drug with the clay in the interaction product MEL1-VHS5. These results corroborated that the MEL drug is present in the interaction product. However, part of the initial MEL has precipitated as a crystal solid mainly in MEL1-VHS1 on the external surface of clay mineral, because cannot precipitate inside the interlayer confined nanospace. Nevertheless, the observation of only one reflection of MEL indicates that is not well crystallized. In MEL1-VHS5, the appearance of a very low intensity single peak and the disappearance of the rest of the peaks of MEL were observed, because the MEL intercalated in the interlayer space loses its crystallinity.
The oriented aggregate X-ray diffraction pattern of MEL1-VHS1 showed a shift of the (001) reflection of VHS after the adsorption of MEL from d(001) = 1.14 nm to 1.55 nm (Figure 1b), which indicates the intercalation of MEL on the interlayer space. However, this reflection peak showed a complex broad shape that could be considered as an overlap of several peaks corresponding to different d(001) spacings from several intercalation grades. Similarly, the increase in the interlayer space of the montmorillonite in the interaction product MEL1-VHS5 was also observed, being the d(001) shift = 1.65 nm slightly higher than above (Figure 1b). This loss of crystallinity of the drug in the interaction product with more clay (MEL1-VHS5) corroborates the X-ray diffraction results, as observed in previous work with other drug [42]. This amorphization is slightly higher in MEL1-VHS5 product. The TGA profile of MEL showed a thermal decomposition of the drug in the 200-575°C range (Figure 3a). The VHS profile showed a weight loss around 2% w/w, corresponding to evaporation of water ( Figure 3b) [42,45]. In the interaction products, MEL-VHS complexes, thermal decomposition of the drug resulted in weight losses in the range between 200 and 700°C, corroborating the adsorption of MEL in both cases. In addition, the TGA confirmed that the drug is in lower proportion in the interaction products MEL1-VHS5 than in MEL1-VHS1 (Figure 3b). The MEL1-VHS1 sample showed a fast weight loss at 260-300ºC similar to the pristine MEL profile. This indicates that part of the adsorbed MEL is in the external surface of the mineral solid. However, the slope of these profiles changes at T>320ºC, indicating the decomposition of MEL intercalated in the interlayer space of the mineral (Figure 3b). This behaviour could show a certain additional stability of MEL inside of mineral nanospaces.

In vitro release studies
In vitro drug release profiles of the pristine MEL and the MEL1-VHS5 interaction product are shown in Figure 4 in percentage of drug released. The results showed that the pure drug was not released or diffused through the membrane under the tested conditions. Similar previous findings indicated that the amount of pure MEL permeated is extremely low [79,80]. Conversely, MEL1-VHS5 interaction product revealed the release and diffusion of the drug, demonstrating a controlled diffusion rate in Franz cells in isotonic saline solution pH 7.4 as a medium simulating the plasma. In the MEL1-VHS5, approximately 50% of the drug was released between 4-5 hours and the 100% of the drug release through the membrane was achieved at around 10 hours. This corresponds to a total permeated amount of drug of 4600 µg/cm 2 . Therefore, the MEL1-VHS5 showed a slow dissolution rate and diffusion of the drug through the skin simulating membrane. This is typical of a prolonged release system, where the drug is gradually released over time at a rate limited by the release system (adsorbed clay system). With this design, the plasma concentration of the drug is prolonged, thus extending the therapeutic effect. Hence, drugclay interaction products in a 1:5 w/w ratio can be new modified release systems of MEL for skin delivery. Then to analyze the drug release, the dissolution data of the PZQ1-VHS5 interaction product was fitted to different kinetic models. In table 1, the adjust correlation coefficient (R 2 adjust), taking into account the same number of data (not data greater than 60% of the drug release to adjust a kinetic following Higuchi and Peppas requirements) and the same number of parameters, was calculated. The Akaike Information Criterion (AIC) was also calculated. R 2 adjust and AIC showed that the first-order model is the most appropriate to described the dissolution kinetic, because it presented the highest R 2 adjust and the lowest AIC value compared with the rest of the proposed models. Therefore, the release rate of the MEL1-VHS5 interaction product is concentration dependent with a first order kinetics [81]. 3.3. Molecular modeling and infrared spectroscopy of the crystal structure of melatonin The molecule of MEL was previously studied with different molecular modeling methods by several authors [1,82,83]. However, thorough studies for the crystal structure of MEL are still missing. This is interesting because the experimental IR spectroscopy of MEL is performed normally in its crystalline phase instead of as isolated molecules. Nevertheless, we start our computational study with the isolated MEL molecule in order to compare and validate different methods prior to our extended study in the crystal lattice.
The main geometrical features of the optimized MEL molecule ( Figure S3) were compared (Table S1). The distances of the main bond lengths and angles of the optimized MEL molecular structure with the different methods showed that the calculated values are consistent with the experimental data. The optimized structure is quasi planar as seen in Figure S3, except in the cases of the Compass FF and the INTERFACE FF. In the former, the H-N-C-C and C-C-C=CN dihedral angles are approximately 90º, and therefore the N-H and C=CN bonds are not coplanar with the CH2-CH2 aliphatic chain. In the later, the H-N-C-C and C-C-C=CN dihedral angles are approximately 0º and 90º, respectively. The observed conformation of the methoxy group with respect to the rest of the molecule is in agreement with DFT results in the literature [84].
The calculations of the crystal structure were performed with a unit cell containing four MEL molecules (194 atoms) arranged in the P21/C space group as experimentally known. The crystal structure is anhydrous. Initially, the MEL crystal structure was optimized at variable volume with the same FF methods previously used in the isolate molecule and with CASTEP quantum mechanics methods (Table S2), relaxing both the atomic positions and the lattice parameters ( Figure 5). The packing energy, calculated at DFT level was -54.77 kcal/mol per unit cell, i.e. -13.7 kcal/mol per molecule. The main bond lengths and angles in the optimized crystal structure of the drug were similar to those described above in all methods, being close to the experimental geometry. In the MEL crystal lattice, the indole rings are parallel but are shifted and no significant π-π interaction is observed. However, the carbonyl group has strong bifurcated H bonds with the  In the Table 2, the crystal lattice parameters of the optimized MEL structure with all atom positions at variable volume and different methods were compared with experimental data. According to these results, the two best methods for reproducing the experimental lattice parameters were INTERFACE among the FF methods, and CASTEP with Tkatchenko-Scheffler dispersion correction among the DFT methods ( Figure 5), whereas the Universal FF is the method that experimental data are worst reproduced. Therefore, we selected the INTERFACE FF and CASTEP with Tkatchenko-Scheffler corrections for performing the IR spectroscopic and adsorption studies. IR frequencies of crystal structure of MEL were calculated (Table 3) and compared with the experimental IR spectrum obtained by us as described in the Supplementary Materials ( Figure 6 and Figure S4). IR frequencies calculated with INTERFACE and DFT with CASTEP were consistent with our experimental results, and some new frequencies have been assigned. No negative frequencies were found, indicating that these structures correspond to minima in the potential energy surface. Nevertheless, the ν(NH) vibration modes calculated with INTERFACE showed high frequencies, whereas the values calculated with DFT fitted better with experiment, though both methods showed similar tendency. These theoretical calculations showed that the stretching of the NH bond, which is close to carbonyl, appears at a higher frequency, while the stretching of the pyrrole NH bond appears at lower frequencies. This assignment is different to that reported previously, [84] where an isolated molecule was considered. In our crystal model the ν(NH) frequencies are closer to experimental values than the previous one [84] because the intermolecular interactions are included in our crystal model. Both NH groups form H bonds with the carbonyl O atoms of the vicinal molecules. The effect of the H bond decreases the frequency of ν(NH) and this effect is higher in the pyrrole NH group, d(NH···OC) = 1.847 Å, than in the amide NH group, d(NH···OC) = 1.847 Å. Hence, ν(NH) of pyrrole group appears at lower frequency than ν(NH) of amide group, the opposite of previously assigned [84]. These bands cannot be assigned only experimentally. Along with these bands, the stretching of the CH bond of the pyrrole ring, asymmetric stretching of the CH3 bonds, bending of the CH bond of the pyrrole ring, bending of the NH bonds, etc., have been newly assigned in this work, because they could not be experimentally assigned properly. The pyrrole ν(CH) mode appears at higher frequencies than the ν(CH) of the benzene aromatic ring. We can appreciate some shoulders at lower frequencies of the ν(NH) band, that can be assigned to this pyrrole ν(CH) mode. The ν(CH) of the rest of CH groups are coupled and the assignment has been tentative.
The ν(C=O) vibration mode calculated with INTERFACE showed higher frequencies than experimental values, probably due to an overestimation in highly polar bonds like in ν(NH). Whereas the ν(C=O) frequencies calculated with DFT fitted better with experiment. The pyrrole δ(NH) modes appear at lower frequencies than the amide δ(NH) modes, though both cases are mixed with some δ(CH) modes.

Molecular modeling of the adsorption of melatonin on montmorillonite
Computationally the adsorption of MEL in the interlayer space of montmorillonite in presence of water molecules was studied. First of all, one MEL molecule was placed in the centre of the interlayer space of VHS in a 3x2x1 supercell, Na6(Al19Mg5)(Si47Al1)O120(OH)24, and twelve molecules of water per supercell were included around the adsorption complex. The adsorption complex (1MEL-VHS) was optimized by using INTERFACE FF (Figure 7) at variable volume yielding a d-spacing d(001) = 12.3 Å, slightly higher than the optimized d(001) = 11.6 Å of the VHS without drug absorbed.
After the optimization, the intercalated molecule adopts a linear and quasi-planar conformation, due to the electrostatic interactions with the mineral surface. The aromatic part of MEL adopts a parallel orientation with respect to the interlayer mineral surface. The aliphatic chain of the molecule presents a certain twist, forming an angle of 107º with respect to the aromatic ring.
The main interactions between MEL and VHS were hydrogens bonds and electrostatic interactions. MEL presented weak hydrogen bonds between its carbonyl O atom and the water H atoms with d(COMEL…Hw) = 2.  The adsorption is energetically favourable as shown by computing ΔEads = (E1MEL-VHS) -(EMEL + EVHS) = -18.6 kcal/mol. However, the d(001) spacing is lower than the experimental one obtained above by us in MEL-VHS hybrid systems. We infer from this fact that either a higher amount of MEL could be adsorbed in the interlayer space of VHS in the experimental data or MEL is adsorbed in a different conformation than that found in our calculations.
In a second calculation, two MEL molecules were adsorbed in the interlayer space of the VHS. The adsorption complex (2MEL-VHS) was optimized with INTERFACE FF at variable volume resulting a interlayer space of d(001) = 13.5 Å. This is larger than the space 1MEL-VHS system. If we consider the two molecules get a quasi-planar conformation, similar to 1MEL-VHS, this interlayer space can be also justified by water molecules intercalated between MEL and the VHS internal surface (Figure 8a). After optimizing this complex, the adsorbed drug molecules are found as a monolayer in the interlayer space of the clay (Figure 8). The main interaction forces between MEL and VHS were also hydrogens bonds and electrostatic interactions as was previously described. The adsorption is again favourable, with ΔEads = -36.3 kcal/mol, approximately twice the value of a single molecule. Nonetheless, the adsorption of two molecules loses approximately 0.9 kcal/mol with re-Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 19 July 2021 doi:10.20944/preprints202107.0395.v1 spect to a linear interaction, which could come from either the repulsive interaction between both molecules, a change of conformation of one or both molecules at once, and/or the interaction with the water molecules between MEL and the interlayer surfaces. Finally, three molecules of MEL drug were intercalated in the VHS interlayer space. This optimized adsorption complex (3MEL-VHS) with INTERFACE at variable volume showd a d-spacing d(001) = 16.0 Å, similar to the experimental one, suggesting that this could be the actual MEL-VHS ratio in the interlayer space. Nonetheless, the three molecules confined in the interlayer space of our model interact producing different conformations than the interaction of two molecules. In effect, in the optimized structure, the MEL molecules are oblique to the interlayer plane as showed in Figure 9. The main interactions between the drug and the clay were hydrogens bonds and electrostatic interactions, as those previously described. The adsorption energy in this case was -48.4 kcal/mol, approximately 7 kcal/mol less than three times the adsorption of a single molecule. In effect, this indicates in the model with three a certain repulsion between them is produced, coming from the lack of room in the interlayer space, and the three interactions previously mentioned in the case of two. Taking the point of view of the adsorption energy per molecule in each of the cases (i.e. ΔEads/n being n = 1, 2, 3) is -18.6, -18.15 and -16.13 kcal/mol for the 1, 2 and 3 MEL/supercell cases respectively, showing an extremely small negative cooperative effect in the adsorption of MEL on VHS, coming from the three effects previously mentioned. 3.5. Ab initio Molecular Dynamics As exposed in the Methodology section, in this ab initio molecular dynamics (MD) we initially placed a MEL molecule in the interlayer space of the montmorillonite model and many water molecules in the interlayer space and outside. Multiple Walkers Metadynamics simulations were run until observing the diffusion of the MEL from the interlayer space to the bulk water reservoir. The different replicas (walkers) of the system sampled many different configurations of MEL in the interlayer space until the diffusion was eventually obtained (Figure 10). The simulation allowed us to reconstruct the Free Energy Surface (FES) associated to the diffusion process from the confined MEL to the solution (Figure 11). MEL looks like to be visiting different minimum positions around a global minimum region being most of them secondary minima, but all of them are in a region between z = 2 and 4 Å, but when MEL goes out from the minima region, FES trajectories ride in a wide FES plateau up to fall into the solution. We can estimate the highest free energy barrier as  20.5 kcal/mol. On the other hand, it was not possible to simulate the re-entry of the MEL from the bulk region to the interlayer space due to the high computational cost associated.  . Free energy surface at 300 K for the release of the MEL from the montmorillonite interlayer space to the bulk water reservoir. The global minimum is seen around x=5.5, z=2.5, where the MEL is in the interlayer space, while the MEL is released to the bulk water region when x>17.
More interesting is the sampling of the different configurations of the adsorbed MEL within the interlayer space. In contrast to what can be seen from the static picture provided by the 0 K geometric optimizations, the finite-temperature dynamics reveals that the MEL adopts many different conformations in the interlayer space. Two representative such configurations consistent with the global minima in the CV space are shown in Figure 12. In order to reach a deeper understanding of the interactions that determine the retention of MEL in the montmorillonite (MON) system, we computed the bonding charge density in these two structures as Δρ = ρfull -ρMEL -ρMON, where ρfull is the electronic density of the full system, ρMON is the electronic density of the system without the MEL molecule, and ρMEL is the electronic density of the isolated MEL in the gas phase. Visualizing the resulting Δρ in real space, we can clearly see which are the regions with electron accumulation/depletion upon the insertion of the MEL in the MON interlayer, thus providing an intuitive picture of the interactions involved. In the configuration shown in Figure  12a, the MEL lies with its indole ring parallel to the montmorillonite surface, on top of two Na + cations, which are introduced in the tetrahedral cavity, and interact considerably not only with the methoxy oxygen but also with two C atoms of the indole ring, which are seen to be polarized. The Na + interacting with the indole group is much deeper in the tetrahedral cavity. The cations seeking property in the tetrahedral cavity make the adsorbed molecules have more room in the interlayer space, and tetrahedral cavity···cation system be a more neutral charged locus. On the other hand, the longest side chain interacts with the topmost montmorillonite layer via the amino group and even some of the CH units, while the keto group accepts two H bonds from the surrounding water molecules.
Another representative configuration (Figure 12b) of the main minima shows the MEL molecule obliquely across the whole interlayer space, with the Na + cations at one side interacting with the ether oxygen and the aliphatic H atoms of the side chain, while at the same time the C-H and N-H units of the indole ring and of the side chain in touch with the topmost montmorillonite surface are seen polarized. Again, the keto group is involved in H bonding with the neighbouring water molecules. We note that in some other configurations, it is possible to see the keto group also interacting with the interlayer Na + cations. It is possible to note, some water molecules inside the tetrahedral cavity with the hydrogen pointing to the basal oxygens, and going up and down because of the dynamics. Figure 12. Representative structures of the MEL in the interlayer space sampled from the AIMD metadynamics simulations at 300 K (side (a) and top (b) views). Colour code as in Figure 10, except for the montmorillonite surfaces, which are shown in silver van der Waals spheres for clarity. The isosurfaces of the bonding charge density Δρ are shown in orange (indicating electron accumulation) and blue (electron depletion) corresponding to the isovalues +0.004 e-/(a0)3 and -0.004 e-/(a0)3, respectively. For clarity, in b we only show the water molecules that are within 4.0 angstroms around the MEL molecule.
For the sake of comparison, we run a regular AIMD simulation of the MEL molecule in bulk water and analyze some of their representative solvation structures ( Figure 13). In this case, the keto oxygen is also involved in different H bonds with its surrounding waters, as it is also the case of the amino group on the side chain and the NH of the indole ring. On the other hand, the interactions of the methoxy oxygen with the solvating waters are weaker than those observed inside the montmorillonite with the Na + cations, and the CH bonds as well as the indole ring itself are seen much less polarized than in the montmorillonite. This weakness in the interactions of the water molecules of the solution with MEL could justify the lack of solubility of MEL in water. However, the special stability of MEL in the interlayer space comes from the interaction with the interlayer environment.
We conclude that the surface charge on the montmorillonite strongly determines the adsorption of MEL in its interlayer space, thanks to the variety of the possible electrostatic interactions between MEL and the interlayer cations and the basal atoms on the montmorillonite surfaces. Thus, it would be interesting to explore the possibility of tuning the retention/release ability of the system by varying the concentration or composition of the interlayer cations (e.g. exchanging Na + by K + ), or by introducing substitutions in the tetrahedral sheets (e.g. Al 3+ by Si 4+ ). Figure 13. Representative structures of the MEL in bulk water sampled from AIMD at 300 K. Colour code as in Figure 10. For clarity, we only show the water molecules that are within 4.0 angstroms around the MEL molecule.

Conclusion
The intercalation of melatonin in the nanospaces of montmorillonite was demonstrated by observations using different experimental techniques. Remarkably, the results showed that the melatonin-clay system gradually releases the drug at a regular rate, reaching the diffusion of 100% of the adsorbed drug in 10 hours. On the contrary, pure melatonin does not permeate the skin-like membrane under these in vitro conditions. Therefore, melatonin-clay systems enhanced the diffusion and release of the drug and can be considered a new dosage form for topical application.
On the other hand, the molecular and crystal structure of melatonin have been studied theoretically by using different techniques based on classical and quantum mechanical methods. Among the different methods based on molecular mechanics, the INTERFACE force field showed the best performance, while among the quantum-mechanical ones it was DFT with Tkatchenko-Scheffler corrections, reproducing the experimental lattice cell parameters of the melatonin crystal structure. These theoretical methods are useful tools for assignment of experimental infrared frequencies of melatonin, which are performed mainly in solid state, where the intermolecular interactions can alter its frequencies. This approach is more realistic than previous works calculating isolated molecules models in gas phase. With these methods we have successfully studied the IR spectra of melatonin, in excellent agreement with the one obtained experimentally in this same work. Our approach has facilitated the interpretation and assignment of new IR bands of melatonin, assigning IR bands that had not been previously described.
Finally, we have explored the structural and energetic features of the melatonin-clay system via static calculations at the classical and ab initio level, and also via finite temper-ature ab initio MD simulations, describing at atomistic level the release process of melatonin from clay minerals and characterizing the most important interactions taking place in the melatonin-montmorillonite system, finding that this adsorption is energetically favourable and is determined by coulombic interactions and hydrogen bonds with the basal oxygen of the tetrahedral sheet of the interlayer space. Interlayer cations take part of these electrostatic interactions, but they have a mobility role by being able to introduce in the tetrahedral cavity and doing room for the molecule to adsorb. The diffusion of MEL in the confined interlayer nanospace and its release to another medium is confirmed by AIMD simulation at room temperature, where the MEL is stable in different intercalation minima, but its budging is able to found different trajectories up to overcome a FES 20 kcal/mol transition state to go to a trajectories flat plateau where its way out of the interlayer space is easy to found. These facts account for the results found experimentally.  Informed Consent Statement: Any research article describing a study involving humans should contain this statement.