Non-linear continuous time random walk models

A standard assumption of continuous time random walk (CTRW) processes is that there are no interactions between the random walkers, such that we obtain the celebrated linear fractional equation either for the probability density function of the walker at a certain position and time, or the mean number of walkers. The question arises how one can extend this equation to the non-linear case, where the random walkers interact. The aim of this work is to take into account this interaction under a mean-field approximation where the statistical properties of the random walker depend on the mean number of walkers. The implementation of these non-linear effects within the CTRW integral equations or fractional equations poses difficulties, leading to the alternative methodology we present in this work. We are concerned with non-linear effects which may either inhibit anomalous effects or induce them where they otherwise would not arise. Inhibition of these effects corresponds to a decrease in the waiting times of the random walkers, be this due to overcrowding, competition between walkers or an inherent carrying capacity of the system. Conversely, induced anomalous effects present longer waiting times and are consistent with symbiotic, collaborative or social walkers, or indirect pinpointing of favourable regions by their attractiveness.


Introduction
An essential assumption for a continuous time random walk (CTRW) model is that the walkers do not interact with each other. This is a premise which in the application to e.g., biological systems may obscure important features. It is well-established in biology that the walkers "sense" the walker population indirectly via a chemotactic substance [1][2][3][4][5][6][7][8] or directly [7][8][9][10][11][12] as a result of interactions or volume exclusion processes [13][14][15]. Inspired by the demands of such model features, we seek to understand the repercussions of including these interactions with the mean population into the already known equations for a CTRW. In the classical formulation where the walkers do not interact, a CTRW is described by two integral equations P (x, t) = Ψ (t)P (x, 0) + t 0 j(x, τ )Ψ (t − τ )dτ, j(x, t) = ψ(t) R 3 dyλ(x − y)P (y, 0) where P (x, t) is the probability density function (PDF) of being at position x ∈ R 3 at time t > 0. Further, j(x, t) is the PDF of having just arrived, ψ(t) is the PDF of the waiting times between jumps and λ(x) is the jump density [16,17]. Consequently, Ψ (t) = ∞ t ψ(u)du is the survival function (probability of no jump occurring during [0, t]) [18]. Note that j(x, t) includes both the recent arrivals at t > 0, as well as the contribution from walkers starting from their initial (t = 0) spatial distribution P (x, 0). Starting from this initial spatial distribution, the position of the random walker is denoted by the random variable X(t).
It is well-known that from these equations we obtain the famous Montroll-Weiss equation which is given in Fourier-Laplace space by where ψ(s) = L{ψ(t)} = ∞ 0 ψ(t)e −st dt is the Laplace transform of the PDF of the waiting times between jumps, and λ(k) = F{λ(x)} = R 3 λ(x)e ik·x dx is the characteristic function of the jump density. Anomalous walks with a heavy-tailed waiting time PDF for large τ follow ψ(τ ) ∼ 1/τ 1+µ with 0 < µ < 1. Thus, in the limit when s → 0 we obtain the known asymptotic approximation ψ(s) = 1 − (sτ 0 ) µ Γ (1 − µ) + O(s) where τ 0 > 0. If the anomalous random movement through space is symmetric such that λ(k) = 1 − σ 2 k 2 + O(k 4 ) (where σ 2 is the variance of the jumps and k 2 = k · k), in the long-time limit we find from (2) the fractional diffusion equation where K µ is the generalised diffusion coefficient, 0 D 1−µ t the Riemann-Liouville derivative [19] and ∇ 2 x the Laplacian operator. This equation describes subdiffusion which is characterised by a mean squared displacement (MSD) X 2 (t) ∼ t µ . Equations (2) and (3) at the time of their introduction constituted breakthroughs in the field, by the argument that the key variable of interest is the waiting time T i the walker spent in a given state i (or location x) before jumping elsewhere. This waiting time could have a heavy tail in its distribution, thus allowing for anomalous effects [16,20]. We now introduce a dependence on the mean walker population which intuitively will either serve to increase or decrease waiting times between jumps.
We term a negative response one wherein a local increase in the mean number of walkers decreases the waiting times. For power-law distributed waiting times this has the potential to temper the anomalous effects. Negative responses arise naturally in systems with e.g. carrying capacities or Allee effects which lead to overcrowding, or competition between walkers [7,9,12]. Conversely, positive responses are ones wherein a local increase in the mean number of walkers increases the waiting times. Hence, systems with shorter waiting times could acquire anomalous behaviour for a sufficiently large increase in the waiting times. Positive responses are consistent with walkers which exhibit social or collaborative behaviour, or are indirect indications of favourable regions as seen by the number of walkers they have already attracted. A positive response need not occur immediately: there may be a small local increase initially before the walkers self-organise to prefer this region.
In nature one has observed both positive and negative responses, corresponding to mutually favourable coexistence or social behaviour (increased waiting times) and hostility or competition (decreased waiting times), respectively. Positive responses (which favour accumulation) have been observed for a variety of bird species, red deer, marmots and squirrels, while negative reponses (which favour dispersals of large groups) have been observed for certain gulls, roe deer, voles and badgers [11]. Reasons for these positive or negative responses depend on the environment the animals move in, the social tendencies they exhibit and the predation to which they are susceptible.
The key challenge of this paper is implementing nonlinear positive and negative responses to the linear CTRW model. If the statistical properties of the random walker depend on the local average population, what is the implication of this dependence on (1)-(3)? Naturally, the waiting time is then affected by the mean number of walkers, but how can these equations be extended to include this effect? Attempts at including non-linear effects into the equations have previously been made in [21][22][23][24][25][26][27]. In the classical case of Markovian transport equations, one can easily extend the description to include non-linear effects. This is done by letting the jump rate depend on the mean number of walkers [12,28]. However, anomalous transport is non-Markovian and the answer to this question becomes less trivial. While elegant in their simplicity, we do not believe equations (1)-(3) constitute the best starting point if one is to consistently introduce these biologically motivated non-linearities into the transport. Instead, we propose another method by which to achieve our goal.
The chosen method by which these walker interactions are introduced is by considering the escape rates T i at which they leave their current position/state i. If a walker leaves its current position i at a constant rate T i = λ, we recover an exponential waiting time PDF ψ i (τ ) = λe −λτ . However, an escape rate T i (τ ) = µi τ +τ0 for constant τ 0 , µ i > 0 results in a power-law PDF ψ i (τ ) = For those readers so inclined this can be interpreted as a semi-Markov process where the escape rate depends on the age.
One can now introduce a dependence of the waiting times on both the age τ and the mean local number of walkers, ρ i (t) at a position i via the escape rate T i = T i (τ, ρ i (t)). We shall show that this dependence on ρ i (t) can both induce and temper the anomalous effects.
In order to describe a negative response the escape rate T i (τ, ρ i (t)) must increase with the number of walkers such that ∂Ti ∂ρi > 0. This is intuitively sound: if a region is above the carrying capacity or overcrowded, the escape rate of the walkers will increase. Conversely, a positive response must decrease T i (τ, ρ i (t)) with the number of walkers such that ∂Ti ∂ρ < 0. This is, the larger the group (and thereby the more favourable the region), the lower the escape rate.
In Section 2 we introduce the key mathematical notation and our approach to the inclusion of non-linearities. Section 3 illustrates the effects of tempered and induced anomalous effects on a simplified two-state system. Finally, Section 4 derives a new tempered fractional differential equation for the distributed case of walkers on a lattice. Section 5 concludes on our findings.

