Two Relations to Estimate Membrane Permeability Using Milestoning.

Prediction of passive permeation rates of solutes across lipid bilayers is important to drug design, toxicology, and other biological processes such as signaling. The inhomogeneous solubility-diffusion (ISD) equation is traditionally used to relate the position-dependent potential of mean force and diffusivity to the permeability coefficient. The ISD equation is derived via the Smoluchowski equation and assumes overdamped system dynamics. It has been suggested that the complex membrane environment may exhibit more complicated damping conditions. Here we derive a variant of the inhomogeneous solubility diffusion equation as a function of the mean first passage time (MFPT) and show how milestoning, a method that can estimate kinetic quantities of interest, can be used to estimate the MFPT of membrane crossing and, by extension, the permeability coefficient. We further describe a second scheme, agnostic to the damping condition, to estimate the permeability coefficient from milestoning results or other methods that compute a probability of membrane crossing. The derived relationships are tested using a one-dimensional Langevin dynamics toy system confirming that the presented theoretical methods can be used to estimate permeabilities given simulation and milestoning results.

Introduction

Passive transport of solutes across membranes is of key importance for processes such as the uptake and excretion of drug compounds. Poor permeability can cause unfavorable pharmacokinetics or bioavailability leading to drug attrition.[1] From this perspective of drug design, methods which provide molecular level insight into the permeation process are of high interest. As such, many physical models of permeability have been developed; these methods are recently reviewed.[2]

Mathematically, the permeability coefficient P of a solute is described by relating the flux per unit area, J, across a bilayer and the concentration gradient, Δu, asFrom this relation, assuming overdamped conditions and accounting for drift caused by the potential of mean force (PMF) we get the inhomogeneous solubility-diffusion equation (ISD)[3−7]where R is the resistivity, β = 1/kBT where T is the temperature, kB is Boltzmann’s constant, W(z) is the PMF, D(z) is the local diffusivity constant and z1, z2 are two positions just outside of opposing sides of the membrane. Estimates of W(z) and D(z) can be predicted using molecular dynamics (MD) simulations.[2,6,7] While the ISD has been applied successfully in many studies,[2,5,8−22] recent work by various groups has suggested that the dynamics of crowded environments, such as inside the bilayer, may exhibit more complicated behaviors.[2,23] The dynamics inside of the bilayer has not, to the best of our knowledge, been probed in detail.[2]

The theory of milestoning is a promising new approach to study molecular processes by employing many unbiased short independent simulations in subregions of phase space to predict thermodynamic and kinetic properties of interest. Salient descriptions of milestoning theory can be found in refs (24−26). Previously, milestoning has been applied to the permeability problem by Cardenas et al., studying the permeation of block tryptophan across a membrane.[27] Cardenas et al. later studied the permeation of small peptide NATA across a DOPC bilayer using milestoning. In their later work, they proposed a scheme to estimate the permeability coefficient.[28,29] This scheme is comprised of three steps involving (1) estimating the diffusivity of the small molecule in bulk, (2) the measurement of the steady state flux from bulk into the membrane, and (3) the estimation of permeant crossing probability. We will show that the additional calculations in step 2 are not necessary and the permeability can be computed directly from the milestoning results and bulk diffusivity.

In this work, we derive two methods for obtaining permeability coefficients based on milestoning results. The first derivation is a new formulation of the ISD relating the mean first passage time (MFPT) to the permeability coefficient P assuming overdamped dynamics. The second derivation makes no assumptions about system damping and relates P to the transition probability of crossing a membrane over diffusing to a previous position in bulk solvent. We evaluate the compatibility of using either approach depending on damping dynamics using a one-dimensional (1D) Langevin dynamics system across various W(z) and D(z) profiles, showing that the new relations yield solutions in good accord with previous results.

Theory

Relating Permeability to Mean First Passage Time via the Smoluchowski Equation

