Understanding Drug Skin Permeation Enhancers Using Molecular Dynamics Simulations

Our skin constitutes an effective permeability barrier that protects the body from exogenous substances but concomitantly severely limits the number of pharmaceutical drugs that can be delivered transdermally. In topical formulation design, chemical permeation enhancers (PEs) are used to increase drug skin permeability. In vitro skin permeability experiments can measure net effects of PEs on transdermal drug transport, but they cannot explain the molecular mechanisms of interactions between drugs, permeation enhancers, and skin structure, which limits the possibility to rationally design better new drug formulations. Here we investigate the effect of the PEs water, lauric acid, geraniol, stearic acid, thymol, ethanol, oleic acid, and eucalyptol on the transdermal transport of metronidazole, caffeine, and naproxen. We use atomistic molecular dynamics (MD) simulations in combination with developed molecular models to calculate the free energy difference between 11 PE-containing formulations and the skin’s barrier structure. We then utilize the results to calculate the final concentration of PEs in skin. We obtain an RMSE of 0.58 log units for calculated partition coefficients from water into the barrier structure. We then use the modified PE-containing barrier structure to calculate the PEs’ permeability enhancement ratios (ERs) on transdermal metronidazole, caffeine, and naproxen transport and compare with the results obtained from in vitro experiments. We show that MD simulations are able to reproduce rankings based on ERs. However, strict quantitative correlation with experimental data needs further refinement, which is complicated by significant deviations between different measurements. Finally, we propose a model for how to use calculations of the potential of mean force of drugs across the skin’s barrier structure in a topical formulation design.


