Diffusion of finite-size particles in channels with random walls
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] D ec Diffusion of finite-size particles in channels with random walls
Maximilian Bauer,
1, 2
Aljaˇz Godec,
1, 3 and Ralf Metzler
1, 4, ∗ Institute of Physics and Astronomy, University of Potsdam, D-14476 Potsdam-Golm, Germany Physics Department, Technical University of Munich, Garching, Germany National Institute of Chemistry, Ljubljana, Slovenia Physics Department, Tampere University of Technology, FI-33101 Tampere, Finland
Diffusion of chemicals or tracer molecules through complex systems containing irregularly shapedchannels 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. Inthis study we use numerical simulations to consider the transport of finite-sized particles throughasymmetrical two-dimensional channels. Additionally, we examine transient binding of the moleculesto the channel walls by applying sticky boundary conditions. With the application of diffusingpathogens in hydrogels in mind, we consider an ensemble of particles diffusing in independentchannels, which are characterised by common structural parameters. We compare our results forthe long-time effective diffusion coefficient with a recent theoretical formula obtained by Dagdugand Pineda [
J. Chem. Phys. , 2012, , 024107].
I. INTRODUCTION
Civil aviation traffic has been broadly democratisedwithin the last few decades, leading to ever-increasingnumbers of both business and pleasure passengers. Thedownside of this increased worldwide connectivity is thevery rapid global spreading of diseases [1]. Simultane-ously, new infectious diseases are emerging, driven byhuman or ecologic causes [2], and microorganisms are de-veloping various forms of multiple drug resistance [3, 4].It is therefore imperative to develop more rapid, mobile,and reliable methods for pathogen detection. One of thepaths being followed is the development of smart hydro-gels, which respond to various stimuli [5–8].The scenario we have in mind is that viral pathogensdiffuse into the hydrogel, where they bind to specificsites and effect a local mechanical deformation. To causemacroscopic effects, the pathogens’ size must be compa-rable to the typical mesh size of the hydrogel. A similarsituation was investigated in a single particle trackingstudy of submicron tracer beads in reconstituted F-actinnetworks, observing anomalous diffusion [9] h r ( t ) i ≃ t α , (1)where the anomalous diffusion exponent α depended onthe ratio of the tracer bead size versus the typical meshsize of the actin network: for relatively large beads α was shown to decrease to zero, while for small beads α converged to one, the case of normal diffusion [9]. Sub-diffusion of the type (1) was observed for similarly-sizedtracer beads in wormlike micellar solutions [10], as well asfor the motion of various objects in the macromolecularlycrowded cytoplasm of living cells such as the infectiouspathway of adeno-associated viruses in living HeLa cells[11], as well as for submicron lipid and insulin granules ∗ Electronic address: [email protected] in living fission yeast and MIN6 insulinoma cells [12, 13].For reviews on anomalous diffusion, see Refs. [14, 15], andfor subdiffusion in crowded systems, see Refs. [16–18].Given these experimental findings it is pertinent to askwhether the motion of pathogens in a hydrogel is equallysubdiffusive. Here we address this question by consider-ing the path of the pathogen through the hydrogel as themotion of a particle of finite size through a tortuous, cor-rugated channel with varying width. The particle getsrepeatedly held up by bottlenecks in the channel, andmay transiently bind to the channel walls, representingthe interaction with specific binding site for the pathogenin the hydrogel. Using an ensemble of channel geometriesin our numerical analysis, we account for the differentpaths the pathogen can take through the hydrogel. Wefind that indeed the particle in their channel performstransient subdiffusion, that we analyse in terms of theanomalous diffusion exponent, the number of successfulmoves with respect to the number of simulations steps,and the effective long-time diffusivity, as function of thecharacteristic channel geometry parameters.The theoretical description of the motion of particles inconfined geometries like channels (2D) or pores (3D) withvarying width has a long-standing history (see [19] andreferences therein). In a seminal work Zwanzig deriveda modified Fick-Jacobs equation which is at the basis ofmost subsequent quasi-one-dimensional descriptions [20], ∂G ( x, t ) ∂t = ∂∂x D ( x ) w ( x ) ∂∂x G ( x, t ) w ( x ) (2)Here G ( x, t ) describes the local concentration of particlesat position x and time t , w ( x ) the width of the channelat position x and most importantly D ( x ) is an effectiveposition-dependent diffusion coefficient. Subsequently,several different forms for D ( x ) were derived and appliedto various systems [20–30]. Taking along only first orderderivatives of the width profile w ( x ), Kalinay and Per-cus [24] (see also Martens et al. [27]) obtained the closedform result D KP ( x ) = D arctan( w ′ ( x ) / w ′ ( x ) / , (3)where D denotes the position-independent diffusion co-efficient in the absence of confinement. We note thatthe presence of an x -dependent diffusivity in free spaceis sufficient to effect various forms of anomalous diffusion[31–33].However, the above forms of D ( x ) are restricted tosymmetric channels, i.e., channels with a straight centre-line. This constraint was removed in an approach byBradley [35], which was subsequently generalised by Dag-dug and Pineda to [36] D DP ( x ) = D arctan (cid:16) y ′ ( x ) + w ′ ( x )2 (cid:17) w ′ ( x ) − arctan (cid:16) y ′ ( x ) − w ′ ( x )2 (cid:17) w ′ ( x ) , (4)where y ( x ) denotes the vertical position of the centre-line at horizontal position x . Note that Eq. (4) gener-alises all the previous results for D ( x ), for instance, oneobtains the result (3) for symmetrical channels by setting y ′ ( x ) = 0 [36].As detailed by Zwanzig [20], in a system with periodicboundary conditions the effective diffusion coefficient inthe long-time regime, D eff , is obtained by using the fol-lowing formula introduced by Lifson and Jackson [37] andgeneralised by Festa and Galleani d’Agliano [38],1 D eff = (cid:28) D ( x ) w ( x ) (cid:29) h w ( x ) i , (5)where h·i denotes the average over one period. Thus,for any channel with width profile w ( x ) and centre-line y ( x ), Eq. (4) can be used to calculate D DP ( x ), which inturn is used in Eq. (5) to calculate the effective diffusioncoefficient D eff , DP .The above theories based on the Fick-Jacobs formal-ism apply to point-like particles, that is, the particle canmove along the channel as long as the width is not equalto zero. For our scenario of pathogens moving in a chan-nel, we argued that the pathogen size is comparable tothe channel width. Thus, the pathogen can only fullycross the channel when the width profile w ( x ) at anyposition is larger than the particle size. Systems contain-ing finite-size particles were indeed studied in literature[39, 40]. Essentially, in all formulae the channel width w ( x ) has to be replaced by an effective channel width.Another modification with respect to the Fick-Jacobsapproach that we consider here concerns the boundaryconditions. Usually, reflecting boundary conditions areused, i.e., when the particle collides with the channelwalls, its perpendicular motion is simply reversed. In the pathogen-hydrogel scenario, the pathogen will (tran-siently) bind to specific binding sites incorporated intothe hydrogel. Here we explicitly include such effects byreactive boundary conditions, due to which the modelparticle will transiently be immobilised by binding to thechannel wall. Moreover, it may immediately rebind tothe channel wall after unbinding. As we will see, this hasa major effect on the particle motion.Finally, the conventional Fick-Jacobs approach de-scribes the motion of particles in a single channel withquenched geometry. However, for the application to thepathogen motion in the hydrogel, we mimic the possi-bility for particles to move on different paths across thehydrogel by sampling an ensemble of tracer particles inan ensemble of channel geometries. This ensemble ofchannels is characterised by a common set of structuralparameters.The present study therefore represents an applicationof the Fick-Jacobs approach to the biophysical problemof pathogen motion in a complex, confining environment.At the same time, however, it significantly generalises theFick-Jacobs model. The paper is structured as follows.In the subsequent section we introduce the details of thenumerical approach used in this study. In section 3 wediscuss how the numerical results are analysed in termsof time and ensemble averaged observables. In section4 we present the detailed results. Section 5 puts ourfindings in perspective with respect to theory by Dagdugand Pineda, before drawing our conclusions in section 6. II. SIMULATION DETAILS
We study the diffusion through a two-dimensionalchannel with periodic boundary conditions in horizontal x -direction. In vertical y -direction the system is limitedby two walls, see Fig. 1. These two walls are describedby N points connected by straight lines, represented bythe blue lines in Fig. 1. For numerical convenience thepoints of the channel wall reside on a lattice with unitlattice constant, which effectively determines the funda-mental length scale of the system. Due to the horizontalperiodic boundary conditions the two leftmost and thetwo rightmost wall points are identical. Their verticaldistance (in y -direction) is denoted by g , see Fig. 1. Thisis one of the fundamental parameters of the system andis called gap opening parameter in the following.We only consider channel configurations withouttouching or overlapping wall. Moreover, we solely allowwall configurations in which the y -coordinates of nearestneighbour points within one wall differ by at most unity.This excludes the occurrence of extremely rugged walls.Thus the size of the system is fully described by the gapopening parameter g and the two ‘displacement vectors’of size ( N − × i th entry of this displacement vector denotes the differencebetween the y -coordinate of the ( i + 1) st and i th wallpoints for each of the two walls ( i is counted from left (cid:0) ✁✂✄ ✂☎ FIG. 1: Schematic of the channel with periodic boundaryconditions. The blue lines depict the channel walls, while thedotted (brown) lines mark the excluded volume for a finite-sized particle. Parameters: gap opening g = 6, particle size s = 1, step size s = 0 .
6, ruggedness parameter M = 6, andlateral channel length N = 12. to right). Due to the constraints mentioned above only0 , ± M describes the number of dis-placements of size ± M is an even number and liesin the interval [0 , N −
1] (for odd N ) or [0 , N −
2] (foreven N ). We only consider configurations in which M isequal for the upper and lower walls. However, this doesnot confine our study to symmetric channels, compareFig. 1.The position of the random walking particle is de-scribed by the position of its centre of mass, illustratedby the green cross in Fig. 1. Its motion is off-lattice. Thisis shown in Fig. 1, where a circle of radius s (the stepsize of the random walk) is drawn around the particle’scurrent position. For each step a random angle with re-spect to the x -axis is drawn, and the particle attempts tomove its centre to the corresponding point on the dottedcircle. To account for the diffusion of finite-sized par-ticles through the channel, before executing a step wecheck whether the distance from the current position tothe wall is sufficient. To this end, the minimal distanceto the wall is calculated for the trial position. Only ifit is larger than the particle size s , the step is actuallyperformed. This accessible space is limited by the twodotted brown lines in Fig. 1. Their vertical distance isthe effective channel width for the finite-sized particle,and it is the quantity to be inserted into Eqs. (4) to (5).If the particle aims to move at a forbidden position, thestep is cancelled and the particle remains at its current position, but time is increased by one unit. This cor-responds to ‘sticky’ boundary conditions, which mimictransient binding to the channel wall.The diffusing particle is initially placed in the middle ofthe channel in both x and y directions. However, if sucha starting position is not possible in the sense describedabove, the given channel configuration is dismissed and anew one chosen. T max random walk steps are performedand 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 configurationsusing the parameters N = 100, T max = 10 , s = 1, and s = 0 . III. EVALUATION PROCEDURE
A quantity of central interest when tracking the mo-tion of single particles is the time-averaged mean squareddisplacement (TA MSD) δ i ( t ) = 1 T max − t T max − t Z dt ′ h x i ( t ′ + t ) − x i ( t ′ ) i (6)for the i th time series x i ( t ) along the horizontal direction.We use the fixed simulation time T max , and in what fol-lows the bar denotes a time average. The TA MSD wassubsequently averaged over all configurations to obtainthe ensemble and time averaged mean squared displace-ment (EATA MSD) D δ ( t ) E = 1 N conf N conf X i =1 δ i ( t ) . (7)Thus, the usual ensemble-averaged MSD is nothing buta special case of the ensemble- and time-averaged MSD,where the point of reference is the starting position. Wealso consider the ensemble averaged mean squared dis-placement (EA MSD) in x -direction, h x ( t ) i = 1 N conf N conf X i =1 h x i ( t ) − x i (0) i , (8)where N conf denotes the number of different configura-tions, in our case N conf = 25 , h·i denotes the ensemble average over channel re-alisations.As transient anomalous diffusion behaviour is not read-ily discernable in a conventional log-log plot of the MSDversus time, we visualise the data in terms of the MSDdivided by time, as function of the logarithm of time,see Refs. [41, 42]. This method emphasises deviationsfrom normal diffusion behaviour: curves with a negativeslope represent subdiffusion. A typical plot is shown inFig. 2. We observe a weaker form of subdiffusion fortimes from roughly 10 to 10 time steps. Around 10 < x ( t ) > /t, < δ ( t ) > /t t
0 20 40 60 80 100 y x FIG. 3: Channel wall configurations characterised by the gapopening 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 region bounded by the red and blue curves are accessiblefor the particle. IV. RESULTSA. Fixed channel wall configuration
Before studying the effect of different parameter sets g and M to characterise the diffusion through this class ofcorrugated channels, we investigate in detail the featuresseen in Fig. 2 from simulations with fixed channel wallconfigurations.We explicitly consider three exemplary configurationsto illustrate the effect of the sticky boundary conditions.These configurations are characterised by the parameterpairs g = 6 and M = 0, g = 4 and M = 30, and g = 6and M = 50. The two configurations with non-flat wallsare shown in Fig. 3. The grey curves denote the actualposition of the channel walls, while the bold blue andred curves mark the region, which is inaccessible for theparticle’s centre.While for a given wall configuration by the choice of n a v g x g=6, M = 0g=6, M =50g=4, M =30 FIG. 4: Mean number n avg of unsuccessful attempts to moveto the next position along the channel as a function of theposition x along the channel, for the three wall configurationscharacterised by g = 6 and M = 0 (green line), g = 6 and M = 50 (blue line), and g = 4 and M = 30 (red line). the gap opening parameter g we make sure that the par-ticle finds sufficient space in the middle of the channelwhere it is initially placed, it is a priori not certain thatthe particle can traverse the entire channel. This is thecase when at some point the upper and lower walls over-lap. Strictly speaking, however, a passage is impossibleonly if there exists an overlap area whose width is atlast of the step size s . Otherwise, due to the finite stepsize the particle can actually ‘tunnel’ through such bot-tlenecks. The wall configuration depicted in the upperpanel of Fig. 3 does not allow such a tunnelling for thegiven parameters g = 6 and M = 50: the red and theblack curve do not overlap. The situation is different inthe lower panel of Fig. 3 with g = 4 and M = 30: thechannel is blocked for the particle at x ≈ × s , such a narrow straightconstitutes a severe entropic bottleneck for the diffus-ing particle: there exists an appreciable possibility thatthe particle repeatedly sticks to the channel walls. Toquantify the influence of the sticky walls we binned thechannels into 99 cells of length 1 and extracted from oursimulations how often the particle unsuccessfully tries tomove to a new position while being in the correspondingbin.We first study the mean number n avg of unsuccessfulattempts for the three above sample configurations asa function of the particle position along the channel inFig. 4. The green curve for the parameters g = 6 and M = 0 shows that for flat walls n avg is approximatelyconstant. The fact that this value is finite is a conse-quence of the sticky boundary condition at the edges ofthe system: namely, move attempts which would end ata position which is forbidden due to the finite size of the particle are not executed. If reflective walls were consid-ered, any step could be executed and then n avg = 0. Inthe present case, the value of n avg depends on the stepsize and the (constant) width of the channel. From Fig. 4we deduce that n avg ≈ . g = 6 and M = 50 showssome variation as function of x : where the channel isnarrow, e.g., around x ≈
20 in the upper panel of Fig. 3, n avg is much higher than at locations where the channelis wider, e.g., in the middle of the channel. Comparingthe minimum and maximum of n avg along the channel,the variation of n avg makes up a factor of approximately7. This effect is much more pronounced in the red curvein Fig. 4 for the parameters g = 4 and M = 30, corre-sponding to the lower panel of Fig. 3: the curve is brokenas two bins of the channel are inaccessible for the parti-cle. In the bin to the right of the channel blockage theaverage number of unsuccessful tries is larger than 1. Inother regions of the channel n avg attains values similarto the ones in the other two configurations. Hence, themean number of unsuccessful motion attempts along thechannel directly reflects the effective channel width andthus the local transport properties.Additional information can be deduced from studyingthe probability distribution p ( n uns ) of the number n uns of unsuccessful attempts in a row shown in Fig. 5, wherewe focus on the most distinct configuration with param-eters g = 4 and M = 50 corresponding to the lowerpanel of Fig. 3. For better visibility we only considerextreme cases: namely, only the two bins with the high-est and the two bins with the lowest mean number ofmotion attempts. In all four cases, the probability tofind higher values of n uns decreases. In regions where thechannel is wide (green and blue symbols in Fig. 5) thisdecay is very fast, such that within our simulation timethere were never more than 19 subsequent unsuccessfulattempts. Otherwise, it becomes obvious that near se-vere bottlenecks the distribution of waiting times is muchmore heavy-tailed (black and red symbols in Fig. 5). Upto 100 unsuccessful attempts in a row are possible, witha probability of about 10 − .With this information, let us go back to the featuresof Fig. 2. According to the Fick-Jacobs theory, wheneverthe width of the channel changes in the form of a bot-tleneck or a bulge, this slows down the diffusion of theparticle [20]: in the case of a bottleneck the particle maybe reflected back, while in the case of a bulge the parti-cle may execute many motion events off the minimal pathalong the channel. The effect of the entropic bottlenecksin our present case is amplified by the presence of thesticky boundary conditions. On all time-scales on whichthe particle interacts with the walls, it is slowed down incomparison to a particle moving in free space. This in-duces the transient subdiffusion in the ensemble and timeaveraged trajectories. On very long time scales, when theconfigurations shown in Fig. 3 are equivalent to a singlestep size in a coarse-grained random walk on the wholeperiodic structure, normal diffusive behaviour is restored, -8 -6 -4 -2
0 20 40 60 80 100 p n uns x=12.5x=13.5x=57.5x=69.5 FIG. 5: Probability distribution for the number of subsequentunsuccessful motion attempts, n uns , obtained with the wallconfiguration corresponding to the lower panel of Fig. 3 withparameters g = 4 and M = 30. We show the statistics forthe bins centred around: x = 12 . x = 13 . x = 57 . x = 69 . x = 57 . x = 69 . but now with a reduced diffusion coefficient. This re-duced coefficient takes into account all the intermediatecontacts with the channel walls. Thus, it is expectedthat more corrugated and/or tighter channels, which im-ply more frequent interaction with the walls should showreduced values of α and a reduced effective diffusion co-efficient on the ensemble level. To study these effects,in the following we systematically study the impact ofthe parameters g and M on the transport through thechannels in ensembles of size 25,000. B. Simulations of channel ensembles
Two main simulations series were performed with thefixed values g = 6 and g = 4 for the gap opening parame-ter and 15 different values of the wall ruggedness param-eter M , spanning the whole possible range [0 , D rel , are depicted in Fig. 6as function of M . Here and in the following, solely theEATA MSD values were used, as they supply the mostextensive data. The results obtained with the EA MSDdata are very similar (data not shown), which is not sur-prising 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 slopeof this decrease is steepest for small M values and thengradually flattens off. Conversely, at fixed values of M the effective diffusion is always substantially faster for g = 6 than g = 4.To study the impact of the gap opening parameter g in more detail, we took five different values at fixed val- D r e l M D r e l g FIG. 6: Normalised effective long-time diffusion coefficient D rel from fitting of the simulations results, as function of thewall ruggedness parameter M . Parameters: gap opening pa-rameter g = 4 (black symbols) and g = 6 (blue symbols). In-set: D rel as function of g for M = 0 (black symbols), M = 10(blue symbols), and M = 50 (green symbols). ues of the wall ruggedness parameter, namely, M = 0, M = 10, and M = 50. Thus, we consider flat chan-nels, slightly corrugated, and heavily corrugated chan-nels. The corresponding results for the fitted values of thenormalised effective diffusion constant D rel are shown inthe inset of Fig. 6. For fixed value of g we see once morethat the diffusion is quickest when the wall is smootheror, i.e., when M is smaller. As was already observed inthe preceding paragraph a wider gap opening at a fixedvalue of M leads to faster diffusion. Thus, wider chan-nels can be traversed quicker.At first sight surprisingly, we observe that even forcompletely flat upper and lower channel walls ( M = 0,black line in the inset of Fig. 6) the diffusion in tighterchannels is slowed down in comparison to the situationin free space. This is not expected in systems with re-flecting boundaries, which are usually described with theFick-Jacobs equation. However, for the finite-size parti-cles studied here it is the result of the sticky boundaryconditions at the channel walls: in a tighter channel theparticle is more often close to the walls and binds tran-siently. Alternatively, the slow-down due to the interac-tion with the walls can be quantified by measuring theanomalous diffusion exponent α in the intermediate timeregime, on time scales 10 to 10 simulations steps. Thefitted values for α are depicted in Fig. 7 as function ofthe ruggedness M . The same trend as for the effectivediffusion coefficient is seen for the anomalous diffusionexponent α of the transient subdiffusive regime. For in-creasing contour lengths of the channel wall the motionis increasingly subdiffusive. As expected, the transientsubdiffusion is heavier for the tighter channel ( g = 4).The dependence of α on the gap opening parameteris shown in the inset of Fig. 7. In this case, only values α M α g FIG. 7: Anomalous diffusion exponent α as function of wallruggedness M from power-law fit of the EATA MSD data.Parameters: g = 4 (blue) and g = 6 (black). Inset: α asfunction of g for M = 10 (black), and M = 50 (blue). α D rel g=4g=6 FIG. 8: Anomalous diffusion exponent α of the transientlysubdiffusive regime as function of the long-time diffusion ex-ponent D rel with gap opening parameters g = 4 (black) and g = 6 (blue). obtained with corrugated channels ( M = 10 and M =50) are shown.[44] Again, the curves for α are similar tothe ones obtained for the effective diffusion coefficient: α is an increasing function of g and M = 10 leads to lesspronounced transient subdiffusion than M = 50. Thisanalogy motivates the study of the relation between α and D rel in more detail.In Fig. 8 we plot all α values for gap opening g = 4and g = 6 as function of the corresponding fitted values of D rel . The results show that there is a strong (nonlinear)correlation between both parameters. For increasing val-ues of D rel the value of α also increases, with decreasingslope. The relation between both parameters is bijec-tive, thus, both quantities are appropriate and sufficientto quantify the slow-down of the motion along the chan- f open M f open g FIG. 9: Fraction f open of unblocked configurations as func-tion of the channel ruggedness M for channel opening g = 4(black) and g = 6 (blue). Inset: f open as function of g for M = 10 (black) and M = 50 (blue). nel. For similar values of D rel the values for α obtainedwith the tighter channels ( g = 4) are slightly larger thanthose for the wider channels ( g = 6). However, this factshould not be overstated: all data sets were fitted in thetime interval 10 to 10 , irrespective of the exact shapesof the curves. It is conceivable that a closer connectionbetween the values for α could have been obtained bychoosing the fitting time window for each curve individ-ually.As mentioned above, in an ensemble of systems withcorrugated boundaries not all channels can be traversedcompletely. If a channel is blocked somewhere, the cor-responding squared displacement of the particle positionhas an upper limit. On an ensemble level these trajecto-ries will reduce the average values of α and D rel . Thus,it is important to extract from our simulations solely theunblocked configurations. The corresponding parameter f open is plotted as function of the ruggedness M for fixedgap openings g = 4 (black symbols) and g = 6 (blue sym-bols) in Fig. 9, and as function of the gap opening g (for M = 10 and M = 50) in the inset of Fig. 9. An in-spection of Fig. 9 shows that this fraction is a decreasingfunction of M and an increasing function of g . Over-all, the curves look similar to the ones of the normalisedeffective diffusion coefficient D rel (compare Fig. 6).To better understand why above a certain threshold ofthe boundary ruggedness parameter M the effective dif-fusion coefficient only decreases slightly (see Fig. 6), it isinstructive to study the average weighted effective width w wgt of the channels. Here, weighted means that for anyblocked channel the width is set to zero. The weightedwidth is plotted as function of M in Fig. 10. For in-creasing, yet small values of M the effective weightedchannel width decreases, until it reaches a shallow min-imum beyond which w wgt levels off to a plateau. Whileheavily rugged channel walls, in principle, can be tighter w w g t M FIG. 10: Weighted effective channel width w wgt as functionof channel ruggedness M for gap opening g = 6 (black lineand symbols) and g = 4 (blue line and symbols). or more spacious than a flat configuration with the samegap opening g , the tighter ones are in constant ‘danger’of precluding the particle passage. Thus, the averagewidth of those traversable channels is larger for morerugged wall configurations. This facilitates the trans-port through these configurations, as the sticky bound-aries are further away. However, as remarked earlier,more rugged walls slow down the diffusion due to the oc-currence of bulges and constrictions [20], such that wehave two opposing effects, which mostly (almost) canceleach other. Consequently, w wgt (and thus D rel ) reachesa plateau above a threshold value of M .Fig. 11 shows the normalised effective diffusion coeffi-cient D rel as function of the weighted channel width w wgt .For better visibility data points with identical M valuesare connected by lines ( M = 0: black line, M = 10:blue line, and M = 50: green line). While for a fixedvalue of M more spacious channels allow faster diffu-sion, the heavy scatter between values of D rel obtainedwith similar values of w wgt (grey symbols) shows thatthe knowledge of the mean channel width of an ensembleof channels alone is insufficient to predict the transportproperties. V. COMPARISON WITH DAGDUG ANDPINEDA’S FORMULA
In order to compare our results obtained from ensem-bles of channel configurations with the result of Dag-dug and Pineda, we make a simplifying assumption: forall unblocked configurations, we determine D DP ( x ) fromEq. 4 and subsequently D eff , DP from Eq. 5. For allblocked configurations the effective diffusivity D eff , DP =0 accounts for the fact that on long time-scales particlesin these configurations do not contribute significantly tothe MSD. Finally, we average over all configurations in D r e l w wgt FIG. 11: Normalised effective diffusion coefficient D rel asfunction 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). the ensemble and normalise through division by the dif-fusion coefficient in free space, D rel , DP = h D eff , DP i /D .We compare these values with our fitted values of the nor-malised diffusion coefficient in the upper panel of Fig. 12,where data points with the same value of the gap opening g are represented in the same colour.This is the central result of our study: As demon-strated in the upper panel of Fig. 12, for fixed values ofthe gap opening g there is a linear relation between D rel and D rel , DP . Most of the data points are located withina 10% confidence interval around Dagdug and Pineda’svalue. More explicitly, a linear fit of the data points ob-tained in a simulation series yields the following results.For g = 4, the fitted relation between the two is D rel =(0 . ± . . ± . D rel , DP , while for g = 6we find D rel = (0 . ± . . ± . D rel , DP .For g = 8, 10, and 12 we did not fit the data as therewere only three values available.The fact that the slope of the fits is somewhat below1 shows that Dagdug and Pineda’s formula, which onlyapplies to systems with perfectly reflecting boundariesslightly overestimates the diffusion coefficient comparedto our system with sticky boundary conditions. Thus, asexpected the additional interaction with the boundariesfurther slows down the diffusion, which is already reducedby the occurrence of entropic bottlenecks. This reason-ing is further substantiated by the observation that intighter channels (with g = 4), where these surface effectsplay a more important role, the slope of the conversionformula is smaller, and the deviation from Dagdug andPineda’s formula is more pronounced. Given the quiteintricate form of the effective channel width (see Figs. 1and 3) it is not feasible to quantify this effect analyti-cally. However, in the future other values of s and s could be considered to study the magnitude of the cor-rection terms numerically. D r e l D rel,DPg= 4g= 6g= 8g=10g=12 0 0.04 0.08 0.12 0.16 0 0.2 0.4 0.6 0.8 1 p D rel,DP,iM =98M =50M =20M = 6 FIG. 12: Upper panel: normalised effective diffusion coef-ficient D rel from our simulations in the sticky, corrugatedchannel as function of the value D rel , DP obtained from theresult (4) of Dagdug and Pineda. Black: gap opening param-eter g = 4, red: g = 6, green: g = 8, blue: g = 10, andcyan: g = 12. The two grey lines mark the range of ± g = 6 in all four cases. M = 6: black, M = 20:red, M = 50: green, M = 98: blue. That deviations of our results from Dagdug andPineda’s formula are also based on the fact that theiranalysis applies to one given channel configuration.Driven by the physical application we here consider anensemble of different channel wall configurations, solelydefined by the fixed macroscopic parameters g and M .Individual configurations may therefore differ consider- ably. This is illustrated in the lower panel of Fig. 12,where for all unblocked configurations the expected ef-fective diffusion coefficient was calculated with Dagdugand Pineda’s formula. For g = 6 and four different valuesof M we see that the distribution of values of D rel , DP isfar from narrow. The values increasingly scatter for morerugged conformations (higher values of M ). Thus, thesemi-quantitative agreement of our data with the theoret-ical prediction on an ensemble level is indeed remarkable. VI. CONCLUSION AND OUTLOOK
We studied the motion of finite-size particles throughrandomly corrugated channels with sticky walls, observ-ing transient anomalous diffusion of the particles in theirpassage of the channel. We also obtained the long-timediffusion coefficient for this motion on time scales overwhich normal Brownian diffusion is restored. The con-trol parameters in our study were the gap opening pa-rameter fixing the distance between the channel walls atthe entrance and exit of the channel, as well as the wallruggedness parameter setting the maximal variation ofthe channel wall configuration. We quantified the depen-dence of the anomalous diffusion exponent and long-timediffusion coefficient as function of the ruggedness and gapopening, and showed that both quantities are in fact cor-related. We especially analysed the blocked channels,which the particle cannot fully traverse. The long-timeeffective diffusion coefficient was shown to agree well withthe prediction for point-like particles in channels withreflecting boundary conditions by Dagdug and Pineda.Translated into the language of pathogens in hydrogels,whose motion we want to mimic in our numerical studyof an ensemble of ‘parallel’ channels with identical gapopening and ruggedness, our study provides statisticalinformation on how many channels are blocked for thepathogen passage.In future studies it might be interesting to replace ourminimal model of interactions with the constituents ofthe hydrogel with a more realistic model. In particular,the introduction of dynamic boundaries, for example byvarying the gap opening during the simulations might beworth considering in order to model the structural changeof the hydrogel due to the binding of pathogens.MB and AG acknowledge funding from the GermanFederal Ministry for Education and research, and RMfrom the Academy of Finland within the FiDiPro scheme. [1] D. Brockmann, L. Hufnagel, and T. Geisel,
Nature , 2006, , 462.[2] K. E. Jones et al.,
Nature , 2008, , 990.[3] A. J. Alanis,
Arch. Med. Res. , 2005, , 697.[4] R. S. Sellar and K. S. Peggs, Brit. J. Haematol. , 2012, , 559.[5] J. Holtz and S. Asher,
Nature , 1997, , 829.[6] T. Miyata, N. Asami, and T. Uragami,
Nature , 1999, , 766.[7] K. Gawel, D. Barriet, M. Sletmoen, and B. T. Stokke, Sensors , 2010, , 4381.[8] E. Wischerhoff, N. Badi, J.-F. Lutz, and A. Laschewsky, Soft Matt. , 2010, , 705.[9] 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, , 178101.[10] J-.H. Jeon, N. Leijnse, L. Oddershede, and R. Metzler, New J. Phys. , 2013, , 045011.[11] G. Seisenberger, M. U. Ried, T. Endreß, H. B¨uning, M.Hallek and C. Br¨auchle, Science , 2001, , 1929.[12] 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, , 048103.[13] 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. USA , 2013, , 4911.[14] R. Metzler and J. Klafter,
Phys. Rep. , 2000, , 1.[15] R. Metzler and J. Klafter,
J. Phys. A , 2004, , R161.[16] M. J. Saxton and K. Jacobson, Ann. Rev. Biophys.Biomol. Struct. , 1997, , 373.[17] E. Barkai, Y. Garini, and R. Metzler, Phys. Today , 2012, (8), 29.[18] F. H¨ofling and T. Franosch, Rep. Prog. Phys. , 2013, ,046602.[19] P. S. Burada, P. H¨anggi, F. Marchesoni, G. Schmid, andP. Talkner, Chemphyschem , 2009, , 45.[20] R. Zwanzig, J. Phys. Chem. , 1992, , 3926.[21] D. Reguera and J. M. Rub´ı, Phys. Rev. E , 2001, ,061106.[22] P. Kalinay and J. K. Percus, J. Chem. Phys. , 2005, ,204701.[23] P. Kalinay and J. K. Percus,
Phys. Rev. E , 2005, ,061203.[24] P. Kalinay and J. K. Percus, Phys. Rev. E , 2006, ,041203.[25] A. M. Berezhkovskii, M. A. Pustovoit, and S. M. Bezrukov, J. Chem. Phys. , 2007, , 134706.[26] P. Kalinay and J. K. Percus,
Phys. Rev. E , 2008, ,021103.[27] S. Martens, G. Schmid, L. Schimansky-Geier, and P.H¨anggi, Phys. Rev. E , 2011, , 051135.[28] A. Berezhkovskii and A. Szabo, J. Chem. Phys. , 2011, , 074108.[29] P. Kalinay,
Phys. Rev. E , 2013, , 032143.[30] S. Martens, A. V. Straube, G. Schmid, L. Schimansky-Geier, and P. H¨anggi, Phys. Rev. Lett. , 2013, ,010601.[31] A. G. Cherstvy, A. V. Chechkin, and R. Metzler,
New J.Phys. , 2013, , 083039.[32] A. G. Cherstvy and R. Metzler, Phys. Chem. Chem.Phys. , 2013, , 20220.[33] A. Fuli´nski, Phys. Rev. E , 2011, , 061140.[34] A. Fuli´nski, J. Chem. Phys. , 2013, , 021101.[35] R. M. Bradley,
Phys. Rev. E , 2009, , 061142.[36] L. Dagdug and I. Pineda, J. Chem. Phys. , 2012, ,024107.[37] S. Lifson and J. L. Jackson,
J. Chem. Phys. , 1962, ,2410.[38] R. Festa and E. D’Agliano, Physica A , 1978, , 229.[39] L. Dagdug, A. M. Berezhkovskii, Y. A. Makhnovskii, andV. Y. Zitserman, J. Chem. Phys. , 2008, , 184706.[40] W. Riefler, G. Schmid, P. S. Burada and P. H¨anggi,
J.Phys. Cond. Mat. , 2010, , 454109.[41] P. A. Netz and T. Dorfm¨uller, J. Chem. Phys. , 1995, ,9074.[42] M. Saxton,
Biophys. J. , 1996, , 1250.[43] V. Tejedor, O. B´enichou, R. Voituriez, R. Jungmann, F.Simmel, C. Selhuber-Unkel, L. B. Oddershede, and R.Metzler, Biophys. J. , 2010, , 1364.[44] In flat channels with M1