Generalized master equation for first-passage problems in partitioned spaces

Motivated by a range of biological applications related to the transport of molecules in cells, we present a modular framework to treat first-passage problems for diffusion in partitioned spaces. The spatial domains can differ with respect to their diffusivity, geometry, and dimensionality, but can also refer to transport modes alternating between diffusive, driven, or anomalous motion. The approach relies on a coarse-graining of the motion by dissecting the trajectories on domain boundaries or when the mode of transport changes, yielding a small set of states. The time evolution of the reduced model follows a generalized master equation, which takes the form of a set of linear integro-differential equations in the occupation probabilities of the states and the corresponding probability fluxes. Further building blocks of the model are partial first-passage time (FPT) densities, which encode the transport behavior in each domain or state. The approach is exemplified and validated for a target search problem with two domains in one- and three-dimensional space, first by exactly reproducing known results for an artificially divided, homogeneous space, and second by considering the situation of domains with distinct diffusivities. Analytical solutions for the FPT densities are given in Laplace domain and are complemented by numerical backtransform yielding FPT densities over many decades in time, confirming that the geometry and heterogeneity of the space can introduce additional characteristic timescales.


Introduction
Cellular spaces and membranes show heterogeneous structures due to compartmentalization and macromolecular crowding, leading to complex diffusive transport of molecules with implications for biochemical reactions [1][2][3].The phenomenon of anomalous (sub-)diffusion plays certainly a prominent role here, which is typically more pronounced for large molecules and long-distance transport and which is likely to subsume a number of physical causes.Yet, already concrete microscopic structures can give rise to non-trivial dynamics, examples being spatial domains of varying diffusivity [4], possibly separated by the nuclear envelope or other diffusion barriers [5], and the nowadays established nano-scale partitioning of the plasma membrane [6][7][8].Another aspect are alternating modes of motion, e.g., stochastic switching between actively directed and Brownian motion [9], or between different dimensionalities of space such as sliding on one-dimensional (1D) DNA strands and 3D diffusion in the nucleoplasm [10]; in the context of nanocatalysts, one finds surface diffusion on a nanoparticle interleaved with 3D diffusion in solution [11].
First-passage times (FPT) are a crucial factor in the chemical kinetics of diffusioninfluenced reactions: they measure the time of first encounter between reaction partners diffusing in space.For high copy numbers of the reactants and under well-mixed conditions, the kinetics is appropriately described in terms of reaction rate constants, which are essentially the reciprocal of the mean FPT.If these requirements are not fulfilled, as in femtolitre-sized cells when signaling molecules are not abundant, the mean FPT is insufficient to model the kinetics and it appears advisable to consider the whole probability distribution of FPTs instead [12,13].The influence of the geometry on FPT distributions and thus the reaction kinetics was appreciated only recently [12][13][14].Already for comparably simple geometries, one observes FPT distributions that are governed by two or more widely distinct time scales, meaning that the FPTs can fluctuate considerably from one molecular trajectory to the other [14,15].
The aim of the present paper is a mathematical framework of first-passage problems for diffusion in a continuous space partitioned into domains.These domains are assumed to be homogeneous regions, which however can differ from each other with respect to their transport properties or even their dimensionality; at a later stage, different chemical reactions may occur within each domain.Exploiting the Markov property of (idealized) diffusion, we coarse-grain the process and dissect a given trajectory into parts whenever a domain boundary is crossed (Fig. 1).Following the dynamics from one crossing to the next, we have recast the diffusion problem as a renewal process.These partial trajectories are fully contained within a single domain, and all possible trajectories within one domain are identified with a coarse-grained state.
Therewith, the original diffusion process in continuous space has been replaced by a continuous-time random walk (CTRW) on the domain states.The waiting times between the jumps correspond to the dwell times (or residence times) on the respective partial trajectories and their distributions are given as solutions to (partial) first-passage problems on each domain.The jumps between domains constitute in general not a Poisson process (as for a Markovian random walk), not even for simple diffusion, and thus introduce a memory into the evolution equations of the occupation probabilities of the states.We will refer to the latter set of equations as the generalized master equation (GME) of the coarse-grained problem.
In our approach, we do explicitly not require any mechanism of barrier crossing between the domains, although our approach allows for such barriers.A computationally motivated example for a situation without barriers is Doi's volume reaction model [16], used in particle-based reaction-diffusion simulations [17,18].In the presence of sufficiently high barriers, each domain forms a metastable set with respect to a molecule's motion, transport on long time scales resembles a Markovian hopping process, and the reaction kinetics in such a partitioned space can be described by a spatio-temporal chemical master equation [19].Here, we lay a basis to go beyond these approximations.
Early studies of CTRWs on a finite state space date back almost half a century ago and include a two-state model for the orientational motion of molecules in dense media [20], later amended by a CTRW in space on top of it [21]; a more recent application are the on and off times in blinking quantum dots [22].In all these examples, the sought quantities were computed using classical renewal equations; introductory texts on this method can be found in Refs.[23,24].
In the present paper, we will adopt an alternative, more heuristic approach put forward by Chechkin et al. [25] in the context of a different problem, where the authors also coined the term "generalized master equation" (GME).This approach was applied successfully in the modeling of reaction-subdiffusion systems [26,27].The charme of the method lies in its clarity as it can be derived from a few basic principles such as .Geometries of the first-passage problem with two domains and sketches of exemplary partial trajectories in a one-dimensional half axis (left) and between concentric spheres in three dimensions (right); the larger sphere with radius  > |r0| is not shown.Space is partitioned at the position/radius  into an inner (red) and an outer (blue) domain.The symbols  denote the dwell time densities of the different types of partial trajectories.
local balance and continuity of probability fluxes.An adaptation to more complex situations is relatively straightforward, whereas in such cases the renewal equation framework would become increasingly convoluted and nearly impossible to handle, let alone the explicit calculation of fluxes onto the target via the partial differential equation (PDE) for the transport propagator.
Here, we will adhere to the very simple case of two domains either with completely the same physical properties (i.e.no inhomogeneity at all) in order to establish the method, or with two different diffusion constants.For the latter case we will compare our results to those in Ref. [15], where the FPT density was deduced from the flux onto the target by solving the PDE for the particle density for the whole, non-uniform system.

