Fluid mechanics of luminal transport in actively contracting endoplasmic reticulum

The Endoplasmic Reticulum (ER), the largest cellular compartment, harbours the machinery for the biogenesis of secretory proteins and lipids, calcium storage/mobilisation, and detoxification. It is shaped as layered membranous sheets interconnected with a network of tubules extending throughout the cell. Understanding the influence of the ER morphology dynamics on molecular transport may offer clues to rationalising neuro-pathologies caused by ER morphogen mutations. It remains unclear, however, how the ER facilitates its intra-luminal mobility and homogenises its content. It has been recently proposed that intra-luminal transport may be enabled by active contractions of ER tubules. To surmount the barriers to empirical studies of the minuscule spatial and temporal scales relevant to ER nanofluidics, here we exploit the principles of viscous fluid dynamics to generate a theoretical physical model emulating in silico the content motion in actively contracting nanoscopic tubular networks. The computational model reveals the luminal particle speeds, and their impact in facilitating active transport, of the active contractile behaviour of the different ER components along various time-space parameters. The results of the model indicate that reproducing transport with velocities similar to those reported experimentally in single particle tracking would require unrealistically high values of tubule contraction site length and rate. Considering further nanofluidic scenarios, we show that width contractions of the ER’s flat domains (perinuclear sheets) generate local flows with only a short-range effect on luminal transport. Only contractions of peripheral sheets can reproduce experimental measurements, provided they are able to contract fast enough.


I. INTRODUCTION
The mammalian endoplasmic reticulum (ER) is the single largest intracellular structure (see sketch in Fig. 1a). The organelle is made up of membranous sheets interconnected with the nuclear envelope and branching out into a planar network of tubules extending throughout the cell periphery (Fig. 1b) [1]. The ER dynamics on a second scale include the cytoskeleton-assisted tubular network restructuring [2,3] and an interconversion between two distinct forms, including a narrower form covered by membrane curvature-promoting proteins [4]. The ER morphology and its dynamics presumably enable and facilitate its functions: the ER is responsible for the production, maturation and quality-controlled folding of secretory and membrane proteins, which constitute approximately a third of the cell's proteome [5]. The organelle's membranes also harbour the lipid biosynthesis machinery, while its lumen stores calcium. The contiguous nature of the ER is believed to ensure an efficient delivery of all these components across the cell periphery.
Timely transport of the content within the ER is therefore integral to the function of the cell. The geometry and dimensions of several cell types with extensive ER-containing projections (e.g. neurons and strocytes) pose a kinetic challenge for material distribution with physiological timing. These considerations predict the need for an active luminal transport to ensure timely material homogenisation across the vast ER. Empirically, the active nature of the ER luminal transport is supported by the sensitivity of GFP mobility (measured by FRAP) to ATP depletion [6,7]. Thus past measurements indicate that the transport of proteins in ER is not consistent with Brownian motion [6,8]. Using single particle tracking of ER luminal markers, an active solute transport in the planar ER network was observed directly [7]. A photoactivation chase technique also measured a superdiffusive behaviour of luminal material spread through the ER network [9].
However, the mechanism for generating ER luminal flows remains unclear. Understanding the mode of material exchange across the organelle is crucial for rationalising the ER shaping defects-related neuronal pathologies [10], identifying factors controlling ER transport and informing the development of ER transport modulation approaches with health benefits. Based on ER marker velocity fluctuations measured in single particle tracking and the detection of transient narrowing points in the tubules by improved super-resolution and electron microscopy [7], it has been postulated that these active flows may result from the stochastic contractility (pinching and un-pinching) of ER tubules at specific locations along their lengths (see sketch in Fig. 1e); other plausible mechanisms for flow generation were also considered. However, measurements for testing this pinching hypothesis are currently inaccessible, due to limitations in space-time resolution of live cell microscopy, i.e. the live cell-compatible super-resolution techniques achieve resolution of ∼80 nm at the relevant speed, while the tubular radius is estimated in the range 30−60 nm [9,11,12]. Improvement or resolution currently can only be achieved trading off speed. To circumvent these experimental difficulties, in the current study we use mathematical modelling to quantitatively analyse the relevant scenarios of actively contractility-driven flows. Specifically, we explore how various sets of spatio-temporal parameters of ER contractility may produce flows facilitating solute transport in quantitative agreement with experimental measurements.
Specifically, assuming the flow is driven by tubule contractions (shown schematically in Fig. 1e), we construct a physical and mathematical model of the ER network and solve it for the flows inside the network (description in Methods §IV). We use our model to carry out numerical simulations to study the motion of Brownian particles carried by these flows and show that the tubule pinching hypothesis is not supported by the results of our model ( §II A), a result independent of the network geometry ( §II B). The failure of active pinching to drive strong flows can be rationalised theoretically by deriving a rigorous upper-bound on the rate of transport induced by a single pinch ( §II C), including possible coordination mechanisms ( §II D). Only by increasing both the length of pinches and their rate can beyond admissible values can we produce particle speeds in agreement with measurements ( §II E).
We then explore in §II F two different hypotheses as possible explanations for the active ER flows, first the pinching of the junctions between tubules (Fig. 1c) and then the pinching of the two types of ER sheets, perinuclear (Fig. 1f) and peripheral (Fig. 1d). We investigate the conditions under which both would be consistent with the experimental measurements of active transport in Ref. [7]. Our results suggest that the biological origin of solute transport in ER networks remains open and call for extensive empiricalexploration of the alternative mechanisms for flow generation.

