Membrane permeability of small molecules from unbiased molecular dynamics simulations

Permeation of substances through membranes is a fundamental cellular process that determines many pivotal properties of cells. Crucially, the permeability coefficient P informs drug efficacy since it describes the rate at which substances translocate into and out of cells and thus their ability to attain molecular targets. While specific molecular species such as ions undergo facilitated permeation through channel proteins, permeation of most molecules occurs in a non-facilitated fashion, i.e., by passive diffusion through lipid bilayers. The corresponding permeability coefficients P are, in practice, often represented by de facto estimates. One example is the water–octanol partition coefficient in Lipinski’s rule of five,¹ which characterizes drug-like molecules. Such solubility-based estimates roughly capture the energetic contributions to the translocation rate. However, they do not include the kinetics, which are much harder to obtain from experiment.² 

In contrast, molecular dynamics (MD) simulations provide a detailed view of the kinetics and are, in principle, suited to quantify P precisely, given an accurate force field.³ Various methods have been used for calculating permeabilities from simulations, most of which rely on enhanced sampling methods and the inhomogeneous solubility diffusion (ISD) model.^2,4,5 This one-dimensional model relates P to the local diffusivity and free energy gradients along a transition coordinate. However, it is not clear a priori to which degree and in which situations the assumptions of the ISD model are justified. Moreover, even for unbiased simulations, calculating P involves technical challenges regarding the system and simulation setup, the choice of the membrane interface, length and time scales for fitting the ISD model, and error estimation. Those difficulties are underappreciated and differ widely for different classes of permeants.

The present work provides guidelines, tools, and theoretical underpinnings for calculating permeability coefficients from unbiased simulations. To cover the most important scenarios, three small molecules are considered: oxygen (hydrophobic), water (hydrophilic), and ethanol (amphiphilic). Despite their fundamental differences in solubility, unbiased simulations on the sub-μs time scale provide sufficient sampling to estimate their permeabilities through phospholipid bilayer membranes. The existence of unbiased estimates for these molecules makes them ideal candidates for scrutinizing the capabilities and limitations of simulation models.

Generally, existing methods for calculating permeability from simulations fall in one of two categories based on whether P is derived from global or local observations; see Fig. 1. The global observables are closely related to the definition of P in a setup with two water compartments that are separated by a membrane. In this setup, P is the proportionality constant between the concentration difference Δc of permeants in the water compartments and the net flux J of the permeant through the membrane.³ In a nonequilibrium setup, these observables can be directly measured to evaluate P via “flux-based counting,”

By invoking Fick’s first law, Eq. (1) can be reformulated in a way that applies also for systems in equilibrium. Here, the global observables are the crossing rate r, i.e., the number of permeant transitions through the membrane per unit time and area irrespective of direction, and the equilibrium concentration $cw$ of the permeant in water,³ and

(see also Sec. S1 of the supplementary material). Albeit sometimes formulated differently, this “transition-based counting” method and its variants have been used in various previous studies.^6–15 Both counting methods are almost model-free as they make no assumptions about the nature of permeant kinetics inside the membrane.

The required global observables J, Δc, and r, however, are generally not available from biased simulations, and full transitions are often too rare to characterize P precisely from unbiased simulations. In order to utilize enhanced sampling techniques, one must switch to a local framework, where the permeability can be derived from the individual rates between microstates as in milestoning^16–18 and Markov state models.¹⁹ The two are closely related as milestoning can be understood as Markov state modeling in trajectory space.²⁰ However, milestoning is also applicable to non-Markovian kinetics.²¹ Both approaches can utilize many short simulations for finding the transition rates between microstates. Moreover, they allow for a balanced sampling of low- and high-energy segments of the transition. The same also applies to transition interface sampling methods.^22,23 Note that the transition-based counting method can be considered as a special case of milestoning with milestones only at the two membrane–water interfaces.²⁰ 

Finally, the most-used class of methods builds on diffusive modeling, predominantly the ISD model. The ISD model characterizes the permeability via the position-dependent diffusion D(z) and free energy F(z) along the transition coordinate z, and the permeability directly follows from solving the Smoluchowski equation under stationary conditions,

Here, F^ref is the reference free energy of the permeant in water, h is the width of the bilayer centered around z = 0, and β = 1/(k_BT), with the Boltzmann constant k_B and temperature T. While the permeants’ center of mass along the bilayer normal z is the natural transition coordinate for small molecules, the path can be more complicated for larger permeants,²⁴ and the ISD model has also been used in this situation.²⁵ 

Applying the ISD model involves practical and theoretical challenges. From a practical point of view, identifying the transition path and the free energy along this path can be difficult. To this end, others have applied enhanced sampling methods such as (replica exchange) umbrella sampling,^4,5,26–28 (transition-tempered) metadynamics,^25,29 and adaptive biasing force methods.¹¹ Identifying suitable reaction coordinates becomes especially difficult when deformations and undulations of the membrane have to be taken into account.³⁰ Moreover, some permeants can undergo dynamic changes in the protonation state as the permeant translocates into the hydrophobic core of the bilayer.²⁸ Once the transition path is identified, the diffusion profile can be computed via autocorrelation functions of force,³¹ position,^5,32,33 or velocity^5,32 that are evaluated at multiple locations along the transition coordinate. Alternatively, a global diffusive model can be fit to optimally explain permeant displacements along z.^11,32,34

From a theoretical perspective, a one-dimensional diffusive model presents an idealized description of the true kinetics. The ISD model assumes that the diffusion constant along the transition path can be defined as a function of the transition coordinate, which is not obvious, e.g., if multiple paths exist. Moreover, memory effects of the medium can corrupt the Markovian assumption that underlies diffusive modeling. Therefore, the global model fitting approach has also been extended to a fractional diffusion equation.³⁴ Taken together, the assumptions underlying Eq. (3) may introduce substantial errors to the permeability estimate.

For the simulations in the present work, estimates from the ISD model are compared to the nearly model-free permeabilities obtained from transition-based counting. Section II recapitulates the counting method, establishes its validity in the context of the ISD model, presents a maximum likelihood approach for globally fitting diffusion and free energy profiles, and introduces two new software packages to facilitate those calculations. Section III describes in detail the calculation of permeabilities for oxygen, water, and ethanol through phospholipid bilayers, studies the influence of various method parameters, and compares results between the two approaches. Section IV uses the findings from these examples to develop best practices and discuss implications for biased approaches. Section V provides conclusions. The second part of this study³⁵ (henceforth referred to as Paper II³⁵) applies the methods described here as well as flux-based counting and a compartmental model to study in detail ethanol permeation through various bilayers.

