Particle-Environment Interactions In Arbitrary Dimensions: A Unifying Analytic Framework To Model Diffusion With Inert Spatial Heterogeneities

Abridged abstract: Inert interactions between randomly moving entities and spatial disorder play a crucial role in quantifying the diffusive properties of a system. These interactions affect only the movement of the entities, and examples range from molecules advancing along dendritic spines to anti-predator displacements of animals due to sparse vegetation. Despite the prevalence of such systems, a general framework to model the movement explicitly in the presence of spatial heterogeneities is missing. Here, we tackle this challenge and develop an analytic theory to model inert particle-environment interactions in domains of arbitrary shape and dimensions. We use a discrete space formulation which allows us to model the interactions between an agent and the environment as perturbed dynamics between lattice sites. Interactions from spatial disorder, such as impenetrable and permeable obstacles or regions of increased or decreased diffusivity, as well as many others, can be modelled using our framework. We provide exact expressions for the generating function of the occupation probability of the diffusing particle and related transport quantities such as first-passage, return and exit probabilities and their respective means. We uncover a surprising property, the disorder indifference phenomenon of the mean first-passage time in the presence of a permeable barrier in quasi-1D systems. We demonstrate the widespread applicability of our formalism by considering three examples that span across scales and disciplines. (1) We explore an enhancement strategy of transdermal drug delivery. (2) We associate the disorder with a decision-making process of an animal to study thigmotaxis. (3) We illustrate the use of spatial heterogeneities to model inert interactions between particles by modelling the search for a promoter region on the DNA by transcription factors during gene transcription.


I. INTRODUCTION
Local interactions between mobile agents or particles and their environmental features plays a crucial role in the dynamics of many systems across disciplines and scales [1][2][3][4][5].When such environmental features are inert heterogeneities, the local interactions only affect the movement dynamics of the agents.A wide array of spatial heterogeneities can be classed as inert, e.g.impenetrable or permeable barriers, areas of reduced or increased mobility, lattice defects such a disclinations, and traps that are reversible.
In some instances the presence of such heterogeneities is by design, e.g. in manufacture engineering where ma-terials are constructed to have specified diffusive characteristics [6,7].In other scenarios spatial heterogeneities occur naturally.In ecology, animals alter their foraging behaviour due to variations in vegetation cover [8,9].In molecular biology, particles undergo fence hindered motion in the lipid bilayer membranes of eukaryotes [10,11], and slow down dramatically when moving within the cell cytoplasm due to exclusion processes [12].While the relationship between mobility and spatial disorder in these and other systems has always been a focus of scientific studies, it is the highly resolved nature of modern observations that has made apparent the need for a general framework to model inert particle-environment interactions.
Investigations on movement dynamics in spatially disordered systems date back as early as the 50's [13][14][15][16][17].Despite such a long history most analyses lack a rigorous quantitative description of the 'microscopy' of the arXiv:2209.09014v1[cond-mat.stat-mech]14 Sep 2022 interaction events between the particle and the environment.In the past, transport in highly disordered media has been studied approximately, linking the Hausdorff-Besicovitch dimension of fractal structures to a diffusion constant via scaling arguments [18].Other approaches have kept the geometry non-fractal utilising random walks on regular lattices, the so-called random walk in random environments model [19][20][21][22][23][24][25].The majority of these latter studies pertain to 1D domains [26], and have used techniques such as the effective medium approximation to find statistical properties of the movement dynamics.It is precisely the absence of explicit spatiotemporal representation of higher dimensional particle dynamics, that has hampered the widespread applicability of these various models to current high fidelity observations.
More recent theoretical applications to movement in disordered environment have focused on the diffusive dynamics in the cell (e.g.see reviews in Refs.[27][28][29]).Many attempts in this area are macroscopic and tackle particle dynamics without representing the local interactions.Some efforts, giving importance to the very slow dynamics that emerge from overcrowding effects, have modelled particle movement via fractional diffusion [30].The relative size of accessible versus inaccessible regions has been accounted for using diffusion on percolation clusters and has highlighted the difference between compact versus non-compact exploration of space [31].Other investigations have put emphasis on the spatiotemporal dynamics of the environment and have developed the so-called diffusing-diffusivity models, where the diffusion strength of the medium itself is a random variable [32][33][34][35].
These approaches have brought important insights and have broadened the tools and techniques with which to study disordered systems.However, they too lack the mechanistic connection between the environmental heterogeneities and the moving particle [30].With the advent of new experimental techniques such as super-resolution microscopy and single particle tracking [36], the need for an explicit consideration of particleenvironment interactions has also emerged in microbiology [37,38].
The challenge in fulfilling this need stems from the symmetry breaking role that disorder plays on the underlying diffusive dynamics.In most instances describing explicitly multiple heterogeneities is an unwieldy boundary value problem.The vast majority of theoretical studies have in fact been limited to highly symmetric scenarios, e.g.spherically symmetric domains with concentric layers of different diffusivity [39][40][41][42] and an array of periodically placed semi-permeable barriers in 1D [43][44][45][46][47].
To bypass this challenge, and to avoid the use of computationally prohibitive stochastic simulations, we propose a unifying analytic framework to model interactions between diffusing agents and spatial disorder.We do so by developing a random walk theory where interactions with heterogeneities are represented as a perturbation of the transition probabilities of a homogeneous lattice.By extending the so-called defect technique [48][49][50], we are able to model explicitly any inert particle-environment interactions in arbitrary dimensions, e.g. the passage through porous or permeable barriers, the movement within regions of altered diffusivity, which we call sticky or slippery sites as well as shortcut jumps to far away locations.
The theory allows us to derive mathematical expressions for the random walker occupation probability, the so-called propagator.The generating function of these propagators are exact and obtained in terms of the occupation probability in the absence of spatial heterogeneities, thereby making our framework modular in its application.Multiple derived quantities, such as first-passage, return and exit probabilities, which in the past were obtained either numerically or known only in asymptotic limits [51], can now be readily computed via the evaluation of certain matrix determinants.
Given the generality of our framework, we have opted to provide three examples of application.The first deals with an extra-cellular process, namely the potential optimisation of transdermal drug delivery [52,53].The second example is the modelling of thigmotaxis, the tendency of insects and other animals to remain preferentially close to physical boundaries whilst moving [54,55].The third application concerns with the search dynamics in a two particle coalescing process that is of relevance to early stages of gene transcription [56,57].
The remainder of the paper is organised as follows.In Section II we introduce the general mathematical formalism via a lattice random walk Master equation, and show how we represent different kinds of heterogeneities.In Section III we solve the Master equation and find the exact propagator.Section IV deals with first-passage statistics and their associated mean, i.e. mean first-passage, mean exit and mean return times.The latter half of the paper, Sections V, VI and VII, are devoted to the three applications mentioned previously, which are transdermal drug delivery, thigmotaxis and gene transcription.Lastly, conclusions and future applications form Section VIII.