II. RESULTS
A. Tubule contractility-driven ER luminal motion yields inadequate transport kinetics To assess the kinetics of particle motion in the lumen of tubular structures in response to their contractility, we generated an in-silico simulation model of the process. The model incorporates local calculations for the low-Reynolds number hydrodynamics of a contracting tubule into a global analysis using Kirchhoff's Laws and standard graph theoretical results of the flows throughout the network geometry (see Methods §IV for details). This allows us to compute the dynamics of particles inside structures such as that the ER accounting of the effects of their dynamic morphology.
We initially implemented the model's numerical simulations of particle transport in a reconstructed ER network of a COS-7 cell [7] (which we label C0, see §IV A) with tubules locally contracting stochastically according to the spatiotemporal parameters suggested by microscopy measurements [7]. We will therein refer to these contractions, illustrated in with an average flow speed of 1.3 µm/s (see also Supplementary Video S1). Further, the resulting instantaneous speed distribution of Brownian particles advected by these flows (methodology in §IV E) is considerably shifted relative to the experimental counterpart [7] towards lower values (Fig. 2b). Similarly, the distribution of average edge traversal speeds (defined precisely in §IV F) from our simulations (solid blue line in Fig. 2c) is lower than the experimentally measured speed distribution by an order of magnitude (Fig. 2c, inset).
Moreover, the results for all measures of transport under pure diffusion, in the absence of pinching-induced flows, were virtually indistinguishable from those where the pinchingdriven flows are included (see Fig. 2b-c). Within the framework of transport theory, this conclusion can be rationalised by estimating the relevant Péclet numbers that measure the relative importance of advection and diffusion. Using the average valueŪ ∼ 1.3 µm/s of the mean flow speeds over time and edges as a velocity scale, we may estimate a mean Péclet number as Pe ∼Ū R/D. Using R = 30 nm and the measured diffusivity D ≈ 0.6 µm 2 s −1 [7], this leads to the estimate Pe ∼ 0.07. Flows affect transport for Péclet numbers above order-one values [13], and therefore the pinching-driven flows have a negligible influence on the transport inside the ER network. In order for fluid motion to have a noticeable effect, one would thus need flows to be either stronger, or sustained in one direction for longer.
The chaotic flows produced by the pinching events with stochastic parameters suggested by resolution-limited microscopy [7] appear too weak to generate fast, non-diffusive edge traversals.
B. Tubule pinching-induced transport is network geometry independent and fails to facilitate luminal homogenisation To estimate the network geometry's contribution to transport simulation outcomes, we compared results across four different ER structures (C1-C4, Fig. S8). The distributions of the average edge traversal speeds appeared insensitive to ER structure variations for both pinching-induced and exclusively diffusional transport. This is reflected in the small deviations from the mean of the data points averaged across the different structures ( Fig. 2de). Further, pinching-induced flows inside a regular honeycomb network (Fig. 2f, inset) with a typical ER edge length (1 µm) appears within a reasonable variance compared to the natural networks. Therefore, this mathematical idealisation of the ER network geometry can be used for exploring the consequences of the network ultrastructure contractility on transport kinetics in a standardised manner.
To test the effectiveness of the particle velocities for facilitating luminal material exchange across the ER, we tracked homogenisation kinetics by measuring intermixing of particles of two distinct colours equally seeded in each half of a honeycomb network at t = 0 ( Fig. 3; see also Supplementary Video S2 for the flows inside such a network driven by pinches with the parameters from [7]). The measure of homogenisation is given by the variance Var(φ(t)) ≡ Var(n b (t) − n r (t)) over twenty regions of the network (Fig. 3a, horizontal lines) of the difference between the numbers n b (blue) and n r (red) of particles of each colour in region; note Var = 0 represents perfect homogeneity. The measures of mixing over time for pure diffusion and active transport with parameters [7] from Fig. 2 show a nearly complete overlap ( Fig. 3e; see also Supplementary Videos S3 and S4), distinct from faster mixing under stronger flows in a network driven by pinches whose lengths are increased to their maximum possible value (i.e. the length of the tubule) and which are in addition sped up by a factor of 10 (see Supplementary Video S5; see also §II E for a discussion of the effects of artificially strong pinch parameters on average edge traversal speeds). This indicates that even if the experimental particle speeds were overestimated, the presumed pinching parameters would be inadequate to facilitate ER luminal material exchange.
C. Theoretical analysis of advection due to a single pinch explains weak pinchinginduced transport The slow luminal transport driven by the pinching-induced flows is intrinsically linked to the volume of fluid expelled by each pinch during a contraction. The fundamental reason underlying the weak pinching-induced transport is that individual pinches are very weak generator of flow; even in the best possible configuration, the volumes of fluid pushed by each pinch are too small to impact luminal transport. Specifically, in Methods ( §IV H 1), we mathematically show that, given a pinch of length 2L, the maximum distance ∆z max a suspended particle may be advected instantaneously by an individual pinch is Using the experimentally-measured average value of the pinch length 2L = 0.14 µm [7], a typical pinch can then propel a particle by a maximum distance of ∆z max ≈ 0.19 µm.
This may, equivalently, be framed in terms of velocities. A transported particle experiences an average velocity during the contraction of at most V max = 8L/3T , where T is the duration of a contraction or a relaxation. Using the pinch length as above and the average values of 2T = 0.213 s [7] this leads to V max ≈ 3.9 µm/s, which is an order of magnitude smaller than the measured typical edge traversal speed of ∼ 45 µm/s, consistent with the order of magnitude difference between the measurements and the predictions of the computational model.