This section first describes the methods used to evaluate the permeabilities from counting (Sec. II A) and the ISD model (Sec. II B). It then points out the relationship between the two approaches (Sec. II C) and presents the analysis software (Sec. II D) and simulation setups (Sec. II E).

Let us briefly revisit the transition-based counting method as derived in our recent review.³ The picture underlying Eq. (2) is that of a quasisteady state, where the dimensions of the water compartments and membrane as well as the water-membrane partitioning of permeants are constant up to statistical fluctuations. If the system contains multiple water compartments, the distribution of permeants between these compartments does not necessarily need to be in equilibrium, but the molecular systems in the present paper contain only a single bilayer and a single water phase, and the rest of this section assumes such a single-bilayer setup. Double bilayers are explicitly considered in Paper II.³⁵ 

Given an MD trajectory of length T_sim, let z denote the bilayer normal with the membrane’s center of mass at z = 0, H be the average length of the unit cell along z, and A_xy be the average cross-sectional area of the membrane. Furthermore, let h be the thickness of the bilayer, i.e., the distance between the water–membrane dividing surfaces. The concrete choice of h is not trivial, and its influence on the permeability coefficient will be discussed later (see Sec. III C 5). To compute P from transition-based counting, two quantities need to be computed: the equilibrium concentration $cw$ and the permeation rate r.

The steady state concentration $cw$ is the average number of permeants in water, divided by the total volume of the water compartment. Given the permeant distribution p (normalized to ∫p dz = 1), this concentration is evaluated as

where N is the total number of permeants in the system and $hw/2$ is the dividing surface. When using this equation, it is important that p(z) is roughly constant at and beyond $hw/2$⁠. Another convenient way to evaluate $cw$ from an equilibrium simulation utilizes the potential of mean force,¹⁵ 

where the reference free energy F^ref is the average value of the potential of mean force in the water phase. In practice, F^ref is computed over a range of bins far away from the membrane, where F(z) is roughly constant.

The rate r is the number of permeation events per unit time and membrane area, n_events/(A_xy · T_sim). In the simplest case, permeation events are defined as full crossings, i.e., transitions from z = −h/2 through the bilayer to z = h/2 or reverse, and the permeability is obtained from Eq. (2). In symmetric setups, statistics can be improved by introducing an additional milestone in the center z = 0. If the kinetics are memory-less (Markovian), a particle in the center has equal probability to exit on the either side. In other words, an escape from the membrane center to either water/membrane dividing surface has a 50% chance to be a crossing and a 50% chance to be a rebound. Moreover, each crossing or rebound is composed of two semipermeation events (one entry and one escape). Therefore, previous work has calculated membrane and channel permeabilities from

where Φ is an integer specified in Table I. Note that these relations hold only in equilibrium for symmetric systems, where transitions between the water phase and the membrane center are mostly Markovian, while counting full crossings is also valid in non-equilibrium and asymmetric setups.

As mentioned earlier, these counting methods can be considered as special cases of milestoning.¹⁶ In milestoning, milestones z = m1, m2, … are placed along the transition coordinate, and the trajectory of a particle is divided into segments based on which a milestone was last visited. Whenever a particle hits a milestone, its segment i is updated. This type of modeling allows the calculation of transition rates between segments and the closely related mean first passage times τ(m_i → m_j) between milestones, which, in turn, provide the permeability; see the  Appendix.

This discussion of the transition-based counting method closes with some remarks regarding its implementation. First, the entire analysis is performed in normalized coordinates z/H with respect to the instantaneous box size in order to handle constant-pressure simulations and constant-volume simulations in the same way. Second, for each permeant in the membrane, the direction from which it entered the membrane must be recorded in order to count crossings. For particles that are in the membrane at the start of the simulation, this information is not known. Therefore, the first escape would not be counted, which can lead to underestimates of the permeability.¹⁵ A heuristic solution is to initialize all permeants in z ∈ (−h/2, 0] as having entered from z = −h/2 and all in z ∈ (0, h/2] as having entered from z = h/2. Third, particle positions are only known at discrete points in time. If the output interval, i.e., the time between stored trajectory frames, is too long, transitions might be missed. In order to allow an output interval as large as possible, jumps over certain compartments are handled in a stepwise approach. For example, the simulation domain is divided into four intervals, water to the left and right of the membrane, z ∈ (−H/2, −h/2] and (h/2, H/2], and a left and a right membrane leaflet, (h/2, 0] and (0, h/2]. Whenever a particle jumps over one of these intervals between two trajectory frames, it could have gone either through the periodic boundary or through the membrane center. In this case, we will always assume that it has gone through the periodic boundary since diffusion is usually faster in the water phase than the membrane. The effects of initialization, output interval, and type of the event on the permeability estimates are discussed in Sec. III.

The ISD model, Eq. (3), requires the position-dependent diffusion and free energy profiles in the bilayer. To extract these profiles from a simulation trajectory, the present work uses a modification of Hummer’s Bayesian analysis.³² 

As in Hummer’s method,³² the transition axis z is divided into bins of width Δz,

Permeant transitions from B_j to B_i over a lag time τ_L are counted and stored as a transition matrix $(Nij)ij$⁠.

To find diffusion and free-energy profiles that optimally explain the simulated transition data, Hummer sampled trial profiles using a Monte Carlo scheme and the log-likelihood,

The present work replaces this Monte Carlo sampling with a maximum likelihood estimation. The matrix R stores transition rates between adjacent bins and is generated from the discrete profiles D_i±0.5 and F_i (the diffusion is defined at bin edges, and the free energy is defined at bin centers). The elements of the rate matrix are

This rate matrix is tridiagonal except for the corner elements, which are adapted to impose periodic boundary conditions. Equation (9) is derived from a finite-volume discretization of the Smoluchowski equation; see Sec. S2 of the supplementary material. To reduce the number of degrees of freedom, the profiles are modeled through basis functions ψ_k. Furthermore, to perform the optimization without constraints, the diffusion is enforced to be positive by writing the logarithm of D as a superposition of basis functions. This gives