■ INTRODUCTION
Transdermal drug delivery enables continuous, noninvasive, and pain-free drug administration that, compared to oral or intravenous drug delivery, is associated with reduced side effects and increased patient compliance. Continuous drug administration for days to weeks, combined with the avoidance of first-pass metabolism in the gastrointestinal tract and liver, would give physicians greater control over the medication and dosage levels. The human skin is, however, an effective barrier to permeation, limiting the number of drugs presently available for transdermal delivery.
The skin's main permeability barrier is located to a stacked lamellar lipid structure 1,2 situated in the intercellular space of the superficialmost part of the skin, the stratum corneum. 3 This lipid structure consists of a unique mixture of ceramides, cholesterol, and free fatty acids 4−6 with a water content of approximately one to two molecules per lipid. 7,8 A number of different models have been proposed for its detailed lipid packing. 8−15 Today, the dominant means of assessing skin permeability is by ex vivo/in vitro testing in diffusion cells. 16 These techniques often use excised skin, separating the donor and receptor chamber. To enhance transdermal drug transport, formulations containing one or more chemical permeation enhancers (PEs) are commonplace. 17 However, the modes of action of PEs are often unknown, which limits their use in commercial applications. Various PE action mechanisms have been proposed (cf. ref 17), including (i) altering the solubility of the drug in the formulation, (ii) modifying the partitioning of the drug into the skin's barrier structure, (iii) changing the thermodynamic activity of the drug, (iv) promoting transport by "dragging" the drug across skin, and (v) perturbing the lipid organization of the skin's barrier structure. In addition, many PEs have complex and nonlinear concentration-dependent permeability-enhancing effects. In order to understand these effects, the partitioning of various formulation excipients into the stratum corneum has been investigated. 18−20 Furthermore, the effect of PEs on the structural properties of the stratum corneum has been studied using an array of different experimental methods, including differential scanning calorim-etry, infrared spectroscopy, small-angle X-ray diffraction, and solid-state NMR spectroscopy. 21 −23 In silico modeling can offer additional insight into the absorption of chemicals into and permeation of chemicals across skin. A detailed understanding of how PEs interact with the skin's barrier structure could enable a tailored design of topical drug delivery vehicles, ensuring that they modify the barrier in regions that present the highest permeation resistance to a specific drug. Quantitative structure−permeability relationship (QSPR) 24−28 methods are often used to predict permeability coefficients for molecules permeating the skin. QSPR models are good at predicting drug permeability coefficients, with reported mean square errors as low as 0.48 cm 2 h −2 (log units) and R 2 = 0.7. 25,28 However, due to the wide range of molecular mechanisms for the actual permeation process, QSPR methods are limited with respect to faithfully modeling the permeability-enhancing effects of PEs, as the methods will need to have been trained on PEs that have both similar structures and similar permeation-enhancing molecular effects. The general prediction performance of QSPR methods is thus reduced for new categories of PEs, and the methods are unable to explain the mechanism of action of PEs.
Molecular dynamics (MD) simulations are not as domainlimited as QSPR methods. In addition, they enable the incorporation of PEs into the skin's barrier structure at their most probable concentrations, locations, and orientations. The atomistic resolution of MD simulations can yield an understanding of how PEs and other formulation excipients interact with the barrier in order to modify the permeation of various drugs as well as potential direct interactions with the barrier lipids. It should be noted that the atomistic models are relatively simplistic, and detailed knowledge about the exact lipid organization in the human skin barrier is still evolving. Furthermore, the currently used models do not contain any representation of corneocytes or the cornified cell envelope. However, the ability of MD simulations to address how various formulation excipients influence skin permeability makes them an attractive alternative or complement to QSPR methods with respect to transdermal drug delivery design. MD simulation modeling of skin lipid bilayers has been performed by several groups, 29−35 often based on hairpin bilayers consisting of ceramides, cholesterol, and free fatty acids. More complex structures incorporating ceramide EOS have also been employed, 36,37 yielding better agreement with experimentally measured permeability values.
Recently, we have presented an improved simulation protocol for transdermal drug transport modeling, which enables improved precision through better sampling of the free energy landscape as well as a better calculation of the local diffusion coefficient of drugs permeating the skin's barrier structure. 38 Unfortunately, good reference in vitro skin permeation data are scarce, largely due to variability of reported values between laboratories (e.g., roughly 1 order of magnitude for testosterone 39,40 ) and/or to differences in experimental setups (e.g., in donor concentration, temperature, skin thickness, donor/receptor fluids) which impede direct comparison. In this work, we selected two recent studies for testing our modeling approach. These studies were chosen to avoid comparing permeability data on the same permeant from different laboratories and to ensure that the reported permeability enhancement was obtained using an identical experimental setup. First we compare our results with those of Pham et al., 23 who reported the effect of saturated solutions of single PEs on skin permeation of metronidazole. Then we focus on selected data from the work of Abd et al., 41 which involves the permeability-enhancing effect on caffeine and naproxen of three formulations based on ethanol and water. This was assumed to pose a greater challenge due to the differing ionization of naproxen in the different formulations but enabled the possibility to compare the effect of similar formulations on different permeants. We also decided to not include comparisons with sodium lauryl sulfate due to expected difficulties of modeling charged molecules. Finally, we outline our basic general approach with respect to MD skin permeability modeling and address the challenges ahead.
For more detailed information about the simulation methods used and the atomistic skin barrier model employed, the reader is referred to previous publications. 8,38,42 ■ METHODS The atomistic model of the human skin's barrier structure presented by Lundborg et al., 8 and originally called 33/33/33/ 75/5/0.3, a was used as the starting structure for the simulations. The starting model had been equilibrated for approximately 270 ns, 250 ns of which was without restraints. 8 All MD simulations in this study were performed using GROMACS 2022, 43−45 using the accelerated weight histogram (AWH) method for alchemical free energy calculations. 42,46,47 A source code modification enabled symmetrizing AWH sampling along a spatial (pulling) reaction coordinate dimension and also changed the number of AWH blocks for autocorrelation analysis from 64 to 128. These changes are available from the GROMACS gitlab repository 48 and are expected to be included in future releases of the code. Symmetrizing effectively improves the convergence of the simulations. We have previously not experienced any significant difference in the results when comparing PMFs calculated with and without symmetrization (see Figure S1 in ref 49).
Van der Waals interactions had a cutoff of 1.2 nm with a smooth force switch from 1.0 to 1.2 nm. The simulations were run without a dispersion correction for the energy and pressure. Coulomb interactions were calculated using PME 50 with a radius of 1.2 nm. Bonds to hydrogen atoms were constrained using the P-LINCS algorithm. 51,52 TIP3P 53 parameters were used for water molecules. For the lipid molecules, the CHARMM36 lipid force field 54,55 was used. Ceramide parameters were modified to more accurately reproduce the ceramide NP crystal structure 56 as described in ref 8. In order to allow a 3 fs integration time step, hydrogen atoms were made 3 times heavier by repartitioning the corresponding mass from their bound heavy atoms. 57,58 The temperature was set to 305.15 K by using a stochastic dynamics integrator (also referred to as a velocity Langevin dynamics integrator) with a time step of 3 fs and a time constant (τ) of 2 ps (corresponding to a friction constant of 0.5 ps −1 ). The pressure was set to 1 atm and controlled using a stochastic cell rescaling barostat 59 with a time constant of 1.0 ps and a compressibility of 4.5 × 10 −5 bar −1 . In all simulations containing the skin barrier system, the pressure coupling was semi-isotropic with no compressibility in the Z dimension to maintain the lateral spacing derived from the original CEMOVIS data.
Along the alchemical free energy dimension, 21 equidistantly distributed λ states were used for decoupling both van der Waals and Coulomb interactions simultaneously. There were 100 steps between each Monte Carlo coordinate sampling along the alchemical free energy λ value and spatial position across the lipid structure. There were 10 samples per update of the f λ bias. The estimated initial error was set to 10 kJ mol −1 . Soft-core transformations 60 with α = 0.5 and σ = 0.3 nm were applied to both the van der Waals and Coulomb interactions of the solute.
Topologies (i.e., inter-and intramolecular interaction parameters) for all permeants and formulation components except water were generated using STaGE, 61 which in turn uses Open Babel 62 and MATCH 63 to generate GROMACS topologies compatible with the CGenFF 64 CHARMM force field.
All images representing molecules were prepared using Tachyon 65 in VMD. 66 Solvation Free Energy Calculations. The solvation free energy was calculated starting with a box large enough to solvate the solute (system sizes for all formulations simulated in this work are shown in Table S5). For each solvation free energy calculation, the solute was inserted into the system at a random position with its interactions to the surroundings turned off, as if in vacuum. The temperature was set to 305.15 K. AWH was used to sample the alchemical free energy λ states 42 to calculate the solvation free energy. The initial AWH histogram size, which determines the initial update size of the free energy, was set indirectly by specifying an estimate of the diffusion constant along the alchemical free energy λ axis in combination with a rough estimate of the initial error. An input diffusion constant of 1 × 10 −3 ps −1 was used for the solvation free energy calculations, which means that it is estimated to take approximately 1 ns to cross the alchemical dimension for one AWH walker. Twenty-four communicating AWH walkers were run in parallel, with the requirement that only simulations that covered the whole alchemical reaction coordinate counted toward the covering check in the initial AWH stage. The simulations were 60 ns long per walker for a total simulation time of 1440 ns per solute.
Permeability Coefficient Calculations. The permeability coefficient, K P , and the permeation resistance, R, were calculated from the obtained potentials of mean force (PMFs) across the skin's barrier structure as follows: 67 where ΔG transfer is the free energy relative to the solvation free energy, β is the thermodynamic β, and D is the local diffusion coefficient across the skin's barrier structure (dz). The PMFs were shifted so that the PMF minimum was never below 0, and following the discussion in ref 38, the factor 30 ± 6 represents an estimate of the number of lipid bilayers present in the stratum corneum based on observations from cryo-EM images of human skin. 14,68 Since the data presented herein mostly focus on relative differences between simulated and in vitro enhancement ratios, the number of assumed bilayers does not influence the final conclusions. The AWH simulation estimates a friction metric tensor g(z), 69 which can be used to derive the diffusion coefficient across the barrier structure, D(z) = g −1 (z). To reduce the noise of the local diffusion coefficient curves, a 0.2 nm wide rolling median filter was applied. When symmetrizing the sampling along the spatial dimension by using the absolute coordinate values and accounting for the AWH bias across the sampling boundaries, there are usually artificial spikes at the edges of the PMFs. These were removed by setting the two lowest points (0.005 and 0.015 nm) and two highest points (5.200 and 5.210 nm) in the PMF to the values of their neighbors (0.025 and 5.19 nm, respectively). These minor adjustments were made merely to avoid strange-looking features in the curves and are far smaller than any statistically significant effects on the calculated permeability coefficients.
At the start of each AWH walker simulation, the permeating molecule was inserted in the lipid barrier structure at a random position with all interactions to its environment turned off. The free energy profile through the skin's barrier structure was calculated with a two-dimensional AWH setup using a harmonic potential to steer the permeant along the direction normal to the barrier plane (also referred to as the Z dimension) and an alchemical free energy reaction coordinate. 42 This allows sampling of the free energy along the permeation direction and also the relative insertion free energy of the permeant from vacuum. In turn, this enables a direct calibration to the solvation free energy since the vacuum state is the same in both cases. The AWH input diffusion constant was set to 3 × 10 −5 nm 2 ps −1 for the spatial pulling dimension and 5 × 10 −5 ps −1 along the alchemical free energy dimension. The input diffusion constant affects only the AWH histogram size but not the computed diffusion coefficient, which is obtained from the AWH friction metric during data analysis. The AWH force constant along the spatial Z dimension (normal to the lamellar stack), steering the permeant relative to the ceramide fatty acid chains using a harmonic pull potential, was set to 25 000 kJ mol −1 nm −2 . The force constant also determines the resolution along the reaction coordinate dimension. For each permeant, three sets of simulations were run with heavy hydrogen atoms (see above) and a 3 fs integration time step. These were run using 24 communicating walkers, each running for 450 ns. The combined PMF was derived from the average of the independent PMFs from the three sets of simulations.
Along the alchemical free energy dimension, it is the end states, i.e., the fully interacting and fully decoupled states, that are of highest interest, as the difference in free energy between them corresponds to the probability of transferring the permeant from vacuum into the skin's barrier structure. Therefore, the target distribution used in these simulations puts more weight on the end states, especially the state with interactions fully turned on. The target distribution along the alchemical free energy dimension that was used is the same as shown in the supplement of ref 38. Along the spatial pulling dimension, the target distribution was uniform.
Skin Partition. From the solvation free energy and the PMF we calculate the partition coefficient of the solute into the skin barrier structure as follows: where ΔG rel_water is the Boltzmann-weighted free energy in the lipid system relative to water, R is the gas constant, and T the temperature. Insertion of Permeation Enhancers. Due to the limited knowledge about the structure of the surface at the point where the applied formulations contact the skin as well as the associated large computational cost, it is unfeasible to perform a direct simulation of the partitioning of PEs into the skin's barrier structure. Instead, we utilize the PMFs obtained for each PE through the reference bilayer structure (without any added excipients) to calculate the total number of PEs to insert from the formulation into the skin system, N mol , as follows: where ΔG transfer is the free energy relative to the solvation free energy in the formulation, β is the thermodynamic β, and the constant C is the number of PEs in the formulation corresponding to the same volume as the simulated skin bilayer. In order to avoid inserting the PEs in regions where they have a low probability of partitioning into, we divided the system into 20 slices along the Z direction (normal to the bilayer) and performed the calculation in eq 3 within each slice. The PEs were then inserted at random positions in the lateral XY plane within their corresponding slice with their interactions with their surroundings turned off. The system was then simulated for 5 ns, during which the PEs had their electrostatic and Lennard-Jones interactions linearly increased until they were fully interacting with the surrounding molecules at the end of the simulation. A position restraint along the Z dimension with a force constant of 1000 kJ mol −1 nm −2 was applied to the PE center of mass in order to prevent the PE from leaving its starting position during the early stages of the simulation. In the systems where multiple PEs are present, this procedure was performed sequentially for each type of PE. After all PEs had been inserted, the systems were equilibrated without restraints for 300 ns. This procedure for the insertion of PEs into the skin bilayer structure is summarized in Figure 1. Deuterium Order Parameters. The deuterium order parameters were calculated using the gmx order command from the GROMACS simulation package, with a selection of the ceramide sphingoid and fatty acid chain carbons, in separate command executions, that exist in all ceramides (fatty acids C2−C20 and sphingoid C4−C18). The first 50 ns of the simulation trajectories was discarded as equilibration.

