Diffusion of finite-size particles in two-dimensional channels with random wall configurations

Abstract

Diffusion of chemicals or tracer molecules through complex systems containing irregularly shaped channels is important in many applications. Most theoretical studies based on the famed Fick–Jacobs equation focus on the idealised case of infinitely small particles and reflecting boundaries. In this study we use numerical simulations to consider the transport of finite-size particles through asymmetrical two-dimensional channels. Additionally, we examine transient binding of the molecules to the channel walls by applying sticky boundary conditions. We consider an ensemble of particles diffusing in independent channels, which are characterised by common structural parameters. We compare our results for the long-time effective diffusion coefficient with a recent theoretical formula obtained by Dagdug and Pineda [J. Chem. Phys., 2012, 137, 024107].

---

I. Introduction

Diffusive transport in structured environments is a ubiquitous feature relevant in a large variety of systems. Inter alia, these range from the dispersion of tracers in the permeable rock or loose materials making up groundwater aquifers,¹ over the diffusion of chemicals in ramified matrices such as porous glasses or zeolites,^2,3 up to the random motion of biomolecules and submicron objects in the crowded environment of living biological cells with their complex scaffolding made up of structural, semiflexible polymers.^4–6 A general theme connecting these systems is the geometrical constraints imposed on the motion of the tracer particles, which may lead to interesting entropic effects.⁷

Important model systems for confined geometries are two- and three-dimensional channels and pores. The quantitative mathematical description of diffusion in such objects with varying width has a long-standing history (see ref. 8 and references therein). In a seminal work Zwanzig derived a modified Fick–Jacobs equation which is at the basis of most subsequent quasi-one-dimensional descriptions,⁹

Here G(x,t) describes the local concentration of particles at position x and time t, w(x) the width of the channel at position x and most importantly D(x) is an effective position-dependent diffusion coefficient. Subsequently, several different forms for D(x) were derived and applied to various systems.^9–19 Taking along only first order derivatives of the width profile w(x), Kalinay and Percus¹³ (see also Martens et al.¹⁶) obtained the closed form resultwhere D₀ denotes the position-independent diffusion coefficient in the absence of confinement.

However, the above form of D(x) is restricted to symmetric channels, i.e., channels with a straight centre-line. This constraint was removed in an approach proposed by Bradley,²⁰ which was subsequently generalised by Dagdug and Pineda to²¹

where y₀(x) denotes the vertical position of the centre-line at horizontal position x. Note that eqn (3) generalises all the previous results for D(x), for instance, one obtains the result (2) for symmetrical channels by setting y₀(x) = 0.²¹

As detailed by Zwanzig,⁹ in a system with periodic boundary conditions the effective diffusion coefficient in the long-time regime, D_eff, is obtained by using the following formula introduced by Lifson and Jackson²² and generalised by Festa and Galleani d'Agliano,²³

where 〈·〉 denotes the average over one period. Thus, for any channel with width profile w(x) and centre-line y₀(x), eqn (3) can be used to calculate D_DP(x), which in turn is used in eqn (4) to calculate the effective diffusion coefficient D_eff,DP.

Here we consider two generalisations of the above scenario. The first concerns the fact that the above theories based on the Fick–Jacobs formalism pertain to point-like particles, that is, the particle can move along the channel as long as the width is not exactly equal to zero. For scenarios, in which the size of the diffusing particle is comparable to the channel width, the particle can only fully cross the channel when the width profile w(x) at any position is larger than the particle size. This is particularly relevant for biological systems, for instance, in membrane channels.²⁴ Systems containing finite-size particles were indeed studied in the literature,^25–27 showing that, essentially, in all formulas the channel width w(x) has to be replaced by an effective channel width.

The second modification with respect to the Fick–Jacobs approach that we consider here concerns the boundary conditions. Usually, reflecting boundary conditions are used, i.e., when the particle collides with the channel walls, its perpendicular motion is simply reversed. In the present scenario, the diffusing particle is allowed to (transiently) bind to the channel walls. Such a behaviour occurs when the diffusing particle has a chemical or physical binding of adhesion propensity to the surface. We explicitly include this effect by reactive boundary conditions, due to which the model particle will be transiently immobilised by binding to the channel wall. Moreover, it may immediately rebind to the channel wall after unbinding. As we will see, this has a major effect on the particle motion.

Finally, we consider quenched, randomised channel walls with periodic boundary conditions. The walls are subject to the constraints of a fixed channel width at the boundary and a second parameter characterising the randomness of the wall configuration. We also consider that the channel in a typical system changes its shape along the trajectory of the particle beyond the length of the unit cell with its periodic boundary conditions, or that several such channels exist and can be traversed by tracers in parallel. To account for these ensemble effects we average the dynamics over an ensemble of tracer particles in an ensemble of channel geometries, all characterised by a common set of the structural parameters.

The assumption of a finite-size particle in a channel with randomised, reactive boundary conditions significantly generalises the Fick–Jacobs model. In particular, we find that the particles in their channels perform transient subdiffusion, which we analyse in terms of the anomalous diffusion exponent, the number of successful moves with respect to the number of simulations steps, and the effective long-time diffusivity, as a function of the characteristic channel geometry parameters.

The paper is structured as follows. In the subsequent section we introduce the details of the numerical approach used in this study. In Section III we discuss how the numerical results are analysed in terms of time and ensemble averaged observables. In Section IV we present the detailed results. Section V puts our findings in perspective with respect to the theory by Dagdug and Pineda, before drawing our conclusions in Section VI.

II. Simulation details