II. MOVEMENT IN HETEROGENEOUS ENVIRONMENTS
We start by defining the dynamics of a Markov lattice random walk on a d-dimensional lattice via where n is a d-dimensional vector and A n,m is the transition probability from site m to site n such that m A m,n = 1 for any site n on the lattice, i.e. with d = 1, A is a probability conserving transition matrix, and when d > 1, A is actually a tensor.For convenience in inverting generating functions, as compared to Laplace inversion, we use a discrete time formulation with the variable t.Changes to a continuous time description is straightforward [58], but is omitted here.We refer to this equation as the homogeneous Master equation and its solution, given a localised initial condition, as the homogeneous propagator.The underlying lattice is referred to as the homogeneous lattice whose size can be finite or infinite.
Since spatially heterogeneous dynamics are defined relative to the homogeneous system, we define heterogeneities as locations or defects where the dynamics are different from the corresponding ones on the homogeneous lattice.Examples of heterogeneities are depicted in Fig. 1.The heterogeneities displayed in Fig. 1 emerge from the modification of the outgoing transitions from one or more sites, hence we refer to these altered transitions as heterogeneous connections.For example, given a partially reflecting barrier in between two neighbouring sites, the jump probability from either of the two sites to the other is reduced, while the probability of staying put at either of the sites is increased.Conversely, by connecting together two non-neighbouring sites, we may wish to reduce the probability of staying put at a given site, whilst adding the possibility of hopping to the site further away.One can represent conveniently these or any other heterogeneity through a modification of the transitions as depicted in Fig. 2. Formally, the outgoing connections of the sites u and v are adjusted by introducing the parameters λ v,u and λ u,v to create two heterogeneous connections.Although we choose to modify transitions in both direction, i.e. from u to v and v to u, this does not have to be the case.Modifications of only outgoing connections are also permitted, e.g.see the dashed arrows connecting the additional sites r and s.
FIG. 2. A schematic representation of the transition probabilities after the introduction of spatial heterogeneity or disorder.The probability of hopping from site u to v is given by A v,u .When λv,u is positive, the probability of jumping from u to v decreases, while the probability of staying put increases.When λv,u is negative, the opposite effect occurs with a decrease in the probability of staying, while increasing the jump probability from u to v. The parameter λu,v affects the transition probability from v to u and the probability of remaining at v in an equivalent manner.
The construction implicitly conserves probability, which can be evinced by picking a defect, e.g.u, and summing over all of the outgoing probabilities.The changes induced by the λ parameters cancel each other out leaving w A w,u , with w representing all the neighbours of u, equal to the homogeneous outgoing probability.To ensure positive probability for a given heterogeneous site u, we have the conditions for all w with a heterogeneous connection in the direction u to w, and which enforces upper and lower bounds on the λ parameters although each one of them can be positive or negative.This formulation allows one to perturb arbitrarily the homogeneous lattice creating any type of probability conserving particle-environment interactions.

A. Quantitative representation of heterogeneities
To understand the practicality of the formalism, we focus on the three specific types of heterogeneities in Fig. 1, namely, barriers (Figs.(1)a and (1)b), long-range connections (Fig. 1c) and sticky sites (Fig. 1d).In the following subsection we present convenient parameterisation for the constant λ's to construct such heterogeneities.

Barriers and Anti-Barriers
With u and v two neighbouring sites, we construct a partially reflecting barrier by having is a measure of the reflectivity of the barrier.When α v , α u = 1 we have an impenetrable barrier (shown in Fig. 3), while with α v , α u = 0 we regain the homogeneous transition.Notice that the barrier does not need to be symmetric, i.e. α u = α v , with the extreme scenario being a barrier with λ v,u = A v,u and λ u,v = 0 yields a one-way barrier or gate.In such a case the movement from v to u is allowed but from u to v is not.
It is also possible to have dynamics opposite to the partially reflecting barrier.In this case, again with u and v two neighbouring sites, one has λ As the probability of jumping to the neighbours increases whilst the probability of staying put decreases, we have chosen the name anti-barrier for this type of heterogeneity.
Example of a reflecting barrier between u and v generated by modifying the transition probabilities (from the left to the right of the schematic).The modified transitions are indicated by coloured arrows.The modification in this case results in an impenetrable barrier between u and v, with αv = αu = 1.

Long Range Connection
When adding an outgoing long-range connection one has to draw the probability from one or more of the existing transitions.Let us consider a site u and a nonneighbouring destination site s, where A s,u = A u,s = 0.One way of introducing the outgoing long range connection is to draw upon the lazy (also called sojourn) probability using λ s,u = −β s A u,u and λ u,s = −β u A s,s , where β u , β s ∈ [0, 1] is the proportion of the lazy probability added to the long range connection, see Fig. 4 for a pictorial representation.
Note this is not the only way; one can also rewire an existing connection from a neighbour to the non-neighbour.In such a case, with v a neighbour of u, we let λ v,u = A v,u and λ s,u = −A v,u .The former removes the possibility of jumping from u to the neighbour v, whilst the latter adds the possibility of hopping from u to the non-neighbour s.
Example of a long range connection obtained by rewiring the lazy probability of u and v to create a long range connection between them (from left to right).In this case βu = βs = 1.

Sticky and Slippery Sites
Adding a partially reflecting barrier between two neighbouring sites naturally increases the probability of staying.By harnessing this property one can use multiple one-way partially reflecting barriers between a site w and all of its k nearest neighbours, r 1 , • • • , r k , giving λ ri ,w = αA ri,w , and λ w,ri = 0 with α ∈ [0, 1] for all i = 1, • • • , k.The result is a sticky site, w, where the probability of staying is increased, whilst the probability of jumping to any of its neighbours is decreased.The introduction of α is used to control and distribute the stickiness equally across the neighbours in a convenient manner.See Fig. 5 for a pictorial representation on a 1D lattice.
Conversely, keeping λ w,ri = 0 and letting λ ri ,w = − β k A w,w with β ∈ [0, 1] for all i = 1, • • • , k yields a slippery site with opposite dynamics.As for the sticky site, the introduction of β is used to control the slippery quality of the site w equally among its neighbours.Note that we have chosen to divide β by k so that Eq. ( 3) is automatically satisfied.