■ RESULTS AND DISCUSSION
Modeling the permeability-enhancing effect of PEs is challenging, partly due to the complex nature of skin. Our approach is based on using MD simulations to calculate the free energy difference between a topical formulation and the skin's barrier environment (Figure 1). This enables us to model the partitioning of formulation components into the barrier structure and then calculate their permeabilityenhancing effect on various pharmaceuticals. Here, three major constraints need to be fulfilled: (1) the calculated permeability coefficients of a permeant in the absence of PEs should correlate with in vitro/in vivo measurements; (2) the molecular force field used in the MD simulations must be able to reproduce the partitioning behavior of all excipients, including the permeant and PEs, between the formulation and the skin's barrier structure; (3) the calculated permeability enhancements should correlate with in vitro/in vivo measurements of a wide range of PEs and permeants.
We have previously shown that our model of the skin's barrier structure, combined with the employed MD simulation techniques, is able to satisfy the first of these constraints, i.e., that the calculated permeability coefficients of single permeants correlate with in vitro/in vivo measurements. 38 To address the second constraint, i.e., to verify that our chosen force field is able to adequately model the partitioning behavior of formulation excipients, we have now calculated the correlation between experimental and calculated excipient skin partition coefficients. Finally, as a first step toward addressing the third constraint, i.e., the modeling of complex topical formulations including several PEs, we here chose to study saturated solutions of a single PE (using in vitro data from Pham et al. 23 ) as well as ethanol/water formulations containing a single additional PE in low concentration (with data from Abd et al. 41 ).  Table 1 as an example. (a) The solvation free energy, ΔG solvation , of the PE in the applied formulation is calculated relative to a vacuum. (b, c) The free energy of insertion, ΔG lipid (with vacuum as the reference state), as well as the pulling free energy through the skin's barrier structure are calculated using the two-dimensional accelerated weight histogram method. (d) By combining the results obtained in (a−c), the free energy difference between the formulation and the different parts of the skin's barrier structure, ΔG transfer , can be calculated. This free energy is then used to calculate the number of molecules to insert along the Z dimension of the system according to eq 3. (e) The final system is equilibrated and can then be used for permeability calculations of selected permeants. Carbon atoms are colored as follows: in nonacyl ceramides, green; in acyl ceramide EOS, light blue; in free fatty acids, orange; in cholesterol, yellow; in geraniol, teal. The rest of the atoms are colored according to atom type (oxygen, red; nitrogen, blue; hydrogen, white). The illustrated skin patch in the image is adopted with permission from refs 38 (CC BY 4.0) and 70 (copyright 2016 C. Wennberg).