D. Marginal addition of coordinated contractility to luminal transport
The theoretical upper bound in the previous section for the maximum luminal transport produceable by an individual pinch is realisable only in the hypothetical situation where the flow generated by the pinch is all directed to one end of the tubule i.e. when the other end is effectively blocked. However, content transport produced during the contraction of a single tubule would then be reversed when the tubule relaxes back to its original state, with a single pinch site expected to exhibit only reciprocal (i.e. time-reversible) motions. Any advection contributing to edge traversals must thus be dominated by non-reciprocal motions of multiple pinches resulting in net displacements of solute particles.
The simplest system capable of producing non-reciprocal motions consists of two pinches, and the optimal sequence of motions to maximise the resulting advective particle displacement is illustrated in Fig. 5. We show in Methods ( §IV H 2) that this is indeed the optimal two-pinch coordination, which results in a time-averaged displacement equal to the upper bound derived in Eq. (1). Since this optimal sequence of motions involves one pinch site starting a pinch halfway through the pinching of the other site, it is reasonable to estimate its duration as 3T , and therefore an average particle speed of 8L/9T ≈ 1.3 µm/s. The low particle speed achievable by the optimal coordination between two pinches suggests that a very high level of coordination among multiple tens of pinches per tubule would be required to reproduce the measured edge traversals.
E. A combination of high frequency and pinch length is required to replicate experimental particle speeds Since the magnitude of ER contractility suggested by microscopy [7] (pinch lengths and frequency) do not explain the measured speeds, we set out to explore different sets of parameters that may generate particle velocities closer to the experimental measurements. The currently achievable imaging spatiotemporal resolution limits the detection of tubular diameter contractility by microscopy. Therefore, it is reasonable to postulate that the relevant parameters may have been underestimated. We simulated ER transport varying individually or in combination the values of pinch duration 2T , waiting time T wait between successive pinches on a tubule, and pinch length 2L.
We first decreased both the original [7] values of T and T wait by the same factor of 1/α, where α ≥ 1, and simulated particle transport in the honeycomb network (Fig. 4a).
In effect, this simply 'fast-forwards' the flows in the original active network by a factor of α, and Brownian particles of the original diffusivity are released into this sped-up flow.
Instantaneous and edge traversal speeds exhibited corresponding increases when we increased the value of α (Fig. 4a). An extreme value of α = 100 produces an average edge traversal speed distribution that peaks at around 8 µm/s (Fig. 4a). Similar results are observed in the C0 network from a COS-7 cell (Fig. 4b). The longer tails of these distributions Next, we attempted to obtain a better fit to experimental data by combining changes in both the pinches' time and geometry parameters. In Fig. 4e (labelled 'pinch L max , 10× rate'), we thus sped up these maximally long pinches by a factor of α = 10, which yielded speeds averaging above 60 µm/s. Notably, the tail of the speeds distribution appeared longer than that seen in the experiments. The sensitivity of particle speed to the changes in contraction site size and α is illustrated by the shifts in the speeds distribution toward lowers values upon reduction of these parameters ( Fig. 4f-i).
F. Luminal transport kinetics derived from contractile ER tubular junctions and sheets.
As shown above, establishing effective transport in a tubular constrictions-driven model based on realistic ER fluid dynamics required a set of questionable assumptions, compelling us to explore alternatives. Thus, we set out to investigate how ER luminal transport would be impacted by the contractility of its structural components with volumes larger than tubules: (i) the tubular junctions ( Fig. 1c), (ii) the perinuclear ER sheets (Fig. 1f) and (iii) the peripheral sheets (Fig. 1d).
First, we run numerical simulations of transport driven by contracting junctions on an ER network (C1 from Fig. S8d; sketch of junctions in Fig. 1c). Since junctions contracting at the tubular pinch temporal parameters measured experimentally yielded inadequate transport, Next, we considered the particle transport generated by two types of ER sheets contracting and relaxing over a duration 2T (see Methods for details). The perinuclear sheets are shaped as contiguous layers of flat cisternae with a luminal volume larger than the tubules branching from these structures (see sketch in Fig. 1f). Accordingly, their contraction with 2T = 5 s and V sheet = 10 µm 3 yielded a mean average traversal speed of 35 µm/s and consistent with the single particle tracking experiments (see Fig. 7b). However, the speeds distributions tailed towards higher values ( Fig. 7a-b), something that has not be observed experimentally so far; it may be that high velocities cannot be recovered in experiments due to the constraints on linkage distance imposed to ensure trajectory fidelity in particle tracking.
Further, our simulations reveal that the flow decays sharply with distance away from the sheet where it originated as it branches out into tubules. This is illustrated in Fig. 7c-d and we further quantify the spatial gradient of the average edge traversal speeds across Cartesian 2D coordinates in Fig. 7e, revealing the stark contrast between the homogeneous profile for contracting junctions (red dotted line) vs sheet-driven transport (blue solid line). The short range of influence afforded by contracting perinuclear ER sheets thus argues against its ability to sustain mixing flows and fast particle transport on the distal tubular network.
Instead, we sought to explore whether the the peripheral sheets (i.e. the smaller flat intertubular ER regions, see sketch in Fig. 1d) may overcome the range limit. The peripheral sheets have volume significantly larger than tubules and junctions, which we estimated at 0.12µm 3 ± 0.04 µm 3 (see Fig. S13 and Methods §IV G 4 for details). This suggests that the peripheral sheets could produce average edge traversal speeds compatible with experimental measurements. Assuming that a sheet typically occupies the area enclosed by a 'triangle' of tubules, we may incorporate the transport inside a sheet-driven network using our model for junctions-contraction, but with junction volumes set to the measured sheet volumes (see §IV G 4 for details) and with either (i) each node expelling one third of the volume expelled by a contracting sheet (since peripheral sheet are in contact with three nodes on average), or (ii) a contracting node expelling the entire volume expelled by a contracting sheet, but with only one-third of the junctions actively contracting at any one time. Further assuming that during each contraction an peripheral sheet expels half its total volume so that in simulation (i) node k expels a volume V k /6 in each contraction and in simulation (ii) each active node k expels a volume V k /3, nodes which contract at rates α −1 = 2.5 and α −1 = 5 times slower than the tubule pinches in Ref. [7] in simulations (i) and (ii) respectively give average edge traversal speeds in the correct range of 40 µm/s ( Fig. 7f-g). The contraction of peripheral sheets may thus be offered as a plausible mechanism for fast luminal transport, provided they are able to contract with sufficiently high rates.

