Sampling efficiency of the counting

method for permeability calculations

estimated with the inhomogeneous
solubility–diffusion model
Cite as: J. Chem. Phys. 154, 054106 (2021); https://doi.org/10.1063/5.0033476
Submitted: 16 October 2020 . Accepted: 08 January 2021 . Published Online: 02 February 2021

Samaneh Davoudi, and An Ghysels

J. Chem. Phys. 154, 054106 (2021); https://doi.org/10.1063/5.0033476 154, 054106

© 2021 Author(s).
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

Sampling efficiency of the counting method

for permeability calculations estimated
with the inhomogeneous solubility–diffusion
model
Cite as: J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476
Submitted: 16 October 2020 • Accepted: 8 January 2021 •
Published Online: 2 February 2021

Samaneh Davoudi and An Ghyselsa)

AFFILIATIONS
IBiTech–BioMMeda Group, Ghent University, Corneel Heymanslaan 10, Block B-entrance 36, 9000 Gent, Belgium

a)
Author to whom correspondence should be addressed: an.ghysels@ugent.be

ABSTRACT
Permeability is a key property in various fields such as membrane technology for chemical separation and transport of substances through
cell membranes. At the molecular scale, the counting method uses the number of membrane crossings in a conventional unbiased molecular
dynamics simulation to predict the permeability. This contribution investigates under which conditions the counting method has insufficient
statistics. An equation is derived for a compartmental model based on the inhomogeneous solubility–diffusion (Smoluchowski) model, giving
insight into how the flux correlates with the solubility of permeants. This equation shows that a membrane crossing is a rare event not only
when the membrane forms a large free energy barrier but also when the membrane forms a deep free energy well that traps permeants. Such
a permeant trap has a high permeability; yet, the counting method suffers from poor statistics. To illustrate this, coarse-grained MD was run
for 16 systems of dipalmitoylphosphatidylcholine bilayer membranes with different permeant types. The composition rule for permeability
is shown to also hold for fluxes, and it is highlighted that the considered thickness of the membrane causes uncertainty in the permeability
calculation of highly permeable membranes. In conclusion, a high permeability in itself is not an effective indicator of the sampling efficiency
of the counting method, and caution should be taken for permeants whose solubility varies greatly over the simulation box. A practical
consequence relevant in, e.g., drug design is that a drug with high membrane permeability might get trapped by membranes thus reducing its
efficacy.
Published under license by AIP Publishing. https://doi.org/10.1063/5.0033476., s

I. INTRODUCTION On the other hand, computational simulations have been

proven to give more insight into permeation at the atomic scale.15–17
Permeability is a key property to describe the kinetics of mass Among simulation methods, molecular dynamics (MD) can lay
transport of a substance through a membrane. It is highly relevant in open more details on the kinetics of the process, the property that
various research and application fields, such as membrane technol- experiments often cannot provide directly, under the condition that
ogy for chemical separation1–4 and passive permeation of substances a proper force field for the molecular interactions is available.18 Dif-
through the cell membrane.5,6 In biology, permeation through the ferent simulation methods have been applied to calculate the perme-
cell membrane is a vital process that is relevant for drug delivery5,6 ability, many of which rely on enhanced sampling methods18–20 and
and the regulation of cell death and proliferation.7 A large number the inhomogeneous solubility–diffusion (ISD) model based on the
of experimental studies have been carried out to provide information Smoluchowski equation.21–28 In the ISD model, the permeability P
on this process during the past decade.8–14 However, these methods is expressed in terms of local diffusivity and free energy gradients.21
are either expensive or not informative enough, as they are usually A limitation of this approach, however, is that it is generally not
unable to provide data on a molecular level. clear a priori to which degree and in which situation the diffusive

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-1

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