Structural density approach
Let us consider a random walk on a regular lattice of spacing 1. Each possible position on the lattice thus takes a value R = i for every site i ∈ Z, and has an associated escape rate T i . The position of a walker at a certain time with respect to an arbitrary reference point is given by X(t) ∈ i. Hence, X(0) = 2 denotes a walker present on the lattice at position 2 at t = 0. The mean-field interaction between the walkers is applied locally to each site i.
The walkers will wait for a certain time (or age) T i at a site before jumping elsewhere. We let this waiting time be a random variable which is reset whenever the walker leaves the site, and is thus independent of previous visits. The jumps between sites happen with rates T i (τ, ρ i (t)), which may differ for each site and depend on both the age spent at the site (τ ) as well as the current mean number of walkers at the site (ρ i (t)). As a result of the latter dependence, the escape rate depends implicitly on time. The key factor which constrains our analysis is the dependence on ρ i (t). Because of this, we cannot ab initio introduce a waiting time PDF with this dependence, but instead use the escape rate from each site on the lattice.
As a consequence of our choice of waiting time T i at a site as our random variable, we consider the average density of walkers at a certain time with a certain residence time (age). We call these the mean structural densities n i (t, τ ). That is, n i (t, τ )∆τ is a measure of the mean walker population and gives the mean number of walkers at site i at time t with ages in the interval (τ, τ + ∆τ ). The total rate of change of n i (t, τ ) must balance with the escape rates from each site, such that we obtain the mean-field equation: This approach has previously been employed to study non-Markovian transport [8,29] and is due to Cox and Miller [30]. Changes in each site are hence purely a result of the walkers which leave. For simplicity, we assume the standard initial condition at t = 0 where all walkers across the lattice have zero age. That is, they are all recent arrivals at their sites, such that where ρ 0 i is the initial spatial distribution of the walkers across the sites. This assumption does not affect the longterm dynamics on the lattice. Note that the dependence T i (τ, ρ i (t)) is purely with other walkers in the same site, such that we have zero-range interactions. Then, as n i is the mean number of walkers at site i with a certain age, the total mean number of walkers at the site is simply the sum over these ages, which yields However, our methodology is flexible and can be extended to account for finite range interactions beyond walkers at the same site. To this end, we consider a mean neighbourhood of walkersρ i (t) defined bȳ where the kernel ζ(i|j) governs the interactions of walkers at site j with walkers at site i; a standard approach when considering finite range interactions [31]. A process interacting with walkers at other sites would thus have an escape rate T i (τ,ρ i (t)). We must also include the boundary conditions which describe what occurs to a walker when it jumps from one site to another. This process is governed by the redistribution kernel κ(j|i), the probability of entering site j if leaving site i. Hence, new arrivals at a certain site are apportioned from those which left other sites immediately before: where κ(j|i) is normalised. This formulation of the transport has the advantage that it makes no assumptions about the particular forms of the escape rates T i , which may include whatever non-linearities are desired. If T i = T i (τ ) and is not a function of ρ i (t), then the above methodology is equivalent to the results of (1) if one lets i → x and κ(i|j) → λ(x − y). The key difference between the two methods is that (1) describes probabilities of a walker being at a certain point in time and space, whereas our method concerns the mean number of walkers. In Section 3 we apply our approach to non-linear interactions occurring between two states, which may be thought of as adjacent sites on the lattice.

