To assess the kinetics of particle motion in the lumen of tubular structures, detailed in Figure 1, in response to their contractility, we generate an in silico simulation model of the process. The model incorporates local calculations for the low Reynolds number hydrodynamics of a contracting tubule, assuming in the first instance the no-slip boundary conditions at the tubule walls (i.e. Poiseuille hydrodynamics), into a global analysis of the flows throughout the network geometry, by using Kirchhoff’s laws and standard graph theoretical results (see Materials and methods for details).

We initially implement the model’s numerical simulations of particle transport in a reconstructed ER network of a COS-7 cell (Holcman et al., 2018) (which we label C0, see Network modelling) with tubules locally contracting stochastically according to the spatiotemporal parameters suggested by microscopy measurements (Holcman et al., 2018). We will therein refer to these contractions, illustrated in Figure 1e, as ‘pinching’, with the relevant parameters being the duration and frequency of pinch events, and the length that the pinch sites occupy along the tubule (for details of the pinching kinematics, see Materials and methods, Pinch modelling). An estimate reveals that these pinches are indeed afforded by biologically realistic forces, of the order of 30 pN (see Materials and methods, Estimate of forces required for pinches).

The fluid flows in the edges of the network (model in Hydrodynamic modelling computed as detailed in Solving the hydrodynamic network model) reveal a rapid direction alternation of luminal currents (on average with a frequency of approximately 50Hz), as reflected in the changes of the axial velocity sign (Figure 2a), with an average flow speed of 1.3µm/s (see also Video 1). Furthermore, the resulting instantaneous speed distribution of Brownian particles advected by these flows (methodology in Simulating particle transport) is considerably shifted relative to the experimental counterpart (Holcman et al., 2018) towards lower values (Figure 2b). Similarly, the distribution of average edge traversal speeds (defined precisely in Data processing: instantaneous speeds and average edge traversal speeds) from our simulations (mean 4.4µm/s, solid blue line in Figure 2c) is lower than the experimentally measured speed distribution (mean 45µm/s; Figure 2c, inset) by an order of magnitude.

Figure 2

Download asset Open asset

Video 1

Download asset

Moreover, the results for all measures of transport under pure diffusion, in the absence of pinching-induced flows, are virtually indistinguishable from those where the pinching-driven flows are included (see Figure 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 $\bar{U}\sim 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 $\mathrm{Pe}\sim \bar{U}R/D$. Using $R=30$ nm and the measured diffusivity $D\approx 0.6$ µm² s⁻¹ (Holcman et al., 2018), this leads to the estimate $\mathrm{Pe}\sim 0.07$. Flows affect transport for Péclet numbers above order-one values (Leal, 2007), 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 (Holcman et al., 2018) appear too weak to generate fast, non-diffusive edge traversals.

Importantly, our conclusions remain unchanged upon relaxing the no-slip boundary conditions, an important point to consider since the tubule membranes themselves could be moving in response to the nanoscale fluid flow. This can be modelled by introducing a finite slip boundary condition to couple the membrane-bound lipids with luminal flows (see Materials and methods, Incorporating slip boundary conditions for details). Our simulations show that slip boundary conditions have virtually no effect on the average edge traversal speeds (Figure 3a). Physically, while a non-zero wall slip does change the shape of the flow profiles in the tubules (Figure 3b), it does not modify the cross-sectionally averaged flow speeds which directly affect global particle transport. These results justify our use of Poiseuille hydrodynamics (no-slip boundary conditions) throughout this work.

Figure 3

Download asset Open asset

To estimate the contribution of the network geometry to the outcomes of our transport simulations, we compare results across four different ER structures, which we label C1–C4. We illustrate these reconstructed networks along with the source data in Figure 9. The distributions of the average edge traversal speeds appear 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 (Figure 2d, e). Furthermore, pinching-induced flows inside a regular honeycomb network (Figure 2f, inset) with a typical ER edge length (1 µm) appear to be 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 track homogenisation kinetics by measuring intermixing of particles of two distinct colours equally seeded in each half of a honeycomb network at $t=0$ (Figure 4; see also Video 2 for the flows inside such a network driven by pinches with the parameters from Holcman et al., 2018). The measure of homogenisation is given by the variance $\mathrm{Var}(\varphi (t))\equiv \mathrm{Var}(n_b(t)-n_r(t))$ over 20 regions of the network (Figure 4a, horizontal lines) of the difference between the numbers $n_b$ (blue) and $n_r$ (red) of particles of each colour in region; note $\mathrm{Var}=0$ represents perfect homogeneity. The measures of mixing over time for pure diffusion and active transport with parameters (Holcman et al., 2018) from Figure 2 show a nearly complete overlap (Figure 4e; see also Video 3 and Video 4), 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 Video 5; see also (A combination of high frequency and pinch length is required to replicate experimental particle speeds) for a discussion of the effects of artificially strong pinch parameters on average edge traversal speeds). This indicates that the presumed biological pinching parameters would be inadequate to facilitate ER luminal material exchange.

