Impact of global structure on diffusive exploration of organelle networks

We investigate diffusive search on planar networks, motivated by tubular organelle networks in cell biology that contain molecules searching for reaction partners and binding sites. Exact calculation of the diffusive mean first-passage time on a spatial network is used to characterize the typical search time as a function of network connectivity. We find that global structural properties — the total edge length and number of loops — are sufficient to largely determine network exploration times for a variety of both synthetic planar networks and organelle morphologies extracted from living cells. For synthetic networks on a lattice, we predict the search time dependence on these global structural parameters by connecting with percolation theory, providing a bridge from irregular real-world networks to a simpler physical model. The dependence of search time on global network structural properties suggests that network architecture can be designed for efficient search without controlling the precise arrangement of connections. Specifically, increasing the number of loops substantially decreases search times, pointing to a potential physical mechanism for regulating reaction rates within organelle network structures.


Model
To explore search efficiency, we analytically calculate the diffusive mean first-passage time (MFPT) 17 between an initial and a target node, given the connectivity and physical length of the network edges. Particle diffusion between nodes is governed by the propagator G ij (t), which gives the probability that a particle starting at node i will be at node j after time t, without passing through the target node. G j ij ∑ is the probability that the particle has never reached a target node. The MFPT to reach the target node k is is the Laplace-transform of G ij (t). Adapting recent work 46 , the propagator is given by ij ij j where I is the identity matrix, P nm  is the Laplace-transform of the flux of particles from node n directly to a connected node m without any intervening steps to other nodes, and  Q n is the Laplace-transformed probability that a particle starting at node n has not arrived at another node. Paths that reach the target are assumed to leave the network entirely, so =  P 0 nk . Equation (2) generalizes earlier work 8 to networks with distinct, non-exponential distributions for diffusion time along each edge. It differs from first-passage time calculations which assume all node-node transitions correspond to identical time steps 13,47 or with infinitesimal time spent on edges 23 , and from the numerical integration previously used to evaluate diffusion on systems of containers connected with tubes 48 .
Inserting Eq. (2) into Eq. (1) gives the MFPT between source node i and target node k. The elements  = P s ( 0) nm correspond to the probability that a particle starting at node n will next step to node m, which depends only on the lengths ℓ nm of the connecting edges: where node m and nodes w are directly connected to node n. Similarly, Q s ( 0) n =  gives the mean first passage time for a particle to arrive at any of the directly connected nodes from node n, with where D is the particle diffusivity and nodes w are directly connected to node n (derivations in Methods).

