1 Introduction↩︎

In 2024, the mean global near-surface temperature exceeded the pre-industrial average by more than 1.5 K,[1] underscoring the increasingly tangible impacts of global warming on human health, behavior, and daily life.[2]–[4] A plausible behavioral response—particularly among vulnerable populations such as older adults and individuals with pre-existing health conditions—is to spend more time indoors to mitigate thermal discomfort during periods of extreme heat. However, empirical evidence remains limited regarding how climate change is reshaping patterns of time use.[5], [6]

Nevertheless, because humans spend most of their time indoors, the impacts of climate change are largely mediated through indoor environments.[7], [8] Accordingly, understanding how rising outdoor temperatures influence indoor air quality has become an increasingly important research priority.[9], [10] Models and simulations are expected to play a central role in evaluating the long-term effects of climate change on indoor environments.[11], [12] However, these modeling efforts rely on accurate physicochemical representations of the complex surface–air chemistry occurring indoors.[13]

In this study, we address this challenge by investigating the permeation of two representative volatile and semivolatile organic compounds (VOCs) through human skin. The compounds of interest—acetone and 6-MHO (Figure 1)—are formed during the ozonolysis of skin lipids, which act as the dominant ozone scavengers in densely occupied indoor environments and thus constitute an important source of secondary indoor pollutants.[14], [15] Kinetic models have been developed to describe the mass transport and chemical transformations of these oxidation products across the skin, surrounding air, and clothing.[16], [17] However, the reliability of such models depends critically on accurate physicochemical parameters describing the partitioning of semivolatile compounds—such as 6-MHO—between the human body and indoor air,[18] as well as on reliable descriptions of their diffusional transport through the skin barrier.

To address this challenge, we recently began investigating permeation energetics using atomistic molecular dynamics (MD) simulations of stratum corneum (SC) lipid membranes (Fig 1)—the skin’s primary barrier to chemical penetration.[19] We found that the convergence of free-energy (FE) profiles describing the translocation of permeants across these membranes was remarkably slow, requiring microseconds of simulation time even when state-of-the-art enhanced sampling techniques were employed. We traced the molecular origin of this convergence bottleneck to the highly ordered structure of SC membranes, which contrasts sharply with the fluid-like organization of conventional lipid bilayers such as those composed of 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC). In addition, we identified stochastic lipid flip–flop events between the two membrane leaflets that induce long-lived structural asymmetries.

In this work, we extend our previous investigations by providing estimates of diffusivity profiles along the membrane normal. Predicting these diffusivities remains a challenging task—particularly when permeation is slow and enhanced sampling is required to overcome large FE barriers. We employ two complementary approaches to compute position-dependent diffusivities from harmonically restrained MD simulations, based on either the velocity autocorrelation function (VACF)[20] or the position autocorrelation function (PACF)[21] of the coordinate describing membrane permeation. As shown by our results—and consistent with the findings of Rowley et al. for a dipalmitoylphosphatidylcholine (DPPC) membrane[22]—estimates obtained from the two approaches differ substantially. Through a rigorous propagator-based analysis, we demonstrate that the resulting diffusivity profiles provide lower and upper bounds to the true profile. We further propagate these results to obtain single-membrane permeabilities and we critically assess how methodological uncertainties affect the accuracy of these predictions. Taken together, our work advances the understanding of molecular permeation across skin lipid membranes and establishes a robust framework for improving predictive models of indoor air chemistry under changing environmental conditions.

2 Methods↩︎

The model systems and simulation protocol have been described in detail in our previous study.[19] Here, we briefly summarize the main aspects relevant to the present work and describe the computation of position-dependent diffusivities and the propagator analysis.

2.1 Model Systems↩︎

We studied three solutes—acetone, 6-MHO, and water—representing two skin-oil oxidation products and a reference compound. These solutes permeate through a model stratum corneum (SC) lipid bilayer, originally proposed by Klauda et al., and representative of the short periodicity phase (SPP)—the tightly packed lamellar arrangement of lipids that forms the primary diffusion barrier in the SC (Figure 1).[23], [24]

2.2 Umbrella Sampling↩︎

The translocation of the solutes across the membrane is described by a collective variable (CV), denoted as \(z\), which measures the center-of-mass distance between each solute and the membrane projected onto the membrane normal (Figure 1).