The physical description of permeability can be drawn as a membrane separating donor and acceptor compartments (Figure ). Assuming that the dynamics are overdamped, permeation across a membrane can be modeled as a one-dimensional process using the 1D Smoluchowski equationwhere u(z,t) is the probability distribution at position z and time t, v(z,t) is the drift velocity, and J(z,t) is the infinitesimal flux. A particle travels through a viscous medium at velocity v = F(z)/γ, where F is the force on the particle and γ is the friction coefficient. Using the Smoluchowski-Einstein relation γ–1 = βD, assuming steady-state conditions, and applying that , we can convert the partial differential equation into an ordinary differential equationIntegrating and multiplying by , we getwhere A1 is a constant. The RHS can be condensed into a single differentialBy integrating, we can obtain a general expression for the steady-state concentrationWe now apply the boundary conditions shown in Figure . Assuming that the acceptor compartment is a sink, u(a) = 0, and that in bulk W(a) = 0Solving for A2

Figure 1

Schematic of membrane permeability showing the concentration of permeant with respect to position. We define that the concentration at z = b is u0 and z = a is 0. Assuming overdamped dynamics along with these boundary conditions, we get the inhomogeneous solubility-diffusion equation and its relation to passage time.

Substituting eq into eq , we get a difference of two integrals. Combining them into a single integral and setting the concentration of the donor compartment at position b to that of bulk u0We can solve for the constant A1Combining eqs , 9, and 7, the concentration u(z) can be expressed according toReturning to the calculation of the stationary flux, substituting eq into eq , and recognizing that the flux is just the negative of our first constant of integration, we get thatSubstituting into eq and assuming that Δu = u0, we get the ISD

Mean First Passage Time

The MFPT ⟨t⟩ of a diffusional encounter can be generally expressed as[30]here V is the total volume of the system in question, m is the total number of isolated sinks, and A is the total area of the ith sink. In our 1D case, the expression is simplifiedTaking the integral of eq and reversing the order of integration in the numerator, we getNoting the equivalence of ∫ ∫f(x)dx dz = ∫ ∫f(x)dx dz for symmetrical functions over z|b < z < a and dividing by J(a) = −A1 we getThe MFPT can be related to the permeability coefficient by integrating eq by partswhere the second term on the RHS has been shown to be equal to the MFPT by others[10,27,31,32] and can be derived by applying Fubini’s theorem. Rewriting in terms of P and ⟨t⟩and rearranging we get the permeability as it relates to ⟨t⟩This is formally equivalent to the ISD and obviates the calculation of the local diffusivity, requiring only the local PMF (or equivalently, the local stationary concentration of the permeant) and the MFPT.

Permeability from Milestoning Using Crossing Probability

Our second permeability derivation in this section is very similar to that used in Brownian dynamics (BD) simulations that make use of the Northrup-Allison-McCammon (NAM) algorithm.[33] As in BD simulations, the b-surface is the starting position of particles in the simulations. Conversely, the q-surface and a-surface are absorbing boundaries during the simulations. The quantity ρ is the proportion of permeants that start at the b-surface and cross the membrane to touch the a-surface and not “escape” to the q-surface. Here, we attempt to be consistent with notation used in the NAM derivation. The only exception is that NAM uses a β to represent the probability of crossing the a-surface, but we use the symbol ρ instead to avoid confusion with the thermodynamic β used in the previous section. This probability ρ can be found using milestoning.

We define one additional surface in bulk solution on the donor side of the membrane: 0-surface (placed at z = 0 for mathematical convenience). Note that a > b > q > 0. A steady-state Smoluchowski problem (Figure ) is set up in the bulk domain between the 0-surface and the b-surface, where the concentration at the 0-surface boundary is held at u0 and the concentration at the b-surface is held at 0. Any flux across the b-surface has a probability equal to ρ that it will be taken across the membrane and 1 – ρ to return to the q-surface. Under this scheme, the permeability from the 0-surface to a-surface can be found according towhere u0 is the difference between the concentration at the 0-surface and the concentration at the a-surface, and J(a) is the flux across the a-surface. Because we assume steady-state conditions, the flux across the a-surface must be equal to the flux across the 0-surface

Figure 2

