Exact non-Markovian permeability from rare event simulations

Permeation of compounds through membranes is important in biological and engineering processes, e.g., drug delivery through lipid bilayers, anesthetics, or chemical reactor design. Simulations at the atomic scale can provide insight in the diffusive pathways and they give estimates of the membrane permeability based on counting membrane transitions or on the inhomogeneous solubility-diffusivity model described by the Smoluchowski equation. For many permeants, permeation through a membrane is too slow to gather sufficient statistics with conventional molecular dynamics simulations, i.e., permeation is a rare event. Recent attempts to improve the description of the dynamics of such rare permeation events have been based on milestoning, which allows the study of processes at timescales beyond those achievable by straightforward molecular dynamics. The approach is not relying on an overdamped description, but, still, it uses a Markovian approximation which is only valid for small permeants that are not disruptive to the membrane structure. To overcome this fundamental limitation, we show here how replica exchange transition interface sampling (RETIS) can effectively be used on this problem by deriving an effective set of equations that relate the outcome of RETIS simulations and the permeability coefficient. In addition, we introduce two new path Monte Carlo (MC) moves specifically for permeation dynamics, that are used in combination with the ordinary path generating moves, which considerably increase the efficiency. The advantage of our method is that it gives exact results, identical to brute force molecular dynamics, but orders of magnitude faster.


I. INTRODUCTION
Permeation of compounds through another medium is essential in both biological and engineering processes. In biology, the permeation of molecules, nutrients or nanoparticles through membranes is an integral part of a functioning cell, and understanding the drug delivery process can aid the design of new cancer drugs or anesthetics [1][2][3][4]. In chemical engineering, the transport of molecules can play a role in the selectivity of the molecules [5,6]. In highly complex chemical systems, it is valuable to unravel the various diffusion pathways, which can aid the optimization of a chemical reactor setup. Despite the existence of spectroscopic methods, such as fluorescence spectroscopy [7,8] or EPR experiments [9], the insight in diffusive transport at the molecular scale is difficult to obtain experimentally. Especially in inhomogeneous media where the diffusion of permeants is a function of the location, experiments usually only provide a global effective diffusion constant, without discerning local differences at the molecular scale. It is in this respect that molecular dynamics (MD) simulations can play a major role. MD creates molecular trajectories and transport properties are directly observable at the molecular scale, such as the permeability P.
The permeability of a membrane is the flux through the membrane as a response to a concentration gradient over the membrane. A first standard approach to derive the permeability P from MD simulations is the counting method, which is based on measuring the rate of membrane transitions per unit of time and area [10][11][12]. Another standard approach is Bayesian analysis (BA) using the Smoluchowski equation, which assumes a position-dependent concentration profile as well as a position-dependent diffusion profile across the membrane [12][13][14][15][16][17]. The Smoluchowski equation is a pure diffusion model, where memory effects are not modeled. The serious limitation of these approaches is that, when the permeability is low, the statistics provided by MD might be insufficient, even when the trajectories are extended up to 1 microsecond of simulation time.
A more recent approach is the extraction of P from milestoning. Milestoning is based on the sampling of many short trajectories released from equilibrium distributions at hypersurfaces (milestones) [18]. Typically, these hypersurfaces are defined as subsets of configuration space with fixed values of the reaction coordinate (RC). The milestoning method then counts from each originating hypersurface how often the left or right hypersurface is hit first and the time that it takes to let this happen. These two properties for the different milestones are then combined such that the dynamics of the system can be described by a Markovian hopping sequence from one surface to the other. Rates, mean-first-passage times and other relevant dynamic and thermodynamic data can then be obtained. The central assumption underlying milestoning is that the set of first hitting points, of the trajectories originating from one surface with a hypersurface left or right from it, is again distributed according to the equilibrium distribution. This leads to conflicting requirements. On the one hand, the milestones need to be set closely as this guarantees the highest efficiency boost compared to brute force MD. On the other hand, they need to be well separated for the central assumption to be sufficiently satisfied. In addition, for the assumption to hold, the choice of the RC is crucial and should ideally coincide with the committor [19]. The committor is generally an unknown coordinate that not only describes the position of the permeant relative to the center of mass of the membrane, but also should include nontrivial membrane deformations. Obtaining approximate forms of the committor is generally an immense task [20].
As a way to overcome the limitation of the Markovian assumption, this work will make use of the transition interface sampling (TIS) framework and of its extension, the replica exchange TIS (RETIS) formalism [21][22][23]. (RE)TIS gives a completely non-Markovian treatment of the interface hopping probabilities, while still being orders of magnitude faster than brute force MD. The overall rate is obtained in a divide-and-conquer mindset, by constructing a series of conditional probabilities to reach the next interface. Here we will show how the permeability, as defined in steady-state nonequilibrium conditions, can be transformed into functional form that depends on the equilibrium path ensemble properties which can be obtained from RETIS (instead of the typical canonical, NVT, or isobaric, NPT, phase space ensembles).
In an extension, we also show how the rates of transition can be defined in a meaningful way for an infinite system or a system with periodic boundaries such that it can be linked to the permeability. When the phase neighboring the membrane is unbounded, particles have a probability to indefinitely diffuse in the opposite direction of the membrane rather than through the membrane. Also when a system has periodic boundary conditions, it is difficult to detect whether particles reach the other side of the membrane through the periodic boundary or through the membrane itself. The permeability formula will be adapted to treat those cases. In addition, two additional moves will be introduced in the MC sampling of the interface ensembles. They will improve the efficiency when the simulation box contains multiple permeant molecules or when the membrane is symmetric with respect to the membrane center.
This article is organized as follows. In Sec. II, we review the existing approaches for obtaining P from MD simulations. In Sec. III, the (RE)TIS formalism is revised and the path ensembles are defined. In Sec. IV, the theoretical derivation of the permeability from (RE)TIS is presented. In Sec. V we derive the needed adaptation to treat systems with periodic boundary conditions. In Sec. VI, the two new MC moves are presented. In Sec. VII, we illustrate the accuracy and effectiveness of our approach by computing the permeability of a few basic model test systems. We end with concluding remarks in Sec. VIII.

A. Direct counting
The permeability is defined as the ratio of the net flux J of particles transiting a membrane in the steady-state regime when imposing a small concentration difference c over the membrane, To compute the permeability, one could consider imposing the concentration difference at the membrane boundaries, as suggested by the definition in Eq. (1). However, at the molecular scale it is troublesome to impose a nonequilibrium steady-state concentration gradient, especially when periodic boundary conditions are applied. A more common approach is therefore to count transitions in both directions through the membrane in an equilibrium simulation [10]. At first sight, Eq. (1) is no longer adequate for the computation of P. Indeed, note that J in Eq. (1) is the result of both crossings from left to right and from right to left, corresponding to a positive flux J + and a negative flux J − . These fluxes are proportional to the concentrations at the left and right hand side of the membrane, respectively, and bear opposite signs. The net flux is, hence, the result that is left after a partial cancellation of J + and J − , and in an equilibrium situation, where c = 0, the cancellation is complete and the net flux is zero. Whereas the permeability in this case is no longer measurable experimentally, in simulations it is still possible since it is relatively easy to trace the J + and J − fluxes individually from the MD trajectories and one can write Here, z = 0 is considered to be the center of the membrane of thickness h with borders at z = ±h/2, and c(z) is the concentration profile across the membrane. The reference concentration is the concentration outside the membrane c ref = c(−h/2) = c(h/2) and corresponds to the concentration of permeants in the bulk liquid.

B. Smoluchowski equation
A second common approach is to run equilibrium MD simulations and to analyze these assuming the validity of the inhomogeneous solution-diffusivity (ISD) model, where transport is modeled by position-dependent Brownian diffusion (diffusion profile D(z)) on a free energy landscape (profile F (z)), as governed by the Smoluchowski equation [15]. The free energy is related to the permeant concentration in equilibrium through the Boltzmann probability, c(z) ∼ exp(−βF ), where β = 1/(k B T ) is the inverse temperature, k B is Boltzmann constant. Given the two profiles F and D, the permeability follows directly from solving the Smoluchowski equation for its steady-state solution in the presence of a fixed concentration difference c over the membrane [13,16], where F ref is the free energy at the reference location, usually at z = − h 2 . The structure of Eq. (3) shows that the free energy profile is the dominant contribution to the permeability. Several cases can now be thought of for F (z): F (z) has a barrier, is flat, or has a well. The last two cases are rather academic, as realistic membranes often form a combination of free energy barriers and wells, e.g., for O 2 permeation through phospholipid bilayers [16]. First, when the permeation implies crossing a high free energy barrier, the integration range in Eq. (3) can readily be reduced from if the integrand can be neglected for the outer regions −h/2 < z < −h /2 and h /2 < z < h/2. Given the exponential dependence on the free energy barrier, the integration boundaries −h /2 and h /2 can often be chosen rather close to the maximum of the free energy barrier as long as F ref is taken in the bulk region. In the direct counting method, the neglect of the outer regions in the integration relates to the fact that nearly all transitions from −h/2 to h/2 contain one and only one single −h /2 to h /2 transition. Second, the permeability in a homogeneous medium where F (z) is a flat profile, equals P = D/h with D the diffusion constant. The permeability halves when doubling the thickness h. This shows that the integration boundaries and reference region need to be chosen with some care whenever the free energy barrier is relatively small. Third, when F shows a free energy well, the permeants may be trapped in the membrane, and the integration boundaries should also be chosen with care.
Besides these conceptual fundamental issues related to the definition of the permeability, there are also practical issues. The danger exists that not all regions inside and outside of the membrane are accurately sampled in the equilibrium MD simulations. Lastly, slow sampling can also originate from a low diffusivity [D(z) in Eq. (3)], e.g., for permeants that are bulky. The challenge of the Smoluchowski approach is to determine the model parameters, i.e., the F (z) and D(z) profiles, which can be extracted a posteriori from the MD trajectories in various ways [10,12,15,16,28,29,34]. When MD is inadequate because of its limited timescale, the approach can be combined with rare events simulation techniques like umbrella sampling [35], adaptive bias force [36][37][38], and biased diffusion [39,40]. All these free energy methods, with the aim to compute the permeability via Eq. (3), have however as a shared limitation that they build on the validity of the Smoluchowski equation. The validity of this equation is questionable for many complex systems in which the transfer of permeants over the barrier involves other types of motion, like a rotation of the permeant or a local stretch of a membrane opening [38,41].