For each solute, we performed 91 independent simulations (windows), in which the solute was harmonically restrained at integer-Å positions along the \(z\) coordinate (umbrella sampling, US). A detailed analysis of the influence of the harmonic force constant was presented by Rowley et al.,[22] and our own related findings on this topic are provided in Section S1 of the Supporting Information (SI). For the results reported in the main text, a force constant of 30  mol Å⁻⁻² is employed. Each simulation is initiated from well-equilibrated configurations obtained in our previous study and run for 15 ns per window.

All simulations are conducted at the physiological skin temperature of \(\require{physics} T = \qty{305.15}{\kelvin}\), maintained by a stochastic velocity-rescaling thermostat[25] with a time constant of 1 ps, and performed using NAMD 3.0 (alpha13).[26], [27]

2.3 Diffusivity Profiles↩︎

Two approaches are used to compute the position-dependent diffusivities, \(D(z)\), along the membrane normal.

The first approach, proposed by Woolf and Roux, relies on the calculation of the velocity autocorrelation function (VACF):[20] \[C_v(t) = \left<\dot{z}(t)\dot{z}(0)\right>,\] followed by the evaluation of a Laplace-frequency–dependent diffusivity, \[D(s) = \frac{ -\hat{C}_v(s)\left<\delta z^2\right>\left<\dot{z}^2\right> }{ \hat{C}_v(s)\!\left[s\!\left<\delta z^2\right> + s^{-1}\!\left<\dot{z}^2\right>\right] - \left<\delta z^2\right>\left<\dot{z}^2\right> }, \label{eq:ds}\tag{1}\] where \(\hat{C}_v(s)\) denotes the Laplace transform of \(C_v(t)\). The local diffusion coefficient is then obtained in the limit \[D(z=\left<z\right>) = \lim_{s \to 0} D(s), \label{eq:ds95extrapol}\tag{2}\] which, however, is mathematically ill-defined and therefore requires an empirical extrapolation procedure.

Hummer later proposed an alternative approach that avoids Laplace-domain extrapolation by using the position autocorrelation function (PACF):[21] \[\begin{align} C_{\delta z}(t) &= \frac{\left<\delta z(t)\delta z(0)\right>}{\left<\delta z^2\right>}, \\ D(z=\left<z\right>) &= \frac{\left<\delta z^2\right>}{\int_0^\infty C_{\delta z}(t) \, dt}. \label{eq:dt} \end{align}\tag{3}\] In practice, however, one must either impose a cutoff time for the integral or assume a functional form to fit the correlation decay—both of which can introduce systematic bias.

The velocities required for computing the VACF are obtained from finite differences of the CV trajectory, which is saved at 2 fs intervals. Position-dependent diffusivities are determined using the ACFCalculator code developed by Rowley et al., which implements a heuristic procedure for the extrapolation in eq. 2 .[22]

We note that alternative strategies for estimating position-dependent diffusivities are available. In particular, Bayesian-inference methods were originally proposed by Hummer[28] and later extended to anisotropic media by Ghysels et al.[29]. They have also been coupled to adaptive biasing force calculations by Comer, Chipot, and co-workers.[30], [31] Complementary formulations have been developed by Rosta and collaborators.[32], [33] We focus here on the VACF/PACF approaches because they are widely used in membrane-permeation studies,[34]–[38] straightforward to implement on the same umbrella-sampling data, and allow direct comparison to prior results for our benchmark system (water).

2.4 Propagator Analysis↩︎

Diffusion in inhomogeneous media can be described by the Smoluchowski equation. In the one-dimensional case relevant to this study, the position- and time-dependent probability density \(p(z, t)\) evolves according to \[\partial_t p(z, t) = \partial_z \left\lbrace D(z) e^{-\beta \Delta F(z)} \, \partial_z \!\left[ e^{\beta \Delta F(z)} p(z, t) \right] \right\rbrace ,\] where \(\beta = 1/(k_\mathrm{B}T)\) is the inverse thermal energy, \(\Delta F(z)\) denotes the FE profile, and \(D(z)\) represents the position-dependent diffusivity.

For given \(\Delta F(z)\) and \(D(z)\) profiles and a delta-function initial condition, \(p(z, t_0) = \delta(z - z_0)\), this equation can be solved numerically using a standard forward-in-time, centered-in-space (FTCS) finite-difference scheme. The resulting solution, \(p(z, t \,|\, z_0, t_0)\), represents the propagator of the diffusion model—that is, the conditional probability of finding a particle initially located at \(z(t_0) = z_0\) at position \(z(t)\) at a later time \(t\). The same quantity can also be obtained directly from molecular dynamics (MD) simulations without assuming diffusive behavior, thereby providing a rigorous test of the accuracy of the proposed \(\Delta F(z)\) and \(D(z)\) profiles.

