Particle invasion, survival, and non-ergodicity in 2D diffusion processes with space-dependent diffusivity

We study the thermal Markovian diffusion of tracer particles in a 2D medium with spatially-varying diffusivity $D(r)$, mimicking recently measured, heterogeneous maps of the apparent diffusion coefficient in biological cells. For this heterogeneous diffusion process (HDP) we analyse the mean squared displacement (MSD) of the tracer particles, the time averaged MSD, the spatial probability density function, and the first passage time dynamics from the cell boundary to the nucleus. Moreover we examine the non-ergodic properties of this process which are important for the correct physical interpretation of time averages of observables obtained from single particle tracking experiments. From extensive computer simulations of the 2D stochastic Langevin equation we present an in-depth study of this HDP. In particular, we find that the MSDs along the radial and azimuthal directions in a circular domain obey anomalous and Brownian scaling, respectively. We demonstrate that the time averaged MSD stays linear as a function of the lag time and the system thus reveals a weak ergodicity breaking. Our results will enable one to rationalise the diffusive motion of larger tracer particles such as viruses or submicron beads in biological cells.


I. INTRODUCTION
For a typical bacterial cell such as E. coli , various proteins, large cellular complexes, nucleic acids, lipids, etc. occupy some 30-40% of the cell volume [1][2][3][4].The implications of this macromolecular crowding on the characteristics of diffusing particles of various sizes are still under debate [5,6].Another source impeding the free diffusion of particles in eukaryotic cells stems from a network of cytoskeletal filaments and internal membranes like the endoplasmic reticulum or the nuclear membrane.Such forms of crowding impair the particle diffusivity inside a cell and may alter the law of diffusion altogether, from Brownian motion to a subdiffusive law.In the latter case, the mean squared displacement (MSD) scales as [7] x 2 (t) with the anomalous diffusion exponent 0 < β < 1. Experimental data are available, inter alia, for the in vivo subdiffusion of proteins [8] and enzymes [9], endogenous submicron particles (lipid and insulin granules) [10][11][12][13][14], viral particles [15], fluorescently labelled gold particles [16], messenger RNA molecules [17], as well as the telomeres of chromosomes [18].In vitro, dense solutions of coil-like polymers, proteins, or worm-like micelles often mimic the effects of molecular crowding which depend on the particle size, the solution viscosity, and the effective medium porosity [19][20][21][22].Similarly, in large scale computer simulations of crowded lipid membranes, subdiffusion is observed for various membrane chemistries [23][24][25].
Measuring the apparent local diffusivity of smaller proteins in bacterial [26] and eukaryotic [27] cells reveals a nontrivial dependence on the position in the cell.One reason for this spatial variation of the diffusivity may be the cells' geometrical shape [26].Thus, certain cell types possess a 'fried egg-shape' (Fig. 1) with a significant variation of the cell thickness from the periphery towards the nucleus.A higher apparent abundance of proteins in the cytosol near the nucleus, interpreted as a higher cytoplasm diffusivity, may simply originate due to the 2D imaging of the fully 3D particle trajectories.Away from the thicker perinuclear region, the cell periphery offers only a thin, nearly 2D domain for the particle diffusion.
Another source for the variations of the local diffusivity is the heterogeneity of the density of the macromolecular crowding in the cytoplasm and nucleoplasm, as well as of the dense cytoskeletal meshwork near the cell periphery, and the accumulation of large cellular organelles in a perinuclear region.How exactly this affects the porosity of the cytoplasm and the diffusivity of tracers of different sizes is not well established [28,29].Specifically, substantial deviations from the Stokes-Einstein law for protein tracers of varying molecular weights (MW) diffusing in the E. Coli cytoplasm were observed and the diffusivity shown to follow the scaling law D ∼ MW −0.7 [30].Small tracer proteins apparently experience a higher porosity near the nucleus of mammalian cells [27], while the diffusion of larger proteins becomes progressively restricted [30].Thus, from a biological perspective, a stochastic model with spatially-varying diffusivity may mimic the effects on the diffusion of tracer particles in the heterogeneous environment of the crowded cellular cytoplasm and will serve as an empirical description of secondary processes such as intracellular, diffusion-controlled reactions.We here study the physical properties of such a heterogenous diffusion process (HDP).
Important clues come from viral particles, a major class of natural diffusers in the bacterial cytoplasm.After internalisation by receptor-driven endocytosis [31], vir- uses often recruit highly processive cellular motor proteins [32,33] which ensure fast and efficient viral transport between the cell periphery and the nucleus, where the viral replication and assembly often occurs [34,35].The intra-cellular dynamics of viruses and their multistep infection pathways, as monitored by single-particle tracking, exhibits some features of anomalous diffusion [15].For instance, the scaling exponent β of the viral motion is shown to depend on the region of the cytoplasm in which the diffusion takes place.
Three different modes of transport for adenoassociated viruses [36,37] were identified in living substrate-adhered HeLa cells [15,38,39].The first is Brownian motion with β=1 albeit with a much smaller diffusion coefficient than in dilute aqueous solution.The second mode is that of subdiffusion with β = 0.5 . . .0.9 and with a broad apparent distribution of diffusivities.The third mode is that of motor-driven transport of viruses via quasi-1D persistent walks along microtubular filaments, mediated by molecular motors driven by energy from ATP conversion.Upon infection, the ballistic, driven motion with β=2 yields an effective drift of virions towards the nucleus.These diffusion-based and active modes of viral transport can interchange.Although the fraction of actively-transported virions is relatively small [15], this 'active pathway' is often vital for a successful viral infection.Small viruses can reach the nucleus solely by thermal diffusion, while larger virions have no chance but rely on the active transport mechanism.Indeed, an accumulation of viruses near the nucleus was shown to be inhibited by microtubuli-depolymerising drugs (e.g., nocodazole) that suppress motor-assisted virus transport [40].
Several recent models of intermittent transport [41,42] were implemented to describe kinetics of viral infection, and search optimisation models with 2D versus 3D intermittent dynamics were developed [43,44].A number of diffusion [45], diffusion-reaction-advection [40], and kinetic transport [46] models were suggested to rationalise the features of intracellular virus trafficking.In particular, the kinetics of spreading of a viral population starting at the cell membrane and the accompanying nucleus invasion times were computed [46] and compared to typical time scales of viral infection recorded experimentally [47,48].In a series of theoretical and computer simulation studies Holcman and colleagues [49][50][51][52][53] modelled the process of viral trafficking as a sequence of alternating Brownian 2D diffusive excursions and ballistic motorpowered propulsions along radially-ordered microtubuli filaments.The dynamical characteristics of viral invasion were computed in such a 2D planar pie-like model.The probability density function (PDF) and the mean time of nucleus invasion by viruses were evaluated [51].
More advanced theoretical models can also include a rate of viral degradation in the cytoplasm [46], the kinetics of viral binding to microtubuli, and some bi-directionality of virus transport by the motors.
Here we consider the passive diffusion of tracer particles of sizes comparable to a virus capsid in a model cell.To construct our model we include the following information.From the viral trajectories reported in Refs.[15,38] we conclude that those particles exhibiting normal diffusion with β = 1 take azimuthal journeys, at about constant separation from the cell nucleus.Such propagation likely takes place in a region of roughly constant diffusivity.In contrast, that part of the viral population that diffuses anomalously mainly travels in the radial direction.We propose below that heterogeneities of the medium during the journey of a particle from the cell membrane to the nucleus gives rise to anomalous (in particular, sub-diffusive) features for this second population of particles.
Recently, extending previous studies [54,55] we examined effects of position-dependent diffusivities on the ensemble and time averaged characteristics for 1D HDPs [56,57].We tested several functional forms for D(x) (power-law, logarithmic, and exponential) to rationalise their implications onto diffusive and ergodic properties of the process.For power-law forms of D(x) we predicted from stochastic simulations and analytical calculations the regimes of sub-and super-diffusive behaviour.The conditions for weak ergodicity breaking were also analysed in details, an important feature when information from single particle tracking studies is evaluated in terms of time averages [6].The diffusion is non-ergodic in 1D due to the heterogeneities of the medium.In general, despite a non-Brownian scaling of the MSD, the time averaged MSD was shown to follow a strictly linear growth with lag time.These features are similar to those for continuous time random walk processes [58].
For the 2D HPDs examined below we also find nonergodic behaviour.In addition, we compute a number of biologically relevant quantities such as the survival probability S(t) of particles in a circular domain for both diffusion from the inside of the cell to the outside and vice versa, the first-passage time dynamics for reaching the domain boundary, the PDF for the spreading of diffusing particles starting at the cell centre, at the cell boundary, and for initially uniformly-distributed walkers.
This paper is organised as follows.In Section II we introduce the basic notations and the main quantities to be analysed.We outline the numerical scheme used in computations as well as the implemented theoretical concepts.In Section III we report the main simulations results and support them by asymptotic analytic calculations.We analyse the effects of the system heterogeneity and polar cell geometry on diffusive, kinetic, and ergodic properties of the HDP.In Section IV the Conclusions are drawn and possible applications of the results are discussed.