Shown here is the milestoning system setup. The stationary concentrations in the bulk u1 and u2 are subject to boundary conditions at z = 0 and z = b. A permeant that hits the b-surface has a probability of ρ to cross the membrane and a probability 1 – ρ to return to the q-surface. These boundaries and domains are not necessarily drawn to scale.

Because the region between the 0-surface and b-surface exists in bulk solution, we assume that the PMF is constant and the dynamics overdamped, and we can use a 1D Smoluchowski equation for a continuum description of this region. Although the q-surface is an absorbing boundary during the simulations and milestoning calculations, the q-surface will become a source boundary in the continuum Smoluchowski calculations for the region where 0 < z < b. Here, it is useful to split this domain into two pieces: 0 < z < q, q < z < b. We then solve a system of two Smoluchowski equations. Domain 1:Domain 2:Domains 1 and 2:where D is the diffusion coefficient of the particle in bulk water. The last equation, eq , assumes that the flux through the 0-surface (J1) is less than the flux going through the b-surface (J2) because some of that flow comes back and appears once again on the q-surface. Because anything that goes through the b-surface has a probability of 1 – ρ to return to the q-surface, the flow out of the q-surface is what flows through the b-surface (J2) multiplied by 1 – ρ. This is equal to the difference between J2 and J1. We find a general solution to eqs and 24d by integrating twicewhere c1, c2, d1, and d2 are constants of integration. By eqs and 24d we know that J = D(−du/dz), and therefore J1 = −Dc1 and J2 = −Dc2 and thus eq simplifies to c1 = ρc2. Upon applying the boundary conditions from eqs and 24f to eq we get that d1 = u0 and d2 = −c1b/ρ. Substituting the above relations and eq into 24g and solving for c1 we getRearranging for the flux, J1, and dividing by u0, we get the permeability from the 0-surface to the a-surfaceNow we remove the contribution to the permeability by the bulk layer, 0 < z < b, to obtain the true permeability across the membrane P. This is done most effectively by working with the resistivity, R = 1/P. Expressing the resistivities across the domain as a sum of resistivities over subdomainsAssuming a constant PMF and diffusivity in bulk, 0 < z < b, the resistivity in that domain can be expressed asobtained by solving the Smoluchowski equation for a flat PMF and diffusion profile. Plugging eq into eq and inverting we get the true membrane permeability as related to the crossing probability ρFor subsequent discussion, we refer to eq as the MFPT in ISD relation (MFPT-ISD), and eq as the permeability based on crossing probability (PBCP).

Methods

Langevin Dynamics Model

To demonstrate the validity of the above relations, we have developed a 1D Langevin dynamics model to generate sample trajectories over a user specified PMF and viscosity profile. All code used has been made available through GitHub at https://github.com/ctlee/langevin-milestoning. The profiles are represented using interpolated cubic Hermite polynomials shown in Figure . Four different systems are considered, flat, small-barrier, urea-like barrier, and codeine-like. The flat system emulated diffusion in bulk with viscosity chosen to match water. For the small-barrier system, an arbitrary barrier height was chosen such that brute force computation could yield successful crossings. Meanwhile the urea- and codeine-like systems are interpolated from profiles from ref (34). To facilitate convergence of calculations, the viscosity of the membrane environment is set to 0.005 kg m–1 s–1. Each of the membrane simulations employ the same viscosity profiles. The hydrodynamic radius and mass of the particle were arbitrarily chosen for the flat and small-barrier cases while the urea and codeine cases use molecular weight and a radius as predicted by software HydroPro.[35] The parameters used in each case are shown in Table . The diffusivity along the profile, D, can be calculated from the Stokes–Einstein relationwhere η is the viscosity, Rhyd is the hydrodynamic radius, and c is the friction coefficient. Notably, new systems can be generated with any user-specified PMF, viscosity, mass, and hydrodynamic radius parameters of interest.

Figure 3

PMF and diffusivity profiles used for each system. Each profile is represented using a piecewise cubic Hermite polynomial interpolated from a set of user specified values.

Table 1

Parameter Inputs to the Langevin Dynamics Simulations for Each Case

	flat	small barrier	urea	codeine