with n_f (n_d) being the number of basis functions for the F profile (ln D profile) and f_k (d_k) being the coefficients that are optimized by the maximum likelihood routine. Two different sets of basis functions are used, cubic splines and trigonometric functions as in the work of Ghysels et al.³⁶ The optimal profiles are found by numerically solving

where $R$ is the space of real numbers. The minimum is not necessarily the same as the average from Monte Carlo sampling since the distribution of sampled f and d may be skewed around the maximum of the log-likelihood. However, both converge to the same global optimum in the limit of an infinitely long, diffusive trajectory, where the Monte Carlo method accepts new samples only in a very close neighborhood of the global optimum.

The optimization problem in Eq. (10) is solved using a robust, derivative-free global optimization method, the Covariance Matrix Adaptation-Evolution Strategy (CMA-ES).³⁷ This method is more directed toward the optimum than random Monte Carlo sampling, which speeds up the search for an optimal set of profiles and thus allows utilization of a higher-dimensional space of basis functions. Another consideration in switching to a maximum-likelihood estimation is that the error estimates from the Bayesian analysis tend to be overly narrow as they only represent the error within the model, which includes a fixed set of basis functions and a single set of transitions.³⁸ Instead, the present work reports error estimates as 95% confidence intervals over multiple independent replicas and multiple sets of basis functions, which turn out considerably larger.

Diffusion profiles D(z) generated from fitting an ISD model to a transition matrix N_ij usually vary between different lag times τ_L^11,36 mostly because the short-time dynamics contain non-Markovian contributions that are not represented in a diffusive model.³⁴ Moreover, the finite-volume discretization of the Smoluchowski equation may lead to non-negligible numerical diffusion; see Sec. III C 1. In the context of the ISD model, these factors usually cause an overestimation of diffusivities, which is most severe for small lag times. Problematically, shifting the focus to long lag times is not a practical solution because the transition matrix becomes increasingly insensitive to local transitions, which, in turn, makes it harder to fit a meaningful diffusion profile.

While this dilemma reflects some of the general theoretical issues with the ISD model, a remedy has been proposed by Ghysels et al.³⁶ They applied the Bayesian analysis to different lag times and then applied linear regression to fit the diffusion coefficients D(z; τ_L) to 1/τ_L → 0. The present work follows the same approach. Concretely, the extrapolation is performed for each replica of a simulation and each set of basis functions. The resulting profiles D(z; τ_L = ∞) are averaged over all replicas and sets of basis functions, and error estimates are calculated as 95% confidence intervals.

As apparent from its derivation, the transition-based counting method does not require diffusive dynamics inside the membrane. However, it is reassuring that Eq. (2) can also be obtained in the framework of the Smoluchowski equation. The proof is given in the  Appendix.

Two software packages were developed to facilitate the permeability calculations. Rickflow³⁹ is a Python package that implements well-validated OpenMM simulation protocols using the CHARMM force field. It also provides analysis classes to compute observables including potentials of mean force, the number of entries, escapes, and crossings of a membrane, characteristic permeation times, and transition matrices. Many common trajectory formats are supported by the way of the MDTraj package,⁴⁰ and analysis routines are vectorized using NumPy.

The maximum-likelihood estimation for the ISD model is implemented in the package diffusioncma.⁴¹ Its core functionality is implemented in C++ and uses the libraries Eigen and libcmaes for linear algebra and optimization routines. A wrapper module transfers the optimization routine into the Python ecosystem and provides many useful utility functions for postprocessing and analyzing the profiles.

Both packages are open source. Examples are given in Sec. S4 of the supplementary material.

Water and oxygen trajectories generated with the CHARMM program^43,44 were taken from previous work.^3,36,42 Trajectory frames were stored at 1 ps intervals.

Ethanol simulations were run using two different approaches, mostly to study the effect of the integration algorithm. First, simulations were run using the setup described in Paper II.³⁵ Systems containing water and lipids were run for 30 ns of NPT equilibration in NAMD.⁴⁵ Ethanol molecules were then added to the water phase using Packmol.⁴⁶ Production simulations were run using a Langevin thermostat at 298.15 K in the isothermal–isobaric ensemble at atmospheric pressure (henceforth denoted as NAMD/LD). The time step was 2 fs. In order to reduce the effect of added friction, the collision frequency was relatively small, 1 ps⁻¹.²⁷ Trajectory frames were written out every 10 ps.

Second, simulations of the same systems were set up in OpenMM version 7.3.1.⁴⁷ using the Rickflow package. Dynamics were run in the isothermal–isobaric ensemble using a Nosé–Hoover chain^48–50 velocity Verlet integrator as implemented in openmmtools⁵¹ and a Monte Carlo barostat⁵² for membranes. The integration time step was 1 fs, and trajectory frames were written out with an output interval of 1 ps. All covalent hydrogen bonds were constrained using the SETTLE and CCMA algorithms.^53,54

The CHARMM36 lipid force field was employed for all simulations.⁵⁵ As advised for membrane simulations, Lennard-Jones interactions were cut off at 12 Å with a force-switching function between 8 Å and 12 Å.⁵⁶ Particle mesh Ewald was used for long range electrostatics.

Table II summarizes the primary simulations presented in this work. Figure 2 shows snapshots of the POPC simulations. Methods for the simulations of O₂ in POPC bilayers and homogeneous solutions (hexadecane and water) discussed in Sec. III C 4 (influence of the thermostat) are contained in Secs. S5 and S6 of the supplementary material, respectively.

Results are presented starting with permeant distributions (Sec. III A), then moving on to permeabilities from counting (Sec. III B), and the ISD model (Sec. III C).

Table III and Fig. 3 show system properties and potentials of mean force (PMF) for the five systems. PMFs were computed directly from the probability distributions from all replicas. The highly symmetric shapes of the PMFs indicate good sampling and convergence. The water and ethanol PMFs are in good agreement with those calculated by enhanced sampling methods using the same force field,^57,58 although some discrepancies are expected due to a different ethanol concentration. We are not aware of equivalent studies of oxygen.

As expected from their solubilities, the distributions for water, oxygen, and ethanol differ greatly. Water has a high barrier at the bilayer center, ∼10 k_BT for DPPC and ∼11 k_BT for POPC. In contrast, oxygen partitions favorably into the bilayer. It has a 1 k_BT barrier around z = 20 Å, which is the average position of the phosphorus atoms, and a minimum in the inter-leaflet space, which is ∼3 k_BT below F^ref. Ethanol has a similar 1 k_BT barrier near the phosphate groups, a minimum below the head groups, and a 5 k_BT barrier in the center.