II. MODEL
Our model cell is a circular disc with a reflecting outer boundary, mimicking the situation that internalised viruses do not leave the cell again.We analyse the following form of the diffusivity that is solely dependent on the radius r away from the cell centre.At small r values the diffusion is fastest, continuously slowing down towards the outer cell region (the 'cell membrane').To avoid divergencies in the discrete simulations scheme implemented below, we regularised D(r) in Eq. ( 2) by introduction of the constant A > 0. At r A, the diffusivity exhibits the power-law scaling D(r) ∼ 1/r 2 .The constant D 0 fixes the units of the diffusivity.The dependence (2) is in qualitative agreement with the experimentally measured trends for the diffusivity of small fluorescently-labelled proteins in the cytoplasm of mammalian NLFK and HeLa cells [27].It also reflects the above observation that azimuthal diffusion is fully Brownian, i.e., in our language, the diffusivity remains constant.The simulations method described below is readily applicable to other D(r) forms.
We characterise the HDP in terms of the ensembleaveraged MSD of particles defined via the PDF P (r, t), r 2 (t) = r 2 P (r, t)2πrdr. ( For a 2D trajectory r(t) = {x(t), y(t)} (r = √ r 2 ) of length T , the time averaged MSD is defined as the sliding average with the lag time ∆, While for an ergodic process for sufficiently long measurement times T the equivalence r 2 (∆) = δ 2 (∆) holds, the behaviour of the two quantities remains different even for T → ∞ in weakly non-ergodic systems [6,58,59].In particular, individual realisations of time averaged quantities becomes irreproducible [6,58,59].The ergodicity breaking parameter EB characterises the deviation of the system from the ergodic behaviour.It contains the second moment of the time averaged MSD and is defined as follows [60,61] EB(∆) = lim For the canonical Brownian motion in 1D (d = 1) one obtains [60] EB This means that EB BM → 0 at ∆/T → 0, and the spread of time averaged MSD traces around the mean computed over N traces, approaches a sharp δ-function shape, i.e., the experiment is fully reproducible [6,58,59].To extract a statistically meaningful spread of δ 2 (∆) values around the mean δ 2 , the condition ∆/T 1 should be satisfied.
We also define the survival probability S(t) of particles in the circular domain when either of the boundaries is considered absorbing, and the particles are released at the opposite boundary.The probability of particles in this scenario is not conserved and S(t) tends to zero as time progresses.The PDF of first passage is then defined as −dS(t)/dt, and the mean first passage time as MFPT = ∞ 0 S(t)dt.We evaluate the statistics of the first arrival times directly from the generated trajectories, r(t).
At every time step in the computer simulations we use the Klimontovich-Hänggi [62] post-point scheme to evaluate the two coupled Langevin equations with independent noise sources, Here, the increments of the Wiener processes for the corresponding coordinate, (W x,i+1 − W x,i ) and (W y,i+1 − W y,i ), each represent a different δ-correlated Gaussian noise with unit variance.Unit time intervals δt separate consecutive iteration steps in the simulations.From N 2D stochastic trajectories {x(t), y(t)} generated for the initial particle position, x(t = 0) = x 0 and y(t = 0) = y 0 , the ensemble and time averaged characteristics of the HDP are evaluated.We note that we could also use the Stratonovich scheme to simulate the process.For the MSD, similar to the 1D case [56], the difference between the two representations occurs only in the prefactor and is of order unity.

A. MSD and time averaged MSD
The computed ensemble and time averaged MSD as well as the mean time averaged MSD are shown in Fig. 2. For the 1D case, the MSD for a diffusivity of the form This asymptote, derived within the Stratonovich scheme in Ref. [56] and shown as the black dashed curves in Fig. 2, is in good agreement with our 2D simulations for the D(r) defined in Eq. ( 2).At proximate initial positions x 0 and y 0 , the deviations from the theoretical MSD asymptote (9) almost vanish after several simulation steps.For more distant initial positions {x 0 , y 0 }, the sub-linear r 2 (t) t 1/2 -scaling is approached somewhat later, giving rise to an initial plateau.
The time averaged MSD trajectories are linear functions of the lag time ∆, their mean scaling as see Fig. 2. The spread (amplitude scatter) of the time averaged MSD traces is very pronounced, with large trajectory-to-trajectory variations, see the red traces in Fig. 2.This indicates an ergodic violation, see below.Also, at shorter T values the spread of individual δ 2 (∆) in the region of ∆/T 1 progressively decreases for larger values of the particle initial position (not shown).Here, we do not quantify the details of the distribution φ(δ 2 / δ 2 ) of individual traces δ 2 .We refer the reader to Refs.[56,57] where this procedure is discussed in detail for the 1D HDP.
We checked that for D(x) = const we obtain the standard 2D result with with only a minute scatter of δ 2 (∆) traces at ∆/T 1.Moreover, we find that indicative of the self-averaging behaviour typical for Brownian motion.
For the choice of D(r) used here, leading to subdiffusion, usually the ratio δ 2 (∆) / r 2 (∆) 1 for not too small values of the initial positions {x 0 , y 0 }, see Fig. 2. As a connection to experiments, let us define the model parameters that can describe the MSD magnitudes measured in the tracking experiments of small adenoassociated viruses [15,38,39], as mentioned in the Introduction.Specifically, for the subdiffusive population of viruses the MSD measured in the cells after the diffusion time of t ≈ 0.32s was ≈ 0.4µm 2 , see Fig. 3G in Ref. [15].The viral diffusivity was D ∼ 0.2µm 2 s −0.6 for their subdiffusive motion with exponent β ≈ 0.6.To get the same MSD value in the same physical time t the diffusion coefficient D 0 in our model would be D 0 ≈ 10µm 4 s −1 .9) and (11).The MSD corresponds to the blue line.Parameters are the same as in Fig. 2, except for x0 = y0 = 1 and N = 40.

