Diffusion of molecules through nanopores under confinement: Time-scale bridging and crowding effects via Markov state model

Diffusion model

In many cases, the transport of molecules though nanopores can be adequately described within the diffusion approximation by considering the molecules as over-damped Brownian particles [19,26,27], and assuming the adiabatic separation of the time scales, i.e., the fast degrees of freedom are in the thermodynamic equilibrium, while few slow ones are evolving quasi-stationary according to a stochastic Markov process. Often, it is sufficient to consider only one slow coordinate, e.g., the one characterizing the position of the molecule along the pore axis or the reaction path. The corresponding one-dimensional Fokker–Planck equation (FPE) for the probability density of molecule’s coordinate, ϕ ( x , t ) , is reduced to the Smoluchowski diffusion-drift equation,

where k stands for Boltzmann’s constant, T is the temperature, and D ( x ) is the position-dependent diffusion coefficient. U ( x ) is called the potential of mean force (PMF) as it generates the mean drift force acting onto the Brownian particle, f ( x ) = − d U ( x ) / d x . The difference between the PMF values at two points, x and x = 0 , equals the minimal thermodynamic (reversible) work required for moving the molecule from point x = 0 to x , U ( x ) = ∫ 0 x f ( x ′ ) d x ′ = W min ( x ) ; here, we set U ( 0 ) = 0 . According to the fundamental principles of thermodynamics [28], the minimal work required to change a parameter of a system while keeping it in the thermodynamic equilibrium, e.g., maintaining constant pressure, p , and temperature, T , is equal to the corresponding difference of the appropriate thermodynamic potential of the total system (particle plus medium), e.g., the Gibbs free energy W min ( x ) = G ( x , p , T ) − G ( x = 0 , p , T ) ≡ Δ G ( x ) ; therefore, U ( x ) = Δ G ( x ) . Equation (1) has a Boltzmann-type equilibrium (zero-flux) solution, determined by the PMF and independent on the diffusion coefficient,

where c 0 is a constant dependent on the boundary conditions (the concentration of the molecules outside the pore).

As far as the PMF, U ( x ) , is used to characterize the distribution of the particles in the thermodynamic equilibrium, the timescales of the fluctuations of the variables are not important. The description is equivalent to the classical Gibbs distribution. When one uses the PMF, U ( x ) , in the diffusion-drift equation to treat the system out of the thermodynamic equilibrium, the adiabatic separation of x from all other degrees of freedom must be assumed, in general. Fortunately, this very strong condition may be softened and the validity of the Smoluchowski equation may be extended if one introduces an appropriate modified position-dependent diffusion coefficient, D ( x ) [29,30, 31,32]. Still, it is worth to remind that diffusion is a time-average concept with respect to the all-atom dynamics. It assumes Maxwell’s distribution in the momentum space, so that the diffusion time scale is much larger than the corresponding momentum relaxation time of the particle, t ≫ τ p = D M / k T , where D is the diffusion constant and M is the mass. For a typical biologically related small molecule ( M = 250 Da , D = 1 nm 2 / ns ), at room temperature, one obtains τ p = 0.1 ps .

In the practically important case of the steady-state, the particle distribution is time independent, ϕ ( x , t ) = ϕ ( x ) , but the diffusion current, I , is nonzero and depends on boundary conditions. In particular, if the molecular bulk concentration on one side of the pore (cis side, at x = 0 ) is c C and on the other side of the pore (trans side, at x = L ) is c T and the diffusion current is determined by the Kramers-type formula [26]:

where S is the cross section of the channel at the entrance of the pore (the active surface of the pore). The average number of (noninteracting) diffusing molecules in the pore reads:

For small numbers, N c ≪ 1 (at low concentrations c C ), this quantity has the meaning of the probability to find a single molecule in the channel (channel occupancy), q = N c .

For the specific cases, equation (3) has been known and practically used for many decades, e.g., the Fick’s law for the diffusive current through membranes driven by the gradient of the concentration, the Goldman–Hodgkin–Katz equation for ionic current through membrane driven by the gradient of the electric potential and ion concentration difference, etc. The formula for the diffusion current is modified in the case when the radiation boundary conditions are more appropriate [33], or for a single particle in a periodic PMF driven by a constant force [19,34,35].