III. DISCUSSION
The motion of solutes in cellular compartments is now understood to be facilitated by active components. This is evident from direct motion measurements and the dependence of motion speed on the availability of ATP-contained energy [6][7][8]14]. The origin of these active driving forces is however challenging to identify. The cytoplasm's currents are often believed to originate from the motion of large complexes such as ribosomes and large vesicle cargo shuttled by the cytoskeleton-mediated motorised transport [14]. In the case of enhanced transport in the lumen of the ER tubules, the contractility of tubules has been suggested as the flow generating mechanism, and indeed such tubule deformations have been observed in microscopy [7]. However, establishing the direct empirical link between tubule contraction and active flows, or experimentally testing other hypotheses for the driving mechanism behind ER solute transport, remain currently unattainable. In this study, we thus propose a physical modelling approach, which provides a platform to explore the nanofluidics behaviour of biological systems such as the ER network. The outcomes of our simulations for a contractile ER argue against the plausibility of local pinch-driven flow; pinches with frequency and size on the order of those estimated by microscopy yield significantly lower speeds than single particle tracking measurements, as well as no enhancement of mixing beyond that from passive diffusion. The deficit stems from the fact that the displaced fluid volume upon local contraction is too small to generate sufficient particle transport.
Given the uncertainty of the empirical measurements for ER tubule deformation, due to the limits in the spatio-temporal resolution of organelle structure imaging, the pinch parameters may have been significantly underestimated. This sanctioned exploring a wider range of spatio-temporal parameters in our simulations, which revealed that a combination of a higher frequency with a much larger pinch length may provide higher particle speeds that are comparable to the single particle tracking measurements. Further, our modelling results suggest the possibility of transport by luminal width contractions of larger volume ER subdomains (which are contiguously interconnected with the network). In that respect, the contraction of peripheral (i.e. inter-tubular) sheets, in particular yielded speed values in a plausible range provided they are able to contrast fast enough. In contrast, the contraction of a large-volume perinuclear sheets led to fast transport but with a limited spatial range and with flows that would not impact transport beyond a few microns into the peripheral network. Similarly, the alternative scenario of tubular junctions' contractility appears unlikely as particle speeds similar to experiments could only appears for junctions pinching at up to 8 times faster than experimentally observed contracting tubules.
The modelling approach in this study, although focused on the ER network, provide a step forward toward understanding intra-organellar fluid dynamics. Our simulation results have been used to rule out several scenarios, which seemed physically intuitive but are nevertheless unable to explain the observed enhanced luminal transport. Furthermore, our results generated a set of potentially testable predictions that can be used to validate or refute each of the envisaged transport mechanisms. For example, as fluid expulsion from peripheral sheets appears to be in broad in agreement with single particle tracking measurements data under even conservative assumptions, future measurements may explore whether an active luminal motion is more pronounced and faster in proximity to peripheral sheets. Moreover, a set of improved spatio-temporal resolution measurements of particle tracking and structural contractions will be needed to complete the physical picture of ER luminal transport.
In conclusion, it is worth emphasising that our in-silico fluid dynamical modelling reveals that for structural fluctuation-based mechanisms to facilitate luminal motion require assumptions currently not supported by empirical data. This warrants explorations of alternative mechanisms for ER luminal transport, for example anomalous diffusion driven by the fluctuations of macromolecular complexes [14] or by osmotic forces, as previously suggested [7].