B. Azimuthal and radial diffusion
We project the increments of the diffusing particles at each simulation step i with particle position r i onto the radial and azimuthal directions and compute the single-step displacements δr i and r i δΦ i .We account for the clock-and anti-clockwise azimuthal rotation of the particle position vector r i .We then restore the corresponding average displacements after t = T /δt simulation steps, computed as the average over all the traces, ρ 2 (t) = ( t i=1 δr i ) 2 and Φ 2 (t) = ( t i=1 r i δΦ i ) 2 .The results of the simulations show that the growth of the radial increments, similarly to the MSD in Eq. ( 9), obeys the subdiffusive law For the azimuthal increments, in contrast, the diffusion is Brownian, Fig. 3, with the scaling

C. Ergodic violation
The simulations show that the ergodicity breaking parameter for short lag times ∆ assumes values close to those for the 1D case with analogous D(x) treated in Ref. [56].The EB values at ∆/T 1 deviate from zero, indicating a weak ergodicity breaking and nonequivalence of ensemble and time averaging for this 2D diffusion process in a heterogeneous environment.Nonhomogeneities in the diffusion coefficient break the ergodicity in the system, see also the discussion in Ref. [63].
For long lag times, when ∆ → T , and for {x 0 , y 0 } values in the high-diffusivity region close to r = 0, the ergodicity breaking parameter approaches 1/2 of the value (6) of the asymptote for 1D Brownian motion obtained in Ref. [60].Such a reciprocal dependence on the space dimension d, has recently also been discovered for multi-dimensional fractional Brownian motion [64].
Clearly, in our system with an inhomogeneous diffusivity, the initial conditions of the diffusing particles affect the magnitude of the time averaged MSD traces and thus the values of the ergodicity breaking parameter.Specifically, as the values of x 0 and y 0 decrease, the value of EB(∆ → 1) decreases, as shown in Fig. 4. The HDP thus becomes more ergodic, and the Brownian asymptote ( 12) is approached at earlier lag times ∆.This is due to the fact that at larger {x 0 , y 0 } the spatial heterogeneities are sampled by the considerably slower walkers to a lesser extent for the same length T .
Note here that the evaluation of the dependence EB(∆) often requires much better statistics than that needed for the MSDs presented in Fig. 2. The reason is the large spread of δ 2 (∆) between trajectories at all lag times ∆, see the red curves in Fig. 2.This scatter has more severe implications on the ergodicity breaking parameter containing the square δ 2 2 , see Eq. ( 5), and involving the averaging over N traces.The dependence of the initial value EB(∆ = 1) on the trajectory length T for different initial conditions is illustrated in Fig. 5. Similar to the 1D situation treated in Refs.[56,57], the variation of EB(T ) depends on how far the system is away from the ergodic state for the imposed initial conditions.For instance, the ergodicity breaking parameter EB(∆ = 1) for x 0 = y 0 = 3 and short traces with T = 10 2 is quite close to the Brownian value given by Eq. ( 12).Conversely, for the same initial conditions but longer trajectories, T = 10 3...5 , the system   is more non-ergodic and the corresponding EB parameter is larger than that for x 0 = y 0 = 0.1 or 1 (Fig. 5).This is the reason why we observe intersection of curves for different {x 0 , y 0 } values shown in Fig. 5.At T → ∞ the ergodicity breaking parameter tends to a universal value.
The data show that system heterogeneities indeed cause a weak ergodicity breaking in the 2D HDP with diffusivity (2).Due to the non-equivalence of the radial and azimuthal diffusion, we predict a direction-dependent ergodicity breaking parameter, see Fig. 6.We observe that the radial azimuthal ergodicity breaking parameters EB ρ and EB Φ become quite close in the limit ∆/T 1.In the limit of long lag times, as ∆ ∼ T , the parameter EB ρ behaves similarly to EB computed from r(t), compare Figs. 4 and 6.In contrast, the azimuthal parameter EB Φ does not approach the Brownian asymptote ( 12) at later stages of the time averaged trajectories.