We study the diffusion through a two-dimensional channel with periodic boundary conditions in the horizontal x-direction. In the vertical y-direction the system is limited by two walls, see Fig. 1. These two walls are described by N points connected by straight lines, represented by the blue lines in Fig. 1. Thus, if one neglects the periodic boundary conditions the channel is of length N − 1. For numerical convenience the points of the channel wall reside on a lattice with unit lattice constant, which effectively determines the fundamental length scale of the system. Due to the horizontal periodic boundary conditions the two leftmost and the two rightmost wall points are identical. Their vertical distance (in the y-direction) is denoted by g, see Fig. 1. This gap opening parameter is one of the fundamental parameters of the system.

Fig. 1 Schematic of the channel with periodic boundary conditions. The blue lines depict the channel walls, while the dotted (brown) lines mark the excluded volume for a finite-size particle. Parameters: gap opening g = 6, particle size s₁ = 1, step size s₂ = 0.6, ruggedness parameter M₁ = 6, and lateral channel length 12.

We solely allow wall configurations in which the y-coordinates of nearest neighbour points within one wall differ by at most unity. This excludes the occurrence of extremely rugged walls. We consider channel walls with different contour lengths. The ruggedness parameter M₁ describes the number of displacements of size ±1 which occur in a wall. Due to the periodic boundary constraint the number of jumps directed upwards must be equal to the one for jumps downwards for both walls. Consequently, M₁ is an even number and lies in the interval [0,N − 1] (for odd N) or [0,N − 2] (for even N). We only consider configurations in which M₁ is equal for the upper and lower walls. However, this does not confine our study to symmetric channels, compare Fig. 1.

Explicitly, to generate a given channel configuration we proceed along the following steps: we start with the upper wall at coordinate (0,g) and then move to the next lattice point 1,g + j, where j ∈ {−1,0,1}, and so on. The value of the vertical offset j is determined by a random number from the unit interval. A non-zero offset is chosen when the random number is smaller or equal to the ratio of the number of remaining vertical displacements (given by M₁ minus the number of already performed nonzero vertical offsets) to the number of remaining intervals in the horizontal direction (N − 1 minus the number of already performed steps). From a second random number we decide whether a vertical offset is taken upwards or downwards. From these configurations we choose those for which the height of the final point at coordinate N is equal to g. This procedure guarantees that exactly the right number M₁ of displacements will be implemented. After this, we further randomise the configuration by multiply interchanging the indices of the individual steps. An analogous procedure is followed to construct the lower channel wall. Finally, only those configurations are accepted for which the upper and lower walls do not touch or overlap.

Thus the width profile of the system is fully described by the gap opening parameter g and the two ‘displacement vectors’ of size (N − 1) for the upper and lower walls. The ith entry of this displacement vector denotes the difference between the y-coordinate of the (i + 1)st and ith wall points for each of the two walls (i is counted from left to right). Due to the constraints mentioned above only 0, ±1 are valid entries, and the sum of all entries per wall must be zero to fulfil the periodic boundary requirements.

The position of the random walking particle is described by the position of its centre of mass, illustrated by the green cross in Fig. 1. Its motion is off-lattice. This is shown in Fig. 1, where a circle of radius s₂ (the step size of the random walk) is drawn around the particle's current position. For each step a random angle with respect to the x-axis is drawn, and the particle attempts to move its centre to the corresponding point on the dotted circle. To account for the diffusion of finite-size particles through the channel, before executing a step we check whether the distance from the current position to the wall is sufficient. To this end, the minimal distance to the wall is calculated for the trial position. Only if it is larger than the particle size s₁, the step is actually performed. This accessible space is limited by the two dotted brown lines in Fig. 1. Their vertical distance is the effective channel width of the finite-size particle, and it is the quantity to be inserted into eqn (3) and (4). If the particle aims to move at a forbidden position, the step is cancelled and the particle remains at its current position, but time is increased by one unit. This corresponds to ‘sticky’ boundary conditions, which mimic transient binding to the channel wall.

The diffusing particle is initially placed in the middle of the channel in both x and y directions. However, if such a starting position is not possible in the sense described above, the given channel configuration is dismissed and a new one chosen. T_max random walk steps are performed and the position in the x direction is traced and analysed. If not stated otherwise, for each parameter set g and M₁ the results were averaged over 25 000 configurations using the parameters N = 100, T_max = 10⁶, s₁ = 1, and s₂ = 0.6.

III. Evaluation procedure

A quantity of central interest when tracking the motion of single particles is the time-averaged mean squared displacement (TA MSD)for the ith time series x_i(t) along the horizontal direction. We use the fixed simulation time T_max, and in what follows the bar denotes a time average. The TA MSD was subsequently averaged over all configurations to obtain the ensemble and time averaged mean squared displacement (EATA MSD) Thus, the usual ensemble-averaged MSD is nothing but a special case of the ensemble- and time-averaged MSD, where the point of reference is the starting position. We also consider the ensemble averaged mean squared displacement (EA MSD) in the x-direction,where N_conf denotes the number of different configurations, in our case N_conf = 25 000. Here and in the following 〈·〉 denotes the ensemble average over channel realisations. Note that the ensemble of trajectories was obtained by generating a single random walk trajectory for each wall configuration, with the exception of the results presented in Section IV A, where the ensemble consisted of trajectories obtained with the same wall configuration.

At first we address whether the observed motion corresponds to normal Brownian diffusion characterised by the mean squared displacement