Formulation of the problem
As a test bed of the proposed framework, we apply the GME method to the calculation of the FPT distribution for target search by free diffusion (i) on a 1D half axis with an absorbing boundary at position  > 0, and (ii) in the 3D space between two concentric spheres of radii  < , where the boundary of the inner sphere is absorbing and the outer one is reflecting (Fig. 1).A domain boundary is placed at the position  >  (1D) and at radius  with  >  >  (3D), respectively, and we refer to the region between  and  as inner domain and to the remaining accessible space as outer domain.The trajectories start in the outer domain at x 0 with |x 0 | > .
(For simplicity, we use the vector notation also in 1D space.)The central quantity of interest is the probability density  FPT () of the random time  until the first encounter with the boundary at .
As outlined above, we dissect the trajectories into partial ones whenever the boundary at |x| =  is crossed.The dwell time on a given partial trajectory (equivalently, in its respective domain) is the time the molecule spends between entrance to and exit from the domain.The probability density of dwell times is an FPT density itself, which is why it will be referred to as a partial FPT density in the following.In the outer domain, there are two types of trajectories: starting at x 0 and starting at the boundary |x| = , both types end at this boundary; the corresponding partial FPT densities are  0 () and  out ().Trajectories in the inner domain always start at |x| = , but either leave the domain at this boundary again or are stopped at |x| =  with partial FPT densities  in () and  in,∅ (), respectively.The specific forms of these dwell time densities for the geometries chosen here are discussed in Sections 3.1 and 4.1.
In a first step, we assume no physical difference between the inner and outer domains so that the boundary at |x| =  is only an artificial one.This allows us to check the proposed GME approach of dissecting trajectories and assembling the overall FPT density  FPT () for reaching the boundary |x| =  from the partial ones.The result must reproduce the known distribution of FPTs for reaching |x| =  directly when starting at x 0 .In Section 5, we elaborate on the situation of different diffusion coefficients in the two domains.

