Parameterization of a coarse-grained model of cholesterol with point-dipole electrostatics

We present a new coarse-grained (CG) model of cholesterol (CHOL) for the electrostatic-based ELBA force field. A distinguishing feature of our CHOL model is that the electrostatics is modeled by an explicit point dipole which interacts through an ideal vacuum permittivity. The CHOL model parameters were optimized in a systematic fashion, reproducing the electrostatic and nonpolar partitioning free energies of CHOL in lipid/water mixtures predicted by full-detailed atomistic molecular dynamics simulations. The CHOL model has been validated by comparison to structural, dynamic and thermodynamic properties with experimental and atomistic simulation reference data. The simulation of binary DPPC/cholesterol mixtures covering the relevant biological content of CHOL in mammalian membranes is shown to correctly predict the main lipid behavior as observed experimentally.

Unfortunately, all-atom (AA) simulations are still far from the temporal and spatial scales of experimental studies.Consequently, supra-atomic coarse-grained (CG) simulations arise as an attractive alternative to explore such systems beyond scales achieved by atomistic simulations.
In the 1980s and 1990s, Ipsen and coworkers [59] and Miao and coworkers [31] developed some of the earliest CG models in an attempt to understand the behavior of cholesterolphospholipid mixtures.Because these models corresponded to two-dimensional model descriptors of phosphatidylcholine-cholesterol systems interacting by a phenomenological Hamiltonian, they allowed phosphatidylcholine-cholesterol phase diagrams to be constructed and coexistence lines between the So (solid, ordered-chain), Lo (liquid, ordered-chain), and Ld (liquid, disorderedchain) phases to be determined.
Murtola and coworkers [60] took a step forward when they proposed a two-dimensional model to represent the dipalmitoylphosphatidylcholine-cholesterol (DPPC/cholesterol) mixture.
These authors obtained the effective interaction potential by the IMC technique from AA-MD simulations, which allowed complex lipid mixtures to be exploited on a large length scale.The original model was further improved in subsequent work [61,62].
In the same spirit, Khelashvili and co-workers [63] used a self-consistent mean-field (SCMF) model based on MD simulations to study a two-dimensional representation of the DPPC/cholesterol model system.Their study enabled complex lipid models to be simulated on large spatial (500 nm in lateral size) and temporal (around dozens of microseconds) scales.De Meyer and Smit also studied both cholesterol and phospholipids by means of dissipative particle dynamics (DPD) via CG models.This particular work provided a new understanding of DMPC/cholesterol mixture phase diagrams [64].