At the lower ethanol concentration, the space below the head group region was equally accessible as the water phase. At the higher concentration, however, the free energy minimum at z ≈ ±12 Å became metastable, and the peak in the center narrower. This suggests that the space below the head groups became increasingly saturated by ethanols, which in effect inserted deeper into the membrane. The generally less pronounced minima and maxima of the profile at higher ethanol concentration also indicate that the structural properties of the membrane were modulated: the bilayer became more flexible. This is reflected by a higher standard deviation of phosphorus atoms along the z coordinate (2.8 instead of 2.5 Å). In comparison, the standard deviation in the water (POPC) system was only 2.1 Å. Table III also shows that a higher ethanol concentration caused a thinning and widening of the bilayer. Concretely, the cross-sectional area increased by ∼18%, and the distance between the phosphate planes decreased by ∼10%, which reflects the membrane-modulating properties of alcohols. A detailed structural analysis of bilayers in response to ethanol is found in Paper II.³⁵ 

The average phosphate distances from the center z_P were also used to define dividing surfaces at z = ±h/2 for computing the permeability 5 Å outside the phosphate planes, h/2 = z_P + 5 Å. This choice allows a consistent definition of dividing surfaces for all systems. It mostly includes the variations in the free energy profiles and leaves enough space in the water compartment to reliably record transitions into and out of the membrane.

The permeant distributions underlying Fig. 3 were then used to compute the steady-state concentrations $cw$ that are required for the counting methods. To this end, it is important to realize that c_w depends on the choice of the dividing surface h_w in Eq. (4). The same applies for Eq. (5), which depends on the region that is used to fit F^ref. Figure 4 shows c_w as calculated from Eq. (4) as a function of the dividing surface. While h_w = h is most natural in theory, the ethanol concentration at (h_w/2)/H = 0.29 was already significantly different than for points further away from the membrane. In fact, choosing h_w = h would artificially increase the permeability estimate by 10% according to Eq. (2). To this end, it is more important to enforce a constant concentration in water than to use the same dividing surfaces for calculating c_w and evaluating permeation events. This can also be seen from Eq. (A3) in the derivation of the counting method from the Smoluchowski model, where it is critical that F(z) be constant throughout the water phase. In contrast, a minor difference between h and h_w only affects the permeability estimate very slightly as will be discussed in Sec. III C 5. We overlooked this nuance in previous work³ and reported c_w = 0.0109/nm³ for oxygen, which resulted from a different choice of h_w. In the present work, all c_w were calculated over the outer 10% of the simulation box, i.e., (h_w/2)/H = 0.45, where the variations in F(z) are negligible.

With T_sim from Table II and A_xy and c_w from Table III, the final variable required for the counting methods in Eq. (6) is the number of event counts. These were computed from the trajectories using the dividing surfaces h/2 from Table III. The numbers of entries (water → center), escapes (center → water), rebounds (water → center → water on the same side), and crossings (water → center → water on the other side) are listed in Table IV. Escapes, rebounds, and crossings were counted using the heuristic memory initialization described in Sec. II A 2 so that all escapes were counted as either rebounds or crossings, i.e., n_escapes = n_rebounds + n_crossings. The number of entries was not affected by this way of initializing the counting method.

In all systems, the entry and escape counts were approximately equal, as expected from Eq. (6). For the water systems, they were even identical, which means that all waters that diffused into the bilayer center also escaped thereafter. According to Eq. (6) and Table I, the number of rebounds and crossings should also be equal (on average) for a diffusive trajectory. This was the case for water and oxygen. However, the ethanol simulations showed a large discrepancy.

The ethanol simulation at lower concentration, “Ethanol (170),” had more than twice as many rebounds as crossings so that permeabilities calculated from Eq. (6) are dramatically different for different types of events. This inconsistency was even more pronounced at the higher ethanol concentration, where the number of rebounds exceeded the number of crossings by a factor of 4. In other words, an ethanol reaching the membrane center was much more likely to return than to move on and permeate to the other side.

To understand this breakdown of the Markovian assumption for ethanol, consider the simulation frames depicted in Fig. 5. Ethanols populate the space below the head groups, where they form hydrogen bonds with lipids and each other. Ethanols that translocate into the center z = 0 still maintain these connections with ethanol molecules in the same leaflet. Permeating to the opposite leaflet is energetically unfavorable since it requires leaving the hydrogen bond network. This effect is especially visible in the snapshot of the high concentration (Fig. 5, right), where an ethanol molecule from the upper leaflet has reached the core of the bilayer, while remaining weakly bound to other ethanol molecules on the same side. At high ethanol concentrations, this effect is stronger, which can be explained from the structural properties reported in Sec. III B. First, the increased flexibility and undulations of the membrane facilitate reaching the center. Second, the increased area of the membrane creates more space for ethanols and, therefore, more potential hydrogen bonding partners. Third, the ethanol saturation of the space below the head group region pushes permeants further into the membrane, which makes brief insertions into the center more likely.

In contrast, oxygen and water did not show such collective behavior. Oxygen molecules do not form hydrogen bonds so that their kinetics in the membrane are much more independent. For water, the energetic barrier is too high and too wide to maintain strong hydrogen bonds with the head groups and water on the same side. As a result, a water molecule in the center can exit the membrane in either direction without energetic penalties, which was reflected by consistent numbers of rebounds and crossings in both water simulations.

In conclusion, the water and oxygen simulations are amenable to all counting methods listed in Table I since the Markovian assumption in the center is satisfied. Ethanol permeabilities at significant concentrations, however, can only be determined accurately by counting crossings, while counting escapes or entries may lead to gross overestimates of the permeability.

Another consequence of such collective behavior is that the permeability becomes a function of the concentration in the two leaflets rather than the mere concentration difference. Since this effect cannot be observed in single-bilayer simulations, this line of research is pursued further in double-bilayer simulations in Paper II.³⁵ In the remainder of the present work, the presented ethanol permeabilities relate to a setup with similar concentrations in the two leaflets. Paper II³⁵ shows how such estimates for symmetric concentrations can be used to estimate permeabilities for asymmetric ethanol distributions via compartmental models.

Before calculating permeabilities from the counts in Table IV, the present section discusses the estimation of uncertainties. We have previously reported uncertainties of the crossing counts using two different approaches: confidence intervals of a Poisson distribution³ and standard errors over multiple replicas.^14,15 The present section investigates whether the counts follow a Poisson distribution and whether escapes and semipermeation events are also amenable to this type of modeling.