Generalized master equation: simplified model
In the two-domain model of the preceding section, the state of the molecule is described by the probabilities  in () and  out () that at time  it is found in the inner and outer domain, respectively.It is amended by the probability  ∅ () that the trajectory was stopped at time  or earlier: the molecule has reached the target and was removed from the system, e.g., by a chemical reaction.The time evolution of the tuple ( in (),  out (),  ∅ ()) is governed by a non-Markovian generalization of the master equation (GME).
The derivation of the GME of our problem follows ideas in Ref. [25] and makes use of two types of conservation laws: (i) probability is conserved locally: net gains or losses within one state determine the change of local probability (ii) probability is conserved in the transitions between the states (continuity of the fluxes).The other key ingredient is an equation accounting for the renewal character of the stochastic process, linking present losses from a state with the gains at a previous time.
The transitions between the states induce probability fluxes  ± in/out , where superscripts denote a gain (+) or loss (−) of probability in each state.The subscript is out for loss from or gain to the out state, and in for overall gains or losses of the in state.Losses from in to the "absorbed" state (∅) are not included in  − in ; these are written as  − in,∅ .Figure 2 shows a transition graph of the model along with the loss fluxes.
For an unbiased diffusion, a partial trajectory ending at the domain boundary |x| =  is continued in either the inner or the outer domain with equal probability 1/2.By the Markov property of diffusion, this applies for partial trajectories irrespective of which side of the boundary they belong to (in or out state).Therefore, continuity of the fluxes in the transitions requires: where the terms on the right represent the incoming arrows of the state given on the left as depicted in Fig. 2.
Local balance in each state demands: that is, the temporal change of probability in a state is the difference between the total gain and loss fluxes.Summing up the three equations shows that the overall probability, including that of the absorbed particles, is conserved as it should: (d/d) [ in () +  out () +  ∅ ()] = 0.The quantity  ∅ () in Eq. (2c) collects the trajectories that have stopped up to time .Thus, the change of  ∅ () is actually equal to the sought FPT density of the target search problem, d ∅ ()/d =  FPT (), and the task is to compute the flux  − in,∅ () onto the boundary |x| = .Finally, we write the following recursive, renewal-like relations for the loss fluxes Invoking the picture of an ensemble of particles, these equations reflect the fact that a loss of particles from a domain can only happen if the particles had been there before, either from the very beginning as encoded in the term  0 () or through previous gains: the integral terms describe particles that entered the domain at an earlier time  −  ′ and leave it at .Note that the fluxes  − in and  − in,∅ are mutually dependent since partial trajectories of the in state can either end at the boundary |x| =  (subscript in) or at |x| =  (subscript in, ∅, absorption), see Fig. 1.Therefore,  in and  in,∅ are no proper FPT densities in the sense that they do not normalize.The full density of dwell times in the state in is given by their sum, which satisfies ]︀ d = 1.This fact expresses the splitting of probabilities to follow one or the other type of partial trajectory (i.e., to survive or to be stopped at the end of the current step).
The system of Eqs.(2a) to (2c) and (3a) to (3c), together with Eqs.(1a) and (1b), forms the GME of the first-passage problem with two domains.Using Eqs.(1a) and (1b), the gain fluxes  + in/out can be expressed in terms of loss fluxes, which leaves us with a system of 6 linear integro-differential equations in 6 variables (3 loss fluxes and 3 occupation probabilities).This type of equation system is conveniently solved in Laplace domain.

Solution in Laplace domain and numerical backtransform
The Laplace transform f := L[ ] of a measurable function  () on R + is defined as [23] f where the frequency  > 0 is the Laplace variable conjugate to .For a convolution , and for the time derivative one has L [︀ ḟ ]︀ () =  L[ ]() −  (0).Laplace transformation of the self-consistent, linear system of integro-differential Eqs.(2a) to (2c) and (3a) to (3c), yields a closed set of linear, algebraic equations in the probabilities and fluxes with the partial FPT densities as coefficients.Substituting  + in/out in Eqs.(3a) to (3c) by means of Eqs.(2a) and (2b), we obtain: where we made use of the initial conditions  out (0) = 1 and  in (0) =  in,∅ (0) = 0. Solving for the loss fluxes, we have: The system of equations in Laplace domain is completed by from Eqs. (2a) to (2c).Solving the linear system for j− in,∅ () provides us with an explicit expression for the sought FPT density in Laplace domain, pFPT () = j− in,∅ (), which is fully specified by the partial FPT densities φ ():  probability density and thus integrable.It allows for the analytic continuation to the upper complex plane and is connected to the Laplace transform by φ() = φ(i).The backtransform is uniquely given by a cosine transform: using that Re φ() is an even function in .For the robust numerical evaluation of the Fourier integral, we used a modified Filon quadrature as developed recently for the back and forth transformation between time correlation functions and their dissipation spectra [28].