III. HETEROGENEOUS PROPAGATOR
We consider an arbitrary collection of heterogeneous connections given by a set of M paired defective sites or defects, S = {{u 1 , v 1 }, • • • , {u M , v M }}.We use u i and v i with subscripts to indicate the two members of the i th pair, while u and v without subscripts refers to a generic pair in S. The pairs are unique, i.e. {u i , v i } = {u j , v j } for any i = j, however, a site can be part of multiple pairs.For example, the set of pairs which represents the schematic depicted in Fig. 2 is S = {{u, v}, {u, r}, {v, s}}, with the sites u and v being part of two pairs while the sites r and s being part of only one pair each.The evolution of the occupation probability is given by the Master equation (for the details see Section I of the Supplementary Materials) where the second summation is over all pairs of heterogeneous connections.When all λ parameters are set equal to zero, Eq. ( 4) reduces to Eq. ( 1) and the occupation probability on the heterogeneous lattice, Φ(n, t), reduces to that of the homogeneous lattice, ϕ(n, t).
One can find the generating function (z-domain) solution of Eq. ( 4) by generalising the so-called defect technique to obtain (see Section I of the Supplementary Materials) where In Eqs. ( 6) and (7) we have used the notation From here onwards we refer to ϕ n0 (n, z) as the homogeneous propagator, which is known explicitly (see Sections V and V B of the Supplementary Materials), while Φ n0 (n, z) is referred to as the heterogeneous propagator.When t = 0, that is z = 0, we have ϕ n0 (n, 0) = δ n0,n , while |H(n, n 0 )| / |H| = 1 and we recover the appropriate initial condition, Φ n0 (n, 0) = δ n0,n .
In general the size of matricies H and H(n, n 0 ) depends on the number of paired defects, M .A ddimensional walk with one sticky (or slippery) site requires two paired defects for each of the d dimensions.However, in this case one can make a simplification and reduce the size of the matrices by a factor of 2d.The simplified matrices H and H(n, n 0 ) are defined, respectively, in Eqs.(S20) and (S21) of the Supplementary Materials.FIG. 6.A snapshot of Φn 0 (n, t) at time t = 100 obtained from Eq. ( 5) with standard numerical methods [59,60].Propagator with different configurations of defects, corresponding with Fig. 1, at t = 100.For all panels, the homogeneous propagator, ϕ n 0 (n, z) is the 2D propagator with reflecting boundaries given in Eq. ( 23) of Ref. [58] (see also Section V B of the Supplementary Materials).The parameters used are: a domain of size N = (10, 10), a localised initial condition with n0 = (6,6).We a use diffusion parameter of value q = (0.2, 0.2), which gives the following transition probabilities: in the bulk of the homogeneous lattice the probability of jumping to one of the four neighbours is A r,s = 0.05 (with r = s), while the probability of staying at the same site is A r,r = 0.8.The reflecting barriers and other heterogeneities are super imposed on top of the probability.For panels (a) and (b), λv,u = A v,u and λu,v = A v,u (for all {u, v} ∈ S) yielding perfectly reflecting barriers.For panel (c) λv,u = − 1 2 A u,u , and With this perturbation, when on one of the defective sites, the probability of staying is reduced to A u,u = A v,v = 0.4, while the probability of jumping to the non-neighbour is increased (from zero) to A u,v = A v,u = 0.4.Lastly, for panel (d), for each of the sticky-sites w with k neighbours r1, • • • , r k we use λr i ,w = 1 4 A r i ,w (see Section II A 3).For convenience we have omitted colour bars for each panel as we are interested only in the relative differences of the occupation probability.
In Fig. 6 we plot a snapshot of Φ n0 (n, t) for the heterogeneities depicted in each of the panels of Fig. 1.In panel (a) the lattice is partitioned by impenetrable barriers represented by the solid white lines.Here, one can observe the lowest probabilities in the top-left corner since the walker has not had the time to travel around the barriers.Panel (b) contains areas enclosed by impenetrable barriers, with occupation probabilities that are always zero.The long range connection shown in panel (c) has enabled the walker to spread further than in other pan-els.Small peaks in the probability can be observed away from the initial condition, in the top-left, bottom-left and bottom right corners.In panel (d) the sticky regions tend to show a higher occupation probability compared to the homogeneous sites.
Note that we have not placed any restriction on whether (the homogeneous propagator) ϕ n0 (n, z) conserves probability or note.When there are fully or partially absorbing sites, one may proceed in two ways.(i) In the first approach one account for the absorbing dynamics by finding the propagator ϕ n0 (n, z) that satisfies appropriate boundary conditions, before adding inert disorder via Eq.( 5).(ii) In the second approach one would take ϕ n0 (n, z) without any absorbing locations, construct Φ n0 (n, z) and then add the absorbing sites using the standard defect technique in the presence of absorbing sites [50].While the choice makes no impact on the final dynamics, depending on the situation one procedure may be more convenient than the other.

IV. FIRST-PASSAGE PROCESSES
An important quantity derived from the propagators is the first-passage statistics to a set of targets.It is relevant to stochastic search in movement ecology [61], swarm robotics [62] and many other areas [63].
The first-passage probability, F n0 (n, t), that is the probability to reach n for the first time at t having started at n 0 , is related to the propagator, Φ n0 (n, t) by the renewal equation.When n = n 0 , the well-known relation in z-domain is given by Having the first-passage probability in closed form allows one to substitute the heterogeneous first-passage probability F n0 (n, t) (or F n0 (n, z)) in place of the homogeneous counterpart in other established contexts where homogeneous space was previously assumed.
One such context is a first-passage in the presence of multiple targets, where one is interested in the probability of being absorbed at any of the targets.We use recent findings [58] to determine the dynamics of a lattice walker to reach either of two sites, n 1 , and n 2 , for the first time at t in the presence of spatial heterogeneities, given by F n0 (n 1 , n 2 , t).The generating function of this probability, given by −1 , (taken from Eq. ( 38) in Ref. [58]), is expressed in terms of the first-passage probabilities to single targets.In Fig. 7 plot the time dependent probability for the heterogeneity examples shown in Fig. 1.The first nonzero probability corresponds to the length of the shortest path to either of the targets, which in the absence of heterogeneities and for the examples in Fig. 1(a), (b) and (d) is 6, whereas for the panel (c) one can reach the target n 1 from n 0 in 4 steps as a result of the nearest long range connection.It is clearly visible in the earlier rise of the curve related to Fig. 1(c).Interestingly, the first-passage probability curves corresponding with excluded regions, shown in Fig. 1(b), and the homogeneous case are almost indistinguishable from each other for 2 decades.While excluding parts of the lattice increases the lengths of some paths to the targets, it also reduces the overall space that can be explored.For the setup chosen, these two effects counteract each other at short and intermediate timescales.
Among all the curves, the case with open partitions related to Fig. 1(a), results in a first-passage probability which is the slowest to rise and with the broadest tail in the distribution.The reasons for such characteristics compared to all other curves is due to the location of the initial site relative to the targets.As the latter ones are partially behind partitions, the more directed paths take more time to reach the targets and the walker remains confined in the region around the initial site for much longer.
The sticky sites in Fig. 1(d) have limited effect on the more directed paths connecting the starting site and the targets.This is why F n0 (n 1 , n 2 , t) in Fig. 7(b) is identical to the homogeneous case at short times.However, stickysites can be both a hindrance or a benefit to the searcher.While it can partially trap the walker and stop it from reaching the target site, it can also stop the walker from exploring regions away from the targets.Since there are sticky sites close to the targets, these two effects counteract one another and we observe marginal difference in the tail of the distribution when compared with the homogeneous curve.