Equation (3) demonstrates that the PMF enters the current in the exponent, while the diffusion coefficient does only linearly: the dependence of the diffusion current on the PMF is stronger. In many cases, the diffusion coefficient can be effectively taken as a constant, D ( x ) ≈ D ∗ . If the PMF creates a large barrier in the constriction region, U ( x c ) ≫ k T , the integral in (3) may be restricted to a short interval, Δ x , around the barrier position, and the diffusion current reads:

This latter concept has been successfully applied recently to model the permeation of antibiotics through outer membrane channels of Gram-negative bacteria [7,36].

The central input quantities of the model, U ( x ) and D ( x ) , may be determined microscopically by using individual trajectories from all-atom simulations including, e.g., the pore, the membrane, the solvent, and the molecules, see, e.g., [37,38, 39,40]. For large molecules and strong interaction, the timescale for the particle to cross the pore increases and plain all-atom simulations may become unable to achieve statistically convergent results. Then, the enhanced sampling methods (e.g., the umbrella sampling [41] or the metadynamics [42,43]) may be utilized. But even when using the enhanced sampling, the all-atom approach is still often computationally demanding as multi-microseconds long trajectories may be required [23,44, 45,46].

Markov state model

The Markov state kinetic model of the pore–molecule interaction, used, e.g., for the analysis of the channel gating in the electophysiology [47], starts with a set of states defined by different ionic conductance levels of the channel (the “bound states”) due to the presence of substrate molecules. The transitions between these states are assumed to be much faster than the corresponding residence times. The ground state is the open-channel state (optionally, several) without any molecule inside. The transitions between the states are described as a discrete Markov process determined by the corresponding rate constants – the target parameters describing the pore–molecule interaction kinetics [22].

Here, we consider a one-binding-site (two-states) Markov model of the channel, the extension to the many-state model is straightforward.

Let p o , p b are the probabilities to have the open channel and a single molecule in the pore, correspondingly. In the Markov state model, one introduces the transition rates, λ o b , λ b o , between the two states of the pore occupation states, and the simultaneous rate equations:

Additional normalization condition, p o + p b = 1 , reflects the physical assumption that only a single molecule a time can stay in the pore. Then, one introduces the kinetic rates, k off C , k off T , for the transition of the particle from the binding site to the “cis” and the “trans” sides of the channel, respectively, so that λ b o = k off C + k off T . Also, one assumes that the “cis” and “trans” compartments are large reservoirs keeping the molar concentrations of the molecules at the sides of the channel invariant, respectively, c C and c T . Thus, one can write down, c C k on C and c T k on T , for the transition of the particle from the “cis” side and from the “trans” side, respectively, onto the binding site in the channel, so that λ o b = c C k on C + c T k on T . Besides, we define the following quantities, which can be determined in the electophysiology:

Below, we summarize necessary results from the two-state Markov model in the steady-state [22,47,48]. If the particles are added on the both sides at the concentrations, respectively, c T and c C , the steady-state probability to find a molecule in the channel reads,

The net current of molecules through the pore reads,

With the aforementioned definition, the current is positive for the net transport from cis to trans, and it is negative for the opposite direction.

From the aforementioned results, and from the general requirements of the zero net particle current, in the case of the symmetric bulk solution on the sides of the channel (equal concentrations, no electric field gradient, etc.), I C → T ( c , c ) = 0 , one arrives at the detailed balance condition for the transition rates,

Thus, one finally obtains,

Equation (13) shows that the molecular pore-translocation current obeys Fick’s law for small concentrations only, c T K T + c C K C ≪ 1 ,

In the opposite case of large concentrations, in contrast with the diffusion approximation, the molecular current saturates and reaches its maximum, I max .

Another important difference between the Markov state model and the diffusion model is that the diffusion molecular current is always symmetric with respect to the side of application of the concentration gradient, while the Markov state translocation current is symmetric only at small concentrations. In general, the translocation current depends on the side of application of the concentration gradient,

In the case of asymmetric solvent on cis and trans sides, due, e.g., to the different solvent activity or an applied electric potential for charged molecules, the definition of the in-rates should be modified, λ o b = c C exp ψ C k T k on C + c T exp ψ T k T k on T . The constant ψ = μ − k T ln ( x ) is the difference between the chemical potential of the particle and the logarithm of its molar fraction; it contains both the activity and the applied electric potential. The detailed balance condition, (12), remains, but the net translocation current becomes ψ -dependent (compared with equation (3) for 1D diffusion),