Regularized GME with an auxiliary boundary
The coefficients in Eqs.(6a) to (6c) become singular if any of the FPT densities φ () = 1, i.e., if   () =  + () is the density of the Dirac measure on the positive reals, R 0 , supported at  = 0. Unfortunately, we are facing this problem for  in () and  out () with the setup of Fig. 1: partial trajectories starting at the boundary |x| =  return to that boundary immediately, almost surely, the starting and ending points are the same.This issue is familiar from the first-return problem for diffusion in the continuum, it does not arise for random walks on a lattice.As a regularization of these singularities, we require a minimum length  > 0 of the partial trajectories in the calculation of partial FPT distributions and we take the limit  → 0 for the final result of the FPT distribution.Technically, we introduce an auxiliary boundary at |x| =  + , with  +  <  in the 3D case (Fig. 3).The boundary is transparent for partial trajectories in the out state, i.e., they start in the outer domain and end at the boundary |x| = .Partial trajectories starting at |x| =  are assigned to the in state, they either end at |x| =  +  or are stopped at the absorbing boundary |x| = .In the former case, the following partial trajectory belongs to the out state as it begins in the outer domain, it ends at |x| = .Thus, the out state is followed by an in state again and so forth.(e.g., diffusion coefficient) along a partial trajectory are those of the assigned state (or domain), which leads to a hybrid character of the region  < |x| <  + .Clearly, this interpretation is an approximation to the original diffusion problem, which is restored in the limit  → 0.
Our definition of the occupation probabilities   () of the states  ∈ {in, out, ∅} remains unchanged.And as before, the probability fluxes  ±  () denote gains (+) and losses (−) of the state ; the loss  − in,∅ from the in state to the absorbed state ∅ is not included in  − in ().The transition graph of the amended two-domain GME is given in Fig. 4. As already discussed and in contrast to the simplified setup, the transitions between the in and out states are strictly alternating, so that the loss fluxes do not split and flux continuity demands: which replaces Eqs.(1a) and (1b).For the local balance in each state, the same Eqs.(2a) to (2c) as previously hold.Also the Eqs.(3a) to (3c) for the fluxes apply without modifications.This set of equations constitutes the GME of the regularized two-domain model.
In the Laplace domain, the three expressions for the loss fluxes, Eqs.(6a) to (6c), carry over as well since their derivation did not rely on Eqs.(1a) and (1b).The linear system is completed by three equations for the occupation probabilities, which follow from Eqs. (2a) to (2c) and (10): Solving the linear system in six variables, we find the desired FPT density pFPT () = j− in,∅ () in terms of the partial FPT densities φ (): which is a main result of this work.Note that the partial FPT densities [except φ0 ()] implicitly depend on the regularization parameter .The complete solution for all probabilities and fluxes is given in Eqs.(A2) of the appendix.

FPT densities for partial trajectories
As mentioned above, the dwell time probability densities   () are themselves FPT densities, namely of the first passage from their entrance to the exit from the respective regions.We will refer to these regions as "outer" and "inner" regions, respectively (see Fig. 3).These partial FPT densities can be obtained from the corresponding partial differential equations (PDEs) for the positional particle (probability) density, with an initial condition and the respective boundary conditions.More precisely, the setup where particles reaching some point in space for the first time are immediately removed from the system, corresponds to the setup in terms of position probability density with an absorbing boundary at that point.The flux into that point x  will be the FPT density to visit that point for the first time, where the particles started at x = x 0 and x Ω is a point at the boundary.
For each region we will solve the respective PDEs in Laplace domain where they take the form of ordinary differential equations linear in x.The frequency  will be kept as a variable, although it takes the role of a parameter during calculations.In a first step, we will solve the pertinent initial-boundary-value-problems more generally for arbitrary starting points  0 or  0 and lower and upper boundaries at  ℓ ,   or  ℓ ,   , respectively and customize them later.A more detailed exposition of the procedure than we are able to give here can be found in [29].
Specializing to 1D spaces, we have the following equation governing the probability to find a particle at some where the right side represents the initial condition of all particles starting at  0 .The fundamental system is , i.e. our solutions will be a linear combination of these two exponentials in  and √︀ /.

Outer region
In the outer region we have an absorbing lower boundary at some  ℓ and a reflecting upper one at   : for which solves Eq. ( 15).The fluxes onto the boundaries give us the FPT densities for a particle to reach the respective boundary.Thus with φ ( ℓ ,   , ;  0 ) = ±   ψ(, ) (+ for the flux onto the lower and − for the flux onto the upper boundary) we have as expected.To see what happens in an infinite system, we let   → ∞: which is the known density for first passage to a point  ℓ on an infinite domain.An expansion in small  yields a non-analytic expression which gives rise to long time tail.(In fact the exponential of a square root of  gives the one sided Lévy stable law   of parameter  = 1/2 in , thus the tail is ∝  −3/2 .)

