Crowded transport within networked representations of complex geometries

Transport in crowded, complex environments occurs across many spatial scales. Geometric restrictions can hinder the motion of individuals and, combined with crowding between individuals, can have drastic effects on global transport phenomena. However, in general, the interplay between crowding and geometry in complex real-life environments is poorly understood. Existing analytical methodologies are not always readily extendable to heterogeneous environments: in these situations predictions of crowded transport behaviour within heterogeneous environments rely on computationally intensive mesh-based approaches. Here, we take a different approach by employing networked representations of complex environments to provide an efficient framework within which the interactions between networked geometry and crowding can be explored. We demonstrate how the framework can be used to: extract detailed information at the level of the whole population or an individual within it; identify the topological features of environments that enable accurate prediction of transport phenomena; and, provide insights into the design of optimal environments.


Introduction
The efficacy of a wide range of cellular processes within living organisms, from protein synthesis [1] to the initiation of a T-cell immune response [2], hinges upon the timely transport of macromolecules through crowded intracellular environments [3]. The motion of individuals within complex environments is central, but not limited to, cell biology, and is important in a wide range of scientific and technological disciplines. Indeed, understanding the roles that both geometry and crowding play in regulating transport processes has immediate and disparate applications across a vast range of spatial and temporal scales, from designing planning algorithms for autonomous robotic motion [4] to utilising the transport of nanoparticles to deliver targeted drug therapies [5]. However, despite the ubiquity of applications, a general means to quantify and characterise the relationships between environmental geometry, crowding and transport phenomena in complex real-life environments remains elusive. Mathematical studies such as this one can play a key role in identifying and quantifying such relationships.
The behaviour of crowded individuals has captured the attention of mathematicians and physicists for decades [6,7,8,9], and the environments in which the individuals are constrained have been found to have large consequences for the emergent transport behaviour. For particles diffusing along an infinite one-dimensional lattice [10], excluded volume interactions between particles significantly hinders the motion of tracer particles, and can induce a qualitative shift in the time-dependence of the mean squared displacement (MSD) from classically diffusive (MSD ∼ t) to anomalously diffusive (MSD ∼ t α , where α < 1). However, for some fractal environments, including diffusion-limited aggregates, the MSD for a tracer particle has the same exponent both in the presence and absence of crowding [11]. Moreover, on comb lattices (one-dimensional backbones with periodic and infinite extrusions representing low-dimensional environments) exclusion interactions along the backbone result in a speed-up of the transport of a tracer particle, and even results in super-diffusive behaviour (MSD ∼ t α , where α > 1) on intermediate timescales [12]. Crowding in complex environments has also been studied on higher-dimensional (greater than or equal to two) regular lattices, such as diffusion in the presence of obstacles [13,14,15] and crowded transport on Manhattan lattices [16] where environmental complexity arises via disordered directionality in the lattice. The above studies make analytical progress in understanding crowded transport behaviour by utilising symmetry and scaling arguments, as well as exploiting the infinite size of these domains. Whilst less common, there are studies that make analytical progress in understanding crowded transport in finite heterogeneous environments. For example the transport of proteins within heterochromatin has been studied using results from transport within fractal environments [17]. Additionally, analytical results from lattice gas cellular automata models have been used to study cell migration in heterogeneous environments, by coupling the automata to a force field representing environmental complexity [18]. However, when interested in exploring microscopic properties of crowded transport in detailed non-fractal real-life environments, such as the intracellular environment, many of these analytical approaches are not readily extendable. In such situations, sophisticated numerical approaches are employed instead.
Significant technological advances in imaging techniques, such as light-sheet microscopy [19] and x-ray tomography [20], have drastically enhanced the availability of high-resolution images of complex environments (see Fig. 1A). Incorporating these detailed geometries into mathematical studies of crowded transport requires discretisation of the geometry into a highresolution mesh, this is a collection of interconnected voxels upon which transport models can be studied numerically. However, numerically integrating equations of motion on these meshes can incur significant computational costs, particularly for stochastic modelling paradigms that require repeated simulations to explore expected transport behaviour [21,22,23]. If these computational costs are too high they present a barrier to studying a large ensemble of geometries, (E-F) Micro-CT scan of the pore space in a Berea sandstone sample, and the corresponding pore network reproduced from [27]. (G) Nanotube vesicle network reproduced from [28]. (H) An electron-microscopy image of three nanotunnels connecting four mitochondria reproduced from [29].
as is necessary to properly characterise how geometry and crowding influences transport behaviour. In addition, it remains unclear how to relate such high-dimensional descriptions of the geometry with statistics of the transport process in order to fully understand the relationship between geometry and crowded transport.
To circumvent the limitations of high-resolution meshes, we require a framework for studying crowded transport in complex environments that avoids traditional spatial discretisation methods. A proposed alternative description of complex geometries is to represent the geometry as a network [24,25,26]. In particular geologists have used these networked representations to study the transport of fluid and sediment through porous media [27]. There are several available algorithms (such as the maximal ball algorithm [25] and SNOW [26]) to extract networked representations that initially dissect complex geometries into several reservoirs (see Fig. 1B).
Connecting these reservoirs are narrow regions of space that are referred to as throats. These throats represent highly restricted regions where individuals experience strong crowding, as seen within heterogeneous porous media [27] (Fig. 1E,F), nanotubes in microfluidic devices [28] (Fig. 1G), and in newly discovered nanotunnels that connect mitochondria [29] (Fig. 1H).
Networks are formed by assigning the reservoirs to be the network nodes, and an edge connects two nodes if the corresponding reservoirs are connected by a throat. The reservoirs are often described as balls with volume equal to the volume of the reservoir, this leads to what is referred to as a balls and sticks network (Fig. 1C). By virtue of the fact that the reservoirs are typically much larger than a voxel in a mesh reconstruction (Fig. 1C), such networks provide a low-dimensional, efficient characterisation of the complex geometry. Indeed, application of a network extraction algorithm (SNOW [26]) to the cardiomyocyte image in Fig. 1A produces a network with 81 reservoirs, whereas the natural Cartesian mesh (where each pixel is represented as a voxel) has 16774 voxels (SI Numerical Methods Section 4.6, Fig. S10).
Beyond dimensionality reduction, networked representations allow us to characterise the role of geometry in regulating transport processes due to the breadth of available topological descriptors that can predict networked transport behaviour [30,31,32]. For example, two classes of summary statistics that are known to accurately predict various emergent transport behaviours are degree-based statistics (these depend solely on the degree distribution of the network, where the degree of a node is equal to the number of edges adjacent to that node) and spectral-based statistics (quantities based on the eigendecomposition of network related matrices, such as the graph Laplacian [33]).
Whilst the majority of the literature on networked transport focusses on single non-crowded individuals, there have been several studies exploring crowded transport processes on networks.
The Totally Asymmetric Simple Exclusion Process (TASEP) is a transport process where particles can move only in one direction along one-dimensional segments and individuals cannot bypass each other. The TASEP has been studied extensively on the one-dimensional line [34], small graphs [35,36] and large networks [37,38,39], where the uni-directionality of the transport results in various behaviours such as shocks and traffic jams [40]. In addition the TASEP has been used to model, for example, protein dynamics along microtubule networks [41,42]. In these models proteins are assumed globally well-mixed within a single large reservoir which is coupled to a filament network upon which the proteins perform asymmetric random walks with exclusion.
In contrast, in this work we are interested in the unbiased transport of individuals that are locally well-mixed within the reservoirs that make up the network. Moreover, the individuals are purely constrained to the network and experience crowding only within geometrically constrained regions (Fig. 1D). As such we build a mathematical framework to study the transport of individuals on general finite networks (in the form of balls and sticks, Fig. 1D) where crowding effects are present only along the edges of the network. We note that there are significantly fewer results concerning symmetric crowded transport on networks compared with the expansive TASEP literature.
We introduce our framework for crowded, networked transport by presenting a hierarchy of diffusive transport models with increasing computational scalability (SI Table 1 provides a summary). The ability of our framework to identify governing principles connecting transport, crowding and geometry is demonstrated through an examination of how crowding and topology affect networked equilibration times [44]. Understanding this relationship is crucial as many results within statistical physics assume an equilibrated population, but the validity of this assumption is rarely addressed [43]. A key result of our work is that heterogeneity in the microscopic structure of complex environments enables low-connectivity networks, as seen in networked descriptions of real-world environments ( Fig. 1A-D), to achieve globally minimal equilibration times. We conclude by extending our framework to provide information on the dynamics of a single motile individual, an extension which opens the door to studying the dynamics of intracellular signalling pathways [45] where the dynamics of individual proteins are known to control a vast range of biological processes, such as T-cell activation and stem cell differentiation.