hypothesis is valid for the system under study. Moreover, access- water layers will be treated as an extra parameter in the discus-
ing the permeability in this model can be challenging regarding the sion. Moreover, this parameter will allow us to discuss the overall
simulation setup and Bayesian analysis (BA).29 sampling efficiency. When the limit of the water layer thickness
Permeability can also be determined by counting crossings is taken to cover the whole simulation box, the associated num-
through the membrane.17,30 In the case of rare events, there might ber of crossings represents crossings through the whole simulation
be an insufficient number of crossings, and the statistics can be poor. box.
To set the mind, imagine aiming for a 10% error on the number of Section II will use the ISD model to derive a theoretical expres-
crossings. When crossings are occurring independently according to sion of how the counting method is affected by the solubility of
a waiting time distribution of a Poisson process, at least 100 events permeants in the membrane. Section III discusses the dependen-
would be needed to achieve the 10% standard error.25 An example cies and limiting cases of the derived equation. It is shown under
of using the counting method is oxygen and water diffusion through which conditions membrane crossings can be considered as a rare
phospholipid membranes.24,25,31–33 event. A composition rule for the flux is derived as well. Section IV
In this paper, the number of crossings expected in a simula- presents the number of crossings of 16 different permeant types
tion will be investigated as a function of the simulation properties. through the dipalmitoylphosphatidylcholine (DPPC) bilayer mem-
Several properties are dictated by the simulation setup, such as the brane, illustrating the derived equation. Section V summarizes our
number of permeants and box dimensions. A straightforward way findings.
to improve the sampling of the unit cell is to increase the system
size. When the cross section doubles and the total number of per-
meants doubles, the concentration in the solvent phase will remain II. THEORY
unaltered, while the number of crossings would double. There are
The Smoluchowski equation describes the transport of perme-
unfortunately limits to this approach because the computational cost
ants through an inhomogeneous medium on a free energy profile
increases faster than linearly with the box size, and in comparison, it
F(z) with a position-dependent diffusivity D(z), where z is the coor-
would be more efficient to double the simulation time T sim . Another
dinate along the membrane normal.21 Assume F ref is the free energy
approach could be to insert more permeants without altering the
at the reference location, usually taken to be in the water phase
simulation box size, meaning that all concentrations—both solvent
on both sides of the membrane. The permeability P of a layer of
and membrane—increase. One should be wary of the clustering of
thickness h follows from solving the Smoluchowski equation under
the permeants, however, and the permeability might become con-
steady-state conditions,21
centration dependent. For instance, oxygen showed no clustering for
oxygen concentrations about ten-fold higher than the oxygen solu-
bility in water,24 while ethanol permeation through phosphatidyl- 1 eβF(z)
= e−βFref ∫ dz, (1)
choline (POPC) membranes differs substantially between low and P h D(z)
high concentrations.29,34
Other properties are inherent to the chosen membrane, solvent,
where the integration runs over a certain thickness h, and
and permeants. This paper considers the solubility of permeants
β = 1/(kB T) is the inverse temperature with T being the temper-
(or, equivalently, the partitioning coefficient K) and the membrane
ature. In a compartmental model, as illustrated earlier for oxygen
thickness hm . It is well known that the membrane forms a high free
transport,24 F and D are piecewise constant functions, and the
energy barrier when a permeant is nearly insoluble in the membrane,
integral becomes a sum over the several compartments,
and the associated low permeability prevents good sampling of the
unit cell. It will now be shown that a high permeability, for perme-
1 1 hi
ants that are highly soluble in the membrane, is another cause for =∑ =∑ , (2)
poor sampling. In contrast to low permeability, a high permeabil- P i P i i K i Di
ity membrane is often regarded as a membrane where the flux is
not rate-limiting, but we will show that high permeability cannot where each layer i is characterized by its thickness hi , diffusivity Di ,
guarantee good sampling of the simulation box either. and free energy F i . The factor K i = exp(−β(F i − F ref )) is the partition-
The membrane thickness follows from the self-aggregation of ing constant in layer i compared to the water phase, and it is equal
phospholipid molecules in the bilayer. It is less discussed in the per- to the concentration ratio between layer i and the water phase, i.e.,
meability literature than the partitioning coefficient. Still, using the K i = ci /cref . The composition rule in Eq. (2) shows that the per-
simplified Overton’s model P = KD/hm ,35 where D is the diffusion meation through a series of layers can be seen as a series of local
constant, it is clearly of equal importance as the partitioning coef- resistances 1/Pi , if the same reference F ref is used for each layer.
ficient. The membrane thickness is however ill-defined since there Let us start with a one-compartmental model for a uniform
is no unique answer as to where an atom ends and another atom membrane of thickness hm , for which the integral in Eq. (1) sim-
begins, and moreover, membranes exhibit thermal undulations with plifies to
lipid protrusions. In experimental work, a roughly estimated mem-
brane thickness is used to estimate the permeability.10,36 In simula- KD
P= . (3)
tion work, it is often derived in relation to the plane of the phosphate hm
groups.29 Hence, another aim of this paper is to showcase the sensi-
tivity of the permeability to the considered thickness. Therefore, in The partitioning coefficient K relates to the free energy difference
addition to the membrane thickness, the thickness of two additional ΔF between the membrane and the water phase, where ΔF is also

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-2

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

referred to as the free energy of transfer to bring a permeant from to the previous equation, P ≈ KD/hm , as if the water layer is unim-
the water phase into the membrane. portant. The permeability is then indeed dominated by the height
This model is extended to the three-compartmental model of the barrier and is a linear function of K. This is expected from
depicted in Fig. 1, where a water layer of thickness d is added on the integral form in Eq. (1), which is dominated by high F values. In
both sides around the membrane. For simplicity, the diffusivity is the limit of large K, at least compared to hm /(2d), the membrane has
assumed to be constant over the whole simulation box. This model a high concentration of permeants, and P becomes approximately
will allow for a more complete discussion of the behavior of the a constant, P ≈ D/(2d), meaning that the permeability of the mem-
counting method. The permeability of this model is obtained by brane is not the limiting factor but rather the resistance coming from
putting the previous permeation resistance hm /(KD) in series with the water layer.
the permeation resistance of the water phase 2d/D, where the par- Given this compartmental model, we will now turn toward the
titioning coefficient of the water phase is simply 1. This gives the practical situation of a simulation setup. In a simulation, P can be
permeability of the membrane plus water layers equal to computed using the counting method, by determining the cross-
ing rate J and the reference concentration cref in the water phase.
D The number ncross of complete transitions through a membrane with
P= . (4) cross section area σ in either direction is counted during a long equi-
hm
2d + librium MD simulation of length T sim . Next, the (bidirectional) flux
K J per unit of time and unit of area in either direction is derived,