Inner region
Here we consider the boundary condition ψ( ℓ , ) = 0 (21) ψ(  , ) = 0 (22) for which the solution to Eq. ( 15) yields Again we compute the fluxes onto the boundaries [cf.Eq. ( 14)]: These fluxes are the FPT densities to the respective boundaries.The splitting probabilities, i.e. the cumulative probabilities to leave the system via the respective boundary, are obtained by sending  → 0, which corresponds to  → ∞ in time domain: Sending the outer boundary to infinity, we end up again with the expected FPT density on an infinite domain.

Scaling form and full solution for the FPT density
Let now, more in accordance with Fig. 1, start the partial trajectory at a small distance  to its end point , thus  0 =  + ,  ℓ =  for φout and  0 = ,   =  +  for φin .The lower boundary for the   in be  and the upper boundary of the  ℓ out be .
Furthermore let us introduce the scaled Laplace variable  :=  2 /.Then, we have All of these expressions for the FPT densities of the partial trajectories attain the same limit for  → 0: lim which permits an analytic inversion of the Laplace transform: with  * := / 2 the scaled time variable conjugate to .The partial FPT densities for the inner and outer domains are depicted in the left panels of Fig. 5.The closer to the absorbing boundary or target a particle starts, the shorter the time of the peak position and the higher the peak value, indicating that the particles are more rapidly absorbed.Thus, moving the starting position closer to the target corresponds to a limiting procedure for the partial FPT densities, lim →0  in/out (, ) =  + (), converging to a singular peak at  = 0: if the particle starts at the position of the target, all probability is absorbed immediately within zero time.Our construction of an auxiliary -shell around the boundary at , see 2.4, circumvents the difficulties arising from such singular FPT densities.At intermediate times, we find the expected  −3/2 power law decays for the particles exploring the outer domain yet without hitting the outer confinement.The FPT densities at large times decays rapidly due to the confinement.
For the inner domain, there is another reflecting boundary at |x| = .A sharp exponential cutoff sets in at times  ≈  2 / (Fig. 5, upper panels), which we attribute to partial trajectories that end as soon as they reach the outer boundary and do not contribute to the FPT density at longer times.The outer domain has a reflecting     27), upper panels] and the outer domain [out(), Eq. ( 28), lower panels] for the 1D problem.Left panels show results for the partial FPT densities in,out(), and right panels the scaled partial FPT densities.The parameters of the geometry are 0 = 10,  = 5, and  = 20, and  =  2 / is the unit of time.The black dotted line denotes the analytical solution in the limit  → 0 [Eq.( 30)].The gray dashed line indicates the asymptotic power law tail  −3/2 . .
boundary at |x| = , which also results in an exponential cutoff of the tail (Fig. 5, lower panels).However, the trajectories are continued and move on in their attempts to reach the target, which they will eventually do, but at a later time.This shifts and smears out the cutoff at large times.Moreover, due to the conservation of probability, the reflected trajectories, which in the case of an unbounded domain would have contributed into the power law tail, are responsible for the small shoulder of the FPT density at large times.The right panels of the figures show the partial FPT densities scaled with respect to the distance of the starting point to the target .Small values of  are equivalent to a large domain size  ≫ , the particle feels the confinement at a later time and the differences between reflecting or absorbing outer boundary vanish.

Full solution for the FPT density
With the partial FPT densities calculated above, we are ready to write up an analytical expression for the Laplace transform of the full FPT density for a particle starting at  0 and being absorbed at  by substituting (12).
As an ultimate test for the validity of the method we compare our GME solution for the FPT density for the homogeneous system (i.e.particles behave the same way in inner and outer region) with the FPT density for particles starting at  0 and ending at first encounter with  (i.e. for non-dissected trajectories).We find perfect agreement between the two of them, Observe that in this case the -dependence vanishes for the GME-solution, taking the limit  → 0 is not necessary anymore.Also the  vanishes from the expressions as it should in the completely homogeneous case.For infinite outer domains  → ∞ we have: In time domain, this limit is which is the expected result for first passage on an infinite domain.As a last step, expression (31) has to be Laplace-inverted which can be done numerically, see Fig. 6.It is indeed qualitatively the same picture as already for the partial FPT densities with reflecting outer boundary: for small times it fits the Lévy-Smirnow law, Eq. (33), at large times we get an exponential cutoff and at times just before the cutoff a small elevation compared to the  −3/2 decay, which collects the probability in the truncated tail.

