Skin permeability prediction with MD simulation sampling spatial and alchemical reaction coordinates

Abstract

A molecular-level understanding of skin permeation may rationalize and streamline product development, and improve quality and control, of transdermal and topical drug delivery systems. It may also facilitate toxicity and safety assessment of cosmetics and skin care products. Here, we present new molecular dynamics simulation approaches that make it possible to efficiently sample the free energy and local diffusion coefficient across the skin's barrier structure to predict skin permeability and the effects of chemical penetration enhancers. In particular, we introduce a new approach to use two-dimensional reaction coordinates in the accelerated weight histogram method, where we combine sampling along spatial coordinates with an alchemical perturbation virtual coordinate. We present predicted properties for 20 permeants, and demonstrate how our approach improves correlation with ex vivo/in vitro skin permeation data. For the compounds included in this study, the obtained log K_Pexp-calc mean square difference was 0.9 cm² h^-2.

Figure 1

A schematic representation of the structure of epidermis and dermis. To the right is shown a snapshot of the atomistic model of the intercellular lipid structure of human stratum corneum that constitutes the skin’s main permeability barrier (6). The arrow indicates the main permeability direction, also referred to as the Z dimension. In the MD simulations the system is in practice repeated infinitely in the three dimensions. The carbon atoms are colored based on the molecule type, where ceramides are green, acyl ceramides (ceramide EOS) are light blue, free fatty acids are orange, and cholesterols are yellow. Hydrogen, nitrogen, and oxygen atoms are colored white, blue, and red, respectively. To see this figure in color, go online.

Figure 2

Calculating the PMF of testosterone through the skin’s barrier structure. (a) On the right, the 2D free energy landscape of testosterone where the alchemical free energy λ state is shown on the y axis—0 is fully interacting and 20 is fully decoupled (as illustrated by the miniatures to the left of the free energy landscape). On the x axis the z coordinate (in nm) through the simulated system is shown. (b) The skin’s barrier system, aligned to illustrate how the PMF is mirrored around 0 (the middle of the system). (c) The PMF of testosterone through the skin’s barrier structure of the fully interacting λ state (0 on the y axis in a). The PMF is from one set of simulations of 24 communicating AWH walkers, for a total of 11 μs. AWH analyses do not give a reliable error estimate from one set of simulations, therefore there are no error bars presented in this plot. For error estimations based on multiple simulations see Fig. 3. To see this figure in color, go online.

Figure 3

Faster convergence with 2D AWH. The PMFs of testosterone are from 2D spatial/alchemical AWH (in the fully interacting alchemical λ state) and 1D AWH with only a spatial reaction coordinate dimension. In both cases the PMFs are symmetrized and only half the PMFs are presented, and the PMFs are calibrated to 0 at the ceramide sphingoid chain interface (5.2 nm) to make comparisons easier. The uncertainties represent one standard error of the mean between the independent sets of simulations. The 11 μs 2D AWH results are from one set of 24 communicating walkers. AWH analyses do not give a reliable error estimate from one set of simulations, therefore there is no error presented for the 11 μs AWH plot (in black). The 1D AWH simulations of 58 μs required approximately the same computation time as 22 μs of 2D AWH. A snapshot of the molecular system is shown below the plot to indicate where the head groups are located (at $\approx$ 3.1–3.3 nm). To see this figure in color, go online.

Figure 4

The correlation between experimental log K_P and those calculated from MD simulations. The blue line shows the linear regression, whereas the dashed black line is the identity line. The two largest outliers, urea and hydromorphone, are marked in red (still included in the linear regression). The error bars represent 1 mean ± SE (standard error of the mean), which is approximated for experimental values. Details about the approximation is presented in the header of Table S1. To see this figure in color, go online.