Filling in d = 0 in Eq. (4) results back in Eq. (3). ncross

Two regimes exist for P, depending on the partitioning coeffi- J= = ∣J← ∣ + ∣J→ ∣, (5)
Tsim σ
cient K of the membrane. In the limit of low K, the concentration in
the membrane is much lower than in the water phase, and P reduces
where J ← and J → are the flux through the membrane in the negative
or positive direction, respectively. The well-known formula17,30 for
P from the counting method is

∣J← ∣ + ∣J→ ∣ J
P= = . (6)
2cref 2cref

The counting method for permeability relies on the number of cross-

ings. The aim is to discuss the number of crossings as a function of
the key properties of the membrane and the simulation setup: (1) the
free energy of transfer ΔF, (2) considered layer thickness h = hm + 2d,
(3) the total number of permeants N t that are present in the simu-
lation box, and (4) the dimensions of the simulation box, i.e., the
height L and cross area σ. For this reason, ncross will be written as a
function of these parameters in order to observe their role. Using the
first equality in Eq. (5), it is clear that a discussion of J contains sim-
ilar information as a discussion of ncross , apart from the role of the
simulation time T sim and cross area σ, and therefore, the discussion
of ncross will be largely based on the discussion of J.
From Eq. (6), it immediately follows

J = 2cref P, (7)

where P is given by Eq. (4). The task at hand is thus to express the
reference concentration located in the water phase, as a function of
the key parameters K and N t . The water concentration cref = cw is
the average number of permeants ⟨N w ⟩ in the water phase divided
by the volume of the water phase, which is σ(L − hm ) according
to Fig. 1. Similarly, the membrane concentration cm is the average
number of permeants ⟨N m ⟩ in the membrane divided by the mem-
FIG. 1. Simulation box with length L and cross area σ containing the membrane brane volume σhm . Note that, while the total number of permeants
(red) and water (blue). The compartmental model considers the membrane (thick- N t remains constant during the simulation, the number of particles
ness hm ) and a water layer on both sides (thickness d each). Diffusivity is assumed in the water and membrane, N w and N m , respectively, will fluctuate
to be constant, and the free energy profile is piecewise constant.
over time. Their averages are given as

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-3

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

⟨Nw ⟩ = cw σ(L − hm ), (8) III. DISCUSSION

The flux in the compartmental model [Eq. (11)] is first dis-
cussed for the two limiting cases where the layer does not comprise
⟨Nm ⟩ = cm σhm = Kcw σhm . (9)
any water (d = 0, so h = hm ) and where all water is included in the
layer (h = L), and next, it is illustrated for the intermediate cases (hm
Using the equality N t = N w + N m = ⟨N w ⟩ + ⟨N m ⟩, the reference < h < L).
concentration becomes For d = 0, only the membrane of thickness hm is considered,
and the flux simplifies. For large K (free energy well) and small K
Nt (free energy barrier), the dependence on K is constant or linear in K,
cref = cw = , (10) respectively,
σ(hm K + L − hm )

and we find J as a function of the key parameters ⎧ Nt D

⎪
⎪ , K≫1
⎪ σh2m
⎪
J(hm ) = ⎨ Nt D (16)
⎪
⎪ K, K ≪ 1.
2Nt D 1 1 ⎪
⎪
J= . (11) ⎩ σh(L − hm )
σ hm K + L − hm hm
2d +
K
When the membrane is highly soluble for the permeants (ΔF > 0),
the number of crossings is independent of K and thus ΔF. However,
This is the equation that will be used for the discussion. The num- when the permeants’ solubility in the membrane is lower than in the
ber of crossings follows immediately by multiplying J with T sim σ water phase (ΔF < 0), the number of crossings will drop exponen-
[Eq. (5)]. tially with ΔF. A very high free energy of transfer will give a very few
In a final step, we leave the case of the compartmental model membrane transitions, with poor statistics for the counting method
and return to the case of general F(z) and D(z) profiles that can as a consequence.
take any shape. It is still assumed that the transport kinetics can For 2d = L − hm , the whole water layer in the simulation box is
be described by the Smoluchowski equation, i.e., the kinetics are included in the counting of membrane transitions. The considered
Markovian. In the general case, the reference concentration is given thickness is h = L, so J(L) represents the flux of complete crossings
by of the simulation box,

