Target finding in fibrous biological environments

We study first-passage time distributions of target finding events through complex environments with elongated fibers distributed with different anisotropies and volume occupation fractions. For isotropic systems and for low densities of aligned fibers, the three-dimensional search is a Poisson process with the first-passage times exponentially distributed with the most probable finding time at zero. At high enough densities of aligned fibers, elongated channels emerge, reducing the dynamics dimensionality to one dimension. We show how the shape and size of the channels modify the behavior of the first-passage time distribution and its short, intermediate, and long time scales. We develop an exactly solvable model for synthetic rectangular channels, which captures the effects of the tortuous local structure of the channels that naturally emerge in our system. Our results shed light on the molecular transport of biomolecules between biological cells in complex fibrous environments.


I. INTRODUCTION
Many biochemical reactions between chemically active molecules involve molecules distant in space, and commonly, at least one molecular species is free and searches for its target molecule. Thus, biochemical reactions depend on how a molecule diffuses toward its target, and the probability associated with the molecules to react once they are in close contact [1][2][3]. The former process depends on the reactant's diffusion coefficient D and the shape and size of the confining region. The latter process depends on an intrinsic reactivity k. Under ideal conditions, reactants in high concentrations are uniformly distributed in space, leading to uniform and independent encounters between molecules. Consequently, diffusion and kinetic controls are each correctly described by a single time-scale, and in particular, the mean reaction time is given as the sum of two time-scales: the mean time for molecular encounter ∼ 1/D, and the mean time for chemical reaction ∼ 1/k [4][5][6]. In recent years it has become more evident that it is necessary to question this simplified description of biochemical reactions and consider more elaborate models. In particular, for many biochemical reactions, the number of reactants can be low, limited to a few copies [7]. For example, gene expression and gene regulation occur at low-copy protein numbers, and their stochastic behavior has been the focus of many studies [8][9][10][11]. Another example is the sensory systems used by swimming bacteria responding to the activationdeactivation of membrane receptors by a limited amount of molecules [12][13][14]. In these cases, it is no longer appropriate to describe reaction rates with the mean time for a molecular encounter or the mean first-passage time (MFPT), but one needs to know the whole distribution of first-passage times (FPTs). It becomes then clear that at low molecular concentrations, the MFPT is not the only relevant time-scale of the reaction process, but the mostprobable FPT (MPFPT) becomes essential too. Hence, different diffusive-controlled events in the same system can vary widely in their time scales since the MPFPT, and the MFPT can differ by orders of magnitude [15,16].
We have recently shown that increased density and alignment of fibers facilitate molecular transport from a source to a target, which may support long-range cell-cell biochemical interactions [37]. In our model, we consider random walks of a diffusing molecule searching for its target within a system with fixed elongated fibers distributed with a nematic order parameter S and taking up a volume fraction φ. Only excluded volume interactions are considered. As the density of fibers increases, the system percolates differently depending on the alignment: for isotropic fibers (S = 0), the system undergoes a drilling percolation transition at φ 3D C = 0.75 [45][46][47]. At fiber occupation fractions of φ ≥ φ 3D C , the molecule gets caged by fibers, impeding target finding. Whereas for aligned fibers (S = 1), the system follows a 2D random site percolation process in the cross section of fiber positioning, with a critical fiber density of φ 2D C = 0.408 [48]. As φ reaches φ 2D C , the components of the diffusion coefficient perpendicular to fiber alignment decay to zero while the parallel component remains unaffected. This effect on the diffusion coefficient results from the emergence of channel-like structures that confine the dynamics to a 1D process, which is more effective than the 3D case. This caged state of the dynamics can be modulated toward the 1D process by continuously increasing fiber alignment S and the fiber volume fraction φ.
In this paper, we study the effect of channel shape and size on the FPT distribution, the MFPT, and the MPFPT of a target-finding process by numerical simulations and by analytically solving for several simplified geometries of channels with square and rectangular cross sections. We show that the channel size as well as its fractal shape influence the FPT distribution. Additionally, we show that the FPT distribution is no longer characterized by a single time-scale, implying that the typical notion of describing molecular reactions as the sum of two MFPTs (1/D and 1/k) is not appropriate.
The paper is organized as follows: in Sec. II, we introduce our model. In Sec. III, we study the FPT problem in the free case (φ = 0) and in the low fiber volume occupation fraction regime, both for aligned (S=1) and for isotropic (S=0) fiber distributions. Next, in Sec. IV we consider the very high volume occupation fraction limit (φ ≈ 1) of aligned fibers (S = 1), such that the dynamics are entirely 1D. Then, in Sec. V we consider high densities (φ > φ 2D C ) of aligned (S = 1) fibers and study the FPT in elongated channels with different shapes and sizes; We numerically analyze the FPT distributions that occur in the disordered channels that naturally emerge in our system of aligned fibers. Additionally, we approximate these complex channels with synthetic square and rectangular channels that allow us to study the FPT distribution analytically. In Sec. VI, we study the FPT problem for intermediate alignment 0 < S < 1, and map these more general cases to our results for perfectly aligned fibers (S = 1), both above and below the percolation threshold. Finally, Sec. VII concludes with a summary and discussion of our work and its implications.