Poisson processes describe independent rare events that are separated by exponentially distributed “interarrival times.” It is thus useful to test whether the interarrival times, in fact, follow an exponential distribution. Therefore, the times at which an event concluded (an entry or escape occurred) were recorded. The times between two consecutive events were analyzed for different event types, and the statistics were fit by exponential distributions. Subsequently, 10 000 random samples from the theoretical distributions were created, and the quantiles between the theoretical and observed interarrival times were compared in a logarithmic Q–Q plot, Fig. 6.

For water in DPPC, the quantiles for entries and escapes fell almost on a straight line, which indicates that the fitted exponential model provides a sound theoretical basis for the observations. In contrast, the statistics for all events (entries and escapes) were less well represented by an exponential distribution. This can be explained from the short residence times of water in the membrane.³ Each entry of a water molecule into the center was closely followed by the same water molecule escaping. In effect, the Q–Q plot takes a nonlinear shape, indicating that entries and escapes were not independent. Therefore, considering both entries and exits does not harm the water permeability estimate but does not help the statistics either. When counting crossings or escapes, a Poisson distribution is appropriate to model uncertainties. When counting semipermeation events, the uncertainties are larger than predicted by Poisson confidence intervals.

For oxygen, the interarrival times considering all semipermeation events were much better described by a Poisson process, as indicated by the more linear Q–Q plot. Small deviations are expected due to the relatively small number of oxygen molecules in system (10), which increases the correlation between entries and escapes. Nevertheless, the distribution of interarrival times was close enough to an exponential distribution to justify Poisson confidence intervals for all three types of counting methods (entries and escapes separately not shown).

Finally, the interarrival times for ethanol were almost perfectly exponential for both concentrations, which is visible in the Q–Q plot as well as the plot of the distributions in Fig. 7. This result is slightly surprising due to the complications for ethanol transitions described in Subsection III B 1. It shows that, despite the strong hydrogen bonding of ethanols in the membrane, crossings through the bilayer midplane were still single molecule. Ethanols did not permeate in pairs or clusters, and even brief insertions into the center z = 0 occurred independently for single molecules. This justifies the use of Poisson confidence intervals to estimate the uncertainty in the number of transitions. It also shows that transitions through the center can be non-Markovian, even when semi-permeation events occur one at a time.

Finally, the above findings can be combined to compute permeabilities and estimate statistical errors. The error estimates were dominated by statistical uncertainties in the counts, and in the case of oxygen, the equilibrium concentration c_w. Since error propagation in Eq. (6) is not easily expressed analytically, 10⁶ values for the counts were sampled from Poisson distributions around the number of independent events (water entries and exit were not considered as independent events). The same number of samples for c_w and A_xy was created from normal distributions with mean and standard deviations according to Table III. These were combined to generate samples of P for each method and permeant and compute 95% quantiles for P.

Table V lists the rates and permeabilities. The permeabilities span a wide range from 1.5 × 10⁻³ cm/s for water in POPC to 20 cm/s for oxygen, as expected from the vastly different solubilities of the considered permeants. As evident from the water and ethanol simulations, P also depends on the lipid type, temperature, and permeant concentration.

The effect of ethanol concentration was especially visible for the OpenMM simulations, while the simulations from NAMD using the Langevin thermostat lacked precision. A detailed study of such concentration-dependent behavior is given in Paper II.³⁵ The rates and permeabilities for oxygen and water were slightly but not significantly different than previously reported,³ which is due to a more accurate calculation of c_w. In addition, the error estimates were more conservative due to the propagation of the uncertainty on c_w into the permeability, which affected the precision especially for oxygen.

The counting of escapes increased precision for water and oxygen. The uncertainty for oxygen permeability was further reduced by counting all semipermeation events (entries and escapes; this reduction is partially artificial due to some dependence between entries and escapes). In contrast, these more precise counting methods led to gross permeability overestimates for ethanol due to the breakdown of the Markovian assumption in the center (see Sec. III B 1). Comparing the two different simulation protocols for ethanol, the crossing rates obtained from Langevin and Nosé–Hoover simulations were consistent. In contrast, the rates of semipermeation events and escapes diverged for the higher concentration. The most important differences between those simulations were the thermostat (Nosé–Hoover vs Langevin), integration time step (1 vs 2 fs), and trajectory saving frequency (every picosecond vs every 10 ps). As shown by an application of the counting methods to trajectories with different output intervals in Sec. III B 4, this discrepancy was caused by the output interval rather than the integration method, at least for a Langevin collision frequency of 1 ps⁻¹.

The frequency of writing out the coordinates in the MD simulation plays an important role for the counting methods. Table VI shows the event counts for different output intervals. Even for an interval of 5 ps, the first discrepancies started to appear with respect to the finest frequency of 1 frame/ps: entries and escapes were missed for all permeants. The situation became more severe for larger intervals. For water, frames written only every 25 ps led to dramatic overcounting of all types of events because water molecules translocated from one leaflet through the periodic boundary into the other leaflet between subsequent frames. These crossings through the boundary could not be distinguished from crossings through the bilayer since permeation through the water phase is fast, and the dividing surfaces h/2 were placed relatively close to the periodic boundary. The situation can be remedied without corrupting the permeability estimate for water by moving the dividing surfaces further into the bilayer (see also Sec. III C 5).

The counting was more stable for oxygen, where even the largest interval of 250 ps did not affect the numbers of events severely due to long residence times in the membrane.⁶⁰ Ethanol crossings through the membrane could also be easily distinguished from crossings through the periodic boundary, mostly due to a larger box size in the z direction. Furthermore, ethanol diffusion in water is slower than for the other two permeants (see Sec. III C 2), which makes leaflet-to-leaflet transitions through the periodic boundary more unlikely. However, dramatic undercounts of rebounds and entries were observed since many insertions of ethanols into the center were too brief to be captured at larger intervals. This also explains the above-described discrepancy between the OpenMM and NAMD counts at the higher ethanol concentration. The NAMD/LD simulations were only saved every 10 ps. The rates from OpenMM simulations at the same output interval were comparable with the ones in Table V: 99 (94–104) nm⁻² μs⁻¹ for semipermeation events and 50 (46–53) nm⁻² μs⁻¹ for escapes.