where K₁ denotes the diffusion coefficient. Deviations from this linear time dependence are termed anomalous diffusion. For a power-law form we distinguish subdiffusion (0 < α < 1) and superdiffusion (α > 1).^28,29 Subdiffusion of submicron tracers is quite widely observed in systems with a dense and highly structured environment, including the motion of tracer beads in reconstituted F-actin networks³⁰ or in wormlike micellar solutions³¹ or the passive, thermally driven motion of labelled messenger RNA in living bacteria cells,³² of the infectious pathway of adeno-associated viruses in living HeLa cells,⁵ and of endogenous granules in living fission yeast and MIN6 insulinoma cells.^33,34 Anomalous diffusion may be ergodic, that is, the time averaged MSD for sufficiently long trajectories converges to the ensemble quantity 〈x²(t)〉. However, there also exists stochastic dynamics for which so-called weak ergodicity breaking occurs, and and 〈x²(t)〉 behave differently.^35–38

As transient anomalous diffusion behaviour is not readily discernable in a conventional log–log plot of the MSD versus time, we visualise the data in terms of the MSD divided by time, as a function of the logarithm of time, see ref. 39 and 40. This method emphasises deviations from normal diffusion behaviour: curves with a negative slope represent subdiffusion. A typical plot is shown in Fig. 2. We observe a weaker form of transient subdiffusion for times up to some 10³ time steps. At longer times, there occurs a crossover to a more pronounced subdiffusive regime, which persists until approximately 10⁵ time steps. After that, terminal normal diffusion appears. Such behaviour could also be interpreted as hindered diffusion or as a process characterised by a time-dependent diffusion coefficient.^40,41 We note that both time and ensemble MSD coincide at longer times, which confirms the ergodic nature of the diffusion process. On shorter time scales the EA MSD curve lies above the EATA MSD curve due to the fact that, by construction, at the beginning of each trajectory the particle is placed in the middle of the channel. At such short times the probability that the particle sticks to the wall is greatly reduced compared to later times, and thus the EA MSD attains larger values than the EATA MSD, which averages the behaviour along the entire time series. The anomalous behaviour displayed in Fig. 2 is one of the major results of this study.

Fig. 2 Transient subdiffusion in an ensemble with gap opening parameter g = 6 and wall ruggedness parameter M₁ = 30. We plot 〈x²〉/t (blue symbols) and (red symbols) as a function of time t. Note the logarithmic abscissa. Inset: the time-local scaling exponent α(t) obtained from EATA MSD data for g = 6, M₁ = 10 (black), 30 (red) and 70 (green). The red line corresponds to the ensemble for which the curves are shown in the main plot.

In the following we study the slowing-down of the particle diffusion in terms of two quantities. First, in the normal diffusive behaviour beyond 10⁵ time steps we fit the EATA MSD data with a linear function in order to obtain the effective long time diffusion coefficient D_eff. This quantity is then compared with the theoretical value D_eff,DP given by eqn (4), since for each channel configuration we calculate D_DP(x) via Dagdug and Pineda's formula (3), where w(x) is replaced by the effective channel width. We normalise the value of the effective long time diffusivity by the corresponding value in the absence of walls, D_rel = D_eff/D₀. Second, we obtain the anomalous diffusion exponent α on time scales ranging from 10³ to 10⁵ time steps. Similar results are obtained by studying the mean maximal excursion of the particle⁴² (data not shown). This somewhat arbitrary choice of the time range in which α is fitted is motivated by the fact that most of the curves showed the more pronounced subdiffusive behaviour in this time regime. We use this averaged scaling exponent as a characteristic of our system and analyse its dependence on other parameters below. Alternatively, a time-local scaling exponent can be defined via. This quantity is plotted in the inset of Fig. 2 for the ensemble of configurations described by the main figure (red symbols) and two other ensembles with g = 6: M₁ = 10 (black symbols) and M₁ = 70 (green symbols). The inset clearly underlines our statement that normal diffusive behaviour is only observed in the short-time limit and in the long-time limit, whereas on intermediate time-scales, transient subdiffusive behaviour of variable values occurs. Thus, the fitted value of α studied in the remainder of this work can be considered as an effective value of the exponent for an intermediate time range.

IV. Results

A. Fixed channel wall configuration

Before studying the effect of different parameter sets g and M₁ to characterise the diffusion through this class of corrugated channels, we investigate in detail the features seen in Fig. 2 from simulations with fixed channel wall configurations.

We explicitly consider three exemplary configurations to illustrate the effect of the sticky boundary conditions. These configurations are characterised by the parameter pairs g = 6 and M₁ = 0, g = 4 and M₁ = 30, and g = 6 and M₁ = 50. The two configurations with non-flat walls are shown in Fig. 3. The grey curves denote the actual position of the channel walls, while the bold blue and red curves mark the region, which is inaccessible for the particle's centre.

Fig. 3 Channel wall configurations characterised by the gap opening parameter g = 6 and the wall ruggedness parameter M₁ = 50 (upper panel) and g = 4, M₁ = 30 (lower panel). The actual upper and lower walls are plotted as grey lines. The regions bounded by the red and blue curves are accessible for the particle.

While for a given wall configuration we make sure that the particle finds sufficient space in the middle of the channel where it is initially placed, it is a priori not certain that the particle can traverse the entire channel. This is the case when at some point the upper and lower walls are sufficiently close so that the finite-size particle cannot pass through such a bottleneck. Strictly speaking such a passage is impossible if there exists a horizontal interval whose width is at least of the step size s₂. Otherwise, due to the finite step size the particle can actually ‘tunnel’ through such bottlenecks. The wall configuration depicted in the upper panel of Fig. 3 does not allow such a tunnelling for the given parameters g = 6 and M₁ = 50: the red and the black curve do not overlap. The situation is different in the lower panel of Fig. 3 with g = 4 and M₁ = 30: the channel is blocked for the particle at x ≈ 12.