Results
Individual crowding combined with geometry-induced crowding drastically slows equili- V is the set of reservoirs and their connectivity is specified by the set of edges E. Each edge represents a narrow channel within the geometry where crowding between individuals is nonnegligible (Fig. 1D). To incorporate crowding, a narrow channel (i, j) ∈ E connecting distinct reservoirs i and j is discretised using a one-dimensional lattice with integer length K (i,j) , where individuals undergo a symmetric random walk and at most one individual may occupy each lattice site (SI Extensions Section 3.1 provides a relaxation of this assumption), this process is known as a Symmetric Simple Exclusion Process (SSEP) [46]. Crowded transport within this networked environment is modelled using a canonical framework for stochastic processes, the continuous-time Markov chain (CTMC). A population of N individuals is distributed on the network and their positions evolve as follows. An individual within a narrow channel lattice site attempts to jump into an adjacent lattice site or reservoir at rate α. If the adjacent site is already occupied a collision event occurs and the jump is aborted (Fig. 2D). However, due to the volumetric differences between reservoirs and narrow channels (Fig. 1B) crowding effects in reservoirs are assumed negligible and so jumps into reservoirs are never aborted. Individuals within each reservoir are assumed well-mixed. Those in reservoir i attempt to jump into the first lattice site of one of their connecting narrow channels at rate γ i (Fig. 2D), this means that is the average time taken for an individual in the i-th reservoir to exit the reservoir.
This exit time depends strongly on the local geometry of each reservoir [47,48] and can be calculated as follows. To calculate the mean reservoir exit time for a well-mixed individual we can first consider narrow exit time problems, for which both analytical and computational approaches are readily available [48,49,50]. Narrow exit times quantify the time taken for a particle to reach a narrow opening (such as the opening of a narrow channel) conditional on the initial position of the individual. Integrating these quantities over the reservoir and normalising by the reservoir volume, will yield mean exit times for well-mixed individuals. However, in the interest of maintaining generality we keep the parameters τ i to be abstract, but note particular reservoir geometries can be incorporated as described above. This CTMC is referred to as the full Markov model (FMM) and a technical description of the FMM is found in SI Models Section 1.1.
Evaluating how geometry combined with crowding can impact transport behaviour requires a suitable summary statistic with which to characterise transport. Intuitively, it can be reasoned that the effects of crowding between individuals along narrow channels should result in a "slowing down" of the transport process as individuals in the channels block the paths of others. The insightful statistic to quantify such an effect is the time taken for an initially unequilibrated population of individuals to become well-mixed. This time is known as the equilibration time and is calculated as the reciprocal of the spectral gap (second smallest eigenvalue in absolute value) of the transition matrix for the CTMC [51]. The average reservoir mean exit time, Whilst the FMM is a conceptually ideal way to describe crowded transport in complex geometries, evaluating the equilibration time is computationally infeasible for all but the simplest networks ( Fig. 2A, inset) as the dimensionality of the transition matrix for the FMM is . Therefore, to further investigate the role of crowding and geometry on transport behaviour it is critical to consider dimensionality reduction techniques.
Scaling to large geometries Networked representations of complex environments (Fig. 1A) may contain on the order of hundreds or even thousands of interconnected reservoirs. In light of such, we require models that scale computationally to large environments whilst incorporating details of their microscopic spatial structure. The high dimensionality of the FMM arises from explicitly modelling the occupancy of every lattice site along every narrow channel.
To make progress we introduce a reduced Markov model (RMM) which, in lieu of considering the dynamics within the narrow channels in detail, allows for direct exchange of individuals between reservoirs (Fig. 2C). This approach is reminiscent of the average current calculations for exclusion processes between two open boundaries [52]. For the RMM let n be the configuration vector, where n i is the number of individuals in the i-th reservoir. Focussing on the high-density regime (SI Models Section 1.2.3 discusses the low-density regime), where crowding effects are most important, the rate at which exchange between two connected reservoirs i and j occurs, denoted k HD i,j ( n), is calculated by considering the dynamics of the interacting individuals along the narrow channels (SI Models Section 1.2). By invoking particle-hole duality [53], it is convenient to consider the dynamics of the vacant sites rather than the individuals explicitly. The resulting expression for k HD i,j ( n) is given by Similar to the FMM, the RMM is a CTMC and the reciprocal of the spectral gap of the transition matrix provides the equilibration time. However, the dimensionality of the RMM is still prohibitively high, O N |V|−1 , and does not scale computationally to large complex environments.
Further dimensionality reduction can be achieved through a continuous mean-field approximation of reservoir occupancy in the RMM 2 . Introducing x, such that x i = n i /N HD is the fraction of the population that lies within the i-th reservoir, and expanding the chemical master equations governing the RMM as a Taylor series provides the corresponding Fokker-Planck equation, a partial differential equation describing the evolution of the probability density for the distribution of individuals (SI Models Section 1.3.1, Eq. [S24]). In equilibrium, the networked distribution of the population is determined by the reservoir geometry and is given by . Localising the Fokker-Planck equation about x * reveals that the population dynamics as the networked distribution equilibrates is governed by an Ornstein-Uhlenbeck (OU) process which is known to be Gaussian [54]. Thus, from an initial configuration of individuals x 0 , the configuration at time t follows a multivariate normal distribution with known mean vector µ(t; x 0 ) given by where I is the identity matrix, and covariance matrix Σ(t; x 0 ) given by where F HD and D HD are drift and diffusion matrices, respectively (see SI Models Section 1.3.3 for a derivation). The entries of the drift matrix F HD are given by for i = j, and F HD i,i = − j =i F HD j,i for 1 ≤ i ≤ |V|, where A i,j are the entries of the network adjacency matrix and k HD i,j (N HD x * j ) is the continuous extension of k HD i,j ( n) defined in Eq. (1).
The entries of the diffusion matrix D HD are for i = j, and D HD i,i = − j =i D HD i,j for 1 ≤ i ≤ |V|. Equation (2) reveals that the equilibration time of the OU process is dictated by the spectral gap of the drift matrix F HD which is in the form of a networked graph Laplacian. The spectral gap of F HD accurately predicts the equilibration time for the FMM ( Fig. 2A) and is inexpensive to compute. In particular, the transition matrix for the FMM in Fig. 2A has a dimension of 6968 whilst the weighted graph Laplacian F HD has a dimension of three. We now exploit this remarkable dimensionality reduction to reveal the fundamental principles that govern topological optimisation of equilibration times in crowded environments.
Equilibration times are highly sensitive to an environments network topology To demonstrate how network topology can affect equilibration times we first consider a toy network consisting of five reservoirs (Fig. 3A inset). The three-dimensional reservoir positions x i are uniformly sampled within a unit sphere. The integer narrow channel lengths between reservoirs i and j are given by K i,j = || x i − x j ||/∆ , where ∆ = 2 × 10 −3 is the width of the narrow channels and || x i − x j || is the Euclidean distance between reservoirs i and j. We now consider all 728 possible connected networks with five reservoirs and their corresponding equilibration times.
As a function of topology, equilibration times, calculated from the spectral gap of F HD , vary over several orders of magnitude (Fig. 3A). The topology that induces the quickest equilibration is the complete network (Fig. 3A,B-(xiv)), because the opportunities for individuals to exchange between connected reservoirs are maximised when every connection is present.
However, a complete network is inappropriate to describe most complex environments due to spatial constraints limiting the connectivity of the reservoirs (Fig. 1A-D). The connectivity of a network can be quantified by its total edge length, the sum of the narrow channel lengths present in the network. Imposing a restriction on the total edge length (Fig. 3A, vertical line) reveals a new non-complete optimal network (Fig. 3A,B-(viii)). By varying the restriction over the range of total edge lengths, as defined by the minimum spanning tree(s) and the complete network (Fig. 3A,B-(i) and (xiv), respectively), a frontier of optimal networks arises (Fig. 3A,B).
Networks that lie on the optimal frontier, which we term the optimal networks, represent envi-ronments in which a population equilibrates efficiently, under a given restriction on the environmental connectivity.
A global envelope of optimal networks To explore properties of the optimal frontier (Fig. 3A) we temporarily assume homogeneous reservoir exit times, an assumption that models environments with regular periodic structure such as synthetic porous nanomaterials [55]. Thus, networks that lie on the optimal frontier depend solely on the ensemble K of all possible narrow channel lengths which, in turn, depends upon the spatial configuration of the reservoirs (SI Numerical Methods Section 4.3.1). Comparing optimal frontiers between distinct ensembles Numerical evidence strongly supports the hypothesis that the rescaled optimal frontier follows a global curve that is independent of K and hence is independent of the spatial configuration of the reservoirs (Fig. S6). Over many distinct ensembles K, the variation in the optimal frontier becomes vanishingly small as the number of reservoirs increases (Fig. S6B-D). The global curve persists even for an ensemble of channel lengths that are intentionally sampled from an extremely heterogeneous distribution to ensure that the curve is not a feature of our sampling procedure (Fig. S6E). The apparent persistence of the globally optimal frontier has