FPT densities for the partial trajectory sections
Analogously to our calculations in Section 3.1, we now will calculate the partial FPT densities from the fluxes onto the boundaries of the respective partial differential Equations.for the radial symmetric setup (see Fig. 1b).Thus in this case, Eq. ( 13) becomes in Laplace domain: By substituting η(, ) := ψ(, )/, the equation for η attains a similar form as the one in the 1D case: where on the right hand side there is again the initial condition (all trajectories start initially from the radius  0 ).Again the fundamental system is

Outer region
For the outer region, we have the transformed boundary conditions (absorbing at  ℓ , reflecting at   ): so that the solution to Eq. ( 35) is given piecewiese for  ≶  0 by With ψ′ = η′ / − η/ (42) and confirms that with probability 1 the particles ultimately reach the boundary  ℓ , as expected for a finite domain with the boundary at  ℓ as the only exit.We may consider the transition to an infinite domain by taking the limit of large   : where again, as in the 1D case, the √ -term generates the long time tail.Note that the limit  → 0 in the infinite domain is less than unity, indicating the transience of diffusion in three dimensions.

Inner region
For the inner region, we have the boundary conditions: We find for the fluxes onto the boundaries: The long time limits give us the splitting probability: and indeed the sum of both is 1 as it should.  → ∞ gives us the same cumulative FP probability for an infinite system as in the case for a reflecting upper boundary [cf.Eq. ( 43)], lim

Scaling form
Again, according to Fig. 1, we let the partial trajectory start at a small distance  to its end point , thus  0 =  + ,  ℓ =  for φout and  0 = ,   =  +  for φin .The inner radius for the  in be  and the outer radius of the  out be .With the scaled Laplace variable  :=  2 / we have The limit for  → 0 is again the same as in the one dimensional case, Eqs. ( 29) and (30).
Figure 7 shows the partial FPT densities for the inner and outer domains (left panels).Similarly to the 1D case we find again the peak at small times, which shifts t/τ τ ϕ in (t)  t/τ τ ϕout(t) to the left and gets the higher (and narrower) the closer to the target the particle starts.Intermediate times are again governed by a  −3/2 power law decay.At large times, the cutoff sets in, again relatively sharply at  ≈  2 / for the inner domain with its absorbing outer boundary at |x| = , and more blurred out near  ≈  2 / for the outer domain with its reflecting boundary at |x| = .In the 3D case, the small shoulder at large times is more pronounced as compared to the 1D case.This is a well-known effect and is due to the compact exploration of space in low dimensions, where the first passage is more dominated by the direct trajectories [14].A smaller large-time shoulder indicates that the geometry plays a minor role in low dimensions than in 3D or higher.The right panels of Fig. 7 demonstrate again data collapse of the partial FPT densities if scaled with respect to the distance  of the starting point to the target.

Full solution for the FPT density
Again we use the partial FPT densities to obtain the Laplace transform of the full FPT densities for a particle starting at  0 and being absorbed at  by substituting into Eq.( 12).The coincidence of our solution for the FPT density for the homogeneous system obtained via GME approach with the time density for first passage to  starting from  0 , again underlines its vanishing dependence on  and , as expected.For an infinite outer domain  → ∞ we have: which becomes in time domain: The numerical Laplace inversion of Eq. (54) yields the final FTP densities  FPT (), (Fig. 8).Here, the difference to the 1D case (Fig. 6) becomes obvious: the plateau region is much more pronounced in the 3 dimensional case.Notice that the limiting FPT density for infinite domains Eq. ( 56) does not normalize due to its prefactor, while the FPT densities in confined domains do.The plateau in the FPT profile for reflecting boundaries does not only compensate for those particles that would have contributed to the power law tail in the case of an unbounded domain, but in addition it compensates also for the particles that would have got lost forever due to the transient character of diffusion in dimensions higher than 2.