D. PDF and spreading of particles
The spreading of particles starting at the cell boundary at r = R is characterised by the PDF shown in Fig. 7.The accumulation of particles near the reflecting outer wall in the region of low diffusivity contributes to the enhanced azimuthal spreading.This spreading remains profound also at later times, because of the Brownian diffusive behaviour in azimuthal direction, as contrasted to subdiffusive spreading in the radial direction, see Fig. 3.The overall trend is similar to the 1D case [56,57], where at long times the particles tend to accumulate in the regions of lower diffusivity.Naturally, the average effective jump length of particles diffusing near r = 0 is larger than in the region of slow diffusion near the cell boundary.As one can see from Fig. 7, a strong azimuthal spread at t = T /30 . . .T /10 turns into a profound invasion of particles over the entire cell at t = T (for trace length T = 10 5 and N = 150 analsed trajectories in this figure).The time evolution of the radial PDF shown in Fig. 8 quantifies the 2D plots in Fig. 7 when particles are initially released at the fringe of the cell.We observe that for longer trajectories the maximum of the PDF, initially localised at r = R, progressively spreads and approaches a universal scaling law given by p(r) r.
We also simulated the diffusive 'focusing' of walkers, that were initially homogeneously distributed in the cell.We find that fast-diffusing particles leave the region near the origin at r = 0 relatively quickly and progressively shift the maximum of the PDF towards the region of slow diffusion near the cell periphery, see the graphs for different times t in Fig. 9.This trend is similar to the 1D situation with power-law diffusivity [56].