II. COMPUTATIONAL MODEL
We study transport in complex environments using a model of particles moving on a 3D simple cubic lattice with periodic boundary conditions in all directions. The lattice is set with a total volume of V = L x × L y × L z sites, and three types of molecules can occupy the different lattice sites: a tracer molecule that is released from a source located at r 0 = (x 0 , y 0 , z 0 ), with x 0 = L x /2, z 0 = L z /2, y 0 = L y /4, a static target placed at r T = (x T , y T , z T ), with x T = L x /2, z T = L z /2, y T = 3L y /4, and elongated fibers, each one running along one of the three principal directions of the lattice and that span the whole system length. Fibers have a thickness of one lattice site and thus occupy a total volume fraction of with M the total number of fibers and p i the probabilities of a fiber to run along the x, y and z axes. We calibrate the probabilities p i to obtain a desired nematic order parameter S = (3 cos 2 θ − 1)/2, with θ the angle between the fiber orientation and the preferred direction of orientation [49], which we arbitrarily set to be the y-axis. We set p x to be equal to p z , thus S = (3p y − 1)/2, and p x = p z = (1 − p y )/2. A system with S = 0 has fibers isotropically distributed, with p x = p y = p z = 1/3, and one with S = 1 has all the fibers aligned along the y-axis, i.e., p x = p z = 0 and p y = 1. The FPT is defined as the time needed for the tracer molecule to reach the target site for the first time. After every finding event, the tracer is placed back to the source location, and then, a new random configuration of fibers is generated, and a new target search process begins.
At each unit of time, the molecule advances to one of its six neighboring sites, provided the desired location is empty of fibers. Otherwise, the move is rejected, and the time increases by a unit of time without the position of the molecule changing. For the case φ = 0, with no fibers, the molecule moves in an empty lattice and diffuses with a diffusion coefficient D 0 = 1/6. As the volume occupation fraction increases toward the percolation threshold φ C , the diffusion coefficient decays algebraically as predicted by percolation theory [48]. Hence, to describe the complex behavior of the diffusion coefficient one needs the critical density φ C and the exponent µ, which controls this algebraic decay [50]. In the isotropic (S = 0) case, diffusion is equally hindered in all directions and the diffusion coefficient decays following the Swiss-cheese model [50], as: Using the critical density for drilling percolation φ 3D C = 0.75, we obtain µ 0 = 2, agreeing with the value reported in [48], see Appendix A. For aligned fibers (S = 1), diffusion is not affected along the y axis, but is hindered along the x − z cross section. Thus, we decouple the effect of fibers in the components of the diffusion coefficient, and see that it decays as: with the critical density for 2D random site percolation φ 2D C = 0.408. We get that for our system µ 1 = 1.3, in agreement with the value reported in [48], see Appendix A.

III. 3D SEARCHING
We start by considering the reference case of target finding without fibers (φ = 0). This system with periodic boundary conditions allows us to solve the FPT distribution, based on the evolution equation for the probability, P (r, t), of the particle to be at position r = (x, y, z) at time t. The FPT distribution F (r T , t) of finding the target located at r T at time t, given that the molecule ( (1-f C 2D )V 2 3 9.8x10 5 9.2x10 5 FIG. 1. MFPT and FPT distributions for target finding for different values of S and φ. A) MFPT as a function of the available volume for S = 0 and S = 1. In the free case φ = 0, the MFPT increases with the volume. For low φ, the MFPT is given by Eq. (8) (lime line). Inset: The target finding dynamics strongly changes as φ approaches φ 2D C (blue line) and φ 3D C (orange line), for S = 1 and S = 0, respectively. B) Normalized FPT distribution for S = 0 and S = 1 for values of φ ≤ φ 3D C and φ ≤ φ 2D C , respectively. Black dashed line shows exponential behavior. C) Normalized FPT distributions for aligned fibers and values of φ > φ 2D C . The red line is Eq. (11) for the 1D FPT distribution. Statistics are performed for 2 × 10 5 finding events.
started at r 0 at t = 0, is related to P (r, t) by: The right hand side of Eq. (3) is the probability of reaching r T for the first time, given that at any previous time t the molecule was at r 0 . After taking the Laplace transform of both sides, we get the relation: By separation of variables, P (r, t) = P y (y, t) P x (x, t) P z (z, t), we obtain the independent probabilities P j (j, t), and find that the MFPT scales with the system's volume, see Appendix B: where is a geometrical prefactor, and Ω 0 (x, y, z) = ω (πx) + ω (πz) + ω (2πy) with the function ω(ψ) = (1 − cos ψ)/3. Here, the periodic boundary conditions ensure that over long times, the tracer molecule is equally likely to be at any lattice site in the system, making the finding events a Poisson process, with an exponential FPT distribution F (t) = 1/ t × exp (−t/ t ) [51]. To test the theoretical prediction of Eq. (5), we plot in Fig. 1A the MFPT for two system sizes with L x = L y = L z = 100 (black dashed line) and L x = L y = L z = 97 (red diamond).
We recover the MFPT dependence on the system's volume and the free diffusion coefficient D 0 , as shown by the green dashed line with φ = 0. Next, we examine the case of low values of fiber occupation fraction φ and their effect on the MFPT. The introduction of fibers has two competing effects on the MFPT. On the one hand, as φ increases, the available volume decreases, reducing the MFPT for target finding. On the other hand, the presence of fibers hinders diffusion, thus increasing the time needed for the tracer to find its target. In our model, we capture the MFPT behavior by using a mean-field approximation for the case φ 1. Specifically, for aligned fibers (S = 1), diffusion is hindered only along the x − z plane, and we approximate by 1 − φ the probability of succeeding to move in this plane. For other values of S, the diffusion is also hindered in the y axis. Hence, we approximate by 1 − φ (p x + p y ) = 1 − φ (p z + p y ) = 1 − φ(2 + S)/3 the probability of succeeding to move in the x − z plane, and by 1 − φ (p x + p z ) = 1 − 2φ(1 − S)/3 the probability of succeeding to move in the y axis. Therefore, the function Ω 0 (x, y, z) used in Eq. (6) for φ = 0 depends on φ and on S via the new function Ω (x, y, z, φ, S) = ]ω (2πy). We suggest that α 0 in Eq. (5) should be replaced by α(φ) which is obtained by substituting this expression for Ω 0 in Eq. (6). After expanding α(φ) to first order in φ, we obtain that α(φ) = α 0 − α 1 φ, with Note that by definition p x + p y + p z = 1, and thus to first order in φ, the MFPT is independent of S. Thus we expect that for low φ the MFPT will be given by To test our predictions, we plot in Fig. 1A the MFPT as a function of (1 − 2φ/3)V for values of fiber volume occupation φ < 0.15. We see that for these low values of φ, fiber alignment S does not affect target-finding (circle and purple data). We observe agreement between our numerical results and Eq. (8).
To characterize the target finding dynamics, we obtain the FPT distributions for values of fiber occupation fractions below the percolation thresholds. Both for S = 0 and for S = 1, the FPT distributions follow an exponential behavior for low φ, while for fiber densities approaching φ C , the FPT distributions exhibit deviations from an exponential distribution, indicating that as the system approaches percolation the FPT distributions are not fully characterized by a single time scale, see Fig. 1B.
Target finding dynamics for fiber densities higher than the critical thresholds are very different if fibers are aligned or isotropically distributed. In the case S = 0, the diffusion coefficient decays to zero for φ > φ 3D C , and thus, the MFPT diverges, and the FPT distribution is no longer defined, as shown in the inset of Fig 1A. For S = 1, the MFPT follows a complex behavior for φ ≥ φ 2D C , exhibiting a sharp decrease around φ 2D C , as shown by the inset in Fig. 1A. Moreover, the FPT distribution for φ ≥ φ 2D C follows a non-monotonic behavior characterized by three different time scales: the MPFPT at short FPTs, the MFPT at intermediate time scales, and the time scale of distribution tail (TSDT), at long FPT, as shown in Fig. 1C. We define the TSDT as the decay rate of the exponential tail of the FPT distributiont. For these cases, channels in the x−z cross section are formed, reducing the dimensionality of the dynamics from 3D to 1D. In this paper, we focus on how the channel structure affects the FPT distribution, and we study in detail this regime in Sec. V below, albeit we consider first the simpler full 1D limit, which is obtained for φ ≈ 1.