However, even if inaccessible regions in a given wall configuration exist but the overlap of the walls stretches over less than the distance 2 × s₂, such a narrow straight constitutes a severe entropic bottleneck for the diffusing particle: there exists an appreciable possibility that the particle repeatedly sticks to the channel walls. To quantify the influence of the sticky walls we binned the channels into 99 cells of length l and extracted from our simulations how often the particle unsuccessfully tries to move to a new position while being in the corresponding bin.

We first study the mean number n_avg of unsuccessful attempts for the three above sample configurations as a function of the particle position along the channel in Fig. 4. The green curve for the parameters g = 6 and M₁ = 0 shows that for flat walls n_avg is approximately constant. The fact that this value is finite is a consequence of the sticky boundary condition at the edges of the system: namely, move attempts which would end at a position which is forbidden due to the finite size of the particle are not executed. If reflective walls or diffusion in free space were considered, any step could be executed and then n_avg = 0. In the present case, the value of n_avg depends on the step size and the (constant) width of the channel. From Fig. 4 we deduce that n_avg ≈ 0.1056.

Fig. 4 Mean number n_avg of unsuccessful attempts to move to the next position along the channel as a function of the position x along the channel, for the three wall configurations characterised by g = 6 and M₁ = 0 (green line), g = 6 and M₁ = 50 (blue line), and g = 4 and M₁ = 30 (red line).

The blue curve in Fig. 4 for g = 6 and M₁ = 50 shows some variation as a function of x: where the channel is narrow, e.g., around x ≈ 20 in the upper panel of Fig. 3, n_avg is much higher than at locations where the channel is wider, e.g., in the middle of the channel. Comparing the minimum and maximum of n_avg along the channel, the variation of n_avg makes up a factor of approximately 7. This effect is much more pronounced in the red curve in Fig. 4 for the parameters g = 4 and M₁ = 30, corresponding to the lower panel of Fig. 3: the curve is broken as two bins of the channel are inaccessible for the particle. In the bin to the right of the channel blockage the average number of unsuccessful tries is larger than 1. In other regions of the channel n_avg attains values similar to the ones in the other two configurations. Hence, the mean number of unsuccessful motion attempts along the channel directly reflects the effective channel width and thus the local transport properties.

Additional information can be deduced from studying the probability distribution p(n_uns) of the number n_uns of unsuccessful motion attempts in a row shown in Fig. 5, where we focus on the most distinct configuration with parameters g = 4 and M₁ = 30 corresponding to the lower panel of Fig. 3. For better visibility we only consider extreme cases: namely, only the two bins with the highest and the two bins with the lowest mean number of motion attempts. In all four cases, the probability to find higher values of n_uns decreases. In regions where the channel is wide (green and blue symbols in Fig. 5) this decay is very fast, such that within our simulation time there were never more than 19 subsequent unsuccessful attempts. Otherwise, it becomes obvious that near severe bottlenecks the distribution of waiting times is much more heavy-tailed (black and red symbols in Fig. 5). Up to 100 unsuccessful attempts in a row are possible, with a probability of about 10⁻⁶.

Fig. 5 Probability distribution for the number of subsequent unsuccessful motion attempts, n_uns, obtained with the wall configuration corresponding to the lower panel of Fig. 3 with parameters g = 4 and M₁ = 30. We show the statistics for the bins centred around: x = 12.5 (black symbols), x = 13.5 (red symbols), x = 57.5 (green symbols), and x = 69.5 (blue symbols). The symbols for x = 57.5 are almost fully covered by those for x = 69.5.

With this information, let us go back to the features of Fig. 2. According to the Fick–Jacobs theory, whenever the width of the channel changes in the form of a bottleneck or a bulge, this slows down the diffusion of the particle:⁹ in the case of a bottleneck the particle may be reflected back, while in the case of a bulge the particle may execute many motion events off the minimal path along the channel. The effect of the entropic bottlenecks in our present case is amplified by the presence of the sticky boundary conditions. On all time-scales on which the particle interacts with the walls, it is slowed down in comparison to a particle moving in free space. At this point, we have to differentiate between relatively flat and rugged channel walls. In the former case it is expected that an effective, constant diffusivity characterises the motion along the entire channel, since there is no typical bulge size. In particular, no transient subdiffusion should occur. In the latter case, however, the combined effect of entropic bottlenecks and sticky boundaries induces the transient subdiffusion in the ensemble and time averaged trajectories. This can also be understood when considering two very different representatives of channels characterised by the maximal value of M₁ (see appendix).

On very long time scales, when the configurations shown in Fig. 3 are equivalent to a single step size in a coarse-grained random walk on the whole periodic structure, normal diffusive behaviour is restored, but now with a reduced diffusion coefficient. This reduced coefficient takes into account all the intermediate contacts with the channel walls. Thus, it is expected that more corrugated and/or tighter channels, which imply more frequent interaction with the walls, should show reduced values of α and a reduced effective diffusion coefficient on the ensemble level. To study these effects, in the following we systematically study the impact of the parameters g and M₁ on the transport through the channels in ensembles of size 25 000, where for each configuration we generate one random walk trajectory.

B. Simulations of channel ensembles