Application: domains with different diffusivity
As an application to diffusion in a non-uniform, piecewise homogeneous medium made of two concentric spherical shells, we modify the above calculations for the 3D case and assign distinct diffusion constants,  in and  out , to the inner and outer domains, respectively; such a geometry was studied in Ref. [15].The modified diffusion constants enter merely the partial FPT densities, so that the same derivations as before go through.In particular, pFPT () obeys Eq. ( 12), but the explicit result will be different from Eq. (54).
In this setup it is also interesting to consider the case where the particle starts in the inner region,  <  0 < .As [15] point out, here the FPT density is governed by a third timescale which manifests in an additional intermediate regime in the FPT density profile.When the particle starts on the inside, the governing equations change slightly due to the different initial condition,  in (0) = 1,  out (0) = 0. Eqs.(10),(2a, 2b, 2c) remain the same, but for the recursive equations for the loss fluxes we get: where  0,in ,  0,∅ () is the dwell time probability density for a particle on an inner partial trajectory starting from  0 and ending at ( + ) or , respectively.In Laplace domain we get the linear system consisting of j− in = φin () for the loss fluxes and for the probability densities of a particle to occupy a state (i.e.be inside of radius , outside or stopped), from which we may extract the FPT density in Laplace domain (full solution see appendix): pFPT () = j− in,∅ = φ0,∅ () + φ0,in () φin,∅ () φout () For both cases of the particle starting from outside or within the radius , the  as calculated in 4 were inserted into Eq.( 12) or Eq. ( 60) with the respective diffusion constants.Finally the resultant expressions were Laplace inverted numerically.
Here, in contrast to the homogeneous case the -and -dependence of the FPT prevails, but the limit  → 0 always exists.We can see from Figs. 9 and 10 that already for moderate  there is hardly an influence on the final results regarding the FPT densities.Obviously the overall time spent in the region (,  + ) is negligible compared to the time spent in the inner and outer region so that the accumulated relative error made by constructing the auxiliary -shell around  remains small.As a comparison and for a qualitative discussion, we have also drawn the respective FPT density for an infinite domain into the FPT plots: Eq. ( 56) with  =  out for the particle starting outside (Fig. 9) and  =  in for the particle starting inside of  (Fig. 10).
In the case of the particles starting outside (Fig. 9), we see either an enhancement or a delay in the left peak indicating the particles that proceed directly to the target , depending on whether the diffusion in the inner region is fast (left panel) or slow     (right panel).Moreover, we have a broader tail in the case where particles are slower in the outer region (left panel).It takes a longer time for the particle to traverse the whole domain before the exponential cutoff due to finiteness of the system comes into effect.In the case of the particles starting from within the radius , we can see from Eq. (60) that the profile is a superposition of the FPT distribution of the direct trajectories φ0,∅ () and another term.For small times, φ0,∅ () is well approximated by the FPT density on an infinite domain.With respect to the broadening of the FPT density for small  out we have the same effect as for particles starting outside, but in addition, an intermediate regime emerges for those particles that initially leave the inner region but then proceed to the target without much exploring the outer region.A full quantitative analysis of all these regimes using a completely different approach can be found in [15].Thus, our approach provides the expected results for the FPT densities regarding the situation in radially symmetric domains with distinct diffusion constants and a central target.

Summary and conclusions
We have introduced a novel framework to stochastic first-passage problems in nonuniform, piecewise homogeneous environments, which is suitable to yield FPT distributions, but also splitting probabilities and total sojourn times in a domain.The central requirement is the Markov property of the underlying transport process, which allowed us to cast the FPT problem into a renewal problem.To this end, we coarsegrained the diffusion trajectories by dissecting them at the crossing points between adjacent domains, for example, between regions of different diffusivities; each partial trajectory is interpreted as a state, labeled by  (Fig. 1).The dwell times on each domain are determined by the dynamics in the domain, e.g., by the time it takes a molecule to leave the domain.The dwell time distributions   () parameterize the coarse-grained model and summarize the geometry, dimensionality, diffusion properties, etc. of the respective domain.They can be obtained from solving a first-passage problem on the (homogeneous) domain, or from simulations or even experiments.In general, the dwell times are not exponentially distributed with the consequences that the coarse-grained stochastic process is not a Markov jump process and that the evolution of the occupation probabilities   () does not obey a classical master equation.The latter is replaced by a generalized master equation (GME), which is a set of linear integro-differential equations for the probabilities   () and their fluxes, thereby accounting for memory effects.
A subtle difficulty of the approach is the singular character of first-return times to the same domain boundary in a continuous space.We showed that this issue can be solved by assigning a finite width  to the domain boundaries where needed (Fig. 3), which regularizes the partial FPT densities; the limit  → 0 is then taken in the final results.
As a test case, we exemplified the GME approach for two-domain models in 1D and 3D with domains of equal diffusivities and showed that the known FPT distributions  FPT () are recovered; here, the dependence on  dropped out.An analytical solution of the GME is readily obtained in frequency domain by a Laplace transform [Eq.( 12)] with explicit results for the overall FPT density [Eqs.(31) and ( 54)].We have complemented these results by a numerical backtransform to the temporal domain using an algorithm [28] that has proven robust for broad FPT densities extending over many decades (Figs. 6 and 8).Based on these results, it was straightforward to address a 3D target search problem with two domains of distinct diffusivities.The obtained FPT densities (Figs. 9 and 10) exhibit a complex structure, characterized by three time scales, and resemble the findings of Ref. [15].
Having the validity of the GME approach corroborated, we are now in a position to investigate first-passage problems in other, more complex scenarios.The modular approach of the method makes it particularly simple to extend the discussed twodomain models to heterogeneous spaces formed by a larger number of domains, or to assemble models for the interplay of different modes of transport, giving scope for a variety of problems of heterogeneous diffusion in space and also in time.On the other hand, it enables us to trace back the origin of certain features that emerge, e.g., in the FPT densities, since the results in frequency domain always algebraic compositions of the FPT densities of the partial trajectories, reflecting the properties of the respective spatial domain (or transport mode).
Manifestations of domain-specific transport behavior that can easily be accounted for in our framework on the level of the partial FPT densities are, e.g., the degradation (or even reaction) of molecules, changing the dimensionality of the space of motion, and different types of trapping subdiffusion.One could also implement crowding effects by using models of anomalous diffusion in certain domains, albeit this often introduces memory on the trajectory level that would be lost when molecules leave a domain.However, such a loss of memory may be justified if the transition between domains is associated with barrier crossing, e.g., due to a cellular membrane.Barrier crossing and directed channeling can be modeled in our framework by altering the flux balance equations and by, e.g., modifying the splitting of probability at the domain interfaces or by adding a bias.We leave these more complex models for future research.In brief, we believe that the presented GME-based framework is a potentially powerful toolbox to address first-passage related questions of heterogeneous diffusion processes for a variety of transport modes.