This section describes the results from the maximum likelihood estimation of the ISD model. The maximum likelihood estimation was performed for various choices of the lag time, number of bins, and basis function space. The application of the CMA-ES optimization method improved performance by 1–2 orders of magnitude compared to the Bayesian analysis based on Monte Carlo sampling as fewer evaluations of the log-likelihood are required for convergence. The bottleneck in both schemes is the matrix exponential in Eq. (8), which is computed via a Padé approximation of complexity $O(nbins3)$⁠. Consequently, discretizations ≫100 bins quickly become computationally infeasible in both approaches. However, high orders of ansatz functions are only practical in the CMA-ES optimization due to slow convergence of the Monte Carlo scheme in high-dimensional parameter spaces.

To analyze the effect of lag time and discretization, Fig. 8 shows different oxygen diffusion profiles obtained over the data from all replicas. The lag-time dependence is depicted in the left plot. At small lag times, the diffusion was overestimated, but the profiles converged smoothly over the whole domain with increasing lag time. For long lag times (400 ps), unreasonable oscillations appeared in the water phase, which suggests that the fit becomes ill-conditioned due to insufficient local information in the transition matrix.

The spatial discretization error was only relevant at small lag times. For τ = 4 ps, the profiles were strongly affected at coarse spatial resolution with 40 bins or less; see Fig. 8 (middle). These errors are not surprising since central-difference schemes as the one used for discretizing the Smoluchowski equation introduce numerical diffusion (see Sec. S2 of the supplementary material), especially when they are applied to discontinuous initial conditions as in the present case. In contrast, the diffusion profiles at τ = ∞ were consistent for different spatial discretizations (see Fig. 8, right).

The rational behind such a fit is that the oxygen kinetics are mostly diffusive beyond 10 ps, while the displacements still contain contributions from the non-diffusive regime. The same problem is encountered in the computation of self-diffusion coefficients from mean squared displacements (MSDs), where the linear fit of MSD vs lag time does, generally, not pass through the origin due to inertial contribution at small lag times. Consequently, the fit to the infinite lag time is strongly encouraged in order to remove non-diffusive short-time behavior and discretization errors from the profiles.

Finally, fits were performed for different basis function spaces to confirm that the profiles were not biased by the choice of ansatz functions. Six different settings were applied:

• (symmetric) cosine series with nine and six coefficients for F and D, respectively,

• (symmetric) cosine series with 11 and eight coefficients for F and D, respectively,

• (asymmetric) cosine/sine series with 18 and 13 coefficients for F and D, respectively,

• (asymmetric) cosine/sine series with 20 and 15 coefficients for F and D, respectively,

• asymmetric B-splines with 18 and 13 coefficients for F and D, respectively, and

• asymmetric B-splines with 20 and 15 coefficients for F and D, respectively.

While previous work³⁶ used symmetric ansatz functions to better exploit the available transition data, the present paper intentionally allows asymmetry, which can be a useful indicator for poor convergence.

The fits were performed individually for each replica and lag times τ = 10 ps, 20 ps, 30 ps, 40 ps, and 50 ps. The diffusion profiles for each set of fits were extrapolated to τ = ∞. The final profiles and uncertainties were obtained as averages and 95% confidence intervals of D(z, τ = ∞) over all replica and basis function spaces.

Figure 9 shows the diffusion profiles fit with this protocol for all permeants. The uncertainties for oxygen and ethanol were low in the membrane, indicating good agreement of the individual fits. In contrast, the uncertainties for water in the POPC membrane spanned one order of magnitude due to poor sampling in the core of the bilayer. The uncertainties for water in DPPC (see Fig. S1 of the supplementary material) were similarly large. Although the average profile agreed well with previous work by Sajadi and Rowley,⁵⁷ an uncertainty this high indicates that the ISD model is not well suited to calculate water permeability from unbiased simulations.

Although the diffusion profiles were well converged for ethanol, the non-Markovian collective behavior in the bilayer (Sec. III B 1) suggests that diffusive modeling of a single ethanol molecule is not suitable, especially at the higher concentration. This is because the Smoluchowski equation models the dynamics as Markovian on all time and length scales. While this type of diffusive modeling is generally infeasible at very small times, the ISD model assumes that there exists a regime in which the dynamics are nearly diffusive. Previous work by Chipot and Comer suggests that small alcohols do not reach this regime within the relevant permeation times even at very low concentration.³⁴ The collective effects at higher concentrations strongly amplify this behavior, and it is highly questionable that diffusive modeling of the permeation is appropriate.

An important consistency check for the ISD fit is the comparison of the free energy profile F(z) coming from the fit with the directly observed potential of mean force. Figure 10 shows the profiles from the maximum likelihood estimation. The profiles for oxygen and water matched the PMFs in Fig. 3 although water uncertainties were high in the center. Agreement was significantly poorer for ethanol despite small error bars. At the lower concentration, the plateau in the center was smoothed out, and at the higher concentration, the barrier was lowered by almost 1 k_BT. These disagreements indicate limited applicability of the ISD fit for the ethanol simulations.

The diffusion and free energy profiles in Figs. 9 and 10 were used to compute permeabilities via the ISD model, Eq. (3). The results are shown in Table VII. While the water permeabilities through POPC and DPPC were consistent with the estimates from counting methods (Table V), their error estimates in the ISD model were much wider, as expected from the large uncertainties in the profiles. The oxygen profiles were significantly better converged in the membrane, which also shows in the permeability estimates. However, the estimates were much more conservative than reported in our previous work.^3,36 This is understandable from significant variation in the reference free energy, which was previously not propagated into the error estimates. In comparison, the counting of semipermeation events led to a more precise estimate than the ISD fit.

For ethanol, the values from the ISD model and counting crossings were inconsistent. The ISD model overestimated P by a factor of 2 for the low concentration and a factor of 4 for the higher concentration. This mismatch supports the above conclusion that ethanol permeation at significant concentration in the membrane is not amenable to diffusive modeling at the individual particle level.

Table VII also shows a statistically significant difference between the Langevin thermostat simulations and the Nosé–Hoover thermostat simulations at high ethanol concentration. Figure 11 shows that this discrepancy comes from inconsistent diffusivity. Although the collision frequency was low (1/ps), the artificial friction added by random thermal noise changed the kinetics, particularly in the water phase. It is well-known that a Langevin thermostat can alter the kinetics in this way. However, the difference was small in the membrane, suggesting that simulations with a 1/ps collision frequency are acceptable for permeation studies of ethanol.