Two main simulation series were performed with the fixed values g = 6 and g = 4 for the gap opening parameter and 15 different values of the wall ruggedness parameter M₁, spanning the whole possible range [0,98]. The fitted values of the normalised effective diffusion coefficient in the long-time regime, D_rel, are depicted in Fig. 6 as a function of M₁. Here and in the following, solely the EATA MSD values were used, as they supply the most extensive data. The results obtained with the EA MSD data are very similar (data not shown), which is not surprising due to the ergodicity of the process at long times. In both cases, an increase of the wall ruggedness (larger M₁ values) effects slower effective diffusion. The slope of this decrease is steepest for small M₁ values and then gradually flattens off. Conversely, at fixed values of M₁ the effective diffusion is always substantially faster for g = 6 than g = 4.

Fig. 6 Normalised effective long-time diffusion coefficient D_rel from fitting of the simulations results, as a function of the wall ruggedness parameter M₁. Parameters: gap opening parameter g = 4 (black symbols) and g = 6 (blue symbols). Inset: D_rel as a function of g for M₁ = 0 (black symbols), M₁ = 10 (blue symbols), and M₁ = 50 (green symbols).

To study the impact of the gap opening parameter g in more detail, we took five different values at fixed values of the wall ruggedness parameter, namely, M₁ = 0, M₁ = 10, and M₁ = 50. Thus, we consider flat channels, slightly corrugated, and heavily corrugated channels. The corresponding results for the fitted values of the normalised effective diffusion constant D_rel are shown in the inset of Fig. 6. For fixed values of g we see once more that the diffusion is quickest when the wall is smoother, i.e., when M₁ is smaller. As was already observed in the preceding paragraph a wider gap opening at a fixed value of M₁ leads to faster diffusion. Thus, wider channels can be traversed quicker.

At first sight surprisingly, we observe that even for completely flat upper and lower channel walls (M₁ = 0, black line in the inset of Fig. 6) the diffusion in tighter channels is slowed down in comparison to the situation in free space. This is not expected in systems with reflecting boundaries, which are usually described with the Fick–Jacobs equation. However, for the finite-size particles studied here it is the result of the sticky boundary conditions at the channel walls: in a tighter channel the particle is more often close to the walls and binds transiently. Alternatively, the slow-down due to the interaction with the walls can be quantified by measuring the anomalous diffusion exponent α in the intermediate time regime, on time scales of 10³ to 10⁵ simulation steps. The fitted values for α are depicted in Fig. 7 as a function of the ruggedness M₁. The same trend as for the effective diffusion coefficient is seen for the anomalous diffusion exponent α of the transient subdiffusive regime. For increasing contour lengths of the channel wall the motion is increasingly subdiffusive. As expected, the transient subdiffusion is heavier for the tighter channel (g = 4).

Fig. 7 Anomalous diffusion exponent α as a function of wall ruggedness M₁ from the power-law fit of the EATA MSD data. Parameters: g = 4 (black) and g = 6 (blue). Inset: α as a function of g for M₁ = 10 (black), and M₁ = 50 (blue).

The dependence of α on the gap opening parameter is shown in the inset of Fig. 7. In this case, only values obtained with corrugated channels (M₁ = 10 and M₁ = 50) are shown. Note that in flat channels with M₁ = 0 no transient subdiffusion occurs, see above. Again, the curves for α are similar to the ones obtained for the effective diffusion coefficient: α is an increasing function of g and M₁ = 10 leads to less pronounced transient subdiffusion than M₁ = 50. This analogy motivates the study of the relation between α and D_rel in more detail.

In Fig. 8 we plot all α values for gap openings g = 4 and g = 6 as a function of the corresponding fitted values of D_rel. The results show that there is a strong (nonlinear) correlation between both parameters. For increasing values of D_rel the value of α also increases, with decreasing slope. The relation between both parameters is bijective, thus, both quantities are appropriate and sufficient to quantify the slow-down of the motion along the channel. For similar values of D_rel the values for α obtained with the tighter channels (g = 4) are slightly larger than those for the wider channels (g = 6). However, this fact should not be overstated: all data sets were fitted in the time interval 10³ to 10⁵, irrespective of the exact shapes of the curves. It is conceivable that a closer connection between the values for α could have been obtained by choosing the fitting time window for each curve individually. We note that this relation between the scaling exponent α and the effective diffusion coefficient was recently studied in molecular dynamics simulations dealing with proton and water mobility in (tetramethyl)urea solutions.⁴³

Fig. 8 Anomalous diffusion exponent α of the transiently subdiffusive regime as a function of the long-time diffusion exponent D_rel with gap opening parameters g = 4 (black) and g = 6 (blue).

As mentioned above, in an ensemble of systems with corrugated boundaries not all channels can be traversed completely. If a channel is blocked somewhere, the corresponding squared displacement of the particle position has an upper limit. On an ensemble level these trajectories will reduce the average values of α and D_rel, since this is one of the main factors contributing to the anomalous features of Fig. 2. Thus, it is important to extract from our simulations solely the unblocked configurations. The corresponding parameter f_open is plotted as a function of the ruggedness M₁ for fixed gap openings g = 4 (black symbols) and g = 6 (blue symbols) in Fig. 9, and as a function of the gap opening g (for M₁ = 10 and M₁ = 50) in the inset of Fig. 9. An inspection of Fig. 9 shows that this fraction is a decreasing function of M₁ and an increasing function of g. Overall, the curves look similar to the ones of the normalised effective diffusion coefficient D_rel (compare Fig. 6).

Fig. 9 Fraction f_open of unblocked configurations as a function of the channel ruggedness M₁ for channel opening g = 4 (black) and g = 6 (blue). Inset: f_open as a function of g for M₁ = 10 (black) and M₁ = 50 (blue).