A. Network modelling
We represent the 'skeleton' of a two-dimensional ER network as a planar graph with each node assigned a position x ∈ R 2 . Given an edge of the network labelled (i, j) and of length |x i − x j | = l, we model the lumen of the corresponding tubule to occupy a cylinder of radius R whose axis lies along the edge and has length l. This assumption avoids the intrinsic difficulty in defining a precise boundary between a tubule and a tubular junction, as well as leading to a simplified model of the intra-nodal dynamics of a solute particle (see below); since the tubules are long compared to the size of the junctions, the impact of their small overlap can be safely neglected.
We mathematically reconstruct a model ER network, which we refer to as C0, using the skeleton image of the COS7 ER network given in Fig. 4

of the Supplementary Material from
Ref. [7]; this network is reproduced in Fig. 8a. We use the multi-point tool in ImageJ [15] to place numbered points at the positions of the nodes on the source skeleton image, and then obtain a list of the indices and position coordinates of each node. The edges are then manually tabulated as a list of pairs of nodes; this gives us all the information required to construct a mathematical graph as shown in Fig. 8b, with the original network superimposed on the mathematical model in Fig. 8c. Note that in what follows we work with the largest connected component of the source skeleton image in order to study transport in a fully connected network. We also use the same procedure to extract the graph structures from microscopy images, of four smaller ER networks which we label C1-C4 (original network with mathematical graph superimposed in Fig. 8d). In Supplementary Fig. 10, we show the distributions of the edge (tubule) lengths in each of the C0-C4 networks, as well as the mean edge lengths and the mean degrees (the degree of a node is the number of edges connected to it). The mean edge lengths are around 1 µm and mean degrees are approximately 3. In order to compare the biological network to an idealised ER system, we also consider below a honeycomb network, i.e. one where every node (apart from those at the boundaries) has a degree of 3, with all edge lengths exactly 1 µm.

B. Pinch modelling
To describe the pinches (Fig. 1e), we consider a model where the kinematics of each pinch is fully prescribed in time. We therefore assume that the active biochemical forces responsible for the deformation of the tubules balance with the elastic resistance of the tubules and with the dissipative forces in the fluid in such a way that the pinches occur as described.
The geometrical model for a pinch is illustrated in Fig. 9. Each tubule is assumed to have a pinch site at its midpoint, with each pinching event happening stochastically and independently of each other. Each pinch is defined geometrically by three parameters from a random distribution (see below): (i) the duration of a pinch 2T , (ii) the time T wait between the end of a pinch and the beginning of a new one on the same site, and (iii) the length of a pinch 2L (see Fig. 9). We assume for simplicity that all pinches are axisymmetric so, using the notation of Fig. 9, each tubule remains a cylinder of time-varying radius denoted by r = a(z, t). Pinches are assumed to have reflectional symmetry about the plane in the crosssection through the centre of the pinch (i.e. Fig. 9 each pinch is characterised by the same length L on either side of it). The total length of a tubule is denoted by l = L 1 + 2L + L 2 , so the portion before the pinch has length L 1 , and that after the pinch is of length L 2 .
We model the geometrical profile of each pinch of length 2L as following a linear radius change (see Fig. 9). Within a pinch located at z = z 0 , the radius of the cylinder is therefore given by The time-varying function b(t) ≤ R is therefore the minimum pinch radius in the centre of In our simulations the stochastic pinch parameters are drawn from the distributions measured in Ref. [7]. The pinch duration 2T is therefore drawn from an exponential distribution with rate parameter λ = ln 10/0.167 s −1 (so the mean value is 2T = 0.0725 s) , the time between pinches T wait from from an exponential distribution with rate parameter λ wait = ln 10/0.851 s −1 (mean value T wait = 0.370 s) and the pinch length 2L from a uniform distribution with mean µ = 0.12 µm and variance σ 2 = (0.06 2 / ln 10) µm 2 ≈ (0.040 µm) 2 .
Throughout we set the tubule radius R = 30 nm and b 0 = 0.01R (recall from Fig. 9).
C. Hydrodynamic modelling

