Using molecular simulation to understand the skin barrier

Skin’s effectiveness as a barrier to permeation of water and other chemicals rests almost entirely in the outermost layer of the epidermis, the stratum corneum (SC), which consists of layers of corneocytes surrounded by highly organized lipid lamellae. As the only continuous path through the SC, transdermal permeation necessarily involves diffusion through these lipid layers. The role of the SC as a protective barrier is supported by its exceptional lipid composition consisting of ceramides (CERs), cholesterol (CHOL), and free fatty acids (FFAs) and the complete absence of phospholipids, which are present in most biological membranes. Molecular simulation, which provides molecular level detail of lipid configurations that can be connected with barrier function, has become a popular tool for studying SC lipid systems. We review this ever-increasing body of literature with the goals of (1) enabling the experimental skin community to understand, interpret and use the information generated from the simulations, (2) providing simulation experts with a solid background in the chemistry of SC lipids including the composition, structure and organization, and barrier function, and (3) presenting a state of the art picture of the field of SC lipid simulations, highlighting the difficulties and best practices for studying these systems, to encourage the generation of robust reproducible studies in the future. This review describes molecular simulation methodology and then critically examines results derived from simulations using atomistic and then coarse-grained models.

. Structure and nomenclature for the 24 CER subclasses identified in the unbound lipids from human SC. Two-tailed CERs are designated as CER ZFAZSB where ZFA and ZSB represent by one or two letters the fatty acid and sphingoid base, respectively. Three-tailed CERs are designated as CER E_S where E designates a non-hydroxy fatty acid ester-linked to the primary hydroxyl of the sphingosine base (S), which is connected by an amide bond to a second fatty acid, either the non-hydroxy (N) or alpha-hydroxy (A) fatty acid. The possible structure shown for dihydroxy sphinganine (T) is from Schmitt and Neubert [2].  Figure S1 (continued). Structure and nomenclature for the 24 CER subclasses identified in the unbound lipids from human SC. Two-tailed CERs are designated as CER ZFAZSB where ZFA and ZSB represent by one or two letters the fatty acid and sphingoid base, respectively. Threetailed CERs are designated as CER E_S where E designates a non-hydroxy fatty acid ester-linked to the primary hydroxyl of the sphingosine base (S), which is connected by an amide bond to a second fatty acid, either the non-hydroxy (N) or alpha-hydroxy (A) fatty acid. The possible structure shown for dihydroxy sphinganine (T) is from Schmitt and Neubert [2].

S2. Ceramide compositions in pig and human SC
Tables S1 -S5 provide the published compositions for CERs in pig and human SC determined using either thin layer chromatography (TLC) or liquid chromatography with mass spectrometry (LC/MS). Compositions in weight % are converted to mole % using the molecular weight (MW) values listed in Table S1, which are copied from Schmitt and Neubert [2]. Source in paper   Fig. 3e of Ishikawa et al. [18] as ng/µg of protein are reported in Table 1 of Kovacik et al. [19] as weight %, which are the numbers listed here. Kawana et al. [3] repeated the numbers from Kovacik et al. [19] in Table S6 without specifying weight %. g NS and NdS were not determined separately in the TLC analyses h AS and AdS were not determined separately in the TLC analyses i CER EOP was not identified in human SC until 2003 [17] j CER EOdS was not identified in human SC until 2011 [20]   e Mole % numbers were calculated from the weight % data using the molecular weight (MW) values listed in Table S1. f These same calculated weight % numbers are listed in Table 1 of Schmitt and Neubert [2]; however, the mole % numbers listed in their Table 1 were calculated incorrectly. g NS and NdS were not determined separately in the TLC analyses h AS and AdS were not determined separately in the TLC analyses i CER EOdS was not identified in human SC until 2011 [20] Fig. 3b of reference [25] and as weight % in Table 1 of reference [6] (although reference [6] did not specify the units, and the amounts for AdS, and EOdS, and possibly EOP, are larger than expected, perhaps due to typographic or copying errors). The mole % numbers from van Smeden et al. [25] are also listed in Table  S6 of Kawana et al. [3] without specifying mole %. Table 1 of Schmitt and Neubert [2] reports weight % numbers derived from the data listed in Table 1 of van Smeden et al. 2014 [6], which they adjusted to 100% after excluding EOdS (1.3%). Schmitt and Neubert [2] calculated the mole % numbers they included in Table 1 of their paper using molecular weight (MW) values listed in their paper (and presented here in Table S1). We recommend using the mole % numbers presented here, which are from Janssens et al. [23] d Although not stated in the paper, the results were reported as relative ceramide abundance (in percentage) of all CER subclasses measured (i.e., mol %); personal communication with the corresponding author K Sandra (email 02 September 2020). These same numbers are also listed in Table S6 of Kawana et al. [3] and in Table 1 of Schmitt and Neubert [2]. Kawana et al. did not specify mole % or weight %. Schmitt and Neubert incorrectly assumed weight % and therefore the mole % numbers listed in Table 1 of their paper are incorrect. Table 1 of Kovacik et al. [19] tabulates these same numbers for only the 12 most abundant CERs (i.e., without OS, OP, OH and NT), which sum to 96.9% (and not 100%); Kovacik et al. [19] incorrectly specify that the results are weight %. e Results were reported as relative ceramide abundance (in percentage) of all CER subclasses measured (i.e., mol %); this totals to 99.4%. f Calculated mole % of the 12 most abundant CERs.