e−βFref Nt
cref = . (12) Nt D 1 1
σ ∫ e−βF(z) dz J(L) = . (17)
σ(L − hm )2 hm 1 hm
L 1+ 1+ K
L − hm K L − hm
P is given by Eq. (1). Using the ⟨. . .⟩ notation for spatial averages in
the z-direction, over h or L, these are rewritten as Here, J(L) can be interpreted as a measure of the sampling efficiency
of permeants over the whole z-direction of the simulation box. It
eβFref is interesting to see that this equation is symmetric for the replace-
P= (13) ment of K by 1/K. Replacing ΔF by −ΔF does not alter the sampling
eβF(z)
h⟨ ⟩ efficiency J(L). In other words, J(L) is an even function of ΔF by
D(z) h the appearance of K in this expression. This is an important realiza-
tion: the flux in a simulation can be low because either a high barrier
and makes it difficult for the permeants to pass through the membrane
or a deep well makes it difficult to leave the membrane once they
e−βFref Nt are trapped. In summary, the larger the difference between the free
cref = . (14) energy minima and maxima across the membrane normal, given by
σL⟨e−βF(z) ⟩L
|ΔF|, the lower is the sampling efficiency because J drops with both
The flux in Eq. (7) now becomes K and 1/K. A high variety in solubility hence leads to lower sampling
efficiency.
We now approximate J(L) to allow comparison with the pre-
2Nt vious approximations of J(hm ) in Eq. (16). For large K (free energy
J= βF(z)
. (15)
e well) and small K (free energy barrier), J(L) is inversely proportional
hLσ⟨ ⟩ ⟨e−βF(z) ⟩
D(z) h L to K or linear in K, respectively,

The above equation is the generalization of the compartmental ⎧

⎪ Nt D 1
⎪
⎪
⎪ , K≫1
model in Eq. (11) to arbitrary free energy and diffusivity profiles. ⎪ σh(L − hm ) K
J(L) = ⎨ Nt D (18)
It gives the flux J as a function of the key parameters in the general ⎪
⎪
⎪ K, K ≪ 1.
⎪
case. ⎩ σh(L − hm )
⎪

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-4

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

When K is small, J(hm ) ≈ J(L). This can be understood because the According to Eq. (7), the permeability and the reference con-
exact location of the boundaries, be it at h = hm or h = L, is relatively centration are the main contributing factors, which are plotted in
unimportant when the barrier is high. When K is large, however, Figs. 2(b) and 2(c), for several values of d. The semi-log plots clearly
we find that J(h) ≫ J(L). When the membrane is highly soluble and show the constant (zero slope), linear (positive slope), and inversely
only few particles are in the water phase, the exact boundary loca- linear (negative slope) dependence on K in the limiting cases K ≫ 0
tion does matter. Omitting the water phase or including a small or and K ≪ 0. The number of crossings is mainly the product of these
larger part has a large effect on the number of crossings. The reason two graphs, and the semi-log plot in Fig. 2(a) is indeed the sum of
is that the permeability of the water phase is in this case the limiting the curves except for a vertical shift. The widespread curves for small
resistance for the overall permeability of the considered layer, hence K (i.e., ΔF ≪ 0) in Figs. 2(a) and 2(b) are a signature of the sensitivity
the sensitivity to the included water layer thickness. to the chosen layer thickness when the membrane is a permeant trap.
For d between 0 and (L − hm )/2, the number of crossings lies Let us focus on the number of complete crossings of the simulation
between J(hm ) and J(L). Figure 2 visualizes the trends for the toy sys- box (black line), which refers to the sampling efficiency. The vertical
tem (Fig. 1) with numerical values that are typical for MD settings, dashed line in Fig. 2(a) indicates 100 crossings corresponding with
for instance for a bilayer consisting of 72 phospholipid molecules.24 a standard error of about 10%. This corresponds to a free energy
The simulation box had sizes L = 70 Å and σ = (50 Å)2 , the mem- difference of about |ΔF| = 5kB T, indicated with vertical lines in the
brane thickness was set to hm = 55 Å, the box contained N t = 10 other plots. This means that the sampling of the complete simulation
permeating particles with diffusivity D = 0.5 Å2 /ps, and the total box becomes poor when |ΔF| is larger than 5kB T and one should
simulation time was set to T sim = 1 μs. be wary of statistics. When ΔF ≫ 0, the statistics are poor because
the membrane forms a high barrier making crossings a rare event.
However, when ΔF ≪ 0, the poor statistics are caused by the low
permeant concentration cw in the water phase [Fig. 2(c)]. When all
permeants are trapped in the membrane, meanwhile not contribut-
ing to the crossings, there is barely any permeant left in the water
phase as a candidate to perform a crossing.
Figure 3 plots the number of crossings as a function of the layer
thickness for several values of K. The curves for low K are almost
unaffected by considering some water layer. In contrast, the curves
for high K indeed decrease much faster than those for low K as those
are more sensitive to the layer thickness. For complete crossings at h
= L, the curve for a given K coincides with the corresponding curve
for 1/K (indicated with black crosses in the graph), illustrating the
symmetry between ΔF and −ΔF.
For all curves in Fig. 3, ncross decreases with increasing thick-
ness h. This is a consequence of the composition rule in Eq. (2) for
permeability. A higher thickness gives a lower permeability accord-
ing to Eq. (2). Using J = 2cref P, a similar composition rule holds for
the flux,

