Impact of global structure on diffusive exploration of organelle networks

TL;DR

This study analyzes how the overall structure of organelle networks influences diffusive search efficiency, revealing that global properties like total edge length and loops significantly impact search times, with implications for cellular function and network design.

Contribution

The paper demonstrates that global structural properties of planar networks largely determine diffusive search times, connecting network topology with physical models and biological regulation.

Findings

Search time depends on total edge length and number of loops.

Increasing loops reduces search times significantly.

Global network properties can be used to optimize search efficiency.

Abstract

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.

Equations

$$⟨T_{ik}⟩={\int }_0^{\infty }t\left[-\frac{\partial }{\partial t}\sum _{j\ne k}{G}_{ij}(t)\right]dt=\sum _{j\ne k}{\hat{G}}_{ij}(s=0).$$

$${\hat{G}}_{ij}(s)=[(\mathbf{I}-\hat{\mathbf{P}}{)}^{-1}{]}_{ij}{\hat{Q}}_j,$$

$${\hat{P}}_{nm}(s=0)=\frac{{\ell }_{nm}^{-1}}{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nw}^{-1}},$$

$${\hat{Q}}_n(s=0)=\frac{1}{2D}\frac{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nm}}{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nm}^{-1}},$$

$$⟨T_{ik}⟩={\left[{\left(\mathbf{I}-\hat{\mathbf{P}}\right)}^{-1}\hat{\vec{Q}}\right]}_{i,s=0}.$$

$$P_1(t)=\frac{1}{3}{f}_{2{\ell }_1}(t)+{\int }_0^{t}d{t}_2{\int }_0^{t_2}d{t}_1$$

$$\frac{1}{3}{f}_{2{\ell }_1}(t-t_2)F_{{\ell }_2}(t_2-t_1)\frac{1}{3}{f}_{2{\ell }_1}(t_1)+\dots .$$

$${\hat{P}}_1$$

$$\displaystyle+\frac{1}{3}\hat{f}_{2\ell_{1}}(\hat{F}_{\ell_{2}}+\hat{F}_{\ell_{3}})\left[\frac{1}{3}\hat{f}_{2\ell_{1}}+\ldots\right]\bigg{]}$$

$$=\frac{1}{3}{\hat{f}}_{2{\ell }_1}\sum _{n=0}^{\infty }{\left[\frac{1}{3}{\hat{f}}_{2{\ell }_1}({\hat{F}}_{{\ell }_2}+{\hat{F}}_{{\ell }_3})\right]}^n$$

$$\displaystyle=\frac{1}{3}\hat{f}_{2\ell_{1}}\bigg{/}\left[1-\frac{1}{3}\hat{f}_{2\ell_{1}}(\hat{F}_{\ell_{2}}+\hat{F}_{\ell_{3}})\right]\ ,$$

$$\hat{f}_{d}(s)=2\sinh\left[\sqrt{s/D}(d-\ell_{1})\right]\Big{/}\left[\sinh(\sqrt{s/D}d)\right]\ ,$$

$$\hat{F}_{d}(s)=\sinh\left[\sqrt{s/D}(d-\ell_{1})\right]\Big{/}\left[\sinh(\sqrt{s/D}d)\right]\ ,$$

$${\hat{P}}_1(s=0)={\ell }_2{\ell }_3/\left({\ell }_1{\ell }_2+{\ell }_2{\ell }_3+{\ell }_1{\ell }_3\right).$$

$$Q_j(t)=1-\sum _{w=1}^3{P}_{jw}(t),$$

$${\hat{Q}}_j(s)=\frac{1}{s}\left[1-\sum _1^3{\hat{P}}_{jw}(s)\right].$$

$$\begin{array}{rl}{\hat{Q}}_j(s=0)& =-\frac{\partial }{\partial s}\sum _1^3{\hat{P}}_{jw}(s=0)\\ & =\frac{1}{2D}\frac{{\ell }_1+{\ell }_2+{\ell }_3}{{\ell }_1^{-1}+{\ell }_2^{-1}+{\ell }_3^{-1}}.\end{array}$$

$${\hat{P}}_{nm}(s=0)=\frac{{\ell }_{nm}^{-1}}{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nw}^{-1}},$$