To better understand why above a certain threshold of the boundary ruggedness parameter M₁ the effective diffusion coefficient only decreases slightly (see Fig. 6), it is instructive to study the average weighted effective width w_wgt of the channels. Here, weighted means that for any blocked channel the width is set to zero. The weighted width is plotted as a function of M₁ in Fig. 10. For increasing, yet small values of M₁ the effective weighted channel width decreases, until it reaches a shallow minimum. For larger values of M₁ it rises to higher values than at the origin, slowly approaching a plateau.

Fig. 10 Weighted effective channel width w_wgt as a function of channel ruggedness M₁ for gap opening g = 6 (blue line and symbols) and g = 4 (black line and symbols).

Due to the procedure to generate the wall configurations, for those geometries with a low value of both g and M₁, there is a substantial probability that the passage is blocked: after a displacement in one wall which decreases the effective channel width to a critical value, there are not many possibilities which could increase it again. At the same time the maximum width is rather restricted due to the low total number of displacements. The situation becomes different above a certain threshold: heavily rugged walls which reach a critical value have many more opportunities to become wider again and they can reach far bigger widths. Thus, the average width of those traversable channels is larger for more rugged wall configurations. This facilitates the transport through these configurations, as the sticky boundaries are further away. However, as remarked earlier, more rugged walls slow down the diffusion due to the occurrence of bulges and constrictions,⁹ such that we have two opposing effects, which mostly (almost) cancel each other. Consequently, D_rel remains nearly constant above a threshold value of M₁ (compare Fig. 6).

Fig. 11 shows the normalised effective diffusion coefficient D_rel as a function of the weighted channel width w_wgt. For better visibility data points with identical M₁ values are connected by lines (M₁ = 0: black line, M₁ = 10: blue line, and M₁ = 50: green line). While for a fixed value of M₁ more spacious channels allow faster diffusion, the heavy scatter between values of D_rel obtained with similar values of w_wgt (grey symbols) shows that the knowledge of the mean channel width of an ensemble of channels alone is insufficient to predict the transport properties. At the same time an inspection of Fig. 6 shows that a single of the parameters g and M₁ is not sufficient to make such a prediction. At least two such parameters are needed.

Fig. 11 Normalised effective diffusion coefficient D_rel as a function of the weighted effective channel width w_wgt for M₁ = 0 (black), M₁ = 10 (blue line), M₁ = 50 (green line) and other values of M₁ (grey symbols).

C. Displacement autocorrelation function

Another possibility to assess the transport properties of our channel systems is to calculate the velocity autocorrelation function. This quantity is particularly interesting since it is long known that in the case of two-dimensional (hard-disk) fluids the diffusion coefficient calculated via the Green–Kubo formula diverges due to the power-law scaling of the autocorrelation function.^44–47

Since, by definition, there is no velocity in our random walk simulation we calculate the displacement autocorrelation function along the channel C_xx,i(t) for the ith time series, which is defined as:

This quantity can then be averaged over different realisations to obtain the ensemble-averaged displacement correlation function, . The data show that this function stays negative for small and intermediate time regimes and slowly levels off to zero. On longer timescales the function fluctuates around zero. Thus, in order to conveniently create a log–log plot in Fig. 12 we show the absolute value of this function for two ensembles characterised by the parameters g = 4, M₁ = 30 (upper panel) and g = 6, M₁ = 50 (lower panel).

Fig. 12 Absolute value of the ensemble-averaged displacement correlation function for two ensembles of size N_conf = 10 000 characterised by g = 4, M₁ = 30 (upper panel) and g = 6, M₁ = 50 (lower panel). The red symbols correspond to the values obtained from our simulations, the dashed lines to a fit consisting of the sum of two stretched exponentials: C_fit(t) = c₁ exp(at^b) + c₂ exp(ct^d). Fit parameters: a = −3.656, b = 0.1781, c = −1.621, d = 0.1479, c₁ = 0.2387, c₂ = 0.005139 (upper panel) and a = −5.033, b = 0.1329, c = −5.601, d = 0.07081, c₁ = 0.9091, c₂ = 0.3142 (lower panel).

For t ≲ 10⁴ time steps both data sets are almost of power-law form, however, on a closer look two shoulders may be discerned around such a power-law trend. The curves can be approximated very well by a fit function, which is the sum of two stretched exponential functions: C_fit(t) = c₁ exp(at^b) + c₂ exp(ct^d). The data are thus consistent with a fast, exponential decay. Note that at times shorter than 10⁵ time steps the particle still has not explored the entire channel. At longer times the correlations will clearly be negligible so that the Green–Kubo integral converges.

V. Comparison with Dagdug and Pineda's formula

In order to compare our results obtained from ensembles of channel configurations with the results of Dagdug and Pineda, we make a simplifying assumption: for all unblocked configurations, we determine D_DP(x) from eqn (3) and subsequently D_eff,DP from eqn (4). For all blocked configurations the effective diffusivity D_eff,DP = 0 accounts for the fact that on long time-scales particles in these configurations do not contribute significantly to the MSD. Finally, we average over all configurations in the ensemble and normalise through division by the diffusion coefficient in free space, D_rel,DP = 〈D_eff,DP〉/D₀. We compare these values with our fitted values of the normalised diffusion coefficient in the upper panel of Fig. 13, where data points with the same value of the gap opening g are represented in the same colour.

Fig. 13 Upper panel: normalised effective diffusion coefficient D_rel from our simulations in the sticky, corrugated channel as a function of the value D_rel,DP obtained from the results (3) of Dagdug and Pineda. Black: gap opening parameter g = 4, red: g = 6, green: g = 8, blue: g = 10, and cyan: g = 12. The two grey lines mark the range of ±10% around the theoretical value. Lower panel: distributions of effective diffusion coefficients for individual unblocked channel configurations obtained with Dagdug and Pineda's formula. Parameters: g = 6 in all four cases. M₁ = 6: black, M₁ = 20: red, M₁ = 50: green, M₁ = 98: blue.