Heterogenous non-optimal networks
Heterogenous optimal networks rescaled total edge length 0 1  The range of equilibration factors of both optimal and non-optimal networks for homogeneous, t(0), and heterogeneous, t(1), reservoir geometries. (E) The bar charts (i)-(x) show the weighted degree distribution for both optimal and non-optimal networks across the range of rescaled total edge lengths. The heterogeneous vector of reservoir mean exit times used in (C) and (E) is τ = (2 1 , 2 2 , . . . , 2 10 ). All data presented in (C)-(E) uses the same 5000 configurations of 10 reservoirs with randomly generated ensembles of narrow channel lengths K (SI Numerical Methods Section 4.3.1).
significant implications for optimal network design; testing the optimality of a proposed network merely requires direct comparison of the rescaled total edge length and equilibration factor (two cheap-to-compute network statistics) to the global curve. Moreover, the global curve provides a benchmark to compare the efficacy of algorithms designed to efficiently construct optimal or close to optimal networks (SI Optimal Networks Section 4.2).
Reservoir heterogeneity leads to minimised equilibration times in cases of restricted connectivity Complex geometries can exhibit a range of microscopic spatial structures (Fig. 1B), and this gives rise to reservoir heterogeneity within the corresponding networks (Fig. 1C). The effects of such heterogeneities can be encapsulated by a vector of distinct reservoir exit times τ .
To meaningfully compare the optimal frontiers that arise from two different vectors of reservoir exit times we require that both vectors have the same ensemble average τ . Such a requirement guarantees that the narrow channel equilibrium occupancy is held constant (SI Models Section 1.1.2 Eq. [13]), and thus changes in optimal equilibration times occur solely due to reservoir heterogeneity rather than a change in crowding effects 1 . To systematically explore how heterogeneity impacts the optimal frontier we introduce a vector of reservoir exit times The new reservoir exit times t (φ) depend upon a heterogeneous yet arbitrary vector of mean exit times τ , and the parameter φ controls the extent of reservoir heterogeneity, where t(0) = ( τ , . . . , τ ) represents homogeneous reservoirs, ensuring t (φ) = τ for all values of φ. For increasing levels of reservoir heterogeneity, which is achieved by increasing φ (Fig. 3C, inset), the shape of the optimal frontier changes significantly, and optimal networks achieve globally minimised equilibration times (close to the equilibration time of the complete network) with significantly reduced connectivity (Fig. 3C). Thus, heterogeneity within the internal structure of complex environments has the potential to facilitate globally efficient equilibration in environments with spatially restricted connectivity.
Optimal networks have distinct topological structure The topological structure of optimal networks can be characterised via a weighted degree distribution that arises from considering the diagonal entries of the weighted graph Laplacian, F HD . For the complete network the i-th diagonal entry of F HD is proportional to R i = j =i τ j K (i,j) −1 which represents the total transition rate out of the i-th reservoir when viewing the weighted graph Laplacian as a rate matrix (SI Optimal Networks Section 2.4). As every connection between reservoirs is present in the complete network we term R i the potential transition rate of reservoir i and we order the reservoirs such that R 1 < · · · < R |V| . The ratio of the i-th diagonal entry of F HD for a network G with adjacency matrix A and the i-th diagonal entry of F HD for the complete network is −1 /R i and represents the fraction of the potential transition rate of reservoir i present in the network G. The weighted degree distribution is defined by are translated such that if w i (G; K, τ ) > 0 then the i-th reservoir has a ratio W i (G; K, τ ) greater than the network average, and is referred to as being over-represented in the network. Similarly, a reservoir with w i (G; K, τ ) < 0 is said to be under-represented. Naively, one might expect reservoirs with the highest potential transition rates to be over-represented in any network that lies on an optimal frontier, as connections between reservoirs with high transition rates should encourage faster equilibration. However, we discover that the structure of optimal networks varies greatly depending on the restriction on the rescaled total edge length.
The structure of networks that lie on the optimal frontier compared with non-optimal networks (randomly sampled connected topologies) is distinct across all rescaled total edge lengths ( Fig. 3D), with the greatest contrast seen for networks with mid-range connectivity (Fig. 3D,E (iii)-(viii)). Moreover, the weighted degree distribution highlights how reservoirs can be connected to achieve optimality. Optimal networks with low connectivity and reservoir homogeneity prefer connections between reservoirs with high potential transition rates (Fig. 3E(i) orange/green), on the other hand, for highly connected optimal networks the relative importance of the reservoirs is reversed (Fig. 3E(x) orange/green). Interestingly, as connectivity varies between the two extremes, optimal networks involve over-representation of reservoirs with both high and low potential transition rates (Fig. 3E,(vi)-(viii) orange/green). This transitional behaviour vanishes if reservoir heterogeneity is sufficiently high, where reservoirs with high potential transition rates are over-represented for all levels of connectivity ( Fig. 3E blue/purple).
Even for highly connected networks, the absence of a single important connection between reservoirs can significantly increase equilibration times by over a factor of four in heterogeneous environments (SI Section 2.5, Fig. S8). Collectively our results demonstrate the ability of the weighted degree distribution, as well as the graph Laplacian F HD , to reveal connections between geometric structure and optimal transport that could not have been identified with traditional modelling approaches.
The dynamics of tagged individuals are highly-sensitive to narrow channel length The detailed dynamics of a tagged individual within a population is of interest across a broad range of disciplines [7,56]. For example, the differentiated fate of a stem cell can hinge upon the spatial and temporal dynamics of a single protein within the crowded intracellular environment [57]. A significant benefit of our framework is that it readily extends to provide information at the level of a single tagged individual. The microscopic dynamics of a tagged individual are identical to the dynamics of any individual within the FMM. Therefore, the transition of a tagged individual between adjacent reservoirs i and j via a connecting narrow channel of length K (i,j) occurs as follows. First a tagged individual in reservoir i will jump into the first lattice site adjacent to reservoir i on the narrow channel connecting the i-th and j-th reservoirs. In the high-density regime, the tagged particle will jump back into reservoir i at rate α or jump into the adjacent (second) lattice site when a background individual is exchanged from the i-th to the j-th reservoirs. The latter jump occurs at a significantly lower rate k HD i,j ( n) (SI Eq. [16]) and almost always the tagged individual returns to the i-th reservoir. For this reason we consider a tagged individual to have 'properly' entered the narrow channel only when it reaches the second lattice site (Fig. 4A(i),(ii)). Then, as background individuals continue to be exchanged between reservoirs i and j, the tagged individual undergoes a random walk along the sites of the narrow channel until either being absorbed back into the i-th reservoir (Fig. 4A(i)), or absorbed into the j-th reservoir (Fig. 4A(iii)). The latter occurs with probability p TI x i , x j ; K (i,j) , where the fractions of the population that occupy the i-th and j-th reservoirs at the moment when the tagged individual first enters the narrow channel ( Fig. 4A(ii)) are denoted x i and x j , respectively. The probability p TI x i , x j ; K (i,j) is given by where the functions f k (x i , x j ) are given by The probability p TI is referred to as the tagged individual crossing probability (see SI Models The successful crossing of a tagged individual that has just entered the narrow channel ( Fig. 4A(ii)) requires a net exchange of several background individuals from the i-th to the j-th reservoir (SI Models Section 1.4.1). The probability, p TI , that this net exchange occurs depends on x i / (x i + x j ), the fractional occupancy of the i-th reservoir relative to the j-th is (Fig. 4B, vertical line), as the population equilibrates there is a bias favouring exchange of background individuals from the i-th to the j-th reservoir which subsequently increases the tagged individual crossing probability. However, the probability p TI becomes incredibly sensitive to the fractional occupancy as the length of the narrow channel increases (Fig. 4B, arrow). Indeed, a successful crossing of a tagged individual can require that the fractional occupancy of the i-th reservoir is significantly above the equilibrium occupancy (Fig. 4B, red curve). Thus, in an equilibrated population, the probability that a tagged individual traverses between two adjacent reservoirs is effectively zero if the narrow channel is too long 3 .
The temporal dynamics of tagged individuals are drastically affected by narrow channel length. For an individual that has just entered the second lattice site along the narrow channel ( Fig. 4A(ii)), the tagged individual mean exit time, m TI x i , x j ; K (i,j) , denotes the mean time taken for the tagged individual to exit the channel at either end, and is given by , and the f k (x i , x j ) are as in Eq. (7)  Section 4.5). As such, many useful statistics describing the dynamics of tagged individuals within complex and crowded environments are immediately available. To highlight the utility of the discrete random walk model we consider a first-passage process [58], reminiscent of the Notch signalling pathway where an intracellular protein traverses between the cellular and nuclear membranes [57]. We adopt a caricature representation of the intracellular geometry in the form of a random geometric network (Fig.5A,B and SI Numerical Methods Section 4.3.2), and consider the first-passage properties of a tagged individual initially within a reservoir adjacent to the cellular membrane (Fig. 5A,B, outer circle) whose position evolves according to the networked discrete random walk (SI Numerical Methods Section 4.5) before terminating at a reservoir adjacent to the nucleus (Fig. 5A,B, inner circle).
By comparing our discrete random walk model to an almost identical model that does not consider crowding effects (SI Models Section 1.4.4, Eqs. [S44,S45]), we discover that crowding effects within the narrow channels drastically alters the paths taken by tagged individuals (Fig. 5). The expected number of times an individual traverses a narrow channel becomes highly sensitive to local network topology when crowding effects are incorporated (Fig. 5A,   B). The negligible chance that a crowded tagged individual traverses a long narrow channel ( Fig. 4B, Fig. S6) significantly widens the distribution of the expected number of crossings ( Fig. 5C and A,B highlighted regions), and tagged individuals follow paths that visit reservoirs connected via short narrow channels significantly more often than longer channels (Fig. S4).
Favouring shorter narrow channels subsequently favours indirect paths (Fig. S4A) resulting in an increase in the total path length of an individual that moves from the cellular to the nuclear membrane (Fig. 5D). The alterations of the dynamics of tagged individuals due to crowding, as detailed above, has important implications for the efficiency of signalling pathways and subsequent downstream processes [59].