E. Survival probability: diffusion from the nucleus to the membrane
After starting at the cell centre and diffusion towards the cell membrane at r = R, the probability of staying in the domain of radius R is described by the survival probability S(t) shown in Fig. 10.The simulations results obey the universal scaling Naturally, for larger cells the diffusing particles start to follow this asymptote at later times, as it takes longer to reach the outer cell border by diffusion.Relatively strong variations of S(t) at later times are due to back-andforce diffusion of individual particles through the outer boundary (which was treated permeable in the algorithm for computing S(t) in Fig. 10).We have checked that the scaling law ( 16) is valid also for other initial positions {x 0 , y 0 } in the cell (results not shown).We also expect Eq. ( 16) to remain valid for other choices of the diffusivity variation D(r) in the cell.
We also computed the distribution of arrival times of particles diffusing from the cell centre to the cell boundary, see Fig. 11.These distributions p(t arr ) reveal a wide spread, particularly at large R values, indicating large trajectory-to-trajectory fluctuations.From these distributions we determine the threshold time t 1/2 at which 50% of the fastest particles reach the outer cell boundary.Such a threshold characteristic is often important for biological problems, e.g., in the dynamics of population spreading or proliferation of viral infections.
The function t 1/2 (R) obtained via the analysis of the histograms presented in Fig. 11 often turns out to be bounded by two asymptotes.The first one is defined via the slowest diffusivity at the cell boundary r = R. Namely, from elementary scaling arguments we can write The second characteristic time scale is defined via the average diffusion coefficient in the domain, namely These asymptotes (respectively, the black and green lines in Fig. 12) indicate the leading-order scaling t 1/2 ∼ R 4 (apart from the logarithmic correction).The time t 1/2 characterises the arrival of the fastest half of a population of diffusing walkers, and it can be related to the effectiveness and reliability of the target search in such a heterogeneous medium.It is particularly important as the arrival time distributions are skewed [66], compare Fig. 11.Consequently, the mean of the distribution and its width are not the best indicators of the arrival statistics [66], and instead t 1/2 should be used.In our 2D bounded domain the first-passage time dynamics and the histograms for p(t arr ) can be fitted, e.g., by a generalised Gaussian distribution [66,67] (not shown).Clearly, for a larger domain size R the width of the distributions of arrival times grows, because of accumulated statistical fluctuations among diffusing particles with longer trajectories, see Figs. 11 and 14.