FIG. 2. Number of crossings ncross (a), permeability P (b), and reference concentra-
tion cw (c) in the model system of Fig. 1 as a function of the free energy difference
ΔF = −k B T ln K, where ΔF > 0 is a barrier and ΔF < 0 is a well. Various water lay-
ers d are considered. Blue dashed line: water slab not taken into account. Black
line crosses: complete simulation box considered; this is an indicator of overall FIG. 3. Number of crossings ncross as a function of the water layer d in the model
sampling. The horizontal dashed line in (a) is at 100 crossings, and the vertical system of Fig. 1, for various free energy barriers (K < 0) and free energy wells
dashed lines in (b) and (c) indicate the corresponding free energy difference. (K > 1). Black crosses: complete simulation box considered.

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-5

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

1 1 created manually. The MARTINI force field was used in all systems,
=∑ , (19)
J i J i mapping atoms of each phospholipid molecule and water molecule
into 12 beads and 1 bead, respectively.40
where J i is the bidirectional flux through layer i [see Eq. (5)]. An Water was modeled by applying the SPC/E model.41 The parti-
equivalent composition rule holds for ncross = T sim σ J. Hence, J cle mesh Ewald approach was employed to calculate the Coulombic
and ncross also decrease when the considered thickness is increased. interactions. Both the Coulombic and the van der Waals interac-
For a symmetric bilayer membrane, it immediately follows from tions were truncated at 1.1 nm, with the potentials shifted to zero
the composition rule that the entrance flux through one leaflet is at the cutoff using the potential modifiers. The neighbor list length
identical to the escape flux through that leaflet, J entr = J esc , and the is 1.4 nm updated using the Verlet neighbor search algorithm. With
complete crossing flux through the entire membrane is half of that, a time step of 30 fs, the leapfrog integrator was used to integrate
J = J entr /2 = J esc /2. Likewise, P = Phalf /2, where Phalf is the perme- the equations of motion. A temperature of 320 K was considered
ability of one membrane leaflet. It is important to note that these for the simulated systems by the coupling of velocity rescale ther-
composition rules are based on the compartmentalized inhomoge- mostats, with coupling constant set to 1.0 ps. At every time step,
neous solubility–diffusion model and are thus only true under the the center of mass motion of the system was removed and peri-
assumption of perfectly diffusive behavior. Non-Markovian kinet- odic boundary conditions were applied in all directions (x, y, z). A
ics might distort these composition rules,37 e.g., when memory Parrinello–Rahman barostat was employed to subject the box vec-
effects are at play for ethanol crossings through a phospholipid tors to semi-isotropic pressure, using a reference pressure of 1 bar,
membrane.29 a coupling parameter of 12 ps, and an isothermal compressibility of
Finally, the symmetry in ΔF is shown in the general case of arbi- 3 × 10−4 bar−1 . After adding the permeants, an energy minimiza-
trary profiles in Eq. (15). When the diffusion profile is taken to be tion was performed. Next, two equilibrations were performed: NVT
constant, D(z) = D, and considering the full simulation box, h = L, simulation of 30 ns and NPT simulation of 30 ns. The production
the flux becomes run was 270 ns of NPT simulation for each system. Coordinates
were stored every 3 ps resulting in 90 000 snapshots in total for each
Nt system.
J(L) = βF(z)
. (20) The number of membrane crossings was determined with
hLσD⟨e ⟩ ⟨e−βF(z) ⟩
L L the Rickflow package code using dividing surfaces at |z| = hm /2
= 2.3 nm.29 The diffusion tensor with diagonal elements Dxx , Dyy ,
This highlights that the same number of crossings will be encoun- Dzz was computed with Gromacs.
tered when a profile F(z) is replaced by −F(z), even though the
permeability changes. While it may be counterintuitive, a profile B. Numerical results
with a high free energy barrier will therefore give the same sampling
The free energy profile was determined from the histogram p(z)
efficiency as a profile with a low free energy barrier. Highly solu-
of the position of the CG beads as F(z) = −kB T ln p(z). Figure 4
ble membranes are thus not necessarily improving the statistics in
presents F(z) for three exemplary bead types. Bead AC2 has a deep
the counting method. This realization is important for drug design,
free energy well; hence, it is highly soluble in the membrane (K ≫ 1).
where substances with high permeability might be trapped by the
membrane thus reducing their efficacy.

IV. ILLUSTRATION WITH COARSE GRAINED

MOLECULAR DYNAMICS
To illustrate the discussion of Sec. III, a membrane of dipalmi-
toylphosphatidylcholine (DPPC) was simulated at the molecular
scale using a coarse grained (CG) model. A range of different per-
meants was considered with both positive ΔF (barrier) and negative
ΔF (well): 16 bead types of the Martini coarse-grained model were
used as the permeating molecule. These beads span a wide range of
free energy profiles F(z) across the membrane.

A. Simulation details
Using the Gromacs-2019.3 package,38 simulations were per-
formed on 16 systems consisting of a phospholipid membrane, CG
beads, and water. Each system contained 20 permeant beads of one
specific CG type, and it was solvated in 2000 water molecules. The FIG. 4. Free energy of three selected bead types that are representative of a free
membrane consists of 128 DPPC lipid molecules in each simulation, energy barrier (P1), a free energy well (AC2), and a mix of both (SN0). The profile
and its initial coordinates were taken from the MARTINI website.39 F(z) = −k B T ln p(z) is plotted with p(z), the histogram of the 20 permeants’ z-
coordinates over all 90 000 snapshots.
The topology file and initial coordinates of the permeant beads were

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-6

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