We compute MD-derived propagators for the critical membrane-core region as follows. A total of 1800 snapshots were selected in which the solute was located near the membrane center (\(|z| < \SI{0.05}{\angstrom}\)). The harmonic restraint was then removed, allowing the solute to freely explore the entire membrane. Unrestrained simulations were run for 15 ns, and propagators \(p(z, t \,|\, z_0 = 0, t_0 = 0)\) for selected lag times \(t\) were estimated from histograms along the \(z\) coordinate.

3 Results and Discussion↩︎

3.1 Diffusivity Profiles↩︎

Diffusivity profiles, \(D(z)\), of the solutes along the transmembrane axis were estimated from the same set of simulations using the PACF and VACF approaches described above. The resulting profiles are shown in Figure 2.

For all three systems, the PACF method yields constant diffusivities in the aqueous regions (\(|z| > \SI{40}{\angstrom}\)) surrounding the membrane, which agree well with the corresponding bulk-phase diffusivities obtained from separate simulations (Section S2). The VACF method produces comparable results in this region for acetone and 6-MHO but deviates by approximately 24 % for water.

In the headgroup region, both methods indicate a local diffusivity bottleneck for all solutes. The most pronounced discrepancy between the two approaches, however, appears at the membrane center, where the estimates differ both quantitatively and qualitatively: the VACF method predicts higher-than-bulk diffusivities, whereas the PACF method yields lower values. Similar observations were reported by Rowley et al.for a DPPC membrane.[22] In the following, we discuss the challenges associated with determining diffusivities using both approaches and the likely origins of these discrepancies. Our analysis focuses primarily on two representative regions of the membrane: the aqueous region, where the two methods show good agreement, and the membrane center, where the discrepancies are most pronounced.

3.2 Challenge 1: Extrapolation of the VACF Results↩︎

Obtaining position-dependent diffusivities from VACFs requires computing a Laplace-frequency–dependent diffusivity, \(D(s)\), for each US window (eq. 1 ), followed by extrapolation to zero frequency (eq. 2 ). However, this limit is ill-defined because of the \(s^{-1}\) dependence of \(D(s)\) and therefore necessitates an empirical extrapolation procedure. Our VACF-based diffusivity profiles were obtained using the automated linear extrapolation method of Rowley et al.[22] Although this approach is relatively robust in the aqueous bulk region, it may not be optimal at the membrane center (Figure 3). Because \(D(s)\) exhibits pronounced curvature in the membrane core, the linear extrapolation is highly sensitive to the placement of the tangent line, yielding diffusivity estimates between roughly 0.1 Å² ps⁻¹ and 0.75 Å² ps⁻¹.

3.3 Challenge 2: Determination of PACF Correlation Times↩︎

The challenges associated with extrapolating \(D(s)\) to zero frequency were recognized by Hummer, who proposed an alternative approach based on determining the PACF correlation time (eq. 3 ).[28] For purely diffusive dynamics, this method is formally equivalent to the VACF approach, but it is considerably simpler to apply in practice.

In Figure 4, we show selected diffusivities, \(D(t)\), as a function of the integration cutoff time used in eq. 3 . The PACF method performs well when the solute resides in the aqueous phase, where \(D(t)\) reaches a clear plateau between approximately 8 ps and 10 ps. Although noise leads to some fluctuations, it does not significantly alter the predicted diffusivity. In contrast, at the membrane center, the decay of the PACF is substantially slower, and no clear plateau in \(D(t)\) is observed even after 10 ps. These observations are consistent with those reported by Rowley et al.for a DPPC membrane.[22]

As a consequence, solute- and region-specific integration cutoffs would be ideal but require tedious manual adjustment. Instead, the diffusivity profiles in Figure 2 were obtained using practical compromise cutoffs: 8 ps for water and 10 ps for acetone and 6-MHO. Although these fixed cutoffs may slightly overestimate PACF-derived diffusivities in certain regions (e.g., the membrane center), our main conclusions are unaffected: PACF-derived profiles still represent lower bounds to the true \(D(z)\).

3.4 Propagator Analysis↩︎

Having obtained two qualitatively distinct diffusivity profiles from the same set of simulations, it is natural to ask which of the two, if either, is correct. To address this question, we computed the propagators of both diffusion models in the critical membrane-center region, \(p(z, t \,|\, z_0 = 0, t_0 = 0)\)—that is, the probability of finding a particle initially located at \(z_0 = 0\) at position \(z\) after a lag time \(t\). The corresponding FE profiles were taken from our previous study,[19] as shown in Figure 5. For the present analysis, we used the umbrella-sampling (US)–based free-energy profiles. A comparison between the model-derived propagators and those obtained directly from molecular dynamics simulations is presented in Figure 6 for two representative lag times.