IV. 1D SEARCHING
In this section we present the liming case of a high density (φ ≈ 1) of aligned fibers (S = 1), for which the diffusing molecule is confined to a 1D line along the direction of the fibers. In our lattice model, the channel is aligned along the y axis and has a cross section equal to one. Due to the periodic boundaries of the system, the topology of the channel can be understood as a ring-like structure with a circumference of L y and a single target that can be reached by the tracer molecule either from the left or from the right side of the ring. The distribution of FPT to the target is related to the survival probability, H(t), that the tracer did not yet reach the target up to time t, by [52]: Due to the periodic boundary conditions, the survival probability is equal to the probability that a tracer diffusing on a finite system of length L y with absorbing boundary conditions, remains in the system: We obtain the probability of finding the tracer at a position y at a time t, P (y, t), by solving the diffusion equation. In Appendix C we show that for the specific case that the initial position of the tracer is equidistant from the two boundaries, i.e. y 0 = L y /2, the FPT distribution is given by: where ω(ψ) is defined above and k n = πn/L y . This FPT distribution has multiple time scales. At long times, the FPT distribution decays exponentially, F (t) ∼ exp −t/t , with the TSDT given by: The MFPT of the distribution given in Eq. (11) is: We also obtain from Eq. (11) the MPFPT of the distribution: with β ≈ 0.411 being the positive solution of the transcendental equation These characteristic time scales all scale with the system size and the diffusion coefficient as L 2 y /D 0 , but exhibit different prefactors: 1/π 2 ≈ 0.1 for the TSDT, 1/8 for the MFPT, and β/(π 2 ) ≈ 0.04 for the MPFPT. The first two time scales are similar, with the TSDT slightly smaller than the MFPT, while the MPFPT is about one order of magnitude smaller than the TSDT and the MFPT. Figure 1C shows the perfect agreement between the analytical expression for the 1D FPT distribution from Eq. (11) with numerical simulations of 1D channels with L y = 100. At long times, the behavior of the FPT distribution is fully described by a simple exponential decay. Importantly, when considering our percolation system, we see that for aligned fibers (S = 1) at high occupation fractions, the FPT distribution approaches the 1D FPT distribution. This results in a reduction of dimensionality in the dynamics of the system and the emergence of narrow, elongated channels, and will be the focus of the following section.  [48]. B) Distribution of the radius of gyration of 4 × 10 5 channels with a fixed value of N = 484. Lower inset: The radius of gyration increases as N 1/d f . Upper inset: The variance of N increases quadratically as a function of the mean channel size N . Statistics are performed for 5 × 10 4 randomly generated fiber configurations for each volume fraction.

V. QUASI-1D SEARCHING
For fiber occupation fractions φ 2D C < φ < 1 and S = 1, elongated channels with complex cross-sectional shape emerge, and as φ approaches φ 2D C from above, the channel structure becomes fractal. We thus focus on how channels change as a function of the volume occupation fraction and plot in Fig. 2A the distribution of the number of lattice sites F N (N ) in the x−z cross-sectional area of the channels for different values of φ. For values of φ just above φ 2D C , the channel size is distributed with an exponential tail but with a clear shoulder. Hence, channels manifest two characteristic size scales, one for large channels and another one for more compact channels. As φ increases, the distributions shift toward smaller values of N ; the distribution becomes narrower, and singleexponentially distributed [53]. In the inset of Fig. 2A, we plot the average channel size in the cross-sectional area N as a function of φ − φ C , together with the known relation from percolation theory [48]. Specifically, we show in Fig. 1C that systems with fiber occupation fractions of φ ≥ 0.6 follow the FPT distribution of the 1D target-finding process, indicating that for narrow channels with values of N ≤ 10, the dynamics are effectively 1D.
To further understand the structure of the channels, we calculate the channel radius of gyration R g , which is defined from: with x i and z i , being the position of a lattice site within the channel's cross section. We now fix a cross-sectional area of the channel to N = 484 sites and get the distribution of the radius of gyration F Rg (R g ), as shown in Additionally, the lower inset in Fig. 2B shows that the radius of gyration scales as the fractal dimension of 2D random percolation [48,54]. Finally, in the upper inset of Fig. 2B we plot the variance of N as a function of N , and observe that it increases following the relation Var(N ) ∼ N 2 .