A. Explicit mean first-passage quantities
The first moment of F n0 (n, t), that is the mean firstpassage time (MFPT), , is given by where F n0→n is the homogeneous MFPT from n 0 to n and the elements of the matrices H, H (1) and H (2)  are defined in terms of homogeneous MFPTs.They are given, respectively, in Eqs.(A3) to (A5) for general heterogeneities, while for the case with only sticky or slippery heterogeneities the matrices are given by Eqs.(S82) to (S84) in the Supplementary Materials.In the coming sections we use the mathfrak notation e.g.F, R and E, for statistics involving the heterogeneous dynamics, while the mathcal notation, e.g.F, R and E, is reserved for the homogeneous counterpart.The dependence on the target at n is only present in the matrices H (1) and H (2) ; the dependence on the initial condition n 0 is only present in H (1) ; and the dependence on the location of the heterogeneities are in all three matrices.FIG. 7. Time-dependent first-passage probability to either of two targets in the presence of heterogeneities.The location of the targets n1 = (4, 2) and n2 = (10, 7) in relation to the initial condition n0 = (6, 6) and visualisation of the heterogeneities present can be seen in the schematic diagram in Fig. 1.We use a homogeneous propagator, ϕ n 0 (n, z) with a reflecting domain of size N = (10, 10) and a diffusion parameter of value q = (0.2, 0.2).The explicit form of ϕ n 0 (n, z) is given by Eq. ( 23) of Ref. [58].The lines are obtained through numerical inversion of the generating function of the first-passage probability to either of two targets (see text), while the corresponding marks-shown only in panel (a)-are obtained through 1.5 × 10 6 stochastic simulations.
The probability distribution of the first-return time is related to the propagator via R(n, z) = 1 − Φ −1 n (n, z), with a mean return time (MRT) given by where R n is the homogeneous mean return time.
When the heterogeneities preserve the symmetric properties of the homogeneous lattice, i.e. the disorder does not add any bias to a diffusive system or remove any bias present in a system with drift, then the ratio is satisfied, for all {u, v} ∈ S, and the heterogeneous system maintains the steady state of the homogeneous system.In this case, H (2) = 0, the MFPT given by Eq. ( 8), can be simplified to F n0→n = F n0→n − 1 + H − H (1) / |H|, while the MRT remains the same as the homogeneous MRT, R n = R n as expected from Kac's lemma [64].(see Appendix A 1 and Section II C of the Supplementary Materials).
In the presence of multiple targets at the outer boundary of the domain, we relate the first-passage probability to any of the targets to a propagator with the appropriate absorbing boundaries.In this case, the first-passage is referred to as the first-exit, and its probability generating function is related to the propagator through the relation E n0 (z) = 1 − (1 − z) S n0 (z), where S n0 (z) is the survival probability.Taking the mean of the distribution (see Appendix A 2) gives where E n0 is the mean exit time starting at n 0 without any heterogeneities, and the matrix S(n 0 ) is given explicitly in Eq. (A11).The presence of one or more absorbing boundaries on the homogeneous propagator ϕ r (s, z) allows for a simple evaluation at z = 1.That is to say ϕ r (s, z = 1) is finite for any r and s in the domain; and therefore H and S(n 0 ) also remains finite and can be easily evaluated.
In Fig. 8 we show the effect of randomly distributed barriers and anti-barriers as a function of the barrier strength in a 2D domain with absorbing boundaries.The M neighbouring defective site pairs are uniformly distributed on the lattice with λ = λ v,u = λ u,v for all {u, v} ∈ S. One can see that for λ > 0, E n0 increases, as the heterogeneous connections behave as a partially reflecting barrier slowing down the walker.Furthermore an increase in the number of heterogeneities resulting in larger exit times.Conversely, when λ < 0 the heterogeneous connections become anti-barriers increasing the probability of jumping across compared to the homogeneous case, which effectively increases the spread of the walker leading to shorter exit times.When the barriers are impenetrable, increasing the number of barriers also increases the likelihood of the walker being trapped and unable to reach the boundary and will cause the MET to diverge.Although we do not study it here, a similar setup could be used to analyse percolation in finite The explicit form of ϕ n 0 (n, z) is given by the z-transform of Eq. ( 23) of Ref. [58].Each curve is obtained using Eq.(10) and performing an ensemble average with 10 2 sample realisations of locations of barriers (λ > 0) or anti-barriers (λ > 0) for each λ.We consider a simple spatial heterogeneity in a 1D domain with a partially reflecting barrier between u and u + 1.To study the dependence of the position and strength of the barrier (or anti-barrier) on the firstpassage dynamics, we first fix the position of target and initial sites with n > n 0 ; assume a reflecting boundary between n = 0 and n = 1; and take λ u,u+1 = λ u+1,u = λ with λ ∈ [−(1 − q), q/2].In this case the first passage probability can be written using the convenient notation where , and with the probability of moving given by q ∈ (0, 1].The homogeneous first-passage probability, F n0 (n, z), can be recovered from Eq. ( 11) by letting λ → 0. When the barrier is to the left of the initial condition, the limit λ → q 2 creates an impenetrable barrier, the behaviour is equivalent to shifting both the target and the initial condition to the left by u giving F n0 (n, z) = F n0−u (n − u, z).Whereas, when n 0 ≤ u < n, the same limit gives F n0 (n, z) = 0 as the walker becomes blocked by the barrier and can never reach the target.In Fig. 9, we plot the time dependence of Eq. ( 11) for the two different scenarios u < n 0 and u ≥ n 0 represented, respectively, by panels (a) and (b).With u < n 0 and λ < q/2, that is the barrier to the left of the initial condition, as one increases u from u = 1 one observes an increase in the modal peak.When the walker is reflected by the permeable barrier, it stops the walker from straying further left and effectively reduces the space that can be explored, increasing the probability of reaching the target at an earlier time.However, if the walker passes through the barrier, the partial reflection dynamics becomes a hindrance: the walker is kept in the range [1, u], causing the probability of reaching the target at long times to increase also.As probability in the tail and the mode increases, the probability conserving F n0 (n, t) demands a reduction at intermediate times, which is clearly visible from the figure.This permeability induced modetail enhancement can also be witnessed by fixing u and changing λ ∈ [0, q/2), and we have also observed the inverse effect, mode-tail compression by having anti-barrier with λ ∈ [q−1, 0].We have chosen not display these latter cases for want of space.Similar features have been observed in a diffusing diffusivity model in Ref. [35], where increases in the probability at short and long timescales were attributed to the dynamic diffusivity.Our findings point to the fact that such richness can also emerge from a static disorder at a single location.
Differently from the case when the barrier is to the left of the initial condition, is the case when u ≥ n 0 .In this scenario, the barrier is always acting to slow the search process down, reducing the probability of reaching the target at early times and increasing the probability at long times as seen by the flattening of the mode and the broadening of the tail, as shown in Fig. 9(b).
Computing the mean via either Eq. ( 11) or from simplifying Eq. ( 8) yields the compact expression where F n0→n = (n − n 0 )(n + n 0 − 1)/q is the 1D homogeneous MFPT for n 0 ≤ n (given by Eq. ( 14) of Ref. [58]).Astonishingly, the mode-tail enhancement present in the time-dependent probability when u < n 0 has no effect on the mean.This is what we have termed the disorder indifference phenomenon.
To explain why there is such an effect of disorder indifference, we split the first passage trajectories into two mutually exclusive subsets: the trajectories that never return to the initial condition before reaching the target site on the right and those that return at least once before reaching the target site.Clearly, the former trajectories are unaffected by the presence of a barrier.The latter trajectories can be affected by the barrier, however, in computing the mean one deals with mean return times which are unaffected from the homogeneous case when λ u,u+1 = λ u+1,u as stated in the previous section (see Appendix B 1 for the mathematical details).
An analogue of this indifference phenomenon was observed in Ref. [40], where the MFPT in a quasi-1D domain in continuous space with two layers of different diffusivity was studied.When the initial condition was in between the interface of the layers and the target, they observed that the MFPT was indifferent to the diffusivity of the media beyond the interface.In that study, the first-passage probability was not considered and the cause of this indifference could not be quantified.However, one can relate the location of the interface of their system with the position of the barrier in ours.Through this relation, we believe that the behaviour observed in Ref. [40], is closely related to the dynamics presented in Fig. 9.
Given a barrier between the initial condition n 0 and the target n, the effect on the MFPT increases linearly as the displacement from the boundary increases.While the effect, which can be to speed up (λ < 0) or to slow down (λ > 0), is due to the disorder, the linear dependence is not.This linear dependence is present in all 1D situations and is proportional to the distance between the initial condition and the reflecting boundary (see Appendix B 2).
To explore the effects of asymmetry in the heterogeneities we consider the MRT with λ u+1,u = λ u,u+1 .In this cases, the steady state is no longer homogeneous, effectively creating an out-of-equilibrium system due to the loss of detailed balance in the Master equation.To illustrate this point, we consider the MRT of a 1D walker within a segment of length N with reflecting boundaries and with a barrier between u and u + 1, (λ u+1,u = λ u,u+1 ).In this case, Eq. ( 9) simplifies to , n < u + 1, One can see that when λ u+1,u = λ u,u+1 , the MRT reduces to N regardless of whether n ≤ u or n > u.In the extreme case, where the barrier is impenetrable in both direction, λ u+1,u = λ u,u+1 = q/2, one can recover the appropriate MRTs when n ≤ u and n > u, which are, respectively, u and N − u (see Section IV B of the Supplementary Materials).
Thus far we have focused on technical development and theoretical insights.As we move forward, the remainder of the article is devoted to practical examples, and is used to demonstrate the applicability of the framework.For practical convenience the details of the modelling set up are given in the appendices and only the results are discussed.