S3. Comparison of simulations for bilayers with different SC lipid compositions
For systems with different lipid compositions, comparisons of area per lipid (APL) or normalized lipid area (NLA) should be performed using data from one study or between studies that used the same force field and computational protocol. Moore et al. [26] and Wang and Klauda [27][28][29] each generated simulation results, summarized in Table S6, of hydrated bilayers with other lipid compositions (i.e., CER NS C16, CER NS C24, and CER AP C24 alone or mixed with different amounts of CHOL and FFA C24) that can be compared. Table S7 compares APL and NLA for simulated bilayers of pure CER and equimolar mixtures of the same CER with CHOL and FFA C24. NLA values in these tables were calculated assuming the effective number of hydrocarbon tails per lipid is one for FFAs, two for CERs, and 1.9 for CHOL as proposed by Shamaprasad et al. [30]. Tables S6 and S7. First, adding CHOL to CER bilayers causes a minimal change in the APL; see results from Moore et al. [26] for CER NS. This is expected given that CER and CHOL have similar cross-sectional areas [31]. Second, as also expected, the addition of FFA, which has a single hydrocarbon tail, to either pure CER or CER-CHOL mixtures decreases the APL significantly with almost no effect on the NLA. Third, changes in the CER NS acyl tail length from C16 to C24 have no effect on the APL whether the bilayer consists of pure CER or a mixture of CER with CHOL or with both CHOL and FFA C24. Fourth, from the Wang and Klauda results listed in Table S7 [27][28][29], CER AP C24 has a larger APL compared with pure CER NS C24 (46.4 Å 2 and 42.8 Å 2 ), whereas the APL for equimolar mixtures of each with CHOL and FFA C24 are nearly identical (32.6 Å 2 and 32.8 Å 2 , respectively), suggesting that the steric hindrance caused by the additional hydroxyls in the CER AP headgroup is mitigated by the presence of CHOL and FFA C24. Moreover, this effect is observed even when smaller amounts of equimolar CHOL and CER are added to the CER (Table  S6); e.g., the APL values for AP and NS are identical for a CER:CHOL:FFA C24 mole ratio of 1:0.5:0.5 and are only slightly different for a mole ratio of 1:0.21:0.21 (37.4 Å 2 and 36.8 Å 2 for AP and NS, respectively). Consistent with the APL observations, the NLA values are essentially the same (~20 Å 2 ) for both CER AP C24 and CER NS C24 with mole ratios of 1:X:X CER:CHOL:FFA where X = 0.2, 0.5 and 1 (Table S6).