C. Path sampling approaches
Path sampling methods seem to provide a natural solution to the permeation problem since they are designed to maintain the natural dynamics of the process as much as possible, while still allowing the sampling of events that happen on long timescales. Among the different path sampling methods, applications on the permeation problem have so far mostly adopted the milestoning method [18]. In this technique, phase space is divided in domains that are separated by interfaces, called milestones. Trajectories are initiated at each milestone and run until they cross another milestone. The statistics over this set of short trajectories give the mean first passage times between pairs of milestones, which are then incorporated in a Markovian rate network model to extract the overall rate. Cardenas and Elber [42] proposed the formula for the permeability where J 1 is the flux of particles hitting the first milestone per area and per time in equilibrium, q is the absolute flux vector of trajectories crossing the milestones, when solving the Markovian rate model for its steady-state solution with specific boundary conditions. Cardenas and Elber applied this to the permeation of a small peptide [42] or water molecule [43] through a 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) phospholipid bilayer, while Fathizadeh and Elber simulated potassium permeation through this DOPC membrane [44].
Recently, Votapka et al. derived an alternative formulation for the permeability based on milestoning [45]. The advantage of milestoning is that the MD trajectories that need to be generated are very short. However, the milestoning relies on the Markovian assumption that the system loses memory when it hits the next interface/milestone [18,46]. This assumption is only correct if the milestones are chosen along the isocommittor lines [19]. Hence, since the milestones are defined by fixed values of the RC, the RC should be the committor [47][48][49][50]. The committor is therefore often considered as the ideal RC since the dynamics projected on this one-dimensional coordinate becomes Markovian and models relying on the Markovian assumption, like Smoluchowski and milestoning, become exact. Note that the committor is in principle defined in phase space though by assuming that the dynamics is overdamped the committor can be defined in configuration space alone.
Hence, if milestoning is applied with this ideal RC (the committor) then trajectories released from the same milestone have the same probability of reaching h/2 before −h/2 regardless the point of origin within that milestone. This RC should account for all relevant rotations of the permeant, collective motions, and deformations of the membrane that could be vital for the permeation process. Including these motions within a single coordinate is highly nontrivial even if the committor is assumed to be only configuration space dependent based on the overdamped approximation. In practice, this is generally not even attempted. Rather, a simple intuitive RC is chosen such as the z coordinate of the permeant that is followed. This pragmatic choice will generally invalidate the Markovian assumption leading to a systematic error that can be mild or severe depending on the system. Another path sampling approach that has similarities with milestoning is the partial path transition interface sampling (PPTIS) [51]. The PPTIS method is a Markovian variant of the TIS method. Still, it uses a less stringent assumption than milestoning by including more memory in the dynamical description. In PPTIS, if the system hits an interface, the chance to move to either its left or right interface depends on the history on the path i.e., from which interface (left or right) it came from. Still, no memory is retained before that point. The systematic error due to the nonideal RC is therefore presumably lower in PPTIS than in milestoning.
Path sampling methods that include the complete history dependence of the dynamics, such as transition path sampling (TPS), transition interface sampling (TIS), forward flux sampling (FFS) [52], adaptive multilevel splitting (AMS) [53], and replica exchange TIS (RETIS), are exact and independent of the RC, which is a big advantage in complex systems. The sampling of complete trajectories will make these methods generally more computationally intensive than milestoning or PPTIS, though a quantitative comparison relies on the tradeoff between systematic and statistical error. FFS and AMS are based on a splitting approach which makes it applicable for nonequilibrium dynamics as well. On the other hand, the lack of backward-in-time integration limits the applicability of splitting methods to stochastic dynamics and leads to a relatively high risk of producing nonrepresentative transition trajectories [54]. The other three methods, TPS, TIS, and RETIS, are all based on a MC sampling procedure in path space. The original TPS approach for rate calculations is no longer being used in practice as TIS is both faster and more accurate than TPS. RETIS is even more efficient than TIS, but requires a more complex implementation. By the emergence of open source path sampling libraries like OPS [55,56] and PyRETIS [57,58] the latter aspect has become less of an issue.
In this paper, we show how RETIS can effectively be used to compute permeability coefficients equally exact as the direct counting approach, but orders of magnitude faster. In the next section, we shortly introduce the RETIS approach, then show how it can be amended for permeability calculations in Sec. IV, and show in Sec. VI two new path MC moves that can further enhance ergodic sampling.

III. REPLICA EXCHANGE TRANSITION INTERFACE SAMPLING
TIS and RETIS are rare event techniques that allow to compute rate constants k when transitions have to overcome a high free energy barrier, and thus transitions are unlikely to be observed in a standard MD simulation. Intuitively the rate constant k can be expressed as a number of transitions, from reactant state to product state, per unit time and per amount of reactants. However, the translation of this phenomenological rate constant into a computational measurable property is far from trivial since it requires a microscopic definition of the reactant state and product state. If a single dividing surface is used to assign whether a molecular system is in the reactant state or product state, e.g., based on a geometric observable being lower or higher than a specific value, it is expected to observe many correlated recrossings of the dividing surface. Only when the system has moved far enough beyond the dividing surface, it can be considered as stabilized to the other state. To avoid the issue of correlated recrossings, it can be generally better to consider two separate dividing surfaces left and right of the barrier, λ A and λ B , respectively [see Fig. 1(a)]. At the left of λ A and at the right of λ B , the system is considered committed to the reactant state and product state respectively. The disadvantage is that the barrier region between the dividing surfaces is not assigned to either state. A transition from reactant state to product state needs to cross the no man's land between A and B, which makes it difficult to assign for each transition a specific point in time at which it takes place, which is essential to avoid overcounting.
The TIS and RETIS approaches circumvent this problem by introducing overall states which are history-dependent state definitions. Using the two dividing surfaces, λ A and λ B , the system is part of stable state A if the value of the chosen RC is below λ A and it is part of stable state B if it is higher 033068-4 than λ B . The overall states are denoted by the curly letters A and B, and they include the stable states A and B, respectively. Further, any instance of the system traveling in no man's land is assigned to overall state A or overall state B based on the stable state that was most recently visited.
The advantage is that the system is always assigned to either overall state A or B, and the overall state regions are rather insensitive to the placement of the dividing surfaces λ A and λ B as long as these are reasonable [23]. In addition, a transition from A to B is well defined without hindrance of recrossings and leads to a microscopically measurable net flux without perturbing the equilibrium conditions. As such, the phenomenological rate constant can be expressed as where h X (t ) equals 1 when the system is in state X at time t and 0 when the system is not in state X at time t. The interpretation of the ensemble average that involves history dependent functions is given in Appendix A.
The above equation counts the rare crossings from left to right through the surface λ B of points that actually come from λ A (when the equations of motion are followed backward in time, λ A is reached without crossing λ B ). Since we can assume equilibrium, the same number of counts per second will be obtained by considering crossings with λ A that after crossing that interface reach λ B without recrossing λ A . This allows us to write the rate k as where f A is the conditional flux (crossings per time) through λ A counting all crossings in the positive direction per time spent in state A, and P A (λ B |λ A ) is the crossing probability, i.e., the chance that a crossing with λ A is followed by a crossing with λ B before a recrossing with λ A occurs. f A can easily be computed with MD, as is done in TIS, since crossing λ A is not a rare event. In RETIS, the flux is computed from the average path lengths of two path ensembles. The crossing probability P A (λ B |λ A ) is generally too low to be computed by straightforward MD, but can be recast into the following factorization by defining a set of nonintersecting interfaces: Here, P A (λ i+1 |λ i ) is the history dependent conditional probability that, given there is a crossing with interface λ i for the first time since last λ A = λ 0 crossing, interface λ i+1 will be crossed before λ A . These conditional probabilities are computed in n − 1 different path sampling simulations. Each simulation samples a different path ensemble using a set of different MC moves to generate new paths in the ensembles. The ensemble [i + ] consists of all possible paths that start at λ A and end at λ A or λ B and have at least one crossing with λ i . The MC approach is tuned such that the same statistical distribution of paths is generated as the distribution of paths that would result if these are cut out from a hypothetical extremely long MD simulation. The fraction of paths in the [i + ] path ensemble that cross λ i+1 equals P A (λ i+1 |λ i ). Performance of TIS is optimal when the number of interfaces and the spacing between the interfaces are tuned such that each conditional probability is around 0.2 [59]. The most important MC move is the so-called shooting move in which a random point of the previous path is picked, its velocities are randomly modified to generate a new phase point, and from this phase point the equations of motion are followed backward and forward in time to generate a new path. The new path will be accepted if it fulfills the requirements for the specific ensemble (like crossing λ i in the [i + ] ensemble) and, depending on the type of velocity randomization procedure, an additional Metropolis acceptance/rejection step will be invoked. In this work, we applied the aimless velocity modification in which velocities are regenerated, independent from the old velocities, from a Maxwell-Boltzmann distribution [60]. Compared to TIS, RETIS has one additional path ensemble called [0 − ]. Like the other ensembles, this ensemble contains paths starting at λ 0 , but from there the paths move in the negative direction away from the barrier. The paths are terminated when they cross λ 0 again. While the flux f A in Eq. (7) is computed with straightforward MD in TIS, in RETIS it is computed from the average path lengths in the [0 − ] and [0 + ] ensembles (see Sec. IV A).
Another difference between RETIS and TIS, is that RETIS employs additional MC moves. Since RETIS is purely based on path sampling simulations, instead of MD and path sampling simulations, replica exchange moves between the different path ensembles (parallel path swapping [22]) can be applied throughout the complete RETIS simulation. The swapping moves enhance the sampling in a similar way as parallel tempering [61], since the [i + ] path ensemble for a given i tends to contain trajectories moving on a higher energy surface than the trajectories of the [(i − 1) + ] ensemble. The higher energy trajectories are less likely to get trapped in specific reaction channels that are separated by free energy barriers orthogonal to the RC λ [62]. As paths between [(i − 1) + ] and [i + ] are sometimes swapped, also the [(i − 1) + ] path ensemble will sample these different reaction channels more easily. The combination of TPS and standard parallel tempering can provide a similar effect [63], though in RETIS there are no additional simulations at alleviated temperatures needed. In contrast, the [(i − 1) + ] ↔ [i + ] swapping move is very inexpensive as it does not require any force evaluations. If the move is accepted, it provides a new path for both the [(i − 1) + ] and the [i + ] path ensemble. In addition, the accepted swapping moves provide generally paths that are more decorrelated from the previous path than a shooting move.
The swapping moves between the [0 − ] and [0 + ] path ensembles are done by exchanging the end and start points of the paths and extending those forward and backward in time respectively. Despite not being a free move, the barrierless diffusion within the reactant well of the [0 − ] paths, followed by exchange moves, will basically feed the [0 + ] ensemble with fresh initializations. This facilitates the decorrelation of the MC sampling even further and orthogonal barriers can be avoided without having to cross them in any of the path ensembles. An additional advantage is that this method works in case where parallel tempering is not effective, that is when barriers are mainly entropic in nature.