Results
Visualizing mean first-passage times. Global mean first-passage time (GMFPT) is defined as the MFPT to a single target node averaged over all possible source nodes 49 . Figure 1a shows GMFPTs for a 'decimated' honeycomb network: a complete honeycomb network is constructed inside a circle of unit radius, and edges are removed while maintaining a single connected component. The radius of the circular domain R sets the lengthscale of the system, and the particle diffusivity D sets the time-scale. Here, and in all subsequent results with synthetic networks, the first passage times are nondimensionalized by R 2 /D. Figure 1b shows GMFPTs for a particle diffusing in an example ER network from a COS-7 cell (see Methods). Both networks in Fig. 1 have higher GMFPT for nodes nearer the network periphery, as compared to centrally located nodes. Better-connected (higher degree) nodes are found more quickly (Fig. 1b inset).
Loops and total edge length constrain TA-GMFPT. Given the substantial GMFPT variation between target nodes on each network, we define a single metric characterizing the efficiency of target search processes on a particular network. Namely, the target-averaged GMFPT (TA-GMFPT) is defined as the GMFPT averaged over all possible target nodes in the network, and we use the TA-GMFPT as a typical 'search time' hereafter.
We investigate the impact of network structure on diffusive search over planar networks with nodes placed homogeneously throughout a circular domain. Decimated honeycomb networks are generated with different node densities and numbers of edges randomly removed, keeping only networks with all nodes connected. These networks all have the same spatial extent (set by domain radius R = 1), but different connectivities and node densities.
The choice of decimated lattice planar network structures, with homogeneously distributed nodes, is motivated by suggestions that yeast mitochondrial networks are evenly spread along the cell surface 32 and ER networks in several adherent cell types span throughout the relatively flat periphery of the cell 31,50 . We choose honeycomb networks because their three-way junction structure matches ER 51 and mitochondrial 32 networks, and the 120° angles between edges at network junctions match the peak angle for ER junctions 31 . This construction enables the generation of a varied family of planar networks that connect well-distributed nodes while retaining some of the geometric and topological features of cellular network structures.
Each network is characterized by the sum of all edge lengths ('total edge length') L and the cyclomatic number 16 , which is the number of elementary cycles in the network, hereafter termed 'loop number' . Loop number is given by Γ = N e − N n + 1, with N e the number of edges and N n the number of nodes. We use the loop number as a simple measure of redundant connectivity. Figure 2a shows mean TA-GMFPT vs. total edge length L and loop number Γ, averaged over many distinct decimated honeycomb networks. Larger L increases search time, by increasing the one-dimensional volume of the search space. Higher Γ substantially decreases search time -for some L values the mean search time varies by more than an order of magnitude over the explored range of Γ.
The TA-GMFPT coefficient of variation c v (ratio of standard deviation to mean) for a given total edge length L and loop number Γ does not exceed 0.3, with typical c v substantially lower (Fig. 2b). For most L values, the c v given both L and Γ is significantly smaller than the c v given L alone (Fig. 2c), demonstrating that both total edge length and loop number are necessary to accurately predict search times on a decimated lattice network. These parameters (L and Γ) incorporate the number of nodes and edges on a network within a fixed spatial region, in a manner that highlights the network density and connectivity, respectively. Times in (a) are nondimensionalized by R 2 /D, where R is the domain radius and D is the particle diffusivity; times in (b) are given for a particle diffusivity of D = 1 μm 2 /s. Inset of (b) shows GMFPT vs. node degree for both networks, normalized by the overall target-averaged GMFPT (TA-GMFPT).
www.nature.com/scientificreports www.nature.com/scientificreports/ To establish the utility of total edge length L and loop number Γ in predicting network search times, we consider several alternate network architectures, each averaged over many distinct individual networks. Decimated Voronoi networks are generated by randomly placing points within a circle with an exclusion radius around each preceding point, constructing a Voronoi tessellation using these points, and removing edges while maintaining all nodes in the single connected component. These decimated Voronoi networks maintain the three-way junction geometry of honeycomb networks, and their mean search times are very similar to decimated honeycomb networks with the same L and Γ (Fig. 2d). We also generate decimated square networks, which have node degrees up to 4. The ratio of search times between these square networks and decimated honeycomb networks, matched by L and Γ, shows greater variation (Fig. 2e). Nonetheless, search times on these square networks are generally within 40% of comparable honeycomb networks -this deviation is small in comparison to the orders of magnitude variation in the TA-GMFPT over the network structures in Fig. 2a. Figure 2 thus highlights the importance of total edge length and loop number in determining diffusive search times over a broad variety of planar network structures with well-distributed nodes.