Calculated Partition Coefficients.
To assess the quality of the automated generation of small-molecule parameters with STaGE 61 (and MATCH 63 ), the octanol−water partition coefficients, log K ow , of 16 different molecules were calculated by combining the results from three different sets of solvation free energy calculations in water and octanol for each solute: ln (10) ow ow (4) where ΔG ow is the difference in average solvation free energy in water compared to octanol, R is the gas constant, and T the temperature. The correlation with experimental values (Table S1 and Figure S1) gave an RMSE of 1.11 log units with R 2 = 0.78, similar to previously reported log K ow values for the CHARMM36 force field. 71 Following this, the partition coefficient from water into the skin's barrier structure, log K w−lip was calculated for all 16 molecules using eq 2, and the correlation with experimentally measured partition coefficients is displayed in Table S1 and Figure 2. The obtained RMSE was 0.58 log units with R 2 = 0.87, indicating that the generated parameters seemed to adequately predict the partitioning behavior of small molecules into the skin's barrier structure, with most of the outliers positioned around log K w−lip < 1. The experimental relationship between log K w−lip and log K ow has previously been observed as 72 = K K log 0.69 log w lip ow (5) which fits the in vitro data used here as well ( Figure S2). For our calculated partition coefficients, the correlation is best fitted using a slope of 0.89.