IV. PERMEABILITY FROM RETIS SIMULATIONS
The rate k in Eq. (6) (with unit 1/time) describes the kinetics of a process, while the permeability P of a membrane (with unit length/time) describes the transport kinetics through the layer, and clearly a link between k and P is expected. However not only the units, but also the formalism for permeability and rates are somewhat different. This paper will derive the correct connection between the TIS framework for rare events and the permeability.
The progression of permeation is most easily measured by the z coordinate orthogonal to the membrane as the RC. Despite this simple order parameter choice, some details ask for attention. In case of periodic boundary conditions, common in MD simulations of membranes or porous catalytic crystals, the question arises how it may be detected whether a molecule crossed the membrane, or whether it simply circled back to the other side of the membrane through the water phase because of the periodic boundary condition in the z-direction. Therefore we first derive the connection between P and k when the solvent phase is bounded on both sides (e.g., by hard walls as in Fig. 1) in Sec. IV A, and second we derive a relation for the unbounded system when either periodic boundary conditions or an infinite particle bath [64] are applied in Sec. IV B.
A. Connecting permeability and rate Figure 1(b) shows a hypothetical long MD equilibrium run. The timescales are not very realistic as we would in practice expect that thousands of crossings with λ A would proceed a transition to state B. Yet, it will be used to show how we can subdivide an ensemble average into different regions A, B or A, B as follows.
The long trajectory can be seen as a series of visited phase points. Here, we assume that the MD equilibrium run is in fact sampling the distribution of interest. The MD integrator is hence coupled to some kind of thermostat or barostat, whenever this distribution is different from the NVE ensemble. The MD integrator gives a trajectory, which is effectively a series of phase points (time slices) because of the discrete time step t in the numerical integration, such as the velocity Verlet integration algorithm [65]. Given ergodicity, the phase points will be visited with the correct relative probability when the trajectory length T goes to infinity. The trajectory phase points can be divided into subsets corresponding to the regions A, M, B in Fig. 1(b) based on the value of the order parameter in each phase point. The ensemble average becomes where M is the membrane region, previously referred to as no man's land in the context of TIS and RETIS. Here, p X is the probability that the system is in state X . Given that the ergodicity makes the time averages respect the relative probabilities, it follows that p X = T X /T , where T X is the total time spent in state X , and T = X T X . In Eq. (9), the notation . . . A refers to the ensemble average over all trajectory phase points associated to A. As mentioned earlier in Sec. III, the trajectory phase points can also be assigned to either overall state A or B. The assignment to A or B is generally not based on the evaluation of the order parameter in a single phase point but rather based on the series of phase points in the trajectory (history dependent). The ensemble average can be divided in two contributions, Let us compare the three different flux terms: k in Eq. (6), f A in Eq. (7) and J + in Eq. (2). Here, f A provides the frequency of state changes from A to M within the overall state ensemble A. J + on the other hand measures all crosses from left to right along the full region M of all permeants, per time and per membrane surface area σ . The rate constant k finally also measures the frequency of full crossings like J + but with the same time normalization T A as for f A . In summary, One could imagine another flux definition, similar to f A , but where the denominator is the total time T of the long equilibrium run rather than the time spent in A, The appearance of the factor p A gives f A the interpretation of a "conditional" flux compared to f . The flux f is however not accessible in a RETIS simulation and will not be further discussed. If all of the N p permeating particles are identical, (13) 033068-6 and we can relate J + with f A and k as follows: where P A (λ B |λ A ) is the previously introduced crossing probability. All terms in Eq. (14) are measurable in a RETIS simulation except p A , the probability of overall state A.

RETIS simulates the
Further, we can use the subdivision in global states in Eq. (10), giving where z t is the z coordinate of the target permeant, and (ρ ref ) A refers to the probability density within the A ensemble. The factor δ(z ref − z t ) B is zero, since z ref is located in stable state A and z t does not lie in A for any trajectory phase point that is assigned to B (see Appendix A). Consequently, we can write Substitution of Eqs. (17) and (14) in Eq. (2) gives This is the first theoretical result connecting P and k. Equation (18) shows that the B path ensemble needs not be simulated to find the permeability. In contrast, the direct counting method of Eq. (2) In practice, one would count the number of steps that the z coordinate is inside the interval, and divide this by z and by the total number of steps that the system is part of A paths, assuming a constant MD integration time step t. If the reference location is in a region where the free energy profile F (z) is flat, then T ref scales linearly with z, and (ρ ref ) A will in principle not be affected by the chosen interval length z.
On the other hand, the conditional flux is the number of crossings through λ A in the positive direction divided by the time spent in T A , where, as said earlier, N [0 + ] is the number of [0 + ] trajectories that can be cut out from the long equilibrium trajectory. Substitution of the two previous equations into the permeability in Eq. (18) makes the T A drop out, and gives a practical expression for P, This is the second expression linking P with RETIS quantities. It gives the insight that P depends on time spent in the reference region, but not explicitly on time spent in the [0 + ] nor [0 − ] ensemble. However, Eq. (21) still refers to the quantities N [0+] and T ref , which are obtained from a long equilibrium simulation. In the next last step, the conversion to path ensemble averages is made. With T X the time spent in a path ensemble X and N X the number of trajectories in the path ensemble X , the average path length is given by τ X = T X /N X . An advantage of τ X is that it is in principle independent of the simulation computer time, i.e., if the number of simulated paths in the ensemble is doubled, the average τ X will not change, which gives τ X an intrinsic meaning. We would like to stress that a path ensemble average like τ X is a property that is averaged over all paths in the path ensemble, and it differs from the common phase space average, denoted as . . . Y , where the average is taken over all phase points within the ensemble Y .
The conversion to path averages is now simply executed by introducing factors For the probability density, this gives while the conditional flux is converted to (see Ref. [22]) This is the third theoretical result expressing P in terms of intrinsic RETIS quantities. Equation (24) is an expression that can be fully computed in a RETIS simulation since P A (λ B |λ A ) is standard output of this approach, while τ ref,[0 − ] can easily be obtained from the paths generated in the [0 − ] ensemble by histogramming the z coordinate. Figure 1 is however not typical for an actual membrane system that has no natural confining energy barrier that prohibits the permeants to drift far away from the membrane. To deal with the situation of an unconfined system, we derive in Sec. IV B how a λ −1 interface can be used to constrain the system in a way that the permeability coefficient can still be determined exactly, as if the system were unrestrained.