$${\hat{Q}}_n(s=0)=\frac{1}{2D}\frac{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nm}}{{\sum }_{w=1}^{\text{deg}(n)}{\ell }_{nm}^{-1}},$$

$$\begin{array}{rl}{c}_{\text{v},i}& =\frac{1}{{⟨m⟩}_i}\sqrt{{⟨{\sigma }^2⟩}_i+{\text{var}}_i(m)},\end{array}$$

$$c_{\text{v},i}=\frac{1}{n_i}\sum _j\frac{{\sigma }_{ij}}{m_{ij}}.$$

$$p=\frac{1}{2}+\frac{1}{2n}+\frac{\Gamma }{2n(n-1)}.$$

$$D\sim$$

$$\sim$$

$$n=\frac{L}{2}\left(1+\sqrt{1-\frac{4(L+\Gamma -1)}{L^2}}\right).$$

$$D\sim$$

$$T\sim \frac{{\ell }^2}{D}N\sim \frac{(n-1{)}^{-2}}{D}{n}^2\sim \frac{1}{D}.$$

$$T\sim \frac{{\ell }^{2/\alpha }}{D}{N}^{\alpha /2}\sim \frac{(n-1{)}^{-2/\alpha }}{D}{\left(n^2\right)}^{\frac{2/\alpha }{2}}\sim \frac{1}{D}.$$

$$S_g\sim N_g^{1+\log 2/\log (m+2)},$$

$$\begin{array}{rl}{T}_g=\Delta t_g{S}_g& \sim {\ell }^{1-\frac{\log 2}{\log (m+2)}}{L}_g^{1+\frac{\log 2}{\log (m+2)}}\\ & \sim L_g\left[2^{-g+\frac{g\log 2}{\log (m+2)}}\right(\frac{m+2}{2}\left){}^{\frac{g\log 2}{\log (m+2)}}\right]\\ & \sim L_g\end{array}$$

Full text

††preprint: APS/123-QED

Network models have been employed to describe and understand a wide variety of phenomena Boccaletti et al. (2006), ranging from transparently physical processes such as flow Durand (2007) and conductivity Cheianov et al. (2007) to the more abstract examples of human physiology Bashan et al. (2012); Ivanov et al. (2016), social interactions Wasserman and Faust (1994), and mortality Farrell et al. (2018). Random walks on networks Masuda et al. (2017) can model the dynamics of epidemic spreading Pu et al. (2015), animal foraging Perna and Latty (2014), brain signaling Chavez et al. (2010), and electron transport Nelson (1999). The rate at which such random walkers find target sites within the networks is known to depend on factors such as dimensionality Burioni and Cassi (2005), target connectivity Hwang et al. (2012), and number of shortest paths passing through the target Noh and Rieger (2004).

In comparison to generalized complex networks, spatial networks have physical constraints that limit connections to spatially proximal nodes Barthelemy (2011). In addition, physical diffusion along network edges gives rise to broadly distributed non-exponential waiting times that depend on edge length Redner (2001), in contrast to classic models of random walks on networks Noh and Rieger (2004); Barthelemy (2011); Masuda et al. (2017). Thus, diffusive search on spatial networks is perhaps better described by physical variables such as fractal dimensionality Ben-Avraham and Havlin (2000); Benichou et al. (2010) or tortuosity Shen and Chen (2007), rather than the number of nodes and edges often used to characterize general network structures. Random walks on spatial networks have similar dynamical properties to those in complex, porous, or crowded media Havlin and Ben-Avraham (1987); Condamin et al. (2007); Benichou and Voituriez (2014). Although there has been progress in calculating how geometry affects diffusive search times on complex domains Condamin et al. (2007); Benichou et al. (2010); Benichou and Voituriez (2014), and understanding the impact of distinct search strategies Chupeau et al. (2015) there is little guidance on how to structure spatial networks to accelerate diffusive search.

Intracellular structures provide a key example where diffusive processes over complex geometries have an important role to play in cellular function. Reticulated organelles, such as the peripheral endoplasmic reticulum (ER) Westrate et al. (2015); Schwarz and Blower (2016) and mitochondria Collins et al. (2002); Rafelski et al. (2012), are composed of membranous tubules enclosing a single connected luminal volume, whose physical structure has recently been mapped in extensive detail Speckner et al. (2018); Harwig et al. (2018). These organelles constitute spatial networks that span throughout the cell interior, comprising hundreds of nodes and edges connected in a highly looped architecture Lin et al. (2014); Viana et al. (ress). Within these networks, proteins and other molecules diffuse to find reaction partners and binding targets. For instance, secretory proteins must encounter an exit site in order to leave the ER Hughes et al. (2009) and DNA-binding proteins must find mitochondrial nucleoids to participate in DNA maintenance and replication Ruhanen et al. (2010).