V. TRANSDERMAL DRUG DELIVERY
In the first application we consider the problem of optimising transdermal drug delivery, that is the transfer of drugs through the skin.One of the challenges of transdermal drug delivery is traversal of the outer-most layer of the epidermis called the stratum corneum (SC) by hydrophilic molecules [65].This layer is made up of dead cells called corneocytes which are arranged in a dense 'brick-and-mortar' like pattern [66].Inspired by some of the recent strategies proposed to enhance drug absorption [67], we consider the use of a micro-needles to pierce first the SC before applying a drug patch.We study the effectiveness of this method by using our modelling framework to represent the SC as heterogeneities on a lattice and modelling the movement of drug molecules as a random walk.10.Representation of the 'brick-and-mortar' arrangement of the corneocytes in the stratum corneum (SC).The red square depicts the starting location of the random walker.The geometry is given by an an absorbing boundary at n1 = 1, a reflecting boundary at n2 = N1, shown as a thick solid black line, and a periodic boundary in the second dimension, depicted as dashed black lines.By using a number of paired defects, one is able to cordon off sites (shaded grey), creating the 'brick-and-mortar' pattern of the SC.The dashed blue rectangle with width w and height h models the destruction of the SC structure via a micro-needle puncture, with h and w representing, respectively, the puncture height and width.
This destruction may open up some of the 'bricks', allowing the walker to easily travel inside.The initial position of the walker is at the centre-top of the puncture, n0 = (N1, N2/2).
We use a homogeneous 2D nearest-neighbour random walker subject to mixed boundary conditions: an absorbing boundary located at n 1 = 1 and a reflecting boundary located at n 1 = N 1 , for the first dimension and a periodic boundary condition on the second dimension.The heterogeneities are impenetrable barriers representing the lipid matrix.These are arranged in a manner to create excluded regions that form the 'brick-and-mortar pattern of the SC, see Fig. 10.The pattern is partially destroyed to represent the needle piercing in a rectangle with height h and width w resulting in an area absent of barriers as shown by the blue dashed rectangle in Fig. 10.
The quantity of interest is the MET with an initial condition starting at the reflecting end of the domain.We plot the MET as a function of the puncture depth and width in Fig. 11.The overarching qualitative changes in the MET can be explained by two competing effects.The first is the breaking of enclosed bricks to create open partitions.The additional sites available for exploration makes the paths to reach the absorbing boundary longer.The second effect is that the puncture allows the walker more direct movement towards the bottom layers leading to smaller MET.The removal of some of the impenetrable barriers allows for more direct paths to the absorbing boundary, which leads to smaller mean exit times.The strength of this effect is dependent on the size of the puncture hw.For small values of h, the first effect has greater influence leading to an increase in the METs.As h is increased the second effect becomes more prominent and drives down the METs resulting in a global maximum.The interplay between the two effects also gives rise to the oscillations.Puncture of a brick layer opens it up, leading to larger exit times as the walker becomes temporarily confined inside a brick.Increasing the puncture height further destroys the brick structure of a layer and allows the walker to traverse the latter via a direct route thereby decreasing the exit times.
The global maximum and the oscillations are only present when the barriers are highly reflecting or impenetrable i.e. 0 λ v,u λ u,v ≤ q/2 for all {u, v} ∈ S. The maximum is lost when the permeability gets larger as the random walker is only partially confined by the barriers, leading to a monotonic decrease in the MET as seen in the inset of Fig. 11.With permeable barriers all the sites are always accessible independently of h and w, puncturing only creates more direct routes to the absorbing boundary leading to smaller exit times.

VI. THIGMOTAXIS
For the second application we look at thigmotaxis, which broadly speaking, is the movement of an organism due to a touch stimulus.We are interested specifically in the tendency of animals to remain close to the walls of an environment, a behaviour that is observed in many species from insects to mammals [68,69].We quantify the thigomotactic tendency by appropriately selecting defects location and λ to represent regions which are more easily accessible when moving in one direction (approaching boundaries) versus another (moving away from boundaries).
Since we are able to construct arbitrary shapes with the formalism, we consider two concentric circles within a square domain.The first is used to restrict the walker to a circular reflecting domain of radius R. The second has a radius r, with r < R, and is used to partition the domain into two regions: an inner region; and an outer region, which is the annulus between r and R, representing the preferred area of occupation.By placing one-way partially reflecting barriers along the radius r, we allow the walker to leave the inner region to enter the outer region without any resistance, while the tendency of remaining in the outer region is controlled by the parameter α i ∈ [0, 1].With α i = 1 the walker never leaves the outer region once it gets there, whereas with α i = 0, the partially reflecting barrier are removed and all areas of the circular domain becomes easily accessible.For a details on the placement of the defects and the construction of the circular domain see Appendix C 1. 13. Saturated mean-squared displacement for a thigmotaxis process for different values of the ratios of r/R.We study the dynamics as a function of the normalised parameter αi ∈ [0, 1] which represents the tendency of the walker to remain close to the boundary.When αi = 0, there are no outer or inner regions, while with αi = 1, the walkers never leave the outer region once they they get there.The saturation MSD is normalised by M, which is the saturation value when αi = 0. Other parameters used are described in in the caption of Fig. 12.
Given these constraints we study the dynamics as a function of the α i .In Fig. 12, we plot the probability Φ n0 (n, t) for different values of t = 500, 750, 1000 and ∞.The walker is initially at the centre of the domain n 0 = (51, 51) and can freely move inside the inner-region and is able to enter into the outer region without any resistance.However, once inside the outer region there is a greater tendency not to leave, due to the high value of α i = 0.95.We observe this effect when going from panel (a) to (d).Initially the separation between the inner and outer regions are barely visible but as time progresses this separation becomes increasingly clear, culminating with a sharp step at the steady-state.In panel (e) we plot a cross section of the probability at n 2 = 51 for the times corresponding with panels (a)-(d).
To examine the system further, we plot in Fig. 13 the mean-squared displacement (MSD) at steady-state, M, as a function of α i , for four different ratios of inner and outer regions.The MSD at steady-state is given by M With the curves normalised to the case where there are no internal barriers, M, i.e. when α i = 0. We find that as we increase α i from zero, for small values of α i , M initially increases logarithmically, while further increases of α i causes M to saturate.The value of saturation is dependent on the ratio of r/R: with a high ratio the outer region is thinner keeping the walker closer to the boundary and yielding greater saturation value, whereas a smaller ratio results in a thicker outer region allowing the walker to remain closer to the initial condition leading to a smaller value of M. Note that the reason for the r/R = 0.92 curve not being on top of the others is due to the discretisation of space when the outer region is very thin.