F. Survival probability: diffusion from the cell membrane to the nucleus
The survival probability for the diffusion from the outer boundary to the cell 'nucleus' exhibits the exponential scaling  in contrast to the t −1/2 law (16) for the opposite direction in the same domain.This exponential scaling is akin to the standard problem of 2D diffusion in a circular domain with a sink [65], see also Ref. [66] for the exponential scaling of S(t) in the 3D case.
The characteristic time t of the decrease of S(t) with time corresponds to the time the particles spend diffusing from the outer to the inner boundary in a medium with average diffusivity.As the radial diffusion is quasi-1D we can write For not too large (R−a) values, when the medium diffusivity varies only moderately in the concentric shell, such an ansatz for t works quite well.These asymptotes are shown as dashed lines in Fig. 13, in comparison to the simulation results for S(t).
Similarly to the results for nucleus-to-membrane diffusion in Fig. 11, we evaluate the distribution of the arrival times from the cell periphery to the nucleus of different sizes, compare Fig. 14.

IV. DISCUSSION AND OUTLOOK
We studied the diffusion of particles in a 2D circular domain with a radially varying diffusivity D(r).We showed that the resulting HDP is weakly non-ergodic in the sense that time and ensemble averages of physical quantities such as the MSD behave differently.This ef-fect was shown to depend on the initial conditions of the diffusive walkers.The diffusion in the direction of the diffusivity gradient was shown to be anomalous, while the azimuthal diffusion occurs in a nearly constant environment and is Brownian.This behaviour is reminiscent of the radial and azimuthal diffusion of viral particles monitored in the bacterial cytoplasm, with purely radially varying diffusivity [15].
Specifically for the evaluation of single particle tracking data, our results for the non-ergodicity imply that (i) the time averages of physical quantities such as the MSD behave differently from their ensemble analogues, and that (ii) individual time averages are not reproducible, i.e., there occurs a major scatter in the amplitudes of these quantities.Both need to be taken into account for a proper physical interpretation of data.
We demonstrated that the diffusion from the domain centre to its boundary (nucleus to membrane) and the reverse process obey entirely different behaviours for the respective survival probabilities.Namely, the S(t) t −1/2 scaling law was found for nucleus-membrane diffusion and the exponential S(t) e −t/t decay was identified for membrane-nucleus diffusion.This latter fact as well as the spreading of particles according to these two scenarios can be rationalised in terms of a domain-averaged diffusion coefficient.
A quantitative understanding and the ability to tune viral diffusion in living cells has enormous potential as a tool to control and hopefully suppress the proliferation of infection.Viral gene delivery carriers [37,68] with a high transfection efficiency actively transported by motors [69] are nowadays extensively used for gene delivery purposes.We note that our model may also be applied to macroscopic systems.Thus, the spatial spreading of epidemics in a population of animals subject to non-homogeneous habital or foraging conditions is another possible area for application for our model.
In the present paper, we focused on the statistical and nonergodic properties of HDPs in circular domains.A mathematical investigation of the process of viral infection in the presence of three inter-connected diffusion pathways (anomalous diffusion, normal diffusion, and active directional transport) is currently under way [70].