Figure 4

Download asset Open asset

Video 2

Download asset

Video 3

Download asset

Video 4

Download asset

Video 5

Download asset

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 generators 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 Materials and methods (Advection due to a single pinch), we mathematically show that, given a pinch of length $2L$, the maximum displacement $\mathrm{\Delta }{z}_{\mathrm{m}\mathrm{a}\mathrm{x}}$ a suspended particle may be advected by an individual pinch is

Using the experimentally measured average value of the pinch length $2L=0.14$ µm (Holcman et al., 2018), a typical pinch can then propel a particle by a maximum distance of $\mathrm{\Delta }{z}_{\max }\approx 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 (Holcman et al., 2018) this leads to the estimate $V_{\max }\approx 3.9$ µm/s, which is an order of magnitude smaller than the measured typical edge traversal speed of $\sim 45$ µm/s, consistent with the order-of-magnitude difference between the measurements and the predictions of the computational model.

Note that a slip boundary condition keeps the volume flux expelled out of (or driven into) the pinching regions unchanged, but decreases the maximum flow speed (see Figure 3). Therefore, a non-zero slip length would only decrease the theoretical maximum advective displacement ${\displaystyle \mathrm{\Delta }{z}_{\mathrm{m}\mathrm{a}\mathrm{x}}}$ a suspended particle can achieve; therefore slip length plays no significant role in our theoretical argument to rationalise the weak pinching-induced transport.

The theoretical upper bound in the previous section for the maximum luminal transport producible 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 that is when the other end is effectively blocked. Moreover, 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 Figure 5. We show in Materials and methods (Extension to nonlinear interactions between two pinches) that this is indeed the optimal two-pinch coordination, which results in a time-averaged displacement equal to the upper bound derived in Equation 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\approx 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.

Figure 5

Download asset Open asset

Since the magnitude of ER contractility suggested by microscopy (Holcman et al., 2018) (pinch lengths and frequency) does 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 simulate 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 decrease both the original (Holcman et al., 2018) values of $T$ and $T_{wait}$ by the same factor of $1/\alpha$, where $\alpha \ge 1$, and simulate particle transport in the honeycomb network (Figure 6a). In effect, this simply ‘fast-forwards’ the flows in the original active network by a factor of $\alpha$, and Brownian particles of the original diffusivity are released into this sped-up flow. Instantaneous and edge traversal speeds exhibit corresponding increases when we increase the value of $\alpha$ (Figure 6a). An extreme value of $\alpha =100$ produces an average edge traversal speed distribution that peaks at around 8 µm/s (Figure 6a). Similar results are observed in the C0 network from a COS-7 cell (Figure 6b). The longer tails of these distributions (compared to those from the honeycomb network) result from the variation in edge lengths in the real network, with shorter edges, across which edge traversals are correspondingly fast, contributing to the tails.

Figure 6

Download asset Open asset

These results suggest that, in order to produce average edge traversal speeds of the same order as the experimental values, we would need an active network sped up by an unrealistic factor considerably greater than 100, probably on the order of $\alpha \sim \mathcal{O}({10}^3)$ or higher and corresponding to multiple thousands of pinches occurring on average per second on each tubule. Similarly, it takes an extreme increase in pinch site size spanning the entire length of an average tubule, only to yield an average edge traversal speed of $\sim 10$ µm/s (Figure 6d).