A. FPT distribution in disordered channels
In Fig. 1C we show that for channels with L y = 100, S = 1, and φ φ 2D C the FPT distribution is not the 3D exponential distribution, nor the 1D distribution. Instead, for these cases, as φ approaches φ 2D C from above and the channels become fractal, the MPFPT becomes more pronounced and shifts to smaller values. Note that for φ = 0.42, the MFPT and the MPFPT differ by almost two orders of magnitude.
An important point to consider is how the channel size modulates the FPT distribution. Thus, we next choose a fiber occupation fraction of φ = 0.45 and randomly select nine channels with different values of N and run our random-walk simulations for each channel configuration, as shown in Fig. 3A. Here, it becomes evident that the channel's cross-sectional area determines the shape of the FPT distribution. For example, for N = 5, the distribution follows the 1D behavior of the situations with high volume occupation fractions. For N = 963, the distribution has a pronounced maximum at the MPFPT at low FPTs, similar to the case φ = 0.42, in which the fiber density is close to the percolation threshold. Simulations show that as the cross-sectional area of the channel increases, the MPFPT monotonically decreases in values F(t In the case N = 26, the channel is narrow, and the FPT distribution is similar to that of the 1D case. The shape of the channel with N = 443 is more complex, exhibiting internal holes and sharp edges, that support relatively directed trajectories of the tracer toward its target, leading to an FPT distribution with a pronounced MPFPT at short time scales.
We now focus on the effect that channel shape has on the FPT distribution in the natural channels obtained from randomly positioning fibers. For this, we fix the cross-sectional area at N = 484, and choose channels with three different shapes, with R g = 12, 14.9, and 19.4. Figure 3B shows that, intriguingly, the FPT distributions for these three channels seem qualitatively similar, despite the difference in their values of R g . This behavior contrasts with the previously observed effects of φ and N . The FPT distributions exhibit a pronounced MPFPT similar to the systems with φ φ 2D C . Quantitatively, the channel with R g = 19.4 displays a more pronounced MPFPT than the channel with R g = 12.
These findings suggest that for the natural channels, the radius of gyration moderately modulates the FPT distribution. For completeness, we plot in black in Fig. 3B the FPT distribution obtained from the whole ensemble of channels with N = 484.

B. Synthetic channels
To further understand the effect of channel shape and size on the FPT distribution, we now consider synthetic channels of predefined shapes and sizes, with square cross section N = L x × L x (Fig. 4A), as well as rectangular channels with different aspect ratios N = L x × L z (Fig.  4B). Note that the lengths L i are the lattice dimensions in all our simulations, whereas L i are the dimensions of the synthetic channels considered here. Thus, in general L i ≤ L i . Specifically, for both channel shapes we choose L y = L y = 100. The radius of gyration for these systems is These channel geometries allow us to solve the FPT distribution, based on the evolution equation for the probability P (r, t) of the particle to be at position r = (x, y, z) at time t. For that, we make use of Eqs. (3) and (4) above, and express P (r, t) = P y (y, t) P x (x, t) P z (z, t).
We separately obtain P x (x, t) and P z (z, t), by implementing reflecting boundary conditions at x = 0, L x and z = 0, L z , and P y (y, t) by taking periodic boundary conditions, i.e., P y (y, t) = P y (y + L y , t). Additionally, for simplicity, we assume that L y is even. After taking the Laplace transform of P (r, t) and using Eq. (4), we ob-tain the Laplace transform of the FPT distribution, or its generating functioñ with and the function ω(ψ) defined in Sec. III above. For the complete derivation, see Appendix B. Similar expressions for the generating function of the FPT distribution were derived in [55] for a d-dimensional system with arbitrary boundary conditions. Here, we concentrate on thoroughly investigating the specific system at hand of rhombic 3D channels.
From the Laplace transform of the FPT distribution, we obtain the MFPT t = nx,ny,nz We note that the sum over n x , n y , n z in Eq. (20) includes all values of n x , n y and n z between 1 and 2L x , L y and 2L z respectively, except for the single point n x = 2L x , n y = L y and n z = 2L z , such that Ω 0 is always positive.
In general, inverting the Laplace transform of the FPT distribution is not trivial. Therefore, we approximate the Laplace transform of the FPT,F , by functions,F m , which agree in the first m terms in their Taylor expansion. For m = 1, 2 these approximations can be inverted explicitly by Higher order approximations are given by where ξ k are the roots of an m'th order polynomial. The polynomial and the coefficients A k are further detailed in Appendix B. We note that the second order approximation, The range of validity of higher order approximations is smaller, as discussed in Appendix B. In particular, we find that if both L x and L z are smaller than L y , then the second order approximation is valid. Figure 4A shows the FPT distributions of our simulations for square channels with different values of N . Additionally, we plot in dashed lines in Fig. 4A our second-order approximation F 2 (t) of the FPT distribution, showing excellent agreement with our simulations. Interestingly, for all considered cases, the shape of the FPT distributions differ from the ones that are naturally obtained from percolation. For small squares, the distribution follows the 1D behavior, as shown above for narrow channels. The MPFPT shifts toward lower FPTs when increasing the square lateral size, but without sharply increasing its peak, in contrast to the natural percolation case. Note that the shape of the distribution gradually changes from the 1D limit to the exponential behavior seen for the cubic case with φ = 0 and L x = L y = L z = 100. We conclude that channel size strongly modulates the magnitude of the MPFPT, in some cases making it more than one order of magnitude smaller than the MFPT. However changing channel size using the simplest square-shaped synthetic channels is not enough in order to capture the qualitative evolution of the FPT distribution seen above for natural channels.
In order to better see the location of the MPFPT and the behavior of the FPT at small times, we compare the FPT distribution from the numerical results to the analytical approximations up to third order. For narrow channels, the second and third order approximations give progressively better results for short times, see Fig. 5. For wider channels, the first order approximation (a simple exponent) agrees very well with the results at times 3 x 3 x 100 longer than the MPFPT, as shown in Fig. 5C and D.
Note that the location of the peak at the MPFPT as predicted by the second order approximation agrees very well with the numerical results, and there is no appreciable improvement given by the third order approximation in this regard, even for the narrower channels shown in Fig. 5A and B.
Next, we fix the cross-sectional area at N = 484 and consider rectangular channels with different aspect ratios, see Figs. 4B and 6. The FPT in the rectangular case 44 × 11 is distributed similarly to the square shape 22 × 22. As the aspect ratio increases to 121 × 4, the FPT distribution exhibits a pronounced peak around the MPFPT, in clear contrast to the shape of the FPT distributions of square channels. Larger values of the aspect ratio further increase the peak of MPFPT and shift its value toward lower values of FPT/MFPT. For the fully elongated channel 484×1, with a width of one lattice site, the MPFPT is almost two orders of magnitude smaller than the MFPT, stressing the strong effect of the channel shape on the FPT distribution. Dashed lines are the second-order approximation F 2 (t) for the FPT distribution of the corresponding channel shape. We note that our solvable model does not entirely capture the complex shape of FPT distributions for elongated channels, as shown in Fig. 6. From our solvable model we can also obtain an approximation for the location of the MPFPT, by looking at the 44 x 11 x 100 From Figs. 4B and 6B we see that although the position of the MPFPT is accurately captured, the magnitude of the peak is not. Figure 7 shows the very good agreement between the approximation for the MPFPT, Eq. (25), and the simulation results. We find that for square cross sections the MPFPT has a maximum value and that it vanishes for L x close to the system's length L y , while the ratio between the MPFPT and the MFPT is a decreasing function of the cross-section size.