Figure 1 :
Figure 1: The variation of local diffusivity in the cytoplasm of a mammalian cell.The FRAP intensity of the cyan colour refers to the local effective porosity of the cytoplasm, scaling with the volume fraction available for protein diffusion.Scale bar is 10 µm.The image is taken from Ref. [27]; courtesy to Jörg Langowski.

2 Figure 2 :
Figure 2: Dependence of the ensemble-averaged MSD (thick blue curve), the mean time averaged MSD (thick blue curve), and the time averages of the MSD for individual trajectories (red curves) on time t or the lag time ∆.The theoretical asymptote (9) for the ensemble-averaged MSD is represented by the dashed black line.Parameters: A = 0.01, the starting positions are x0 = y0 = 0.1, 1, and 3 from the graphs from top to bottom.The number of traces for the averaging is N = 300, the length of each trajectories being T = 10 5 in units of the simulations time step.The simulation time for each choice of the starting conditions is ∼2.5 days on a standard 3 GHz working station.

Figure 3 :
Figure3: Anomalous behaviour of radial increments (orange) and Brownian diffusion of azimuthal particle increments (green line).The dashed asymptotes are given by Eqs.(9) and(11).The MSD corresponds to the blue line.Parameters are the same as in Fig.2, except for x0 = y0 = 1 and N = 40.

Figure 4 :
Figure 4: Ergodicity breaking parameter as a function of lag time ∆ for different initial conditions.The Brownian asymptote (12) for the 2D case is shown by the dashed line.Parameter are the same as in Fig. 2.

1 Figure 5 :
Figure 5: Dependence of EB(∆ = 1) on the trace length T .At least N =300 trajectories were used to compute each point in the graph.Parameters are the same as in Fig. 2.

Figure 6 :
Figure 6: Direction-dependent ergodicity breaking parameters.The Brownian asymptote (12) is represented by the dashed line.Parameter are the same as in Fig. 2. The colour coding corresponds to Fig. 3.

Figure 7 :
Figure 7: Series of PDFs in our 2D 'cell' for the temporal spreading of random walkers starting at the cell boundary at x0 = y0 = R/ √ 2. N = 150 trajectories of T = 10 5 time steps were analysed.The cell radius is R = 5, and the times t of the snapshots are indicated in the panels.The dark spot in the centre of each graph is due to the faster diffusion at r = 0 and a finite grid for sampling and projecting r(t) traces.

Figure 8 :Figure 9 :Figure 10 :
Figure 8: Dependence of the radial distribution function p(r) for particle invasion into a circular nucleus-free domain 0 < r < R starting from the cell boundary.Parameters are the same as in Fig. 7.

Figure 11 :
Figure11: Distributions of arrival times to the cell boundary for diffusion of particles initially released in the cell centre, plotted for varying cell radius R. In the inset we demonstrate the scaling t −3/2 expected from the survival probability in Fig.10.The colour scheme and parameters are the same as in Fig.10.

Figure 12 :
Figure 12: Scaling of t 1/2 for diffusion from centre to cell boundary with cell size R.The asymptotes correspond to Eqs. (17) (top) and (19) (bottom line).Parameters are the same as in Fig. 11.

Figure 13 :Figure 14 :
Figure 13: Exponential decay of the survival probability for the diffusion of particles starting at the cell membrane.The dashed lines represent Eq. (20).The radii a of the inner absorbing boundary are indicated in the graph, other parameters are the same as in Fig.10.