Discussion
Our results highlight how networks provide an effective framework through which to reveal and quantify the combined influences of geometry and crowding on the transport properties of individuals confined within finite complex environments. Moreover, detailed network analysis uncovered global relationships between population-level transport behaviour and optimal networked topologies, as well as the salient network features responsible for optimisation. Our framework provides new insights into the interplay between geometry, crowding and transport, with such insights being of relevance for a range of geometrically-regulated transport processes.
The work presented in this paper focusses on modelling the temporal fluctuations of particles within the reservoirs of a network. This was achieved through averaging the effects of crowding within narrow channels by assuming that the narrow channel occupancy is held constant.
However, previous work has shown that crowding effects can give rise to strongly fluctuating densities within narrow channels [52]. As such, it would be interesting to extend our frame-work to allow for both the occupancy of reservoirs and narrow channels to fluctuate and explore whether our conclusions about optimal networked topologies change. This would offer a means to test whether the global envelope of optimal networks persists, providing further evidence of a possible universal relationship between rescaled total edge length and optimal equilibration.
Beyond identifying global connections between crowding, geometry and transport, our framework provides a highly efficient computational tool to perform more focussed investigations. For example, cardiomyocyte cells have an intracellular environment consisting of mitochondrial and myofibril filaments. Patients suffering from diabetic cardiomyopathy have been shown to have highly clustered and disordered mitochondrial distributions [60,61]. It has been hypothesised that these alterations to the intracellular geometry help regulate the transport of essential metabolites and increase the overall energy supply to the cell in an effort to maintain regular heart function [60]. A networked modelling approach would offer a sophisticated investigation into the functional role of the observed intracellular restructuring. Analysing the structure of networks extracted from light-sheet microscopy images (SI Fig. S10) would quantify the structural differences between healthy and diabetic cardiomyocytes and studying the transport of metabolites on these networks would provide a means to test whether these changes in intracellular geometry can improve cellular bionergetics. Moreover, the computational efficiency seen within our framework significantly increases the number of experimental images that can be studied compared with traditional approaches that are constrained by using high-resolution meshes. Studying a large number of distinct cellular geometries will be critical to developing a deeper understanding of how intracellular restructuring may regulate cellular bioenergetics.
Our framework offers numerous opportunities for generalisation. The RMM is formulated as a chemical reaction network [62] and thus can immediately support reactions between individuals, where subsequent coarse-graining will result in a Fokker-Planck equation capable of investigating geometry-controlled kinetics [63]. Included within the SI is a generalisation of our framework that supports active transport (SI Extensions Section 3.2). We study a Partially Asymmetric Simple Exclusion Process (PASEP) that allows individuals to undergo bidirectional motion along narrow channels but with a bias in one direction. Extending to active transport opens the door to investigating the geometric influences on highly crowded transport phenomena across a wide array of spatial scales, from mRNA translocation along intracellular microtubule networks [64], to molecular trafficking between cells connected via cytonemes [65] or plasmodesmata [66], to the transport of sediment subjected to flows within porous media [67].
Through reconceptualising how we model crowded, geometrically-constrained transport, we stand to gain significant new insights within fields such as optimal synthetic design [68] and molecular cell biology [69], to name a few. Furthermore, this work presents a versatile framework that paves the way for furthering our understanding of the fundamental connections between crowding, networked geometry and transport, beyond what has already been achieved under more traditional modelling paradigms.