Search on cell biology networks.
We also analyze search times on intracellular reticulated organelle network structures. Figure 3a shows an example fluorescent image of a yeast mitochondrial network 32 . Figure 3b shows TA-GMFPT from 350 mitochondrial networks 32 , which exhibit features similar to the honeycomb networks in Fig. 2a: approximate prediction of TA-GMFPT by total edge length L and loop number Γ, and a substantial decrease in search time as Γ increases and L decreases. Mitochondrial networks from wild-type cells and mutant cells with mitochondrial fission and fusion proteins knocked out occupy distinct regions of the Γ vs. L plane. However, for given values of these two structural parameters, the two network types (wild-type and mutant) exhibit similar search times. The fluorescent image in Fig. 3d shows an example ER network. We calculated the TA-GMFPT for regions of 103 such peripheral ER networks (Fig. 3e). Although ER network structures are restricted to relatively high looping number for each total edge length, the search times vary similarly to honeycomb networks (Fig. 2a) and mitochondrial networks (Fig. 3a).
While the synthetic honeycomb, Voronoi, and square networks in Fig. 2 are planar, mitochondrial and ER networks exist in three-dimensional intracellular volumes. However, imaging of peripheral ER networks in COS7 cells indicates that these structures are relatively flat, with rarely observed crossing of tubules outside the typical 3-way junction nodes 42,45,50 . Three dimensional deformation in the paths of individual edges would reduce the overall effective diffusion coefficient in a planar projection 52 , but would not substantially alter the global search trends described here. Mitochondrial networks in budding yeast cells tend to remain at the cell surface, with little incursion into the three-dimensional bulk of the cell 32 . These networks are thus essentially confined to a two-dimensional manifold in the shape of a spherical shell. For simplicity we also approximate them as planar, neglecting the large-scale curvature of the spherical surface. www.nature.com/scientificreports www.nature.com/scientificreports/ Investigation of how network structural characteristics facilitate searches must account for the expected search time increase as domain size increases. To compare the mitochondrial and ER networks to idealized lattice-like structures, the cellular networks are scaled to the same physical area as the unit circle containing honeycomb networks. Three-dimensional mitochondrial network coordinates extracted from imaging are projected onto a spherical surface. The network area is estimated using a convex hull of both ER nodes on the plane and projected mitochondrial nodes on the sphere (details in Methods). This allows comparison of search times between networks with the same spatial extent but different node density and connectivity.
Figures 3c,f plot the ratio of each mitochondrial and ER network search time, respectively, to the honeycomb network search time at the corresponding total edge length and loop number. The ratios for both mitochondrial and ER networks are near unity, suggesting these organelle network structures have similar search characteristics and dependence on total edge length and loop number as synthetic lattice-like networks.
Dependence of search times on network morphology. Using decimated honeycomb and square lattice network structures (Fig. 2), we explore the dependence of search time on total edge length L and loop number Γ (Fig. 4). Increasing loop number at a constant total edge length corresponds to networks with less dense nodes that are more completely connected (Fig. 4a). The search time varies distinctly for low vs. high loop numbers Γ. Search time weakly depends on Γ for low Γ (Fig. 4b,c), indicating that adding a few loops, within a primarily tree-like structure, will not substantially affect the search process. In contrast, search time steeply decreases with rising Γ for higher Γ values (Fig. 4b,c). This suggests that once a threshold number of loops is reached, further added loops can significantly decrease search time. High search time variability in Fig. 2b,d,e aligns with the neighborhood of these thresholds in Fig. 4b,c, suggesting that at the threshold where loop number begins to perturb global transport, the precise arrangement of the loops can have a substantial impact on search time.
By contrast, increasing total edge length L at a constant loop number corresponds to denser, less well connected networks (Fig. 4d). The dependence of search time on total edge length L becomes more steep as loop number increases (Fig. 4e,f).
The decimated honeycomb and square lattice networks resemble percolation systems, where the fraction of bonds retained is above the critical percolation value p c . Random walks in such systems are effectively diffusive above a certain correlation length, and for planar networks should have the same scaling properties as two-dimensional diffusion 53,54 . In particular, we expect the search time to be largely independent of node density and to scale as T ~ D −1 , where D is the effective diffusivity (see Methods for details). Near the percolation threshold, where p is the fraction of lattice bonds remaining and μ ≃ 1.30 for two-dimensional lattices 54 . www.nature.com/scientificreports www.nature.com/scientificreports/ By treating our synthetic networks as large clusters in a two-dimensional system approaching percolation, we derive the expected dependence of search times on total edge length and loop number (see Methods section). Namely, when loop number is very low (Γ ≪ L), then search times are expected to be independent of Γ but to scale with total edge length as T ~ L μ . Both of these expected relationships are consistent with the observed dependence of search times on L and Γ for synthetic networks with low loop numbers (Fig. 4).
We note that the L 1.30 dependence is intermediate between two extreme cases of loopless networks within a fixed-area domain. One extreme includes linear structures that snake through the domain without branching, or comb-like networks with a single backbone connecting many individual branches, which both exhibit MFPTs scaling as ~L 217 . The other extreme is self-similar tree-like networks 17,55 , which have MFPTs that vary as ~ L when scaled down to unit physical extent (see Methods).
For high loop numbers (Γ ≫ L), search times in a cluster close to the percolation transition are expected to depend on both total edge length and loop number as T ~ L 2μ Γ −μ (see Methods). In Fig. 4b,c, this T ~ Γ −1.30 dependence is consistent with search times for low total edge lengths, but does not entirely explain search time behavior at the highest total edge lengths and high loop numbers. Similarly, in Fig. 4e,f, the predicted T ~ L 2.60 dependence is consistent with search times for intermediate loop numbers, but does not entirely explain search time behavior for the highest loop numbers. The transition from intermediate connectivity to poor connectivity (Γ ≲ L) is evident in the top right corner of Fig. 4e,f, where the search time dependence on total edge length begins to shift from T ~ L 2μ to T ~ L μ . We note that networks with the lowest total edge length for a given loop number (or highest loop number for a given total edge length) correspond to the most fully connected lattices, which are far from the percolation transition. We would thus expect the aforementioned scaling relationships to break down in this regime, as is seen for the data points with highest L and Γ in Fig. 4c, and for the lowest L, highest Γ in Fig. 4e,f.