Hadley and McCabe proposed a new cholesterol model for CG simulations based on the
IBI technique [65].This model correctly described the cholesterol crystalline solid state approximately 180 times faster than analogous AA simulations.The authors also explored the self-assembly of a set of skin lipids and investigated how the cholesterol content affected these systems [66].
Izvekov and Voth [54] introduced a new view of CG based on the FM technique.They mapped dimyristoylphosphatidylcholine (DMPC), cholesterol, and water molecules into fewer site bead representations-thirteen beads, four (or seven) beads, and one bead, respectively.CG-MD simulations of a hydrated DMPC bilayer containing cholesterol molecules were 50-100 times faster than AA-MD simulations of a system of equivalent size.
Marrink and co-workers [36] proposed a new version of the popular MARTINI Force Field (so-called MARTINI-FF 2.0).They suggested a protocol to model ring structures with specific applications to cholesterol.Therefore, the cholesterol molecule was modeled with a total of eight beads-six beads for the rigid sterol body and two beads for the hydrocarbon tail.These authors also simulated mixtures containing DPPC and cholesterol; they covered a broad range of cholesterol content, from 0 mol% to 60 mol%.MARTINI-FF development stimulated subsequent research into biomembrane mimetic systems containing cholesterol [36,[67][68][69][70][71][72][73][74][75][76][77][78][79][80][81][82][83].Moreover, Daily and co-workers carried out a systematic optimization of structural parameters based on atomistic reference simulations, to improve MARTINI-FF for lipids [84].This revised model predicted cholesterol ordering and thickening effects on the palmitoyloleoylphosphatidylcholine (POPC) bilayer better as compared to the original MARTINI-FF.
Shen and co-workers [85] recently proposed a new CG model for cholesterol based on the CAVS-FF, in which two opposite charges embedded into different sites describe the cholesterol electrostatics.This model correctly predicts the effects of cholesterol concentration on membrane dipole potential.
Here, we present a point-dipole CG model for cholesterol (designated CHOL hereafter) that is consistent with ELBA-FF (ELectrostatic BAsed Force Field) [37].We also provide an ELBA-FF extension for the DPPC lipid model, derived from such a set of optimized parameters available for the POPC lipid model [86].A key feature of ELBA-FF is that it describes electrostatic interactions through charges and dipoles rather than using the standard practice of most CG FFs which employ point charge interactions shielded by an effective dielectric constant or incorporates it all into a short range potential.Thus, ELBA-FF is especially suitable to model molecules with different polar groups because it provides a realistic mapping of their polar moieties through explicit point dipoles.In particular, ELBA-FF allows functional mapping of CHOL polar moieties through point dipoles.
This paper is organized as follows.The first section (Computational Methodology) reports the MD simulation details and analysis protocols.The subsequent sections address CHOL model parameterization based on physicochemical properties obtained from experimental and atomistic reference data.We then compare our predictions for the CHOL and DPPC self-diffusion coefficients and the free energy of CHOL desorption membrane partitioning in pure and CHOLrich DPPC bilayers to the AA and CG reference simulations.

Coarse-grained simulations
MD simulations were carried out with LAMMPS (http://lammps.sandia.gov)version February 16, 2016 [87].Initial coordinates for the binary DPPC bilayer systems at 10, 20, 30, 40, and 50 mol% CHOL content (amounting to a total of 128 lipids) were generated with the MOLTEMPLATE program [88].Two water layers were placed at the top and bottom of the lipid bilayer systems, which corresponded to a hydration level equal to 50 water molecules per lipid molecule, in agreement with typical experimental conditions.Additionally, a simulation box containing one CHOL molecule solvated by 8000 ELBA water molecules and another simulation box containing 1000 CHOL molecules were also constructed.
VMD [89] and MEMPLUGIN [90] were used to visualize and to analyze the MD trajectories.Shifted-force interactions for the Lennard-Jones and electrostatic potentials were used with truncation at a cutoff radius of 12 Å.Integration time steps of 10 fs and 8 fs were employed in the absence and presence of CHOL molecules, respectively.The velocity-Verlet integrator [91] was used for all the CG MD simulations.The simulation temperature was maintained at 303 K by using a Langevin thermostat with a damping coefficient of 1 ps -1 .Semi-isotropic pressure coupling was used to control the pressure at 1 atm [92] with an oscillation period of 1 ps and a compressibility of 4.6 x 10 -5 atm -1 .After initial energy minimization, followed by an equilibration step for 0.5 microseconds, production runs reached 1.6 microseconds for all of the studied systems.
Original ELBA-FF parameters for water and lipids were taken from Orsi and Essex [37].The new set of ELBA-FF parameters for CHOL is available in the Supplementary Information (SI).

Atomistic simulations
AA MD simulations were carried out with the NAMD 2.9 program [93] using CHARMM36-FF for lipids [94,95] and the TIP3P water model [96].The Particle Mesh Ewald (PME) method [97] was employed for long-range electrostatic interactions, and 1-4 non-bonded interactions were turned on.The initial configuration for the hydrated DPPC bilayer system containing a CHARMM36-based CHOL molecule [94,95] was constructed with CHARMM-GUI Membrane Builder [98].In addition, a simulation box containing 32154 TIP3P water molecules solvating a CHARMM36 CHOL molecule [94,95] was generated by using Solvate tools of the VMD program [89].Electrostatic and Lennard-Jones potentials were used with a cut-off radius of 12 Å.The SETTLE algorithm [99] was used for rigid water molecules, whereas Newton's equations of motion were solved by means of the MTS integrator Verlet-I [100] with a 2 fs time step.

Alchemical free energy -Thermodynamic Integration (TI)
The electrostatic term of the solvation free energy was calculated with the aid of Thermodynamic Integration (TI), as implemented in both LAMMPS [87] and NAMD [93] packages.In particular, the electrostatic free energy was estimated through the following equation: , where   is the electrostatic potential energy, and λ is the coupling parameter with limits of integration between 0 and 1.The ⟨  ⁄ ⟩  term is the average derivative of the electrostatic potential energy as a function of λ.
Numerical integration through a cubic spline interpolation over ten equally spaced λ points was used to compute the electrostatic free energy.The CHOL molecule was harmonically restrained by a force constant equal to 100 kcal/mol Å 2 in the mass center of both bulk water and the DPPC bilayer.This allowed us to evaluate the electrostatic free energy difference of CHOL desolvation in the membrane.
The statistical uncertainty of the means was estimated by block averaging over the entire λ interval.Standard errors were plotted as a function of the block length.The plateau value calculated by averaging values in the limit of infinite block length gave an estimate of the 'true' standard error [101,102].

Mean force potential -Adaptive Biased Force (ABF)
Free energy profiles were calculated by the Adaptive Biasing Force (ABF) approach [103][104][105] implemented in both the LAMMPS [87] and NAMD [93] programs.By definition, the average force exerted on the transition coordinate, , can be related to the free energy first derivative according to: where |J| is the Jacobian determinant to transform the Cartesian coordinates to generalized coordinates, and ⟨  ⟩  denotes the cumulative average of the instantaneous force   in bins with the size of .In the ABF approach, the force applied to overcome energetic barriers along the reaction coordinate can be defined as: where ⟨  ⟩  denotes the cumulative mean of   , calculated for each step of MD simulations.  represents the biasing force that ideally reaches zero in its convergence.
For ABF calculations, the reaction coordinate was discretized in consecutive windows with size equal to 5Å.A "non-charged" version (only Lennard-Jones interactions are active) of ELBAbased and CHARMM36-based CHOL molecules [94,95] was dragged from bulk water toward the bilayer center in their respective systems.A constant force equal to 10 kcal/mol Å 2 in increments  of 0.1 Å (ABF-AA simulations) and 1.0 Å (ABF-CG simulations) was employed.
The statistical error of the mean force for each particular bin can be obtained from the standard error of the mean according to the following relation [106]: where 〈  2 〉  represents the variance of the instantaneous force collected to bin i , and   is the independent samples.The number of independent samples   can be obtained by the product between the total number of simulation steps,   , and the MD time step, , divided by the autocorrelation time,   , of the instantaneous force,   =      ⁄ .In particular,   values for each specific window were obtained by non-linear fitting of an exponential function  −  ⁄ .The Bienaymé formula [107] was used to calculate the error of the sum of mean forces for points in different windows:

Thermodynamic calculations
We obtained the total free energy difference of CHOL partitioning in both pure and CHOL-rich DPPC bilayers by employing the thermodynamic cycle shown in Scheme 1.
Scheme 1. Thermodynamic cycle employed to calculate the total free energy of CHOL partitioning in lipid bilayers: We assumed that the initial (n) and final (m) states are the CHOL molecule harmonically restrained in the middle of bulk water and different depths in the bilayer, respectively.

Partial volume and area per lipid
Determination of the partial area and volume per lipid in complex lipid mixtures by MD simulations has received considerable attention in recent years [27,108].Atom-based tessellation methods have been widely applied to obtain the partial area per lipid in the perpendicular plane of bilayers composed of complex lipid mixtures.Unfortunately, standard tessellation analysis for molecular area calculations in lipid bilayers depends on an arbitrary choice of which normal coordinate position of the bilayer plane should be tessellated.
The methodological approach taken in this study was based on the Hofsäß method [24], which provides a useful relation between bilayer thickness and molecular volumes to obtain lipid partial areas in binary mixtures.
A straightforward strategy to split the total area per lipid a into its components aA and aB, assuming a binary mixture of lipids A and B, is given by: where x is the molar fraction of lipid B.
The Hofsäß method [24] can be used to split the total area per lipid into partial areas.
Subsequently, Edholm and Nagle proposed a variation of this method to improve the analysis of the volumetric part [109].The central premise of the Hofsäß method [24] is to use the mathematical relation between area and volume to define a standard thickness ℎ() for the lipid bilayer as follows: where  is the total volume of the simulation box,   is the total number of water molecules,   is the volume of a water molecule and  is the total area.
For the specific case of a binary lipid mixture of lipids A and B, the partial area for A and B as function of CHOL  can be expressed as: For convenience, one can define a function () that represents the volume ratio between   and   , and then, one can write: The standard analysis of the partial area per lipid does not take account of the free volume related to the lipid packing in bilayers.The Voronoi tessellation [110] is especially well-suited to avoid this problem.In short, this technique consists of defining several cells surrounding a particular set of particles getting the space closest to them.Thus, we have therefore made use of Voronoi tessellation via Voro++ program [111] in LAMMPS to capture the effective volume   and   occupied by the lipids in the simulation box.

Cholesterol tilt angle distribution
The tilt angle distribution was obtained by the polar angle (given in degrees) between a vector joining specific beads in the CHOL and DPPC molecules and the outward vector normal to the bilayer within the range of 0° to180° (see Section 4 in SI).

Lipid lateral self-diffusion coefficient
Lipid lateral self-diffusion coefficients were obtained from the center-of-mass displacement of DPPC and CHOL molecules in the two-dimensional space parallel to the bilayer.The mean square displacement (MSD) as function of time is given by which is related to the self-diffusion coefficient  by where  is the number of lipids, and   ( 0 ) and   ( +  0 ) are the position vectors of the lipid  center-of-mass in the bilayer plane at the starting measurement time  0 and at the limit measurement time  +  0 , respectively .The bracket notation was satisfied with the linear fit of lipid MSD as a function of time in a normal (linear) diffusional regime.Any "drifting effect" in the center-of-mass of the analyzed lipid groups was removed before the displacement of each of their particles was calculated.The lipid diffusion average and its uncertainties were calculated by block averaging [101,102].

Parametrization protocol and topology
Initially, an appropriate topology for the CHOL molecule was defined.The CG mapping scheme for the CT3-Me2b CHOL model designed by Daily and coworkers [84] for the MARTINI-FF [36] was chosen.Figure 1 shows the molecular structures of CHOL at different levels of resolution.In particular, the pseudo-ring portion of the CHOL model was treated as a rigid body (Figure 1, beads 1 to 5).In this case, such an approximation is possible due to small distortions in the rigid CHOL ring structure.Therefore, numerical instabilities in the time propagation of motion equations were avoided by eliminating high frequency vibrational modes and improper dihedral potentials after the equilibration phase.
The parameterization protocol is similar to that used in MARTINI-FF [36].Electrostatic and non-polar components of solvation free energies obtained from CHARMM36-FF [94,95].MD simulations were used to set the Lennard-Jones (LJ) epsilon parameters of CHOL-water and CHOL-DPPC cross-interactions.
Regarding the other LJ parameters, default values from the ELBA-FF [37] for the ester and non-polar chemical blocks were used as starting point to describe both the hydroxyl-like and the non-polar CHOL beads, respectively.Similarly, the dipole moment present in the hydroxyl-like bead was refined by comparison with electrostatic free energy results from CHARMM36-FF [94,95] MD simulations (Table 1).
Parameterization of DPPC-CHOL cross-interactions and CHOL dipole adjustment considered simulation data from free energy calculations of CHOL partitioning in pure DPPC bilayer.Furthermore, CHOL-CHOL cross-interactions were adjusted on the basis of biophysical properties of high-cholesterol DPPC/CHOL (50:50 mol%) bilayer simulations.In addition, intermediate CHOL content values were predicted by using the optimized set of parameters that best described the structural and energetic properties of both systems cited above.

Optimizing CHOL model LJ cross-interactions
Kessel and co-workers [112] emphasized that a single CHOL molecule desolvation free energy and specific location within the membrane are mainly driven by non-polar free energy contributions from DPPC-CHOL interactions.The authors also pointed out that CHOL distribution within a bilayer is determined by the elastic response of their neighboring lipids.In addition, some authors such as Jo and Khelashvili [63,113] argued that CHOL orientational entropy as well as CHOL-CHOL and CHOL-phospholipid interactions play a crucial role in the free energy associated with the CHOL tilt modulus in heterogeneous lipid environments.
Here, particular attention was paid to calibrating DPPC-CHOL and CHOL-CHOL crossinteractions due to their direct relation with CHOL behavior within membranes.To optimize CHOL-CHOL and DPPC-CHOL cross-interactions, the original set of parameters for the polar and non-polar LJ sites available in ELBA-FF [37] were used as starting point.Default ELBA-FF [37] parameters were maintained for cross-interactions between ELBA water and polar and nonpolar DPPC and CHOL ELBA CG sites during parameterization.
First, the potential of the mean force (PMF) related to the free energy that was necessary to transfer an "uncharged" CHOL molecule from bulk water to the middle of the DPPC bilayer was obtained via CHARMM36-FF and ELBA-FF MD simulations.LJ epsilon values of CHOL-DPPC cross-interactions were scaled up by multiplicative factors (1.0, 1.05, 1.10, 1.15 and 1.20) over default parameter values used in ELBA-FF (the effect over PMFs can be viewed in Figure 2).
Therefore, the cross-interaction between CHOL and DPPC was adjusted by seeking agreement between the AA and CG profiles.Similarly, the CHOL-CHOL self-interaction was also adjusted by employing multiplicative factors (in this case 1.0, 1.05, 1.10, and 1.15) over the default set of ELBA-FF parameters.Our aim was to choose the set of parameters that best described CHOL tilts and the partial area per CHOL and DPPC as compared to atomistic reference data.

Optimizing CHOL model dipole moment
The magnitude of the dipole moment embedded in bead 1 (see Figure 1, item a) was fine-tuned to match the reference value found for the total free energy of the CHARMM36 CHOL [94,95] partitioning from bulk water into the hydrocarbon core of the DPPC bilayer [96].The optimal dipole (4.92 D) was obtained by scaling up the dipole vector that describes the ester group (2.00 D) in ELBA-FF [37] (for more details, see Section 2 in SI).
Table 1 summarizes AA and CG predictions for the Electrostatic and LJ free energy contributions for moving the CHOL molecule from bulk water to the middle of the DPPC bilayer.
Table 1 also provides a comparison between the CHOL model predictions and those obtained from the GROMOS87, CAVS, and MARTINI FFs [85,114] for the total free energy of CHOL partitioning in pure DPPC bilayers.

Structural properties
A characteristic effect of the presence of CHOL in bilayers consisting of mixed lipids is the socalled "cholesterol condensing effect", characterized by a sharp decrease in the phospholipid bilayer area (and free volume).Theoretical models such as the Superlattice Model [115,116] and the Umbrella Model [117,118] were proposed to explain such effect.The Superlattice Model successfully predicts highly regular distributions by symmetrical geometric assumptions of the system, as confirmed by experimental techniques [119][120][121].This model also suggests a link with the long-range repulsive interaction owing to the cross-sectional area difference between sterol and lipids.In turn, the Umbrella Model postulates that the shape complementarity of phospholipids "covering" the CHOL molecule avoids the free energy penalties of the CHOL molecules interacting with interfacial water in phospholipid bilayers.The latter model plays a crucial role in explaining the "cholesterol condensing effect" in phospholipids [117,118].
Regarding DPPC/CHOL mixtures, previous predictions obtained by mean-field calculations and MC/MD simulations demonstrated that the partial area per DPPC decreases with increasing CHOL content [122,123].Hofsäß and coworkers [24] and Edholm and coworkers [109] had already reported that CHOL addition to DPPC affects the partial CHOL area: areas decrease from 29 to 27 Å² and from 28 to 23 Å², respectively.Furthermore, Chiu and coworkers [23] found partial area per CHOL equal to 22.3 Å² in high-cholesterol regime.The discrepancy between the AA predictions for the partial CHOL area (Figure 3) results from their different forms to incorporate the local perturbation of CHOL in the partial area calculations.
In general, our CG DPPC/CHOL simulations correctly describe the reduction in the partial and total area per CHOL and DPPC lipids within the range of 0 mol% to 50 mol% CHOL.For the DPPC/CHOL mixture in the range of 0 mol% to 30 mol%, the partial area per DPPC decreases significantly, being almost constant for CHOL contents greater than 30 mol%.Such reductions in the partial lipid area and the consequent contraction in the bilayer plane correspond to the aforementioned "condensing effect of cholesterol" concept, which has been widely discussed in the literature.[124][125][126][127][128][129][130].
Concerning the area per lipid of the pure DPPC bilayer, experimental data have suggested a value equal to 62.9 Å² [131], whereas recent atomistic simulations have inferred values ranging from 59.1 to 62.9 Å² [94,132,133].Our prediction for the pure DPPC-ELBA bilayer model is 58.9 Å².
Table 2 highlights the average values of tilt angles Ɵ (in degrees) and distances along the  The average distance from CHOL molecules and CHOL-hydroxyl (OH) groups are, respectively, 11.6 Å and 16.8 Å to the bilayer center in the DPPC/CHOL (60:40 mol%) system.
These results are consistent with data acquired by experimental and simulation studies of binary DPPC/cholesterol bilayers containing equivalent CHOL amount [135,137].
The CHOL model reasonably predicts the tilt angle for CHOL molecules in DPPC bilayers with 40 and 50 mol% CHOL content, with average values equal to 18.4° and 17.7°, respectively.
Experimental prediction for CHOL tilts in the DPPC/CHOL mixture with CHOL at 40 mol% has afforded an average value equal to 16° [136].The CHOL model has been shown slightly higher tilt average values as compared to values obtained by AA simulations.Furthermore, DPPC hydrocarbon chains in DPPC/CHOL mixtures containing 40 mol% CHOL has shown an average value equal to 20.3°, which well agrees with AA predictions [136].
We also compared experimental volumetric data for pure and hydrated CHOL between the CHOL model and equivalent data from literature.Renshaw and coworkers [138] observed two values of molecular CHOL volume in bulk water: 627 Å³ and 606 Å³ for the anhydrous and monohydrated states, respectively.Additionally, an experimental bulk density equal to 1.1 g/cm³ has been found for pure CHOL at 298.15 K. Our simulation predictions for the hydrated CHOL molecular volume and for the bulk CHOL density are 613.1 Å³ and 1.2 g/cm³, respectively.These results are in line with the results cited above and contribute to validating our parameterization.

Dynamic properties
The dynamic properties of membrane lipids are intrinsically related to lipid distribution in cells, flip-flop diffusion, lipid raft formation, and membrane permeability.
The free-volume model [139] assumes that lateral lipid diffusion in bilayers is a result of the random occurrence of free volume portions of which a lipid in its vicinity takes place.On the other hand, Falck and co-workers [140] postulate that lateral lipid diffusion is a continuous and correlated motion of lipids.
Particularly, the free volume theory successfully explains why lipid diffusion decreases with increasing CHOL content [141].However, Falck and co-workers argue that more realistic models should be used to predict lipid diffusion in lipid bilayers quantitatively [108].
Subsequently, Almeida and coworkers [142] clarified the validity of the free volume model to describe experimental data for lipid diffusion in PC/CHOL mixtures quantitatively.
Scherfeld and co-workers [152] verified an increase in the lipid diffusion rates for CHOL contents ranging from 33 mol% to 67 mol% in Giant Unilamellar Vesicles (GUVs) of DPPC/CHOL mixture by Fluorescence Correlation Spectroscopy (FCS).Several experimental results corroborated that the lateral PC lipid self-diffusion coefficient diminishes upon increasing CHOL content [143,147,148,150].Moreover, Kuo and co-workers [158] noticed that less than 10 mol% CHOL content accelerates lipid diffusion in DPPC multilayers by a pulsed gradient proton Nuclear Magnetic Resonance (NMR) spin-echo study.The same authors found that the diffusion coefficient decreases in high CHOL concentrations.On the basis of the pulsed field gradient NMR technique, Filippov and co-workers [151] inferred a similar increase in lateral diffusion at small CHOL content in the DOPC/CHOL system.They also observed that lipid diffusion decreases substantially when CHOL concentration increases from 5 mol% to 15 mol% but remains constant thereafter.
Experimental methods like pulsed NMR and FRAP work in the millisecond timescale, whereas contemporary MD simulations usually do not reach such long timescale.Thus, direct comparison of experimental and simulation results concerning lipid diffusion should take this into account.
Most MD simulation studies agree that lipid diffusion coefficients decrease with increasing CHOL content in PC/CHOL mixtures [24,108,159].
Remarkably, our CG simulation results (Figure 6) correctly describe most of the experimental observations cited above.Compared to CHOL, the DPPC diffusion rate is slightly lower for CHOL concentrations up to 30 mol%.As for CHOL concentrations greater than 30 mol%, the DPPC and CHOL diffusion rates are not significantly different.In addition, our DPPC/CHOL model at low CHOL content come up with diffusion rates around an order of magnitude faster than those found at CHOL-rich regime, whereas there is no significant difference in MARTINI predictions [74,160].the data derived from atomistic simulations are indicated in blue [108] and green [24] colors, respectively.For the sake of clarity, the upper inset shows the CHOL and DPPC diffusion rates on a logarithmic scale as function of CHOL concentration.
The fact that CHOL and DPPC have similar diffusion rates at high CHOL content can be interpreted on the basis of the hypothesis of Hofsäß and co-workers [24], who argued that DPPC molecules hold on to almost the same vicinity of CHOL molecules along time, thereby leading to mutual diffusion between them.

Thermodynamic properties
Free energy calculations via MD simulations have been employed to understand how CHOL affects organic solute transport through bilayers.Jedlovszky and Mezei [161] explored how CHOL impacts the transport of small solutes through the membrane, to observe that CHOL hydroxyl groups are located within the membrane plays a pivotal role in the solvation free energy of polar and apolar solutes.
Despite the vast literature on CHOL modulation of biophysical properties and small solute transport in biomembranes, few studies have provided a quantitative picture of the thermodynamic process involved in CHOL partitioning in biomembranes.
Jo and co-workers [113] suggested that an unfavorable enthalpic component due to CHOL  The CHOL model provides a fair prediction for the free energy of CHOL partitioning in membranes consisting of pure and cholesterol-rich DPPC bilayers as compared to previous literature simulation data.

Figure 1 .
Figure 1.Molecular structures of CHOL at the AA and CG resolutions.(a) AA model of the

Figure 2 .
Figure 2. Potential of mean force for a "non-charged" CHOL molecule, related to the reversible

Figure 3
Figure3shows the total and partial area per DPPC and CHOL for different CHOL mole

Figure 3 .
Figure 3.The total (diamond) and partial molecular area of DPPC (triangle) and CHOL (square) bilayer normal D Z (in Å) for the CHOL molecules and DPPC hydrocarbon tails in DPPC/CHOL bilayers containing 40 and 50 mol% CHOL, respectively.The first parameter is closely related to CHOL ordering abilities and orientation within the bilayer.The latter parameter helps to validate the CHOL location normal to the bilayer.

Figure 6 .
Figure 6.Lateral self-diffusion coefficient for DPPC and CHOL as a function of different CHOL

(
hydroxyl) polar group desolvation is balanced by a favorable entropic factor arising from enhanced lipid dynamics promoted by CHOL.Moreover, Bennett and coworkers[114] found a total free energy for moving a single CHOL molecule from bulk water to its equilibrium position in DPPC/CHOL bilayers (60/40 mol%) equal to -21.3 kcal/mol and -20.6 kcal/mol when they used AA (GROMOS87-FF) and CG (MARTINI-FF) models, respectively.These authors also found a favorable free energy difference for the CHOL partitioning process in the middle of DPPC/CHOL bilayers (60/40 mol%) -around -12.7 kcal/mol and -17.0 kcal/mol -by using the GROMOS87-FF and MARTINI-FF, respectively.Additionally, the PMF calculated by Díaz-Tejada and coworkers[80] for lateral CHOL partitioning from a pure DPPC moiety to another moiety containing DPPC at 50% CHOL content gave a total free energy difference of around 3.6 kcal/mol.To calculate the total free energy of CHOL partitioning in DPPC with CHOL at 40 mol% content, we first estimate the Electrostatic and LJ free energy difference of CHOL desolvation at its equilibrium position in the bilayer, to achieve ∆   ( → ) = −0.8± 0.4 kcal/mol and ∆   ( → ) = −21.4± 2.3 kcal/mol, respectively.

Figure 7 .Table 3 . 5 CONCLUSION
Figure 7. Permeation free energy profile for a "non-charged" CHOL molecule going from bulk

Table 1 .
Electrostatic and LJ contributions for the total free energy of CHOL partitioning in pure denotes the total free energy cost for moving CHOL from bulk water to the bilayer center.

Table 2 .
Average values of tilt angles (Ɵ) for CHOL molecules and DPPC hydrocarbon tails, and distances from bilayer center (  ) to pseudo-atom groups of CHOL molecules in DPPC/CHOL bilayers at two different CHOL concentration (40 and 50 mol% CHOL).