Bead P1 has a high free energy barrier and is insoluble in the mem- the MD data could be obtained for a small water layer of thickness
brane (K ≪ 1). Bead SN0 shows a mix of the free energy well and d = 0.1 nm, which was added to the plot with a gray open sym-
barrier, where the solubility depends on the exact location in the bol. The spread of the vertical lines however indicates that almost
membrane. The other bead types had similar profiles with different any value can be the outcome by tuning the parameter d. This
elevations. indicates again that the number of crossings, and hence the per-
Figure 5 shows the number of detected crossings vs ΔF for the meability, is very sensitive to the considered thickness. One should
16 beads (full symbols), where ΔF is computed as the free energy therefore be cautious when interpreting and citing experimental or
difference between the membrane center and the membrane bor- computational results for highly permeable membranes.
der. The graph follows the shape of the curves in Fig. 2(a), where For Eq. (15), the integral is evaluated using the measured
ncross generally drops with |ΔF|. These MD values are compared to free energy profile F(z), the same constant diffusivity, and d ≈ 0
the model values of ncross (open symbols) based on the flux in the [not exact because of binning in p(z)], giving the values indi-
compartmental model [Eq. (11)] or the more general formulation cated with a colored open symbol. Compared to the compartmental
of the Smoluchwoski model [Eq. (15)]. The model parameters were model, this more general model uses a position-dependent solu-
inferred from the molecular dynamics simulations. bility, which allows us to describe more details of the membrane.
For Eq. (11), K is computed as exp(−βΔF), hm = 4.6 nm cor- Despite the crudeness in the parameter estimations, a good agree-
responding to the membrane thickness, the average box length L = ment is found between these model values and the numerical MD
10 nm, and T sim = 270 ns. The normal diffusivity Dzz describes the data in Fig. 5. The permeants in the present simulations are of course
diffusion normal to the membrane surface, but the measured Dzz only a single bead without internal degrees of freedom, and the
already incorporates the effects of the barrier or well. Therefore, the membrane contains only one phospholipid type. For more realis-
lateral diffusivity D∥ = (Dxx + Dyy )/2 was taken as the diffusivity D tic membranes and realistic permeants, the compartmental model
in the model since diffusion parallel to the membrane is not hin- or the more general formulation [Eqs. (11) and (15)] might not
dered by free energy barriers/wells. The water thickness d was varied have a sufficient level of detail to match the MD results. An exten-
from 0 to (L − hm )/2, giving a range of crossings that is indicated sion of the model in Fig. 1 could for instance be a model that
as a dashed vertical line on the graph. The lower end of the line has both a free energy well and barrier. However, the complexity
corresponds to the number of crossings through the whole simu- of more realistic simulations not only comes from the inhomo-
lation box (h = L), while the higher end corresponds to the number geneity in F and D but there might also be memory effects that
of crossings through the membrane only. For ΔF ≪ 0, the number make the kinetics non-Markovian, while the Smoluchowski equa-
of crossings is highly sensitive to the considered membrane thick- tion is based on the Markovian assumption. Another challenging
ness. For ΔF ≫ 0, the number of crossings is no longer sensitive to situation is when the permeation is accompanied by hysteresis,
the water layer thickness (lines mostly hidden by symbols), in accor- which is caused by not observing energy barriers in some orthog-
dance with Fig. 2(a). It was found that a decent correspondence with onal degree(s) of freedom (orthogonal to the permeation direction
z). This is for instance the case when the permeant has an inter-
nal degree of freedom that is changed during the permeation, e.g.,
a dihedral angle in an oligopeptide, or an orientational degree of
freedom, e.g., the ring plane in ibuprofen. A two-dimensional or
higher-dimensional diffusive equation should then be constructed,
which captures all reaction coordinates that are relevant for the
permeation.
In these more advanced approaches, based on more elabo-
rate multi-dimensional models with kernels, the free energy would
still be one of the main parameters. The present analysis based on
a Smoluchowski equation can therefore be seen as a decent first
approximation to gain insight into the sampling efficiency. The
presented curves clearly illustrate that a higher variation in the
free energy profile results in lower sampling efficiency. Both deep
wells and high barriers cause a dramatic drop in crossings, and the
counting method for permeability will perform inaccurately in those
cases.
When the counting method suffers from poor statistics and
running much longer MD simulations is not expected to be a
sufficient measure, a rare event method for rate calculations can
FIG. 5. Number of crossings ncross as a function of ΔF for the 16 permeant types. be considered, such as replica exchange transition interface sam-
Filled symbols: value as measured in the MD simulations with color according to pling (RETIS)42–44 or milestoning.45 In recent work, the connection
the bead type, i.e., well (blue), barrier (red), and mix (orange). Gray open sym- between RETIS and the permeability was derived,37 and this method
bols: value based on the compartmental model [Eq. (11)] with d = 0.1 nm. Vertical can be a valid alternative for the counting method.46 Another train
dashed lines show the range between setting d = 0 (top of the line) and d = (L of thought is the use of biased simulations, where a large vari-
− hm )/2 (bottom of the line). Colored open symbols: value based on the general ety of free energy is canceled out artificially by adding a bias
Smoluchowski model [Eq. (15)].
potential.20,47,48

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-7

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