Interestingly, the maximum current at high concentration is given by the same equations (16).

Bridging 1D diffusion and the Markov state models

To connect the kinetic rates of the Markov state model with the PMF and the diffusion coefficient, one can consider the low concentration limit, where the both models are linear with the concentration. Similar arguments have been used by other authors, see, e.g., [24,25]. We also assume that the solute at concentration c is added from one side only, and the PMF has same value at the ends, U ( 0 ) = U ( L ) = 0 . The linear-regime probability to find a particle in the channel when added only from the trans and only from the cis side, respectively, can be defined as follows:

The linear-regime particle current through the channel when added to one side (the same for trans and cis) at concentration gradient c , reads,

Quantities q T ( c ) , q C ( c ) and I D ( c ) can be calculated within the diffusion model by using equations (3) and (4). Thus, by defining the following quantities:

satisfying R Q = T 0 + T 1 , one finds,

The aforementioned formulas connect the Markov state model with the diffusion approximation. T 0 , T 1 quantities are equal to the mean first passage times defined in refs [24,25].

The translocation current and the channel occupancy at higher concentrations can be calculated by using Markov state model results, equations (10) and (13). In particular, if the molecules are added from one side only, one finds,

It is worth to note again that the two-state Markov state is a time average concept with respect to the continuous diffusion approximation. It assumes that the average time spent by the molecule in each of the states is much larger than the average diffusion time required to cross the barrier between the states. The latter can be estimated as τ D = Δ l 2 / D , where Δ l is the characteristic width of the barrier. Therefore, the conditions of the applicability of the 1DD-2MS bridging discussed earlier read,

Only if the aforementioned conditions are satisfied, one can use the obtained kinetic rates ( k on , k off ) to trace the nonequilibrium kinetics of the molecule–pore interaction, e.g., the opening and closing time series at it is observed by the electrophysiology.

However, equations (31) for the steady-state tanslocation current are also valid at much weaker conditions. It assumes two main approximations. First, the current, e.g., from cis to trans, is presentable in the form (compared with equation (20)),

where p 1 C is the steady-state probability to find a single particle in the channel when added from cis side, and λ off T is the rate of exit from the pore to trans side. This representation is not limited to a 2MS. It is also valid, for example, for any number of Markov states inside the pore, with only one restriction – the exit rates from all the Markov states to trans side equal either zero or λ off T . The second condition of equation (31) is that not more than one particle at a time can occupy the pore, again with no specific restriction to the 2MS model. If the solute molecules are noninteracting in the solution (ideal solution), the probability to find n molecules in the channel can be given by the Poisson distribution:

where A is the normalization constant, ∑ n p n = 1 , and q is a parameter. If the molecules do not interact also in the pore, n runs from 0 to ∞ , A = 1 and q is the average number of molecules in the pore. In the case of no more than one molecule in the pore, then n = 0 , 1 and A = e q / ( 1 + q ) , and one arrives at equation (30) for p 1 C and equation (31) for the translocation current. Methodologically, this approximation is similar to the introduction of the exclusion volume into the equation of state of ideal gas. It explains two main qualitative effects of the interaction of the molecules in the pore – the saturation of the translocation current and its asymmetry with respect to the side of addition, both can be observed and studied.

By following the latter approach, one can formally extend the model of at maximum N noninteracting molecules in the pore,

Analytic single barrier/binding site model

The barrier (binding site) has the width Δ l and is positioned at a distance x 0 from the left entry of the pore (cis compartment) and x 1 from the right entry side of the pore (trans compartment). Let us suggest that the diffusion coefficient, D , is constant along the channel and the PMF is a rectangular function of magnitude Δ U and width Δ l , located at x = x 0 ,

Here, rect ( x ) = θ ( x ) − θ ( x − 1 ) is the rectangular functions on interval [ 0 , 1 ] ; θ ( x ) is the Heaviside step function.

The linear diffusion current due to the gradient of concentration c then reads

The molecular current is independent of the rectangle position, x 0 . For the large barrier, Δ U ≫ k T ln ( L / Δ l ) , I D decrease exponentially with its magnitude, while for the deep well (binding site), Δ U ≪ − k T , the diffusion current becomes independent on its magnitude, respectively,

Thus, the presence of a barrier for the diffusing molecule decreases the diffusion current, while the binding site increases the diffusion current effectively by reducing the channel length to L − Δ l .