At short lag times (\(t = \SI{2}{\pico\second}\)), the VACF-based model shows quantitative agreement with the simulation data over six orders of magnitude, whereas the PACF approach fails to do so. This finding is consistent with the results reported by Rowley et al.,[22] who performed a simplified analysis by approximating the propagator as a normal distribution—an assumption that holds only for lag times long enough to reproduce diffusive dynamics yet short enough that the effects of inhomogeneity remain negligible.

However, it is not obvious that the lag time of 2 ps chosen here (or the 5 ps used by Rowley et al.) is sufficiently long for the diffusive approximation to hold. Indeed, the failure of the PACF to decay within such short timescales may also be interpreted as a breakdown of this assumption. Subdiffusive behavior has, in fact, been reported for lag times up to at least 64 ps for methanol permeation across a POPC membrane.[31]

We therefore extended our propagator analysis to selected lag times of up to \(\SI{500}{\pico\second}\). For a range of lag times between 50 ps and 500 ps, which we consider to lie within the truly diffusive regime, the model propagators exhibit the behavior shown for the representative case of \(t = \SI{100}{\pico\second}\) in Figure 6. The complete set of results is provided in Section S3.

At these timescales, neither the VACF- nor the PACF-based models reproduce the simulation data quantitatively. However, the VACF results consistently overestimate the diffusive spread (i.e., the diffusivities are too large), whereas the PACF results underestimate it (i.e., the diffusivities are too small). We therefore interpret the two profiles as providing upper and lower bounds, respectively, to the true diffusivity profile, which must lie between them. In the following, these bounds are used to assess the sensitivity of derived quantities, specifically the transmembrane permeabilities.

3.5 Impact on Permeabilities↩︎

Within the inhomogeneous solubility–diffusion (ISD) model,[39], [40] the local resistance of a membrane to permeation, \(R(z)\), is expressed in terms of the diffusivity and FE profiles as \[R(z) = \frac{e^{\beta \Delta F(z)}}{D(z)} \, . \label{eq:resistance}\tag{4}\] The overall permeability of the membrane, \(P\), is then obtained by integrating the resistance over the membrane width, \(h\), according to \[\frac{1}{P} = \int_{-h/2}^{+h/2} R(z) \, dz . \label{eqn::perm}\tag{5}\]

Using the FE profiles for the translocation of all solutes across the membrane (Figure 5), we computed the membrane permeabilities according to eqs 4 and 5 . The resulting permeabilities are shown in Figure 7, and the corresponding numerical values are provided in the Supporting Information (Section S4). Because of the exponential dependence on the free energy in eq 4 , even small variations in \(\Delta F(z)\) are amplified, leading to pronounced differences in the resulting permeabilities, which follow the order 6-MHO \(>\) acetone \(\gg\) water.

We have also shown that the FE profiles obtained from US and WTM differ slightly (Figure 5), even after extensive sampling, although these differences remain within the limits of statistical uncertainty.[19] With the diffusivity profiles now available, we can propagate the different FE profiles to evaluate their impact on the resulting transport coefficients. We find good agreement between the results obtained from both methods for water and 6-MHO; for acetone, however, noticeable deviations are observed. Nevertheless, the permeability coefficients for all solutes differ by less than an order of magnitude, indicating that statistical variations between FE methods likely have little to no impact on kinetic model predictions.

More striking are the differences in permeabilities when comparing VACF- and PACF-derived diffusivities. Because the VACF-derived diffusivities are consistently higher than those obtained from the PACF method, the same trend is reflected in the corresponding permeabilities. In absolute terms, deviations reach approximately one order of magnitude. Whether such differences are significant for kinetic models is highly problem specific. In many cases, these discrepancies are likely acceptable, although the final assessment should be guided by a dedicated sensitivity analysis.[17]

Finally, we examined the influence of the integration limits used in eq. 5 . Typically, \(h\) is defined as the full thickness of the membrane. However, for heterogeneous, multicomponent lipid membranes such as the SC system investigated here, determining precise interfacial boundaries (\(\pm h/2\)) is nontrivial. We therefore considered three plausible definitions for \(h/2\), based on the membrane density profile (Section S4): (i) the position where the water density begins to decrease from its bulk value (\(\require{physics} h_1/2 = \qty{32.0}{\angstrom}\)); (ii) the position where the water density is reduced to half of its bulk value (\(\require{physics} h_2/2 = \qty{25.5}{\angstrom}\)); and (iii) the position where the water density approaches zero (\(\require{physics} h_3/2 = \qty{17.0}{\angstrom}\)).

