We present a statistical model for solving and predicting the transport of large molecules through small flexible channels. The average radius of the channel and the average radius of the molecule are the only two quantities determining the steric part of the potential of mean force for the translocation, in the case of a small rigid particle and a large rigid channel: the barrier is completely entropic and is described by the Fick-Jacobs model. However, the flexibility of the channel’s cross section and that of the molecule’s size have a significant effect on transport, especially when a large molecule goes through a narrow channel. In this case, the steric barrier changes its statistical nature becoming enthalpic, and we predict a strong temperature enhancement of the diffusion current through the channel. The flexibility is described in terms of the equilibrium fluctuations of the channel and of the molecule. The model is compared with the all-atom MD simulations of the transport of hard spheres of various radii and of drug molecules through a biological nanochannel. For the case of Gaussian fluctuations, we derived a simple analytical expression for the steric barrier, which can be quantified using average size and fluctuations of the channel and of the molecule.

Passive transport of molecules through nanosized channels in porous media is of fundamental importance for applications in materials science,^{1} nanotechnologies,^{2–4} and biology.^{5–10} The main key features of the nanoworld and then of biological systems are random fluctuations, dictating flexibility and diffusion. Thus, understanding the passive transport of flexible molecules through flexible channels is a priority for the development of biomimetic applications and in nanotechnologies.

In the present work, we consider the diffusion of a flexible particle through a channel with a narrow and flexible constriction region so that the average equilibrium size of the molecule may be larger than the average size of the pore. This is a typical situation in biology, where the opening of a large pore in membranes contrasts with maintaining the specificity of compartments (membrane potential and concentration of particular substrates). But since the pore can be expanded and the molecule can be compressed, the latter still can be translocated over the constriction barrier. The passive transport may often be adequately described within the diffusion approximation by considering the molecules as overdamped Brownian particles.^{11–13} To our knowledge, there are no existing implicit models to tackle the problem of flexibility without performing all-atom simulations to determine the central input quantities of the transport process, the potential of mean force (PMF), U(x), and the position-dependent diffusion coefficient, D(x), x being the position along the pore axis. Even using enhanced sampling methods (e.g., the umbrella sampling^{14} or the metadynamics^{15}), the simulations of the whole system—the particle, the channel, and the medium—are often computationally demanding as multimicrosecond long trajectories may be required.^{16–18} Our aim is to connect analytically the PMF of the particle in the pore with its geometric and flexibility parameters. Here, we limit our analysis to the steric contribution to the molecule-pore interaction free energy which originates from the fact that atoms have finite size and cannot overlap. As we will see, the steric effect has significant contribution to the PMF, creating a barrier for the diffusing particles.

We consider a channel having a circular cross section and a spheroidal particle with its axis always aligned along the channel (Fig. 1). Furthermore, we select a minimal set of variables characterizing the configuration of the particle in the channel—{*x*, *y*, *z*} are the Cartesian coordinates of the center of the particle; *r*_{c} is the radius of the channel at the position of the particle; and *r*_{m} is the radius of particle’s spheroid orthogonal to the axis. The minimum work, *W*_{min} = *W*_{min}(*x*, *y*, *z*, *r*_{c}, *r*_{m}), required to move the particle from a certain point {*x*_{0}, *y*_{0}, *z*_{0}} (e.g., *x*_{0} = *y*_{0} = *z*_{0} = 0) and equilibrium values of pore’s/particle’s radii, to point {*x*, *y*, *z*} with radii *r*_{c}, *r*_{m} determines the equilibrium probability density for the parameters,^{19} $\psi (x,y,z,rc,rm)\u221dexp\u2212Wmin/kT$; here *kT*, stands for the thermal energy. Then, the equilibrium probability density for *x* reads

On the other hand, *ϕ*_{eq}(*x*) may be obtained as the equilibrium (zero-flux) solution to the Smoluchowski diffusion-drift equation, determined by the PMF and independent of the diffusion coefficient

here, *c*_{0} is a constant dependent on the boundary conditions. Therefore, one may introduce a physical model for *W*_{min} and then use Eqs. (1) and (2) to calculate the PMF.

By assuming constant pressure and temperature, *W*_{min} may be formally calculated as the free energy change along the thermodynamic cycle,