Next, we attempt to obtain a better fit to experimental data by combining changes in both the pinches’ time and geometry parameters. In Figure 6e, f, these maximally long pinches are sped up by a factor of $\alpha =5$ and 10, respectively, which yields speeds averaging around 30 µm/s and above 60 µm/s, respectively. Notably, the tail of the speeds distribution appear longer than that seen in the experiments.

As shown above, establishing effective transport in a tubular constriction-driven model based on realistic ER fluid dynamics requires 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: (1) the tubular junctions (Figure 1c), (2) the perinuclear ER sheets (Figure 1f), and (3) the peripheral sheets (Figure 1d).

First, we run numerical simulations of transport driven by contracting junctions on an ER network (C1 from Figure 9d; sketch of junctions in Figure 1c). Since junctions contracting at the tubular pinch temporal parameters measured experimentally yielded inadequate transport, we consider contractions/relaxations with duration $2T$ exponentially distributed with a mean of ${\alpha }^{-1}$ times the original value in Holcman et al., 2018 and the waiting time $T_{wait}$ between subsequent contractions/relaxations exponentially distributed with a mean ${\beta }^{-1}$ times the original value (i.e. values of $\alpha >1$ and $\beta >1$ reflecting faster and more frequent pinches). Naturally, the results depend on the choice of the volume $\mathrm{\Delta }V$ expelled by a junction during each contraction. Using fluorescence microscopy images, we may estimate the volume in junctions (see Materials and methods, Experimental estimates of junction volumes). We assume that the junction volumes are drawn from a normal distribution with the same mean and SD as our data set that is the distribution N(0.0045 µm³, 0.0021 µm³). Volumes greater than the maximum value in our experimental estimate (0.0081 µm³) or less than the minimum value (0.0020 µm³) are rejected.

We show the results in Figure 7a. The thick solid line illustrates a set of values of $(\alpha ,\beta )$ which produces distributions of average edge traversal speeds with means similar to the experimental values in Holcman et al., 2018, and the dashed lines 1 SD away. We are interested in values of $(\alpha ,\beta )$ close to unity because this corresponds to junctions which pinch with similar pinch durations and frequencies as the experimentally observed tubule pinches and are therefore biologically plausible; $(\alpha ,\beta )\approx$ (2.5,2) is the closest pair on this line to unity. However, only for $\alpha \approx 1$ do we obtain average edge traversal speed distributions of reasonable shapes (i.e. approximately Gaussian), which would require a large $\beta =5$ to match experimental results quantitatively.

Figure 7

Download asset Open asset

Furthermore, the assumption that the entire junction volume is expelled during a contraction may not be realistic, as it would require extreme bending/extension of tubule walls and the ability of the entire junction to empty and fill out during each contraction. Reducing the volume $\mathrm{\Delta }V$ expelled in each contraction to two-thirds of the estimated distribution of the junction’s volume (results shown in Figure 7b) or to half (Figure 7c) causes the average edge traversal speeds to drop considerably. The lower the proportion of the junction volume expelled during a contraction, the faster the pinches are required to be (i.e. large value of $\alpha$) and the larger the frequencies of the pinch events (large $\beta$) in order to produce reasonably high average edge traversal speeds.

