Molecular Dynamics Simulation of Small Molecules Interacting with Biological Membranes

Cell membranes protect and compartmentalise cells and their organelles. The semi-permeable nature of these membranes controls the exchange of solutes across their structure. Characterising the interaction of small molecules with biological membranes is critical to understanding of physiological processes, drug action and permeation, and many biotechnological applications. This review provides an overview of how molecular simulations are used to study the interaction of small molecules with biological membranes, with a particular focus on the interactions of water, organic compounds, drugs and short peptides with models of plasma cell membrane and stratum corneum lipid bilayers. This review will not delve on other types of membranes which might have different composition and arrangement, such as thylakoid or mitochondrial membranes. The application of unbiased molecular dynamics simulations and enhanced sampling methods such as umbrella sampling, metadynamics and replica exchange are described using key examples. This review demonstrates how state-of-the-art molecular simulations have been used successfully to describe the mechanism of binding and permeation of small molecules with biological membranes, as well as associated changes to the structure and dynamics of these membranes. The review concludes with an outlook on future directions in this field.


Introduction
Biological membranes are semi-permeable structures present in all living organisms. Membranes protect the cell from the external environment, compartmentalise cells and their organelles and control the selective transport of molecules in and out of these compartments. Biological membranes, and cell membranes in particular, are commonly described using the fluid mosaic model [1] where the membrane consists of a lipid bilayer into which membrane proteins are embedded (Figure 1). Biological membranes fulfil a wide range of functions and are highly complex structures [2]. A typical bacterial or mammalian membrane can consist of several hundreds of different types of lipids and include hundreds of different proteins [3,4]. The lipid bilayer of cell membranes is primarily made up of glycerophospholipids, sphingolipids and sterols, varying in proportion depending on the type of cell or the specific function of the membrane [4]. The lipid bilayer has two main functions. Firstly, lipids are critical for the function and structural stability of membrane proteins, which make up nearly a third of the human cell membrane proteome and account for approximately 60% of known drug targets [6][7][8]. Secondly, the bilayer controls the permeation of small molecules across the membrane and thus controls what enters the cell or cell organelles. In addition, the membrane mediates the lateral diffusion of small molecules bound to the membrane and can thus affect the mechanism of action of molecules that act as ligands for membrane proteins. It is thus not surprising that the interaction of small molecules with membranes affects the pharmacokinetics, bio-availability and mechanism of action of endogenous substances such as neurotransmitters as well as drugs. [9][10][11]. The permeation of small molecules is also critical for assessing the toxicity of exogenous particles such as diesel soot or silica dust [12], the use of nanoparticles for imaging, biosensing and therapeutic applications [13] and the use of organic molecules in the cryopreservation of plant germplasm [14,15]. Consequently, understanding small molecule-membrane interactions (SMMIs) and, in particular, the ability to predict the binding affinity and permeation coefficients of small molecules is an active area of research and aids our understanding of physiological processes and facilitates pharmaceutical development and a diverse range of biotechnological applications.
In this paper, we present an in-depth review on the use of molecular dynamics (MD) simulation methods to investigate SMMIs. In particular, we focus on studies aimed at understanding the interaction of water, small molecule drugs and drug-like compounds, short peptides as well as cryoprotective agents with model membranes mimicking the plasma cell membrane of mammalian or Gram-negative cells and the stratum corneum, the outermost layer of the epidermis. Validated force fields describing phospholipid bilayers with and without sterols and/or sphingolipids as well as force fields for ceramides and free fatty acids have been available for some time. As a result, there are numerus studies investigating the SMMIs with these types of membranes. In contrast, simulating bilayers containing lipopolysaccharide, a key component of the outer membrane of Gram-negative bacteria, was achieved for the first time in 2001 [16]. Subsequent studies have mostly focused on parameter optimisation and validation of lipid A and lipopolysaccharide, and including a variety of oligosaccharide and lipid combinations [17][18][19][20][21][22][23][24]. At the time of this review and to the best of our knowledge, only three simulation studies of SMMIs with bilayers containing either LPS or Lipid A/ lipopolysaccharide have been reported [25][26][27][28].
Further, this review does not cover the interaction of larger molecules with biological membranes, which have been discussed in a number of recent reviews of the membrane-binding properties of anti-microbial [29][30][31], amyloids [32,33] and cell-penetrating peptides [31,34,35]. Finally, challenges and issues that are common to all biomolecular simulations [36] such as the choice for force field [37], the assessment of convergence [38] or the validation of properties predicted from simulations [39] are not discussed in detail unless particularly relevant simulating SMMIs.
The remainder of this review is organised as follows. First, the fundamental concepts of model membranes, surface binding and permeation are introduced, followed by a brief summary of the experimental techniques most commonly used to study biological membranes and their interactions with small molecules. The main part of this paper reviews the use of unbiased MD simulation as well as various enhanced sampling methods to characterise the binding and permeation of small molecules through biological membranes. For each method, a brief definition of important concepts and the most relevant equations are given together with references to appropriate papers that provide a more detailed description of the underlying theory, and some of the advantages and limitations of each methods are discussed. While this review groups studies by the methods used, the main focus is on describing the molecular insight that can be gained from the simulations. The review concludes with a summary and an outlook on future directions in this field.