Two different systems were chosen to study the thermostat effect in more detail, water and oxygen transitions in a POPC bilayer and oxygen diffusion in water and hexadecane. Details of the simulations are provided in Secs. S5 and S6 of the supplementary material. As expected, the number of transitions vary among the replicates; Table VIII lists the results. Beginning with the bilayer, the total number of transitions N_tot is similar for NVE and Nosé–Hoover (NH) for both the water and oxygen permeabilities, lending strong support to the notion that a rigorous extended system thermostat and barostat can yield reliable dynamics. N_tot for the Langevin thermostat with γ = 1/ps (LD1) is comparable to NH for water, as shown above for ethanol, but not for oxygen; that is, the 95% CI for LD1 does not overlap with the mean of the NH. N_tot is only approximately half that for NH and NVE for the Langevin thermostat with γ = 5/ps (LD5) and clearly statistically different.

Table IX shows that diffusivities of oxygen in homogeneous media are similarly perturbed by using the Langevin thermostat. Oxygen diffusion was decreased by ∼20% at γ = 1/ps and by more than 50% at γ = 5/ps. Hence, the simulations of these simple systems yield the same conclusions as those in the far more computationally demanding POPC bilayer. While this correspondence will not necessarily hold for other permeants, it is an easy and worthwhile check. In any case, permeabilities calculated for other solutes using NAMD/LD should be validated against more rigorous thermostats.

As a final comparison between counting methods and the ISD model, Fig. 12 shows the permeabilities from both methods with respect to the dividing surface h/2. For the water systems, the value of h/2 did not affect the permeabilities in either method. While the uncertainties especially from the ISD model were large, results from both methods were consistent for all dividing surfaces. The resistivity is purely dominated by the large free energy barrier in the membrane. This allows for flexibility in the choice of the dividing surface. In practice, one can pick a dividing surface so that crossings through the membrane and through the periodic boundary are easily distinguished in order to allow larger output intervals of the trajectory.

For ethanol, the situation is similar since P is also dominated by the free energy barrier in the membrane center. Figure 12 shows that the overestimates from the ISD model were similar for all dividing surfaces, which again underlines the unfitness of the ISD model to determine ethanol permeability at significant concentrations.

Oxygen, in contrast, has a free energy minimum in the center. Its resistivity is dominated by the local free energy maxima around the phosphate planes (∼±20 Å). Therefore, the permeability was very sensitive to the choice of the dividing surface. This was reflected in estimates coming from both methods, which matched for all choices of h/2.

This paper has described in detail the computation of permeability estimates for water, oxygen, and ethanol using counting methods and the ISD model. On the way, many ways to go wrong were identified, the most important of which are reiterated in the following.

The utility of counting methods is, of course, limited by the number of permeant crossings. This is optimally in the range of 100 for statistical precision, although as few 10 might be acceptable for exploratory studies. The systems presented here all contained numerous transitions. Microsecond simulations of oxygen and water yielded an average of 34 and 119 transitions, respectively (Table VIII). Analysis of the last 200 ns from simulations of ethanol at six different concentrations in three different bilayers showed between 13 and 923 transitions in 200 ns (Tables S2–S4 in Paper II³⁵).

Now, consider a solute with a permeability the same as water, but at a concentration of 11 mM (four molecules in a simulation cell containing 20 000 waters). Hence, c_w is lowered by a factor of 5000 (55M/0.011M) compared to water. Rewriting Eq. (2) as r = 2c_w · P, the rate of crossings observed in the simulation is similarly reduced by a factor of 5000 and well out of the realm of counting methods with current computer resources. Hence, both the solute concentration and the anticipated permeability must be included when assessing the applicability of counting methods. The same considerations hold for estimating the PMF of a permeant from conventional simulations as opposed to enhanced sampling. Sampling for PMF includes unsuccessful transitions so is a little better than for rates.

A separate issue involves too many permeants in the system. The number of permeants has to be carefully considered, especially for membrane-modulating permeants such as alcohols or general anesthetics, where the permeability coefficient is sensitive to concentration. While it is encouraging that bilayer crossings in all considered simulations were single molecule, typical atomistic simulation setups will often have a higher permeant concentration than corresponding experiments. In this case, simulation results have to be scrutinized in order to rule out undesirable collective behavior.

The number of lipid and water molecules has to be carefully considered as well. Large bilayers are prone to undulations, in which case the definition of a suitable transition coordinate can become highly non-trivial.³ When counting crossings in a single-bilayer setup, it is also important that the water compartment is large enough to place a dividing surface away from the membrane and still be able to distinguish crossings through the bilayer from crossings through the periodic boundary in the water phase.

Due to increasing computational power, simulation frames are often written at increasingly large intervals to save disk space. However, permeability calculations require a detailed monitoring of permeant kinetics. The recording of entries and escapes is particularly sensitive to the output interval. It is beneficial if the transition coordinate is known in advance and the analysis can be performed on-the-fly. To this end, the Rickflow package provides an OpenMM reporter for on-the-fly recording of transition matrices and event counts.

While counting entries and escapes can improve statistics over counting crossings, it can also introduce errors. First, it generally requires a fine output spacing in order to not miss brief insertions of permeants into the membrane center. Second, when the kinetics are non-Markovian (as for ethanol in the present study), the counting of semipermeation events is no longer a reliable method for calculating permeability. Counting crossings is generally more robust. In order to not underestimate the actual number of events for lipid–soluble permeants, the permeation event counter has to be initialized with memory in order to incorporate the first escapes into the permeability calculation.

The mismatch between permeabilities from simulation and experimental studies is well-documented³ for water,¹⁴ ethanol,¹¹ and oxygen.^13,36 Currently, one of the major challenges in the force field community is the development of models that perform well in various molecular environments. Current additive force fields do not accurately predict the partitioning of solutes between the water and organic phases;⁶¹ polarizable force fields will ultimately be required to obtain agreement with the experiment for many permeants.¹⁴ Additionally, permeability can often only be measured indirectly in experimental studies, and a one-to-one correspondence between simulation and experiment is often not easily established.