3
V. CONCLUSION Y.-x. Shen, W. Song, D. R. Barden, T. Ren, C. Lang, H. Feroz, C. B. Henderson,
P. O. Saboe, D. Tsai, H. Yan, P. J. Butler, G. C. Bazan, W. A. Phillip, R. J. Hickey,
This paper investigated the efficiency of counting membrane P. S. Cremer, H. Vashisth, and M. Kumar, Nat. Commun. 9, 3304 (2018).
crossings as a method to determine the permeability. A theoretical 4
V. A. Poteryaeva, M. A. Bubenchikov, A. M. Bubenchikov, and A. V. Lun-Fu, Sci.
relation for the flux was derived based on the ISD model using the Rep. 10, 15631 (2020).
compartmentalized toy system in Fig. 1. An additional water layer 5
A. Mandal, V. Agrahari, V. Khurana, D. Pal, and A. K. Mitra, Expert Opin. Drug
was included in the model to show the sensitivity to the considered Delivery 14, 385 (2017).
6
thickness, which causes a large uncertainty for high permeability B. J. Bennion, N. A. Be, M. W. McNerney, V. Lao, E. M. Carlson, C. A. Valdez,
membranes. It was shown that a high permeability in itself is not M. A. Malfatti, H. A. Enright, T. H. Nguyen, F. C. Lightstone, and T. S. Carpenter,
J. Phys. Chem. B 121, 5228 (2017).
a qualified indicator of whether the sampling efficiency of the sim- 7
L. Galluzzi, I. Vitale, S. A. Aaronson, J. M. Abrams, D. Adam, P. Agostinis, E. S.
ulation box will be extensive enough for the counting method to
Alnemri, L. Altucci, I. Amelio, D. W. Andrews, and M. Annicchiarico-Petruzzelli,
be reliable, and caution should be taken for permeants whose sol- Cell Death Differ. 25, 486 (2018).
ubility varies greatly over the simulation box. A practical effect of 8
W. K. Subczynski, M. Pasenkiewicz-Gierula, R. N. McElhaney, J. S. Hyde, and
this realization in, e.g., drug design is that substances with high per- A. Kusumi, Biochemistry 42, 3939 (2003).
meability might be trapped by the membrane, thus reducing their 9
J. Widomska, M. Raguz, and W. K. Subczynski, Biochim. Biophys. Acta,
efficacy. Biomembr. 1768, 2635 (2007).
10
Moreover, the composition rule for permeabilities was shown M. Möller, J. Lancaster, and A. Denicola, Curr. Top. Membr. 61, 23 (2007).
11
to hold for fluxes as well. The limiting behavior of the number of J. Godyń, D. Gucwa, T. Kobrlova, M. Novak, O. Soukup, B. Malawska, and
crossings was determined for a high barrier membrane and a per- M. Bajda, Talanta 217, 121023 (2020).
12
meant trap membrane, and the dependence on the considered water J. W. Choi, S. H. Bae, M. Kwak, T. G. Lee, M. B. Heo, and D. W. Lee, Sens.
Actuators, B 321, 128624 (2020).
layer was discussed. It was also shown that the sampling efficiency 13
P. Bacchin, D. Snisarenko, D. Stamatialis, P. Aimar, and C. Causserand,
of the complete simulation box remains unaffected when the free J. Membr. Sci. 614, 118485 (2020).
energy profile F(z) is inverted to −F(z), if D is a constant. These theo- 14
C. Lin, T.-C. Kuo, J.-C. Lin, Y.-C. Ho, and F.-L. Mi, Int. J. Biol. Macromol. 160,
retical considerations were all made under the assumption of purely 558 (2020).
diffusive kinetics since the ISD model is based on the Smoluchowski 15
E. Awoonor-Williams and C. N. Rowley, Biochim. Biophys. Acta, Biomembr.
equation for Markovian kinetics. 1858, 1627 (2016).
16
Besides the theoretical assessment of the counting method, 16 N. Pokhrel and L. Maibaum, J. Chem. Theory Comput. 14, 1762 (2018).
17
systems were simulated by coarse-grained MD to illustrate the the- R. M. Venable, A. Krämer, and R. W. Pastor, Chem. Rev. 119, 5954 (2019).
18
ory. The permeants’ crossings were detected and the membrane’s W. Shinoda, Biochim. Biophys. Acta, Biomembr. 1858, 2254 (2016).
19
free energy profile was characterized as a free energy well, a barrier, C. Hoffmann, A. Centi, R. Menichetti, and T. Bereau, Sci. Data 7, 1–7
or a mix of these two. In both the cases of a high free energy barrier (2020).
20
and a deep free energy well, it can be concluded that the permeabil- C. T. Lee, J. Comer, C. Herndon, N. Leung, A. Pavlova, R. V. Swift, C. Tung,
ity calculation by counting crossings will face failure because of the C. N. Rowley, R. E. Amaro, C. Chipot, Y. Wang, and J. C. Gumbart, J. Chem. Inf.
Model. 56, 721 (2016).
insufficient number of crossing during the simulation. To overcome 21
S.-J. Marrink and H. J. C. Berendsen, J. Phys. Chem. 98, 4155 (1994).
this problem, one should resort either to running longer MD simu- 22
S. J. Marrink and H. J. C. Berendsen, J. Phys. Chem. 100, 16729 (1996).
lations to observe more crossings using larger simulation boxes with 23
K. Gaalswyk, E. Awoonor-Williams, and C. N. Rowley, J. Chem. Theory
more permeants or implementing other simulation methods that are Comput. 12, 5609 (2016).
applicable for rare events, e.g., RETIS. 24
A. Ghysels, R. M. Venable, R. W. Pastor, and G. Hummer, J. Chem. Theory
Comput. 13, 2962 (2017).
25
O. De Vos, R. M. Venable, T. Van Hecke, G. Hummer, R. W. Pastor, and
A. Ghysels, J. Chem. Theory Comput. 14, 3811 (2018).
ACKNOWLEDGMENTS 26
H. A. L. Filipe, M. Javanainen, A. Salvador, A. M. Galvão, I. Vattulainen, L. M. S.
The computational resources (Stevin Supercomputer Infras- Loura, and M. J. Moreno, J. Chem. Theory Comput. 14, 3840 (2018).
27
tructure) and services used in this work were provided by the R. J. Ferreira and P. M. Kasson, ACS Infect. Dis. 5, 2096 (2019).
28
VSC (Flemish Supercomputer Center), funded by Ghent University, H. Yang, M. Zhou, H. Li, T. Wei, C. Tang, Y. Zhou, and X. Long, ACS Omega
FWO, and the Flemish Government—department EWI. 5, 4798 (2020).
29
A. Krämer, A. Ghysels, E. Wang, R. M. Venable, J. B. Klauda, B. R. Brooks, and
R. W. Pastor, J. Chem. Phys. 153, 124107 (2020).
30
R. J. Dotson, C. R. Smith, K. Bueche, G. Angles, and S. C. Pias, Biophys. J. 112,
DATA AVAILABILITY 2336 (2017).
31
O. De Vos, T. Van Hecke, and A. Ghysels, Adv. Exp. Med. Biol. 1072, 399
The data that support the findings of this study are available (2018).
from the corresponding author upon reasonable request. 32
A. Ghysels, A. Krämer, R. M. Venable, W. E. Teague, E. Lyman, K. Gawrisch,
and R. W. Pastor, Nat. Commun. 10, 5616 (2019).
33
R. J. Dotson and S. C. Pias, Adv. Exp. Med. Biol. 1072, 405 (2018).
34
M. Ghorbani, E. Wang, A. Krämer, and J. B. Klauda, J. Chem. Phys. 153, 125101
REFERENCES
(2020).
1 35
B. Smit and T. L. M. Maesen, Chem. Rev. 108, 4125 (2008). Q. Al-Awqati, Nat. Cell Biol. 1, E201 (1999).
2 36
M. Tagliazucchi, Chemically Modified Nanopores and Nanochannels (William M. N. Möller, Q. Li, M. Chinnaraj, H. C. Cheung, J. R. Lancaster, Jr., and
Andrew Publishing, 2017), Chap. 1, pp. 1–25. A. Denicola, Biochim. Biophys. Acta, Biomembr. 1858, 2923 (2016).

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-8

