Statistical Mechanical Theories of Membrane Permeability

Abstract

A popular theoretical framework to compute the permeability coefficient of a molecule is provided by the classic Smoluchowski-Kramers treatment of the steady-state diffusive flux across a free-energy barrier. Within this framework, commonly termed "inhomogeneous solubility-diffusion" (ISD), the permeability, P, is expressed in closed form in terms of the potential of mean force and position-dependent diffusivity of the molecule of interest along the membrane normal. In principle, both quantities can be calculated from all-atom MD simulations. Although several methods exist for calculating the position-dependent diffusivity, each of these is at best an estimate. In addition, the ISD model does not account for memory effects along the chosen reaction coordinate. For these reasons, it is important to seek alternative theoretical formulations to determine the permeability coefficient that are able to account for the factors ignored by the ISD approximation. Using Green-Kubo linear response theory, we establish the familiar constitutive relation between the flux density across the membrane and the difference in the concentration of a permeant molecule, j = PΔC. On this basis, we derive a time-correlation function expression for the nonequilibrium flux across a membrane that is reminiscent of the transmission coefficient in the reactive flux formalism treatment of transition rates. An analysis based on the transition path theory framework is exploited to derive alternative expressions for the permeability coefficient. The different strategies are illustrated with stochastic simulations based on the generalized Langevin equation in addition to unbiased molecular dynamics simulations of water permeation of a lipid bilayer.

Figure 1:

Physical situation of a membrane separating a bulk region $a$ (blue) and a bulk region $b$ (yellow) containing a permeant molecule. The indicator function $H_a\left(H_b\right)$ is equal to 1 if the permeant molecule $i$ is on side $a\left(b\right)$, and equal to zero otherwise. By construction, all space is covered by the two indicator functions, with $H_a+H_b=1$. The concentration of permeant molecule at the positions $z_a$ and $z_b$ in the bulk region away from the membrane serves as a reference to establish the constitutive relation $j=P\Delta C$ of membrane permeation.

Figure 2:

Part of a ${\text{BD}}_e$ (Brownian dynamics with exponential memory) trajectory of a narrow barrier, showing color-coded crossings between $a$ (red) and $b$ (blue). The two regions are defined by the edge of the boundary regions $z<z_a$ for the side $a$ and $z>z_b$ for side $b$, with $z_a=-35\AA$, and and $z_b=+35\AA$.

Figure 3:

MSD diffusivity computed by ${\text{BD}}_e$ simulations (solid green line) on a flat free energy landscape $\left(W\right(z)=0)$ with various values of $K_e$ and $D_e$. In each case, $D_z$ is $0.1{\AA }^2/\text{ps}$. The corresponding analytical diffusivity computed with eq 26 is shown for comparison (dotted black line). The limits $D_z$ (dashed red line) and $D_{\text{eff}}$ (dashed blue line) are shown. A. $D_e=0.1{\AA }^2/\text{ps},K_e=0.05\text{kcal}/\text{mol}{\AA }^2$, B. $D_e=0.2{\AA }^2/\text{ps},K_e=0.05\text{kcal}/\text{mol}{\AA }^2$, C. $D_e=0.1{\AA }^2/\text{ps},K_e=0.10\text{kcal}/\text{mol}{\AA }^2$.

Figure 4:

Effective MSD diffusivity as a function of the lag time computed analytically from eq 26 with an exponential memory for different barriers and, therefore, dynamical regimes. In each case, the instantaneous diffusivity $\left(D_z\right)$ is set to $0.1{\AA }^2/\text{ps}$ (red dashed line), the minimum effective diffusivity is $0.05{\AA }^2/\text{ps}$ (blue dashed line), and the memory constant $\left(K_e\right)$ is 0.05 kcal/mol ${\AA }^2$. A barrier with a negative value of $W^{″}$ corresponds to an inverted parabola (black line) and a positive value corresponds to a harmonic well (purple line), while zero is a flat landscape (green line).

Figure 5:

PMF, $W\left(z\right)$, (left) and analytical MSD diffusivity (right) for demonstrating the Grote-Hynes correction to the ISD permeability of an inverted harmonic barrier with exponential memory. The PMF has a curvature $\left|W^{″}\right|=0.2\text{kcal}/\text{mol}{\AA }^2$, and it reaches 0 where $z=\pm 6$. The MSD corresponds to parameters: $D_z=0.1,D_e=0.01{\AA }^2/\text{ps},K_e=0.4\text{kcal}/\text{mol}{\AA }^2$.

Figure 6:

PMF, $W\left(z\right)$, and approximate committor, $q\left(z\right)$, for the narrow (left side) and flat-top (right side) barriers used in the ${\text{BD}}_e$ and BD simulations.

Figure 7:

Total number of crossing events from ${\text{BD}}_e$ (solid lines) and BD (dashed lines) simulations of each barrier (narrow = blue, flat-top = green). The mean and standard deviation of total crossings per replicate from ensembles of 50 replicates of $50\mu \text{s}$ each are 42.12 ± 6.4 for the BD narrow, 27.20 ± 4.9 for the ${\text{BD}}_e$ narrow, 43.94 ± 6.5 for the BD flat-top, and 24.48 ± 5.6 for the ${\text{BD}}_e$ flat-top.

Figure 8:

Top row: Heaviside function ($C_{ab}\left(t\right)$, black) and committor time-correlation function ($C_{qq}\left(t\right)$, blue) for BD (dotted) and ${\text{BD}}_e$ (solid) simulations of the narrow (left) and flat-top (right) barriers. Bottom row: Permeability over a range of lag times for the ${\text{BD}}_e$ simulations only. Permeability corresponding to the ISD model (green) computed with eq 1, “Heaviside” (black) with eq 9, “committor” (blue) with eq 19, “short traj” (dashed blue) with eq 36, and “counting” (red) with eq 13. The counting permeability does not depend on the lag time. The ISD permeability is the same for each barrier by construction. The permeability determined from the Heaviside time-correlation function reaches approximately the same plateau after 2000 – 3000 ps (shown in section S3 of the SI).

Figure 9:

A. PMF of water across a DLPC bilayer, computed from the marginal probability distribution from unbiased MD with the Drude2023 force field. B. Approximate committor of water, which is “unwrapped”, or extended, to cover several periodic images of the membrane system. The committor is calculated from the PMF assuming that the diffusivity is independent of $z$. C. Committor $\left(C_{qq}\left(t\right)\right)$ and Heaviside $\left(C_{ab}\left(t\right)\right)$ time-correlation functions, averaged from each unwrapped water trajectory, mapped onto the unwrapped committor and Heaviside function, respectively. D. Permeability of DLPC to water based on the ISD model (green) with eq 1, “Heaviside” (black) with eq 9, “committor” (blue) with eq 19, and “counting” (red) with eq 13. The ISD permeability depends on the position-dependent diffusivity $D\left(z\right)$, which was determined with biased US simulations. The lag time dependence of the ISD is related to the upper limit of the integral of the position autocorrelation function (see section S4 of the SI). The dashed green line corresponds to an ISD estimate using the $D\left(z\right)$ profile of Figure S4. The counting permeability is based on a mean and standard deviation of 52.0 ± 8.0 crossings per replicate, where three replicates were simulated for 400 ns each.

Figure 10:

Successful water crossings beginning at the time the molecule left region $a$ (red) or $b$ (blue) and ending at the time at which it first touches the opposite region. The trajectories are colored based on the initial region. The dashed black line is the center of each image of the membrane (see Figure 9B.). The mean reactive transit time is interpreted as the mean duration of the crossings, which is 2.93 ns.