To investigate the impact of these definitions, we show the local resistance \(R(z)\) in Figure 8. The area under each curve corresponds to the inverse permeability, \(P^{-1}\). Within the membrane interior, this area is several orders of magnitude larger than the small additional contributions from the interfacial regions. We therefore conclude that different reasonable integration limits have little to no effect on the computed permeabilities.

4 Conclusions↩︎

We used MD simulations to quantify position-dependent diffusivities, \(D(z)\), for water, acetone, and 6-MHO across a representative SC lipid membrane and combined these results with previously reported FE profiles to compute transmembrane permeabilities, \(P\). By applying both PACF- and VACF-based analyses to the same set of US trajectories and validating them through a propagator comparison, we demonstrated that the two \(D(z)\) estimates systematically bracket the true diffusivity profile in the diffusive regime, with PACF providing a lower bound and VACF an upper bound. Permeabilities derived from these bounds consequently span approximately one order of magnitude, whereas differences arising from the choice of FE method (US vs.WTM) become negligible once propagated to transport parameters. Consistent with the exponential dependence of the ISD framework on free energy, the ordering of permeabilities follows the underlying FE landscape rather than the local diffusivities—6-MHO \(>\) acetone \(\gg\) water—highlighting that permeation across the membrane is governed primarily by energetic barriers rather than by differences in molecular mobility.

Direct comparison of permeabilities computed from our MD simulations to experimental values (in vitro[41], [42] or in vivo[43]) is not straightforward because our model represents a simplified SC membrane composed of an equimolar mixture of ceramides, free fatty acids, and cholesterol. In contrast, native human SC exhibits substantial lipid heterogeneity, encompassing a broad spectrum of ceramide subclasses and chain-length distributions.[44] Nevertheless, the water permeability obtained using PACF-derived diffusivities is within the same order of magnitude as the experimental estimate of Gupta et al.(\(7.08 \pm 1.18 \times 10^{-5}\) cm s⁻¹ ),[35] with residual differences most likely reflecting variations between the GROMOS and CHARMM force fields.

Permeability coefficients are often estimated empirically using quantitative structure–activity relationship (QSAR) and continuum diffusion models,[45]–[49] which correlate measured permeabilities with molecular descriptors such as molar mass, partition coefficients, and bulk diffusion constants. Although such models are powerful for interpolation within their training domains, they do not explicitly capture the molecular-scale interactions and dynamic heterogeneity that govern transport across the SC lipid matrix—especially when applied to molecules that fall outside their original parameterization range. In contrast, MD simulations capture these effects directly by resolving how solutes interact with and move through the heterogeneous lipid environment. For example, despite its larger molar mass relative to water and acetone, 6-MHO exhibits the highest permeability because of its comparatively low free-energy barriers to permeation (Figure 5), underscoring the dominant role of thermodynamic driving forces—rather than molecular size—in determining SC transport.

With respect to indoor air chemistry, several kinetic models—such as KM-BL, KM-FILM, and KM-SUB-Skin-Clothing[17]—have been developed to describe the dynamics of skin-oil oxidation products across different interfaces. However, these models critically depend on physicochemical parameters that are difficult or impossible to measure experimentally, underscoring the value of molecular-scale simulations like those presented here. For example, Lakey et al.demonstrated that the performance of their KM-FILM model is highly sensitive to the assumed diffusion constants of surface films,[50] illustrating how improved microscopic input can directly enhance macroscopic predictive accuracy.

Taken together, our results provide mechanistic and quantitative constraints—via upper and lower bounds on \(D(z)\) and permeabilities—that can be directly used to refine existing indoor-air kinetic models. Beyond these applications, the present framework establishes a foundation for systematically linking atomistic transport processes to macroscopic exposure and air-quality models. Future directions include (i) extending the approach to multilamellar SC architectures and anisotropic membrane transport; (ii) assessing the influence of force-field and water-model choice; (iii) developing more robust and transferable diffusivity estimators; and (iv) performing targeted experimental validation for VOCs of indoor relevance. Collectively, these efforts will help bridge the scale gap between molecular transport and human exposure in complex indoor environments.

References↩︎