The channel occupancy for the cis molecule addition, for the trans addition, and for the symmetric addition reads, respectively:

Here, the distance to the barrier/well from the trans side, x 1 ≡ L − x 0 − Δ l is introduced for increasing symmetricity of the formulas.

1DD-2MS mapping, large barrier:

Finally, the translocation current through the high barrier reads for the cis and the trans addition, respectively,

At high concentrations, c ≫ 1 / ( S x 0 ) and c ≫ 1 / ( S x 1 ) , respectively, the currents reach their maximum values,

The saturated translocation current is larger when molecules are added from the side to which the barrier is closer (the barrier at the entry side).

The mapping onto the 2MS model is formal in this case, as molecules stay either on the left side of the barrier or on the right one (two different states). Thus, the obtained transition rate constants do not have real physical meaning. Nevertheless, the steady-state current, equations (46) and (47), can be used for the calculations as far as the assumption of a single (at maximum) molecule in the pore holds.

1DD-2MS mapping, strong binding:

The translocation current through a deep potential well reads, respectively,

At the high concentrations, the currents reach their maximum values,

In contrast with the barrier, the saturated translocation current through the well is larger when the molecules are added from the side opposite to the well (the binding site next to the exit).

In the case of strong binding site nonadjacent to the channel ends, i.e., x 0 , 1 ≠ 0 and U ≪ k T ln ( 2 Δ l x 0 , 1 / ( x 0 + x 1 ) 2 ) , the 2MS model parameters,

have physical meaning. Indeed, in this case, the binding site is single and the 1DD-2MS bridging conditions, equations (32)–(34), can be satisfied for a range of concentrations from zero till the saturation.

Analytic one-barrier-one-binding site model

Let us suggest that the PMF is a sum of nonoverlapping rectangular barrier of magnitude Δ U + and width Δ l + , located at x = x + , and rectangular well of magnitude − Δ U − and width Δ l − located at x = x − (here, the quantities, Δ U ± are positive numbers, see Figure 1c),

The diffusion coefficient, D , is constant along the channel. The diffusion current due to the gradient of concentration c then reads

The channel occupancy for the cis molecule addition, for the trans addition and for the symmetric addition, reads, respectively,

Here, the following “distances” to the barrier/well from the trans and cis sides are introduced for more symmetric formulas, x 0 + = x + ; x 1 + = L − x + − Δ l + − Δ l − ; x 0 − = x − − Δ l + ; x 1 − = L − x − − Δ l − As we can see in .

Figure 1

Schematic diagrams of a pore with (a) a single barrier PMF, (b) a single binding site PMF, (c) a barrier and a binding site PMF, and (d) a two-barriers-one-binding-site PMF. The dotted lines represent the pore length.

Now, we will assume that the barrier and the binding site are strong, i.e., Δ U ± ≫ k T ln ( L / Δ l ± ) . Then, one obtains,

That is, in the presence of a strong barrier, the diffusion current is independent on the binding site. The translocation current reads,

At high concentrations, currents reach their maximum values,

As, I C → T > I T → C , we conclude that by combining a barrier with a well, the higher saturated translocation current is obtained when the barrier is closer to the entry side (where the concentration is higher) and the well on the exit side, following the addition rule of the single barrier and single-well behavior.

Markov chain model

The elementary cases discussed in the two previous subsections allows one to write down the transition rates for a more general situation of a PMF presented as a sequence of binding sites, separated by large rectangular barriers. The molecule moves in the pore by jumping between neighboring binding sites. We number the binding sites from 1 (the first site on the cis site) to N (the last one on the trans side), so that U k − , Δ l k − are the corresponding energy and the width, respectively. The barrier to the right of the binding site k has energy U k + and width Δ l k + , so that the first barrier (to enter the pore from cis) has number 0 and the last one (to exit to trans) has number N . Then, the transition rates from site k to sites k ± 1 read, respectively,

Here, D k is the diffusion constant at barrier k ; λ N N + 1 = k off T , λ 10 = k off C . Also,

This mapping of the diffusion onto the Markov chain model assumes that the barriers are large enough for the Markov timescale conditions, like those given by equations (33)–(32), being satisfied.

By using the transition rates obtained by the Markov chain mapping, one can build a kinetic Monte-Carlo scheme [49] to simulate the translocation kinetics. This approach makes it possible to limit the occupancy for each binding site (not the whole pore), thus enhancing the possibility to study molecular crowding in the channel.