C. Channel structure modulates the characteristic time scales
We now examine the effect of channel size and shape on the three different time scales MFPT, MPFPT, and TSDT. First, we plot in Fig. 8A the MFPT for different channel geometries as a function of N . We observe that for small cross-sectional areas of N < 10 (10% of the channel's length L y = 100), the MFPT effectively follows a 1D dynamics. In this regime, channels with a square or elongated shape or obtained from our percolation model, exhibit the same MFPT. Here, the MFPT scales with the source-to-target distance squared, L 2 y . As N increases, the 3D shape of the channels starts to affect the dynamics of the tracer, and the MFPT rapidly increases in all considered channel shapes. Therefore, modulation of the ratio between N and L y controls the reduction of dimensionality in the target-finding dynamics. Remarkably, for values of N ≥ 10, channel shape influences the MFPT. Linear channels with a width of one lattice site, exhibit the highest MFPT. On the contrary, for a given value of N , square channels have the lowest MFPT. Note that for values of N > 200, square channels behave as the 3D system characterized by Eq. (5). Our solvable model captures the MFPT for these two extreme channel shapes, as shown by the blue and purple dashed lines. Interestingly, the ensemble of natural channels quantitatively behaves more similarly to the synthetic linear channels than to the square channels. We note that for this case, we use the average channel size obtained for each value of φ, instead of a fixed value of N . The MFPTs from these channels are located between the square and linear channels.
Motivated by the MFPT dependence on channel structure, we study how the MPFPT and the TSDT correlate with the MFPT. In general, we see that for all channel shapes, the TSDT is highly correlated with the MFPT, see Fig. 8B upper panel, and therefore, with the cross-sectional area N of the channel, see Appendix D. Moreover, we see that the TSDT-MFPT correlation is affected by the channel shape. Specifically, we see that for square channels, the MFPT converges toward the TSDT, increasing with N . Instead, for linear channels with N ≥ 10, the TSDT becomes larger than the MFPT and rapidly grows with N . Note the excellent agreement between our solvable model with our simulations. For channels originated from our percolation model, the TSDT follows a similar behavior of the linear channels.
The short-time behavior of the FPT distribution, which is characterized by the MPFPT, is also channel shape-dependent. We plot in the lower panel of finds its target in a relatively direct manner. We consider square channels with values of N ≤ 70 2 and see that the MPFPT increases with the MFPT. As the channel's cross-sectional area increases, the trajectories become less directed, and the MPFPT increases. Importantly, we previously showed in Fig. 7 that the MPFPT sharply decreases as the system dimensions L x and L z approach L y and the lattice becomes cubic. In this limit, the MPFPT corresponds to perfectly directed trajectories toward the target. Still, the probability of such events is very low, and the particle needs to scan, on average, the systems volume to find its target. For linear channels, the MPFPT starts increasing for low values the MFPT but then saturates and remains constant as the MFPT increases. Here, the narrow structure of the channel ensures that the directed trajectories toward the target are similar, despite the differences in the channel length. In case the tracer leaves the vicinity of the target and diffuses toward the edges of the channel, the MFPT increases, and in particular, the TSDT happens to dominate the long-time behavior for the FPT distribution. Similar to the behavior of the MFPT and TSDT, the MPFPT for channels from our percolation model is quantitatively similar to the linear channels, indicating that the local fractal structure of the channels confines the directed trajectories toward the target.
We now fix a value of N = 484 and quantify how the channel structure, characterized through R g , modulates the different time scales. We take channels with varying values of R g = 12, 14.9, and 19.4, from the natural percolation model and compare their time scales with the square and linear channels. We show in the upper panel of Fig. 8C that the MFPT increases with R g , i.e., channel elongation. As R g increases, the time needed for the tracer to come back to the vicinity of its target increases, spending time in regions of the channel where the target is not present, and thus increasing the MFPT. Remarkably, although all cases have the same channel volume V = N × L y , channel shape modulates the MFPT, implying that the MFPT no longer defines the diffusion-controlled process, affecting the typical notion of reaction dynamics. Similarly, the TSDT increases with the channel elongation, however the MPFPT is highest for square channels. Specifically, the MPFPT is similar for the three considered natural channels and the linear channel. These results support the notion of local channel compactness supporting the fast trajectories toward the target.

VI. INTERMEDIATE FIBER ALIGNMENT
Lastly, we consider an intermediate value of S = 0.5 of the fiber nematic order parameter and numerically obtain the FPT distributions for different values of φ in channels with L y = 100, as shown in Fig. 9. We note that all the systems with S = 1 undergo a drilling percolation transition at some S-dependent critical density φ 3D C (S). Specifically for the case S = 0.5, the critical density is φ 3D C ≈ 0.83 [37]. For values of φ > φ 3D C the FPT diverges. For φ = 0.3, the FPT distribution follows the expected exponential behavior. At high values of φ = 0.65, the FPT distribution qualitatively follows the shape of the distribution for the 1D case. For the intermediate case φ = 0.5, the FPT distribution exhibits a pronounced MPFPT at t/ t ≈ 0.025. Remarkably, for increasing fiber densities, the FPT distributions for the S = 0.5 system behave qualitatively similar to the case with S = 1, indicating that we can map the S = 0.5 system onto the S = 1 case. This mapping consists of identifying the fraction of fibers that run into the x − z plane in the case with S = 0.5, and then, considering a new system of aligned fibers S = 1 with such areal fiber density. Specifically, we find the areal densityφ = M/(L x L z ), by using p x = p z = (1 − S)/3 and p y = (2S + 1)/3 in the previously given relation for φ, and solving the third- order equation: The areal fiber density on the x − z plane is the real solution of Eq. (26) times p y , i.e.,φ(2S + 1)/3. Note that we consider a single value ofφ since L x = L y = L z . For example, for S = 0.5 and φ = 0.5, the corresponding mapping is a system with S = 1 and φ = 0.387. Figure 9 shows the impressive agreement of the FPT distributions for the considered mapped systems. Therefore, our model can be applied to relevant biological systems, which in general have nematic order parameters 0 < S < 1.

VII. DISCUSSION
We demonstrated how the FPT distribution of a target-finding process is affected by the size and shape of the domain in which the process takes place. In our lattice model, before the system reaches percolation, the FPT process is characterized by an exponential distribution, and the MFPT scales with the available volume. As φ approaches the critical densities of the system φ 3D C and φ 2D C , for S = 0 and S = 1, respectively, we observed significant deviations from exponential behavior. In the presence of isotropic fibers, once the system reaches percolation at φ 3D C , the FPT diverges. On the contrary, when fibers are fully aligned, the dynamics are richer, and FPT distributions with different shapes emerge. Thus, fiber alignment and fiber density are essential in the understanding of the target-finding process. Additionally, we showed that for intermediate values of fiber densities φ φ 2D C of aligned fibers, complex FPT distributions are obtained. These distributions are characterized by three time scales: the MPFPT, the MFPT, and the TSTD. We saw that by modulating the channel size and its shape, the time scales were strongly affected. For small cross-sectional areas of N ≤ 10, channel shape did not affect the FPT distribution much, and the characteristic time scales remained invariant. In contrast, for channels with larger cross-sectional areas, the FPT distribution changed with channel shape. Specifically, we saw that linear channels exhibited higher MFPTs than the compact square channel shapes. We also showed that the long-time behavior of the FPT distribution, characterized by the TSTD, is correlated with the MFPT and less sensitive to channel shape. Moreover, we numerically demonstrated that the short-time behavior of the FPT distribution is very sensitive to channel shape. We saw that for linear channels, the MPFPT effectively remains constant as the MFPT increases. These results indicate that the confined structure of the linear channels supports relatively directed trajectories of the tracer toward its target. Such directed trajectories remain unaffected as the channel aspect ratio increases since, in those finding events, the tracer does not escape from the vicinity of its target, and the edges of the channels are not explored. Contrarily, for square channels, the MPFPT increases with the MFPT for channels with cross-sectional areas of N ≤ 70. Here, as N increases along the x − z plane, the directed trajectories defocus, leading to larger MPFPT and MFPT. When considering the ensemble of natural channels from our percolation model, we observed that the FPT distributions and its characteristic time scales qualitatively behave more similarly to elongated channels. These observations indicate that the complex fractal shape of the channels support directed trajectories to the target, in a similar way as the linear channels do.
In general, biological systems have a nematic order parameter that lies between the extreme values of S = 0 and S = 1. Specifically, the analysis of collagen fiber alignment from ECM porcine urinary bladder obtained an average fiber alignment of S ≈ 0.54 [56]. Also, in artificial ECM environments, such as collagen or fibrin hydrogels, the band area between communicating cells reaches values of around S = 0.5 to 0.7 [37]. We, therefore, studied the intermediate cases with fiber alignments 0 < S < 1 and showed that they are accurately captured via a mapping onto the case S = 1. The mapping was achieved by obtaining the occupation fraction of the fibers that are aligned along the preferred direction, and then setting a new system with only those fibers, in such a way that the nematic order parameter is S = 1. By carefully obtaining the critical density of the system with 0 < S < 1, we showed that bellow percolation, those two systems exhibit very similar FPT distributions and char-acteristic time scales. Thus, the understanding gained for the aligned case (S = 1) is applicable for more realistic systems with intermediate values of the nematic order parameter S.
Our findings show that the classical description of diffusion control, given only in terms of the MFPT, is not accurate for biochemical reactions in environments with complex shapes and at a low number of molecules. Instead, the MPFPT and the TSDT should also be considered. We observed that the shape of narrow channels modifies the FPT distribution, modulating the magnitude of the MPFPT and affecting the MFPT. Interestingly, the MFPT is different for channels with the same volume but different shapes. Consequently, the law of mass action, which states that the reaction rates are directly proportional to each of the reactant concentrations [57], does not apply to our case with channels that support pronounced MPFPTs. Therefore, for sensory systems in cells that respond to low molecular concentration, two first-passage events will be characterized by very different reaction times [12,16]. A natural extension of our work is to consider the impact of attractive non-specific interactions between the tracer molecule and the elongated obstacles on the characteristic time scales of the target finding process [58].
Previous studies have shown that cells continuously remodel the ECM structure by applying forces to the fiber elements [25][26][27], and by degrading or generating new ECM fibers [28,29]. As this remodeling takes place, transport of molecules can be affected, leading to the possibility of biochemical-mechanical signaling feedback [37,38]. Our work provides a theoretical basis for such experiments, with a deeper understanding of how fiber remodeling impact molecular transport. Moreover, due to the complex structure of the ECM in bacterial biofilms, our work can provide further understanding of quorum sensing mechanisms and signal transduction in bacterial populations [59]. We plot in Fig. A1 the diffusion coefficient D(φ) as a function of fiber density φ. For isotropic distribution of fibers (S=0), D(φ) follows Eq. (1) with φ 3D C = 0.75 and µ 0 = 2, and for aligned fibers (S=1), D(φ) follows Eq.
(2) with φ 2D C = 0.408 and µ 1 = 1.3. In the inset we plot the algebraic scaling of the diffusion coefficient for S = 0 and S = 1. A1. Fiber density modulates the diffusion coefficient. The ratio D(φ)/D0 decreases as a function φ. Inset: loglog plot of D(φ)/D0 × (1 − φ) as a function of (1 − φ/φC ). For S = 1, we consider only diffusion along the x − z plane. Therefore we remove a value of 1/3, which accounts to the contribution in the diffusion coefficient along the y axis.

Appendix B: 3D synthetic channels
In this appendix we solve the discrete diffusion equation for a particle moving on the cubic lattice in a 3D channel and from it we infer the FPT distribution. In subsection B 1 we derive the solution for the diffusion equation. In subsection B 2 we derive the Laplace transform of the FPT distribution. The moments of the FPT distribution are derived in subsection B 3. The approximations of the FPT distribution are discussed in subsection B 4, and in subsection B 5 we derive the asymptotic expression for the MFPT.
Consider a particle moving on the cubic lattice in a 3D channel. In the x−z plane, the channel has a rectangular cross-section of size L x × L z , with reflecting boundary conditions, while in the y axis it is periodic with length L y . We are interested in the distribution of the FPT for the particle to first reach a specific target position r T = (x 0 , , z 0 ), which differs from its initial location, r 0 = (x 0 , 0, z 0 ), only in the y coordinate. In what follows we assume for simplicity that L y is even and that = L y /2, such that with the periodic boundary conditions, r 0 and r T are the farthest away possible along the y-axis. We derive the expression for general values of x 0 , z 0 , but later on concentrate on the specific case where x 0 and z 0 are in the middle of the channel.
The FPT distribution, F (r T , t) is related to the probability to find the particle at location r at time t, P (r, t) by The first term on the right hand side of Eq. (B1) is the probability that at time t = 0 the particle is already at location r T , and the second term is the probability that the last time before the particle reached site r T it was at site r 0 at time t , and then in the time interval t − t it reached site r T once. Taking the Laplace transform of both sides yields and therefore, since r T = r 0 , Furthermore, we can decompose the 3D motion of the particle into independent motions along the three axes, such that P (r, t) = P y (y, t) P x (x, t) P z (z, t) , where P j (j, t) is the probability for a 1D walker to be at j at time t, for j = x, y, z. In what follows we derive the 1D probabilities P j (j, t) and later transform the probability P (r, t) to Laplace spaceP (r, s) in order to obtaiñ F (r T , s).

Derivation of the 1D probabilities Pj (j, t)
The 1D probabilities P j (j, t) for j = x, y, z evolve according to the discrete-space diffusion equation, where we dropped the explicit dependence of P j on t for brevity. The prefactor of 3 on the left hand side comes from the 1D motion accounting for one third of the total 3D motion of the particle, which moves with rate τ −1 =1. On a coarse-grained level, the rate τ −1 is related to the diffusion coefficient D 0 by Assuming a solution of the form For j = y we impose periodic boundary conditions P y (y) = P y (y + L y ), while for j = x, z we impose reflecting boundary conditions These boundary conditions set restrictions on the allowed values of k Using the initial condition yields for j = y P y (r y , t) = 1 L y Ly t e 2iπny/Ly , and for j = x, z Setting j = j 0 in Eq. (B14) yields Note that ω(0) = 0, while for k > 0 ω(k) > 0.

FPT distribution in Laplace space
The Laplace transform of the FPT distribution, or its generating function is therefore Ly ny=1 Ly ny=1 with g (n, x, L) = cos 2 πn (2x − 1) 2L . (B18) We now assume that x 0 and z 0 are in the center of the cross-section, such that if L x (L z ) is even then . Under this assumption the functions g (n, x, L) for even and odd values of L are given by g even (n, L) = 1 2 1 + (−1) n cos πn L , We compare in Fig. B1 the analytical result for the Laplace transform of the FPT distribution,F (s), with the numerical results. We find that at small values of s the agreement is excellent, while at higher values of s where the agreement is lacking, the value ofF (s) itself is extremely small. The difference is due to limited statistics in the simulations. In order to see that this is indeed the reason, we ran more simulations and saw that the numerical results approach the analytical results as the number of realizations increases.

Moments
From the Laplace transform of the FPT distribution we can find all its moments by t n = (−1) n ∂ nF (r T , s) ∂s n s=0 . (B20) The first moment, or MFPT is t = nx,ny,nz where Ω 0 (x, y, z) = ω (πx) + ω (πz) + ω (2πy) , and the sum over n x , n y , n z includes all values of n x , n y and n z between 1 and 2L x , L y and 2L z respectively, except for the single point n x = 2L x , n y = L y and n z = 2L z , such that Ω 0 is always positive. This special point is excluded since when taking the derivatives ofF (s), there are contributions to the sums from both the nominator and denominator of Eq. (B17), and since in this special point Ω 0 = 0 the contributions cancel each other. The second moment is We note that sinceF (s) is regular around s = 0, all the moments exist.

Approximations for the FPT distribution
Since inverting the full form of the Laplace transform of the FPT distribution is not practical, we consider approximations that can be inverted. In order to find the m'th order approximation, we first expandF (r T , s) to m'th order in s In the next step we approximateF m (r T , s) bỹ such that when Eq. (B26) is expanded to m'th order in s we retrieve Eq. (B25). For k = 0, 1, 2, 3 we find that Inverting these approximations for the Laplace transform yields successive approximations for F (r T , t). In general, the approximated FPT distribution is given by where ξ j are the roots of the polynomial m k=0 a k ξ k . For m = 1, 2, it is possible to find an analytical expression for the roots ξ j , such that Eq. (B28) is explicitlỹ For m ≥ 3 we can find the coefficients a k exactly for any system size, and solve the resulting polynomial numerically.
Note that these approximations are valid only if the real part of all the roots ξ j is negative. For the first order approximation, m = 1, the only root is − t −1 which is negative. The second order approximation is valid only if For the higher order approximations, it is straightforward to check that the m'th order approximation is valid only if a k ≥ 0 for all k ≤ m. Therefore, the range of validity for each successive approximation is smaller than the previous one. Figure B2 shows the values of L x and L z for which the second order approximation is valid for different values of L y . First, we see that the validity depends on the ratios L x /L y and L z /L y . We observe numerically that the approximation is valid when where γ depends on L y and is very close to 1/3. In the case L x = L z , Eq. (B31) reduces to L x /L y < γ + 1/ √ 2. In order to evaluate γ for larger systems, we concentrate on the case L x = L z , and find the largest L x for which the approximation is valid. Figure B3 shows the ratio L x /L y between the largest L x for which the approximation is valid and the system's length L y . It appears to converge to a value slightly below 1.05, i.e. γ converges to a value slightly below 0. 35. In the region where it is valid, the second order approximation, F 2 (t), has a single maximum at the MPFPT t * = 2 t 2 − t 2 Increasing the size of the cross section, L x × L z , or decreasing the length of the channel, L y , decreases the value of t * . Note that although t * = 0 only in the case of an infinite system, its approximation, Eq. (B32), reaches 0 at the edge of validity for the second order approximation when 2 t 2 = t 2 .

Asymptotic expressions for large Lx, Ly, Lz
In this section we derive asymptotic expressions for t when L x , L y and L z are very large. Changing the sums over n x , n y and n z in Eq. (B21) to integrals over x = n x /L x , y = n y /L y and z = n z /L z , and approximating the function g(n, L) as 1/2 yields where we have used Eq. (B6) and In order to see how good this approximation is we use the Euler-Maclaurin formula for the p'th order approximation of a sum where B 2k are the Bernoulli numbers, f (k) (x) is the k'th derivative of f (x), and R p is an error term. Using Eq. (B35) for each of the three sums in Eq. (B21), we find that the error term is of order L −(2p+1) i for i = x, y, z. Next, note that the summand is such that the function f (x) in Eq. (B35) satisfies f (x) = f (L − x). Therefore, the second and third terms on the right hand side of Eq.
(B35) are identically zero. Hence, we conclude that the difference between the exact MFPT, Eq. (B21), and the approximation, Eq. (B34), is smaller than L −p x,y,z for all p. This means that the error is exponentially small in L x,y,z .

Appendix C: 1D searching
In this appendix, we obtain the FPT distribution of a 1D random walk problem by solving the discrete diffusion equation. Due to the periodic boundary conditions of our model, the topology of the channel can be understood as a ring-like structure with a circumference of L y and a single target that can be reached by the tracer molecule either from the left or from the right side. Let P (y, t) denote the probability of finding the particle at a time t at a position y, given that it has not been absorbed either by the target at the left edge, nor by the one at the right edge of the channel. Then, P (y, t) obeys the discrete diffusion equation, Eq. (B5), with the boundary conditions P (0, t) = P (L y , t) = 0, Imposing the boundary and initial conditions on the general solution, Eq. (B7), yields P (y, t) = 2 L y Ly−1 n=1 e −ω(kn)t sin (k n y 0 ) sin (k n y) . (C2) The first passage time to the target is related to the survival probability that the tracer did not yet reach the target at time t, H(t), by [52]: The survival probability is equal to the probability that the tracer diffusing on the channel of length L y with absorbing boundary conditions, remains in the system: Thus, after assuming that the particle starts its diffusive process form the middle of the channel, i.e., y 0 = L y /2, we get that the FPT distribution is given by: The MFPT is given by t = The same result can be obtained by setting L x = L z = 1 in Eq. (B21). The MPFPT is found by solving ∂F (t) ∂t t=t * = 2 L y Ly/2 m=1 ω 2 (k 2m−1 ) e −ω(k2m−1)t * × (−1) m+1 cot π(2m − 1) 2L y = 0, (C7) which in the limit L y → ∞ yields Eq. (14).

Appendix D: Evaluation of the TSDT
The Laplace transform of the FPT distribution,F (s), is finite for all values of s with a non-negative real part, and all its poles have a negative real part. Since the number of poles is finite, the asymptotic behaviour of the FPT distribution is exponential exp −t/t . The TSDT, t, is given byt whereŝ is the pole ofF (s) with the smallest real part (in absolute value). The poles ofF (s) are found by equating its denominator to zero, and thusŝ is the smallest root of G(s), defined by G(s) = For large enough L y , the sum can be approximated by the MFPT, and thust ≈ t .