Published under license by AIP Publishing
 The Journal
ARTICLE scitation.org/journal/jcp
of Chemical Physics

37 43
A. Ghysels, S. Roet, S. Davoudi, and T. van Erp, “Exact non-Markovian perme- A. Lervik, E. Riccardi, and T. S. van Erp, J. Comput. Chem. 38, 2439 (2017).
44
ability from rare event simulations,” Phys. Rev. X (unpublished) (2021). E. Riccardi, A. Lervik, S. Roet, O. Aarøen, and T. S. Erp, J. Comput. Chem. 41,
38
M. J. Abraham, T. Murtola, R. Schulz, S. Páll, J. C. Smith, B. Hess, and E. Lindahl, 370 (2020).
45
SoftwareX 1-2, 19 (2015). L. W. Votapka, C. T. Lee, and R. E. Amaro, J. Phys. Chem. B 120, 8606 (2016).
39 46
S. J. Marrink, A. H. De Vries, and A. E. Mark, J. Phys. Chem. B 108, 750 (2004). E. Riccardi, A. Krämer, T. S. van Erp, and A. Ghysels, J. Phys. Chem. B 125(1),
40
S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman, and A. H. De Vries, 193–201 (2021).
47
J. Phys. Chem. B 111, 7812 (2007). E. Darve and A. Pohorille, J. Chem. Phys. 115, 9169 (2001).
41 48
P. Mark and L. Nilsson, J. Phys. Chem. A 105, 9954 (2001). J. Comer, K. Schulten, and C. Chipot, J. Chem. Theory Comput. 10, 554
42
T. S. van Erp, Phys. Rev. Lett. 98, 268301 (2007). (2014).

J. Chem. Phys. 154, 054106 (2021); doi: 10.1063/5.0033476 154, 054106-9

Published under license by AIP Publishing