Construction of Skin Barrier Model Systems with Inserted Permeation
Enhancers. The skin barrier model systems with inserted PEs are shown in Figure 3, and the calculated deuterium order parameters for the ceramides in the skin's barrier structure are shown in Figure S3. We also calculated the number of hydrogen bonds between all the barrier lipids in each system; these are shown in Table S7. All excipients lower the number of hydrogen bonds compared to the pure skin system by approximately 15−30%. We have not observed any clear correlation between the amount of broken hydrogen bonds and permeability enhancement ratios for the permeants studied in this work.
The potential of mean force for each PE across the skin's barrier structure is shown in Figure S4. The position along the lipid bilayer normal at which each PE is inserted is primarily decided by the lowest points along the PMF. Water and ethanol prefer to insert in the ceramide headgroup regions of the bilayer structure ( Figure S4, z ∼ 3.2 nm). Eucalyptol, geraniol, and thymol prefer insertion close to the interface between either the sphingoid or fatty acid chains of the ceramides ( Figure S4, z < 1 nm or z > 4.5 nm). Lauric acid, oleic acid, and stearic acid prefer to be inserted in the sphingoid or fatty acid chain regions of the ceramides ( Figure  S4, 1.5 < z < 2.5 nm or 4 < z < 4.5 nm).
Since the number of molecules to insert scales exponentially with the free energy difference between the formulation and the most favorable position in the lipid bilayer structure, we have obtained very high insertion numbers for some systems. In order to limit this effect, we currently use the calculated number of waters inserted from a pure water solution (see Table S2) as a guideline for the upper bound of the total mass of PEs to insert in the headgroup region. The calculated amount of water to insert (100 molecules) results in a total water content (190 molecules) that is roughly 2.1 times the water mass already present in the skin model (90 molecules), and we currently limit the insertion such that the total mass of PEs (including inserted additional water) in the headgroup region does not go above 2.5 times the mass of water already present in the system. The effect of this limit on the total number of molecules inserted in each simulated system is shown in Table S2 (note that only water and ethanol partition into the headgroup region�the other PEs are not affected). For the sphingoid and fatty acid tail regions, we have not found a suitable reference for an upper limit of insertion, and therefore, we do not limit the insertion in these regions.
Metronidazole Skin Permeation Calculations. The calculated permeabilities for metronidazole, as well as comparisons with the experimental values obtained by Pham et al., 23 are displayed in Table 1. The skin barrier systems with inserted PEs are shown in the left column of Figure 3 (panels a−c), and the calculated permeation resistance profiles across the bilayer structure are shown in Figures S5−S7.
Although calculated values for individual permeation enhancers in some cases deviate significantly from the corresponding measurements, the relative enhancement ratio (ER) compared reasonably well, with a slight underestimation of the absolute permeation enhancement. For stearic acid, the concentration in the formulation was so low (Table S2) that a negligible amount was predicted to partition into the skin's barrier system during the permeation measurement of metronidazole. The calculated crossing time, from a pure water solution, for stearic acid was also 4 orders of magnitude longer than that of metronidazole. Thus, the system was assumed to be identical to that of the pure water phase since the stearic acid should have almost no influence on the solubility in the formulation or the permeation across the stratum corneum.
Lauric acid displayed a relatively large uncertainty in its PMF at the lowest points (see Figure S4   . Skin partition coefficients. Correlation between calculated and experimentally measured partition coefficients from water into the barrier structure. Dashed lines above and below the identity line correspond to ±0.5 log unit. Data are also available in Table S1. Experimental data were taken from refs 18−20.

Journal of Chemical Information and Modeling pubs.acs.org/jcim
Article on the PMF, 16 molecules were to be inserted into the bilayer structure, but this resulted in almost no difference in observed permeability, in disagreement with the experimental data.
Therefore, we decided to increase the amount of inserted lauric acid molecules in a stepwise manner, increasing the amount by 50% at each step. Snapshots of the equilibrated systems are shown in Figure S8. From the series of simulations with increasing amounts of lauric acid, we observe a fast increase of the permeability enhancement already when inserting 50% additional lauric acid (24 inserted), as shown in Table 1. When the amount is increased further, no statistically significant extra enhancement is observed (see Table S3). From the reference permeation resistance profile for metronidazole (e.g., black curve in Figure S5), we can see the that the main permeability barrier is situated in the ceramide sphingoid chain region (4 nm < z < 5 nm), with the headgroup region (z ∼ 3.2 nm) posing the least resistance followed by a second barrier in the ceramide fatty acid chain region (z ∼ 2.1 nm). When inserting lauric acid into the skin bilayer ( Figure S5), we see a reduction of the barrier in the ceramide sphingoid chain region of the bilayer stucture with increasing lauric acid content. The barrier in the ceramide fatty acid chain region increases with 16 inserted lauric acids but is otherwise relatively unchanged compared to the reference.
For geraniol we found multiple values reported in the literature regarding its solubility in pure water: 100 mg L −1 (at 298 K), 73,74 229 mg L −1 (no temperature reported), 75 and 686 mg L −1 (at 293 K). 76 For completeness, we decided to calculate the effect of three different formulations of geraniol: 100, 200, and 686 mg L −1 (all permeability values are reported in Table S3). The solubility series reported recently by Martins et al. 74 indicates that the most likely correct water solubility of geraniol at a temperature of 305 K would be roughly 200 mg L −1 , which is the formulation used to calculate the data we report in Table 1. From Table S3 we see that the permeability enhancement mimics the solubilities, with the highest solubility having the highest permeation enhancement. The geraniol molecules prefer to aggregate in the fatty acid chain region of the skin's barrier structure ( Figure S9) but there is no apparent effect on the permeation resistance in this region for the two lower-solubility calculations ( Figure S6). Instead, there is a clear reduction of the barrier capacity in the ceramide sphingoid chain region, most likely related to the decreased lipid ordering in this part of the system ( Figure S3a).
For thymol, when using the saturated concentration of 900 mg L −1 , the partitioning calculations resulted in a local concentration (in the most populated insertion slice in the lipid bilayer structure) of thymol that would be higher than that of a pure thymol solution, and the resulting amount of inserted thymol was roughly equal to the mass of the whole bilayer (see Table S2 and Figure S10). Since this concentration would be unrealistic, we also performed calculations with lower concentrations of thymol in the skin's barrier structure. The concentration of pure liquid thymol is 3.9 molecules nm −3 (density 0.97 g cm −3 ) and we constructed two systems in which the highest local concentration was either 2.3 molecules nm −3 or 4.7 molecules nm −3 . These local concentrations corresponded to a calculated thymol concentration in the formulation of either 100 or 200 mg L −1 , which is how they are denoted in Table S2, S3, S5 and Figures S3, S7, S10. The data presented in Table 1 are based on the system corresponding to a thymol concentration of 200 mg L −1 in the formulation. In the calculated permeation resistance profiles ( Figure S7) we see a large reduction of the permeability barrier for all the simulated thymol systems, with an almost complete removal of the barrier for the system with a thymol concentration of 900 mg L −1 in the formulation (as expected given the unrealistic thymol concentration in the system). The in vitro measurement never reached steady state, indicating a very high  permeability enhancement for thymol, which is also reflected in the calculated values. Caffeine and Naproxen Skin Permeation Calculations. The calculated permeabilities for caffeine and naproxen, as well as comparisons with the experimental values obtained by Abd et al., 41 are displayed in Table 2. The skin barrier model systems with inserted PEs are shown in the right column of Figure 3 (panels d−f), and the calculated permeation resistance profiles across the lipid bilayer structure are shown in Figure S11. When calculating the partitioning from the formulations into the skin's barrier structure, there is a need to account for ethanol evaporation in the formulations. A correct estimate of ethanol evaporation during formulation preparation, subsequent formulation skin application, and final permeability measurement (even within covered donor chambers some evaporation occurs) is difficult to assess. 77,78 For these calculations, we have assumed an ethanol evaporation of 50% from the formulations, but a more indepth evaluation is envisaged for future studies. Furthermore, since ionization of naproxen (pK a = 4.19) varies between simulated formulations (pH between 2.9 and 4.8 41 ), the naproxen permeability coefficients were shifted based on the formulation pH according to pH p a aq (6) where K P e is the effective permeability coefficient, K P 0 is the calculated permeability coefficient, and pK a aq is the aqueous pK a of naproxen. This shift results in a better agreement in absolute permeation coefficients as well as in in vitro enhancement ratios for naproxen (the unmodified and pH-corrected permeability coefficients can be seen in Table S4). Similarly to metronidazole, the reference resistance profile for caffeine (black curve in Figure S11) has its main permeability barrier in the ceramide sphingoid chain region (z ∼ 4.8 nm), with a minimum in the headgroup region followed by a second barrier in the ceramide fatty acid chain region (1 nm < z < 2 nm). For naproxen (note that the reference simulation is from an ethanol/water formulation), a major permeability barrier can be seen in the ceramide sphingoid chain region around z ∼ 4.8 nm, similar to caffeine. Contrary to caffeine, naproxen then experiences a second barrier close to the headgroup region (z ∼ 3.5 nm) and possibly a minor third barrier around z ∼ 1.3 nm in the ceramide fatty acid region.
The addition of 60% ethanol to a water solution lowers the main permeability barrier of caffeine in the ceramide sphingoid chain region ( Figure S11) but has almost no effect in other parts of the system. For both caffeine and naproxen we see an increased permeation resistance in the fatty acid chain region (z ∼ 1.8 for caffeine and 1.5 < z < 3 nm for naproxen) in the systems with inserted eucalyptol, while only a small or no reduction of the barrier capacity is seen in the ceramide sphingoid chain region. Since the permeability enhancement of eucalyptol was much lower than expected, we also performed calculations with twice the amount of eucalyptol inserted ( Figure S12), which in both cases resulted in a higher permeability enhancement and thus better agreement with experimental data.
The calculated permeability of naproxen through the systems with inserted oleic acid displayed a very low increase in permeability (similar to metronidazole with lauric acid as a PE) compared to the in vitro data. Therefore, we performed a second iteration of the skin partitioning calculations, using the first constructed version of the system ( Figure S13a) to recalculate the PMF of oleic acid. The PMFs from both iterations are shown in Figure S14, and the final inserted amount of oleic acid was 54 molecules (Figure S13c). In a similar fashion to the simulations with lauric acid, we decided to increase the amount of inserted oleic acid in steps, going from 18 to 36 and then finally to 54 inserted molecules. The permeation resistance profiles ( Figure S11) of caffeine and naproxen calculated across the different systems with oleic acid show a decrease of the permeation resistance in the ceramide sphingoid chain region (4.7 nm < z < 5.2 nm) for all systems. The permeation resistance in the ceramide fatty acid chain region for caffeine appears relatively unchanged, except for the system with 18 inserted oleic acid molecules, in which the barrier capacity decreases slightly. For naproxen the barriers close to the ceramide headgroup region and in the ceramide fatty acid chain region both increase with the insertion of 18 oleic acid molecules, but for 36 and 54 inserted molecules the barrier in the ceramide headgroup region is reduced and the barrier in the ceramide fatty acid chain region is relatively unchanged.
For the calculated permeability values, starting with caffeine, we observe a relatively good qualitative agreement between the in vitro measurements and the calculated ERs. However, the large increase in the in vitro permeability enhancement from both oleic acid and eucalyptol is greatly reduced in the calculated data. It should be noted that the caffeine permeability coefficient reported by ) is lower than those reported elsewhere (e.g., −3.6 79 and −3.2 80 ). Unfortunately, caffeine has been reported to be one of the more problematic small molecules to parametrize correctly in existing force fields, 81,82 and the CHARMM36 force field parameters used in these simulations overestimate the log K ow of caffeine by 2 log units (see Table S1). This might be an explanation for the poor permeation enhancement obtained in the simulations since a lipophilic molecule experiences lower permeation barriers across the skin's lipid bilayer structure. Furthermore, our obtained standard errors in the calculated permeability coefficients for caffeine are so large that they cannot be deemed statistically different, just higher than that of the reference.
For naproxen, we see relatively good agreement between the calculated permeability coefficients and ERs compared to the in vitro data. When compared to the ERs based on in vitro permeability coefficients, the calculated ERs for oleic acid and eucalyptol are on the lower side although still within the same order of magnitude (and relative difference) as the in vitro data. The ERs based on in vitro maximum flux compare very well with the calculated data with an additional amount of inserted permeability enhancers. This trend is consistent for both PEs, i.e., we see a low permeability enhancement using the amount of inserted molecules suggested by eq 3, but the simulations with an increased amount have a better correlation with in vitro data. This indicates that the current methodology might underestimate the amount of PEs that should be inserted into the system.
Flux Optimization Using the Potential of Mean Force. The flux of a permeant across the skin's barrier, J, is calculated as Journal of Chemical Information and Modeling pubs.acs.org/jcim Article where C v is the permeant concentration in the formulation and K P is the permeability coefficient. Since the PMF across the skin's barrier structure is calibrated against the solvation free energy of the permeant in the formulation (ΔG transfer (z) in eq 1), it should be possible to use the minimum in the calculated PMFs to optimize the permeant flux from a specific formulation. For optimal flux, the permeant free energy difference between a formulation and the skin's barrier structure (i.e., the difference between the solvation free energy in the formulation and the minimum of the PMF) should be close to zero, as this enables a high solubility of the permeant in the formulation (thus maximizing its concentration) while maintaining a good partitioning of the permeant into the skin's barrier structure. Thus, depending on where along the y axis (see Figure S15) the PMF minimum of a permeant is located, we would argue that two different methods to modify the permeant's flux across the barrier exist. In the first case, when the PMF minimum is below zero, the solubility of the permeant in the formulation can be increased without modifying the K P (see eq 1 and the paragraph just after the equation). In such cases, one could modify the formulation using inert excipients that have no impact on the skin's barrier structure (e.g., long-chain PEG or Brij 98) in order to increase the concentration of the permeant in the formulation and thus increase its flux across the barrier. In the second case, when the PMF minimum is above zero, an increased permeant solubility in the formulation would increase the PMF and thus result in a reduced K P . In this situation, there would be no positive effect on the permeant flux by increasing the permeant's solubility in the formulation, since an increased formulation solubility would result in a similar increase in the relative permeation barriers across the skin's barrier structure. Instead, focus should be on modifying the formulation using PEs that reduce the permeant's permeation barriers across the barrier structure.
As an example, naproxen has a reported maximal flux of 0.6 μg cm −2 h −1 from water, 83 and its PMF minimum from a water solution in Figure S15 is −5.4 ± 1.0 kJ mol −1 . If one were to add, e.g., Brij 98 to the water solution (increasing the naproxen solubility and thus maximizing its concentration), raising the PMF minimum to zero, one would expect the maximum flux to be ∼10 times higher, which is similar to what Abd et al. reported for their measurements of naproxen permeability when varying the solubility using different inert control vehicles. 41 For their PEG6000/water formulations, they saw a plateau of the permeability enhancement at a solubility of 45 mg cm −3 , at which they reached a maximal flux of roughly 7 μg cm −2 h −1 , thus obtaining a maximal enhancement ratio (ER) of approximately 10 compared to naproxen from pure water. For caffeine their data show a similar concentration-dependent behavior of the flux (an initial increase from 2.2 μg cm −2 h −1 until a plateau in flux values is reached around 6.6 μg cm −2 h −1 ), 41 but our calculated value for the PMF minimum is roughly −13 ± 1.0 kJ mol −1 , which is a clear overestimation (corresponding to a possibility to increase the flux roughly 180 times) of the partitioning compared to the in vitro data. Given that the force field parameters for caffeine appear to be slightly too lipophilic (log K ow in Table S1), we would expect the real partitioning free energy to be much higher.
The authors are aware of the fact that this idea goes against conventional theory regarding permeability from saturated solutions (that maximum flux is independent of the solubility in saturated inert solutions that do not affect the skin barrier) but still believe that it would be of interest to test this in controlled steady-state skin permeability experiments in the future.