The first term stands for the contribution due to the change in the free pore’s radius at point *x* from its equilibrium value to *r*_{c}; it determines the equilibrium probability density of the pore radius fluctuations, $fc(rc,x)\u221dexp\u2212\Delta Gc/kT$. The second term is the contribution from the change in the particle’s radius from its equilibrium value to *r*_{m}, outside the pore, and its probability density reads $fm(rm)\u221dexp\u2212\Delta Gm/kT$. The last term in (3) contains all the effects of the particle-pore interaction, and it vanishes by definition when particle’s coordinate *x* is outside the pore.

By assuming the hard-wall repulsion between the channel and the particle, we define the steric interaction potential as follows: it equals zero when the channel and the particle do not overlap; it is positively infinite otherwise. For the coordinate system chosen in Fig. 1, one obtains

has the meaning of the average available cross-sectional area for the particle inside the pore at coordinate *x*. Then, using Eq. (2), the PMF may be calculated as

Equations (5) and (6) give the steric contribution to the potential of mean force of a molecule diffusing in the channel taking into account both their average size and the flexibility. This is the main result of this work. The input quantities, i.e., probability densities *f*_{c}(*r*_{c}, *x*) and *f*_{m}(*r*_{m}), may be calculated by studying the equilibrium fluctuations of *r*_{c}(*x*) and of *r*_{m} in all-atom MD simulations, separately, for the channel and for the molecule. As the fluctuations are not necessarily small, advanced sampling techniques may be used to accelerate the statistical convergence.

The major assumption behind the present model consists in using the hard-wall repulsion approximation (4) for the interaction free energy. Therefore, the effect of the molecule/pore steric interaction consists in the limitation of the configuration space. The selection of the collective variables characterizing channel/particle deformations is another important approximation of the model. We have used the minimal set of variables in our example, thus neglecting, e.g., the free energy related to the possible rotation of the long axis of particle’s spheroid with respect to the axis of diffusion. 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 (including those we used to calculate *U*) 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 position-dependent diffusion coefficient, *D*(*x*),^{20–23}

The diffusion coefficient may depend on the size of the channel “directly,” e.g., due to the limitation of the configuration space for the coordinates of the diffusing particle,^{20–22} and also “indirectly,” e.g., via the modified hydrodynamic interactions in the confined solvent which affect particle’s mobility.^{24–26}

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, when the diffusion is driven by the concentration difference, Δ*c*, at the boundaries of the channel of length *L*, the diffusion current is determined by the Kramers-type formula^{11,20}