Rhyd (Å)	5	5	2.86	4.32
mass (g/mol)	20	20	60	299
T (K)	298	298	298	298
brute dt (s)	1 × 10–12	1 × 10–12	1 × 10–12	1 × 10–12
milestone dt (s)	1 × 10–15	1 × 10–15	1 × 10–15	1 × 10–15

We integrate the Langevin equation using a velocity–Verlet-like integrator by Grønbech-Jensen–Farago[36] to generate sample trajectorieswhere q is the position, t is the current time, c is the friction coefficient, b and a are unitless constants, f is the force, v is the velocity, and m is the particle mass. This integrator was derived by integrating the noise term at a small time step dt. The new term, Γ, is Gaussian distributed and has the following propertiesThe quality of random number generators (RNG) to sample Γ is shown in Figure S1. In this work, the aforementioned drag coefficient, c, is estimated by Stokes lawwhere Fd is the frictional force. For additional information regarding the implementation details of the Langevin dynamics engine, please refer to the SI.

Using the above dynamics engine, sample trajectories under various conditions can be generated. Here the permeability is calculated via the following several methods: (1) theoretical ISD, (2) theoretical MFPT-ISD, (3) brute forced MFPT-ISD, (4) milestoning MFPT-ISD, (5) brute PBCP, and (6) milestoning PBCP. Methods 1 and 2 are directly evaluated numerically from the interpolated profiles using the relations in the Theory section. Methods 3–6 employ brute force or milestoning calculations to generate sample trajectories that are postprocessed to yield the permeability coefficient. The strategy for brute force and milestoning simulations are described below.

Brute Force Sampling and Statistics

Brute force attempts to sample the membrane crossing event via a single continuous trajectory were attempted as follows: Trajectories to sample ⟨t⟩ and ρ were initiated at one end of the membrane z = −25 Å, the a-surface. Trajectories to sample the MFPT employed reflecting boundary conditions at the start to prevent the particle from diffusing away into bulk. Effectively trajectories, which passed the reflecting boundary at the a-surface, were placed back across the boundary and the velocity was reversed. In contrast, the simulations to obtain ρ have an absorbing boundary placed at z = −26 Å, corresponding to the so-called q-surface. Trajectories that crossed the q-surface were halted and counted as noncrossing events.

For all brute force trajectories, an absorbing boundary was placed at z = 25 Å. Upon crossing this terminal boundary, the simulations are halted and either the crossing time or the event recorded for the MFTP, and ρ sampling cases, respectively. The crossing probability was found according towhere n and n are the number of b- or q-surface crossing events observed. The FPT distribution often takes the form of a Poisson or exponential distribution. We attempt to quantify the quality of the MFPT estimate by computing the root-mean-square error (RMSE) of the population and confidence interval by bootstrapping. The RMSE is defined aswhere t is the time that trajectory i took before touching the absorbing boundary and N is the total number of observations. We also report the 95% confidence intervals of both the MFPT and the crossing probability. This is estimated by bootstrapping using the bias-corrected and accelerated method and percentile method for the MFPT and crossing probabilities respectively with 10 000 resamples.[37,38]

Milestoning

In this section, we explain our experimental approach and include a brief description of milestoning theory. Milestones are spaced at 1 Å increments across the bilayer ([−25,25]), and in the case of ρ estimation, an additional milestone was placed at z = −26 Å. Trajectories were initiated with a position starting at each milestone, M, with a random velocity sampled from the Maxwell–Boltzmann distribution. Two stages are needed to obtain milestoning statistics. According to formal milestoning theory, we must determine the first hitting point distribution (FHPD) on M. To find the FHPD, trajectories are started at milestone M and run until they cross either another milestone (M, M) or the originating milestone (M). If the trajectory is self-crossing then the simulation is halted and statistics discarded. Conversely, if the trajectory crosses an adjacent milestone first then this trajectory is a member of the FHPD on milestone M. Because of the microscopic time-reversibility of classical mechanics, the reverse of this trajectory is a valid simulation of the crossing from (M, M → M). Thus, we term this stage the “reverse” phase. Next we begin the “forward” phase where each member of the FHPD is restarted at M but with the initial velocity vectors reversed. These trajectories propagate until they cross an adjacent milestone, M or M. Self-crossing events of the starting milestone M are ignored in the forward phase.