Model membrane systems
Due to the high complexity of biological membranes, their interactions with other molecules are commonly studied using model membranes. They are usually composed by one or more of the most abundant lipids found in the membrane of interest such that they reproduce the physico-chemical properties critical to the process one aims to study (e.g. fluidity or surface charge). Besides reducing the complexity, model membranes provide the ability to systematically investigate the effect of different environmental factors on SMMIs (e.g. lipid composition, pH, ionic strength, level of hydration or temperature). By controlling lipid composition and environmental factors, model membranes also enable a more direct comparison of the surface binding or permeation of a set of small molecules with varying physicochemical properties.
The most commonly used model membrane systems include micelles, liposomes, supported monolayers and bilayers (Figure 2) [40,41]. The choice of model system and its composition depends on the set of relevant physicochemical properties that are to be reproduced, the process that is under study (binding or permeation) and the experimental technique(s) used to monitor the progress of said process [2,[41][42][43]. The structures of some of the most common lipids used to mimic mammalian plasma cell membranes and the extracellular lipid environment of the stratum corneum of the skin are shown in Figure 3. For mammalian plasma membranes, the most commonly used lipids are neutral (zwitterionic) phospholipids like 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC), 1,2dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) or 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) [43]. Similarly, to mimic bacterial membranes the neutral lipid 1-Palmitoyl-2-oleoyl-snglycero-3-phosphoethanolamine (POPE) is used. To add surface charge, lipids like 1-palmitoyl-2oleoyl-sn-glycero-3-phospho-L-serine (POPS, mainly used in mammalian membranes) or 1palmitoyl-2-oleoyl-sn-glycero-3-phospho-(1'-rac-glycerol) (POPG, mainly used in bacterial membranes) [42,44] are added. Other lipids commonly added to mimic the plasma membrane include sphingolipids, phosphatidylinositol or cholesterol. Model membranes to mimic the stratum corneum are usually composed of the sphingosine-based ceramides, free fatty acids and cholesterol [45,46]. The thermophysical phase of the lipid bilayer used as a model system is also of importance. Lipid bilayers can be found in a variety of phases including sub-gel, gel, rippled and liquid crystalline (also called fluid phase). Phase transitions represent rearrangements in membrane structure leading to changes in the stability of bilayers (Figure 4). Each phase has a characteristic molecular arrangement described best with a set of structural properties like the area per lipid (APL), membrane thickness and lipid tail order parameter. For example, in the fluid phase lipid tails have high mobility and low lipid tail order. In contrast, in the gel phase lipid tail mobility is lower and is accompanied by an increase in lipid tail order. Phase transitions can be induced by changes in environmental conditions such as temperature, pH, ionic strength and hydration. In model membranes, transition from the fluid to the gel phase can be induced with either a reduction in temperature (thermotropic phase transition), an increase in pressure (barotropic phase transition), a reduction in hydration or a decrease in pH [42,47,48].
The interaction of small molecules with membranes can both induce phase transitions in model membranes as well as shift in the fluid-to-gel transition temperature (Tmelt). For example, in the context of cryobiology, a reduction in the Tmelt reflects stabilisation of the liquid crystalline phase of plasma membranes and the retention of their biological function at lower temperatures. A reduction in hydration upon desiccation leads to the opposite effect: an increase in the Tmelt and the stabilisation of the gel phase. During cryopreservation, where both desiccation and liquid-nitrogen temperatures lead to severe membrane damage, addition of sugars and other non-penetrating cryoprotective agents is aimed at stabilising the fluid phase of cell membranes and protect them from damage upon thawing [15]. Cell membranes exist mostly in the fluid phase, which is imperative for their integrity and function, making this the most relevant state in the study of SMMIs. The gel phase is, nonetheless, also found under physiological conditions. For example, dehydrated stacked lipid bilayers in the gel state form the extracellular environment the stratum corneum, which acts as barrier to prevent exogenous compounds entering the body and to prevent water loss [49,50]. As opposed to phospholipids, ceramides are found naturally in the gel state at physiological temperature. The study of the permeation of small molecules across these lipid layers in the gel state is thus of particular interest in the context of transdermal drug delivery [51][52][53][54].

Membrane binding and permeation
SMMIs can be classified into surface binding (or adsorption) and permeation across the membrane.
Surface binding describes the interactions of the small molecule at the water-lipid interface and the binding affinity associated with the molecule moving from a fully solvated state to a membrane-bound state. In contrast, permeation describes the flux of molecules from a solution environment on one side of the membrane, through the hydrophobic core to the other side of the membrane.
Both surface binding and permeation are initiated by the diffusion of the small molecule from the bulk solvent to the membrane surface. This process can be driven by electrostatic and/or hydrophobic interactions. The specific mechanism of interaction depends on the size, shape and physico-chemical properties of the small molecule, the composition of the membrane and environmental factors such as temperature, pH and ionic strengths. Seelig has described in detail the driving forces that govern the adsorption and insertion of peptides to membranes [55]. Much of this theoretical framework is also applicable to small molecules that do not penetrate into the bilayer (i.e. interactions that are restricted to surface binding). In this theoretical framework, the model used to describe adsorption and thus binding affinity depends on the chemical nature of the molecules. For neutral molecules, adsorption is described using 'hydrophobic partitioning', where the concentration of a molecule in the bulk and at the membrane is considered to be at equilibrium, and the binding affinity is obtained from a simple partitioning principle. For charged molecules 'hydrophobic partitioning with an electrostatic correction' is used, whereby the electrostatic attraction or repulsion of a molecule at the membrane surface increases or decreases its surface concentration with respect to the bulk. An alternative to partitioning models is to calculate the binding affinity of a molecule using adsorption isotherms such as the Langmuir isotherm.
Similar to surface adsorption, there are different conceptual and theoretical frameworks to describe permeation. In early studies of permeation the membrane was considered as a homogenous oil slab and permeation data was often correlated to water/oil partition coefficients (referred to as the Meyer-Overton rule) [56]. However, this model underestimates the effect of the interfacial region and boundary effects. Thus, in subsequent models the permeation process was conceptualised using the homogeneous solubility-diffusion model, in which three successive steps are considered: (i) binding of the molecule to the water-lipid interface and subsequent insertion into the lipid core of the membrane, (ii) passive diffusion through the membrane core, and (iii) unbinding of the molecule at the water-lipid interface on the other side of the membrane. While this model treats the water-lipid interface explicitly, the membrane core is still considered to be homogenous.
Instead, the bilayer interior shows distinct regions that exert different resistance to a permeant. This led to the development of the inhomogeneous solubility-diffusion model, which is based on model that divides the lipid bilayer into four distinctive regions that describe a single leaflet in the lipid bilayer: (1) the low head group density, (2) the high head group density, (3) the high tail density, and (4) the low tail density region [57]. Each region poses a different resistance to the permeant ( Figure 5). Figure 5. The four-region model for a phospholipid bilayer in full hydration. Numbers correspond to the regions described in the main text: (1) low head group density, (2) high head group density, (3) high tail density, and (4) low tail density.
As for surface binding, permeation across a lipid bilayer is influenced by the physicochemical properties of the permeant molecule such as its molecular weight, shape, hydrophobicity/lipophilicity and net charge [58]. For molecules with molecular weight smaller than 50 it has been shown that permeation is mostly dependent on the size of the permeant molecule.
For larger molecules, permeation is more likely a function of volume and shape (if corrected for hydrophobicity) [59]. Hydrophobicity/lipophilicity and net charge affect permeation because they dictate the energy barrier experienced by the permeant when moving across the different areas of the membrane. If the molecule is too hydrophobic, it will become trapped in the lipid core reducing its overall permeability [60]. Similarly, permeability is reduced if the molecule bears hydrophilic or charged functional groups as the preferential interactions with the polar head groups at the waterlipid interface will result in a large energy barrier to penetration beyond the interfacial region [61].
Consequently, there is great interest in developing permeability models that can predict permeation by considering the physicochemical properties of the solute.

Experimental characterisation of SMMIs
Various experimental biophysical techniques are used for characterising SMMIs using model cell membranes. These methods provide information about the changes in structure, dynamics and stability of membranes upon their interaction with small molecules, as well the orientation and location. In addition, these techniques are used to investigate the thermodynamics of SMMIs. In the following, a brief overview of experimental techniques to study SMMIs is provided. For details of a particular method the reader is referred to the references cited.
Nuclear magnetic resonance (NMR) spectroscopy is one of the most common methods to characterise the structural and dynamical properties of lipid bilayers [5,62]. As illustrated in Figure   6, solid-state NMR (ssNMR) can be used to study the dynamics of lipids in the fluid-phase over a wide range of time scales, providing information about the conformation and orientation of individual lipids and their diffusion, as well as collective motions such as membrane deformations.
By comparing these data in the absence and presence of a membrane-active molecule, ssNMR can be used to describe the effect of SSMIs on lipids or membrane structure. NMR can also provide information to construct phase transition diagrams of model cell membranes [63], which can subsequently be used to study the effect of SMMIs on phase transition temperatures [64]. In addition, NMR can provide structural insight into the conformation of small molecules when bound at the membrane interface [65][66][67][68]. For example, NMR has been used to determine the preferred position and orientation of drug-like molecules, water, ethanol and flavonoids in lipid bilayers [67,69,70] or to determine the high-resolution 3D structure of small molecule-lipid complexes ( Figure 6) [68]. Other spectroscopy techniques used to study membranes and SMMIs include Föster resonance energy transfer (FRET), electron paramagnetic resonance (EPR), fluorescence correlation spectroscopy (FCS) and Fourier transform infrared (FTIR) spectroscopy. FRET is particularly useful to monitor the spatial organisation and distribution of lipids including the formation of micro domains in membranes [71], and to detect the intercalation of small molecules into lipid bilayers. Whilst FCS and EPR are less commonly used to directly study SSMIs, both techniques provide information on lipid dynamics and can thus be used to study effect of SMMIs on the membrane. For example, FCS can be used to characterise diffusion processes in membranes by 'tracking' individual lipids both in 'simple' model membranes as well as complex lipid mixtures or membranes with high heterogeneity [72]. EPR is used to determine the mobility and order of lipids and thus provides a measure of membrane fluidity. EPR can also be used to provide information on the water accessibly as a function of membrane depth [73,74]. FRET and FCS require fluorescence labelling while EPR relies on paramagnetic labels. In contrast, FTIR does not require any labels. The hydrogen bonding pattern in the lipid head groups as well as the ordering of the acyl chains gives rise to characteristic IR absorption bands that are sensitive to lipid packing. FTIR is thus well suited to study both the overall structure and organisation of membranes as well as the interaction of small molecules with specific parts of the lipid molecule [74]. For example, a number of studies have used FTIR to detect the interaction of steroid hormones [74] or small-molecule drugs [75] with lipid bilayers and liposomes.
In addition to spectroscopy techniques, neutron and x-ray scattering methods are used to probe the structural properties of model lipid bilayers. These methods provide dispersion patterns that can differentiate the lateral packing arrangement of lipids characteristic to a specific phase. In addition, data from scattering experiments can be used to determine bilayer thickness [76,77]. As illustrated in Figure 7, neutron scattering density profiles from lipid bilayers (in the absence and presence of small molecules) can provide information on how SMMIs affect lipid phases and water distribution in membranes. Similarly, scattering and diffraction techniques have been employed to determine the position of sugars, alcohols and drug molecules in membranes, the effect of SMMIs on the structure and bending modulus of membranes, the lipid tilt angle or lipid diffusion [44,[78][79][80][81][82][83][84].
Scattering techniques are also useful to study the effect of different levels of hydration (i.e. the water: lipid ratio) on membrane structure and SMMIs [85]. The aforementioned experimental techniques are often combined to gain understanding of the molecular mechanism of action of drugs and peptides upon interaction with lipid membranes [74,75,86,87] .
Besides the structure and dynamics of lipid bilayers, characterising of the thermodynamics and kinetics of SMMIs are critical to understand surface binding or permeation. These have been studied using a range of techniques including differential scanning calorimetry, isothermal titration calorimetry (ITC) and surface plasmon resonance (SPR). Differential scanning calorimetry is commonly used to studies to determine their phase transition temperatures and associated activation energies of model membranes including any changes to these properties upon interaction with small molecules [44,47,74,75]. ITC can be used to measure the amount of thermal energy absorbed or released as a result of membrane binding, thus allowing the determination of binding constants (free energy of binding), enthalpies and stoichiometry of SMMIs [75,88,89]. SPR is used to determine the kinetics of binding and unbinding of small molecules to lipid bilayers, yielding corresponding association and dissociation rate constants, as well as the overall binding affinity constant [90][91][92]. As the method is based on a change in mass upon membrane binding, monitoring the association and dissociation of molecules with low molecular weight can be challenging.
However, improvements in sensitivity have enhanced the use of SPR for SMMIs [93].   The methods mentioned above provide information on the effect of SMMIs on membrane structure and dynamics as well as the thermodynamics and kinetics for the binding of a small molecule to the membrane surface. Different methods are required to study the permeation of small molecules.
The parallel artificial membrane permeation assay (PAMPA) measures the flux of a molecule from one donor aqueous compartment to the next through a single lipid bilayer barrier. The method is based on the detection of changes in concentration of the molecule in both compartments, using techniques such as UV-visible spectroscopy or mass spectrometry. PAMPA is especially useful for measuring the permeability of drug-like compounds across human skin model membranes [95,96] and the human blood brain barrier [97,98]. Despite PAMPA being routinely used in the pharmaceutical industry for the automated, high-throughput measurement of permeation coefficients of drug candidates, this approach does not provide information about the underlying biophysical mechanism of permeation [99].
3 Molecular simulation approaches to the study of SMMIs MD simulations are a computational method that can complement experiments by providing both atomistic level detail and high temporal resolution. As illustrated by the many studies discussed in this review, MD simulations have been used extensively to study SMMIs, both for rationalising experimental data and for predicting changes to membrane properties such as APL, membrane thickness, order parameters and lateral diffusion coefficient. Simulations have also provided insight into the effect of lipid composition on these properties.

Common quantities in molecular dynamics simulations
By linking ensemble averages of molecular properties to bulk thermodynamic properties, MD simulations can also be used to calculate binding constants, permeation, diffusion and partition coefficients. However, obtaining these properties from MD simulations relies on sampling representative molecular configurations to calculate free energy differences. The free energy F of a system is related to its partition function Z by: Here β = 1/kBT, kB is the Boltzmann constant, T is the absolute temperature, and E(r) is the energy of the system in configuration r and N is the number of degrees of freedom of the system.
In simulations of SMMI, one is commonly interested in the difference in the free energy between an initial and a final state. For surface binding, this usually corresponds to the free energy difference between the small molecule in solution (unbound state) and the small molecule bound to the waterlipid interface. Since the location of energy barriers and minima in SMMIs are usually not known a priori, the process of going from the unbound to the bound state has to be exhaustively described along every point of its trajectory. This can be done by using a reaction coordinate (RC), ξ. In SMMIs, the RC is usually the centre-of-mass (COM) distance between the membrane and the interacting molecule. The probability distribution along ξ integrated over every degree of freedom except ξ is given as [100] ( Substituting equation (3) in (2) expressed as a function of the reaction coordinate ξ enables the calculation of F(ξ), which is also commonly referred to as the potential of mean force (PMF). Once the free energy is known, the binding constant (and thus binding affinity) can be obtained upon integration of the PMF over ξ using the relation As described in previous sections, SMMI depend on lipid composition and environmental factors.
Consequently, changes to these can affect the PMF and the values of the free energy and the properties calculated from it.
In contrast to the binding constant for surface binding, the permeation coefficient P of a molecule is related not only to the ξ-dependent free energy (of permeation), but also to the local diffusion constant of the molecule in the membrane, D(ξ), as defined by Where ri(t) is the time dependent position of the molecule and ri(0) is the position at time t = 0.
The Green-Kubo relations use the autocorrelation function (ACF) of the force, velocity or position of the diffusing particle. The underlying assumption is that the diffusing particle undergoes a random Hummer's equation was derived to directly use data from umbrella sampling simulations that employ a harmonic restraining potential on the COM of the diffusing particle [103] (See umbrella sampling section below for more details). The Bayesian inference method [103] is more generic in its derivation and requires fewer assumptions, which makes it compatible with different biasing schemes, including time-dependent biasing.
A few studies have tested the consistency of these approaches. Lee  These studies also aimed to assess the effect of solute binding on the integrity of the membrane upon interaction in terms of order parameters, bilayer thickness or APL. Based on these information, the molecular mechanism of action of solutes ranging from disinfectants and drug molecules to cryoprotective agent and environmental pollutants. For example, Mukhopadhyay et al. [106] conducted a simulation of 15 molecules of pentachlorophenol interacting with POPC and POPE lipid bilayers. In both types of bilayers, pentachlorophenol was found to interact at the interface with its aromatic ring placed almost parallel to the axis of the lipid tails, preferentially below the carbonyl of the head groups. In this position, the pentachlorophenol molecules form hydrogen bonds with both interfacial water molecules and the oxygen atoms in the lipid head groups. Pentachlorophenol marginally affected the structure of the bilayer showing a small increment in tail order parameters and bilayer thickness. This study helped rationalise the observed enthalpy-driven partitioning of nonpolar molecules into bilayers. Similarly, a study by Markiewvicz et al. [78] showed that the binding of the non-steroidal anti-inflammatory drugs ketoprofen, aspirin and piroxicam (featuring aromatic rings as well as polar groups with hydrogen bonding capabilities) to a POPC/cholesterol bilayer had little to no effect on the bilayer structure. The negligible effect of these drug molecules on the order parameters and APL has been predicted by simulations as well as measured by X-ray diffraction. The difference between the drugs was mainly in their ability form hydrogen bonds with interfacial water, surface bound Na + ions and the polar lipid headgroups resulting in changes in the water-lipid interface.
Unbiased MD simulations also provided the first molecular picture of how cryoprotective agents interact with membranes and demonstrated the role of lipid thinning in the mechanism of action.
Simulations successfully reproduced the lateral expansion measured as an increase in the APL of DPPC upon interaction with aqueous solutions of short chain aliphatic alcohols (ethanol and propanol) in a concentration-dependent manner [107]. Longer aliphatic alcohols (8 to 18 carbons) also showed a concentration-dependent increase in APL and a decrease in bilayer thickness [108].
Similarly, results from Kondela et al. [109] showed that the insertion of alcohol molecules into the lipid bilayer created more space in the head groups region, allowing water to find new interaction sites in that region. In particular, intercalation of methanol and ethylene glycol in the upper region of the hydrophobic core of the lipid bilayer caused the lateral expansion (i.e. increase in APL) and increased lipid tail disorder, in what is known as the 'alcohol-like' mechanism [110]. MD simulation estimates allowed ranking polyols on their effects on DPPC bilayers. Compared to glycerol and ethylene glycol, propylene glycol induced the strongest decrease of tail order parameter [111].
DMSO is a solvent commonly used as a potent cryoprotective agent (alone or in combination with alcohols) showing marked effects on the structure of lipid bilayers with similar trends as alcohols.
Lin et al. [112] simulated a DMPC bilayer in the presence of DMSO at low concentration. The insertion of DMSO molecules beneath the lipid head groups caused thinning of the lipid bilayer and increased APL as expected. These structural changes resulted in increased permeation of water across the bilayer, which is sign of impaired function and can be deleterious for a cell. In addition, DMSO induced a change in the orientation of the lipid head group (the P-N dipole vector), thus changing the charge-density distribution of the head group region. DMSO induces very drastic structural changes above critical concentrations, including pore formation and, ultimately membrane destruction (Figure 9) [113]. Similar effects required higher concentration of alcohols as shown in a comparative study of the effects of individual cryoprotective agents on the structure and integrity of DPPC bilayer [114]. Whilst many cryoprotective agents cause marked structural changes on bilayer structure, unbiased MD simulation underlines the stabilising effect of sugars such as trehalose, maltose, sucrose or glucose on phospholipid bilayers such as DPPC, DPPE and POPC [115][116][117][118][119]. The APL and lipid order parameters were found not to be affected by the presence of the sugar molecules, suggesting that the overall structure of the bilayer remains largely undisturbed. It was further observed that the sugar molecules accumulate at the membrane surface where they interact with the head groups, resulting in a significant degree of water replacement. Whilst not affecting the structure, simulations of a DOPC bilayer, however, showed that sucrose and to a lesser extent trehalose decrease the lateral diffusion of phospholipids. This was confirmed experimentally by fluorescence correlation spectroscopy [120]. Compared to trehalose, sucrose was seen to interact simultaneously with more head groups , resulting in a larger decrease in lateral diffusion [120]. This suggests that sugars may act by controlling the mobility of head groups, preventing the thermal softening of their interactions and the resulting increased dynamics of the lipid tails. Overall these studies suggest that disaccharide molecules are able to intercalate within the phospholipid head groups of the lipid membranes, simultaneously forming a stabilising network of specific hydrogen bonds [115][116][117][118][119].
This stabilising effect can also prevent thermal disruption, as shown by simulations of sugar molecules on DPPC lipid bilayers. It was shown that at temperatures high enough to cause thermal disruption, structural damage to the bilayer is prevented to a significant extent by the interaction of sugar molecules with the phospholipid head groups [119]. Recent MD simulations of the interaction of sucrose with DOPC bilayers at different levels of hydration have shown that sugar molecules interact preferentially with the lipid bilayer interface as described above at low concentration and/or high hydration, but at high concentration or low hydration they preferentially accumulate in the inter-bilayer solution as the bilayer interface becomes saturated with sugar molecules [121].
All of the studies described above investigated the interaction of molecules with membranes in the fluid phase. In addition, unbiased MD simulations were also used to study SMMI for membranes in the gel phase and the effect of SMMI on phase transition and separation of ordered and disordered phases. A biologically relevant example of a membrane in the gel-phase is the stratum corneum (SC). and happens spontaneously to bilayers composed by lipids featuring asymmetric tails (as will be discussed further below). In palmitate bilayers interdigitation is of particular interest because the lipid order resembles the gel phase (in which lipid tails are highly ordered and tightly packed), with the difference that the lipid tails from opposite leaflets of come in full lateral contact so that it is visually hard to distinguish one leaflet from the other [124]. It was observed that methanol induces the G-to-F transition at high temperatures, and the F-to-G transition at low methanol concentrations and low temperature.
Interestingly, the F-to-I transition was induced at high concentrations and low temperatures. This clearly demonstrated that at high concentrations of alcohol, the membrane actually favours the I but not the G phase. However, in the earlier experiment that originally investigated the biphasic effect, the G and I phases could not be differentiated. The simulations thus showed that the process of interaction is better described by a triphasic rather than a biphasic effect. The interdigitated phase is not exclusive to palmitate bilayers; in fact, earlier works had already observed the concentration-dependent increased interdigitated phase for a DPPC bilayer exposed to increasing concentrations of ethanol [126].
Apart from altering phase transitions, SMMI can influence domain formation lipid (de)mixing. A study combining coarse grained MD simulation and confocal fluorescence microscopy investigated the effect of saccharides on mixed component membranes. In a membrane composed of two phases, a liquid ordered (Lo) DPPC rich phase and a liquid disordered (Ld) DPPC rich phase, the addition of non-reducing disaccharides such as trehalose and sucrose was observed to disrupt the separation of these domains and promote lipid remixing, while the addition of monosaccharides (glucose and fructose) or reducing disaccharides (maltose, palatinose and gentiobiose) did not [127].
This difference was attributed to the interaction of bulky disaccharides with local membrane surface defects, which are more abundant in the Ld phase. These interactions are believed to lower the free energy of the mixed system in the Ld phase compared to that of the Lo phase, providing a driving force for the mixing of lipids [127].
Taking into consideration the evidence towards distinctive regions, the seminal work by Marrink and Berendsen [57,59] was one of the first studies to use MD simulation to calculate the permeability of small molecules across a lipid bilayer. They reported a series of short simulations where a small molecule (water, molecular oxygen or ammonia) was inserted and constrained at different positions along the bilayer normal, referred to as the z-constraint algorithm. The free energy profile was later rendered from these constrained simulations. Their theoretical framework for the inhomogeneous solubility-diffusion model of permeation was used to calculate permeation from said free energy profile and the position dependent diffusion coefficient. The widespread fourregion model described in Section 1.2 was derived in their work, and together with their theoretical framework lies at the core of most MD permeation studies. For example, Bemporad et al. [131,132] used this approach to predict the permeation across a DPPC bilayer of a number of small organic molecules: acetamide, acetic acid, benzene, ethane, methanol, methyl acetate and methylamine.
As shown in Figure 11, marked differences in the free energy profiles of permeation were observed, which relate to the chemical nature of the molecules. Hydrophilic compounds exhibit larger free energy barriers in the low density tail region, whereas for hydrophobic compounds the barriers are located closer to the high density head group region (regions 4 and 2 in Figure 5, respectively). It was also observed that the size of the permeant is of lesser importance for diffusion within the bilayer than for diffusion in water. However, the predicted permeability does indeed change with molecule size and this arises from the partition (insertion) into the bilayer. This revealed that permeation coefficient estimates depend to a higher extent on partition (insertion) into the bilayer than on the actual diffusion through the bilayer lipid core, suggesting that partition is the rate limiting step of permeation. Findings regarding the behaviour of larger, more complex molecules such as the adrenoceptor antagonists alprenolol, atenolol and pindolol (featuring both hydrophobic and hydrophilic portions) upon insertion into a DPPC bilayer support that suggestion [133]. Substantial increases in computer power allow now the simulation of larger systems on times scales where the spontaneous permeation of water can easily be observed, thus avoiding the need for imposing constraints. Nevertheless, simulation timescales substantially increase for describing the permeation of larger, bioactive molecules, rendering unbiased MD simulations unfeasible [134,135].
In addition, simulating more realistic lipid composition (with a variety of lipids) requires larger simulation systems thus significantly increasing the computational power required for simulating such timescales. Furthermore, the increased heterogeneity inside bilayer composed of multiple types of lipids requires considerations beyond a simple one-dimensional permeation path. These reasons highlight the use of sophisticated enhanced sampling methods to describe such complex processes using MD simulations.
Enhanced sampling methods to study small molecule -membrane interactions The need for enhanced sampling methods is not unique to simulations of SMMI but a challenge Multiple enhanced sampling methods have been developed over the last few decades to address this sampling problem [36,[136][137][138][139]. In the context of SMMIs these methods can be divided into two groups: i) those that represent SMMIs using a selected subset of representative degrees of freedom to enhance sampling along them, and ii) those which enhance the exploration of phase space along all possible degrees of freedom. The first group of methods are efficient at enhancing sampling once information is available about the degrees of freedom that contribute the most to the free energy of interaction and permeation. Examples of these methods include umbrella sampling and metadynamics. The second group has the advantage that no prior knowledge of this kind is needed, but these methods can be considerably more computationally expensive. Examples of this second group of methods includes parallel tempering and replica-exchange MD. In the following we discuss the most commonly used and widely validated methods of both approaches.
Umbrella sampling Umbrella sampling (US) is one of the most common enhanced sampling methods used in MD simulation and is indeed the most extensively used to characterise SMMIs and obtain related free energy profiles. The method requires a RC that represents the progression of the SMMI, i.e. the transition from state A to B (e.g. unbound to membrane-bound). In US simulations sampling along the RC is enhanced by applying an additional energy term (the biasing potential) that restrains the location of the small molecule at a given point. In some cases, sampling is enhanced along more than one RC, referred to as multi-dimensional US [140]. The typical approach to US consists of setting up a number of simulations (referred to as 'windows'). In each window the biasing potential is centred at a discrete value along the RC (Figure 12). The result is a series of N partially overlapping histograms, where N is the number of windows. Each histogram provides a biased probability distribution along the RC. An alternative is to keep the biasing potential constant and slowly vary the value of the RC, in a protocol referred to as steered MD or pulling simulations [100]. To obtain the free energy, Ai(ξ), where ξ is the reaction coordinate, the unbiased probability distribution of the system, ( ) , is calculated as [100] ( ) = Where E(r) is the energy of the system with coordinates r, β=1/kT and δ[ξ'(r)-ξ] is the number of states in the simulation that have a value of ξ'. Assuming ergodicity in the system, the biased probability distribution, ( ) , obtained from the simulation will be equal to With ωi(ξ'(r)) being the biasing potential applied on the system. Since integration of the numerator is performed over all degrees of freedom except for ξ and the bias depends only on ξ, using Eq. 8 and Eq. 9 we can derive A more thorough derivation of Eq. 10 can be found in the in-depth review of US by Kästner [100].
This derivation is analytical and only assumes that enough sampling is provided for every value of ξ.
As ( ) is obtained from the US simulation and the biasing potential is given analytically, the free energy, Ai(ξ) can be derived as Here Fi is a constant that can be estimated by approaches such as the weighted histogram analysis method (WHAM) [141] when combining results from different windows. However, WHAM requires the windows to have sufficient overlap in their distributions [142]. An alternative to WHAM is umbrella integration, which instead relies on the calculation of the derivative dln ( ) /dξ and does not strictly require an overlap in the energy distributions [143].
The most common RC used in simulations of SMMIs is the COM distance between the molecule of interest and the membrane, as illustrated in Figure 12. This focuses the sampling effort to only one RC and thus simplifies the N-dimensional landscape to one dimension. However, use of this RC will fail to capture the effects of relevant degrees of freedom (DOFs) if these have slower relaxation times than the time dedicated to sample each window. The lack of sampling of these DOFs could affect the accuracy of the free energy calculation, especially if they are responsible for configurational transitions of the small molecule or the reorganisation of the lipids around it [144].
These slow-varying variables are usually referred to as 'hidden variables' that govern the exploration of the configurational landscape [136]. In any case, if the configurational ensemble is incomplete, i.e. it is insufficiently sampled, then the calculation of the free energy will be inaccurate.
One drawback of using the COM distance as the only RC is that it assumes that the COM corresponds to a fixed point in the 3D space of the simulation box. In the case of bilayers it is possible for the COM to move due to bilayer undulation. For this reason, Nitschke et al. concluded that for an accurate calculation of a RC using the COM distance it is advisable to consider only a cylindrical region of the membrane centred on the solute under study [145], rather than the entire membrane.
A recently developed alternative CV called 'distance from contour' can be defined in terms of the distance to both the closest and furthest edges of the membrane, which accounts for local fluctuations in the shape of the membrane due to the SMMI [146]. Additionally, there is a tendency to assume a symmetric behaviour of the SMMI along the COM distance across the entire membrane, and thus many studies only calculate the free energy for windows from solution to the COM of the lipid bilayer (i.e. across one leaflet). Lee et al. [99] conducted a study of the permeation of small hydrophobic molecules (urea, benzoic acid and codeine) across a DMPC bilayer in the fluid phase and showed the occurrence of asymmetric free energy profiles. This is likely due to insufficient sampling along the RC and stresses the importance of extending the windows to cover the entire translocation process (i.e. across both leaflets).
Other authors have used additional descriptors of the process of permeation in conjunction with the COM distance. This is especially relevant when the molecules under study are large in size and more complex in structure, such that other DOF are needed to fully describe the SMMI. These other DOF usually include geometric parameters such as the orientation/rotation of the permeant and chemical parameters such as dihedral angles and formation internal H-bond [104,147]. As more of these molecular descriptors are included in the analysis of SMMI the more reliable the estimates and the more thorough the insight into the mechanism of interaction. However, it also means that more dimensions need to be explored, which translates into additional use of computational resources.

MD simulations of small molecule-membrane interactions using umbrella sampling and related methods
There are many interesting studies using US for the characterisation of SMMIs in the literature as it is a well-established method and is implemented in most MD simulation packages. In this section we present studies selected to highlight the insight gained and their agreement with experiments.
We also comment on the known drawbacks of the method and optional ways to overcome them.
Similar to the previous section on unbiased MD, estimations obtained from US can be validated by comparison to experimental structural information or thermodynamics data. For example, the interaction of tricyclo-decane derivatives with a POPC bilayer was simulated for 15 ns at several COM distances. The preferential location of tricyclo-decane derivatives in the interfacial region and their orientation (in terms of tilt angle with respect to the normal vector of the membrane) were in agreement with NMR studies [148]. And an example of using thermodynamic properties to validate and compare the SMMI of molecules with varying physico-chemical properties is the study by Tieleman and collaborators [149]. their side chains remained hydrated as they penetrated into the acyl chains. In these cases the free energy of partitioning was dominated by the energetics of these water defects and no longer reflected a simple partitioning between the aqueous interface and the lipid region of the bilayer.
Based on the comparison of the PMFs for the charged and uncharged side chains, it was predicted that Lys, Glu and Asp become uncharged before entering the lipid bilayer core, whilst Arg can either be charged or uncharged. Based on these insights authors stressed that lipid bilayers should not be treated as a generic low-dielectric environment slab.
Monitoring properties related to degrees of freedom of the solute provides deeper insight of the mechanism of interaction. One example is the simulation of the permeation of the small polar drug piracetam across a DOPC bilayer conducted by Ribeiro and co-workers [147]. They described the internal flexibility and orientation of piracetam by measuring several angles. Their simulations show that the drug experiences a decrease in internal flexibility and rotational freedom when entering the head group region of the lipids before adopting an energetically favorable combination of internal dihedral angles at the centre of the bilayer. Interestingly, symmetrical directionality in the interaction of the drug with the head groups was observed, such that piracetam entered in a particular orientation and then flipped to leave the bilayer in the opposite orientation. Interestingly, their findings rationalise the empirical rule that the fewer rotatable bonds in a drug molecule, the higher its bioavailability (Veber's rules) [150]. This is because their simulations showed that if a drug molecule has fewer rotatable bonds, it will lose less entropy when it enters the restrictive membrane environment, which favours its permeation.
Similarly, Lee et al. [104] conducted US simulations of the permeation of three capped aromatic amino acids across a DOPC bilayer and monitored their freedom to rotate and to change conformations at selected depths. The molecules studied are N-acetylated forms of phenylalanine (Ac-Phe-NH₂, or NAFA), tryptophan (Ac-Trp-NH₂, or NATA), and tyrosine (Ac-Tyr-NH₂, or NAYA), which can be treated as dipeptides. In the surface of the bilayer dipeptide-membrane interactions are similar for the three molecules: the preferred orientation is with their side chains towards the membrane interior and their backbone towards the solution. However, after insertion into the head group region and lipid core the conformational space available to the dipeptides is severely reduced.
A lower and broader minimum for NAFA (the most hydrophobic) in the PMF shape can be rationalised in terms of increased sidechain conformational freedom compared to NATA and NAYA, as illustrated in Figure 13. The PMF for NAFA also explains its faster translocation times (in experiments conducted at physiological conditions) as it displays a lower free energy barrier at the bilayer centre. Noteworthy, the rotational freedom of all dipeptides was regained at the bilayer centre. A number of studies showed that the protonation state of a molecule can affect the free energy of interaction and therefore affect the rate of permeation, as experimentally shown for NAFA by Lee et al. [104]. To study this effect Ma et al. [69] analysed the PMF of permeation of the tricyclic drugs amitriptyline and clozapine in different protonation states across different lipid bilayers (DPPC, POPC and DMPC). Their results described the spontaneous insertion and preferred location of these molecules. The charged forms of the drugs showed preferential binding at the water-lipid interface, with a higher permeation free energy barrier compared to that of their neutral forms. Another study comparing the permeation of the three physiologically relevant forms of warfarin (one ionised and two neutral forms) across a DOPC bilayer showed that the ionised form was the slowest in penetrating into the head group region [151]. Additionally, Boggara and Krishnamoorti [152] calculated the PMF and permeation coefficients of the neutral and ionised forms of ibuprofen and aspirin across a DPPC bilayer. Their location in the bilayer, corresponding to the energy minimum in the PMF, was affected by both the structure and charged state of the drug molecule. The charged molecules were transported across the membrane with their polar groups fully hydrated, resulting in deformation of the bilayer. In contrast, the neutral forms permeate dehydrated, resulting in little perturbation of the bilayer. Although these prominent deformations provide a plausible molecular scenario that explains differences in permeation, they could arguably be the result of an inability of the force field to reproduce charge redistributions occurring in the permeant molecule upon a change in the dielectric environment. The inclusion of polarisation effects in the study of SMMI brings additional computational costs and is an area still in development [37]. Nonetheless it has been proven to significantly change the estimations of PMF towards more accurate values [153,154]. The increased free energy barrier to permeation helped to rationalise the higher inwards flux of ROS reported for plasma membranes of cancer cells, given their reduced cholesterol content [157]. In contrast, presence of 10% cholesterol lowered the free energy of partitioning of four drugs doxorubicin, paclitaxel, nicardipine, and morphine across a POPC bilayer. The extent of the effect is dependent on the nature of the substrate: For morphine the barrier was lowered by ~30 kJ/mol whilst for paclitaxel the reduction was <5 kJ/mol [158].
Changes in membrane structure and dynamics become even more apparent membranes in gel phase. Permeation coefficients across lipid bilayers in the gel state are always significantly lower in comparison to those in the fluid phase. This is a consequence of the higher density and molecular order in the lipid core of the bilayer, which arises from the decrease in the APL and the accompanying increase in bilayer thickness. The lower permeation coefficients across lipid membranes in the gel state highlights the need for enhanced sampling methods to better describe SMMI with slower dynamics. Paloncyova et al. [159] compared the permeation coefficients of three small molecules across a bilayer made of pure ceramide 2 (CER2, also known as CER NS) in the gel phase and a DOPC bilayer in the fluid phase. The chosen permeants cover a range of lipophilicities: p-aminobenzoic acid and two esters of it that differ in tail length (ethyl and butyl). Comparison of the predicted PMFs showed higher energy barriers for p-aminobenzoic acid in both bilayers in comparison to its esters, but they were much higher in the gel-phase bilayer than in the fluid phase ( Figure 14). This is likely due to the higher density and lipid tail order of the CER2 bilayer that the permeants need to disrupt. The trends in the free energy barriers in the PMFs were found to be conserved between the fluid-and gel-phase bilayers and the relative ordering of free energy differences for the three molecules is in qualitative agreement with experimental permeation data. As expected, the more hydrophobic alkyl esters of p-aminobenzoic acid are more permeable, exhibiting lower insertion free energies. This was consistent with the deeper insertion of the esters into the bilayer observed in a separate unbiased MD simulation of the same three molecules.
Simulations of the permeation of small molecules across a bilayer in the gel phase is particularly relevant to studying the SC, where the lipid gel phase is necessary for the skin to act as a lipid barrier.

MD simulations of single-component ceramide bilayers showed that the high lipid tail order is what
gives rise to the gel-phase and barrier properties of the SC [160][161][162][163][164]. Indeed, simulations of lipid bilayers containing mixtures of free fatty acids and cholesterol without ceramide [165] cannot reproduce the distinct properties of the SC. Only when simulating ternary mixtures containing ceramide (CER), cholesterol (CHOL) and free fatty acids (FFAs) to represent the most abundant lipid species in SC it is possible to better reproduce the properties of the SC [166][167][168][169][170][171][172]. It has also been shown that to reproduce the properties of the SC the asymmetry in the tail length of ceramides is important [173,174], which has also been shown experimentally [175,176]. A shift towards shorter alkyl tail lengths in CER2 results in an 'imbalanced' SC as found in skin conditions such as psoriasis and dermatitis, causing increased permeability (i.e. the SC loses its barrier properties) [177]. The presence of shorter alkyl tails have also been associated with a decrease in the barrier properties of pure CER2 bilayers due to the increased free space inside the lipid core and a reduction in membrane thickness, which leads to an increase in the predicted permeability coefficient of water by one order of magnitude [161]. Longer alkyl tails are pivotal for the interdigitating region in the SC because in the absence of different tail lengths interdigitation between lipids from opposing leaflets cannot occur.
A number of studies have used either US or the z-constrained algorithm to calculate permeability coefficients of water and small molecules across the lipid models of the SC [53,[178][179][180]. Das et al. [178] employed the z-constraint algorithm to obtain the PMF for the permeation of water across a SC model consisting of ceramides, cholesterol and free fatty acids in a 2:2:1 molar ratio. The free energy barrier for permeation for water was predicted to be located in the densest and most ordered region of the lipid tails just below the head groups. The predicted permeation coefficients were nonetheless significantly lower than those measured experimentally using skin samples. In a similar study, Del Regno and Notman [179] calculated the permeation of water across the SC   (Figure 16). Other non-hairpin conformations were also observed in a model of the SC that better represents the alkyl tail polydispersity of the three main CER subclasses and FFAs [181]. Very recently, a new SC model [172] was developed to account for the experimentally determined low hydration of the stacked lamellar structure of the SC [49,182]. The model was further optimised to match the regular patterns of stacked bilayers observed in cryoelectron microscopy studies of human SC. The most radical feature of this model is that CER molecules are modelled in a splayed conformation, with a ~180° angle between the sphingosine and alkyl lipid tails (Figure 17). This new structural model of the SC has been used to predict the permeation coefficients of small polar molecules (water, ethanol and DMSO) and alicyclic molecules (benzene, codeine, naproxen, nicotine and testosterone) [54]. The PMFs and diffusivity profiles were computed using a non-equilibrium free energy method in which the permeant molecules are slowly "grown" into the bilayer at various depths and then 'pulled' in both forward and reverse directions along a RC that is normal to the bilayer using a large umbrella potential. In Figure 17 (bottom panel) the PMFs for water and testosterone are shown as representative profiles for the polar and alicyclic compounds, respectively, illustrating the difference in the location of the free energy barriers and minima. As expected, all calculated permeation coefficients (except for ethanol) were found to be lower than the values determined experimentally using skin samples (in vitro) and ex vivo in diffusion cells where the SC is fully hydrated. The authors suggested that prolonged contact of these skin samples with water results in an enhancement of permeability, explaining the discrepancy between simulation and experiment. Apart from the ability of this membrane model to describe the structure and barrier properties of the SC, there are a number of other issues that can affect the reliable and accurate prediction of permeation coefficients using US and related methods. These include insufficient sampling of the orientation of the molecule in the membrane and the bias introduced by the initial configuration of the membrane. As seen in the examples outlined above, this can be addressed by biasing other degrees of freedom to obtain multi-dimensional PMFs. Furthermore, as described in a recent paper by Neale and Regis [136], many studies have assumed that the properties of interest converge rapidly within short simulation times, often avoiding an in-depth assessment of convergence for each window or the entire PMF. However, convergence may require simulations in the microsecond scale [99,136,183]. Similarly, Lee et al. [99] also stressed the importance of allowing sufficiently long equilibration times (50 -100 ns) for each window and suggested using sampling windows across the whole depth of the bilayer and "symmetrising" the data afterwards.

Metadynamics
Besides US, metadynamics (MetaD) is one of the most extensively used enhanced sampling method for biomolecular simulations. MetaD accelerates rare events and thus allows calculation of the free energy surface (FES) in feasible simulation times [184,185]. In this method, a history-dependent potential along one or more collective variables (CVs) is used to prevent the system from resampling previously visited configurations. The bias usually takes the form of a sum of repulsive Gaussian energy functions. Like for RCs, the CVs are assumed to describe the largest degrees of freedom of the process of interest such that all significant energy barriers in the FES can be sampled.
The biasing potential as a function of the CVs is defined as Here ( ) describes the CV as a function of the coordinates of the system, ( ) = ( ( )) is the value of the CV at time t, is the Gaussian height, is the Gaussian width and is the frequency with which Gaussian functions are added. For sufficiently long simulation times, the biasing potentials added enable the system to exit every energy minima and sample neighboring ones, leading eventually to the free diffusion of the system in configurational space (Figure 18).
From the resulting biasing potential the underlying, unbiased free energy surface can be obtained as F(S)= -VG(S) + C, where C is an arbitrary additive constant [184]. In a single MetaD simulation, the FES does not converge to specific values, but rather fluctuates around the correct values, leading to a statistical error proportional to the square root of the Gaussian deposition rate [185]. Because of this, it can be non-trivial to decide when to end a simulation. On the other hand, if the simulation is run for too long, the continued addition of the biasing potential can push the system into exploring non-relevant and physically meaningless parts of the FES. This issue is mostly circumvented by using well-tempered metadynamics (WT-MetaD).
Here the biasing potential is modified such that the height of the repulsive Gaussian functions is decreased over time in a manner that is proportional to the overall biasing potential already added at a given point of CV space [185].
Simulating permeation using metadynamics Ghaemi et al. used bias-exchange metadynamics (BE-MetaD) to describe the permeation of ethanol across a POPC bilayer [186]. This variation of MetaD allows the efficient calculation of a multidimensional FES as a function of several putative CVs and is particularly suitable to describe complex processes such as permeation. Whilst the system is biased along a large number of CVs, a subset of relevant CVs sufficient to describe the process can be chosen afterwards. In this study, seven CVs were used which describe the position of ethanol in the membrane along the bilayer normal, the torsional angles of the ethanol molecule, and the contacts of ethanol with water, the lipid head groups and the lipid tails. The resulting FES was then combined with a position-dependent diffusion matrix to perform a long kinetic Monte Carlo simulation to produce a realistic trajectory of the permeation process on arbitrarily long time scales. However, the permeation coefficient calculated from this kinetic model was significantly higher than those obtained from experimental data for the permeation of ethanol in POPC liposomes [187]. Ghaemi et al. subsequently extended this approach to more complex permeants: efavirenz and etravirine, which are two lipophilic anti-HIV drugs [60]. In this case more than 10 CVs were used, which can be grouped as follows: i) one CV for the position of the drug along the COM distance between solute and membrane, ii) three CVs to describe the interaction of the drug molecule with its surroundings, as defined by the number of contacts with water, phospholipid head groups and lipid tails; and iii) several CVs describing the internal conformational degrees of freedom in the drug molecules defined by dihedral angles. concentration, but a much smaller increase comparing the 20 to the 40% system [189]. This study along with some others [190] seem to suggest that the cholesterol could affect the permeation of hydrophobic and hydrophilic molecules differently. Bochicchio et al. [191] compared the effectiveness of MetaD, WT-MetaD and US to calculate free energy of partitioning of dimers of polyethylene and polypropylene into a POPC bilayer. For both molecules, the free energy profiles of the partitioning obtained with the three methods were in close agreement. However, both MetaD approaches were more efficient than US in terms of total computational resources and had a lower statistical error associated with the free energy. Results of Jambeck et al [192] demonstrate how WT-MetaD can provide a complete description of the translocation process of three common drug compounds (aspirin, diclofenac and ibuprofen) across a DOPC bilayer. Upon insertion, the change in dielectric media triggers a high-energy transition in the conformation of these molecules (in terms of dihedral angle distributions) resulting in the formation of an internal hydrogen bond. Considering the low probability of observing such highenergy transitions, US estimates will be likely be subject to incomplete sampling. Their results for ibuprofen are within 1 kBT of the experimental binding energy.
Voth and collaborators compared MetaD, WT-MetaD and the recently developed transitiontempered MetaD approach in terms of the accuracy of the PMF of ethanol and trimethoprim [193,194]. Transition-tempered MetaD is a variation of WT-MetaD specifically designed to study state to state transitions in which the states of interest are known ahead of time, but the transition mechanisms are not. The variation lies in the way the height of the repulsive Gaussian function is tempered, which in this method becomes a function of a dynamically defined percolation path between the states under investigation. In other words, the tempering of the Gaussian functions, and hence the decreasing speed of convergence, starts to operate only after a minimum free energy path that goes from state A to state B has been sampled. All four methods were used to calculate the free energy profile for the permeation of ethanol and the anti-bacterial drug trimethoprim across a DMPC bilayer. Figure 20 compares the PMFs obtained for ethanol by these three MetaD approaches with respect to the US approach used as a reference. In the early stages of the simulations (120 ns) transition-tempered MetaD shows the best agreement whereas WT-MetaD required 220 ns to reach agreement. Figure 20. Permeation of ethanol across a DMPC bilayer. Z is the centre-of-mass distance between the molecule and the membrane, whereas Θ is an orientation parameter that defines the angle between the vector that connects the oxygen to its adjacent carbon and the normal vector to the membrane. Left panel: Definition of two CVs used to describe the process and enhance sampling, namely the COM distance between ethanol and the bilayer, and the orientation of ethanol (O-C vector) with respect to the bilayer normal. Right panel: Differences in free energy profiles obtained from different MetaD methods during the early stages of the simulation, with US taken as the reference. Adapted from [193], Copyright 2016, American Chemical Society.
Replica exchange MD simulations Replica exchange MD (REMD) is a well-established approach to enhancing configurational sampling [195]. The method consists of simulating several copies (replicas) of the same system in different perturbed conditions whilst the unperturbed copy acts as the reference system from which the properties of interest are calculated. The perturbed condition is chosen such that each replica can explore a larger volume of phase space. During the simulations, neighbouring replicas are allowed to swap their configurations in a manner that satisfies the detailed balance condition. In this case, the probability of exchange between replica m and a replica n is given by: Here ( ) is the probability of finding a state with coordinates in replica m, and ( → ) is the transition probability. The most common approach is to allow the system to exchange configurations with a certain probability according to the classic Metropolis criterion, such that ∆ = ( − )( − ) with βm = 1/kBTm and Em being the energy of the system m at the moment of exchange [196] with probabilities In this way configurations derived from regions of phase space usually not accessible by unbiased MD simulations can gradually make their way to the unperturbed original system, which would normally be the biologically or chemically relevant replica. The exploration of higher energy configurational states is usually achieved in one of two ways: i) by simulating each replica at increasingly higher temperatures, in which case the method is referred to as 'parallel tempering' or temperature-REMD (T-REMD); or ii) the Hamiltonian of the system is perturbed in a similar incremental fashion, a class of techniques generally known as 'Hamiltonian tempering'.

Parallel tempering
In parallel tempering the variable used to perturb the system is the temperature. Low temperature replicas allow an accurate sampling of energy minima at physically relevant condition but may become kinetically trapped in such minima during the timescale of a typical simulation [196]. Higher temperature replicas have higher kinetic and potential energies and are thus used to sample larger volumes of phase space. The main drawback of parallel tempering lies in its high computational cost due to the large number of replicas needed. This becomes particularly problematic for simulations of large systems with explicit solvent (like those used for studying SMMIs). The large number of replicas are required because of the exchange acceptance ratio described in Eq. 14. In parallel tempering the probability of exchange is exponentially related to the difference in the total energy of the system (including solvent molecules), and for this probability of exchange to be acceptable (i.e. not too small) the two replicas need to exhibit a non-zero overlap in their potential energy distributions [196]. Large systems usually have a large number of solvent molecules, which usually translates into a very narrow potential energy distribution. Consequently, many replicas are required to correctly sample a large span of temperatures. This method, although historically important and rigorous, has seen little application in SMMI simulations, mainly due to the fact that at high temperatures the membrane incurs the risk of become unstable. Nonetheless a number of derived methods have been developed and which can be used to investigate SMMIs, namely Hamiltonian tempering simulations.

Hamiltonian tempering
Hamiltonian tempering simulations constitute a general class of methods that include a variety of different techniques, each differing in the way the Hamiltonian of the system is perturbed [197][198][199][200][201][202][203]. The general principle of these methods is that instead of each replica having a different temperature, as in the case of parallel tempering, each replica in a Hamiltonian tempering simulation uses a perturbed Hamiltonian, which is usually a linear combination of the starting Hamiltonian HA and the target Hamiltonian HB. The aim of the Hamiltonian perturbation is the same as the scaling of the temperature in parallel tempering: to explore higher energy regions and thus sample larger volumes of phase space. HB is usually chosen to flatten energy barriers on some or all degrees of freedom, allowing thermal fluctuations at the given simulation temperature to overcome the energy barriers that impede the sampling of phase space. Exchange between replicas and the subsequent statistical analysis is almost the same as for parallel tempering, but with some important added advantages. The first one is that directly modifying the Hamiltonian of the system allows the selective perturbation of subsets of atoms or molecules in the system, such as solutes or a part of a protein, leaving the solvent unperturbed. The second advantage is that the exchange probability between neighbouring replicas can be modified such that it is based only on the difference in energy between the relevant portions of the system [200]. For example, in the popular replica exchange with solute tempering (REST) method, leaving out the energy terms for the solvent from the functional form of the exchange probability reduces the computational effort by an order of magnitude, decreasing the number of replicas needed and making this method feasible for larger systems [198]. Figure 21 schematically illustrates the differences in the potential energy distribution across replicas between classical parallel tempering and REST. The efficiency of this class of methods in regard to SMMIs has not been fully explored, probably due to the lack of easily accessible and flexible implementations.
Neale and colleagues [204] used a replica exchange simulation exchanging the umbrella sampling potential between neighbouring windows in a n-propylguanidinium-POPC bilayer system. A particular, non-synchronous variation of replica exchange was used to enhance the convergence of the calculated PMF, which was called virtual replica exchange. This method consists of the asynchronous exchange of structures and Hamiltonians between replicas and a target sampling ensemble based on previously saved instantaneous energy values of the target ensemble [205]. This method allowed to increase the statistical convergence of the calculated PMF by a factor of three.
Even in the case of virtual replica exchange, the PMF was observed to remain non converged even after using 1.2 µs of simulation per window, probably because of large barriers in orthogonal degrees of freedom in the centre of the membrane, which gives rise to correlation times of up to 10 µs [204]. Hamiltonian tempering can achieve a wider distribution of the energy, such that a good overlap is obtained with a smaller number of replicas. The distribution is larger because the solvent-solvent term, which has a very narrow distribution around its mean value due to the large number of molecules, is cancelled out. Bottom panel: In parallel tempering, because of the narrow distribution of the solvent-solvent energy, the distributions of all replicas are much narrower because of the solvent-solvent term; therefore, a greater number of replicas is required to achieve a relatively small overlap in the energy distributions.
Other ways to enhance sampling We have discussed so far the most widely used and validated enhanced sampling methods to characterise SMMIs and calculate binding affinities and permeation coefficients. Although these methods are regarded as the most reliable approaches to addressing sampling problems, other lesser-known methods can be used effectively in specific situations. We summarise them briefly below.
Multi-scale methods use coarse-grained (CG) models to reduce the overall number of pair-wise calculations, thus allowing longer simulation times whilst still describing the region of interest with an all-atom (AA) representation for better accuracy. For SMMIs, this means having a CG representation of the bilayer and solvent and an atomistic representation of the solute. Studies using this approach have covered a range of SMMIs [56,206,207]. This approach relies on force field parameterisation that is compatible with CG-AA interactions. Another way to make use the computational speed-up afforded by coarse graining is to run independent CG simulations and systematically switch the entire system to an all-atom representation, in what is called the multiscale approach to conformations sampling [208]. This approach relies on accurate back mapping to AA representation and the premise of truly uncorrelated CG. For a concise summary of the multi-scale simulation of biomembranes the reader is referred to the review by Pluhackova and Böckmann [209].
Alternative membrane models can be used instead of an all-atomistic or coarse-grained description.
One of these alternatives is to describe the membrane as a spatial continuum, i.e. as a region that has a different dielectric constant to the solvent part of the system. This approach is utilised by several implicit membrane models [210][211][212], which has mostly been used to model the structure and/or folding of membrane-bound proteins [213,214]. The use of implicit solvent models greatly reduces the number of atoms of the system and is thus computationally less expensive than allatomistic or coarse-grained descriptions of a membrane. However, by design these approaches lack the ability to explicitly deal with specific lipid-molecule interactions, some of which are often crucial for an accurate description of the thermodynamics and kinetics of the system. Therefore, before using any of these methods it is important to assess carefully the properties of the system that are of interest and how their simulation could be affected by neglecting specific lipid-molecule interactions. A related approach is the use of a hybrid of atomistic and simplified continuum models, such as the highly mobile membrane-mimetic model [215]. This model has been specifically developed to accelerate the lateral diffusion of lipids in membranes, a quantity which is known to be highly underestimated by classical atomistic simulations [215]. The highly mobile membranemimetic model divides the membrane into a core region, made up of an organic solvent, encaged by two surface regions composed of short tailed ordered lipids that emulate the typical interfacial region formed by the lipid head groups. This type of approach has been successful at describing surface peptide-membrane interactions and predicting free energies of adsorption [216,217], but the method has limited capacity to accurately calculate permeation coefficients or, generally, describe SMMIs where a molecule penetrates into the membrane core [218].
Another method that can be used for simulating the permeation of small molecules is the milestoning approach [219]. In this method, instead of defining RCs or CVs, a reduced reaction space is defined to describe the permeation process, which requires prior sampling of the system. The reduced reaction space is defined by a selected group of states that globally describe the process.
These are denoted anchors and are in effect points in the N-dimensional phase space. Milestones are set between anchor domains (spaces between anchors) to differentiate when the system has made a transition from one anchor to the next (therefore, milestones are hypersurfaces in phase space). The method is based on storing the length of time that the system takes to transition between anchors, i.e. to go from milestone Mij to Mjk following a short simulation trajectory that, when played backwards from Mjk arrives at a different milestone than Mij itself, i.e. the particle crosses three different milestones (Figure 22). When this is done multiple times with independent trajectories, a flux q can be calculated for each milestone and a kinetic network emerges from the transitions between milestones. Similarly to the probability function in US, the permeability coefficient can be calculated from the kinetic network. The set of anchor domains (system states) provides a view of the process in a reduced reaction space and thus its choice requires prior knowledge of the system undergoing said process [220]. Thorough testing should be conducted for CV candidates and their different relaxation times, such that use of computing time is targeted to CVs that affect the accuracy of the estimates (usually slow DOFs). We discussed how, even for molecules as small as dipeptides there are around a dozen internal DOFs describing the conformational space and requiring sampling. These are potential hidden variables of the permeant alone! With the software and computing power currently available it is easier to simulate larger permeants interacting with bilayers. Beginners in this field will hopefully be able to recognise some considerations when undertaking the simulation of SMMIs and their subsequent data analyses. Increased awareness of the effects of insufficient sampling in the study of SMMIs will also enable them to critically approach the wealth of literature in this field.
With the examples presented in this review we have highlighted the challenges and drawbacks of the use of different enhanced sampling methods for the study of SMMIs. We recognise, however, that there is not a single way to approach the complex issues involved in the simulation of SMMIs, with the most efficient approaches being system-dependent.
Estimates of the amount of computational resources required to reach convergence (e.g. of the free energy calculations) are, unfortunately, system-dependent. Therefore, it is difficult to provide a comparison of efficiency across different methods. The most effective CV for sampling a SMMI process is usually not known a priori and usually requires a substantial amount of testing before assessment of convergence can be carried out. Once the best CVs have been determined, methods like MetaD can be quick to reach convergence. In comparison, REMD methods require minimal prior testing, but depending on the size of the system, can rapidly become prohibitive. This is due to the large number of replicas necessary to link the ground state with the high energy states (particularly for larger systems) or because of an excessively large amount of conformational states available at the higher energy states. An ideal enhanced sampling method for a given system should be chosen carefully by taking into account both experimental and mechanistic information available a priori and indeed a good dose of physical intuition.
It is also important to stress that achieving convergence in the calculation of free energy profiles of permeation may require simulations in the microsecond scale. Despite these challenges, insights gained from molecular simulations have increased our understanding of the process of permeation of small molecules and have indeed been instrumental in the development of the inhomogeneous solubility-diffusion model. Future work in the simulation of SMMIs will likely remain focused on improvements of enhanced sampling methods as well as the move towards simulations of more complex membranes that are a more realistic representation of actual biological membranes. This will require in turn the development of force fields for a larger range of lipid classes that can represent their properties accurately in more diverse membrane environments.

Conflict of interests
There are no conflicts to declare.

Acknowledgements
CM and LRP acknowledge the award of a Curtin International Research Scholarship.