Discussion
We have investigated the characteristics that control diffusive search time on planar networks connecting homogeneously distributed nodes over a compact domain. To this end, we employ an exact calculation of mean first-passage time on a spatial network (Eqs.(1)-(4)) based on network connectivity and edge lengths.
Search times in a complex medium are known to depend not only on the spatial structure of the domain but also on the dynamic nature of the search process, with the dimensionality of the walk (defined by www.nature.com/scientificreports www.nature.com/scientificreports/ determining a phase transition between compact and non-compact search processes 19,22,23,56 . In the analysis presented here, changing the dimensionality for particle dynamics would alter how splitting probabilities and waiting times for node-to-node transitions depend on edge length (Eqs. (3), (4)), as well as likely modifying the scaling behavior of search times with network structure near the percolation transition (as discussed in Methods). Such extension to different dynamical processes is outside the scope of this work and is left as a fruitful area for further study. Throughout our calculations, particle motion along network edges is assumed to be purely diffusive (d w = 2). This assumption of diffusive transport is consistent with past analyses of large-scale particle spreading on both ER and mitochondrial structures, as measured by fluorescence recovery after photobleaching 44,[57][58][59] . More recent single-particle tracking studies in the ER indicate that membrane proteins move diffusively, while luminal proteins may in some cases be driven by random processive flows along the network edges 45 . Subdiffusive behavior (d w > 2), commonly attributed to fractional Brownian motion in a viscoelastic medium, has been observed for a variety of intracellular particles of size comparable to organelles or RNA-protein complexes (r ≳ 50 nm) [60][61][62] . However, smaller particles such as individual proteins (~5 nm) often exhibit diffusion-like motion 63 . These observations motivate our choice to focus on diffusive exploration for proteins in the ER and mitochondrial networks.
We assess typical search time on each network by averaging all combinations of source and target nodes. Diffusive search time on networks with homogeneously distributed nodes, including both synthetic networks and those from intracellular structures, is found to be largely predicted by simple geometric characteristics: total edge length and loop number (Figs. 2, 3), which characterize network density and connectivity. Increasing loop number substantially decreases diffusive search time, while increasing total edge length can steeply increase the search time (Figs. 2-4). Search times on ER and mitochondrial networks are comparable to those computed for idealized planar lattice structures with equivalent loop number and edge length (Fig. 3c,f), emphasizing the sufficiency of these two global structural parameters for determining diffusive search efficiency on real-world networks.
Using percolation theory to predict diffusivity, which is inversely related to the search time, largely describes the dependence of search time on total edge length and loop number for networks constructed from decimated planar lattices. This link to percolation theory highlights the importance of network connectivity, in the form of the bond fraction p, for determining search times. Although the bond fraction is well-defined for idealized regular planar lattices, it is not a meaningful parameter for realistic cellular networks, such as ER and mitochondrial structures. For these and other off-lattice networks, we show that the simply measurable parameters of total edge length and loop number can be used to predict diffusive search behavior. We have outlined how these network parameters can be connected to an effective bond fraction and thus to the wide range of results available for percolation systems.
Typical mitochondrial and ER networks have many loops, which accelerate search (Fig. 3b,e). These search-accelerating loops align with Murray's law for vasculature radius 64,65 or the balance of competing constraints in fungi 66,67 , suggesting biology is capable of optimizing transport networks and that networks with many loops may have been partly selected for efficient diffusive transport.
Mitochondrial and endoplasmic reticulum networks extend through much of the cellular volume to interact and form direct contacts with other organelles [68][69][70] , providing connections between various subcellular systems. A well-connected reticulated network structure enables proteins within these organelles to explore many contact sites without requiring export or relying on the slower dynamics of the organelles themselves. Mitochondrial and endoplasmic reticulum network morphology can vary based on cell specialization and the functional state of the organelle 71,72 , possibly modulating signals carried by diffusion through mitochondrial and endoplasmic reticulum networks. Our results indicate that higher loop numbers in the network would decrease diffusive search times and thereby may speed transmission of diffusive signals to organelle contact sites. This potential for organelle functional and signaling modulation through network structural properties is similar to the influence of topology and significance of time delay for organismal states described by network physiology 4,5,73 .
Our finding that loops speed diffusive search points towards key structural criteria for spatial networks whose function relies on efficient diffusive transport, including the intracellular networks studied here. Earlier work on random walks in complex networks showed that tree networks maximize the TA-GMFPT (i.e., lead to the slowest search) 23 , which suggested that some loops may lead to more efficient diffusive search.
By contrast to networks designed for diffusive transport, the optimal spatial network structure for potential-driven flow in a variety of scenarios is a loopless tree 2,74-77 . However, loops can assist network flow-based transport outside of steady state, providing resiliency to damage and fluctuations [77][78][79] . Mitochondrial 43 and ER networks 29 are very dynamic, and resiliency to edge removal may be another benefit of the many loops in these cell biology networks. Although we do not include edge or loop production cost 78,79 , our analysis can establish the utility of these network components for improving transport.
The connection between network structure and diffusive distribution efficiency indicates a potential link between architecture and functionality for cellular organelles such as the ER and mitochondria. The dependence of search efficiency on global structural properties of the network suggests that cells may be able to regulate biochemical kinetics without precise local arrangement of network connections.

Methods
Deriving P and Q. Equations (1) and (2) give the MFPT between source node i and target node k, , we consider, as an example, a particle at a degree-three node 0 with edges of length ℓ 1 ≤ ℓ 2 ≤ ℓ 3 connecting to other nodes (1,2,3). Trajectories that reach node 1 before nodes 2 or 3 can be constructed from excursions a distance ℓ 1 from the initial node, with first-passage time distribution For a diffusing particle that starts at position ℓ 1 on an interval with absorbing boundaries at x = 0 and x = d, the function f d (t) gives the total flux out of the interval at time t and the function F d (t) gives the flux at the x = 0 boundary. The first term of Eq. (6) represents a trajectory that first reaches a distance ℓ 1 from the initial node when it arrives at node 1. The second term of Eq. (6) is for a particle that reaches ℓ 1 from the initial node along the edge to node 2, returns to the initial node without first reaching node 2, and then diffuses to node 1. A Laplace transform t → s converts the convolutions over sequential steps into products, giving  More generally,  P nm and Q n  are where node m and nodes w are directly connected to node n.
network and removing a set of edges, subject to the condition that all initial nodes remain attached to all other nodes in a single connected component when edges are removed. Many network variations can be constructed from one complete network by varying the number and identity of removed edges. For honeycomb and square networks, initial complete networks are constructed as a lattice within a circle of radius one, with nearest neighbors connected by an edge. The lattice size is varied to obtain complete networks with different node densities. 6.2 × 10 4 honeycomb networks are generated for data in Fig. 2a, 3.6 × 10 4 Voronoi networks for Fig. 2d, and 6.5 × 10 4 square networks for Fig. 2e.
To construct a Voronoi network, we first randomly place points within a circle of radius one, subject to the condition that each subsequent point cannot be within an exclusion radius of all preceding points. When nodes can no longer be placed (as the entire circle is blocked with the exclusion radius of at least one point), a Voronoi tesselation is constructed around these points. The boundaries of the Voronoi tesselation cells form the network. The exclusion radius around the initial points is varied to obtain networks with different node densities.
Mitochondrial networks. Spatial coordinates and network connections for 350 mitochondrial networks from Saccharomyces cerevisiae budding yeast cells, obtained using Mitograph software, were generously provided by Matheus Viana and Susanne Rafelski 32 . The networks we analyze include wild-type cells and Δdnm1Δfzo1 mutant cells lacking proteins for mitochondrial fission and fusion. For each cell, we used the largest connected component.
Mitochondrial networks have relatively few nodes and edges in comparison to the synthetic networks. Nodes with degree two were added along edges to ensure individual edge lengths were approximately homogeneous, facilitating comparison with decimated lattice networks. Specifically, sufficient nodes were added to make all individual node-to-node edges shorter than the shortest full edge in the original network, and shorter than a 1 μm length ceiling. This procedure does not change the geometry or topology of the original network, but does redefine the set of target nodes used for the calculation of the TA-GMFPT.
Endoplasmic reticulum networks. COS-7 cells were purchased from ATCC (Catalog # ATCC-CRL1651) and were grown in Dulbecco's modified Eagle medium (DMEM) supplemented with 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin (P/S). Prior to imaging experiments, COS-7 cells were seeded in 6-well, plastic bottom dishes at 1 × 10 5 cells/mL about 18 hours prior to transfection. Plasmid transfections were performed as described previously 80 . For all imaging experiments, the ER was fluorescently labeled with 0.2 μg KDEL venus transfected into each well of a 6-well dish 81 . Live cells were imaged at 37 °C in Fluorobrite imaging media (Invitrogen) supplemented with 10% FBS. Confocal Z-stack images of the peripheral ER were collected using Micromanager Imaging Software with a step size of 0.2 μm. All images were acquired on an inverted fluorescent microscope (TE-2000-U; Nikon) equipped with a Yokogawa spinning-disk confocal system (CSU-Xm2; Yokogawa CSU X1) 80 . Images were taken with a 100 × NA 1.4 oil objective on an electron-multiplying charge-coupled device (CCD) camera 50 × 50 (Andor). Images were acquired with Micromanager Imaging Software and then analyzed, merged and contrasted using Fiji (ImageJ) 82 .
A large continuous region of the peripheral endoplasmic reticulum network was selected from each image. The endoplasmic reticulum from this region of each image was skeletonized, and node and edge data from the skeleton extracted, using Fiji (ImageJ). Node and edge data was analyzed to extract a network structure, assuming nodes within 0.001 μm of one another are the same node. For each cell, the largest connected component was used. We obtained 103 ER networks. Fig. 2c was obtained.

Coefficient of variation. This section describes how
In Fig. 2b, network structures are sorted into bins according to their total edge length (bin size of 2) and loop number (bin size of 4). For each bin there is a mean search time m ij (mean of the TA-GMFPTs for all networks falling into the bin) and a variance of the search time ij 2 σ , where i, j indicate the bin indices for total edge length and loop number, respectively.
In the Fig. 2c the red curve labeled 'Total edge length only' is the coefficient of variation over all loop numbers given a total edge length bin i. This coefficient depends both on the variance within individual bins and the overall variability from bin to bin. It is given by www.nature.com/scientificreports www.nature.com/scientificreports/ Network size scaling. To make a direct comparison between search times for synthetic networks constrained to a circle of radius one, and search times for networks from cell biology, we scale lengths in the cellular networks such that the effective area spanned by the network matches the synthetic network area of π (circle of radius one). Search times are scaled by the length scaling factor squared, as diffusive processes in one dimension occur in a time proportional to length squared.
Three-dimensional points along the largest connected component of each mitochondrial network skeleton are projected onto a sphere, whose center and radius are set to minimize the mean square residual of network points from the surface of that sphere. A convex hull of points is then constructed from these projected positions on the sphere, using the convhulln routine in Matlab, yielding a set of triangles. Triangles are rejected if their center is more than 0.3 μm from the sphere surface or the orientation of their normal vector is more than 40° from the radial direction. This procedure effectively removes triangles spanning across large sphere regions not covered by the mitochondrial network. The areas of the remaining triangles are summed and used as an effective area spanned by the mitochondrial network.
For the ER structures, a convex hull is found from the two-dimensional points along the largest connected component of each network. The total area of the convex hull is then used for the effective area of the endoplasmic reticulum network.
Approximating search times with percolation theory. We consider a fully connected n × n square lattice of network nodes on a unit square, giving N = n 2 total nodes connected to nearest neighbors by edges of length ℓ = 1/(n − 1). The complete lattice has E nn 2 ( 1) max = − total edges. Edges are removed from the network until a number Γ of loops remain, without disconnecting any nodes. The number of edges remaining is E = n 2 − 1 + Γ. The fraction of edges that remain is p = E/E max , p n nn 1 2 The critical bond probability for percolation on a square lattice is p c = 1/2 54 , such that our lattice has p > p c . For p > p c but remaining near p c , the diffusivity depends on the bond probability as D p p ( ) c − μ , where μ ≃ 1.30 54 . The diffusivity on the network is Note that for a fixed loop number Γ, increasing the total edge length corresponds to increasing the lattice density n, and hence decreasing the effective diffusivity. Specifically, the total edge length for the network is given by L = ℓE = n + 1 + Γ∕(n − 1). We can then express the lattice density in terms of our control parameters L and Γ according to, Assuming a dense lattice system with L ≫ 1 and Γ ≪ L 2 , we get the scaling n ~ L. This can be plugged into the diffusivity (Eq. (18b)) to show For a highly disconnected system with low loop number (Γ ≪ L), the diffusivity scales as D ~ Γ 0 L −μ . In the opposite extreme of high loop number (L ≪ Γ ≪ L 2 ), the diffusivity scales as D ~ Γ μ L −2μ . Random walks on the largest connected component of a planar lattice above the percolation transition are expected to show the universal scaling behavior associated with diffusion in two dimensions, at sufficiently large length scales 54 . The target-site search time for a two-dimensional random walk with unit time steps is known to scale with the number of sites T step ~ N, neglecting a logarithmic correction term 19,83 . For our system, the particles diffuse with diffusivity D along edges of length ℓ, so that the characteristic time to traverse each edge scales as ℓ 2 ∕D. Consequently, the overall expected search time is We thus expect the search time to vary as T ~ L μ for nearly loop-less networks, and T ~ Γ −μ L 2μ for networks with high loop numbers.
A similar argument can be used to relate the first passage time T and the effective diffusivity D for a compact search process where the underlying particle dynamics is subdiffusive, with The coefficient D in this case characterizes the large-scale spreading of particles over the network structure. The dependence of D on network connectivity in a cluster near percolation (i.e.: the scaling exponent μ) is likely to be altered for subdiffusive motion. However, we do not address this behavior here, focusing instead on diffusive search processes that are expected to be relevant for a variety of proteins in the ER and mitochondrial networks.
Search times on simple loop-less network topologies. Self-similar, hierarchically branched tree networks are constructed iteratively by attaching m additional branches to the center of each branch in an existing tree. The number of steps S to find a central target on a tree generated by g iterations scales as g g m 1 log 2/log ( 2) + + where N g = (m + 2) g + 1 is the number of nodes in the tree 55 . The number of tree edges is given by K g = (m + 2) g , so K g ~ N g .
If such a hierarchical tree network is constrained to a domain of unit radius, the edge lengths of the tree must become shorter with each iteration, scaling as ℓ g = 2 −g . The total edge length will then be L g = K g ℓ g ~ [(m + 2)/2] g . The time required to diffuse across each edge is  (  2) log 2 log (  2) Consequently, the time for diffusive search over a fractal tree network scaled to fit within a domain of fixed spatial extent should scale as T ~ L, as indicated in the main text.

Data availability
Datasets generated and analysed during the current study are available from the corresponding author upon request. Software for computing mean first passage times is available in a GitHub repository at: https://github. com/lenafabr/networkMFPT.