The forward phase provides all the transition and incubation time statistics that will be used in the milestoning analysis. As the milestones are crossed, each trajectory is counted, and the time that each trajectory took is also tracked. A transition kernel matrix K is then generated, whose elements K represent the transition probability that a system starting at milestone M will transition to milestone Mwhere n represents the total number of trajectories that cross from milestone i to milestone j, where N is the index of the final milestone. Similarly, a lifetime or incubation time vector τ is also generated, whose elements τ represent the average time that a system spends in milestone iwhere t is the time that trajectory l takes to reach an adjacent milestone, and n is the total number of trajectories starting from milestone i.

Using K and τ, we can compute thermodynamic quantities such as the stationary flux (qstat) and the stationary probability (pstat), as well as the MFPT (⟨t⟩), a kinetic quantity. First, we compute the stationary probability using the following eigenvalue equation

Note that before solving eq , we typically have to modify K to have the correct boundary conditions. For our membrane system, because the two terminal milestones at each end of the membrane have no additional milestones beyond them, it does not make sense to simulate from them. We can fill in reasonable transition probabilities and lifetimes using the following scheme: Because each terminal milestone only has one adjacent milestone, we can set the transition probability to transition to the singular neighbor or the other side of the membrane with equal probability. These are periodic boundary conditions. That is, for a set of milestones M = M1, M2, ..., Mandwhere N is the index of the last milestone and i is a milestone index. Splitting the transition probabilities equally makes sense because these milestones are in the bulk solvent, where the PMF is flat. For the lagtime vector τ, given that the terminal milestones are in bulk solvent, the lagtime is defined as followswhere Δx is the Cartesian distance to the adjacent milestone, and D is the diffusivity coefficient in bulk. With these changes, we can estimate the stationary probability according towhere we take the element-wise product over all milestone indeces i.

The stationary probability distribution, pstat, can be used to obtain an estimate of the underlying PMF, W, at milestone i accordinglyWe use this calculated PMF to estimate the integral in eq by performing a numerical integration.

To find the MFPT needed to estimate the permeability in eq , we must modify the transition kernel K once again. This time, we must make milestone M a “draining” state by setting all K = 0, ∀k. We also must modify transitions at milestone M1 to transition only to its neighbor, milestone M2:To get the MFPT, we then solve the following equationwhere p0 is a starting distribution of probabilities along the milestone. We choose p0, to be 1 if i = 1 and choose p0, to be 0 otherwise.

Finally, we obtain the crossing probability ρ for eq by further modifying the transition kernel K. Recalling that the computation of ρ requires an additional milestone at q-surface, which we assign index i = 0, we simulate starting from the other interior milestones. The first and last milestones become “sink” states that capture probability and never let it escape; thus no simulations start from sink states, they only end on them. Sink states are defined by settingandOnce again solving eq we obtain a stationary flux vector, qstat, which represents the proportion of systems that would have partitioned between the sink states at the ends of the membrane. The elements of the vector qstat, must be normalized such that all its elements sum to one for this calculation, however, because K is a Markov matrix, the eigenvector should already be normalized. Because milestone M represents the a-surface, ρ in eq can be extracted from the last value in qstat

The error is sampled according to the scheme found in ref (39). For all milestoning error values, we resample 1000 times and report the standard deviation of each measurable.

Results and Discussion

Below we report the computed permeability coefficients and MFPTs according to the six methods: (1) the ISD, (2) ISD from MFPT, (3) ISD with brute forced MFPT, (4) milestoning estimated ⟨t⟩ in ISD, (5) brute force with no assumptions, and (6) milestoning with no assumptions in Table . Note again that the brute-forced MFPT is performed with a reflecting boundary condition at the starting position while no such condition is applied for the brute force crossing probabilities. This explains why we observed no crossings in the brute force MFPT urea system.