Hydrodynamics of network
We assume that the fluid occupying the ER network is Newtonian. Flow inside the network occurs at low Reynolds number, which can be justified in the following fashion.
From experiments in Ref. [7], we know that velocity scale relevant to ER flows is of the order U ∼ 10 −5 ms −1 . With a typical tubule radius R ≈ 10 −8 m and a kinematic viscosity of at least that of water ν ≈ 10 −6 m 2 s −1 , we obtain a Reynolds number of the order Re = U R/ν ∼ 10 −7 , so the flow is indeed Stokesian and inertial effects in the fluid can be safely neglected.
At a typical instance in time, a tubule undergoing contraction (or relaxation) causes a net volume flux to exit (or enter) the corresponding tubule. In the context of our graph theoretical model, we therefore model each pinch site as a 'pinch node' generating a net hydrodynamic source/sink whenever the tubule contracts/relaxes.
However, since it is not guaranteed that the total volume created by all pinching events always adds up to zero, we need a mechanism for the corresponding net volume to exit, or enter, the network. We achieve this through a number of 'exit nodes' that allow mass to be globally conserved. We locate the exit nodes at the exterior of the network and choose their number randomly. From a hydrodynamic standpoint, we impose the pressure condition p = p 0 at each exit node to model their connection to a large fluid reservoir.
We numerically tested the robustness of our results to the details of the exit nodes by where µ is the dynamic viscosity of the Newtonian fluid. Note that when Q > 0 our notation leads to ∆p < 0, meaning that the pressure decreases across the length of the channel.
The result in Eq. (3) is valid for a straight (i.e. not pinching) tubule, and we need to generalise it to the case of a pinching tubule. Consider first a more general axisymmetric pipe whose radius a(z, t) varies with axial position z and time. We may use the longwavelength (lubrication) solution to Stokes' equations for the streamwise velocity u(z, r, t) and flux Q(z, t) inside such a pipe into which a flux Q 1 (t) enters at z = 0 [17], where The equality in Eq. (5) can be derived using an intuitive mass conservation argument, independently of the inspired ansatz in Eq. (4). Conservation of mass inside a small section [z, z + δz] of the cylinder requires Q(z) − Q(z + δz) = 2πaȧδz + O(δz 2 ). Considering the limit δz → 0 and integrating the resulting expression for ∂Q/∂z from 0 to z yields Eq. (5).
In order for the no-slip boundary conditions at r = a(z, t) to be satisfied, and also to satisfy the incompressibility condition ∇ · u = 0, the radial component of the velocity, v(z, r, t), is necessarily non-zero and given by [17] v(z, r, t) = ∂a(z, t) ∂t r a(z, t) 2 − r a(z, t)  Fig. 9, we may substitute the linear shape functions into Eq. (5) and obtain Note that the final expression may be rearranged as which may be interpreted as the instantaneous volume source/sink during a contraction/relaxation at a pinch site.
b. Pressure drop. We next need to compute the pressure drop in the pinches. We integrate the z-component of the Stokes equation along 0 ≤ z ≤ L 1 + 2L + L 2 , and use the solution for u, to obtain (12b) Here, the second term on the right-hand side of Eq. (12a) has vanished because ∂ z u ∝ ∂ z (Q/a 2 ) but Q and a are approximately constant at the entrance and exits of the tubule when the pinch site is sufficiently far from the ends of the tubule so that the flow is fullydeveloped there.
Using the expression for Q 2 , an integration yields and, by symmetry, Subtracting these two results (and noting the integrals cancel out by symmetry) we obtain the modified Hagen-Poiseuille expression for a pinching tubule as Note that this relationship is linear in each of Q 1 and Q 4 , and is to be solved alongside Eq. (10) to relate the flow rates and pressure drops to the change in size of the pinches.
Importantly, the classical Hagen-Poiseuille law is recovered as b → R since Eq. (15) becomes in that limit which agrees with Eq. (3) when taking l = L 1 + 2L + L 2 .

D. Solving the hydrodynamic network model
The incorporation of pinches as 'dummy nodes' into a graph theoretical framework of §IV A along with Eqs. (10) and (15) for the necessary pinch-related quantities allow us to reduce the problem of determining the time-dependent flows in an active pinching network into the simpler problem of solving at each instant for the instantaneous fluxes inside a 'passive' network with newly added nodes, appropriate sources/sinks, and modified pressure drops. Note that since the flows at these sub-cellular scales are inertialess (i.e. Stokes flows), we are able to effectively decouple time from our problem and solve the problem in the quasi-steady limit.
For each edge (i.e. tubule) (i, j) in the network, we define Q ij to be the flow rate from node i to node j, with the sign convention such that flow is from i to j if Q ij > 0. For mathematical convenience, we define Q ij = 0 in all cases where (i, j) is not an edge in the graph. The goal is to solve for the values of the Q ij 's corresponding to each edge.
After the incorporation of the dummy nodes, we denote by N the number of nodes and E the number of edges. We label the nodes such that {1, . . . , M } denotes the M exit nodes.
Let q i be the source or sink carried by the i th node (so that q i = 0 if i is a normal node, q i is as specified by the RHS of Eq. (10) if i is a pinch node, and q i is a quantity to be determined if i is an exit node.). Our M + E independent variables are therefore {q i |i = 1, . . . , M } i.e. the sources/sinks carried by the exit nodes, and the Q ij 's corresponding to each edge. To obtain their values, we employ the viscous hydraulic analogues of Kirchoff's Laws.

Kirchhoff 's First Law (K1)
The first equation is that mass is conserved at each junction i.e. for each node i we have which gives us therefore N equations. Note that these equations together imply global conservation of mass, N i=1 q i = 0.

