Microscopic Interactions of Melatonin, Serotonin and Tryptophan with Zwitterionic Phospholipid Membranes

The interactions at the atomic level between small molecules and the main components of cellular plasma membranes are crucial for elucidating the mechanisms allowing for the entrance of such small species inside the cell. We have performed molecular dynamics and metadynamics simulations of tryptophan, serotonin, and melatonin at the interface of zwitterionic phospholipid bilayers. In this work, we will review recent computer simulation developments and report microscopic properties, such as the area per lipid and thickness of the membranes, atomic radial distribution functions, angular orientations, and free energy landscapes of small molecule binding to the membrane. Cholesterol affects the behaviour of the small molecules, which are mainly buried in the interfacial regions. We have observed a competition between the binding of small molecules to phospholipids and cholesterol through lipidic hydrogen-bonds. Free energy barriers that are associated to translational and orientational changes of melatonin have been found to be between 10–20 kJ/mol for distances of 1 nm between melatonin and the center of the membrane. Corresponding barriers for tryptophan and serotonin that are obtained from reversible work methods are of the order of 10 kJ/mol and reveal strong hydrogen bonding between such species and specific phospholipid sites. The diffusion of tryptophan and melatonin is of the order of 10−7 cm2/s for the cholesterol-free and cholesterol-rich setups.