Table 2

Computed MFPT and Permeability Coefficients According to Each Methoda

	flat	small barrier	urea	codeine
theoretical MFPT (s)	9.19 × 10–9	1.11 × 10–8	1.99 × 10–3	2.89 × 10–7
brute forced MFPT (s)	9 ± 7 × 10–9 (5100)	1.1 ± 1.0 × 10–8 (2500)	3 ± 3 × 10–4 (1694)	10 ± 9 × 10–7 (1600)
milestoning MFPT (s)	1.15 ± 0.05 × 10–8	1.35 ± 0.04 × 10–8	1.9 ± 0.3 × 10–3	4.0 ± 0.3 × 10–7
crossing probability (ρ)
brute force	2.7 ± 0.8 × 10–2 (1722)	1.0 ± 0.4 × 10–2 (2500)	N/A (1627)b	7 ± 1 × 10–2 (1340)
milestoning	1.7 ± 0.1 × 10–2	9.8 ± 0.5 × 10–3	3.1 ± 0.4 × 10–8	2.2 ± 0.3 × 10–2
permeability log(cm/s)
(1) ISD	1.43	1.14	–4.33	1.44
(2) MFPT-ISD	1.43	1.14	–4.33	1.44
(3) brute force MFPT-ISD	1.43 ± 0.01	1.14 ± 0.01	–3.49 ± 0.02	0.91 ± 0.02
(4) milestoning MFPT-ISD	1.33 ± 0.03	1.06 ± 0.02	–4.3 ± 0.09	1.35 ± 0.06
(5) brute force PBCP	1.6 ± 0.1	1.1 ± 0.2		2.04 ± 0.09
(6) milestoning PBCP	1.37 ± 0.03	1.13 ± 0.02	–4.13 ± 0.06	1.54 ± 0.04

Reported error bars for the MFPT corresponds to the RMSE while the error for the crossing probability is from bootstrapping. Error in the brute force log permeability is estimated by the mean log(P) difference of the high and low confidence intervals.

We observed no brute force crossing transitions for the urea system. The values in parentheses represent the number of forward-phase trajectories sampled in each case.

The distribution of the mean first passage times is shown in Figure . Because of the highly tailed nature of the MFPT, it is difficult to estimate it from brute force simulations. Milestoning greatly aids in the computation of the first passage times by reducing the total distance between points of interest. Using clever statistics, milestoning theory is capable of reconstructing the MFPT from many short simulations. For all cases, we obtain MFPTs in good agreement with the theoretical results. Deviations in both the brute force and milestoning cases are likely caused by insufficient sampling. Sampling convergence for the milestoning can be visualized by deviations in the milestoning-derived PMF from the actual PMF.

Figure 4

Shown here are the distributions of the brute force MFPTs for systems (a) flat, (b) small-barrier, (c) urea, and (d) codeine. Because of presence of large barriers in the urea and codeine cases, the sampling is limited. The red points and error bars indicate the mean and boostrapped 95% confidence intervals of the MFPT.

Reconstructed PMF profiles from milestoning are shown in Figure . Here we see that the estimated potentials of mean force are relatively converged. Discrepancies are more prevalent along −18 < z < 18 where the membrane exists. We suspect that the increased viscosity leads to additional self-crossing events during the reverse phase thus hindering the overall sampling. The total amount of collected statistics per milestone for each system is shown in Figure S2. With milestoning, we can easily adaptively sample regions of poor convergence.

Figure 5

PMF profiles of each system: (a) flat, (b) small-barrier, (c) urea, and (d) codeine as estimated by milestoning, dotted red. The actual PMF is shown in black.

By comparing and utilizing the two equations derived, eqs and 30, the difference in the derivations must be considered for proper analysis and application. Equation was derived using the Smoluchowski equation, which assumes that the underlying system dynamics is overdamped. This assumption is also used in the widely employed ISD, eq . Our relation of the MFPT to the permeability does not require the estimation of a position dependent diffusion tensor. While there are many methods for calculation of the local diffusivity, these methods suffer largely from convergence issues.[34,40] The theory of milestoning does not compute diffusion coefficient directly, but rather estimates rates of passage, probability, and flux. Using eq , the MFPT from milestoning can be used to compute permeability coefficients along with the stationary probability. To the best of our knowledge, this formulation of the ISD using the MFPT has not been shown previously.