This is another major result of our study: as demonstrated in the upper panel of Fig. 13, for fixed values of the gap opening g there is a linear relation between D_rel and D_rel,DP. Most of the data points are located within a 10% confidence interval around Dagdug and Pineda's values. More explicitly, a linear fit of the data points obtained in a simulation series yields the following results. For g = 4, the fitted relation between the two is D_rel = (0.018 ± 0.004) + (0.853 ± 0.009)D_rel,DP, while for g = 6 we find D_rel = (0.032 ± 0.095) + (0.903 ± 0.007)D_rel,DP. For g = 8, 10, and 12 we did not fit the data as there were only three values available.

The fact that the slope of the fits is somewhat below 1 shows that Dagdug and Pineda's formula, which only applies to systems with perfectly reflecting boundaries, slightly overestimates the diffusion coefficient compared to our system with sticky boundary conditions. Thus, as expected the additional interaction with the boundaries further slows down the diffusion, which is already reduced by the occurrence of entropic bottlenecks. This reasoning is further substantiated by the observation that in tighter channels (with g = 4), where these surface effects play a more important role, the slope of the conversion formula is smaller, and the deviation from Dagdug and Pineda's formula is more pronounced. Given the quite intricate form of the effective channel width (see Fig. 1 and 3) it is not feasible to quantify this effect analytically. However, in the future other values of s₁ and s₂ could be considered to study the magnitude of the correction terms numerically.

Deviations of our results from Dagdug and Pineda's formula are also based on the fact that their analysis applies to one given channel configuration. Driven by the physical application we here consider an ensemble of different channel wall configurations, solely defined by the fixed macroscopic parameters g and M₁. Individual configurations may therefore differ considerably, even if they share the values of g and M₁. This is illustrated in the lower panel of Fig. 13, where for all unblocked configurations the expected effective diffusion coefficient was calculated with Dagdug and Pineda's formula. For g = 6 and four different values of M₁ we see that the distribution of values of D_rel,DP is far from narrow. The values increasingly scatter for more rugged conformations (higher values of M₁). Thus, the semi-quantitative agreement of our data with the theoretical prediction on an ensemble level is indeed remarkable.

VI. Conclusion and outlook

We studied the motion of finite-size particles through randomly corrugated channels with sticky walls, observing transient anomalous diffusion of the particles in their passage of the channel. We also obtained the long-time diffusion coefficient for this motion on time scales over which normal Brownian diffusion is restored. The control parameters in our study were the gap opening parameter fixing the distance between the channel walls at the entrance and the exit of the channel, as well as the wall ruggedness parameter setting the maximal variation of the channel wall configuration. We quantified the dependence of the anomalous diffusion exponent and the long-time diffusion coefficient as a function of the ruggedness and gap opening, and showed that both quantities are in fact correlated. We especially analysed the blocked channels, which the particle cannot fully traverse. The long-time effective diffusion coefficient was shown to agree well with the prediction for point-like particles in channels with reflecting boundary conditions by Dagdug and Pineda.

It will be interesting to study different models for the construction of the random channel walls, for instance, by adding correlations in the growth algorithm. Moreover it will be of relevance to replace the single particle scenario and consider an ensemble of diffusing particles with excluded volume.

Appendix: V-shaped and saw-tooth channels

In order to illustrate that two channels with the same values of g and M₁ can have very different transport properties, we consider here two configurations with the maximal value of M₁ = 98. One is a V-shaped channel for which both upper and lower boundaries move upwards in the left half of the channel and back downwards in the right half of the channel. Technically this means that for both walls every vertical offset is plus one in the left half and minus one in the right half. The other channel has walls with a saw-tooth shape. Again the values of the offsets are equal for the upper and lower wall, however, this time the offsets alternate between +1 and −1. In Fig. 14, which is analogous to Fig. 2, we plot the EATA MSD divided by time of 10 000 trajectories in these two channels.

Fig. 14 Transport in two sample configurations with M₁ = 98 and g = 4. We plot 〈δ²〉/t as a function of time t for a V-shaped channel (red curve) and a saw-tooth channel (black curve). The results are averaged over N_conf = 10 000 configurations.

For the V-shaped channel (red curve), already on short time scales the particle often encounters the channel walls such that the time-local diffusion coefficient plotted here attains values considerably smaller than that in free space (0.18). The corresponding values decrease further until around 50 time steps after which the value of the diffusion coefficient remains nearly constant.

Even though the words saw-tooth seem to imply a rather rugged wall shape, the bulges in the effective width are rather small, such that the particle experiences a nearly flat channel, beyond the time scale for crossing a single tooth. Accordingly, the diffusion coefficient at the shortest timescales is very close to the value in free space and after a small decrease it is nearly constant on all timescales. As for really flat channels, only the effective long-time diffusion coefficient is reduced, but no transient subdiffusion is observed.

Thus, these two channels constitute extreme cases of the distribution shown in the lower panel of Fig. 13, where it can be seen that for M₁ = 98 the distribution of effective diffusivities is very broad.

Acknowledgements

MB and AG acknowledge funding from the German Federal Ministry for Education and Research, and RM from the Academy of Finland within the FiDiPro scheme.

References