The structure of these living networks is heavily regulated and likely functionally important English and Voeltz (2013a); Viana et al. (ress). Mitochondrial network structure changes Rube and van der Bliek (2004) during the cell cycle Margineantu et al. (2002), differentiation Shin et al. (2016), and disease Willems et al. (2009), suggesting mitochondrial morphology plays a role in physiological functions such as ATP production Chiaradonna et al. (2006); Ghosh et al. (2018). ER structure varies with cell specialization Schwarz and Blower (2016) and with mutations in morphogenic proteins associated with human pathologies Chen et al. (2013); Westrate et al. (2015). Prior work analyzed the structure and morphogenesis of mitochondrial Viana et al. (ress); Zamponi et al. (2018) and ER networks Lin et al. (2014), and mapped some basic parameters of molecular diffusion within these networks Dayel et al. (1999); Holcman et al. (2018); Viana et al. (ress). However the connection between network morphology and search efficiency has not been systematically addressed.

We investigate diffusive search on two network types. Firstly, we construct a variety of synthetic planar networks, with nodes on a lattice or homogeneously scattered and connected into a single component with varying arrangements of edges. Secondly, we use spatial networks extracted from imaging of yeast mitochondria and mammalian ER. We find that typical search times on these biological structures are largely predicted by simple global structural parameters: the total edge length and loop number, which encompass network density and connectivity.

I Model

To explore search efficiency, we analytically calculate the diffusive mean first-passage time (MFPT) Redner (2001) 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. $\sum_{j}G_{ij}$ is the probability that the particle has never reached a target node. The MFPT to reach the target node $k$ is

[TABLE]

where $\hat{G}_{ij}(s)$ is the Laplace-transform of $G_{ij}(t)$. Adapting recent work Koslover and Spakowitz (2012), the propagator is given by

[TABLE]