As aforementioned, the assumption that the dynamics inside of the bilayer are overdamped is largely untested. Several instances of anomalous diffusion have been observed in crowded environments among many others.[23,41−44] Similarly, the dynamics of water near hydrophobic surfaces is often complicated.[45] Thus, we present a formula relating the permeability to the crossing probabilities, ρ, which makes no assumptions about the dampedness. While the Smoluchowski equation is used in the derivation of eq , it was only used in the regions containing bulk solvent. Remaining transport properties are encapsulated by the probability ρ, which, when calculated from milestoning, makes no assumptions about the dynamics.

Our scheme is similar to that proposed by Cardenas et al.[28,29] excepting that we obviate the need to compute the stationary flux across the first milestone by Langevin dynamics simulations. In eq , the flux across the first milestone is captured by our system of two Smoluchowski equations. Therefore, we eliminate a potentially cumbersome, albeit inexpensive, step in the estimation of the permeability coefficient.

We also demonstrated the utility and accuracy of both of the derived equations using a 1D Langevin dynamics model with various PMF and diffusion profiles. In the case of the flat profile, the results of the simulations, both brute-force and milestoning, can be directly compared with theoretically derived results. For all cases, the theoretical ISD and MFPT-ISD were the same. This is not suprising, given that both are derived from equivalent expressions. Comparing the brute forced MFPT-ISD to the theoretical for the flat profile, we see that they are also comparable within error. In contrast, the flat milestoning PBCP yields a permeability coefficient that was smaller than predicted by the ISD. This discrepancy could be an example of how dampedness plays a role in the difference between the true and predicted permeability coefficients. Because the molecule simulated was not perfectly overdamped, there may be artifacts from conservation of momentum in Langevin dynamics, the formula using ρ may have captured the influence of that effect. Considering that an overdamped integrator would likely have yielded permeabilities that were more consistent between the two equations, one may question why we used the more general Langevin integrator in our test simulations. While an overdamped simulation integrator would have allowed us to apply eq more accurately, we wanted to closely reproduce the situation that would be encountered in practical application of these equations, for instance, during an all-atom MD simulation with milestones across a bilayer, where an assumption of overdamped diffusion would not necessarily hold. For this reason, we also wanted to observe how different the predicted permeabilities would be for the same trajectories using the two different formulas. We anticipate that additional investigations of various permeabilities using all-atom molecular dynamics with milestoning may yield helpful insight about the structural characteristics of the bilayer by comparing the permeabilities estimated using eqs and 30.

The two druglike molecules that we simulated, urea and codeine, provide insight about the potential usefulness of using milestoning to approach practical permeability problems. We obtain quantitatively accurate results for both urea and codeine. Unfortunately, no crossing events for urea were observed by brute force despite substantial efforts. The effect of increased sampling on transition probabilities is explored using bootstrapping in Figure S3. Because of the presence of large transition barriers, the brute force calculation of permeability is intractable, excepting for potentially the use of hyperspecialized hardware such as Anton.[46] Milestoning theory provides a framework for the estimation of membrane permeability using the MFPT in combination with transition probabilities via either the MFPT-ISD or the PBCP.

Conclusions

We have derived two useful new relations for the calculation of permeability: first, a relation of the ISD to the MFPT assuming overdamped dynamics; second, we derive the PBCP, which is agnostic to the underlying diffusional dynamics. Using simulations of our toy systems, we have demonstrated the utility and validity of each method. By comparing results obtained by the two methods, one can gain insight into the deviations of the permeability coefficient when the overdamped assumption is not valid. We anticipate that both approaches will be useful to the computational biophysical community when estimating permeability coefficients across membranes, particularly in the case when milestoning is being used in combination with molecular simulation to determine permeabilities.