The Langevin thermostat is known to alter the kinetics by introducing artificial thermal noise on the atomistic level. The introduction of friction decreases the diffusivity. For ethanol, the effect on the diffusivity inside a bilayer is small but noticeable when using a 1/ps collision frequency. The same holds for oxygen, where the permeation rate is lowered by ∼15% due to a 1/ps collision frequency. For water, the number of crossings was too low to indicate a significant discrepancy. Use of larger collision frequencies is not advised; γ = 5/ps already lowered the permeation rate by 50% for both water and oxygen. This damping of the kinetics has to be individually considered for each permeant and was not present in simulations with a Nosé–Hoover thermostat. Therefore, permeation studies should ideally use thermostats that better retain the kinetics, such as Nosé–Hoover or Lowe–Andersen.⁶² 

Error estimates from the Bayesian analysis do not necessarily represent the sampling error, as shown by the detailed analysis in this work. As a consequence, previously reported error estimates were one order of magnitude too narrow. Another factor that can lead to underestimation of errors is the use of smooth, symmetric basis functions. The present work introduced software to perform the fit more efficiently, which enables the use of higher-dimensional basis function spaces. It is advised to perform the ISD fit over multiple replicas of a simulation and to use multiple basis function spaces in order to provide conservative error estimates.

The lag time dependence of diffusion profiles has been discussed in the previous literature.^11,34,36 The present work shows that diffusion profiles at lag times as long as 50 ps are affected by non-diffusive effects coming from shorter time scales, which lead to overestimates of the diffusivity and thus the permeability. The fit to infinite lag time is a practical way to remove such contributions from the profiles and restore consistency with model-free estimates.

ISD-based approaches are often used for permeability calculations from biased simulations. However, it is important to recognize that the ISD model presents an idealized description of the true kinetics. Diffusive models are, by definition, memory-less on all scales. In contrast, milestoning also captures non-Markovian dynamics and might be more suited as a general tool for simulations of permeation.

This work is a simulation study of water, oxygen, and ethanol permeation through lipid bilayers. The permeability coefficients were obtained from unbiased simulations by using counting methods and the ISD model. To this end, two software packages are introduced to facilitate such calculations. The Python package Rickflow implements simulation workflows and trajectory analysis code to count permeation events and record transitions. The Python/C++ code diffusioncma implements a maximum-likelihood estimation for fitting an ISD model. This approach is closely related to a previously published Bayesian analysis,^32,36 but speeds the fitting up by one to two orders of magnitude.

The comparison between permeability estimates from both types of approaches showed good agreement for oxygen and water, although water permeabilities obtained from the ISD model carried large uncertainties. In contrast, ethanol as a membrane-modulating permeant was not amenable to diffusive modeling due to collective effects in the membrane. Consequently, the ISD model overestimated permeabilities by factors of 2–4. This shows that the ISD model should be used with caution, especially for permeants that interact strongly with one another and the membrane. Furthermore, due to the collective behavior, the permeability becomes a function of not only the concentration difference but also the individual ethanol concentrations in the two leaflets. This line of research is further pursued in Paper II.³⁵ 

See the supplementary material for details about the consistency between Fick’s first law and transition-based counting, discretization of the Smoluchowski equation, water diffusion profiles in POPC and DPPC, code examples, computation of diffusion constants in homogeneous media, and supplementary microsecond POPC simulations.

The data that support the findings of this study are available from the corresponding author upon reasonable request.

The authors would like to thank Andrew C. Simmonett, Mahdi Ghorbani, Attila Szabo, and Sally Pias for helpful discussions and John Legato and Rubén Meana-Pañeda for technical assistance. A.K., E.W., R.M.V., B.R.B., and R.W.P. acknowledge support from the Intramural Research Program of the NIH, the National Heart, Lung, and Blood Institute, and the high performance computational capabilities at the National Institutes of Health, Bethesda, MD (NHLBI LoBoS and Biowulf clusters). Ethanol simulations in NAMD used the Extreme Science and Engineering Discovery Environment (XSEDE) allocation under Grant No. MCB-100139 on Stampede2 and Comet. E.W. and J.B.K. acknowledge funding from the NSF (Grant No. CBET-1604576).

This section derives the transition-based counting method from the Smoluchowski equation. Consider the mean first passage time (MFPT), a quantity that is related to both the permeability and the number of crossings. The MFPT considered here is the average time for a particle to travel from b to a in a domain that has a reflective boundary condition at z = q and an absorbing boundary condition at z = a, where q ≤ b = −a ≤ 0 ≤ a (Fig. 13). For a diffusive trajectory, the MFPT can be expressed as⁶³ 

It can be assumed that F^ref = 0 is the free energy outside the membrane (|z| > b) without loss of generality.

To relate the MFPT to the number of crossings in a molecular simulation, the system with a reflective boundary condition is transformed into a periodic system in three steps. For clarity, we prefer a graphical over a formal transformation; see Fig. 13. First, consider the above-mentioned setup with one reflective boundary (S1), where one particle starts off at z = b and is reset to this initial position, as it reaches the opposite dividing surface z = a, in order to repeat the experiment. In this setup, the number of crossings in a simulation of length T_sim is given as T_sim/τ_[q,a](b). Now, consider a second setup (S2), where a reflective boundary is symmetrically introduced at z = −q. In the absence of an adsorbing boundary in S2, the particle does not need to be reset to its initial condition. This change does not affect the number of crossings. Furthermore, replacing the reflective by periodic boundary conditions (S3) also does not affect the number of crossings because of symmetry. Finally, in a periodic system with N independent permeants, one expects

The next step is to relate the MFPT to the permeability P. By splitting the inner integral over the interval [q, z] in Eq. (A1) into two integrals over the intervals [q, b] and [b, z],

where the remaining inner integral was simplified using F^ref = 0 in the interval [q, b]. In the case of a bilayer with a = h/2 and q = b = −h/2, i.e., when the permeants are restrained to the membrane, τ_[b,a](b) is also related to the permeability when the dynamics are diffusive,¹⁷ 

Combining Eqs. (A3) and (A4) yields

since the free energy vanishes outside the membrane. The latter integral over the whole membrane relates to the concentration in water through Eq. (5). By inserting Eqs. (5) and (A2) into Eq. (A5), the permeability turns out as

as in Eq. (2). While Eq. (2) was derived from fairly general assumptions,³ it is now rigorously confirmed in this much more specific setting of a diffusive equilibrium trajectory. The two considered connections with Fick’s first law (supplementary material, Sec. S1) and diffusive dynamics (here) confirm the validity of the transition-based counting method in two extreme situations: the former holds in a situation where all permeants are in the source compartment, and the latter is associated with an equilibrated system.