Kirchhoff 's Second Law (K2)
The second equation is a consistency of pressure, namely that the pressure change around any cycle (i.e. closed loop) of the network is zero. Therefore, in a given cycle Note that the pressure change across a node is negligible.
The K2 statement in Eq. (18) applies to all cycles in the graph, which would give us more equations than we need since the vectors of coefficients in Eq. (18) are linearly dependent.
Instead, we need a minimal set of linearly independent K2 equations, corresponding to the cycles in a cycle basis of the graph, and we need only apply K2 to these cycles.
To construct a cycle basis of the graph G, we use standard results from graph theory [18].
We first construct a spanning tree T , defined as a connected subgraph which contains all the nodes of G and no cycle, as shown in Supplementary Fig. 11a on an example. Any tree T has N − 1 edges. Therefore there are E − (N − 1) edges in the graph G but not in the spanning tree T ; for each such edge e, we denote by C e the unique cycle in the graph created by adding the edge e to the tree T (see Fig. 11b). The set of all such cycles C e is then a cycle basis of G, and thus there are E − (N − 1) cycles in this set. There are therefore E − N + 1 independent cycles in the cycle basis.
To compute the spanning tree T we use a breadth-first search (BFS) algorithm [19]. We start from an arbitrary node and explore its neighbouring nodes. We add to T any previously unexplored node (and the corresponding edge) which does not result in the creation of a cycle in T . We then repeat this (in an arbitrary order) on the neighbours of the previous generation of nodes added to T , until no more nodes are left to explore. This is illustrated on an example in Supplementary Fig. 11c.
A similar algorithm is also used to compute the cycles C e in the cycle basis. This time, however, denoting e = (i, j), the algorithm is started from i and set to terminate as soon as j is visited, yielding a path in the tree from i to j, which, together with the original edge (i, j), completes a cycle.

Pressure boundary conditions at exit nodes
Here the sum i→j is defined to be over any path P ij from i to j. Thanks to Kirchhoff's second law, this quantity is path-independent. This indeed gives us by the flow inside each tubule. Let x(t) denote the position of a Brownian particle in a tubule and x n the finite-difference approximation of x(n∆t), where ∆t is a discrete time step. The displacement of the particle at each time step can be obtained approximately using an explicit first-order Euler scheme where U is the instantaneous flow velocity, and the random noise term X(∆t) is drawn from a zero-mean Gaussian with variance X(∆t)X(∆t) = 2D∆tI [20], where D is the Brownian diffusivity of the particle. In our simulations, we take the diffusion constant to be the mean intranode diffusivity measured in Ref. [7], D ≈ 0.6 µm 2 s −1 . We include interactions between particles and walls by assuming that particles perfectly reflect off walls (i.e. elastic collisions).
The particles are modelled as rigid spheres of diameter 5 nm, and the size of the particle matters only during elastic collisions with the walls of the tubules. As relevant in the limit of low volume fraction, we neglect hydrodynamic interactions between particles and perform ensemble averaging of the trajectories of many independent particles.
When a particle enters a node, we model its dynamics as follows. We consider a particle to have entered a node only if it has reached the end of a tubule, say of length l, at which instance we assign the particle to the node-point (i.e. the single point associated with the node in the graph description of the ER network). Although nodes contain a three-dimensional volume, their typical nodal length scale is of the order R l, and thus approximating them by point nodes is appropriate on the scale of the whole network. To decide towards which of the connected tubules the particle leaves the node, we estimate the values of the Péclet number Pe in each of the tubules. We define a local Péclet Pe i = U i R/D where U i is the mean flow velocity through tubule i, with U i > 0 for flow out of the node and U i ≤ 0 otherwise. We then assume that the particle enters a neighbouring tubule i connected to the node with a probability proportional to max(Pe i + 1, 0). This ensures that we have the expected behaviour in both limits of Pe: at high (positive) Péclet numbers, the probability is proportional to the flow speeds in each of the connected tubules, while at low Péclet the exit of the node is limited by diffusion and thus the exit is equally likely in each tubule.
During each simulation, we compute the edge traversal speeds as follows. A particle is defined to traverse an edge (i, j) if it travels from node (i.e. junction) i to node j, or from j to i, without visiting i or j in between. The corresponding edge traversal time is then the time between the arrival at the target node and the most recent departure from the node of origin. The edge traversal speed is naturally defined as the length of the tubule (i, j) divided by the edge traversal time.
The average edge traversal speed associated with an edge (i, j) is then defined as the mean over all edge traversal events across (i, j) of the edge traversal speeds.
In addition we also compute for each particle the 'instantaneous' speeds defined by V n = |X(t n+1 ) − X(t n )|/∆t, where t n = n∆t with ∆t = 18 ms, which is the same temporal resolution as in the particle tracking carried out in Ref. [7].

G. Modelling of alternative flow generation mechanisms
We have described in detail our modelling of a network driven by the pinching of tubules.
In this study we explore two other flow generation mechanisms: the contraction of tubular junctions and the contraction of ER sheets. Our model for pinching tubules is readily generalised to account for these mechanisms, as we describe below.

Experimental estimates of junction volumes
From fluorescence microscopy images of ER networks, we measure the fluorescence intensity of junctions (i.e. the number of pixels in a junction multiplied by the mean intensity per pixel). We can then translate this intensity into an estimate of junction volume, assuming that intensity is directly proportional to volume. In order to calibrate the fluorescence intensity, we use our measurements for tubules. Specifically, we use the measured intensities of tubules and their known volumes, obtained by measurements of tubule lengths and the assumption that they are cylinders of radius 30 nm, in order to determine the proportionality constant to convert between pixel intensity and volume. Our new measurement of the intensities of the junctions then allows us to obtain estimates of their volumes; we obtain twelve values, ranging between 0.0020 µm 3 and 0.0081 µm 3 , with a mean of 0.045 µm 3 and a standard deviation of 0.0021 µm 3 .