The formula for the diffusion current is modified in the case when the radiation boundary conditions are more appropriate,^{27} or for a single particle in a periodic PMF driven by a constant force.^{12,28,29} Equation (8) demonstrates that the PMF enters the flux in the exponent, while the diffusion coefficient only linearly: the dependence of the diffusion current on the PMF is stronger. In many cases, the diffusion coefficient may be probably taken as an effective constant, *D*(*x*) ≈ *D*^{*}. If the PMF creates a large barrier in the constriction region, *U*(*x*_{c}) ≫ *kT*, the integral in (8) may be restricted to a short interval, Δ*x*, around the barrier position, and the diffusion current reads *I* = Δ*cA*(0)*D*^{*} exp(−*U*(*x*_{c})/*kT*)/Δ*x*.^{30}

To illustrate the results, we introduce the two corresponding average values *R*_{c}(*x*) and *R*_{m}, the corresponding standard deviations *σ*_{c}(*x*) and *σ*_{m}, and the difference in radii *R*_{cm}(*x*) ≡ *R*_{c}(*x*) − *R*_{m} with its standard deviation $\sigma cm2(x)\u2261\sigma c2(x)+\sigma m2$. If the channel is wider than the molecule, *R*_{cm}(*x*) > 0, and when both are stiff, *R*_{cm}(*x*)/*σ*_{cm}(*x*) ≫ 1, then distributions *f*_{c}(*r*_{c}, *x*) and *f*_{m}(*r*_{m}) may be approximated by the Dirac delta functions and one obtains $A(x)\u2248\pi Rcm2(x)$. The PMF profile, $U(x)\u2212U(0)\u2248\u22122kT\u2061lnRcm(x)/Rcm(0)$, is determined by the difference in the pore and the particle equilibrium radii, in agreement with the Fick-Jacobs approximation.^{20–22,29,31}

Now, we suppose that distributions *f*_{c}(*r*_{c}, *x*) and *f*_{m}(*r*_{m}) are normal, and the average radii are significantly larger than the corresponding standard deviations. Then, distribution *f*_{cm}(*t*, *x*) is also a Gaussian having mean value *R*_{cm}(*x*) and variance $\sigma cm2(x)$. If the particle barely fits the pore at the constriction region so that |*R*_{cm}(*x*)|/*σ*_{cm}(*x*) ≪ 1, then one finds $A(x)\u2248\pi \sigma cm2(x)/2$, i.e., the PMF is completely determined by the pore/particle flexibility in terms of equilibrium fluctuations. In the limit of a narrow channel and a thick molecule, *R*_{cm}(*x*) < 0 and |*R*_{cm}(*x*)|/*σ*_{cm}(*x*) ≫ 1, one may limit the integration in (5) to the range from 0 to $\u2248\sigma cm2(x)/|Rcm(x)|$ and arrive at the following estimate:

Therefore, the average available area is exponentially small in this case and depends on both the average radii and fluctuations. By inserting (9) into (6) and neglecting the logarithmic terms of |*R*_{cm}(*x*)|/*σ*_{cm} compared with the quadratic ones, one obtains the estimate of the steric barrier

We consider the effect of fluctuations of both the pore and the particle on the diffusion current [Eq. (8)]. We again assume Gaussian distributions for the average radii *R*_{m} and *R*_{c} (now, independent of *x*), respectively, with variance $\sigma m2$ and $\sigma c2$; *D* is constant all along the channel. From Eq. (8), one finds that the ratio of the diffusion current through this pore to the current of pointlike particles through a rigid channel of radius *R*_{c} reads $I/I0=A/(\pi Rc2)$, where $I0=\Delta cD\pi Rc2/L$ and *A* is given by Eq. (5). This dimensionless diffusion current is reported in Fig. 2 vs the channel/particle variance. One can see that while for small particles (*R*_{m}/*R*_{c} = 0.1), the diffusion current virtually does not depend on fluctuations, at *R*_{m}/*R*_{c} = 1, the current increases linearly with the variance, $\sigma cm2$. For the particles larger on average than the pore, *R*_{m}/*R*_{c} = 1.5, the current exponentially grows with $\sigma cm2$. Equation (10) also shows that, in contrast with the small-particle limit where the PMF profile has logarithmic dependence on *R*_{cm}(*x*), a much stronger power-law dependence on |*R*_{cm}(*x*)|/*σ*_{cm} holds in the opposite limit. One may also conclude that the PMF changes its statistical origin with the increase in the relative molecule/pore size from purely entropic in the Fick-Jacobs limit^{20} to essentially enthalpic in the opposite case of a large molecule in a small pore. Indeed, as *R*_{cm}(*x*) may be considered independent of or weakly dependent on the temperature, the PMF in the FJ limit is proportional to temperature. According to the fluctuation-dissipation theorem,^{19} the fluctuations of the molecular size are also proportional to *T*, $\sigma m2=\alpha kT$; here, *α* = |*dR*_{m}/*df*| determines the compressibility of the molecule under external force *f* in the limit *f* → 0. Analogously, we can conclude for the channel radius. If *α* does not depend (or weakly depends) on the temperature, the steric barrier (10) becomes temperature-independent [*S*(*x*) = −*∂U*/*∂T* ≈ 0], i.e., has enthalpic property. The latter assumption seems reasonable in the case of large deformations and for a stiff molecule/channel, when the deformations are close to elastic and the work is transformed into the internal potential energy of the molecule/channel. As an important exception, the PMF for translocation of a Gaussian polymer through a rigid pore, calculated within the continuous random walk approximation,^{32–34} is always entropic, as in the FJ case. The Gaussian polymer model, however, does not take into account the finite volume and the flexibility of the monomers, and therefore, it is not applicable when the monomer size is close to or larger than the pore. To quantify the temperature effect on the diffusion current, we define the dimensionless temperature sensitivity coefficient of the channel as *δI*/*I* = *η*_{T}*δT*/*T*. By using the arguments discussed above, one obtains $\eta T=\u2202A\u2202\sigma cm2\sigma cm2A$. On the inset graph in Fig. 2, we show the dependence of *η*_{T} vs the channel/particle fluctuations. While there is almost no temperature dependence (*η*_{T} ≪ 1) for small particles, one would expect a linear increase (*η*_{T} ≈ 1) in the diffusion current with *T* at *R*_{m}/*R*_{c} = 1 and strong temperature enhancement of the diffusion current for *R*_{m}/*R*_{c} > 1 [$\eta T\u224852+Rcm22\sigma cm2$, derived from Eq. (9)].

We apply the developed theory to calculate the PMF for translocation of particles through a biological nanopore, the water-filled *β*-barrel protein OmpF from *Escherichia coli*. This porin transports passively polar nutrients through the outer membrane and is considered the main gate for polar antibiotics. Therefore, it is used as a model system for transport in membrane biophysics.^{8,30} The channel is ≈60 Å long and has a quite rigid hourglass internal shape, with an average radius of ≈17 Å at the entrance and of ≈2.9 Å in the constriction region at *x* = −2 Å; see the inset of Fig. 3. The porin was inserted into a lipid bilayer membrane, solvated in a box of water (size 64 × 64 × 80 Å^{3}) and equilibrated at 300 K and at 1 bar pressure by using all-atom MD simulations, as described earlier.^{35} We calculated (the method is described in Ref. 30) the pore radius profile, *r*_{c}(*x*), and the probability density histogram, *f*_{c}(*r*_{c}, *x*), averaging every 50 ps over a 300 ns trajectory performed in the NVT ensemble, inset of Fig. 3.

By using the calculated *f*_{c}(*r*_{c}, *x*) for the channel and Eqs. (5) and (6), we computed the steric PMF for hard spheres of radii from 2.5 Å to 4 Å; the results are shown in Fig. 3. Besides, we calculated the total PMF by including the spheres into 5–10 *μ*s MD simulations and by using metadynamics to increase sampling; the hard spheres were modeled as van der Waals atoms with *R* = *R*_{m} and *ε* = 0.2 kcal/mol; the position of the center of mass along the channel axis, *x*, and its distance from the axis, $y2+z2$, where biased in metadynamics runs. We see that the steric effect dominates in the constriction region although the total PMF includes the cooperative effects of the solvent, e.g., the hydrophobic interaction. The model describes the PMF profile fairly well, including the height, position, and shape of the central barrier.

Furthermore, we calculated the steric contributions [using Eqs. (5) and (6); at T = 300 K] to PMF for two antibiotics, norfloxacin (NOR, zwitterionic at normal pH) and carbenicillin (BAR, charged −2*e* at normal pH) in OmpF, and compared those with the total PMFs obtained by using all-atom simulations with metadynamics; see Fig. 4. The details of the all-atom simulations for NOR were presented previously.^{35} The procedure for BAR was similar: we used 6 walkers for a total simulation time of 3.4 *μ*s. The distribution of minimal radius of the molecular spheroid, *f*_{m}(*r*), was obtained calculating the minimal projection area in a preliminary all-atom simulation of the molecule in a box of water, as described in Ref. 30. The data presented in the insets of Fig. 4 clearly show that, on average, both molecules are larger than the pore at the constriction region, and only the flexibility makes the translocation possible: we are in the regime described by the dotted curves of Fig. 2, or *R*_{m}/*R*_{c} > 1. Although BAR is notably larger in average than NOR, it is also more flexible. As a result, the steric barriers for the two molecules look similar. The model also reasonably agrees with the total PMF profile in the constriction region, demonstrating that the steric contribution dominates the other interactions (electrostatic, hydrophobic, etc.), similar to the case of the hard spheres. The electrostatic interactions are relative weak for different reasons. For NOR, the large dipole moment does not align to the transversal electric field of OmpF in the central region.^{35} On the other hand, for BAR, the compensation of the steric barrier due to the alignment of the small dipole moment is balanced by its charges (−2*e*) that interact unfavorably with the negative charges of the central region.

Although the model is presented here in a simplistic form (we assume one-dimensional diffusion, circular channel, spheroidal particle, etc.), the suggested way of reasoning allows one to make a straightforward extension to the case of several collective variables, *x* = {*x*_{1}, *x*_{2}, …}. The two antibiotics considered, as well as many other relevant molecules, have average size larger than the OmpF in the constriction region. The presented method can be applied to other problems of diffusive transport, apart from the discussed case when antibiotic molecules pass through small porins,^{30} e.g., transport of flexible solutes through smaller rigid pores^{36} and uptake of hydrophobic molecules into the cellular membrane via the spontaneous opening in the lateral part of a membrane protein.^{37} Eventually, the strong dependence of the steric barrier on the particle size and flexibility discovered here could potentially be employed in the design of new molecular filters.

The research leading to these results was conducted as part of the Translocation consortium (http://www.translocation.com) and received support from the IMI Joint Undertaking under Grant Agreement No. 115525, resources which are composed of financial contribution from the European Union’s seventh framework program (No. FP7/2007-2013) and EFPIA companies’ kind contribution. M.C. thanks the MIUR PRIN Project No. 2015795S5W/005.