VII. TWO PARTICLE COALESCING PROCESS
In this final example, we demonstrate the use of our frame to model certain inert interactions between particles.The interactions we consider are partial mutual exclusion and reversible binding, both of which play an important role in coalescing dynamics.
14. Schematic representation of a two-particle coalescing process, modelled as one dimensional interacting random walkers.Panel (a) depicts the dynamics of particles A and B where q ∈ (0, 1] is the probability of moving at each time step.The combined dynamics of A and B can be represented via one next-nearest random walker in a 2D domain.This abstract domain is depicted in panel (b).The red circles represent locations where the two particles are on different sites, while the blue circles along the right diagonal are locations where they are co-located.In this space the interaction of A and B is modelled with partially reflecting barriers.The placement of these barriers are illustrated by the solid, dashed, and blue-dashed lines, the precise locations and permeability are given in Appendix C 2. The solid black lines are heterogeneities used to model the binding interactions, while the dashed black lines are used to control the unbinding interactions.The movement of C is represented by the 2D random walker moving along the diagonal.Its movement is slowed down, relative to A and B, through the placement of partially reflecting barriers along the diagonal depicted by the dashed blue lines.The resulting movement dynamics of the complex C is shown in panel (c), where αc ∈ [0, 1] represents the degree with which the movement of the complex C is slowed down.
Coalescing processes are ubiquitous in biology and chemistry; they consist of two or more entities that interact to bind and form a new one with different movement characteristics.An example of a coalescing process is the search of a promoter region on DNA by transcription factors.These movement dynamics alternates between periods of 3D search in the cytoplasm and periods of restricted search along the 1D DNA [70].We use our framework to study a system of relevance to the latter scenario: a first-passage process of two interacting particle in 1D.
We consider two particles labelled A and B that move independently on a 1D lattice with reflecting boundary conditions (see Fig. 14 for a schematic representation of the process).Their combined dynamics is described by a two dimensional next-nearest propagator ϕ n0 (n, t), with n 0 = (n 01 , n 02 ) and n = (n 1 , n 2 ).It represents the probability that the particle A and B are located, respectively, on the site n 1 and n 2 at time t given that they started, respectively, on n 01 , and n 02 .Two particles instantaneously form the complex C, namely when they encounter each other, that is when n = (m, m) for The interactions between particles is modelled through the placement of heterogeneities on the combined 2D lattice, yielding three control parameters, α e ∈ [0, 1], α u ∈ [0, 1], and α c ∈ [0, 1] (see Appendix C 2 for details regarding the placement of the defects).These parameters are used to constrain, respectively, the binding events via mutual exclusion of A and B, the unbinding events of C and the mobility of C. The parameter α u is proportional to the unbinding probabilities, while α e is proportional to the mutual exclusion probability.When α u = 1 and α e = 0, there is no interaction between the two particles.The other extreme represents strong interaction: when α e = 1 there is mutual exclusion, whereas α u = 0 results in a binding that is irreversible.The parameter α c ∈ [0, 1] represents the fraction of the movement probability of complex C relative to the movement probability of the constituent particles A and B. When α c = 1 there is no slowing down, while α c = 0 results in an immobile C.
In Fig. 15 we plot the log ratios of the MFPT, F n0→n for both particles to reach a site at the same time, compared with the 2D homogeneous next-nearest neighbour analogue, F n0→n .The latter corresponds with the case when α e = 0 and α u = α c = 1.The panels (a)-(d) depict F n0→n for increasing values of α e .The smallest ratios are observed in the upper left quadrant, which corresponds with high cohesiveness of the complex C and with only a slight reduction to its mobility, given respectively by, low values of α u and high values of α c .Within this parameter region, once the two particles bind they rarely separate, consequently the search in 2D reduces to a search in 1D with fewer sites to explore leading to smaller F n0→n .
When there is no exclusion interaction, i.e. panel (a), the dynamics of a similar model was explored in Ref. [57].In their analysis using asymptotics and simulations, equivalent features were observed.The most prominent feature of those and our observations is the minimisation of the MFPT for a slow moving C. In this regime, it is more favourable to have an intermediate unbinding probability, allowing the two particles to travel independently towards the target before recombining and hitting the target.
The ability to explore easily the parameter space of the model allows us to analyse the MFPT for different values of α e .By comparing the four panels we observe that as α e increases, the overall magnitude of the MFPT ratio decreases.This is explained by the fact that for small and intermediate values of α e the 2D walker is partially restricted to the upper or lower triangular regions of the domain, thereby reducing the overall exploratory space resulting in shorter search times.However, if α e is increased further, i.e. when 0 α e < 1, the particles will rarely coalesce, and the MFPT increases.In other words, shorter MFPTs can be achieved by having particles that mutually exclude one another with some probability.The ratio of the MFPT (F n 0 →n ) of the coalescing system compared to the MFPT (Fn 0 →n) of a homogeneous 2D next-nearest neighbour walker as a function of the heterogeneity strength parameters (see Fig. 14 for detailed a description of the parameters involved).The reactive site is located at n = (100, 100) and the two particles are initially maximal distance away from each other, i.e. n0 = (1, 100) and a target location n = (100, 100), on combined 2D domain of size N = (100, 100) with diffusion parameter q = 2/3.

VIII. CONCLUSION
We have introduced an analytical framework to model explicitly any inert particle-environment interactions.We have constructed the discrete Master equation that describes the spatiotemporal dynamics of diffusing particles in disordered environments by representing the interactions as perturbed transition dynamics between lattice sites.To solve this Master equation we have extended the defect technique to yield the generating function of the propagator in closed form.Using the propagator, we have derived useful quantities in the context of transport processes, namely, first-passage, return and exit probabilities and their respective means.We have also uncovered the existence of a disorder indifference phenomenon of the mean first-passage time in quasi 1D systems.
In light of the relevance of our framework to empirical scenarios, we have chosen three examples to demonstrate its applicability of our theory.In the first example we consider transdermal drug delivery, an intercelluar trans-port process, where we represent the 'brick-and-mortar' structure of the stratum corneum with the placement of reflecting and partially reflecting barriers.This representation allows us to study the effect that piercing has on the traversal time of a drug molecule.In the second example, we have examined the effect that an animal's thigomotactic response has on the mean squared displacement at log times.Lastly, in our third example, we have highlighted the ability of our formalism to study inert interactions between particles.We transformed these interactions and the ensuing dynamics into a single particle moving and interacting with quenched disorder in a higher dimensional space.The setup allows us to model analytically the search statistics in a two particle coalescing process, akin to the search of binding sites on the DNA by multiple transcription factors.
The strength of our result is in deriving the propagator in the presence of spatial heterogeneities, Φ n0 (n, z), as a function of the homogeneous propagator, i.e. the propagator in the absence of heterogeneities, ϕ n0 (n, z).This modularity allows one to change the movement dynamics by selecting different forms of ϕ n0 (n, z).In place of the diffusive propagator one may employ a biased lattice random walk [71], or a walk in different topologies such as triangular lattices [72,73], Bethe lattices [74,75] or more generally a network [76].
The modularity carries through to the heterogeneous propagator.This means that in situations where homogeneous space is assumed, one can relax this assumption and replace the homogeneous propagator, ϕ n0 (n, z) with the heterogeneous counterpart Φ n0 (n, z).We have demonstrated this aspect by studying the first-passage probability to either of two targets using results previously derived considering a homogeneous lattice.Further theoretical exploration could include the analysis of cover time statistics [77,78], transmission dynamics [79,80], resetting walks [81][82][83], mortal walks [84], or random walks with internal degrees of freedom [85].
Directions for future applications span across spatial and temporal scales: the role of a building geometry or floor plan on infection dynamics in hospital wards and supermarkets [86][87][88]; the prediction of search pattern behaviour of animals in different types of vegetation cover [89,90]; the heat transfer through layers of skin with differing thermal properties [91]; and the influence of topological defects on the diffusive properties in crystals [92,93] and territorial systems [94][95][96].
We conclude by drawing the reader's attention to the following.As experimental technologies continue to evolve, observations of the dynamics of particleenvironment interactions are increasing in number and resolution.The detailed description of the environment that these technologies bring presents a unique opportunity to rethink modelling techniques, moving away from macroscopic paradigms to a more microscopic prescription.We believe that the mathematical framework we have introduced to quantify the particle-environment interactions will play a crucial role in connecting the mi-croscopic dynamics to the macroscopic patterns observed across a vast array of systems.

Mean Return Time
Through the renewal equation we also have the return probability relation By noticing the identical structures of Eqs.(A1) and (A7), one can use a similar procedure to the one used to derive the MFPT (see Section II C of the Supplementary Materials) to show the mean return time (MRT) to be 2) .(A8)