where $\bm{I}$ is the identity matrix, $\hat{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 $\hat{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 $\hat{P}_{nk}=0$. Eq. 2 generalizes earlier work Masuda et al. (2017) 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 Burioni and Cassi (2005); Maier and Brockmann (2017) or with infinitesimal time spent on edges Benichou and Voituriez (2014), and from the numerical integration previously used to evaluate diffusion on systems of containers connected with tubes Lizana and Konkoli (2005).

Inserting Eq. 2 into Eq. 1 gives the MFPT between source node $i$ and target node $k$. The elements $\hat{P}_{nm}(s=0)$ correspond to the probability that a particle starting at node $n$ will next step to node $m$, which depends only on the lengths $\ell_{nm}$ of the connecting edges:

[TABLE]

where node $m$ and nodes $w$ are directly connected to node $n$. Similarly, $\hat{Q}_{n}(s=0)$ gives the mean first passage time for a particle to arrive at any of the directly connected nodes from node $n$, with

[TABLE]

where $D$ is the particle diffusivity and nodes $w$ are directly connected to node $n$ (derivations in Methods).

II Results

II.1 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 Tejedor et al. (2009). 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 length-scale 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).

II.2 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 Viana et al. (ress) and ER networks in several adherent cell types span throughout the relatively flat periphery of the cell Lin et al. (2014); Friedman and Voeltz (2011). We choose honeycomb networks because their three-way junction structure matches ER Shemesh et al. (2014) and mitochondrial Viana et al. (ress) networks, and the $120^{\circ}$ angles between edges at network junctions match the peak angle for ER junctions Lin et al. (2014). 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 Barthelemy (2011), which is the number of elementary cycles in the network, hereafter termed ‘loop number’. Loop number is given by $\Gamma=N_{\text{e}}-N_{\text{n}}+1$, with $N_{\text{e}}$ the number of edges and $N_{\text{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 $\Gamma$, averaged over many distinct decimated honeycomb networks. Larger $L$ increases search time, by increasing the one-dimensional volume of the search space. Higher $\Gamma$ 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 $\Gamma$.

The TA-GMFPT coefficient of variation $c_{\text{v}}$ (ratio of standard deviation to mean) for a given total edge length $L$ and loop number $\Gamma$ does not exceed 0.3, with typical $c_{\text{v}}$ substantially lower (Fig. 2b). For most $L$ values, the $c_{\text{v}}$ given both $L$ and $\Gamma$ is significantly smaller than the $c_{\text{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 $\Gamma$) 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.

To establish the utility of total edge length $L$ and loop number $\Gamma$ 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 $\Gamma$ (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 $\Gamma$, 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.

II.3 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 Viana et al. (ress). Figure 3b shows TA-GMFPT from 350 mitochondrial networks Viana et al. (ress), which exhibit features similar to the honeycomb networks in Fig. 2a: approximate prediction of TA-GMFPT by total edge length $L$ and loop number $\Gamma$, and a substantial decrease in search time as $\Gamma$ 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 $\Gamma$ 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 Figure 3d shows an example ER network. We calculated the TA-GMFPT for regions of 103 such peripheral ER networks (Figure 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 Friedman and Voeltz (2011); Chen et al. (2013); Holcman et al. (2018). Three dimensional deformation in the paths of individual edges would reduce the overall effective diffusion coefficient in a planar projection Adler et al. (2019), 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 Viana et al. (ress). 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.

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 and 3f 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.

II.4 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 $\Gamma$ (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 $\Gamma$. Search time weakly depends on $\Gamma$ for low $\Gamma$ (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 $\Gamma$ for higher $\Gamma$ 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_{\text{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 Gefen and Aharony (1983); Stauffer and Aharony (1994). In particular, we expect the search time to be largely independent of node density and to scale as $T\sim D^{-1}$, where $D$ is the effective diffusivity (see Methods for details). Near the percolation threshold, $D\sim(p-p_{\text{c}})^{\mu}$, where $p$ is the fraction of lattice bonds remaining and $\mu\simeq 1.30$ for two-dimensional lattices Stauffer and Aharony (1994).

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 ($\Gamma\ll L$), then search times are expected to be independent of $\Gamma$ but to scale with total edge length as $T\sim L^{\mu}$. Both of these expected relationships are consistent with the observed dependence of search times on $L$ and $\Gamma$ 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 $\sim L^{2}$ Redner (2001). The other extreme is self-similar tree-like networks Redner (2001); Lin et al. (2010), which have MFPTs that vary as $\sim L$ when scaled down to unit physical extent (see Methods).

For high loop numbers ($\Gamma\gg 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\sim L^{2\mu}\Gamma^{-\mu}$ (see Methods). In Fig. 4b,c, this $T\sim\Gamma^{-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\sim 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 ($\Gamma\lesssim 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\sim L^{2\mu}$ to $T\sim L^{\mu}$. 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 $\Gamma$ in Fig. 4c, and for the lowest $L$, highest $\Gamma$ in Fig.  4e,f.

III 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 $\left<x^{2}\right>\sim t^{2/d_{\text{w}}}$) determining a phase transition between compact and non-compact search processes Carretero-Campos et al. (2012); Condamin et al. (2007); Benichou et al. (2010); Benichou and Voituriez (2014). 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_{\text{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 Partikian et al. (1998); Dayel et al. (1999); Sbalzarini et al. (2005); Dieteren et al. (2011). 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 Holcman et al. (2018). Subdiffusive behavior ($d_{\text{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\gtrsim 50$nm) Yamada et al. (2000); Tolić-Nørrelykke et al. (2004); Lampo et al. (2017). However, smaller particles such as individual proteins ($\sim 5$nm) often exhibit diffusion-like motion Etoc et al. (2018). 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 (Figs. 3b,e). These search-accelerating loops align with Murray’s law for vasculature radius Sherman (1981); McCulloh et al. (2003) or the balance of competing constraints in fungi Tero et al. (2009); Heaton et al. (2012), 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 Klecker et al. (2014); Phillips and Voeltz (2016); Valm et al. (2017), 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 Benard and Rossignol (2008); Schuck et al. (2009), 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 Bashan et al. (2012); Bartsch et al. (2015); Ivanov et al. (2016).

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) Benichou and Voituriez (2014), 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 Banavar et al. (2000); Durand (2007); Bohn and Magnasco (2007); Bernot et al. (2009); Hu and Cai (2013). However, loops can assist network flow-based transport outside of steady state, providing resiliency to damage and fluctuations Katifori et al. (2010); Corson (2010); Hu and Cai (2013). Mitochondrial Zamponi et al. (2018) and ER networks Speckner et al. (2018) 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 Katifori et al. (2010); Corson (2010), 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.

IV Methods

IV.1 Deriving P and Q

Equations 1 and 2 give the MFPT between source node $i$ and target node $k$,

[TABLE]

To find $\hat{P}_{ij}$ and $\hat{Q}_{j}$, we consider, as an example, a particle at a degree-three node with edges of length $\ell_{1}\leq\ell_{2}\leq\ell_{3}$ connecting to other nodes. Trajectories that reach node 1 before nodes 2 or 3 can be constructed from excursions a distance $\ell_{1}$ from the initial node, with first-passage time distribution

[TABLE]

For a diffusing particle that starts at position $\ell_{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. IV.1 represents a trajectory that first reaches a distance $\ell_{1}$ from the initial node when it arrives at node 1. The second term of Eq. IV.1 is for a particle that reaches $\ell_{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\rightarrow s$ converts the convolutions over sequential steps into products, giving

[TABLE]

with the last line from the infinite summation of a geometric series.

The Laplace-transforms of $f_{d}(t)$ and $F_{d}(t)$ are Redner (2001)

[TABLE]

where $D$ is the particle diffusivity. Inserting Eq. 8 into Eq. 7c and taking $s\to 0$,

[TABLE]

[TABLE]

where the sum is over the nodes directly connected to node $j$.

[TABLE]

Expanding for small $s$ gives,

[TABLE]

More generally, $\hat{P}_{ij}$ and $\hat{Q}_{j}$ are

[TABLE]

[TABLE]

where node $m$ and nodes $w$ are directly connected to node $n$.

IV.2 Generating networks

IV.2.1 Synthetic networks

We generate ‘decimated’ networks by constructing a complete 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\times 10^{4}$ honeycomb networks are generated for data in Fig. 2a, $3.6\times 10^{4}$ Voronoi networks for Fig. 2d, and $6.5\times 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.

IV.2.2 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 Viana et al. (ress). The networks we analyze include wild-type cells and $\Delta dnm1\Delta 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\mu$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.

IV.2.3 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 a 6-well, plastic bottom dishes at $1\times 10^{5}$ cells/mL about 18 hours prior to transfection. Plasmid transfections were performed as described previously Hoyer et al. (2018). For all imaging experiments, the ER was fluorescently labeled with 0.2 $\mu\text{g}$ KDEL venus transfected into each well of a 6-well dish English and Voeltz (2013b). 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 $\mu\text{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) Hoyer et al. (2018). Images were taken with a 100$\times$ NA 1.4 oil objective on an electron-multiplying charge-coupled device (CCD) camera 50$\times$50 (Andor). Images were acquired with Micromanager Imaging Software and then analyzed, merged and contrasted using Fiji (ImageJ) Schindelin et al. (2012).

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$\mu$m of one another are the same node. For each cell, the largest connected component was used. We obtained 103 ER networks.

IV.3 Coefficient of variation

This section describes how Fig. 2c was obtained.

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 $\sigma_{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

[TABLE]

where $\left<m\right>_{i}=\frac{1}{n_{i}}\sum_{j}m_{ij}$ is the search time averaged over all bins with a given edge length, $\left<\sigma^{2}\right>_{i}=\frac{1}{n_{i}}\sum_{j}\sigma_{ij}^{2}$ is the average of variance within each bin, and $\text{var}_{i}(m)=\frac{1}{n_{i}}\sum_{j}m_{ij}^{2}-\left<m\right>_{i}^{2}$ is the variance of mean search times across all bins for the given edge length. For each edge length, the averages are done over $n_{i}$ bins containing at least 10 networks.

In the Fig. 2c the blue curve labeled ‘Total edge length and Loop number’ gives the average of the coefficients of variance for each individual bin fixing both edge length and loop number. The average is carried out over all bins corresponding to a particular total edge length:

[TABLE]

IV.4 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 $\pi$ (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\mu\text{m}$ from the sphere surface or the orientation of their normal vector is more than $40^{\circ}$ 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.

IV.5 Approximating search times with percolation theory

We consider a fully connected $n\times n$ square lattice of network nodes on a unit square, giving $N=n^{2}$ total nodes connected to nearest neighbors by edges of length $\ell=1/(n-1)$. The complete lattice has $E_{\text{max}}=2n(n-1)$ total edges. Edges are removed from the network until a number $\Gamma$ of loops remain, without disconnecting any nodes. The number of edges remaining is $E=n^{2}-1+\Gamma$. The fraction of edges that remain is $p=E/E_{\text{max}}$,

[TABLE]

The critical bond probability for percolation on a square lattice is $p_{\text{c}}=1/2$ Stauffer and Aharony (1994), such that our lattice has $p>p_{\text{c}}$. For $p>p_{\text{c}}$ but remaining near $p_{\text{c}}$, the diffusivity depends on the bond probability as $D\sim(p-p_{\text{c}})^{\mu}$, where $\mu\simeq 1.30$ Stauffer and Aharony (1994). The diffusivity on the network is

[TABLE]

Note that for a fixed loop number $\Gamma$, 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=\ell E=n+1+\Gamma/(n-1)$. We can then express the lattice density in terms of our control parameters $L$ and $\Gamma$ according to,

[TABLE]

Assuming a dense lattice system with $L\gg 1$ and $\Gamma\ll L^{2}$, we get the scaling $n\sim L$. This can be plugged into the diffusivity (Eq. 18b) to show

[TABLE]

For a highly disconnected system with low loop number ($\Gamma\ll L$), the diffusivity scales as $D\sim\Gamma^{0}L^{-\mu}$. In the opposite extreme of high loop number ($L\ll\Gamma\ll L^{2}$), the diffusivity scales as $D\sim\Gamma^{\mu}L^{-2\mu}$.

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 Stauffer and Aharony (1994). 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_{\text{step}}\sim N$, neglecting a logarithmic correction term Condamin et al. (2008); Benichou et al. (2010). For our system, the particles diffuse with diffusivity $D$ along edges of length $\ell$, so that the characteristic time to traverse each edge scales as $\ell^{2}/D$. Consequently, the overall expected search time is

[TABLE]

We thus expect the search time to vary as $T\sim L^{\mu}$ for nearly loop-less networks, and $T\sim\Gamma^{-\mu}L^{2\mu}$ 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 $\left<x^{2}\right>\sim t^{\alpha}$ on each edge. Namely, if the dimensionality of the walk ($d_{\text{w}}=2/\alpha$) is less than the dimensionality of the domain ($d_{\text{f}}$) then we have $T_{\text{step}}\sim N^{d_{\text{w}}/d_{\text{f}}}$ Benichou et al. (2010). For a planar network above the percolation transition, $d_{\text{f}}=2$ and Eq. 21 can be modified to

[TABLE]

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 $\mu$) 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.

IV.6 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

[TABLE]

where $N_{g}=(m+2)^{g}+1$ is the number of nodes in the tree Lin et al. (2010). The number of tree edges is given by $K_{g}=(m+2)^{g}$, so $K_{g}\sim 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 $\ell_{g}=2^{-g}$. The total edge length will then be $L_{g}=K_{g}\ell_{g}\sim[(m+2)/2]^{g}$. The time required to diffuse across each edge is $\Delta t_{g}\sim\ell_{g}^{2}$. Overall, the total time for diffusive search to the target will scale as

[TABLE]

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\sim L$, as indicated in the main text.

V Additional information

Competing interests: The authors declare no competing interests.

Author contributions: A.I.B. and E.F.K. designed the study, developed the model, performed the calculations, and revised the manuscript. A.I.B. analyzed the data and wrote the initial manuscript draft. L.M.W. cultured and imaged cells. All authors discussed the results and commented on the manuscript.

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.

Acknowledgements.

Acknowledgments: This work was supported by the Alfred P. Sloan Foundation (EFK), the Hellman Fellows Fund (EFK) and the National Institutes of Health (F32GM116371 to LMW and GM083977 to Gia Voeltz). The authors thank Gia Voeltz (University of Colorado, Boulder) for facilities and Gia Voeltz, Matheus Viana, Susanne Rafelski (Allen Institute), Saurabh Mogre, and Anamika Agrawal (UCSD Physics) for useful discussion and feedback.