Mathematical modelling of contractions of tubular junctions
To include the contribution of tubular junctions into the model (Fig. 1c), we assume that in addition to the pinch sites along tubules, each tubular junction pinches independently of other tubules and other tubular junctions. Given a junction, we assume that it expels the same volume ∆V of fluid during a contraction for all its pinches (and takes in the same volume when it relaxes). Each pinch is assumed to create a sinusoidal flow source, so a pinch lasting for a time 2T produces a flow rate S(t) = ∆V sin(πt/T )π/2T at a time t measured from the beginning of the pinch, where the numerical factors in S(t) are chosen such that the volume expelled during a contraction is indeed T 0 S(t) dt = ∆V . We accommodate the flows from the junctions mathematically by modifying our K1 equations, Eq. (17), to allow the normal nodes to also carry non-zero sources (as opposed to just the pinch nodes, as was the case before); the other equations in the model remain unchanged.

Mathematical modelling of contractions of perinuclear sheets
In the tubules-pinching model, we included M exit nodes located towards the exterior of the network and through which flow could enter and exit the system in order to conserve mass. To account for the connection to a perinuclear sheet (Fig. 1f), we now assign randomly a number of these exit nodes, denoted M 2 < M , to be 'sheet nodes' i.e. nodes which are directly connected to a peripheral sheet, so that a number M 1 = M − M 2 > 0 of exit nodes remain. This is illustrated in Supplementary Fig. 12, where we show the C1 network from Similarly to the mathematical model with tubular pinches only, our independent variables are the E tubule fluxes, the M 1 sources at the exit nodes, and the M 2 sources at the sheet nodes. As before, the K1 (Eq. 17) and K2 (Eq. 18) equations give us E + 1 equations. The requirement that the exit nodes are at the same (exit) pressure give us M 1 2 = M 1 − 1 equations, namely the equality ∆P ij = 0 for i, j exit nodes. We make the additional assumption that the sheet nodes are all at the same mechanical pressure (i.e. that they are connected to the same reservoir), which provides an additional M 2 − 1 equations. Requiring that the sources at the sheet nodes sum to the prescribed flow rate S sheet yields one additional equation, so we again have a system of E + M independent linear equations.

Experimental estimates of volumes of peripheral sheets
To model the contraction of peripheral sheets (Fig. 1d), we need to estimate their contained volumes. Using the open-source image analysis software Fiji [21], we identify nine regions roughly occupied by peripheral sheets in a microscopy image of an ER network; this is illustrated in yellow in Supplementary Fig. 13. We then measure their areas (in µm 2 ) and convert them to volumes by multiplication with the diameter of a tubule; we take the diameter to be 60 nm i.e. we assume that the effective thickness of a sheet is equal to the tubular diameter. From our nine data values, we finally obtain a mean sheet volume of 0.12 µm 3 and a standard deviation of 0.04 µm 3 .
H. Derivation of theoretical bounds for active flows driven by pinching tubules

Advection due to a single pinch
In this section we calculate an upper bound for the axial distance ∆z a particle can be advected by the flow produced by an individual pinch. Recall the formula for the volume 'source' due to a pinch, Among all possible ways for a particle to be transported by a pinch, an upper bound on the transport distance ∆z ≤ ∆z max can be reached if all of the following conditions are satisfied: (i) All of the source flows to one side of the pinch (i.e. there is no leakage on the other side); (ii) The particle travels outside the pinching region (axial flows within the pinching region produce smaller advective displacements than those outside, as may be verified numerically); (iii) The particle travels along the centreline of the tubule (i.e. at twice the cross-section averaged flow velocity, a standard Poiseuille result); (iv) The minimum pinch radius b 0 (recall from Fig. 9) is 0.
Under conditions (i) and (ii), the cross-section averaged speed corresponding to maximal transport generated from a single pinch can then be computed using Eq. (21) as Using condition (iii), the position of the particle along the centreline z(t) satisfies then the ordinary differential equationż Integrating this equation from t = 0 (start of contraction, with b(0) = R) to T (end of contraction, with b(T ) = b 0 ) then leads to The value of ∆z is maximized when b 0 = 0 (i.e. when condition (iv) holds), yielding the upper bound

Extension to nonlinear interactions between two pinches
An isolated pinch is only capable of reciprocal motions. The simplest system capable of producing non-reciprocal motions is illustrated in Fig. 5 and consists of two pinch sites arranged in series near the midpoint of a long horizontal tubule.
We may calculate an upper bound on the net particle displacement that can be achieved after both pinch sites pinch exactly once. We make the assumptions (ii)-(iv) as above, and in addition, allow pinches to spend extended amounts of time in their completely closed state; any deviation from these assumptions will result in net transport that is further reduced.    while that after has length L 2 ; the pinch is symmetric and has a length 2L. Q 1 , Q 2 , Q 3 , and Q 4 denote the volume fluxes through the tubule in the four different regions as indicated.
The pressures at the end of the tubule are p = p 1 at z = 0 and p = p 2 at z = L 1 + 2L + L 2 . formed by adding an edge e ∈ G\T to T . c. A breadth-first search (BFS) starting at the rightmost node; the graph is explored in the order red, green, blue. and M 2 = 9 peripheral sheet nodes (blue asterisks).