Non-linear switching between two states
The aim of this section is to demonstrate the general theory for a lattice simplified into two states which can be thought of as neighbouring sites. We believe the essential behaviour is preserved in this illustration and that key results can be generalised to the full system.
Let us consider two states denoted by i ∈ {1, 2} with escape rates T i (τ, ρ i (t)). We assume there are two different, and independent, processes which affect the escape rate: the waiting time since last arrival in the state, and the mean number of walkers in the state. We shall discuss two types of non-linearities resulting from considering the mean population in a state: a positive response (leading to reduced escape rates from the state -joint presence of walkers is desirable) or a negative response (increasing the escape rate from the state in order to avoid overcrowding). These both introduce dependencies on the mean-field population of walkers in different ways, as shall be illustrated in the respective escape rates.
In a two-state model all walkers which leave one state must enter the other, where this process is called an "event". Using the method of characteristics (see Appendix A), we can now solve (4); in doing so we assume that the age spent in a state is less than the total time elapsed. This is consistent with the initial condition given by (5). The general solution to (4) is hence given by where τ < t and we can identify a non-linear survival function Ψ i [t, ω|ρ i ] = e − t ω Ti(u−ω,ρi(u))du which is the probability of remaining at site i for a duration of time t − ω in the interval [ω, t). Note that this probability depends on the mean number of walkers in the state. Hence, (9) is the statement that the current mean structural density in a state consists of the earlier arrivals which remained (aged) until now. Using the relation ψ i [t, By substitution of (9) into (6) and using (5), we can find the mean number of walkers in a state which cannot be easily treated without specifying the escape rate. However, in the case when the rate does not depend on ρ i , this simplifies to the convolution of n i (t, 0) and Ψ i (t, ω), and there is thus a correspondence between this equation and the first part of (1). However, (11) considers the mean number of walkers instead of the PDF of position in space of a single walker. It is also useful to introduce a switching term I i (t) for each state defined as (12) where in the second line we have used (10). The switching term gives a flux at which the walkers leave the state, taking into account their ages. Since the walkers' ages are constrained by (5) we need only consider waiting times τ < t. From (8) we further know that in a two-state system I 1 [t|ρ 1 ] = n 2 (t, 0) and vice versa. We can hence again observe that in the case where the escape rates do not depend on ρ i (t), (12) can be simplified in analogy to the second part of (1). We are interested in the total mean number of walkers in each state, and how this quantity is related to the switching terms. This can be obtained by differentiating (6) and using (4) such that where in the last line we have used (12). For state 1 we can then use (8) to write dρ1 dt = I 2 − I 1 and dρ2 dt = I 1 − I 2 . These equations often constitute the starting point of many random walks, with the explicit forms of I 1,2 [t|ρ 1,2 ] to be determined from the details of the random walk [17,[32][33][34]. However, the above method has the benefit of providing a systematic method by which to obtain these expressions. Here we have assumed a constant number of walkers N = i ρ i (t), such that we can write We remind the reader that the trivial case of constant escape rates T i = λ i yields I i [t|ρ i ] = λ i ρ i (t), such that we obtain the well-known steady state distributions of the walkers ρ st 1 = λ 2 /(λ 1 + λ 2 ) and ρ st 2 = λ 1 /(λ 1 + λ 2 ) (where st indicates the value at steady state).
We have demonstrated the equivalence between the two methods in the case when there is no dependence on the mean number of walkers. The generalised expressions, written via Ψ i [t, ω|ρ i ], do not always allow for further analysis of the behaviour of the system as they cannot generally be expressed as the convolution of two quantities. However, in the following case we demonstrate a case where this is indeed possible by the separation of the ageing and non-linear effects. In particular, we consider an escape rate with a negative response to a local increase in walkers.

Non-linear tempering
Let us first consider the case when the escape rate from each state is given by where α i > 0 is a positive function of the mean population of walkers in the state. If α i is an increasing function of ρ i (t) such that ∂αi ∂ρi > 0, then the escape rate (in addition to the anomalous effects resulting from the age dependence) increases with the number of walkers present in the state. An example of such a form is given by α i = βρ k i , where β, k > 0. However, there are no restrictions on the functional form of α i (ρ i (t)). Naturally, for α i = 0 escape events are anomalous with no tempering, and in the constant case ∂αi ∂ρi = 0 there is a constant tempering where the precise number of walkers in the state has no effect. This may well be a suitable approximation if α i is almost unchanged when the walker population varies (very weak dependence on the number of walkers).
The more walkers in the state, the lower the waiting times in said state, which can lead to a tempering of the anomalous effects in the case when µ i < 1. Physical processes which are subject to such conditions often have saturation limits or carrying capacities above which an increase in the number of walkers is unfavourable and there is a trend to decrease the population in the state. Examples include ion concentrations present in certain transporter channels [35] or regulators in the enzymes which facilitate protein folding [36]. Another interpretation of such rates is that rather than functioning as an internal regular (e.g. a carrying capacity), α i simply contains geometrical or physical effects due to the finite size, number of binding sites, binding radius, etc. of the state. When filled to this limit, any new additions are very weakly bound to the state and so the mean escape rate from the state increases.
For the escape rate given in (14), the survival probability takes the form This expression clearly has two factors: the standard anomalous survival function and a non-linear term which depends explicitly on the prehistory via the integral over α i (ρ i ), This second term includes the interactions of the walkers with the mean population in each state and changes the qualitative behaviour of the system. The reason for this is clear: α i introduces decaying exponential terms which eliminate the heavy tails of Ψ i (τ ) at sufficiently long times. We thus have a tempering of the anomalous behaviour of the system. It is convenient to introduce a tempering functional f i [t|ρ i ] to be Hence, in addition to the non-Markovian effects arising from the age-dependence of Ψ i (τ ), there are additional non-linear effects resulting from the entire prehistory of interactions. The form expressed in (15) is not conducive to further analysis besides the observation that the anomalous waiting times have a linear PDF ψ i (τ ) = µi τ0+τ Ψ i (τ ). In order to utilise the convolution properties of Laplace transforms, it is consequently helpful to rewrite (15) as We observe that net survival requires not leaving due to ageing in a certain time interval, and not leaving due to non-linear escape rate contributions in the interval [t − τ, t) resulting from interactions with ρ i . If we define a memory kernel in Laplace space which is valid for any PDF, we thus obtain (see Appendix B) a non-linear switching term. This obeys where it is important that the additional escape rate α i not only leads to a separate flux, but also substantially affects the flux resulting from the escape rate µi τ0+τ via the tempering functional of (16). So even if these escape rates are assumed to be independent of each other, they still couple in the total switching between the two states. The switching term depends on the ages and populations of walkers at all previous times due to its non-Markovian structure, which makes it a highly non-trivial expression for the movement of the walkers between the two states. Similar expressions for the switching terms have been obtained before in the context of proliferating glioma cells [37], reactions in spiny dendrites [38], and persistent random walks with death [39]. Exponential tempering factors have previously been found for subdiffusive reaction-transport processes [40,41]. For state 1 we can then use (13) and (19) to write If the process is purely anomalous with T i (τ ) = µi τ0+τ where 0 < µ i < 1, then in the long time limit s → 0 we obtain the result that Consequently, the memory kernel is given by . This expression leads to the fractional derivative imposed by the Riemann-Liouville operator which in Laplace space obeys L t { 0 D 1−µi t [ρ i (t)]}(s) = s 1−µi ρ i (s) as s → 0 [19]. If both states are anomalous with with different values of µ 1,2 , then the walkers will slowly accumulate in the state with the smallest value of µ i . If µ 2 < µ 1 , at long times one thus observes ρ 2 (t) → N, ρ 1 (t) → 0. If one anomalous process is tempered, effectively increasing the anomalous exponent in that state, this could drastically change the long-term behaviour of the entire system. Unsurprisingly, in a system with one anomalous state and one with a constant escape rate, the walkers will aggregate in the anomalous state. This was shown by Shushin [42] and was recently studied in the context of human residences [43], as well as how aggregation is affected by degradation effects [44].
We shall now consider the effects of introducing an anomalous waiting time PDF into a tempered two-state system. As shown in Appendix B, we obtain the following expression for the switching term It is clear that this expression contains a non-linear tempered fractional operator acting on ρ i (t). The linear case of f i (t) = γ i t has previously been studied in [45,46]. Implementing these switching terms we find that the mean rate of change of the population of walkers obeys and dρ2 dt + dρ1 dt = 0. We thus have two non-linear expressions for the mean number of walkers entering and leaving each state. The fractional derivative is tempered by the non-linear factor e −fi [t|ρi] , which prevents the slow but sure aggregation one would expect in the state with the smallest anomalous exponent.
If one assumes that in the long-time limit a steady state is reached where dρ i /dt = 0, then the mean number of walkers are constants ρ st i ( st denotes values at steady state). Then, where we have used the result that as t → ∞ e −f st [44]. This expression simplifies to where we have introduced the effective escape rates which are only valid as t → ∞ for a system in steady state. Note that these rates λ * i depend on the stationary number of walkers in each state. Typically, τ 0 1, such that λ . If α st i < 1 the rate of movement between the two states λ * i increases. Otherwise, if α st i > 1 the rate λ * i decreases. Equation (27) thus describes the effective rates of movement of walkers between the two states. By manipulation of (26) we find that where is is important to note that these are non-linear equations in ρ st i due to the implicit dependences on α i (ρ i ). The effective rates of movement from one state to the other therefore depend on the chosen form of tempering, but do not necessarily lead to the long-time aggregation of all walkers in one state. However, over shorter time scales the anomalous aggregation effects may still dominate the dynamics.
Furthermore, from (9) we find that in the steady state the mean structural density n st i obeys where from (8) n st . Using the result from (28) it follows that as t → ∞. Note that this expression is also non-linear in the expressions of the effective rates λ * i . If we consider only one of the two states to be tempered, such that state i = 2 has a a constant escape rate k 2 , then (24) becomes We can rewrite (4) for state i = 1 to be In analogy to before, we consider the effects at steady state where the walkers are still ageing but do not change with time. Letting ∂n st 1 /∂t = 0, (32) thus becomes where we have used the result that at steady state α st 1 = α 1 (ρ st 1 ) > 0 is constant as the mean number of walkers no longer changes. The solution to this equation is given by which corresponds to the new arrivals k 2 ρ st 2 and a survival function Ψ 1 (τ ) = τ µ1 0 e −α st 1 τ /(τ 0 + τ ) µ1 . The probability of longer waiting times is still high, but no longer diverges as a result of the tempering induced by the exponential cut-off e −α st 1 τ . The total mean number of walkers in the state can be found by integration: (35) where Γ (a, x) = ∞ x t a−1 e −t dt is the incomplete Gamma function [19]. By definition, we know that the mean waiting time in a state is given by T i = ∞ 0 Ψ i (τ )dτ , which is exactly the above expression. Hence, we can write where . Similarly, T 2 = 1/k 2 . We underscore that this result is only valid when α 1 (ρ st 1 ) > 0, as the integral otherwise diverges for µ 1 < 1. Since the total number of walkers is preserved, we find that We have thus illustrated the various ways in which nonlinearities can be introduced such that anomalous effects are tempered, as illustrated by our key results (28) and (30). That is, the introduction of nonlinearities in the dynamics has suppressed the anomalous behaviour one might otherwise expect, allowing for the existence of a steady-state solution. A crucial point of this as shown in (27) is that the resulting rates change the form of this steady-state solution by their dependence on the steadystate population. However, anomalous behaviour is not completely absent as the rates still depend on the anomalous exponents µ i . We turn now to the case when these non-linearities might induce anomalous effects instead. This is considered for the case wherein an increase in the number of walkers in a state decreases the escape rate.

Self-organised anomaly (SOA)
The previously examined case lends itself most naturally to a tempering of the anomalous behaviour such that the escape rate from the state does not decrease indefinitely. We are now interested in the opposite effect where anomalous behaviour arises in the case when initially µ i > 1. If we consider a modified anomalous exponent which depends on the mean number of walkers of the form µ i (ρ i (t)) and let this be a decreasing function of ρ i (t), we can obtain a trapping state without a priori defining it to be so. In particular, let the escape rate from state i = 1 be and the escape rate from state i = 2 a constant Our reason for this choice is evident: we know that in the long-time limit heavy-tailed waiting times with µ 1 > 1 are analogous to a constant escape rate, and thus for ease of illustration of the phenomenon, one state is fixed at a constant escape rate. Using (13), we find that where again dρ2 dt + dρ1 dt = 0 by conservation of the total number of walkers. In the long-time limit we expect to reach the steady state of the system wherein the movement between states is entirely balanced. In this case, there are no variations with time (∂n st 1 /∂t = 0), although the walkers in each state still age. Let us write the stationary distribution n st 1 (τ ) which from (4) follows Note that since in the steady state ∂n st 1 /∂t = 0, it follows that ρ st 1 is constant such that µ 1 (ρ st 1 ) is also constant. By substitution of (38) and integration we obtain the formula (where we have used the result that n st which depends on the stationary total number of walkers in state i = 1. We then obtain the ρ st 1 -dependent power law survival function such that n st 1 (τ ) = k 2 ρ st 2 Ψ 1 (τ, ρ st 1 ). This is a fairly intuitive result: the stationary distribution of walkers in the state corresponds to the number of walkers which entered the state (k 2 ρ st 2 ), modulated by the probability of surviving in the state for a variety of ages (Ψ 1 (τ, ρ st 1 )). In order to find the total number of walkers in the state we must integrate over all ages using (6) such that It is imperative to note that in the above expression, the result changes depending on the value of µ 1 (ρ st 1 ) when the above integral diverges as the first moment or mean waiting time T 1 → ∞. However, when µ 1 (ρ st 1 ) > 1 we can solve (40) to obtain where at steady state the number of walkers in state i = 2 (ρ st 2 ) can be determined from the constant total number of walkers N in the states. However, this result ceases to hold as soon as µ 1 (ρ 1 ) < 1.
A key feature of the behaviour observed in the states is one of self-organised aggregation. State i = 1 need not start in an anomalous configuration where µ 1 (ρ 1 (t)) < 1; this can occur slowly as the state gradually gains a larger population and eventually serves as a gathering point for walkers. As such, the aggregation is a natural phenomenon and not one we need to introduce artificially. Once there is a sufficient amount of walkers in the state the value of µ will decrease and the steady-state solution calculated in (45) breaks down. That is, the system self-organises between two states: one in which there is a steady state distribution and another anomalous state where steady state is never reached and aggregation ρ 1 (t) → N, ρ 2 (t) → 0 occurs instead. In conjunction with (42), (45) constitutes a key result of this paper where we demonstrate the emergence of anomalous effects from the inclusion of mean-field interactions. Some work has previously been carried out on the self-organising phenomenon in anomalous transport [21,22]. The marked difference in the behaviour of the system depending on the value of µ 1 (ρ 1 (t)) can be regarded as a self-organising phase transition, an already well-documented phenomenon [47][48][49].

Non-linear distributed case with non-local interactions
The aim of this section is to extend the previously considered two-state system to the distributed case. Secondly, we incorporate non-local interactions between the walkers across the entire space of the system. We consider a one-dimensional space R 1 (which can be easily generalised to higher dimensions) where sites are separated by a physical distance a. That is, a site at x will have neighbours at x ± a. As this spacing decreases, we approximate the dynamics of an analogous continuous system. Escapes from a site at x occur with rates T with equal probability to move in either direction. The walker is thus equally likely to jump left or right. However, the interactions between walkers need not be constrained only to those present at the current lattice site but may be anywhere on the lattice. This motion is encompassed by the equation n(x, t, τ ) is analogous to the quantity described in (4), where each site is designated a spatial position instead of an index i. We describe the non-local interactions between the walkers viaρ[x|ρ] defined as where G(x − y) is the appropriate kernel which governs the interactions between walkers at different positions on the lattice. Naturally, if G(x − y) = δ(x − y) we recover an exclusive dependence on the local number of walkers at the current site. The mean total population of walkers at a certain position in analogy to (11) follows but the survival function must now vary with space andρ instead of ρ. The switching term is set to be which is analogous to (12) but with different spatial dependence. By the same method employed in (13), we find that Recent arrivals n(x, t, 0) at a certain position must come from the adjacent neighbours, such that n(x, t, 0) = Notice that the term in the bracket is the second spatial derivative in discrete space, such that we can write If one assumes escape rates which are qualitatively similar to those present in anomalous tempering (see (14)), Here, the interaction of the walker with the mean neighbourhoodρ[x|ρ] thus affects the escape rate from a site, where any additional features of the interactions can be included into the kernel G(x − y) defined in (47). Previous works have considered the consequences of spatially varying anomalous exponents µ(x) (see [33,[50][51][52]). By the same methods as before, the survival function is given by In analogy to before we can identify a spatially dependent , and a non-linear term which depends explicitly on the prehistory via the integral overρ. As before, it is sensible to introduce a tempering functional which we notice is spatio-temporal by the dependence on both prehistory and walkers in other lattice sites. This spatio-temporal fractional tempering appears in the total survival probability which equals where Ψ (x, τ ) = . It follows that the mean total population of walkers is given by Using (49) and the same methodology employed for (B.2), we find that where we have used ψ(x, τ ) = µ(x) τ0+τ Ψ (x, τ ) for an anomalous waiting time PDF ψ. By rearranging and taking the Laplace transform analogously to the result of (19), where the memory kernel K(x, τ ) is defined in Laplace space via the expression By substitution of (59) into (52), we obtain the key result which is only valid in the case when 0 < µ(x) < 1. In the case whenρ(x, t) = 0, we find the simpler form with a spatially dependent diffusion coefficient )] −1 similar to the result derived in [33]. Note that this is equivalent to (3) when µ does not vary with x.

Discussion and conclusion
In this work we have introduced a general methodology for the inclusion of mean-field interactions where the statistical properties of a random walker depend on the average walker population ρ, either locally or non-locally. The result is a significant qualitative change in the behaviour of the systems as compared to standard results without this mean-field dependence. In the past tempering in such systems has appeared as an imposed exponential factor with the fractional derivative. In our work we do not modify the residence time distribution, but rather obtain tempering of long residence times as an emergent phenomenon resulting from non-linear interactions. We thus manage to obtain a natural tempering from a mean-field dependent escape rate which does not suppose linearity of walker interactions. We have illustrated how escape rates which depend on the mean number of walkers can increase or decrease the waiting time PDFs. This seemingly innocuous effect has the ability to suppress the anomalous behaviour of the system by decreasing the waiting times of the random walker via well-established population effects such as volume exclusion, carrying capacities, etc. (see (28) and (30)). As a result of this tempering, the dynamics are instead governed by non-linear effective rates like (27) which depend on the average walker population. The implicit dependence of (27) on the steady-state population changes this steady-state population from what it otherwise may have been. Crucially, it is the presence of the tempering which allows this steady state to be obtained, and its anomalous behaviour is not completely absent as the rates still depend on the anomalous exponents µ i .
We have also shown that the inclusion of mean-field effects generates anomalous behaviour by an increase in the waiting times in a certain state. This is a direct result of these non-linear population effects and does not arise in the linear case of no population dependence. This generation of anomalous behaviour corresponds to selforganising behaviour resulting from the interaction of a random walker with the mean population, and may arise even for initial values of µ > 1. An increased mean number of walkers in a certain state thus has be ability to further the attractiveness of the state, effectively decreasing µ(ρ i ) to fall below one (see (42) and (45)).
We have also shown that non-local interactions can be incorporated into the distributed case of a series of walkers arranged on a lattice in R 1 . We have derived a new fractional diffusion equation (61) subject to spatio-temporal tempering which lends itself to generalisation in higher dimensions.
Future work will explore the repercussions of these nonlinearities on the resulting fractional diffusion equations. S. Fedotov gratefully acknowledges the support of EPSRC grant EP/N018060/1.

Author contribution statement
The paper has been jointly developed, where both authors have contributed to the derivation of key results and development of ideas. The article has been drafted by H.

Stage.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Appendix A: Method of characteristics
This appendix demonstrates how to calculate the survival function of random walkers subject to non-linear mean-field interactions using the method of characteristics. Starting from (4) and considering movement along a characteristic curve s, we find dn i ds = −T i (τ (s), ρ i (t(s))) n i , (A.1) such that we can find dt/ds = 1 and dτ /ds = 1. Since we assume τ < t, it follows that t − t 0 = τ where t 0 > 0 is the time at which the current residence in the state started. The characteristic is chosen to coincide with the age such that τ | s=0 = 0 and t| s=0 = t 0 = t − τ . Integrating the above equation we find Applying the above values and a change of variable u = t(s), it follows that: where using t 0 = t − τ yields (9).

Appendix B: Derivation of a non-linear switching term
This appendix demonstrates how to calculate the switching term from a given survival function (17) of the random walk. By substitution of (17) into (11), we find where we have used (6), (12) and (17). One can rewrite this such that Note that this is expressed as a modified Riemann-Liouville operator, such that we obtain (23).