Next, we consider the particle transport generated by two types of ER sheets contracting and relaxing over a duration $2T$ (see Materials and 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 Figure 1f). Accordingly, their contraction with $2T=5$ s and $V_{sheet}$ = 10 µm³ yields a mean average traversal speed of 35 µm/s consistent with the single-particle tracking experiments (see Figure 8b). However, the speeds distribution tails towards higher values (Figure 8a, 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.

Figure 8

Download asset Open asset

Furthermore, 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 Figure 8c, d and we further quantify the spatial gradient of the average edge traversal speeds across Cartesian 2D coordinates in Figure 8e, 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 explore whether the peripheral sheets (i.e. the smaller flat inter-tubular ER regions, see sketch in Figure 1d) may overcome the range limit. The peripheral sheets have volumes significantly larger than tubules and junctions, which we estimated at 0.12 µm ± 0.04 µm³ (see Figure 13 and Materials and methods, Experimental estimates of volumes of peripheral sheets 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 contracting junctions (also referred to as ‘nodes’ in our graph theoretical methodology), but with junction volumes set to the measured sheet volumes (see Experimental estimates of volumes of peripheral sheets for details) and with either (1) each node expelling one third of the volume expelled by a contracting sheet (since a peripheral sheet is in contact with three nodes on average), or (2) 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. If we further assume that during each contraction a peripheral sheet expels half its total volume so that in simulation (1) node $k$ expels a volume $V_k/6$ in each contraction and in simulation (2) each active node $k$ expels a volume $V_k/2$, then nodes which contract at rates ${\alpha }^{-1}=2.5$ and ${\alpha }^{-1}=5$ times slower than the tubule pinches in Holcman et al., 2018 in simulations (1) and (2), respectively, give average edge traversal speeds in the correct range of 40 µm/s (Figure 8f, 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.

Furthermore, by relating the work done by a peripheral sheet contraction to the dissipation due to the flows induced inside the sheet and in the rest of the network, the energetic cost of one such contraction may be estimated to be of the order of 1000 molecules of ATP (see Materials and methods, Energetic cost estimate for contracting peripheral sheets for details). For reference, a kinesin motor protein uses one molecule of ATP to move 8 nm (Coy et al., 1999; Schnitzer and Block, 1997), whereas a muscle fibre consumes hundreds of thousands of ATP per second (Barclay, 2017). These peripheral sheet contractions may be directly ATP-driven, or they might come as a result of other mechanical processes in the cell, similar to mechanisms generating cytoplasmic mixing flow (Klughammer et al., 2018; Guo et al., 2014; Lin et al., 2016).

We represent the ‘skeleton’ of a two-dimensional ER network as a planar graph with each node assigned a position $𝐱\in {ℝ}^2$. Given an edge of the network labelled $(i,j)$ and of length $|{𝐱}_i-{𝐱}_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 Supplementary file 3 of Holcman et al., 2018; this network is reproduced in Figure 9a. We use the multi-point tool in ImageJ (Schneider et al., 2012) 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 Figure 9b, with the original network superimposed on the mathematical model in Figure 9c. 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 Figure 9d). In Figure 9e, we show the distributions of the edge (tubule) lengths in each of the C0–C4 networks, as well as the mean edge lengths, the mean degrees (the degree of a node is the number of edges connected to it), and the number of nodes. 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 a honeycomb network, that is one where every node (apart from those at the boundaries) has a degree of 3, with all edge lengths exactly 1 µm.

Figure 9

Download asset Open asset

To describe the pinches (Figure 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 Figure 10. Each tubule is assumed to have a pinch site at its midpoint (i.e. with $L_1=L_2$); we make this simplifying assumption in our simulations, after verifying that more general pinch locations ($L_1\ne L_2$) have virtually no effect on particle transport. The pinching events occur at the pinch sites stochastically and independently of each other. Each pinch is defined geometrically by three parameters from a random distribution (see below): (1) the duration of a pinch $2T$, (2) the time $T_{wait}$ between the end of a pinch and the beginning of a new one on the same site, and (3) the length of a pinch $2L$ (see Figure 10). We assume for simplicity that all pinches are axisymmetric so that using the notation of Figure 10, 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 cross-section through the centre of the pinch (i.e. in Figure 10 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$.

Figure 10

Download asset Open asset

We model the geometrical profile of each pinch of length $2L$ as following a linear radius change (see Figure 10). Within a pinch located at $z=z_0$, the radius of the cylinder is therefore given by

The time-varying function $b(t)\le R$ is therefore the minimum pinch radius in the centre of the tubule. As the simplest modelling choice, we assume in our simulations that $b$ changes in time sinusoidally and thus will use $b(t)=(R+b_0)/2+(R-b_0)\cos (\pi t/T)/2$, where $t$ is time after the pinch begins and $b_0$ is the value of $b$ halfway in time through the pinch. Choosing a smooth time variation for the function $b(t)$ will ensure the continuity of fluxes in time (see Equation 10 later). We have verified that changing the pinch shape to a smoother geometry has essentially no impact on the fluid volume expelled/taken in, and therefore no significant effect on the flows/transport.

In our simulations, the stochastic pinch parameters are drawn from the distributions measured in Holcman et al., 2018. The pinch duration $2T$ is therefore drawn from an exponential distribution with rate parameter $\lambda =\ln 10/0.167$ s⁻¹ (so the mean value is $\overline{2T}=0.0725$ s), the time between pinches $T_{wait}$ from an exponential distribution with rate parameter ${\lambda }_{wait}=\ln 10/0.851$ s⁻¹ (mean value $\overline{T_{wait}}=0.370$ s) and the pinch length $2L$ from a uniform distribution with mean µm and variance ${\displaystyle {\sigma }^2=({0.06}^2/\ln 10)\mu {\mathrm{m}}^2\approx (0.040\mu \mathrm{m}{)}^2}$. Throughout we set the tubule radius $R=30$ nm and $b_0=0.01R$ (recall from Figure 10).

The incorporation of pinches as ‘dummy nodes’ into the graph theoretical framework of Network modelling along with Equations 10 and 14 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 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,\dots ,M\}$ denotes the $M$ exit nodes. Let $q_i$ be the source or sink carried by the $i^{\mathrm{th}}$ node (so that $q_i=0$ if $i$ is a normal node, $q_i$ is as specified by the RHS of Equation 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,\dots ,M\}$ (that is, 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 Kirchhoff’s laws.

Kirchhoff’s second law (K2)

Request a detailed protocol

The second equation is a statement of 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 $C=\{v_1,v_2,\dots ,v_n,v_{n+1}=v_1\}$, if $\mathrm{\Delta }{p}_{v_i{v}_{i+1}}$ denotes the pressure change from node $v_i$ to node $v_{i+1}$, we necessarily have

(18) $\begin{array}{ccc}& \sum _{i=1}^n\mathrm{\Delta }{p}_{v_i{v}_{i+1}}=0.& \end{array}$

Note that the pressure change across a node is negligible.

The K2 statement in Equation 18 applies to all cycles in the graph, which would give us more equations than we need since the vectors of coefficients in Equation 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 (Wilson, 2015). 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 Figure 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 Figure 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.

Figure 11

Download asset Open asset

To compute the spanning tree $T$ we use a breadth-first search algorithm (Goodrich and Tamassia, 2015). 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 Figure 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.

With our solution for the flows in the active network at each instant of time, we now proceed to track the motion of Brownian particles inside the network using a discretisation of their stochastic equations of motion, as a model for the transport of proteins in the ER network.

We use the simplest approach where we superimpose Brownian motion onto advection by the flow inside each tubule. Let $𝐱(t)$ denote the position of a Brownian particle in a tubule and ${𝐱}_n$ the finite-difference approximation of $𝐱(n\mathrm{\Delta }t)$, where $\mathrm{\Delta }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 $𝐔$ is the instantaneous flow velocity, and the random noise term $𝐗(\mathrm{\Delta }t)$ is drawn from a zero-mean Gaussian with variance $⟨𝐗(\mathrm{\Delta }t)𝐗(\mathrm{\Delta }t)⟩=2D\mathrm{\Delta }t𝐈$, 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 Holcman et al., 2018, $D\approx 0.6$ µm² s⁻¹. 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\ll 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 $\mathrm{Pe}$ in each of the tubules. We define a local Péclet ${\mathrm{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\le 0$ otherwise. We then assume that the particle enters a neighbouring tubule $i$ connected to the node with a probability proportional to $\max ({\mathrm{Pe}}_i+\mathrm{1,0})$. This ensures that we have the expected behaviour in both limits of $\mathrm{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=|𝐗(t_{n+1})-𝐗(t_n)|/\mathrm{\Delta }t$, where $t_n=n\mathrm{\Delta }t$ with $\mathrm{\Delta }t=18$ ms, which is the same temporal resolution as in the particle tracking carried out in Holcman et al., 2018.

The methodology we have detailed thus far assumes no-slip boundary conditions at the tubule walls for the fluid flow. However, the membrane-bound lipids themselves could also flow in response to the nanoscale luminal flows. This may be modelled by introducing a finite slip boundary condition on the tubule wall. The slip boundary conditions with a slip length $\lambda \ge 0$ at the wall $r=a(z,t)$ are given by

and $(𝐮-\frac{\partial a}{\partial t}{𝐞}_r)\cdot 𝐧=0$, with $𝐧=\frac{\partial a}{\partial z}{𝐞}_z-{𝐞}_r$, which simplifies to

We may then derive the long-wavelength solution for the flow field inside an axisymmetric deforming tubule as follows. Using an ansatz for the axial component $u$ that is motivated by the uniform-radius Poiseuille flow with slip,

we may solve the incompressibility condition $\frac{1}{r}\frac{\partial }{\partial r}(rv)+\frac{\partial u}{\partial z}=0$. Regularity at $r=0$ constrains the integration constant to be zero, yielding

which automatically satisfies the boundary condition in Equation 22. Note that the mass conservation equations are not affected by the introduction of a slip length.

Using the new solution for $u$, the modified Hagen–Poiseuille expression with slip may be derived as before to be

which does recover the no-slip result as $\lambda \to 0$.

These results may then be used to simulate particle transport with slip boundary conditions. In Figure 3a, we plot the distributions of average edge traveral speeds obtained from simulations of a C1 network pinching with the original pinch parameters from Holcman et al., 2018 for four different values of the boundary slip length. In Figure 3b, we further display the profiles of the longitudinal flow (Equation 23) for different slip lengths with the volume flux fixed.

In this section, we derive an order-of-magnitude estimate for the forces required to pinch an ER tubule.

We first estimate the difference in the membrane’s elastic energy, $\mathrm{\Delta }E=E_{\text{pinched}}-E_{\text{unpinched}}$, between the pinched ($E_{\text{pinched}}$) and unpinched configurations ($E_{\text{unpinched}}$). In the absence of spontaneous curvature, the Helfrich free energy $h$ per unit area of a membrane is given by

where $k_c$ and $\bar{k}$ are bending rigidities, $H$ is the mean curvature, and $K$ is the Gaussian curvature (Helfrich, 1973). The mean and Gaussian curvatures may be expressed in terms of the principal curvatures ${\kappa }_1,{\kappa }_2$ as $H=({\kappa }_1+{\kappa }_2)/2$ and $K={\kappa }_1{\kappa }_2$. We take ${\kappa }_1$ and ${\kappa }_2$ to be the principal curvatures in the directions normal and parallel, respectively, to the tubule’s longitudinal axis.

The dominant contribution to $\mathrm{\Delta }E$ is from $E_{\text{pinched}}$, specifically from the region near the centre of the pinch, taken to be at $z=0$, where the tubule radius is smallest. We have ${\kappa }_1=1/a(z,t)$ and ${\kappa }_2=\mathcal{O}(R/L^2),$ yielding $H=1/2a+\mathcal{O}(R/L^2)$ and $K=\mathcal{O}(R/a{L}^2)$. The dominant contribution to the membrane energy $E=2\pi {\int }_{-L}^{L}ha\mathrm{d}z$ is therefore from $H$, in a neighbourhood of $z=0$. Since real pinches do not have kinks, a linear term, for the purposes of estimating membrane energy, is unphysical, and we expect $a(z,t_0)=b_0+\mathcal{O}(z^2)$ near $z=0$, where $b_0$ is the pinch radius in the centre of the pinch in the maximally contracted state and $t_0$ denotes a time at which the tubule is in a pinched state. We may therefore write $a(z)=b_0+R{z}^2/L^2$ to obtain an order-of-magnitude estimate of the membrane energy in the pinched state.

Evaluating the integral for $E_{\text{pinched}}$ then yields the leading-order estimate

We may take $k_c,\bar{k}\sim 50k_{B}T$ ($k_B$ is the Boltzmann constant times and $T$ the room temperature) (Faizi et al., 2019), with estimate values ${\displaystyle R=30\mathrm{n}\mathrm{m}}$ and ${\displaystyle L=70\mathrm{n}\mathrm{m}}$. The minimum pinch radius $b_0$ may be estimated to be 10 nm allowing for membrane thickness and incomplete squeeze, yielding

The hydrodynamic contribution to the energy expenditure during a pinch may be calculated as the sum of the dissipation inside the pinch and in the rest of the network (see Energetic cost estimate for contracting peripheral sheets for an analogous calculation for a contracting peripheral sheet). We find that the hydrodynamic contribution is negligible compared to the elastic component of the work done to pinch a tubule.

An estimate for the force required to pinch the tubule may then finally be obtained as ${\displaystyle F\sim \mathrm{\Delta }E/R\sim 30\mathrm{p}\mathrm{N}}$.

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.

Mathematical modelling of contractions of perinuclear sheets

Request a detailed protocol

In the tubule-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 (Figure 1f), we now assign a number of these exit nodes, denoted $M_2<M$, to be ‘sheet nodes’, i.e. nodes which are directly connected to a perinuclear sheet, so that a number $M_1=M-M_2>0$ of exit nodes remain. This is illustrated in Figure 12, where we show the C1 network from Figure 9 with both sheet nodes (blue asterisks) and exit nodes (red squares).

Figure 12

Download asset Open asset

A sheet contraction + relaxation lasting a time $2T$ produces a total source $S_{\mathrm{sheet}}(t)=V_{\mathrm{sheet}}\pi \sin (\pi t/T)/2T$ at a time $t$ from the beginning of the pinch, where again, the integral of $S_{\mathrm{sheet}}$ over a contraction gives a volume $V_{\mathrm{sheet}}$ of fluid expelled by the sheet.

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 (Equation 17) and K2 (Equation 18) equations give us $E+1$ equations. The requirement (analogous to Equation 19) that the exit nodes are at the same mechanical pressure gives us ${\displaystyle M_1-1}$ equations. 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_{\mathrm{sheet}}$ yields one additional equation, so we again have a system of $E+M$ independent linear equations.

Experimental estimates of volumes of peripheral sheets

Request a detailed protocol

To model the contraction of peripheral sheets (Figure 1d), we need to estimate their contained volumes. Using the open-source image analysis software Fiji (Schindelin et al., 2012), we identify nine regions roughly occupied by peripheral sheets in a microscopy image of an ER network; this is illustrated in yellow in Figure 13. We then measure their areas (in µm²) and convert them to volumes by multiplication with the diameter of a tubule, taken as 60 nm, assuming 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³ and an SD of 0.04 µm³.

Figure 13

Download asset Open asset

The dominant energy expenditure in deforming a pinching tubule is in creating the large membrane curvatures at the narrowest sections of the pinch, and much less work is done against the small volumes of fluid displaced (Estimate of forces required for pinches). A contracting peripheral sheet, however, displaces a relatively large volume of fluid without attaining the extreme curvatures required in a pinching tubule. We therefore expect, intuitively, that the hydrodynamic contribution will dominate the energy expenditure budget. We now explicitly show this by estimating both elastic and hydrodynamic contributions.

To derive an order-of-magnitude estimate of the energetic cost, we consider an idealisation of a peripheral sheet consisting, in the relaxed state, of two parallel circular membranes $r<R$ located at $z=\pm D$, and in the fully contracted state, of two paraboloids at $z=\pm D{r}^2/R^2$ (see Figure 14). The membranes deform as paraboloids between these two states, and we denote the ‘vertical’ distance between the two membranes as $d(r,t)$. We take $R=0.8$ µm, the value consistent with a sheet thickness $2D=60$ nm and the mean sheet volume (in the fully relaxed state) of $2\pi R^{2}D=0.12$ µm³ (see estimation in Experimental estimates of volumes of peripheral sheets).

Figure 14

Download asset Open asset

To estimate the elastic contribution, we present an argument similar to the one carried out in Estimate of forces required for pinches. In the fully contracted state, the principal curvatures of one membrane scales as $\sim 2D/R^2$, so the Gaussian curvature $K\sim 4D^2/R^4$ and the mean curvature $H\sim 2D/R^2$, yielding the Helfrich energy density $h\sim k[2H^2+K]\sim 12k{D}^2/R^4$, where $k$ is the typical bending rigidity. The total bending energy of the fully contracted membrane is then ${\displaystyle E\sim \pi R^{2}h\sim 2\times {10}^{-20}\mathsf{J}}$. In the relaxed state, the membranes are flat and have zero bending energy. Thus, the total energy required to contract the two membranes can be estimated as

To estimate the hydrodynamic contribution, we first note that the work done instantaneously by the contracting sheet against the fluid, of dynamic viscosity μ, say, is the sum of dissipation rate in the sheet itself and the dissipation rate in the rest of the ER network outside the sheet due to the sheet-induced flows. Any work done against the fluid outside the ER network is neglected since flows decay over short length scales across the network.

We first consider the contribution from within the sheet. We may use lubrication theory to estimate the flows inside a contracting peripheral sheet since the sheet is fairly flat. The leading-order lubrication flow between the contracting membranes (located at $z=\pm d(r,t)$) is in the radial direction, and given by $u_r=-\frac{1}{\mu }\frac{\partial p}{\partial r}(d^2-z^2)$, where the unknown radial pressure gradient $\frac{\partial p}{\partial r}$ may be calculated from mass conservation $2\pi r{\int }_{-d}^d{u}_r\mathrm{d}z=-\frac{\partial }{\partial t}({\int }_0^{r}4\pi r^{\prime }d(r^{\prime },t)\mathrm{d}{r}^{\prime })$, thus yielding a radial velocity

The dominant contribution to the rate-of-strain tensor is $e_{rz}\approx \frac{1}{2}\frac{\partial u_r}{\partial z}\sim \frac{3R}{2DT}$, where $T\sim 0.5$ s is the contraction duration (in the main text we considered a value of $T$ of 2.5 or 5 times larger than the experimentally measured pinch duration; here we take 5). Therefore, the total work done over a contraction scales as the total dissipation rate times the contraction duration $T$, and we scale

Taking the viscosity of the intra-luminal fluid to be 10 times that of water, this is computed to be

The dissipation rate in the network outside the sheet may be estimated by calculating the total dissipation rate in an idealised network consisting of ‘generations’ of tubules, each of length 1 µm, with the first generation consisting of three tubules connected to the peripheral sheet, and each tubule in the $i^{\mathrm{th}}$ generation branching out into three tubules of the $(i+1{)}^{\mathrm{th}}$ generation ad infinitum. The dissipation rate due to a Poiseuille flow of flux $q$ inside a tubule of length $l$ and radius $R$ may be calculated to be

Denoting by $Q$ the total volume flux created by the contracting sheet, we may calculate the dissipation rate in the network of tubules by summing the above expression across all the tubules as follows. The $i^{\text{th}}$ generation of tubules is comprised of $3^i$ tubules each carrying a flux of $Q/3^i$. Summing the total dissipation rates across all generations gives

We may scale $Q$ as $Q\sim V_{\text{s}}/2T$, recalling our assumption that each peripheral sheet contraction expels half of the volume $V_s$ contained in the sheet (which is consistent with the paraboloidal membrane profile we have taken for the fully contracted sheet). The total dissipation over the contraction duration $T$ in this network of tubules is therefore

Scaling the sheet volume with its mean value $V_s\sim 0.12$ µm³, we may compute this to be

The total work done during a contraction is therefore ${\displaystyle W=W_{\text{inside}}+W_{\text{outside}}\sim {10}^{-16}\mathsf{J}}$. Given that the energy released by the hydrolysis of one ATP molecule is of the order of ${\displaystyle {10}^{-19}\mathsf{J}}$ (Bray, 2001), we thus estimate that each contraction of a peripheral sheet would require on the order of 1000 molecules of ATP.

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

EL and PHH are supported by the European Research Council under the European Union’s Horizon 2020 research and innovation program (Grant No. 682754 to EL) and by the Cambridge Trust. EA is supported by the UK Dementia Research Institute [award number UK DRI-2004] which receives its funding from UK DRI Ltd, funded by the UK Medical Research Council, Alzheimer’s Society and Alzheimer’s Research UK.

1. Preprint posted: October 21, 2023
2. Sent for peer review: December 7, 2023
3. Reviewed Preprint version 1: April 10, 2024
4. Reviewed Preprint version 2: October 4, 2024
5. Version of Record published: December 13, 2024

You can cite all versions using the DOI https://doi.org/10.7554/eLife.93518. This DOI represents all versions, and will always resolve to the latest one.

© 2024, Htet et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.