Equations (63) and (64) can be interpreted as a diffusion analog of the classical dynamic transition-state theory [26], where D / ( Δ l + Δ l − ) ratio replaces the frequency prefactor.

Analytic two-barriers-one-binding site model

A more realistic example is provided by the two-barrier-one-binding-site model, introduced earlier to model the transport of metabolites with binding in specific porins [22,48].

Here, Markov chain mapping formulas, (63)–(66), give (the diffusion coefficient is assumed constant along the pore for simplicity),

Then, one obtains,

The translocation current is given by equations (31). The maximum cis/trans asymmetry,

is determined by the magnitude and the width of the barriers.

This model contains two previous examples as special cases: (i) a symmetric case, where the barriers from the two sides are the same; (ii) a strongly asymmetric case, where the barriers are different.

Thus, in the asymmetric case, ( U 1 + − U 0 + ) ≫ k T , the results and consequences are not different from the case of a single barrier with a binding site. When the barrier heights differs less than k T , we can consider the symmetric case.

Numerical validation

In order to show the potency of the method introduced, we considered a realistic model system, the transport of a beta-lactamase inhibitor (tazobactam) through the main porins expressed in the outer membrane of Escherichia coli, OmpF and OmpC. Molecules as tazobactam are used in combination with beta-lactam antibiotics because they bind strongly to beta-lactamase enzyme inhibiting them from modifying chemically antibiotics and thus revert their antibacterial efficacy [50,51]. Today, there is a wide interest in developing such inhibitors to counteract the spreading of resistant pathogens, and one common problem is to predict their transport through porins. By using multiple walkers, metadynamics simulations we reconstructed the one-dimensional free energy surface for the interaction of tazobactam with the two porins, as done for other inhibitors recently [52]. As we can see in Figure 2, the two FES profiles are a complex combination of local minima and barriers. In particular, we note a main central barrier and some minima on the positive region z (extracellular region), thus reminding the simple model of a barrier and a binding site. By calculating the diffusion current with equation (3) and thus applying equation (31), we obtained the total currents as function of gradient concentration, as reported in Figure 3.

Figure 2

Free energy surface profiles vs tazobactam position along the diffusion axis inside the OmpF and OmpC porins.

Figure 3

Current of tazobactam through the OmpF and OmpC porins versus the concentration of the molecule being added to one side only: Cis to Trans (left panel) and Trans to Cis (right panel) currents.

A comparison of the two currents at constant concentration is not sufficient, as, for example, for the cis to trans current in the left panel of Figure 3. Though the barrier is lower in OmpF, than OmpC, and correspondingly larger the current, in the case of OmpF the saturation occurs before, and at concentrations higher than 10 mM, the current through OmpF becomes lower.

It is also interesting to calculate the rate of the cis → trans and trans → cis currents, as provided in Figure 4. As we illustrated with in the aforementioned simple examples, the two current are symmetric at low concentration, while at high concentration, there is an asymmetry modulated by the position of minima and barrier, see equations (57), (58), and (31):

• The current trans → cis is always superior to the cis →trans, because the highest barrier is on the trans side.

• The presence of minima is responsible for the saturation of the current at high concentration, and the change in current is more drastic on the cis → trans current because the minima are on the cis side.

Figure 4

Ratio between the Cis to Trans and Trans to Cis current of tazobactam (being added to one side, Cis or Trans, respectively) through the OmpF and OmpC porins.

Experimental validation

Our model can be employed to guide and plan experiments till the interpretation of results. A quick look at recent published data shows that kanamycin, a positively charged molecule, is transported through both OmpF/OmpC [53]. At high negative external potential applied, the association rate in OmpF saturates, contrarily to OmpC. Because kanamycin is a molecule charged positive, a high negative potential applied means to force a large association rate toward the pore, corresponding to the case described here as an increase in the gradient concentration. Interestingly, the obtained results are described by the top panel of Figure 3. Increasing more and more the external potential we would see higher association rate through OmpC than OmpF.

Further, the behavior described in Figure 4, an asymmetric transport upon asymmetric addition of molecules in solution, was measured for maltoheaxose through the maltoporin trimeric channel [21]. This sugar molecules has higher association rate when added on the trans side, exactly as for tazobactam through OmpF/OmpC. In that case, the asymmetry is visible already at 0.01 mM concentration, because of the high specificity of the maltoporin channel for sugar molecules.