B. The λ −1 interface
The previous section linked the rate constant to a permeability, though the system depicted in Fig. 1 is not exemplary for a permeation system where we imagine a membrane in a near-infinite solvent. The rate constant k that follows a single particle will naturally decrease with the amount of solvent added to the model system since the target particle will spend more time away from the membrane. The permeability P is insensitive to this since it does not depend on the size of the simulation box. Increasing the solvent while maintaining constant concentration automatically implies an increase in the number of permeants. This increase cancels the decrease of the rate for permeation of the individual particles. Still, as RETIS computes rates rather than permeation directly, some confinement is required in practice. This confinement is achieved by introducing an extra interface λ −1 . With the introduction of this interface the overall state A is reduced to the λ > λ −1 region, and is denoted A . The same can be done at the product region side, but this is not essential for the A → B rate calculation. The new state division is depicted in Fig. 2, where the outer regions are here called freeze-time zones which will not be accessed in this adaptation of the RETIS algorithm. A prime is added to show that the states A , B , A , and B contain fewer phase points than A, B, A, and B, respectively. The ensemble average in Eq. (10) is updated to reflect the additional freeze-time zones,  Fig. 1(b)]. Whenever the system enters a freeze-time zone, the time T is stopped and continued whenever it exits this zone. The trajectories that can be cut out from the [λ −1 , λ 0 ] interval constitute the [0 − ] path ensemble. This ensemble is different from the [0 − ] ensemble since its paths can start and end at λ −1 in addition to λ 0 . Figure 2 shows again a hypothetical unrestrained equilibrium MD run just like in Fig. 1(b) with a free energy surface that is flat at either side of the membrane. Conceptually, the existence of an equilibrium in such an infinite system is not obvious as the ergodicity hypothesis implies that for instance ensemble averages are identical to time averages of an infinite equilibrium run. However, if the partition function diverges, even an infinitely long equilibrium run might not visit all the relevant phase space regions. This issue can be solved conceptually by taking the infinite limits for space and time in a controlled way. That is, we consider the potential of Fig. 2 as a special case of the potential shown in Fig. 1, but with the vertical-like increase of the free energy occurring at z = −W and +W . By letting T → ∞ and W → ∞ such that W 2 /(T D) → 0 with D the average diffusion constant, it can be shown that the ergodicity hypothesis holds. While this solves the conceptual problem whether we can assume equilibrium statistics, it does not solve the practical problem that the rate k is still zero in this limit.
The additional λ −1 interface, however, ensures that the overall state A becomes finite and that anything that happens in the freeze-time zone can be ignored as if the stopwatch, measuring T , is paused and statistics are not updated each time that the trajectory enters that region. Now that the rate k and the flux f A can be computed from this long equilibrium run using the new definition for the overall state A , the same counting strategy applies: it is a number of crossings/transitions divided by the time spent in overall state A . The only difference is that, besides time spent in overall state B , also the freeze-time zone is ignored in the 033068-8 normalization. The path ensemble [0 − ] is however changed to [0 − ], which statistically corresponds to the same ensemble one would obtain by cutting out the trajectory segments between λ −1 and λ 0 from the equilibrium run. Specifically, it is no longer valid that N [0 − ] and N [0 + ] are equal. We will therefore introduce the parameter ξ to measure the mismatch between number of paths in these two path ensembles, Here N →R,[0 − ] is the number of [0 − ] paths ending at the right side (at λ 0 ) that can be cut out of the long equilibrium run. It is Fig. 2. It follows that ξ < 1.
Similarly, h →R is the characteristic function that for each path provides the output 1 if the path ends at the right and 0 if the path ends at the left, irrespective of the starting point. In order to reformulate the expressions for f A and P, we first update the expressions for T A and T ref to link them to average path lengths τ X , which are intrinsic quantities of a given path ensemble.
The time T A can be written to be proportional to The time T A is smaller than T A since A comprises fewer phase points than A (the clock is stopped). The time T ref spent in the reference interval is not affected by the presence of the λ −1 as the reference interval is located between λ −1 and For P in Eq.
The last expression is the central expression that links the permeability with thermodynamic averages that can be computed in a RETIS path sampling simulation. The parameter ξ is which leads to an alternative expression for the permeability Equation (32) combines path ensemble averages (ξ, (30) is only based on path ensemble averages. All quantities in Eqs. (30) and (32) are intrinsic, meaning they are expressed as averages over the paths. This means that the same expression is applicable when doubling the number of paths, i.e., running the MC steps in the path ensembles longer. This will not affect the absolute value of any of the terms provided in Eq. (30), but will naturally increase the accuracy of their numerical estimates.

V. RC FOR PERMEATION WITH PERIODIC BOUNDARY CONDITIONS
In the previous section, we solved the problem to link the permeability with the rate constant in an infinite system by introducing the interface λ −1 and the rate constant k , which is nonzero unlike k. In most practical simulations the infinite system is represented by a system with periodic boundary conditions (PBC) which poses the need to properly determine the relative position of the permeant with respect to the center of the membrane, which defines the RC, λ. The simple minimum-image convention will generally not work since the RETIS trajectories can span the full [λ −1 , λ B ] region and a permeant located at λ −1 might actually be closer to the left periodic image of the membrane than to the membrane in the central image.
In order to properly deal with PBC, we first map the z i coordinates of all particles i in the system within the [−L z /2, L z /2] interval, where z = 0 is matched by convention at the center of mass of the membrane and L z is the box length in the z dimension, Here, z i is the z coordinate of particle i provided by the molecular simulation program. By this operation, the center of mass of the membrane is set at 0, while the center of the solvent slab is at ±L z /2. Here we assume that original coordinates z j of the membrane particles j are constructed such that the center of mass at the z axis can be computed without the need to add or subtract L z or multiples of it to any of the membrane particles. As the membrane is stable, we can assume that the above remains valid during the full simulation. That is, membrane particles might move over large distances only if all membrane particles move in cohort. The RC is given by the relative position z i of the tagged permeant i (target) with respect to the membrane. This implies that if we want to compute the permeability of oxygen through a membrane and our model system contains N p 033068-9 oxygen molecules, one of those will be selected and considered as our target permeant.
In some cases, it can be advantageous to select a collective RC such as the maximum value of the z coordinates of all permeants. This has the advantage that the rate increases which makes it less of a rare event and is therefore easier to compute. In addition, the collective RC facilitates the decorrelation of the sampling since the target permeant defining the RC can switch during the simulation. This strategy was for instance applied to study water dissociation where the RC was defined as the largest OH bond in the system [66]. Also in a recent paper on permeation by some of us [67] such an approach was applied to compute escape rates of permeants being trapped in a membrane.
The reason that we nevertheless choose here a RC based on a single target permeant is because the permeation problem is in some applications less of a rare event than for example water dissociation. This implies that for a system with many permeants there is almost always one of them in the membrane region. In fact, the membrane often does not correspond to a single peaked free energy barrier, but may have a well in the middle where permeants get temporary trapped. This makes a collective RC impracticable for describing stable state A. In addition, there are technical and theoretical problems associated to such a collective RC for the calculation of τ ref,[0 − ] that is needed for Eq. (30). The implementation of the collective RC would also be more cumbersome with periodic boundaries than when the RC just depends on a single target permeant. Instead, we can recover the advantages of the collective RC by the introduction of the new MC moves discussed in Sec. VI.
Still, also in the target permeant approach, care has to be taken with periodic boundaries when the relative position of the target permeant with respect to the membrane has to be determined. As a start, we obey the convention λ −1 <λ 0 < λ n .
In addition, the RC should change continuously along the sequence of time slices of a complete path, i.e., it should not suddenly jump by a value equal to L z which could lead to untrue transitions between states. Trajectories end and start with a time slice outside the [λ −1 , λ n ] region and the RC of these points define the states at which they start and end, so the 'jump-free' interval needs to be extended slightly beyond the [λ −1 , λ n ] interval. Based on the above, the safest option is to allow the jump in the RC to occur in the mid-point between λ n and the periodic image of λ −1 (at λ −1 + L z ). This yields for the final RC where z ∈ [−L z /2, L z /2] is the z coordinate of the target permeant following the convention of Eq. (33).
The RC as a function of the target permeant's position within a periodic system is shown in Fig. 3. For NPT simulations with fluctuating box dimensions it might be convenient to define the RC relative to the box size along z: λ NPT = λ/L z where λ is still defined by Eq. (34). The only difference is that L z is the instantaneous box length, which is a variable instead of a constant.

VI. NEW MC MOVES IN PATH SPACE
The choice to select a single target permeant instead of a more collective RC and allowing the target permeant to cross the membrane in just one direction (left to right) can lead to a somewhat restrictive sampling speed in comparison with a more collective RC. In the next two sections, we show how to remove these restrictions by introducing two new MC moves for the [0 − ] ensemble without the need to alter the definition of the RC or the setup of interfaces.

A. Target swap move
As discussed above, the RC is determined by the z coordinate of a single target permeant. The other permeants basically contribute to the environment around the target permeant like any of the other nonpermeating particles in the system. The occurrences of these particles crossing the membrane do not have to be counted as they are part of the natural fluctuations in the environment.
Especially when the membrane is not uniform containing different channels through which the permeants could transfer, it would still be advantageous to utilize the contributions of all the permeants. Some regions in the membrane that are easier to penetrate could be blocked by a nontagged permeant. Waiting for a swap through diffusion of both permeants within the bulk might take a long time. In order to speed up this process, we design a MC move in path space that allows for a swap without diffusion, but by simply reassigning the target.
The target swap move is explained in Fig. 4 and a stepwise description of the algorithm is given below.
(1) Assume the old path (upper panel in Fig. 4) has length L (o) (including start and end points) and is represented by time frames numbered from 1 to L (o) .
(2) For each frame ( (4) Select the permeant/time-frame combination corresponding to the i-th count at step 2. We will call this the new target permeant and the frame index at which this count was registered is from now on called j.
(5) Starting from time slice j, the path is followed (as detailed below) backward until we detect a frame in which the new target has a position outside the [λ −1 : λ 0 ] interval. We call the number of backward steps n b .
(6) Starting from time slice j, the path is followed (as detailed below) forward until we detect a frame in which the new target has a position outside the [λ −1 , λ 0 ] interval. We call the number of forward steps n f . (7) The new path length is n b + n f + 1 which starts with frame index j − n b (which can be negative) and ends with frame index j + n f . Renumber all frame indices by adding − j + n b + 1 to each frame index. The new indices now run from 1 to n b + n f + 1.
Then, using Metropolis-Hastings rule [68,69], the acceptance probability can be written as where p(o) and p(n) are the probabilities of the old and new path, respectively, and P gen (X → X ) is the generation probability to generate path X from X . As the identity of the target has no effect on the path probabilities, the probabilities p(o) and p(n) are essentially the same except for stochastic force terms related to extending or shortening of the path. However, these terms cancel in Eq. (37) as they are also part of the generation probabilities [54]. The only remaining terms that do not cancel are, hence, the selection probabilities for selecting the permeant/time slice.

B. Mirror move
In the case that the membrane is symmetric, transitions through the membrane from left to right and from right to left are statistically indistinguishable within an equilibrium sampling. It is then favorable to count transitions in both directions [10,12]. For instance, the direct counting method described in Sec. II uses this strategy to improve statistics.
One obvious way to include two-directional transitions could be achieved by defining the RC as the absolute distance |z| between the target permeant and the center of the solvent slab (see also Ref. [67]). However, as explained in Sec. V, this can lead to an overlap in the z-coordinate space. While this could still be solved by letting the RC value depend on the history of the path, i.e., the solvent slab's periodic image to be considered is determined by the minimum distance image at the start of the path, this becomes problematic when the path's history is not yet fully determined. For instance, this is the case when a shooting move is carried out or when the target swap move implies that some backward in time integration is required.
The mirror move in path space (see Fig. 5) achieves the same versatility of the two-directional approach in the counting method without having the problems discussed above. The mirror move in [0 − ], which is always accepted, mirrors the whole system with respect to the membrane center. The z coordinates of every particle are mirrored and the z-component of the velocities are multiplied with −1. Because of the periodicity, this mirror move is equivalent to mirroring the whole system with respect to the midpoint between the membrane and its periodic image, which is, loosely speaking, the midpoint of the solvent slab. By construction, this midpoint of the solvent slab lies in the middle of the [λ −1 , λ 0 ] interval related The mirror move swaps the roles of the λ 0 and λ −1 interfaces. This requires that λ 0 and λ −1 are placed at the same distance from the mirror plane. To achieve this, consider the membrane's center-of-mass position at z = z mem and its left periodic image at z = z mem − L z . The distance between λ 0 and membrane should equal the distance between λ −1 and this left periodic image membrane. In other words, λ −1 should be placed such that Equivalently, the interfaces λ −1 and λ 0 should have the same distance to the midpoint of the solvent slab at z = z mem − L z /2. Since we applied the periodic coordinates z as defined in Eq. (33), we have z mem = 0 and therefore λ −1 = −(λ 0 + L z ). Given the positioning of λ −1 and λ 0 according to Eq. (38), the mirror move simply mirrors the z coordinates of every particle in the system at every time slice with respect to the 033068-12 plane z = −L z /2. Hence, it proceeds according to the following steps.
(1) Let k run over all time slices of the path. This includes the start and end point of the path and might include stored time slices (see Sec. VI A).
(3) Consider the z coordinate of permeant j and its velocity component at discrete time step k: z j (k) and v z, j (k).
(5) Accept the new path. If the old path started/ended at λ −1 then the new path will start/end at λ 0 and vice versa.
For code-technical reasons we implemented a slightly different approach where we did not alter the coordinates or velocities, but instead the definition of the RC function [Eq. (34) with z replaced by −z]. The system was hence assigned an additional flag which indicates whether Eq. (34) has to be used with the plain z coordinate of the target permeant or with −z. This pragmatic choice made it easier to use PYRETIS [57,58] with external MD engines, as these might have very different ways of altering the coordinates and velocities. The flag is also exchanged in the replica exchange moves of the RETIS algorithm.
The new moves are only implemented for the [0 − ] ensemble. For the mirror move, this is because this is the only ensemble where paths can start at both the left or the right hand side. In addition, the target swap move is not expected to give a high acceptance for the [i + ] ensemble when crossing λ i is a rare event.

VII. NUMERICAL RESULTS
The theoretical derivation for the permeability calculation from RETIS has been implemented in the python based opensource code PYRETIS [57,58]. First, a one-dimensional toy system was constructed where a Langevin, Brownian, or deterministic Newtonian particle permeates through a medium with or without barrier. For some limiting cases, an analytical expression for the permeability is available, which can serve as a validation of the new RETIS permeability formula. Second, a two-dimensional membrane was simulated with periodic boundary conditions, where permeants can pass the membrane through two different permeation channels. This last system is used to illustrate efficiency of the new Monte Carlo moves.

A. One-dimensional system setup
For simplicity, the membrane is located symmetrically around z = 0, in the region |z| < a, with h = 2a as the membrane height. The effect of the membrane is modeled by an external cosine-shaped potential that acts on a single permeant particle, Here, V 0 is barrier height. This membrane model ensures that the force on the particle is continuous at the membrane borders z = ±a. The harmonic potential 1 2 k harm (|z| − b) 2 is added only to allow the system also to be studied by a reference simulation without λ −1 interface. In most simulation setups, the k harm parameter is set to zero, which reflects the real physical situation for a permeation system that is unbounded at either side of the membrane. The dynamics of the permeant is governed by either Langevin dynamics, Brownian motion, or deterministic dynamics.
In the Langevin dynamics, the permeant experiences both friction, inertia effects, and random collisions with a degree that is controlled by the friction parameter γ . The friction constant γ (unit 1/time) of the particle relates to the particle's diffusion constant as D = k B T /(mγ ). Note that some other  textbooks use a friction coefficient with unit mass/time. This refers to an alternative definition of the friction coefficient γ = mγ . The two other types of dynamics can be viewed as limiting cases of the Langevin dynamics.
A Brownian particle propagates in discrete steps without memory under influence of a random Gaussian force and the force −dV/dz exerted by the potential. It can be considered as the overdamped limit of the Langevin particle, when γ → ∞.
A Newtonian particle on the other hand propagates deterministically in time according to the equations of Newton. It has inertia but undergoes no friction nor random collisions, and energy is conserved. In the RETIS algorithm, the effect of temperature is then only present in the MC moves when a new constant-energy path is created from an old constant-energy path. The detailed balance MC procedure allows the energy to change between paths such that the overall path ensembles are canonically distributed. This simulation set up reflects a system that is so weakly coupled to a thermostat that the dynamics of a single crossing basically occurs at a constant energy (NVE, microcanonical). However, at the much longer timescale between crossing events, the energy could change. The Newtonian particle can be seen as the low friction limit of the Langevin particle, when γ → 0.
The driver for solving the equations of motion was the internal MD-engine of PYRETIS in our toy system [57]. It is also possible to let PYRETIS manage the path ensembles book keeping, while it calls simulation programs, such as GROMACS, OPENMM, or LAMMPS, to drive the molecular dynamics [58]. Throughout Sec. VII, reduced units are used in which mass m, the Boltzmann constant k B , and temperature T are equal to unity (m = k B = T = 1 in reduced dimensionless units). Potentially, several physically realistic systems could be mapped on the model presented in Eq. (39) by tuning appropriate units of energy, length and mass. We examined the model system Eq. (39) using the following parameters: a = 0.1 (h = 0.2), k harm = 100 or 0, m = 1, T = 1, and integration time step t = 0.002. The friction coefficient γ had values 0.1, 5, 10, 20, 40, 60, 80, or 100. The Brownian dynamics propagates through configuration space at discrete steps with a displacement that is governed by a step-size parameter. This step-size parameter can be associated to a t in an equivalent Langevin simulation for a given mass and friction. In our case, the step-size parameter was set such that the time between steps was the same as t of the corresponding Langevin simulation with γ = 100.
In the RETIS simulations, the number of cycles was set to 20 000 (this is the number of MC moves in the path ensembles), the RETIS swapping move frequency was 0.1, the shooting frequency 0.45, and the time reversal move frequency 0.45. The position z was used as the order parameter λ. Three interfaces were used, λ 0 , λ 1 , λ 2 , which were chosen at λ = −0.1, 0, and 0.1. The additional interface λ −1 was chosen at λ = −0.2 and varied to a few other locations in the calculations of Table I The reported error is based on block averaging and error propagation rules assuming independence of the different path ensemble simulations. As the latter assumption is not fully valid due to the replica exchange moves (Ref. [22]), we also estimated the error on P via 10 independent realizations, shown in the last column. . The time spent per path per length would become larger when the friction increases. The crossing probability is lower than for the deterministic particle, which makes the Langevin permeability lower than the deterministic value 0.399 (see Appendix B). Figure 6 compares the computed permeability for the potential Eq. (39) with V 0 = 0 (flat), 0.5, and 1 (cosine barrier membrane) for the three types of particle dynamics. In the Appendix we derived analytical expressions for P based on the Smoluchowski equation and Kramers' expression for Brownian dynamics and Langevin dynamics, respectively. These theoretical curves are shown in the same graphs. The analytical result for deterministic dynamics can be obtained by taking the limit of Kramers' expression for γ → 0. The validity of these theoretical results relies on different kind of approximations. The Smoluchovski expression [Eq. (B1)] is reliable for high friction and low to high barriers, while Kramers' theoretical result [Eq. (B6)] is the reliable reference for high barriers and low to high friction. There is henceforth a blind spot in system parameter space: dynamics with low friction and low barriers is poorly described by both theories.
Indeed, consider the flat potential membrane (V 0 = 0) in Fig. 6(a). It shows good agreement between the theoretical Smoluchowski curve, P = D/h, and the simulated Langevin results for the large γ values. Also the computed values for Brownian and Langevin at γ = 100 agree, as Brownian dynamics can be seen as the high friction limit of Langevin dynamics. However, for low friction, the Smoluchowski curve and the numerical Langevin results deviate. The Kramers prediction of the permeability is zero for any γ > 0 and thus also fails to approximate the Langevin numerical results. Only the limiting case γ → 0 with γ 2 /V 0 → 0 has a nonzero solution equal to √ k B T /(2π m). The Kramers curve in Fig. 6(a) actually shows the theoretical results for V 0 = 10 −4 to visualize this limiting case. The deterministic √ k B T /(2π m) limit to the permeability is indicated by the dashed horizontal line. Our numerical data on deterministic dynamics and the Langevin result with γ = 0.1 agree with this limit.  6(b) and 6(c) show the numerical and theoretical curves for systems with a membrane barrier, V 0 = 0.5 and 1.0, respectively. For these membranes, the two theories agree in the large friction regime. In the low friction regime, the Smoluchowski expression is not a good approximation of the Langevin dynamics. The numerical Langevin simulations agree with the Kramers' curve for all values of γ . The deterministic and Brownian dynamics simulations also agree with Kramers' expression in the limiting γ = 0 and γ = 100, respectively.

C. Two-channel membrane system setup
A two-dimensional system is constructed that mimics a membrane barrier through which particles can permeate through two competing pathways. It could for instance represent a membrane with two transmembrane protein channels. Three noninteracting Langevin particles are subjected to the potential V (y, z) = e −cz 2 V 1 + A + A sin 2π y L y + B + B cos 4π y L y , The membrane is located in the center of the unit cell around z = 0, while V (y, z) is approximately zero far away from the center due to the factor e −cz 2 (see Fig. 7). Periodic boundary conditions are applied, where the system is made periodic in the z direction with a period [−L z /2, L z /2] and in the y direction the period is L y . Particles can permeate the membrane through two channels: one channel at about y = −0.25L y with barrier height V 1 and another channel at about y = 0.25L y with barrier height V 2 . The maximum barrier height is V max .
Reduced units are used as in the one-dimensional case. The parameters in our simulations were V 1 = 10, V 2 = 11, V max = 20, c = 1, and L z = L y = 6. The PYRETIS simulations were run with three Langevin particles with settings t = 0.02, γ = 5, T = 1, and m = 1.
The order parameter is the reduced z coordinate of the target permeant: λ = z j if permeant j is tagged ( j = 1, 2, 3). (1) TIS: In TIS, there are no swapping moves between the path ensembles. The MC moves are the shooting move and time reversal move, with equal frequency 0.5. Differently to standard TIS, the sampling of state A was done using the [0 − ] path ensemble simulation and not via MD.
(2) RETIS: In standard RETIS, the swapping move of paths between the path ensembles is also allowed as an MC move. The swapping move frequency was 0.5, shooting move frequency 0.25, and time reversal move frequency 0.25.

033068-15
TABLE II. Numerical results (reduced units) for permeability P of 3 Langevin particles permeating through a two-channel membrane. Standard error from block averaging and error propagation between brackets. Each of the three particles could be selected as the target permeant when performing a target swap move. The swapping move frequency was 0.5, time reversal move frequency 0.25, mirror move frequency 0.05, target swap move frequency 0.05, and shooting move frequency 0.15.
After an equilibration run of about 1600 MC moves, the analysis was performed based on a production run of 35 000 MC moves. Table II shows the permeability together with the calculated variables that enter in Eq. (30), ξ , τ ref,[0 − ] / z, and P A (λ B |λ A ). Note that the RETIS results on ξ and τ ref,[0 − ] / z are somewhat off compared to TIS and RETIS*. Naturally, for this symmetric system ξ = 0.5 is the exact result which agrees with TIS and RETIS* while RETIS is 8% too high. This is due to the relatively wide region between λ 0 and λ −1 and the lower frequency of shooting moves in the RETIS simulation compared to TIS. In RETIS, 50% of the moves are swapping moves (replica exchange moves between path ensembles), which are very useful to improve the sampling of the barrier region, but not necessarily help the exploration of the water phase. Both the swapping and time-reversal moves are unable to generate a λ −1 → λ −1 path from a λ 0 → λ 0 path and vice versa. In RETIS*, the target swap move and the mirror move repair this weakness even if these moves only represent 10% of the executed MC cycles.

D. Two-channel membrane: analysis
The crossing probabilities P A (λ B |λ A ) and the permeabilities in Table II give quantitative good agreement within about 10%. Given a barrier of at least 10 k B T this is a notable result. However, even more difficult challenges for a simulation method in this system, are (i) the ability to sample transitions through both channels, and (ii) to achieve this with the correct ratio. Since the channels' barriers only differ by 1 k B T , both permeation routes are competing, but successful permeation transitions are expected to proceed via the lowest barrier channel in about 73% of the cases. Getting this ratio right is extremely challenging for any rare event method. Figure 8 shows the distribution P(y * ) of the orthogonal coordinate at the first crossing with λ i , y * , for different path ensembles, [i + ], i = 0, 1, . . . , 10. The crossing point y * of a path is indicative of the channel visited by that path. The distributions show that for TIS, all y * crossing points in the ensembles [6 + ] and higher are in the V 2 channel, while for RETIS and RETIS* the other channel is visited as well in all ensembles. This clearly demonstrates the deficit of the shooting move. The chance for this move to generate an acceptable path is highest when the shooting point is chosen on the barrier region, close to λ i for ensemble [i + ]. However, switching between channels can practically only occur if the shooting is initiated from the well region. The fact that TIS got stuck in the high-energy channel, rather than the low-energy channel is purely accidental reflecting the memory of the initial path that was used to bootstrap the simulation. The TIS result (Table II) is lower than the other values, as one might expect based on its bias towards the high-barrier channel. Yet, since the crossing probability up to λ 6 is based on the progress through both channels, the TIS permeability is still rather close to the RETIS and RETIS* results.
Provided ergodic sampling, TIS and RETIS should be capable to sample nontrivial multiple-channel systems where splitting based methods, like FFS and AMS, would fail. An example of such a case is a system with two channels in which the lowest-barrier channel goes initially much steeper uphill than the channel with the higher barrier [62]. However, as is clear from Fig. 8, the TIS simulation is not ergodic since it is not able to switch between channels for ensembles [6 + ] and higher with just the shooting move. For this academic model, this aspect could be repaired using nonlocal shooting moves in which not only the velocities, but also the configuration point is changed by a nonlocal displacement. Such a move, however, would have vanishingly low acceptance in a realistic condensed matter system as nearly every attempt will lead to a molecular overlap.
The RETIS and RETIS* simulations, however, are able to sample both channels and get the ratio between low-and highbarrier pathways at least qualitatively correct. The replica exchange swapping moves allow the exchange of paths be- If we examine the height of the distributions in Fig. 8, we see that both RETIS and RETIS* predict that the majority of transitions will pass via the low-barrier channel. The RETIS simulation, however, seems to overestimate the preference of the V 1 channel, especially when the [10 + ] ensemble is considered. Integration of exp(−βV (y, λ 10 )) over y along positive and negative values indicates a 2.54 higher probability to be in the negative y-range. Even if these relative probabilities deviate a bit from the distribution of first crossing points y * , this deviation is expected to be marginal for this model system.
To further analyze the effectiveness of the MC schemes we analyzed the number of channels switches observed in the path ensemble simulations. For path ensemble [i + ] each path is assigned to channel V 1 , channel V 2 , or neither of the two, based on the first crossing point y * with λ i . It is assigned to belong to the V 1 -channel if −2.5 < y * < −0.5 and to the V 2 channel if 0.5 < y * < 2.5. A channel switch is counted when the MC move produces a V 2 channel path while the V 1 channel was most recently visited, and vice versa. Figure 9(a) shows the number of channel switches that are calculated via this approach for different path ensembles. The TIS results are magnified by a factor 20 for visualization as this approach shows dramatically less channel switches than RETIS and RETIS*. This shows that replica exchange (swaps between path ensembles) is absolutely necessary for efficient sampling. From [0 + ] to [5 +  The difference in channel switches seems negligible between RETIS and RETIS* up to λ i = −1. After that point, RETIS* seems to produce significantly more switches. This is remarkable since the extra moves, the mirror and target swap move, are only executed in the [0 − ] ensemble and its effect on the [10 + ] path ensemble is only indirect via the replica exchange moves. It requires at least 10 path ensemble swaps to process any information from [0 − ] up to [10 + ]. Still, the effect is most noticeable for the last seven path ensembles, but hardly before. Our conjecture is that the channel switches due to the mirror move and certainly due to the target swap move are more effective in decorrelating the ensemble. A channel switch is likely more effective if it enters the new channel along its central line and indeed this happens more often with the target swap move than with the shooting move. Also, the number of channel switches does not tell the full story. If two path ensemble [i + ] and [(i + 1) + ] are at different channels and then make a lot of successful replica exchange moves solely between each other, this will yield a lot of channel switches. However, the effectiveness in decorrelating the sampling will be modest.
To examine this decorrelation, Fig. 9(b) plots the ratio of the number of paths in the V 1 and V 2 channels as a function of the number of MC moves in the [10 + ] ensemble. The TIS curve is a flat zero line as it is stuck in the V 2 channel for the full simulation. Comparing RETIS and RETIS*, we see that RETIS* converges much faster and approaches the predicted value of 2.54 that was obtained from the numerical integration. Figure 9(c) shows a zoom of the curve together with the index for the target permeant. These running averages for RETIS and RETIS* have a typical sawtooth shape though the latter is able to flip much more frequently the slope of the curve. The zoom in Fig. 9(c) shows that such a flip often coincides with a change of the index for the target permeant. This indicates that the target swap move has a strong influence in the overall sampling even if it is only applied in the [0 − ] ensemble.

VIII. CONCLUSION
In this work, we derived a formula for the permeability based on path sampling quantities that can be determined in a RETIS simulation. As the idealized permeation model represents a membrane inside an infinite solution, the RETIS path ensembles require an adaptation via the introduction of an additional interface prior to λ 0 , called λ −1 , and a newly defined path ensemble [0 − ] that replaces the [0 − ] path ensemble. The resulting approach is exact and does not depend on the positions of the interfaces including λ −1 . Their positions are therefore set to optimize efficiency.  In addition to this theoretical derivation, we introduce a few algorithmic developments such as a consistent way to define the reaction coordinate for permeability whenever periodic boundary conditions are applied and two additional MC moves that mainly operate in the new [0 − ] path ensemble. One of these new MC moves is the mirror move which can be applied whenever the membrane is symmetric. The other MC move is the target swap move and can be used when more than one permeant is present in the simulation model system.
Our new theoretical formulation and algorithmic developments have been implemented in the open-source PYRETIS code [57,58], and it was successfully tested on a onedimensional Langevin system for which analytical results exist. After this, a challenging two-dimensional model membrane with two competing permeation channels was simulated to test the effectiveness of the new MC moves. These simulations show that the replica exchange moves are essential to simulate this system as the plain TIS method gets trapped inside a single channel. The inclusion of the two new MC moves considerably improves the sampling efficiency even further, as is clear when inspecting the relative transmission through the two channels. This noticeable difference is surprising given the fact that the new moves only operate in the [0 − ] ensemble and it takes at least 10 replica exchange moves to transfer the effect of these moves up to the last path ensemble [10 + ]. Still, the direct relation between the improved efficiency and the new MC moves was demonstrated by a correlation between the channel-switches and changes of the target permeant's identity.
The theoretical derivation in this paper is valid for all kinds of microscopically reversible dynamics (e.g., deterministic Newtonian dynamics, Langevin, Brownian, Nosé-Hoover, etc.). Besides the standard ergodicity hypothesis, it does not rely on any further assumption nor approximations. This implies that our approach will in principle give the same value for the permeability as the direct counting method based on brute force simulations, but orders of magnitude faster.
Our approach has the great advantage that the Markovian assumption of memoryless hopping between interfaces (see Sec. I) is not needed like the approaches based on milestoning [42][43][44][45]. Especially for large molecules the permeation process is often driven by nontrivial membrane fluctuations such that the projected dynamics on a one-dimensional coordinate gets a memory dependent character. Since RETIS is inherently non-Markovian in its description, it allows a much broader range of applications. An interesting route of thought could be the combination of our facilitating RETIS framework with the high-throughput methods that efficiently scan chemical space [70] On the other hand, a milestoning type approach avoids the creation of full transition trajectories which can be computationally demanding when the transit time through the membrane is long. In this case, PPTIS [51] could be an interesting alternative. PPTIS avoids the sampling of complete transition paths like milestoning, but still maintains some of the history dependence. Alternatively, the exact non-Markovian character could be kept by alternating between short and long paths by means of stone skipping/web throwing [71]. Both PPTIS and stone skipping/web throwing can straightforwardly be implemented in our theoretical framework.
In conclusion, our permeability method presents a modelfree approach for the computation of permeability and it is expected to become a valid standard method when membrane crossings are rare events.

APPENDIX A: ENSEMBLE AVERAGES IN TRAJECTORY SPACE
At several instances in our article [e.g., Eqs. (6), (10), and (16)], we refer to phase space ensemble averages with the remark that these should actually be viewed as an average over "trajectory phase points." Even if this point has mainly conceptual importance, we will outline here its mathematical interpretation since it is yet underreported in literature. One convenient way is to refer to path space ensemble averages instead of phase space ensemble averages [72]. Here, the path X = {x 0 , x 1 , . . . , x L } can be viewed as a "chain of states" [73] with x i the phase point that is visited after i MD steps, at time t = i t with t the time step. From this, the path probability follows as where ρ(x 0 ) is the probability density of the initial state (x 0 at t = 0) of the path, and p(x i → x i+1 ) are the single time step transition probabilities. The latter are dependent on the type of dynamics. Mostly, we assume that ρ(·) is the equilibrium phase space density given by the Boltzmann distribution: ρ(x) ∝ exp(−βE (x)) with E (x) the total energy of phase point x. Actually, while P[X ] and p(x i → x i+1 ) are commonly referred to as a type of probabilities, it would have been more accurate to call these probability densities as well. Now, by expressing an observable f as a functional of X , the path ensemble average can be formally written as As we assume microscopically time-reversible dynamics, we can write [72] ρ( where x refers to the momenta-reversed phase point: if x = (r, v) with r the configuration and v the particles' velocities, then x = (r, −v). Here, it is also assumed that ρ(x) = ρ(x) for any phase point x. Applying Eq. (A3) multiple times on Eq. (A1), allows us to write alternative expressions for the path probability [59]: In the TIS and RETIS theoretical framework, the path concept is extended by including time slices before x 0 . In this view, x 0 is considered the present state, the principle phase point, while x i is a state in the future or in the past whenever i is, respectively, positive or negative. Hence, for a path (A4) Using Eq. (A4), we can in principle redefine all phase space ensemble averages, in which one integrates over x, as path ensemble averages in which one integrates over x 0 and additional phase points x i =0 with both positive and negative index via Eq. (A2) with dX = L i=−M dx i . Whenever the value of f is instantaneously available from the present phase point x 0 , the integrals over the additional phase points can be ignored, since they are unity.
Another way to generalize the ensemble average, which in some cases could be arguably more intuitive, can be derived from another perspective on the path object as stated by Crooks and Chandler [73]: "A stochastic trajectory can be defined by the chain of states that the system visits, but it can also be represented by the initial state and the set of random numbers, the noise history, that was used to generate the trajectory." [73]. As an example they show that the probability density of a one-dimensional Brownian dynamics path X consisting of L time slices can be written as where each ξ i is a Gaussian random number of zero mean and variance, the stochastic force acting at time between t = (i − 1) t to t = i t. Here, we deliberately shifted the indexing of the noise terms from 1 to L instead of the original [73] indexing from 0 to L − 1. The reason becomes clear when we introduce Eq. (A8).
Since the phase point of the system at t = t, x 1 , is fully determined by the first phase point and first stochastic noise term, we can write x 1 = φ(x 0 , ξ 1 ) with φ being the MD timestep integrator. Likewise, x 2 = φ(x 1 , ξ 2 ) = φ(φ(x 0 , ξ 1 ), ξ 2 ) etc. It is thus apparent that, when we add to x 0 all the information of the random noise sequence ξ 1 , ξ 2 , . . . , ξ L to make an "extended phase point" or "trajectory phase point," x = {x 0 , ξ 1 , . . . , ξ L } = {x 0 , ξ L }, basically every property of the system between t = 0 and L t becomes a function ofx, as if the dynamics would be deterministic.
Henceforth, for a general type of stochastic dynamics that proceeds via random noises that are drawn from a distribution p ξ (·), we can define phase space density of an extended phase pointx as So equivalently to Eq. (A2), by expressing an observable f as a function of an extended phasepointx, we write its ensemble average as This automatically becomes a standard phase space ensemble average with x instead ofx when f is not noise-dependent since all integrals over dξ i become 1. While the concept of an extended phase point is generally not explicitly referred to, it is often implicitly used. For instance, time-correlation functions are often casually introduced as C(t ) = a(0)b(t ) without being specific about the noise dependence. Based on Eqs. (A5) and (A6), we can rigorously define the ensemble average as an integral over extended phase space: where a and b are functions of the phase point of the system at the time under consideration, at t = 0 and L t, respectively. The absolute timescale is irrelevant here since we generally assume we are at an equilibrium distribution at t = 0 and the dynamics conserves this distribution, i.e a(0)b(t ) = a(t )b(t + t ) or any arbitrary t . Hence, the correlation function C(t ) becomes an ensemble average a(x)b(x; t ) where one just integrates overx and b is parametrically dependent on t in addition to its dependence onx. Comparing Eqs. (A1) and (A5), it is apparent that p(x i → x i+1 ) = p ξ (ξ i+1 ) for stochastic dynamics with ξ i+1 being the noise that forces the dynamics to produce x i+1 from x i ; x i+1 = φ(x i , ξ i+1 ). For deterministic dynamics, we can write p(x i → x i+1 ) = δ(x i+1 − φ(x i )). In addition, it is clear that the path interpretation and the extended phase point interpretation are equivalent; if one knows the initial phase point and the noise sequence, one knows the path X and vice versa.
As stated before, the TIS and RETIS theoretical framework requires the description of phase points before x 0 . This means that the "noise history" term by Crooks and Chandler to denote ξ L + = {ξ 1 , ξ 2 , . . . , ξ L } is now recoined as noise future while ξ M − = {ξ −1 , ξ −2 , . . . , ξ −M } is the actual noise history or noise past.
Then, equivalent to Eq. (A4), we can define the phase space density ofx = {ξ −M , . . . , ξ −1 , x 0 , ξ 1 , . . . , ξ L } as where the noise terms have a slightly different interpretations depending on the index being positive or negative. For i > 0, ξ i is the noise needed for φ to produce x i given x i−1 , while ξ −i is the noise needed for φ to producex i givenx −i+1 : φ(x i−1 , ξ i ) = x i and φ(x i−1 , ξ i ) =x i . Hence, the history of the path X follows from the negative noise terms as: , etc, again showing that there is a one-to-one relation between X andx. Based on Eq. (A8), it is now possible to define the probability of overall state A as: such that h −n A is a function of (x 0 , ξ −1 , ξ −2 . . . ξ −n ) which is a part ofx. The product term is simply 1 if none of the points In principle, we consider the path or random noise sequence to extend to infinite in both time directions [L → ∞, M → ∞ in Eq. (A8)]. However, as h −n A (x) does not depend on ξ i with i > −1 nor i < n, many noise integrals are simply 1 and therefore with dξ n − = dξ −1 dξ −2 · · · dξ −n . Now, suppose f is a function of phase space: f = f (x). Then, the ensemble average of f does not require the integration of any noise terms though the conditional ensemble average f A does as The division used in Eq. (10) can now be understood with this extended phase point picture in mind as If we consider the third line of Eq. (A15) separately, we can identify this, for a given phase point x 0 , as the chance that the stochastic dynamics needed exactly n steps backward in time to move outside no man's land, i.e., enter either stable state A or B. Since we assume that no point x 0 can be trapped into no man's land forever, the sum over n of this probability equals 1. Hence, Note that the we can use integration over x [Eq. (A13)] or x 0 [Eq. (A16)] interchangeably since ρ(·) refers to the equilibrium phase density that is time-invariant. So indeed, . . . = p A . . . A + p B . . . B like stated in Eq. (10). As a special case, we can take f (x; z) = δ(z − z t ) with z t being the z coordinate of a specific particle (the target permeant). Here, z t is a part of the system's phase point x that on its turn can be viewed as x 0 which a part ofx (recall that the phase space density is time-invariant). In addition, z is a parameter that specifies a reference region in configuration space. Note that the parametric dependence of f on z is not vanishing when taking the ensemble average since it is not a part ofx and therefore not integrated out. Therefore, we can write All ensemble averages in Eq. (A17) are in principle integrals overx, though in the first line an integral over configuration space would be sufficient since the integrals over momenta and noise terms are unity. In the last line of Eq. (A17), the integrals need to be carried out on the principle phase point x 0 and the backward noise terms ξ −1 , ξ −2 , . . . The integrals over the forward noises ξ 1 , ξ 2 , . . . are still unity. The delta-function is only nonzero whenever z t (x 0 ) equals z. This implies that if z is inside stable state A, then the product δ(z − z t (x 0 ))h B (x) is by definition zero; if z t (x 0 ) = z ∈ A then x 0 ∈ A and, therefore, h B (x) = 0. This explains the last equality of Eq. (16) Finally, for the rate in Eq. (6), the product of h A (0) and h B ( t ) should be evaluated. By defining the principle phase point to be x 0 , we havex = (. . . , ξ −2 , ξ −1 , x 0 , ξ 1 , . . .) and x 1 = φ(x 0 , ξ 1 ) with φ the t time step integrator. Further, the product can only be nonzero if both h A (x) and h B (x 1 ) are equal to 1, and h B ( t ) can be replaced by h B (φ(x 0 , ξ 1 )) in the product, (x 0 , ξ 1 )).
The ensemble average in Eq. (6) can be written as 033068-20 where noise history and noise future appear within one equation. The limit t → 0 only exists formally since t will be taken equal to the typical MD step in any practical case.

APPENDIX B: THEORETICAL PERMEABILITY FOR 1D TOY SYSTEM
In the 1D toy system, we can express the permeability in an analytical shape for the following three situations: for a Brownian particle (based on the Smoluchowski equation), for a Langevin particle crossing a high barrier (based on the Kramers equation), and for a deterministic Newtonian particle.
The Smoluchowski equation for a Brownian particle leads to the permeability expression in Eq. (3). In the onedimensional case, the absence of any other degrees of freedom implies that the free energy F (z) and the potential energy V (z) are the same. Hence, for overdamped dynamics, we can derive a theoretical expression for P by inserting V (z) of Eq. (39) as F (z) into Eq. (3), Here, h = 2a and I 0 (x) = (1/π ) π 0 exp(cos θ )dθ is the 0th order modified Bessel function of the first kind. When V 0 = 0, then I 0 (0) = 1, and the resulting flat potential yields P = D/h. For large V 0 , the cosine barrier can be approximated by a second order Taylor expansion about z = 0, F (z) = V (z) ≈ V 0 (1 − ( πz h ) 2 ) and it can be assumed that exp(−βV (z)) rapidly decays when moving away from the membrane. This can be inserted into Eq. (3) and the integration boundaries can be moved from ±h/2 to ±∞. Solving the resulting Gaussian integral yields an approximation of P for large V 0 , An alternative approach to the Smoluchowski approach is to use Kramer's relation for the rate constant k instead. The permeability P is then obtained via Eq. (18) by first computing the rate k while assuming a hard wall at z = −W that can be taken to infinite. Using a harmonic approximation and a high barrier assumption, this rate constant k can be written as [74] where κ is the transmission coefficient that can be approximated using Kramers' relation Here, ω + is the frequency associated to the curvature at the top of the barrier: ω + = √ k + /m with V (z) ≈ V 0 − 1 2 k + z 2 . From the above Taylor expansion, we have k + = 2V 0 π 2 /h 2 and ω + = (π/h) √ 2V 0 /m. where we assumed that overall state A condition is statistically equivalent to the condition z < 0 for this case, which is a valid assumption for a high barrier. Inserting Eq. (B4) in Eq. (B3) and inserting Eqs. (B3) and (B5) in Eq. (18) gives the permeability for a high barrier, In the high friction limit where γ ω + , Eq. (B4) reduces to κ = ω + /γ , and P for large V 0 becomes Since D = k B T /(mγ ), this equation is equal to the Smoluchowski equation Eq. (B2) within the harmonic approximation for the high barrier. A Langevin particle with high friction is indeed well described by the overdamped dynamics of a Brownian particle.
In the low friction limit γ ω + , Eq. (B4) reduces to κ = 1, and P in Eq. (B6) becomes, for any V 0 , This friction-less limit is exactly the permeability of the deterministic particle. It can also be obtained from Eq. (24). Here, ξ = 1/2, since a deterministic particle in a flat free energy region either moves to the right (velocity positive), either to the left (velocity negative), which have equal Boltzmann probability. A particle moving to the right will reach the barrier top with a probability exp(−βV 0 ) and it will not recross, and therefore P Let us recap the case of a flat potential membrane (V 0 = 0). For the Brownian particle, the permeability is P = D/h. If the particle has low friction γ → 0, then D → ∞, and the permeability diverges. For the Langevin particle, the high friction limit of P in Eq. (B7) based on Kramer's relation vanishes when V 0 = 0, which is not an adequate approximation of a flat potential's permeability. Nevertheless, again considering a Langevin particle and Kramer's equation, the low friction limit in Eq. (B8) converges to P = √ k B T /(2π m), which is finite.
In conclusion, the two theoretical expressions for the one-dimensional case, Eq. (B1) based Smoluchowski and Eq. (B6) based on Kramers, use respectively an overdamped assumption or a harmonic approximation to describe the top of the barrier. For high friction and low barriers, Eq. (B1) will be more accurate than Eq. (B6). For high barriers and low friction Eq. (B6) will prevail over Eq. (B1). In the case that both the friction and the barrier is high, both converge to the same value. In the case that both the friction and the barrier is low, neither Eq. (B1) nor Eq. (B6) will be accurate descriptions of a Langevin particle.