S4. Analysis of SC permeation predictions from Gajula et al.
Gajula et al. [32] assumed, like many others before them, that permeation through the SC can be represented as Fickian diffusion across a brick-and-mortar structure in which the corneocytes are the bricks and the lipid layers surrounding the corneocytes are the mortar. Gajula et al. assumed further that the corneocytes were impermeable and, although not stated explicitly, diffusion in the lipid layers is isotropic [32]. Mathematically, this scenario is described by Eq. (S.1) for two dimensions 2 2 where C is the concentration and Dlip is the diffusion coefficient of the permeant in the lipids, t is time, and x and z are the coordinate directions parallel and perpendicular to the SC surface, respectively. In their calculations, Gajula et al. used the diffusion coefficient calculated from molecular simulations for Dlip [32].
Because the corneocytes have zero permeability, there is no flux in the direction normal to the surface of each corneocyte. Thus, at the surface of all corneocytes. As a result, a numerical solution such as the finite element method is required to solve Eq. (S.1).
However, Kushner et al. [33] showed that this two-dimensional description of diffusion in the SC lipids surrounding impermeable corneocytes could be represented by the following onedimensional expression which combines Eqs. (S.1) and (S.3), for their assumed corneocyte-lipid geometry (see Figure  S2 and Table S8). As a result, their calculations for transport across the SC have accounted for the impermeable corneocytes and extended lipid pathway twice. Figure S3 shows the cumulative mass transfer per area (Q) versus time they calculated compared with experiments for caffeine, fentanyl and naphthol (reproduced by digitizing the plots presented in their Figure 6 [32]). Figure S2. Schematic diagram of the brick-and-mortar configuration Gajula et al. [32] assumed in their calculations (see Table S8 for parameter descriptions and numerical values).
We determined the Deff values in the Gajula et al. [32] calculations by fitting their Q versus t curves to the following expression  Table S8. Descriptions, numerical values, and defining equations for the parameters used by Gajula et al. [32] in their calculations of Q versus t (see Figure S2 a ω = 1 when the offset is symmetrical (i.e., the diffusion path around the left and right sides of a corneocyte are equal) and ω → ∞ when corneocytes are completely aligned. In Figure S2, ω = 3.5.  Table S9. Table S9 lists the tlag and Jss (Jss /C0 for naphthol) values determined by fitting the Q versus t curves from Gajula et al. [32] presented in Figure S3. The Q versus t curves calculated using these tlag and Jss (or Jss/C0) values in Eq. (S.6) closely match the curves from Gajula et al. (see dashed lines compared with solid lines in Figure S3). Table S9 also lists Deff calculated from tlag and the ratio of Deff with the Dlip values from Gajula et al., [32] which range from about 50,000 to 320,000. These are significantly larger than the expected values of 230 to 515 calculated from the tortuosity parameters for the assumed SC brick-and-mortar configuration presented in Figure S2 (see Table S9). However, as shown in Table S9, the Deff /Dlip values from Gajula et al. are within 25% of the square of the expected Deff /Dlip, which is consistent with a calculation that adjusted Dlip for the SC configuration twice. Had Dlip been correctly adjusted for the SC configuration only once, the resulting Deff (which range from 6.0 × 10 -13 m 2 s -1 for caffeine to 9.2 × 10 -13 m 2 s -1 for naphthol) would have been much too large to describe the experimental permeation data; i.e., tlag is less than 1 min and Jss is 200 to 600-fold larger than experimentally observed values (see the 'expected' values of tlag and Jss (or Jss /C0) in Table S9). Expected Deff (m 2 s -1 ) f 6.0 × 10 -13 7.1 × 10 -13 9.2 × 10 -13 Expected tlag (s) g 48 40 31 Expected Jss (µg cm -2 h -1 ) g 1100 3550 Expected Jss / Co (cm h -1 ) g 0.14 0.  Figure S3) to Eq. (S.6); for naphthol Jss /C0 was used instead because C0 was not specified. c Calculated from tlag using Eq. (S.7) d From Table S8