1. L. W. Gelhar, C. Welty and K. R. Rehfeldt, Water Resour. Res., 1992, 28, 1955; B. Berkowitz, A. Cortis, M. Dentz and H. Scher, Rev. Geophys., 2006, 44, RG2003.
2. R. Kimmich, NMR: Tomography, Diffusometry, Relaxometry, Springer, Berlin, 2011.
3. J. Kärger and D. M. Ruthven, Diffusion in Zeolites and Other Microporous Solids, Wiley, New York, 1992.
4. T. Kühn, T. O. Ihalainen, J. Hyväluoma, N. Dross, S. F. Willman, J. Langowski, M. Vihinen-Ranta and J. Timonen, PLoS One, 2011, 6, e22962.
5. G. Seisenberger, M. U. Ried, T. Endreß, H. Büning, M. Hallek and C. Bräuchle, Science, 2001, 294, 1929.
6. A. G. Cherstvy, A. V. Chechkin and R. Metzler, Soft Matter, 2014, 10, 1591.
7. L. Liu, P. Li and S. A. Asher, Nature, 1999, 397, 141.
8. P. S. Burada, P. Hänggi, F. Marchesoni, G. Schmid and P. Talkner, ChemPhysChem, 2009, 10, 45.
9. R. Zwanzig, J. Phys. Chem., 1992, 96, 3926.
10. D. Reguera and J. M. Rub, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2001, 64, 061106.
11. P. Kalinay and J. K. Percus, J. Chem. Phys., 2005, 122, 204701.
12. P. Kalinay and J. K. Percus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2005, 72, 061203.
13. P. Kalinay and J. K. Percus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2006, 74, 041203.
14. A. M. Berezhkovskii, M. A. Pustovoit and S. M. Bezrukov, J. Chem. Phys., 2007, 126, 134706.
15. P. Kalinay and J. K. Percus, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 021103.
16. S. Martens, G. Schmid, L. Schimansky-Geier and P. Hänggi, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 051135.
17. A. Berezhkovskii and A. Szabo, J. Chem. Phys., 2011, 135, 074108.
18. P. Kalinay, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2013, 87, 032143.
19. S. Martens, A. V. Straube, G. Schmid, L. Schimansky-Geier and P. Hänggi, Phys. Rev. Lett., 2013, 110, 010601.
20. R. M. Bradley, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 80, 061142.
21. L. Dagdug and I. Pineda, J. Chem. Phys., 2012, 137, 024107.
22. S. Lifson and J. L. Jackson, J. Chem. Phys., 1962, 36, 2410.
23. R. Festa and E. D'Agliano, Physica A, 1978, 90, 229.
24. J. J. Kasianowicz, E. Brandin, D. Branton and D. W. Deamer, Proc. Natl. Acad. Sci. U. S. A., 1996, 93, 13770; L. Brun, M. Pastorizo-Gallego, G. Oukhaled, J. Math, L. Bacri, L. Auvray and J. Pelta, Phys. Rev. Lett., 2008, 100, 158302.
25. L. Dagdug, A. M. Berezhkovskii, Y. A. Makhnovskii and V. Y. Zitserman, J. Chem. Phys., 2008, 129, 184706.
26. W. Riefler, G. Schmid, P. S. Burada and P. Hänggi, J. Phys.: Condens. Matter, 2010, 22, 454109.
27. M. L. Henle, B. DiDonna, C. D. Santangelo and A. Gopinathan, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2008, 78, 031118.
28. R. Metzler and J. Klafter, Phys. Rep., 2000, 339, 1.
29. R. Metzler and J. Klafter, J. Phys. A: Math. Gen., 2004, 37, R161.
30. I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch and D. A. Weitz, Phys. Rev. Lett., 2004, 92, 178101.
31. J.-H. Jeon, N. Leijnse, L. Oddershede and R. Metzler, New J. Phys., 2013, 15, 045011.
32. I. Golding and E. C. Cox, Phys. Rev. Lett., 2006, 96, 098102.
33. J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede and R. Metzler, Phys. Rev. Lett., 2011, 106, 048103.
34. S. M. A. Tabei, S. Burov, H. Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L. H. Philipson, A. R. Dinner and N. F. Scherer, Proc. Natl. Acad. Sci. U. S. A., 2013, 110, 4911.
35. E. Barkai, Y. Garini and R. Metzler, Phys. Today, 2012, 65(8), 29.
36. S. Burov, J.-H. Jeon, R. Metzler and E. Barkai, Phys. Chem. Chem. Phys., 2011, 13, 1800.
37. A. G. Cherstvy, A. V. Chechkin and R. Metzler, New J. Phys., 2013, 15, 083039.
38. A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys., 2013, 15, 20220.
39. P. A. Netz and T. Dorfmüller, J. Chem. Phys., 1995, 103, 9074.
40. M. Saxton, Biophys. J., 1996, 70, 1250.
41. A. M. Berezhkovskii, L. Dagdug and S. M. Bezrukov, Biophys. J., 2014, L09.
42. V. Tejedor, O. Bénichou, R. Voituriez, R. Jungmann, F. Simmel, C. Selhuber-Unkel, L. B. Oddershede and R. Metzler, Biophys. J., 2010, 98, 1364.
43. J. Xu, T. Yamashita, N. Agmon and G. A. Voth, J. Phys. Chem. B, 2013, 117, 15426.
44. B. J. Alder and T. E. Wainwright, Phys. Rev. A: At., Mol., Opt. Phys., 1970, 1, 18.
45. M. H. Ernst, E. H. Hauge and J. M. J. van Leeuwen, Phys. Rev. Lett., 1970, 25, 1254.
46. J. R. Dorfman and E. G. D. Cohen, Phys. Rev. Lett., 1970, 25, 1257.
47. B. Lin, S. A. Rice and D. A. Weitz, Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top., 1995, 51, 423.

---