In Silico Prediction of Stratum Corneum Partition Coefficients via COSMOmic and Molecular Dynamics Simulations

Stratum corneum (SC) is the main barrier of human skin where the inter-corneocytes lipids provide the main pathway for transdermal permeation of functional actives of skin care and health. Molecular dynamics (MD) has been increasingly used to simulate the SC lipid bilayer structure so that the barrier property and its affecting factors can be elucidated. Among reported MD simulation studies, solute partition in the SC lipids, an important parameter affecting SC permeability, has received limited attention. In this work, we combine MD simulation with COSMOmic to predict the partition coefficients of dermatologically relevant solutes in SC lipid bilayer. Firstly, we run MD simulations to obtain equilibrated SC lipid bilayers with different lipid types, compositions, and structures. Then, the simulated SC lipid bilayer structures are fed to COSMOmic to calculate the partition coefficients of the solutes. The results show that lipid types and bilayer geometries play a minor role in the predicted partition coefficients. For the more lipophilic solutes, the predicted results of solute partition in SC lipid bilayers agree well with reported experimental values of solute partition in extracted SC lipids. For the more hydrophilic molecules, there is a systematical underprediction. Nevertheless, the MD/COSMOmic approach correctly reproduces the phenomenological correlation between the SC lipid/water partition coefficients and the octanol/water partition coefficients. Overall, the results show that the MD/COSMOmic approach is a fast and valid method for predicting solute partitioning into SC lipids and hence supporting the assessment of percutaneous absorption of skin care ingredients, dermatological drugs as well as environmental pollutants.


SASA and individual density analysis
Distributions of SASA from figures S2 to S7 are fitted with one (CHOL and LIGNP) or two (CERs) gaussians. In a few cases, CERs distributions showed only one peak and have been therefore fitted with only one gaussian. The fit has been performed via python3 with the lmfit 1 package. The gaussian has three independent parameters for the fit, that is the mean μ, the standard deviation σ, and the amplitude A. For two peaks, the sum of the gaussians has the following form Where the subscripts 1 and 2 indicate the two gaussians. For the one-gaussian fit, only the first term is considered (A 2 = 0). Table S1 summarizes all the fit values obtained. For the cases where an average is performed, reported as , the average is calculated by weighting the two averages μ i μ for the corresponding amplitude A i .

Partition coefficient calculation
As reported in the main text, the lipid/water partition coefficient can be calculated 2 as K ext where n is the total number of slices of the system, V(z i ) is the volume of the i th slice, ΔG(z i ) is the free energy at z i , R is the universal gas constant, T is the absolute temperature, n w (z i ) is the number of water molecules contained in the i th slice, and n w,tot is the total number of water molecules in the system. Lastly, x w (z i ) = n w (z i )/n w,tot with n w (z i ) being the number of water molecules contained in the i th slice and n w,tot being the total number of water molecules in the system. Thus, x w (z i ) is fraction of water molecules present in the slice i, and it's equal to the value at z = z i of the normalized probability distribution histogram for water.
in Eq. (2) is an extensive quantity, meaning that its value depends on the system size, and x w (z n ) (3) where N A is the Avogadro number, n lip is the total number of lipids, and M lip is the molecular weight of the lipids.
Since the simulations presented in this work have all rectangular boxes and the system is split along the z axis in slices of equal thickness, it follows that all the volumes V(z i ) are independent from i, that is, they have the same value δV. Then, Eq. (2), when converted to [L/kg] via Eq. (3) and by exploiting the symmetry , becomes (red open triangles). RMSE is calculated only for log K ow ≥ 1.

Fit of experimental log K lip data
As reported in the main text, the experimental SC lipid/water partition coefficients log K lip have been compiled together. The data comes from Wang et al. 5 (and references therein) and Ellison et al. 6 and is fitted against the following linear relationship   Figure S10 and Figure S11 show the residuals of the predictions plotted against the residuals of the log K ow , fitted against a line without intercept. In the main paper, the predicted log K lip values have been fit to the linear relationship in Eq.(5) both for the case of experimental and predicted log K ow . The results consist in a series of angular coefficients β for each system. In the following, angular coefficients obtained from fits against experimental log K ow ( ) are denoted as β exp , log K exp ow while those obtained from fits against COSMOtherm-predicted log K ow ( ) are denoted as log K pre ow β pre . The residuals (RES) are defined as  difference between predicted (pre) and experimental (exp) log K lip (blue open circles, first columns in Figure S10 and Figure S11)

Residuals analysis
 difference between predicted log K lip and fit against experimental log K ow (orange open squares, second columns in Figure S10 and Figure S11)  difference between predicted log K lip and fit against predicted log K ow (green open triangles, third columns in Figure S10 and Figure S11) log K pre lip -β pre ⋅ log K pre ow (8)  difference between predicted and experimental log K ow (the abscissa in all subplots in Figure S10 and Figure S11) column, green open triangles), as a function of the residual between predicted and experimental log K ow (as in Eq.(9)) for the hairpin systems (HPm2, HPm3, and HPm2t).

Summary of fit results
Figure S12 Summary of all the results for the linear fit between predicted log K lip and experimental (blue) and predicted (red) log K ow .