Mean Exit Times
The first-exit probability is given by , where S n0 (z) is the survival probability given by Substituting Eq. ( 5) into Eq.(A9) and evaluating the sum in n and simplifying the summation over k one finds where S n0 (z) = n ϕ n0 (n, z) is the homogeneous survival probability, and where the elements of H are given in Eq. ( 6) and By taking the mean of the first-exit distribution, i.e. gives Eq. (10).The mathematical details to derive Eq. ( 10) are given in Appendix A 2 while simple expressions of the 1D problem are given in Section IV of the Supplementary materials.
Appendix B: First-passage quantities in 1D systems

The MFPT disorder indifference phenomenon
We start with a heterogeneous lattice reflecting boundary between n = 0 and n = 1, and a partially reflecting barrier between u and u + 1, with u < n 0 < n as depicted in Fig. B1.The trajectories that contribute to the firstpassage probability can be split into mutually exclusive sets based on the number of return visits, m, to the initial site n 0 .We now formally represent the first-passage probability in terms of a set of mutually exclusive independent events.Let us define F n0 (n, t; m = 0) as the first-passage probability to reach n for the first time at t having started at n 0 and having never returned to the initial site.Clearly, the trajectories that make up F n0 (n, t; m = 0) (coloured blue in Fig. B1), can never be affected by the presence of the barrier as they never move towards the barrier.The trajectories that could be affected by the presence of the barrier are those that return at-least once to the initial site before reaching n.The first-passage probability to visit n and having visited the initial site m times is constructed through the convolution (dashed trajectories in Fig. B1) and where The function h n0 (n, t) represents the probability of returning without visiting the target and is constructed by considering the probability of returning to n 0 and subtracting those that reach n without returning to n 0 at some prior time and subsequently reaching n 0 from n.In z domain the relation can be written more conveniently as Notice that the only term with the dependence on the barrier on the RHS of Eq. (B3) is R(n 0 , z), and in the absence of the barrier R(n 0 , z) reduces to R(n 0 , z).The full first-passage probability for the system with the barrier, which can be written as the sum of the mutually exclusive probabilities giving The relation given by Eq. ( B4) is an alternative method of constructing the first-passage probability i.e. its not one of the standard approaches which are through the survival probability or the ratio of propagators in z-domain.
To confirm the normalisation of the RHS of Eq. (B4) consider the following.By definition F n0 (n, t; m = 0) is not normalised over t, hence, F n0 (n, z = 1; m = 0) = p where 0 < p < 1.Since all other terms in the RHS of Eq. (B4), namely, F n (n 0 , z) and R(n 0 , z) are normalised over time, we find that at z = 1 the RHS becomes ∞ m=0 p(1 − p) m = 1.Differentiating Eq. (B4) with respect to z and taking the limit z → 1, we obtain the first-passage time When the barrier is such that λ u,u+1 = λ u+1,u = λ , the mean return time is equal to the reciprocal of the steady state value, and R n0 becomes R n0 , i.e. the mean return time in the absence of the barrier.To find p in Eq. (B5) explicitly, we first construct F n0 (n, z; m = 0) in terms of known quantities using the approach presented in Ref. [58] to construct timedependent splitting probabilities.We write the two relations by considering the two splitting separately: the first-passage probability of reaching the target n and never returning to the initial condition F n0 (n, t; m = 0); and the first-return probability to n 0 and having never reached the target site n.In time domain they are written via a convolution and are, respectively, where R(n 0 , t; n) is the probability of returning to the site n 0 at t and having never visited the target n.One can take the z-transform and solve for F n0 (n, z; m = 0) and R(n 0 , t; n) giving, respectively, Evaluating Eqs.(B8) and (B9) at z = 1 gives, respectively, the fraction of all the first-passage trajectories that reach the target without returning to n 0 and the fraction of all trajectories that return to n 0 without ever reaching n.Using de L'Hôpital's rule once in Eq. (B8) we find Inserting Eq. (B10) in Eq. (B5) one finds that F n0→n = F n0→n .

MFPT linear dependence on disorder location
To understand the linear dependence in u, with n 0 ≤ u < n, present in Eq. ( 12), we consider building up first passage probability by convolution in time to go from n 0 to u first, then from u to u + 1 and then from u + 1 to n.In z domain one has where the first term on the RHS has no dependence on the barrier as it is after the absorbing site u, while the other two terms are dependent on the barrier.Computing the mean of Eq. (B11), we obtain where we have substituted F u+1→n = F u+1→n using the justification presented in the previous section.By using the relation F n0→n = F n0→s + F s→n with n 0 < s < n, one can rewrite Eq. (B12) to give In the diffusive case, the MFPT to a neighbouring site, F u→u+1 , is always proportional to twice the distance between u and the reflecting boundary to the left, i.e.
with s ≤ u being the position of the reflecting boundary.By simplifying the general MFPT given in Eq. ( 8), we find which is analogous to Eq. (B14) but with the multiplicative (time rescale) factor increased to 2(q−2λ) −1 .Letting s = 1 we obtain Eq. ( 12).Two sets of defects must be placed, one set along a circle radius R to create a (circular) reflecting domain, while the second is used to divide this domain into two different regions (see Section VI) and is placed along a circle of radius r.To place defects on either circle one must first know which sites are within which circle.To determine this use the Euclidean distance as a heuristic, with the site n = (n 1 , n 2 ) being part of the circular domain if and only 1/2 with the size of the bounding square domain given by N = (2R + 1, 2R + 1).
Similarly, a site is part of the inner region if and only if h(n 1 , n 2 ) ≤ r, while the outer region is given by r < h(n 1 , n 2 ) ≤ R. Given these site partitions, one can define two sets of defects, S d , and S i describing, respectively, the impenetrable barriers to restrict the walker to a circular domain, and partially-reflecting inner barriers.In both cases u represents sites inside the circle of defects while v represents sites outside.For all {u, v} ∈ S i we have λ v,u = A v,u and λ u,v is irrelevant as the walker initially starts inside the circular domain.For all {u, v} ∈ S d , we let λ v,u = 0 providing no resistance for the walker to enter the outer-region and λ with α i ∈ [0, 1].

Two particle coalescing process
The interactions that need to modelled are binding and unbinding.Binding can occur via two distinct events.The first is when two particles are located on neighbouring sites with n = (m + 1, m) and at the following time step, one of the particles remains at the same site while the second particle jumps onto the site occupied by the first resulting in n = (m + 1, m + 1) or n = (m, m).The second possible event occurs when the two particles are located two sites apart i.e. n = (m − 1, m + 1) and at the following time step they both jump towards each other landing on n = (m, m).The reverse of these two events gives rise to unbinding of the complex C.These transitions can be modified by placing paired defects of the forms: 4 α e .Intuitively, the movement of the coalesced is slowed as it more massive.To encode this detail we interpret jumps along the leading diagonal n = (m, m) for all 1 ≤ m ≤ N as the jumps made by the coalesced particle C, and we slow its movement by placing paired defects of the form u 4 (1 − α c ).