Figure 1
Figure1.Geometries of the first-passage problem with two domains and sketches of exemplary partial trajectories in a one-dimensional half axis (left) and between concentric spheres in three dimensions (right); the larger sphere with radius  > |r0| is not shown.Space is partitioned at the position/radius  into an inner (red) and an outer (blue) domain.The symbols  denote the dwell time densities of the different types of partial trajectories.

Figure 2 .
Figure 2. Graph of the simplified two-domain model (without an auxiliary shell).Arrows indicate the transitions between the states and are annotated by the corresponding loss fluxes  − from each state.

Figure 3 .
Figure 3. Geometries and exemplary partial trajectories of the first-passage problem with two domains as in Fig. 1, amended by an auxiliary boundary at |x| =  + .Note the hybrid character of the region  < |x| <  + : partial trajectories passing this region belong to either the in or the out state, depending on the domain where they start.

Figure 4 .
Figure 4. Graph of the two-domain model with an auxiliary shell.Arrows indicate the transitions between the states and are annotated by the corresponding loss fluxes  − from each state.

Figure 5 .
Figure 5. Numerical Laplace backtransforms (symbols) of the partial FPT densities in the inner domain [in(), Eq. (27), upper panels] and the outer domain [out(), Eq. (28), lower panels] for the 1D problem.Left panels show results for the partial FPT densities in,out(), and right panels the scaled partial FPT densities.The parameters of the geometry are 0 = 10,  = 5, and  = 20, and  =  2 / is the unit of time.The black dotted line denotes the analytical solution in the limit  → 0 [Eq.(30)].The gray dashed line indicates the asymptotic power law tail  −3/2 . .

Figure 6 .
Figure 6.FPT density as obtained by numerical Laplace inversion of Eq. (31) for different domain sizes  (in units of R), shown as symbols in different colors.Further parameters are 0/ = 10 and / = 5, and  =  2 / is the unit of time.The black dotted line denotes the analytical solution in the limit  → ∞ [Eq.(33)]

Figure 7 .
Figure 7. Numerical Laplace backtransforms (symbols) of the partial FPT densities in the inner shell [in(), Eq. (52), upper panels] and the outer shell [out(), Eq. (53), lower panels] for the radial problem in 3D.Left panels show results for the partial FPT densities in,out(), and right panels the scaled partial FPT densities.The parameters of the geometry are 0 = 10,  = 5, and  = 20, and  =  2 / is the unit of time.The black dotted line denotes the analytical solution in the limit  → 0 [Eq.(30)].The gray dashed line indicates the asymptotic power law tail  −3/2 .

Figure 8 .
Figure 8. FPT densities as obtained by numerical Laplace inversion of Eq. (54) for different domain sizes  (in units of ), shown as symbols of different color.Further parameters are 0/ = 10, / = 5, and  =  2 / is the unit of time.The black dotted line denotes the analytical solution in the limit  → ∞, Eq. (56).