CONCLUSIONS
Understanding how chemical permeation enhancers (PEs) impact skin to improve transdermal delivery of drugs has great implications for clinical medicine, as replacing oral or intravenous drug administration with topical drug administration would give physicians better control over medication and dosage levels. Due to the inability of current in vitro methods to predict the effects of PEs on transdermal drug transport, complementary in silico methods are requested. We have presented our current approach using molecular dynamics (MD) simulations to predict the effect of PEs on drug permeation through skin. The atomistic resolution of MD simulations allows a detailed description of the permeability barriers in skin and enables investigation of how PEs modify these barriers at different locations across the skin's barrier structure. Qualitative agreement between calculated and in vitro-determined permeability coefficients is good, and MD simulations are able to reproduce permeability enhancement ratio (ER) rankings. This indicates that the fundamental models, simulations, and sampling algorithms used are sound, and with smaller statistical errors we are able to identify both outliers and systematic limitations of the models.
Still, the quantiative correlation of the calculated enhancement ratios with in vitro data needs further refinement, with future work focusing on modeling formulation excipient evaporation, permeant ionization in relation to pH, permeation speed of PEs in relation to permeant permeation speed, and skin partitioning of formulation excipients. Some excipients or permeants, such as caffeine, may be difficult to model using the CHARMM36 molecular force field. However, a preliminary screen of molecular properties (e.g., octanol−water partition coefficients) might give an indication of such cases.
Calculated permeant potentials of mean force across the skin's barrier structure can be used to determine whether modification of a permeant's permeation barriers in skin with PEs or modification of its solubility in the formulation will be most efficient to increase the flux of a permeant across skin.
Future refinement of the calculation of formulation excipient partitioning into skin might be called for, as our current setup appears to underestimate the excipient amounts inserted. Improvements regarding the incorporation of long-chain fatty acids, or other lipids, into the skin's barrier structure may also be requested. Although the use of the pH − pK a difference to shift permeability coefficients appears to work well for charged molecules, a correct treatment of ionization states and salts might require a more advanced simulation protocol such as constant-pH simulations. 84

■ DATA AND SOFTWARE AVAILABILITY
The input data and parameters are available for download at 10.5281/zenodo.7620240. coefficients; molecular compositions of simulated formulations; calculated octanol−water partition coefficients; calculated deuterium order parameters; potentials of mean force for each permeability enhancer; potentials of mean force and permeation resistance profiles for naproxen, caffeine, and metronidazole; and skin barrier structures after insertion of permeation enhancers (PDF) ■ ■ ADDITIONAL NOTE a Relative composition in mol % ceramides/mol % cholesterol/ mol % free fatty acids/relative amount of cholesterol on ceramide sphingoid side/mol % acyl ceramide EOS (included in the relative ceramide concentration)/water molecules per lipid (not included in the mol % concentrations of the lipids).