A. Sticky and Slippery Heterogeneities
We start with set of defective sites S = {w 1 , • • • , w L } and use the notation w (ls) i and w (rs) i representing, respectively, the left and right neighbours of w i in the s th -dimension, given w , with d the lattice dimension.A schematic representation of the jump probabilities on a defective site, w is given in Fig. A1 , from which it is clear that to ensure positive probabilities one must have w w (ls)  w (rs) A w (rs) ,w − λ w (rs) ,w A w (ls) ,w − λ w (ls) ,w A w,w + This formal solution is as special case of Eq. (S8) where M = 2Ld where each of the sites in S bears two paired defects for each of d dimensions.However, by noticing that the incoming connections of the sticky sites are left unmodified, i.e. λ w,w (ls) = λ w,w (rs ) = 0 one can to simplify Eq. (S8) to Eq. (S16) thereby reducing the number unknowns by a factor of 2d.
To find the full solution we let n = w k and solve the simultaneous equations where, in this case, the matrix H is simplified to and Y is the same as H but with the k th column replaced by ϕ n0 (w 1 , z), • • • , ϕ n0 (w M , z) T .Substituting Eq. (S19) into Eq.(S16) and summing over k gives the full solution in Eq. ( 5) where H(n, n 0 ) i,j = H i,j − Q wj (n, z) ϕ n0 (w i , z). (S21) where the factor 6 comes from repeated application of the product rule.To carry out the summation in Eq. (S50) we employ the following property of determinants.Given two matrices, A and B of size M × M , where B = ab T and where a and b are two column vectors, the following relation holds when A ( ) is the same as A, but with the th column replaced by the th column of B. Comparing Eqs.(S50) and (S53) with Eq. (S54), we observe that A → J − J (2) , J ( ) = (−1) M (2M + 1)! 3 J − J (2) − J − J (2) + 3J (1) .(S59) The factor 1/3 can be taken into the determinants using the property given in Eq. (S24).For the first determinant on the right hand side (RHS) of Eq. (S59) we can equate A → J , from Eq. (S43) we can equate to derive the relation Similarly, for the second determinant on the RHS of Eq. (S59), we can equate A → J , and using Eq.(S58) we identify J ( ) = (−1) M (2M + 1)! J − 1 /3J (2) − J − 1 /3J (2) + J (1) .(S66) Lastly, employing the property with A → J , 3a with Eq. (S60), 3b with Eq. (S61) and 3c with Eq. (S63), we can simplify Eq. (S66) to Putting it all together gives the final mean first-passage time with M paired defects J − J (2)  (S69) If we divide through all the terms by (u,v)∈S R u R v , we obtain where H (1)  i,j = H (2)  i,j = To build the first-passage probability with sticky and slippery heterogeneities (see Section I A and Fig. A1) we use the propagator given in Eq. ( 5), where the matrices H and H(n, n 0 ) are given, respectively, by Eqs.(S20) and (S21).As the procedure is similar to the one used to derive the general MFPT given by Eq. (S70) (and also Eq. ( 8) in the main text), we outline only the key steps.
Starting from

FIG. 1 .
FIG.1.Examples of the spatial heterogeneities within a square lattice of width 10 with reflecting boundaries.Panel (a) represents an open partitions with the solid black lines indicating impenetrable barriers.When these barriers enclose a region, some space becomes inaccessible indicated by the sites coloured dark grey in panel (b).Panel (c) shows a lattice where three pairs of non-neighbouring sites have a long range connection, i.e. transitions from a dashed site include the nearest-neighbours as well as the site connected via the dashed line.Panel (d) is an example of where the diffusivity of the striped sites is smaller than the regular (non-striped) sites.The central site flagged by a red square and the two sites flagged by a blue diamond are, respectively, the initial condition and the absorbing targets for use in later sections.

FIG. 5 .
FIG.5.Example of a sticky site on a 1D lattice generated by reducing all of the outgoing probability to the neighbours as shown by the thinner arrows, whilst increasing the staying probability of w as shown by the thicker self loop (from left to right).

FIG. 8 .
FIG.8.The ratio of the heterogeneous mean exit time E n 0 to the homogeneous mean exit time E n 0 for randomly distributed barriers and anti-barriers.We use a homogeneous propagator with a domain of size N = (51, 51) with absorbing boundary conditions, an initial condition at the centre of the domain n0 = (26, 26) and a diffusion parameter of q = (0.8, 0.8).The explicit form of ϕ n 0 (n, z) is given by the z-transform of Eq. (23) of Ref.[58].Each curve is obtained using Eq.(10) and performing an ensemble average with 10 2 sample realisations of locations of barriers (λ > 0) or anti-barriers (λ > 0) for each λ.
multidimensional domains.B. First-passage processes in 1D with a single barrier and the phenomenon of disorder indifference of the MFPT

FnFIG. 9 .
FIG.9.The time dependent first passage probability given through the numerical inversion of Eq.(11).Panel (a) represents the scenario u < n0 < n, while panel (b) is the case when n0 ≤ u < n.The values of the other parameters are: λ = 0.975 • q 2 ; a diffusion parameter of value q = 2 3 ; the initial condition n0 = 8; the target site n = 15; and a reflecting boundary between n = 0 and n = 1.The arrows indicate the MFPTs: in panel (a) all of the curves have the same MFPT of F n 0 →n = 231 (disorder indifference), whereas in panel (b) the two arrows indicate the minimum and maximum MFPTs of the curves which are F n 0 →n = 1167 and F n 0 →n = 1869, and obtained, respectively, when u = 8 and u = 14

FIG. 12 .
FIG.12.The propagator Φn 0 (n, t) at different moments in time, t = 500, 750, 1000, ∞ where the walker is initially at the centre of the domain, n0 =(51,51).When inside the inner region the walker can freely enter the outer region without any resistance, that is λv,u = 0 and when in the outer region the probability to move inward is modified via λu,v = αiA v,u .Other parameters used are the diffusivity of value q = (0.8, 0.8) and a square domain of size N = (101, 101), (see Appendix C 1 for details on the placement of defects) FIG.13.Saturated mean-squared displacement for a thigmotaxis process for different values of the ratios of r/R.We study the dynamics as a function of the normalised parameter αi ∈ [0, 1] which represents the tendency of the walker to remain close to the boundary.When αi = 0, there are no outer or inner regions, while with αi = 1, the walkers never leave the outer region once they they get there.The saturation MSD is normalised by M, which is the saturation value when αi = 0. Other parameters used are described in in the caption of Fig.12.
FIG.15.The ratio of the MFPT (F n 0 →n ) of the coalescing system compared to the MFPT (Fn 0 →n) of a homogeneous 2D next-nearest neighbour walker as a function of the heterogeneity strength parameters (see Fig.14for detailed a description of the parameters involved).The reactive site is located at n = (100, 100) and the two particles are initially maximal distance away from each other, i.e. n0 = (1, 100) and a target location n = (100, 100), on combined 2D domain of size N = (100, 100) with diffusion parameter q = 2/3.From panel (a) to (d) we have, respectively, the parameters αe = 0, 0.5, 0.75 and 0.875.
FIG.B1.A schematic representation of a one dimensional heterogeneous lattice with a reflecting boundary to the left (vertical line) and with a permeable barrier between the sites u and u + 1 represented by the shaded rectangle.The firstpassage event can be split into mutually exclusive events represented by arrows of different colours.The blue arrows represent trajectories that never return to the initial site, while the black ones represent trajectories that return m times before reaching n.The green arrow represents first-passage trajectories that reach n0 having starting at n.The solid arrows represent trajectories that are unaffected by the presence of the partially reflecting barrier between u and u + 1, while the trajectories that are affected are represented by dashed arrows.Note that this schematic depicts the case when 1 ≤ u ≤ n0 −1

d
FIG. A1.A schematic representationshowing the modified transition probabilities of a sticky (or slippery) heterogeneity.We only the dynamics in the s th dimension, but the same is present for the other dimensions.λw (rs ) ,w ≤ A w (rs) ,w and λ w (ls ) ,w ≤ A w (ls) ,w for all s = 1, • • • , d, and 0 ≤ A w,w + d s=1 λ w (rs ) ,w + λ w (ls) ,w for all w ∈ S .These conditions are a recast of the one given by Eqs.(2) and (3) in the main text.The full dynamics is described by the Master equation