{'fraction_non_alphanumeric': 0.028233374133006937, 'fraction_numerical': 0.01379028967768258, 'mean_word_length': 4.607992678462477, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Simulation studies of stratum corneum lipid mixtures

We present atomistic molecular dynamics results for fully hydrated bilayers composed of ceramide NS-24:0, free fatty acid 24:0 and cholesterol, to address the effect of the different components in the stratum corneum (the outermost layer of skin) lipid matrix on its structural properties. Bilayers containing ceramide molecules show higher in-plane density and hence lower rate of passive transport compared to phospholipid bilayers. At physiological temperatures, for all composition ratios explored, the lipids are in a gel phase with ordered lipid tails. However, the large asymmetry in the lengths of the two tails of the ceramide molecule leads to a fluid like environment at the bilayer mid-plane. The lateral pressure profiles show large local variations across the bilayer for pure ceramide or any of the two component mixtures. Close to the skin composition ratio, the lateral pressure fluctuations are greatly suppressed, the ceramide tails from the two leaflets interdigitate significantly, the depression in local density at the inter-leaflet region is lowered, and the bilayer have lowered elastic moduli. This indicates that the observed composition ratio in the stratum corneum lipid layer is responsible for both the good barrier properties and the stability of the lipid structure against mechanical stresses.


I. INTRODUCTION
Stratum corneum (SC), the outer layer of the skin [1], provides the main barrier against water loss [2] and invasion by foreign pathogens. During their life-cycle, cells formed in the basal layer of the epidermis change their shape and composition. They progressively occupy outer layers of the epidermis, until they are peeled off from the outer layer. SC is often viewed as a bricks and mortar structure [3] with corneocytes, the keratin filled nonviable disc like cells, arranged like bricks in a lipid mixture forming the mortar phase [4]. The three main components of the SC lipid matrix are a family of ceramide sphingolipids (CER), cholesterol (CHOL) and free fatty acid (FFA) [5,6]. Selective inhibition of any one of the ceramide [7], cholesterol [8] or FFA [9] is known to compromise the barrier function of the skin. However, how the three components affect the lipid matrix properties at the molecular level is not known. To our knowledge, there have been only a few previous attempts at atomistic modeling of the SC lipid layer. Most notably, Höltje et al. [10] simulated a two component mixture of fatty acid and cholesterol. Pandit and Scott [11] simulated a bilayer composed of symmetric CER NS 16:0 molecules. Notman et al. [12] used atomistic simulation to investigate the effect of DMSO molecules on a hydrated bilayer composed of CER NS 24:0. But a systematic study of the effects of the three components is lacking and this study is specifically concerned with understanding the interplay of the three components which endows the skin with an almost contradictory combination of pliability and an extremely high penetration barrier.
There are at least 9 different classes of ceramide found in human SC, with minor modifications in the head group region and the addition of an esterified fatty acid in the case of Ceramide 1. All the ceramides are conspicuous by having a large asymmetry in the length of the two tails and a large polydispersity in the fatty acid tail lengths [13]. Similar polydispersity is found also in the length of the free fatty acids [14]. Realistic representation of such a complex collection of molecules with atomic details is beyond current computational capabilities. Instead we choose just one representative ceramide, ceramide NS (also referred to as ceramide 2), with an asymmetric but monodisperse tail length. Ceramide NS is the most abundant species among the ceramide family. Its fatty acid tail is chosen to be 24:0, guided by the relative abundance of the different tail lengths in human SC [13]; while its sphingosine motif is chosen to have 18 carbons. Similarly we choose only FFA 24:0 because it is the most abundant free fatty acid found in SC lipid layer [14]. Between the corneocytes, the lipid matrix shows regular electron density variations, similar to lipid multilayers. This is not necessarily the only possible arrangement in the SC lipid matrix. In vitro experiments show the possibility of asymmetric leaflets [15] and multiple layer thicknesses, with indications that ceramide 1 connects different bilayers [16]. The lipids are predominantly in a gel phase, possibly a single continuous gel phase [17] or with fluid regions [18]. In the skin, the corneocytes are arranged in clusters with the lipid matrix extending through the full depth of the SC at intervals [19]. These regions show a much lower permeation barrier than the layers between corneocytes [19] and it is not known whether or not the lipids there are arranged in a multilayer structure. There is a hydration gradient across the SC, with the average water content around 30% by weight [20,21]. How the water molecules are arranged in the SC is not completely known.
With this uncertainty about the arrangement of the lipids in SC, and with the computational limitations im-posed by detailed atomistic simulations, we study lipid mixtures arranged in symmetric bilayer structures in excess water. An atomistic potential is chosen in this study to have a faithful representation of the lipids involved. In the absence of prior atomistic simulations, or, experimental results on hydrated SC lipid bilayers, a coarse-grained description is not appropriate for this system. The water is required to stabilize the bilayer structure, and can be viewed as a replacement for the layering field imposed by the flat corneocytes. The lateral arrangement of the components and the properties of such a hydrated bilayer will probably be representative of the lipid matrix confined between the corneocytes.
The rest of the paper is organized as follows: In the section Simulation Details we outline the force field and protocols used for the simulations, and define the various properties measured. Next we provide the results for a wide range of compositions, concentrating mainly at temperature T=340K, which is much higher than the skin temperature. However, much of the results were found to remain qualitatively the same even at the skin temperature ≃ 300K. We present the results for the higher temperature because we are more confident of the equilibration at 340K than at skin temperature T≃300K. Finally we summarize the main findings and present an outlook for further studies on SC lipid structures.

A. Force field
The interaction parameters used in the simulation are based on the united atom OPLS force field [22] with modifications for the nonpolar CH 2 /CH 3 groups [23] that accurately reproduce experimental quantities for lipid molecules [24]. The polar hydrogens were included explicitly and the partial charges for the headgroup atoms were selected to conform with molecules having similar structures and simulated using the same force field in the literature. The dihedral potentials in the hydrocarbon tails were described by the Ryckaert-Bellemans term [25].
The force field used for the ceramide molecules has been used previously [12]. The topologies of the polar part of the fatty acid and the cholesterol are the same as used in [10]. The water is modeled with the SPC potential [26].

B. Simulations
We used extended ensemble molecular dynamics simulation at constant temperature and pressure (NPT) ensemble with GROMACS molecular dynamics package [28,29,30]. Temperature was controlled by Nosé-Hoover thermostats [31,32] coupled separately to the lipid and water molecules with time constant 5ps. Pressure was FIG. 1: Skeletal representation of the atomic arrangements of the ceramide (CER NS 24:0), the free fatty acid (FFA 24:0) and cholesterol (CHOL). Partial charges used in the simulation are shown in brackets. Inset: (a) V-shaped and (b) hairpin arrangements of the two tails of a CER molecule. [27].
controlled by an anisotropic Parrrinello-Rahman barostat [33,34] with time constant 5ps and compressibility 4.5 × 10 −5 /bar. The off-diagonal terms of the compressibility matrix were set to zero to preserve the orthogonal shape of the simulation box. Standard periodic boundary conditions were applied in all three directions. A simple group-based cut-off was used for calculating electrostatic interaction. The cut-off length was chosen to be 1.2 nm for both the Van der Waals and electrostatic potentials. Because the lipid molecules considered in this study do not fully ionize, the dominant electrostatic interaction is the dipole-dipole interaction, with the dipoles made up of the partial changes separated by fixed bond lengths. This is quite different from the case of phospholipids, where fully ionized phosphate group introduces a large dipole moment from well separated N − and P + charges. As a specific example, the C25=O26 atoms on CER molecule (see fig. 1) form a dipole from equal charges of magnitude 0.28 e separated by a fixed distance of 0.123 nm. This creates a dipole moment of 1.65 D, which is roughly 15 times smaller than the dipole moment at the head groups in DPPC bilayer [35]. Noting that the dipolar interaction falls of as 1/r 3 and that for DPPC bilayer effect of electrostatic interaction with the periodic images become negligible for system sizes ∼ 1024 lipid molecules [35], we estimate that for the dipole moments involved in our system, the effect of electrostatic interaction with the periodic images become neglibile for systems containing ∼ 168 lipids. We have carried out simulations with pure CER molecules and 2:2:1 composition system both with group-based cut-off and particle mesh Ewald summation (PME) for several system sizes. The results are presented in the appendix and validates the independence of the results presented on long range electrostatic interaction and system size. Because the use of PME requires more than twice the computational cost as compared to using cut-off and it has no effect on the validity of the results presented here, we restrict ourselves to group-based cut-off for handling electrostatic interaction in the simulations presented in the main paper.
All bond lengths were constrained using the SHAKE algorithm [36]. The analytic SETTLE algorithm was used to handle the rigid SPC water molecules [37]. We used a timestep of 2 fs for T ≤ 340 K and 1 fs for high temperature simulations. During the the production run, configurations at intervals of 0.5 ps were stored for further analysis.
The SC lipid matrix (in mice) shows a pH of ∼ 6 [38]. The ionization state of FFA is known to modify the structure of CHOL-FFA (C16:0) mixtures [39]. For ceramidecholesterol mixtures containing palmitic acid (C16:0) and oleic acid (C18:1), reported values of the pK a lie in the range of 6.2-7.3 [40]. Typically longer tailed FFA are expected to have higher pK a ; e.g., for hexacosanoic acid (C26:0) in an egg-phosphatidylcholine bilayer pK a ∼ 7.4. This suggests that a fraction of the FFA should ionize at skin pH. However, the experimental evidence for the effect of pH on SC lipid mixture is not clear. Surface force and AFM measurements on SC lipid mixtures containing free fatty acid seem to show no detectable effect in a wide range of pH (between 3.0 and 7.0) [41]. The apparent insensitivity to pH for SC lipid mixture may be because the bilayers containing long tail ceramides remain in a dense gel phase and do not allow much structural freedom to the FFA molecules. With this background, we chose to simulate only un-ionized FFA -not having free charges allow fast calculation of the forces with cut-off as alluded before.

C. Initial structures
In all but one of the six different crystal structures formed under varying conditions, the sphingosine and the fatty acid tails of CER NS 24:0 arrange themselves with a large opening angle V-shaped structure [27]. Only low temperature crystallization from solution leads to a hairpin arrangement of the tails (inset of Fig. 1). To ensure a hairpin arrangement of the ceramides, we start with a multilayer anhydride system with the ceramides placed in a slightly expanded hexagonal lattice with a hairpin structure at 300K, and relax the configuration with a 20 ns NPT MD run. A bilayer from the multilayer was then placed in a larger box and the box was filled with water. We ran the system for 5 ns with the ceramide molecules fixed (to hydrate the bilayer properly), and then for 2 ns with only the terminal methyl group on all the lipid tails frozen, so that the rest of the molecules reorient themselves to accommodate the water environment. This configuration was used as the initial configuration for the pure CER system. A further 10 ns NPT MD run was used at each temperature before the production runs.  For mixtures containing CER molecules, we transform the required number of CER molecules, chosen at random but maintaining the same composition in the two leaflets, to either CHOL or FFA. To convert to FFA molecules we simply separate the two tails of CER as two separate molecules and slowly grow the shorter chain derived from the sphingosine motif to the required length. FFA molecules are significantly more mobile as compared to CER molecules. Thus, although in the initial configuration the two FFA molecules generated from the same CER molecule are next to each other, they do not introduce any artificial correlation after the equilibration step. To convert to CHOL we map certain atoms of the CER to CHOL. Repeated short NVT MD runs were performed while keeping the CHOL molecules frozen, and after each such short run, the atoms of the CHOL were displaced by a small amount until they reached the equilibrium positions on the molecule. The mixtures were adequately equilibrated, typically by energy minimization followed by a series of short NVT simulations and finally with at least 10 ns of an NPT MD run. To find a typical relaxation time, we quenched a CER bilayer from 360 K to 300 K and found that the area/lipid and the bilayer thickness re-equilibrate to 300 K values in ∼ 1 ns time-scale. For pure FFA or CHOL bilayers, we simulated multilayers and placed one bilayer in water in the same way as in the case of CER. Table I shows the number of molecules used, as well as the corresponding molar ratios, rounded to closest integer ratios. Some of the composition ratios explored in this study may not be achievable in experiments as hydrated bilayers. However, once prepared, the bilayers are kinetically stable and allow us to isolate the effects of the different components by studying the extreme compositions. In the rest of the paper the different compositions are referred to by the corresponding molar ratios.

III. MEASURED QUANTITIES
In this section we define the different quantities measured from the simulations. For all composition ratios and temperatures, we use averages over 10 ns runs. We estimate the error bars for the measured quantities from the variance of intermediate averages over 2 ns windows.
We define the bilayer normal direction to be the z direction. The orthogonal box shape along with the small system size ensures that this direction remains the same as the z direction of the simulation box. All z-dependent quantities reported are averaged over the lateral direction (x − y plane) in the entire simulation box.

Structural properties
At the molecular length scale, the water-lipid interface has a finite width. To assign an unique value to the bilayer thickness 2d, we calculate the density of the water molecules as a function of z and define the position of the interface between the lipid and water as the z at which the water density decays to 1 e of the bulk water density. We use this criteria, as opposed to the more usual criteria of the Gibbs dividing surface, because we consider multi-component lipids in this study and the different components have significantly different densities at the boundary. Concentrating on the approximately exponentially decaying water density gives a simple and unique method to assign a value to the bilayer thickness. The average lipid layer densityρ L is computed by assuming that all the lipid mass is homogeneously distributed between the two lipid-water interfaces. The local density of the lipid molecules depends on z with a minimum ρ min L between the two leaflets.
The asymmetry in the lengths of the two hydrocarbon tails of CER can lead to significant interdigitation. We define the following dimensionless overlap parameter as a quantitative measure of interdigitation, where ρ t (z) and ρ b (z) are the densities at z for the CER from the top and bottom layers respectively. ρ ov (z) = 1 if half of the density at z is from the top layer and the other half is from the bottom layer CER molecules. If only the top or bottom layer of CER is at z, ρ ov (z) = 0.
Integrating over z, we define a single length scale λ ov ≡ L 0 ρ ov (z)dz to compare the amount of interdigitation for different compositions. The integration is carried over the whole box, since if no CER is present, then ρ ov (z) = 0 and there is no contribution to λ ov .
The area compressibility κ A of the bilayer is related to the area fluctuation in the NPT ensemble by [42] where k B is the Boltzmann's constant and T is the temperature. The angular brackets refer to averages over time for the area A and square of the area A 2 .
The system sizes investigated in this work are too small to calculate the bending modulus by analyzing the undulation modes [43,44]. For fully saturated fluid membranes, the bending modulus κ is approximately related to the area compressibility modulus κ A and the thickness of the hydrocarbon tail region 2d c through where the constant c e is estimated to be 24 from a theory based on polymer brushes [45]. This simple prescription is found to work well for a number of different lipid bilayers [45,46]. In the present simulation, the head groups of the lipid molecules are much smaller than that of phospholipids and we use the bilayer thickness 2d as a good estimate of the hydrocarbon thickness 2d c .

Tail order parameter
The orientational (nematic) order of the tails is probed through an order parameter P 2 , defined by the largest eigenvalue of the second rank tensor whereû iα is the cartesian component α of some specific orientation vector on the lipid molecule i, and the < > denotes average over time and the lipid molecules. For fluid bilayers, P 2 is related to the NMR deuterium order parameter S CD through the relation P 2 = −2S CD [47]. The eigenvector corresponding to the largest eigenvalue of Q defines the average orientation. Different carbon atoms along the tail are expected to have different amount of ordering. To probe the local ordering, we define a local orientation P 2 (n) at carbon atom n on the lipid molecule i by the vector directionû iα (n) between carbon atoms n − 1 and n + 1.

Local pressure
The anisotropic arrangement of lipid molecules in a bilayer leads to an anisotropic local pressure profile, and hence local stress variations. While macroscopic pressure is a well defined quantity in molecular simulations, its microscopic description is not unique [48]. We use the formalism of [49] to define the local pressure tensor at a given height z as (5) Here, v i is the velocity of particle i, F ij is the force on particle i due to particle j and r ij is the relative position vector of particle i from particle j. V slice is the volume of a thin slice with thickness ∆z and < > denotes an average over time. The first sum runs over the particles in the slice centered at z. The function f determines the contribution from the virial (the second sum) to the current slice. f is unity when both particles are in the current slice. If one or both of the particles are outside the slice (but the shortest distance between the two particles goes through the slice); f is respectively chosen to be ∆z/|z i − z j | and ∆z/ (2|z i − z j |). Because constraints in the simulation (SHAKE algorithm) transfer some of the kinetic contributions to the constraint potential, we need to consider both the sums in Eq. 5 explicitly. We calculate the local pressure profiles by starting with stored configurations separated by 2 ns and re-evolving the configurations for 200 ps.
Of particular interest for lipid bilayers is the difference between the lateral and the normal pressure δP (z) = P LAT (z) − P zz (z), where P LAT (z) = 1 2 [P xx (z) + P yy (z)]. The surface tension γ of the bilayer is related to δP (z) through [49] Anisotropic box fluctuations ensure that the bulk water is isotropic, i.e. the pressure tensor has equal diagonal components and zero off-diagonal components. This in turn leads to a zero average surface tension for the lipid bilayer. But locally δP (z) goes through a number of maxima and minima. We define a microscopic stress bȳ Here, the integration is over the bilayer thickness and the angular brackets reflect time averages.

IV. RESULTS AND DISCUSSION
A. Pure Ceramide bilayers   2 shows snapshots of the CER bilayer at (a) 300K and at (b) 360K. The tails of the CER align along the z direction with strong nematic ordering. The terminal methyl groups at the bilayer midplane show a liquid-like disordered arrangement. With increasing temperature, the disordered region at the bilayer midplane gradually thickens. The tail order parameters at 300K and 360K are shown in Fig. 3. The terminal groups (smallest and largest atom indices) have small ordering. P 2 also decreases close to the head group (atom indices 16 and 25). The C-C bonds close to the head group orient at an angle to the vertical direction. Fig. 4 (a) shows a top view of one leaflet. The figure shows that the molecules arrange in a zigzag fashion with rows of molecules having the head group arrangement in alternate (orthogonal) directions. Thus the tail carbon bonds close to the head group point in one of these two orthogonal directions, reducing the overall average of P 2 . In the well-ordered layer immediately beneath the head groups, the tails have fairly well-defined hexatic order with very few defects (Fig. 4b). Moreover, the orientations of the planes containing the lipid tails (Fig. 4  c ) suggests that the phase of this ordered layer is intermediate between crystal and rotator phases (familiar from studies of alkanes of different lengths [50]). While the two dimensional slices show strong sixfold order, the molecules undergo relatively fast slithering motion of the tails (Fig. 2 c). The coupling between the two leaflets is weak, allowing the two leaflets to oscillate about each other by more than the inter-lipid distance in roughly 10 ns. Diffusion of the center of masses of the molecules in the x − y plane is much slower because they can only diffuse at the defect sites. Fig. 5 shows the density of the CER and water molecules across the bilayer at 300K (open symbols) and 360K (closed symbols). The CER density shows a peak at water-lipid interface from the close packing of head group atoms. There is an almost constant density shoulder from the dense packing of the hydrocarbon tail atoms, followed by a dip at the bilayer midplane, which corresponds to the amorphous inner layer due to the asymmetric ceramide tails. The arrow indicates where the water density falls below 1/e of the bulk water density at 300K, which is considered to be the lipid-water interface to calculate the bilayer thickness.   6 shows the area per lipid and bilayer thickness as a function of temperature. The area/lipid at 300K is 0.389 nm 2 , which is close to the X-ray result ∼ 0.4 nm 2 for CER bilayer with molecules in hairpin arrangement [27]. With increasing temperature, the area/lipid gradually increases, while the increasingly disordered terminal tail atoms make the thickness decrease. The average density and area compressibility (Fig. 7) both decrease with temperature. Even at 360K, the area compressibility is an order of magnitude larger (∼ 4000 dyne/cm) than in most phospholipid fluid bilayers. The results do not show a sharp transition. Instead the disordered inter-leaflet region expands smoothly with temperature. The long chains involved in CER probably will lead to a broad gradual softening like a waxy material, instead of a sharp transition. Experimentally, main transition temperatures in ceramide systems have been found, for multilayer systems, to be of order 394 − 420 K [27]. One expects a lower temperature for a single hydrated membrane in solution. For all the composition ratios investigated the temperature behavior in the simulated range of temperatures is gradual and similar to that of CER system. (Fig. 8 shows results for the composition CER:CHOL:FFA=2:2:1). For presentation of the results in this section, we concentrate on 340K. A summary of the different measured quantities is included in Table II. Fig. 9 shows snapshots (a) for pure CER, (b) equimolar CER-FFA, (c) equimolar CER-CHOL and (d) 2:2:1 mixture of CER, CHOL and FFA. The tails of CER molecules (Fig. 9a) retain substantial ordering. Long chain FFA molecules (Fig. 9b) fall in registry with the CER molecules. The slightly longer length of the FFA molecules is accommodated by partially increasing the tail order in the CER molecules. The tails in the leaflets arrange in a slightly tilted orientation with respect to the layer normal. Thus, even though the amount of ordering is the same on both leaflets, the particular orientation of Fig. 9 highlights the ordering in the lower leaflets at the cost of obscuring the order in the upper leaflets.
The head groups of the CHOL (Fig. 9c) mostly stay at the water-lipid boundary. The shorter length of the CHOL molecules squash the bilayer, with the central re-composition 2dρL ρ mid L λov κA (10 3 ) κ (10 −11 )ǭP (molar ratio) (nm) (g/cc) (g/cc) (nm) (dyn/cm) (erg)  gion comprising primarily CER tails. The ordering in the CER tails is decreased compared to pure CER bilayers and the tails from the two leaflets overlap strongly. The molecular arrangements in the ternary systems (Fig. 9d) are somewhat intermediate to the binary CR2-FA and CR2-CHOL mixtures. Single chain FFA molecules are more flexible than CER molecules, which induces more FFA atoms in the midplane disordered phase and thus increases the tail order of CER molecules to some degree compared to CER-CHOL mixtures. To investigate the arrangement of the different components further, in Fig. 10 we plot the density of the lipid components along the z direction for different compositions. Fig. 10a shows the density profile of the CER atoms. Both pure CER and CER-FFA systems have a high density near the head groups, followed by a constant density region due to the ordered tail atoms and then a region of lower density at the midplane, which covers almost 2 nm. For CER-CHOL and three component mixtures, the density from the tail atoms of the CER molecules remains almost constant throughout the bilayer. In fact, for CER-CHOL mixtures, there is an increase in the local CER density at the bilayer midplane due to increased interdigitation of the molecules from the two leaflets. The plateau densities in Fig. 10 are fixed by the relative abundance of the different molecules in the mixtures. Fig. 10b shows the density profile of the FFA atoms. For pure FFA and CER-FFA mixtures, the density profiles are qualitatively similar to that of pure CER bilayers, except the region of low density inter-leaflet space is narrower. In the 3-component mixtures, the density minimum is replaced by a local density maximum because of inter-leaflet overlaps. The density profile of the CHOL molecules (Fig. 10c) does not have the constant density hydrocarbon regions present for the FFA and CER molecules. The profile for the pure CHOL bilayer is not symmetric and the interleaflet density quite high, signifying significant transfer between the two leaflets during the simulation time scale. Fig. 10d shows the total lipid densities. The effect of fatty acid is to increase the bilayer thickness and reduce the density by a small amount. CHOL, when present, reduces the bilayer thickness and increases the density at the tail region.
The overlap of the lipids from the two leaflets (partial interdigitation) is expected to increase the inter-leaflet friction and couple the dynamics of the two leaflets. Fig. 11 shows a surface plot of λ ov , a measure of interdigitation, as a function of composition. The three component mixture is characterized by large λ ov . One of the SC lipids, not considered in this study, ceramide 1 (ceramide EOS) contains an additional long-chain ω-hydroxy acid (with number of carbon atoms > 30) linked to the fatty acid tail. This effectively makes the length of the long tail of ceramide 1 nearly double to that of the other members of the ceramide family. The presence of ceramide 1 will probably introduce significantly more interdigitation and inter-leaflet coupling.   12 shows the tail order parameter of CER atoms at 340K for different composition ratios. The presence of CHOL reduces the nematic order, while FFA increases the order. The presence of either CER or CHOL reduces the order of the FFA atoms compared to a pure FFA bilayer (Fig. 13). This differs from phospholipid membranes, where the planar shape of cholesterol typically increases the nematic order of the phospholipid tails. When mixed with the long ceramides, the shorter cholesterol molecules tend to encourage a thinner membrane by disordering the longer ceramide tails so that they can fill the space around the cholesterol. This also accounts for the increased overlap of the CER tails upon adding cholesterol, as depicted in Fig. 11.  The bending modulus κ, calculated from the polymer brush theory [45] using κ A and the bilayer thickness, behaves in a similar fashion. Close to the skin composition CER:FFA:CHOL=2:2:1 the bilayer becomes softer, with comparatively smaller κ A and κ. The absolute magnitude of the elastic constants remain much higher than in the fluid phase of phospholipid bilayers.   15 shows the difference between the lateral and the normal pressure δP (z) along the bilayer normal direction. Inside the bilayers, pure CER (Fig. 15a) shows large variations in δP (z) as a function of z. Similar but less pronounced variations in δP (z) are observed in phospholipid membranes. For example, Lindahl and Olle [44] found a maximum variation of |δP (z)| < 500 bar for a DPPC bilayer at 323K, which is 5 times smaller than the present case. The addition of CHOL (Fig. 15b) introduces additional peaks in δP (z). FFA reduces the peak height marginally (Fig. 15c). The three component mixtures (15d) show much less pronounced variations in δP (z) than the pure ceramide or any of the 2 component mixtures. A large δP can be interpreted as a local moment acting on the molecules, making a deviation from the flat interface of bilayer more likely. The magnitude of δP also can be viewed as a local stressǭ P . Fig. 16 shows a ternary plot ofǭ P defined through Eq. 7. The vertical line shows the position of composition 2:2:1, whereǭ P shows a minimum. To gain some insight into the energy scale, we note that 10 3 bar ∼ 0.7k B T per methyl group, while the trans-gauche energy difference is ∼ 2k B T per dihedral bond. Thus, in the absence of prohibitive energy barriers, the system can introduce gauche defects to reduce the local variation of δP .

V. CONCLUSIONS
We have presented molecular dynamics simulation results for various composition ratios of the three main constituent components of stratum corneum lipid layers, namely, CER NS 24:0, FFA 24:0 and CHOL. The long asymmetric tails of the CER molecules form a dense bilayer phase in water, with considerable interdigitation of the tails from the two leaflets and strong nematic order. CHOL, being rigid and smaller in length compared to CER, acts as a molecular clamp, squashing the bilayer and increasing the in-plane density. The resulting structure shows large lateral stress variations, which is relieved by the presence of FFA. The three component mixtures are characterized by a higher density and comparatively smaller area compressibilities and bending moduli.
The lipids in real skin stratum corneum have large polydispersity in tail length and have different head groups. This work does not address the effect of polydispersity. We also study bilayers in excess water, which is quite different from limited water environment in the skin. The strong local pressure fluctuations in the bilayer suggests that the bilayer may not be the most stable structure and it will be interesting to simulate the lipids in limited water conditions. The hydrocarbon densities in the bilayers were found to be quite large compared to phospholipids. This is probably responsible for three orders of magnitude smaller permeability of water through stratum corneum as compared to plasma membranes.

APPENDIX A: EQUILIBRATION TIME
To probe the typical time it takes for the system to adapt the equilibrium conformation for a given temperature, we used an equilibrated configuration of 128 CER lipids (1:0:0 composition) at 360 K and monitored the instantaneous area/lipid as a function of time (small black circles in Fig 17) after abruptly changing the desired thermostat temperature to 300 K. We also show the area/lipid calculated from separate runs at 360 K and at 300 K as shaded big circles with error bars at the two extreme ends of the data. As can be seen from the plot, the area/lipid adapts to the new temperature with a time scale of around 1 ns (position of the arrow). In the paper we use equilibration times of 10 ns for each 10K temperature shift -which gives ample time for the molecules to reorganize themselves into the equilibrium conformation corresponding to the set temperature.

APPENDIX B: EFFECT OF LONG RANGE ELECTROSTATICS AND SYSTEM SIZE
To explore the dependence of measured quantities on system size and run times, we have carried out some additional simulations on pure CER bilayers (at 300 K) and bilayers with 2:2:1 composition ratio of CER, CHOL, FFA (at 340 K). These simulations use much larger number of lipid molecules and much longer run times than the simulations reported in the paper. For these simulations, we use both the group-based cut-off and PME methods to probe the effect of long-range electrostatics. In simulation with PME, the grid spacing was chosen to be 0.1 nm and fourth order polynomial interpolation was used. A summary of the number of molecules used and timescales probed is given in Table. III. To get the initial configuration for a larger system, we use the already equilibrated configuration at the required temperature from the simulations presented in the main paper and replicate it in the x and y directions. To reduce disk space requirements, in most cases of these set of simulations, we save configurations only at intervals of 0.2 ns (as opposed to 0.5 ps in the main paper). Thus we expect larger statistical errors in the data reported in this appendix as compared to the main paper.  In Table IV, we report the energies (normalized by the number of lipid molecules in the system) for 2:2:1 system calculated separately with group based cut-off and particle mesh Ewald summation (PME). The two important observations from the energy values are that (i) the difference between the group based cut-off and PME schemes in the total energy is about 2% and (ii) the normalized energies are independent of the system sizes. Both of these two observations show that the interactions with the periodic images do not contribute significantly in the electrostatic energy. Table IV also shows the area compressibility and the average bilayer thickness calculated from these simulations. The effect of the system size and handling of electrostatics do not affect the values significantly.  IV: Short range (SR) and long range (LR) electrostatic energies, total energy, area compressibility (κA) and bilayer thickness (2d) for bilayers containing 2:2:1 molar ratio of CER2, CHOL and FFA at 340K. Statistical errors of the mean values are indicated in the brackets as the uncertainty on the last digit. ( * Energies are in MJ/mol and are normalized by the number of lipid molecules in the system. † κA is in 10 3 dyn/cm.) Table V shows the effect of system size and method of calculating electrostatic interaction for pure CER bilayers at 300K. The conclusions about insensitivity of the results on both the system size and the long range electrostatic interaction in the case of 2:2:1 lipid mixture remains equally valid for the case of CER bilayers.