Introduction
The cell membrane plays a central role in the control of the exchange of key elements (nutrients, wastes, drugs, and heat as the most relevant) between the exterior of a cell and its cytoplasm. Lipids, proteins, and cholesterol (CHOL) are among the main components of human cell membranes. Phospholipids are usually formed by two leaflets of amphiphilic lipids that are divided into a hydrophilic head and one or two hydrophobic tails. Such lipids can self-assemble by hydrophobicity [1,2]. Lipid bilayer membranes formed by dipalmitoylphosphatidylcholine (DPPC, C 40 H 80 NO 8 P) and dimyristoylphosphatidylcholine (DMPC, C 36 H 72 NO 8 P) are of great interest to computational studies [3][4][5][6] because of the abundance of experimental data [7][8][9][10][11] and their usability as model systems [12,13]. In the last decades, cell membrane systems have been extensively studied in regards to their association with drugs and small molecules [14][15][16][17][18]. Among these, we will focus our attention to three species deeply related, belonging to the indole group: melatonin (MEL, C 13 H 16 N 2 O 2 ), its precursor serotonin (SER, C 10 H 12 N 2 O), and its precursor tryptophan (TRP, C 11 H 12 N 2 O 2 ), given their interest as pharmaceutical compounds with a wide variety of applications, in particular a large number that is related to sleep disorders [19][20][21]. The relationships between the three components have been reviewed a significant number of times in different contexts [22][23][24][25]. The sequence of transformations in vivo is as follows: first, the hydroxylation of L-tryptophan by tryptophan hydroxylase forms 5-hydroxytryptophan; when followed by the decarboxylation by the aromatic aminoacid decarboxylase, we get serotonin; have employed two types of computational methods, molecular dynamics (MD) and well tempered metadynamics (WTM). MD is a classical simulation tool that is able to generate a bundle of Newtonian trajectories, one for each single particle of the system, at the atomic level. As atoms interact through pairwise force fields, their trajectories (composed of positions and linear momenta of all particles) are deployed and stored at regular time intervals, in order to be analysed using tools from Statistical Mechanics [68]. MD is a versatile method that is able to successfully reproduce a large number of microscopic properties of a wide variety of systems, from simple atomic liquids, such as argon [69] to molecular liquids as water [70,71], aqueous solutions at interfaces [72][73][74][75][76], up to complex biophysical systems like DNA [77][78][79] or model cell membranes [6,[80][81][82][83][84]. In order to handle the problem of computing free energy landscapes in multidimensional systems, different classes of methods have been proposed, such as quantum mechanics/molecular mechanics [85], transition path sampling [86][87][88][89][90][91][92], adaptive biasing force [93], umbrella sampling methods [94,95], density functional theory molecular dynamics [96], or calculations of potentials of mean force [97] based on reversible work methods [98]. In this work, we have employed reversible work and WTM, a method that is able to efficiently explore free energy surfaces of complex systems while using multiple reaction coordinates what has been revealed to be very successful [99] for a wide variety of complex systems [100][101][102][103][104]. Section 3 reports the technical characteristics of all simulations.

Structural Properties of the Membranes
The structural characteristics of the membranes and the local distributions of atomic species are the first group of properties to be analysed. To do so, we have sketched the detailed atomic structures of the three small molecules, the two phospholipids composing the membranes (DMPC, DPPC) and CHOL in Figure 1. There, the highlighted sites of TRP are the zwitterions 'H1', sharing a positive charge between the three hydrogens that are bound to 'N1'; 'H2', bound to 'N2' and the zwitterions 'O1' and 'O2', bound to 'C1' and sharing a negative charge. In SER, we highlight 'H1', 'H4', and 'O', and, for MEL, we will keep special attention into 'O1', 'O2', 'H15', and 'H16'. Finally, the sites 'O1' and 'O2' sharing the negative charge and oxygen atoms 'O6', 'O8' will be considered for DPPC and DMPC and the hydroxilic OH pair for CHOL.
A common test in computer simulations of cell membranes is the comparison of the area per lipid and thickness of the membrane with experimental data from scattering density profiles [105]. We have monitored the surface area per lipid A when considering the total surface along the XY plane (plane parallel to the bilayer surface) that is divided by the number of lipids N l plus the number of cholesterol molecules N chol. in one lamellar layer [106], as defined in Equation (1): where L x and L y are the length of the simulation box along X-axis and Y-axis, respectively. Z-axis is the (instantaneous) normal direction to the surface of the bilayer, set along plane XY. Fluctuations in the thickness of the membrane are related to the effect of cholesterol on the rigidity of the membrane and its capability to allow the passing of species in and out of the cell. In this work, we defined the thickness ∆z as the distance between the phosphate groups of the lipids at the two sides of the membrane. The area per lipid and thickness along the last 500 ns of each simulation have been computed (see Figure 2) and Table 1 reports the average values.  Figure 1. Atomic structures of melatonin, serotonin, tryptophan, dipalmitoylphosphatidylcholine (DPPC), dimyristoylphosphatidylcholine (DMPC), and cholesterol. Backbone hydrogens are not explicitly shown. The highlighted labels will be referred in the text.  The area per lipid decreases as cholesterol concentration increases: this is a well known trend, as will see below. We obtained a value of around 0.61 nm 2 for a cholesterol-free system and smaller values down to 0.40-0.42 nm 2 when cholesterol has been incorporated. In all cases, the area per lipid is practically independent of the small molecule that is imbedded in the membrane, and it has little influence of the main type of phospholipid. These results are in excellent agreement with experimental data [109,110], where the value for pure DMPC is of about 0.6 nm 2 at 303 K. Further, and according to Nagle et al. [1], values of A of pure DMPC membranes can be obtained from multiple methods (neutron scattering, X-ray and NMR) and they have been reported to be between 0.59 and 0.62 nm 2 at the liquid crystal phase. In the case of DPPC, the best estimations were of between 0.48 and 0.52 nm 2 in the gel phase (293 K) and 0.64 nm 2 in the liquid phase. These results are also in overall good agreement with other computational data in a wide variety of thermodynamical conditions [11,109,[111][112][113], where the values for pure DPPC ranged between 0.50 and 0.63 nm 2 and the trend of decreasing areas for increasing cholesterol percentages was clearly reported. The huge change that is produced at 30% cholesterol concentration is consistent with the fact that phosphatidylcholine membranes experience a phase transition liquid disordered (cholesterol-free system) to liquid ordered phase (systems of cholesterol 30% and 50%) [114,115].
The thickness of the membranes are in good agreement with those that were reported by Kucerka et al. [110] by means of X-ray and neutron scattering. The reported value was of 3.67 nm at 303 K for the DMPC membrane at 0% cholesterol. From the results that are reported in Table 1, we obtain values around 3.5-4 nm for pure bilayers and of 4.4-4.9 nm when cholesterol is considered. We observe a tendency to larger bilayer thickness for increasing cholesterol concentration. Given the reduction of the area per lipid at high cholesterol percentages, we can conclude that cholesterol favours the compression of the structure of the bilayer membrane. This feature can increase the rigidity of the membrane and, by extending the tails of the lipids, give larger bilayer thickness. Such an increase of the rigidity of the membrane was observed from both experimental and computational sides [60,61] in MEL-CHOL mixtures nearby phosphatidylcholine bilayers. According to these studies, the effect of MEL reducing the thickness of the membrane and enhancing its fluidity was partially compensated by the condensating effect of cholesterol.

Preferential Localisations of the Small Molecules at the Interfaces of Phospholipid Membranes: Atomic Radial Distribution Functions
Each of the three small molecules considered in the present work has been simulated for long MD trajectories of hundreds of nanoseconds. We have monitored their positions and velocities and obtained structural, energetic, and dynamical information. In this section, we will focus our attention on the local structure of the probes when embedded in the membrane. As a general fact, we have observed that all three selected species show a strong tendency to be continuously adsorbed at the interface of the membrane during long periods of the order of 10 ns. In the remaining time, the small molecules move away to be solvated by the ionic solution surrounding the membrane. As an example, in Figure 3 we report the evolution in time (window of 60 ns) of the position of the center of mass of melatonin when adsorbed at a DMPC-cholesterol membrane.
In Figure 3, we can observe that the influence of cholesterol is of paramount importance: when the concentration of cholesterol in the membrane reaches 50% of all lipids, MEL can easily shift between the interface of the membrane and the solvating aqueous ionic solution, but, at lower concentrations, the small molecule is likely inside the membrane during the whole time span considered. This indicates that moderate changes of cholesterol concentration may induce some specific organic probes to retreat from the inside of cellular membranes to the outer regions and remain outside the cell. This might have strong implications in melatonin delivery. The relationships between MEL and CHOL and their interactions have been studied since long time ago [40,60,61], but the knowledge of their effects are still quite elusive.  Because the computation of radial distribution functions (RDF) is the best way to investigate atom-atom local structures, we have computed a series of specific RDF in order to have an overview for TRP and MEL. We define the RDF for an atomic pair composed by particles '1' and '2' as g 12 (r), and it is given by: where n 2 (r) is the number of atoms of species '2' surrounding a given atom of species '1' inside a spherical shell of width ∆r. V stands for the total volume and N 2 is the total number of particles of species '2'. In the case of TRP, we have considered the partial RDF that is reported in Figure 4, whereas, for MEL, we will analyse the RDF that is presented in Figure 5.
From the data that are reported in Figure 4, we can observe that TRP stays bound to the inner part of the membrane during long periods of time, according to the time scale of our simulations, in good agreement with the results indicating that, in a cholesterol-free DOPC bilayer membrane, TRP is preferentially located in the interfacial region [55]. In this work, we have observed that the average continuous lifetime of TRP at the interface of the DPPC bilayer is of the order of 10 ns (data not shown). From Figure 4, hydrogen bond (HB) connections between sites 'H1' and 'H2' of TRP and DPPC sites 'O2' and 'O8' (labels according Figure 1) have been found. HB are very short, since the maxima of the g(r) related to 'H1' hydrogens in TRP are located around 1.7 Å for 'O2' and around 1.75 Å for 'O8' of DPPC. Accordingly, the presence of cholesterol reduces 'H1-O2' binding, but enhances the 'H1-O8' one. As a general fact, the presence of cholesterol increases the length of HB, but also making such bonds stronger. This indicates that the influence of the cholesterol in the TRP-DPPC binding is a major effect. Interestingly, 'H2'-DPPC binding was observed in all of the analysed setups, but with maxima found at larger distances (1.9-2.0 Å). Again, the presence of cholesterol showed a major influence on the characteristics of hydrogen bonding. Figure 5 reports the structural results for MEL. All RDF show fluctuating profiles, especially at distances that are larger than r = 3 Å and beyond (higher order coordination shells). There is a clear first coordination shell in all cases, located around 1.8-2.0 Å, due to HB between MEL and the remaining species, such as in the TRP case. The largest maximum of all RDF is the one for MEL-CHOL association (not shown here), which is centered at 1.9 Å when the concentration of cholesterol is of 30%. Choi et al. reported interactions of MEL-CHOL in DPPC bilayers [61], but at finite melatonin concentration. In the remaining cases, HB lengths are around 1.9 Å and they were between both 'H15' and 'H16' of MEL and DMPC sites 'O1' (or 'O2', since both of the sites are sharing the negative charge of the zwitterion). In a similar fashion, MEL can also form HB between both 'H15' and 'H16' with the DMPC's sites 'O6' (or 'O8') for all three percentages of cholesterol. The present findings are in good agreement with the experimental data from Severcan et al. [50] that were obtained by Fourier transform infrared spectroscopy, who observed hydrogen bonding connections between the N-H group of the furanose ring of MEL ('H16' in this work) and the carbonyl (C=O) and phosphate (PO 4 ) groups in DPPC membranes. Our results indicate HB between 'H15' and 'H16' of MEL with DMPC's phosphate group ('O1') as well as with the more internal carbonyl groups ('O6' and 'O8'). The previously unobserved hydrogen bonds of 'H15' with the two well known acceptor groups in the phosphatidylcholines indicated above are responsible for the absorption of MEL into the membrane deeper than TRP, with the two selected donors ('H15' and 'H16'), together with the 'H15-O Chol.' bridges. These findings are in excellent agreement with the results that were reported by Drolle et al. [60] by means of small angle neutron diffraction and MD simulations.

Orientational Distributions of Melatonin
Several previous studies have shown that the orientations of drugs on membranes significantly impact their function in cells [116][117][118][119][120]. In the present work, we have computed the principal orientations of MEL through the definition of three different dihedral angles. We have observed that, in all cases, two preferential orientations arise, since the averaged angular distributions of MEL are centered around two well defined angular values, which we call "folded" and "extended" configurations of MEL, found at all cholesterol concentrations. The dihedral angle, which has a better distribution regarding its fluctuations around mean values, is the angle Ψ that is represented in Figure 6. The torsional angle considered here is related to the nitrogen atom labelled 'N1' in Figure 1, namely the nitrogen chemically bound to the hydrogen labelled 'H15'. For this angle Ψ, after analysing 100 ns of equilibrated trajectories (production runs), we found averaged values that corresponded to 81 ± 10 • , (folded) and 170 ± 23 • , (extended), as shown in Figure 7. Further, from the distributions reported there we can observe that Ψ is neatly defined and it reaches nearly the same mean value, regardless of the concentration of cholesterol of the system. Interestingly, we find that the extended configuration of MEL is most favoured in the case of the highest concentration 50% (green triangles), which suggests that introducing cholesterol into the system could help MEL change from its folded to its extended configuration more easily through hydrogen-bonding between MEL-DMPC and MEL-cholesterol. In addition, according to this, Ψ is an excellent candidate for being used as a collective variable in metadynamics calculations [121,122] of free energy landscapes for MEL binding in biomembranes, as we report in Section 2.4.

Free Energy Profiles of Small Molecules and Free Energy Hypersurfaces of Melatonin Binding
Once we have established preferential locations and angular distributions of the small molecules, if assuming some of these coordinates as good candidates for collective variables, we are ready to use the WTM technique to obtain precise, quantitative values of the free energy barriers that need to be surmounted by the small molecules to move throughout the system, mainly exchanging positions between the interfacial regions and the bulk like aqueous regions of the system. Computationally speaking, WTM is a very expensive method that requires very long trajectories, so that the target subsystem, i.e., the small molecule, can move in the full configurational space, visiting regions of low energy with high probability as well as regions of high energy, being very unlikely to be accessed. In this work, we will complement WTM with a much simpler technique, based in the knowledge of RDF described above, namely the computation of the reversible work that is needed for the target to move between selected regions, being indexed by one dimensional coordinate, such as a radial distance r. The theory has been nicely described in chapter 7 of Ref. [123]. It states that we can obtain W 12 (r) i.e., the reversible work (sometimes also known as potential of mean force, PMF) that is required to move two tagged particles from infinite separation to a relative separation r from: where β = 1/(k B T) is the Boltzmann factor, k B the Boltzmann constant and T is the temperature. W 12 (r) can be understood as the relative Helmholtz (canonical ensemble) or Gibbs (isothermal-isobaric ensemble) free energy that is associated to atomic pairing. Figure 8 reports the W(r) found for the three small molecules and Table 2 reports the quantitative estimations of the main energy barriers. Reversible work calculations can only give us a rough approach to the size of real barriers, since it is based in the use of the interparticle distance r as the only reaction coordinate, which is known to produce some underestimation [97]. However, given that accurate reaction coordinates are usually unknown, very hard to obtain, and multidimensional, W(r) is a reasonable way to estimate the order of magnitude of the free energy barriers. The sata reported in Table 2 and Figure 8 reveal to us that the highest barrier corresponds to the pairing of 'H1' of TRP with 'O2' of DPPC. We have found that all small molecules are able to establish HB with 'O2' and also with the site 'O8' of DPPC, with the latter being located deeper in the membrane (see Figure 1). Conversely, we did not find bindings between 'H15' site of MEL and 'O2' site of DPPC. The position of maxima of the first barrier are mostly centered around 2.45 Å for small molecule-'O2' binding, whereas barriers of ligand 'H4' of serotonin that are associated to the 'O8' sites are centered around 2.75 Å. In the case of SER, only a first minimum is clearly found, which indicates that SER is normally bound to the plasma membrane and it does not move to the solvent bulk.
The binding of 'O2' in DPPC to TRP is located at 1.75 Å, corresponding to the first minimum of the PMF between TRP and DPPC, which is of the order of the typical HB distance in water. Nevertheless, stable positions for 'O8' sites of DPPC are found between 1.7 and 2 Å, a remarkable wider distance. We can compare these values with the barrier for TRP (attached to a polyleucine α-helix) inside a DPPC membrane of 12 kJ/mol [52] or the barrier of the order of 16 kJ/mol found for TRP in a dioleoylphosphatidylcholine bilayer membrane [55]. Finally, the agreement of the barriers reported in the present work ( Table 2) with other neurotransmitters, such as glycine, acetylcholine, or glutamate, of around 2-5 kJ/mol when it is located close to the lipid glycerol backbone [124], is also quite remarkable. One way of getting much more precise free energy estimations is through methods operating with multidimensional reaction coordinates. One of best methods is well tempered metadynamics, although it is a very expensive computational tool, as we will explain in Section 3.2. As a specific example, we have applied WTM to the calculation of the hypersurface of free energy for the system that is composed by MEL and DMPC, at the three cholesterol concentrations of 0%, 30%, and 50% that are described in Section 3.1. The WTM specifications have been reported with full details in Section 3.2. We need to define several specific collective variables (CV) that are able to meaningfully describe characteristic configurations of MEL in order to compute the three sets of two dimensional (2D) well tempered metadynamics simulations. The results shown Figure A2 give us an indication of the convergence of WTM. To achieve this goal, we had to run trajectories of 1400 ns. These trajectories followed from the MD production runs that were employed to obtain structural and dynamical information. Figure 9 shows the resulting 2D free energy surfaces (FES) of MEL bound to DMPC membranes and they correspond to Gibbs free energy calculations. Each state has been indexed by two CV: (1) the z distance between the center of mass of MEL and the center of the membrane (z = 0); (2) the torsional angle Ψ defined and analysed in Section 2.3. The inspection of Figure 9 shows that regions with clear minima are present in the FES in all cases. The main features are the global minima of the FES located between z ∈ [0.7, 3] nm and around two distinctive values for the dihedral angle, namely those that are around |Ψ| ∼ [70 • , 180 • ]. Such orientations are in excellent agreement with the average values of Ψ that were obtained from ordinary MD simulations (see Figure 7) that correspond to the two "folded" and "extended" geometries of melatonin previously reported.
The 2D free energy landscapes reveal that the most favourable stable states of melatonin binding to the membrane (basins A,B,C,D) correspond to z-distances around 0.8 nm at the cholesterol-free system, whereas such a distance tends to significantly increase around to 1.3 nm for the 3% cholesterol concentration and up to 2.3 nm when cholesterol reaches 50%. As a general fact, the 2D surfaces that are shown in Figure 9 correspond to contour plots with values being referred to a global zero. The zero of each 2D plot has been set at the highest free energy value of all, in our case corresponding to locations at the computed maxima of the coordinate z.
According to the CV1, MEL is preferentially located at the interface of the DMPCcholesterol bilayer (regions with 0.8 < z < 3.0 nm). The locations of MEL outside the interface and far enough of lipid headgroups (z > 4.3 nm) show very larger free energies and they cannot be considered to be stable states of the system. Those regions will be considered as the "bulk", i.e., the region containing the electrolyte solution surrounding the membrane. When considering the information revealed by CV2, we can distinguish two sets of minima: (1) for |Ψ| = 67 • (basins B and C) and (2) for |Ψ| = 180 • (basins A and D. These minima are related to the two preferential configurations of MEL close to a DMPC-cholesterol bilayer (folded, extended) indicated above around 80 and 170 • (see Section 2.3).
We collected the data extracted from Figure 9 to estimate the main free energy barriers for the main configurational changes on MEL in the quantitative side. Table 3 reports the values. Our findings have revealed a rather wide range of absolute free energies, which are in good agreement with the range that was reported by Jämbeck and Lyubartsev [125] for small molecules (ibuprofen, aspirin, and diclofenac) at the surroundings of lipid bilayers, of the order of free energy ranges up to 70 kJ/mol and barriers around 40 kJ/mol. For the sake of comparison, we should remark that the barriers of 2-10 kJ/mol reported in Table 2 obtained from the PMF of Figure 8 were related to the formation and breaking of HB, when the small molecules were located inside the interfacial region, regardless of its orientation. However, free energy barriers of orientational changes or those that are related to large displacements of MEL to the center of the membrane or to the extracellular bulk are much larger. For instance (see Table 3), the free energy that is required to exchange between folded and extended MEL configurations is very stable, around 15-20 kJ/mol for all cholesterol concentrations. Florio et al. [126] using a combination of several fluorescence and spectroscopic techniques, found the conformational preferences of an isolated MEL molecule under molecular beams. These authors found MEL three trans and two cis conformers showing free energy gaps of approximately 12.5 kJ/mol, in quantitative agreement with the values reported here for orientational changes. However, the barrier that is to be surmounted by MEL to move from the interface of the membrane to the extracellular fluid is strongly dependent on cholesterol concentration. We observed that it decreases with larger amounts of cholesterol, between 25 kJ/mol at the cholesterol-free case to around 10 kJ/mol for the 50% concentration. Finally, the probability for MEL to access the central, hydrophobic regions of the membrane is scarce, since it will require surpassing free energy barriers of more than 40 kJ/mol. This will make it very difficult to observe transmembrane crossings in the simulated scale of 1 µs. In a recent work conducted by Wang and coworkers [127], small solutes, such as glycerol, caffeine, isopropanol, or ethosuximide, were simulated nearby a model cell membrane. These authors found that, in order to observe transmembrane crossings of such small solutes in the time length of a simulation at the atomic level of description, they needed to run trajectories of 10 µs at low temperatures (310 to 330 K) or, alternatively, raise the temperatures to more than 400 K (for simulation times of 1 µs). In our case, we did not record any transition of MEL between the two sides of the DMPC membrane along the simulated trajectory.

Diffusion Coefficients of Small Molecules: Tryptophan and Melatonin
Microscopic translational dynamics of tryptophan and melatonin have been considered. We have evaluated the mean square displacement (MSD) of the carbon 'C2' in TRP (see Figure 1) and of the center of mass of MEL. From the long time slopes of both MSD, we obtained the corresponding self diffusion coefficients D through the Einstein formula of Brownian motion: where r i (t) is the instantaneous position of particle i. In this general procedure, the spatial dimension of the diffusion regions d is considered. TRP and MEL showed lateral like diffusion (d = 2). Table 4 summarises the results. The main finding is that self-diffusion coefficients D for TRP (Table 4) are between 3-14 × 10 −7 cm 2 /s, which, even within the same order of magnitude, are significantly larger than those of the diffusion of DMPC molecules [6] (0.6 × 10 −7 in the absence of cholesterol). The main trend is that D increases for rising cholesterol concentrations. Overall, we find that the mobility of TRP is significantly higher than that of DMPC. Nevertheless, the effect of temperature is remarkable here, since, at complementary simulations at 310 K, TRP diffusion was of about 2 × 10 −7 , i.e., being significantly slower given the gel like state of the membrane in such a case.
In the case of MEL, D also shows a tendency to increase when cholesterol is mixed with DMPC, regardless of its concentration. In Table 4, at 30% cholesterol, the value of D for MEL is six times larger than the value of D of DMPC molecules in pure DMPC bilayer membrane systems [6]. This fact would suggest that its mechanisms of diffusion may be similar to those of an individual particle (such as in Fickian diffusion) and qualitatively different of those of lipids, whose diffusion was observed to occur in a sort of collective way, being associated in local groups of a few units (around 5-10 units) [6].

Molecular Dynamics
We have performed seven independent series of MD simulations for TRP, SER, and MEL in different environments (DMPC, DPPC, and different concentrations of CHOL, namely: 0%, 30%, and 50% for TRP and MEL, whereas, for SER, only the cholesterolfree membrane was simulated. Each system contains a total of 204 lipid and/or cholesterol molecules that were fully solvated by ∼5000-10,000 TIP3P water molecules and 17-21 sodium chloride pairs at the human body concentration (0.15 M), yielding a system size of about 40,000-60,000 atoms. Table 5 sumarises the characteristics of all simulations. Our MD inputs were created with the CHARMM-GUI web-based tool [128]. All of the systems were simulated at the the isobaric-isothermal ensemble. i.e., at constant number of particles (N), pressure (P), and temperature (T) conditions, with equilibration periods for all simulations being more than 200 ns. After equilibration, we recorded statistically meaningful trajectories of more than 600 ns. A typical size of the system was of 80 Å× 80 Å× 81 Å. The simulation time step was of 2 fs in all cases. Given its ability to reproduce area per lipid of DMPC and DPPC in excellent agreement with experimental data, the CHARMM36 force field [129,130] was used. All of the bonds involving hydrogens were fixed to constant length, allowing for fluctuations of bond distances and all sorts of angles for the remaining atoms. Van der Waals interactions were cut off at 12 Å with a smooth switching function starting at 10 Å. Long ranged electrostatic forces were taken into account by means of the particle mesh Ewald method [131], with a grid space of about 1 Å and updated every time step. The periodic boundary conditions were considered in each spatial direction.

Well Tempered Metadynamics
As we pointed out above, obtaining free energy profiles and estimating the height of the main barriers between stable states is a very difficult task in condensed matter systems [132]. In the present work, in Section 2.4 we have presented two possible pathways to do the job: (1) using a direct method that is based on the reversible work theorem, but knowing that it is, at its best, a first approach to the real barriers and (2) employing a more sophisticated tool, called "metadynamics", which, given a well chosen set of a few reaction coordinates (the collective variables), is able to provide a much more exact picture of the free energy hypersurface. Huber et al. [133] and Grubmüller [134] initially proposed the method and it was developed later on by Laio and Parrinello [99,121] as a method to explore multidimensional free energy surfaces as a function of a a priori unknown CV. Given some deficiencies of the original method, well tempered metadynamics [122,135] was introduced. In the present work, we have run 1.4 µs well-tempered metadynamics simulations in order to obtain Gibbs free energies of the binding states of MEL at phospholipid membrane surfaces that were made by DMPC lipids and CHOL in sodium chloride aqueous solution. Starting from the long trajectories generated by unbiased MD simulations for MEL-DMPC, we could make a reliable guess of two potentially appropriate CV. All of the metadynamics simulations were carried out by means of the PLUMED2 plugin [136,137] within the joint GROMACS/2018.3-plumed tool and they were performed in the NPT ensemble. Table 6 reports the particular details of the WTM simulations. The ussual periodic boundary conditions in all directions of space were considered.

Conclusions
The interactions of some small molecules with human cells are undoubtedly a relevant field of research. In particular, the hormone melatonin has an important role in the treatment of a wide variety of diseases and problems that are related to sleep. It works as a regulator of circadian rhythms and as an antioxidative. Further, its precursor serotonin is a neurotransmitter playing a key role in a variety of physiological processes and in the regulation of mood and cognitive learning. Serotonin is synthesised by the body from its precursor, the essential aminoacid tryptophan. Tryptophan is a zwitterion, with a protonated amino group (NH 3+ ) and a deprotonated carboxylic acid (COO − ), and it is used as an antidepressant. In the present work, we are reviewing a series of MD and WTM simulations of different lipid bilayer membranes in an aqueous ionic solution of NaCl with embedded small molecules. The calculations have been performed using the CHARMM36 force field. Among them, cholesterol at two concentrations (30% and 50%) has been considered together with the cholesterol-free reference systems in order to explore the influence of CHOL concentrations on the properties of the small molecules.
In a preliminary study on the adsorption of tryptophan at a DPPC bilayer membrane at 310.15 K (gel phase) [138], we observed a strong first coordination shell for TRP-water and TRP-DPPC pairs. In this study, we focussed in the liquid phase and only found relevant changes in the local structure and dynamics of TRP for cholesterol concentrations above 30%. TRP-DPPC binding involved coordination shells for the different oxygen sites of DPPC that are able to associate ('O2' and 'O8') versus the two tagged hydrogens ('H1' and 'H2') in TRP. Additionally, the distribution functions of TRP-CHOL revealed very stable hydrogen bonding. TRP is able to establish strong interactions with all solvating particles (water, DPPC, and CHOL), including a sort of double bridge between DPPC and cholesterol species. Typical HB distances have been found to be around 1.7-2.0 Å, which is in good agreement with experimental data [139]. Finally, the self diffusion coefficients of TRP are of the order of 10 −7 cm 2 /s, being strongly dependent of cholesterol's concentration.
In the case of melatonin, we have simulated its behavior when embedded in a cholesterol rich DMPC membrane at 303 K and 1 atm. Our interest was firstly focused on the local structure and angular distributions of MEL. In a similar fashion as in the case of TRP, strong hydrogen bonds between MEL-DMPC and MEL-CHOL have been found. The most important structures of MEL have been observed for two angular configurations: "folded" and "extended". Using a particular dihedral angle (Ψ), we observed two preferential values. The angle Ψ was revealed to potentially act as a reliable reaction coordinate, since two neat angular distributions around ∼81 • and ∼170 • were clearly distinguished. We also observed that introducing cholesterol into the system can favour MEL to exchange between extended and folded configurations. Again, the self diffusion coefficient of MEL was found to be of the order of 10 −7 cm 2 /s, although, in this case, with a very mild dependence on cholesterol's concentration.
Free energy barriers for serotonin, melatonin, and tryptophan at 323.15 K and 1 atm have been analysed using the reversible work theorem, which provides us a simple way to estimate the height of the barriers that are related to the interatomic distances. These features will be directly related to the formation and breaking of hydrogen-bonds. We have found marked first and second coordination shells that correspond to two minima of the PMF, with energy barriers for TRP-DPPC of the order of 10 kJ/mol. Most remarkable have been the binding between hydrogen 'H1' of TRP and oxygens 'O2' and 'O8' of DPPC. In the case of serotonin, we have found it to be a molecule strongly anchored at the membrane unlike to be solvated by water. Interestingly, melatonin has revealed to be able to interact both with water and DPPC, still showing moderately strong free energy barriers. In order to get more precise information, we have conducted well-tempered metadynamics simulations, and obtained 2D free energy landscapes for MEL binding to the DMPC-CHOL membranes. Two CVs have been considered: a dihedral angle Ψ and the distance z between the center of mass of MEL and the center of the lipid bilayer (set at z = 0). From our results, we have found that MEL is usually bound to the external side of the membrane, at distances z ∼1-2 nm and in two main configurations with Ψ = 70 • (folded) and 180 • (extended), with an energetic cost for the exchange between the two conformations of about 15-20 kJ/mol. After CHOL is introduced into the system, it pushes MEL to escape outside the interfacial region of the membrane and move away until it is fully solvated by the aqueous ionic solution. The energetic cost for MEL to leave the interface of the membrane towards the water bulk (z-distances around 4 nm) has been estimated at ∼10-25 kJ/mol. A very uncommon situation in our simulations was that of MEL accessing the center of the membrane, an energetically expensive process (free energy barriers of 40 kJ/mol). We believe that the findings that are presented in this work could be of practical use in two ways: (1) for the design of new reaction coordinates in similar systems of small molecules of biochemical interest, such as amino acids, neurotransmitters, drugs, or hormones and (2) from a more general perspective, to contribute the unveiling of the microscopic interactions of small molecules with cell membranes and the key role that is played by cholesterol in the properties of such molecules. All of this can lead to advances in the research of new pharmaceutical compounds and to a better understanding of the currently available ones.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: Here we can see how full convergence is obtained for the two relevant RDF, related to TRP-water and TRP-DPPC binding. After 50 ns (initial setup + full equilibration) and 100 ns (50 ns of production) the results are qualitatively similar but not fully converged yet, but when we extend to 120 ns, the results are virtually the same as those of 100 ns. Remarkably, the TRP-water structures quickly converged, because of the large amount of statistics (+5000 waters, stable hydration shells of TRP). This show that 100 ns is a reasonable length to achieve fully equilibrated results in the present case.

Appendix A.2. Convergence of WTM Simulations
Finally, to further evaluate the convergence of the metadynamics simulations, we reported the time cumulative average of 1D free energy profiles as defined in a previous work (see Formula S2 in Ref. [140]), i.e., averaging the two leaflets and projecting onto (integrating out) the alternative CV in a range larger than 1 microsecond in all cases. We have taken the case of 30% concentration of CHOL as an example. From the results of Figure A2, we can see that after long cumulative time lengths, the differences between a profile and the one immediately before are very small (up to 2.5 kJ/mol) and lead us to fully converged free energies for the two CV considered in the present study. Free energy [KJ/mol] t=500 ns t=800 ns t=1000 ns t=1100 ns t=1200 ns t=1300 ns t=1350 ns t=1400 ns