Active elastohydrodynamics of vesicles in narrow, blind constrictions

TL;DR

This paper investigates the active transport of vesicles through narrow, blind constrictions using simulations and theory, revealing how geometry, forces, and hydrodynamics influence transit times and dynamics.

Contribution

It introduces a combined simulation and semi-analytical approach to understand vesicle transport in constricted geometries, emphasizing the role of hydrodynamics and motor forces.

Findings

Transport into blind ends forms lubrication layers that impede motion.

Motor and fluid forces create multistable vesicle dynamics.

Hydrodynamic effects significantly influence vesicle transit kinetics.

Abstract

Fluid-resistance limited transport of vesicles through narrow constrictions is a recurring theme in many biological and engineering applications. Inspired by the motor-driven movement of soft membrane-bound vesicles into closed neuronal dendritic spines, here we study this problem using a combination of passive three-dimensional simulations and a simplified semi-analytical theory for active transport of vesicles that are forced through such constrictions by molecular motors. We show that the motion of these objects is characterized by two dimensionless quantities related to the geometry and the strength of forcing relative to the vesicle elasticity. We use numerical simulations to characterize the transit time for a vesicle forced by fluid pressure through a constriction in a channel, and find that relative to an open channel, transport into a blind end leads to the formation of an effective lubrication layer that strongly impedes motion. When the fluid pressure forcing is complemented by forces due to molecular motors that are responsible for vesicle trafficking into dendritic spines, we find that the competition between motor forcing and fluid drag results in multistable dynamics reminiscent of the real system. Our study highlights the role of non-local hydrodynamic effects in determining the kinetics of vesicular transport in constricted geometries.

Tables

Table 1. Table 1: Lattice Boltzmann parameters.

Symbol	Definition	Value in lattice units	Value in SI units
$\delta t$	Time step	1	4.4 $\times {10}^{-10}$	s
$\delta x$	Lattice constant	1	$2.3\times {10}^{-8}$	µm
$\mu$	Fluid dynamic viscosity	1/6	$2\times {10}^{-4}$	$\text{Pa}\cdot \text{s}$
$\rho$	Fluid mass density	1	$1\times {10}^3$	kg/m3
$F_{\text{ex}}$	Applied force	1 – 5	(1.5 – 7.5) $\times {10}^{-9}$	N
${\kappa }_{\alpha }$	Area dilation modulus	1	0.06	N/m
${\kappa }_s$	Shear modulus	0.015	$1\times {10}^{-3}$	N/m
${\kappa }_B$	Bending modulus	0.018	$5\times {10}^{-19}$	J
$R_c$	Channel radius	30 – 40	($6.9$ – $9.2)\times 10{}^{-1}$	µm
$R_p$	Vesicle radius	15 – 32	($3.5$ – $9.4)\times 10{}^{-1}$	µm
$L$	Domain length	384	8.9	µm

Table 2. Table 2: Lubrication model parameters.

Symbol	Definition	Value	Units
$\mu$	Fluid viscosity	$1.2\times {10}^{-3}$	$\text{Pa}\cdot \text{s}$
$\rho$	Fluid density	$1\times {10}^3$	$\text{kg}/{\text{m}}^3$
$F_{\text{ex}}$	Applied force	$40$ – $200$	pN
$C$	Compliance	$5\times {10}^{-9}$	$\text{m}/\text{Pa}$
$R_c$	Channel radius	$1.22$ – $2.15$	µm
$R_p$	Vesicle radius	$0.96$ – $1.5$	µm
$l_{\text{narrow}}$	Narrow channel length	$2.5$	µm
$l_{\text{trans}}$	Transition length	$2.5$	µm

Equations

$$f_i\left(\vec{r}+{\vec{e}}_i\delta t,t+\delta t\right)-f_i\left(\vec{r},t\right)={\Omega }_i,$$

$$f_i^{eq}(\rho ,\vec{u})={\omega }_i\rho \left[1+\frac{{\vec{c}}_i\cdot \vec{u}}{c_s^2}+\frac{({\vec{c}}_i\cdot \vec{u}{)}^2}{2c_s^4}-\frac{\vec{u}\cdot \vec{u}}{c_s^2}\right],$$

$$E_S+E_B=\int {ϵ}_{s}dA+\frac{{\kappa }_B}{2}\int H^{2}dA,$$

$$u(z)=\frac{1}{2\mu }\frac{\partial p}{\partial z}r(r-h)+\frac{U}{h}r,$$

$$Q=2\pi R_c\left(\frac{-h^3}{12\mu }\frac{\partial p}{\partial z}+\frac{1}{2}Uh\right)=\text{constant}.$$

$$\frac{p(z)-p_0}{6\mu }=U{\int }_{Z-R_p}^z\frac{1}{h^2(s)}\mathrm{d}s-\frac{2Q}{2\pi R_c}{\int }_{Z-R_p}^z\frac{1}{h^3(s)}\mathrm{d}s.$$

$$Q=2\pi R\left(U\int_{Z-R_{p}}^{Z+R_{p}}\frac{1}{h^{2}(s)}\mathrm{d}s-\frac{F}{6\pi R^{2}\mu}\right)\Bigg{/}\left(2\int_{Z-R_{p}}^{Z+R_{p}}\frac{1}{h^{3}(s)}\mathrm{d}s\right).$$

$$Q=-\pi R_c^{2}U,$$

$$h(z)={\tilde{R}}_c(z)-\sqrt{R_p^2-(z-Z{)}^2}+C(p(z)-p_0),$$

$$\frac{p^{i+1}(z)-p_0}{6\mu }$$

$$h^{i+1}(z)$$

$$Q^{i+1}$$

$$U^{i+1}$$

$$\frac{F}{6\pi \mu R_p^2}=U\left({\int }_{Z-R_p}^{Z+R_p}\frac{1}{h^2(s)}\mathrm{d}s+R_c{\int }_{Z-R_p}^{Z+R_p}\frac{1}{h^3(s)}\mathrm{d}s\right).$$

$$\begin{array}{rl}{\int }_{Z-R_p}^{Z+R_p}\frac{1}{h^n(s)}\mathrm{d}s& ={\int }_{-R_p}^{R_p}\frac{1}{(R_c-R_p\sqrt{1-(s/R_p{)}^2}{)}^n}\mathrm{d}s\\ & \approx {\int }_{-R_p}^{R_p}\frac{1}{(R_c-R_p(1-(s/R_p{)}^2/2){)}^n}\mathrm{d}s\\ & \approx \frac{1}{h_0^n}{\int }_{-R_p}^{R_p}\frac{1}{(1+(R_p/2h_0)(s/R_p{)}^2{)}^n}\mathrm{d}s\\ & \approx \frac{R^{1/2}}{h_0^{n-1/2}}{\int }_{-R_p/\sqrt{h_0{R}_p}}^{R_p/\sqrt{h_0{R}_p}}\frac{1}{(1+y^2/2{)}^n}\mathrm{d}y\approx \frac{R^{1/2}}{h_0^{n-1/2}},\end{array}$$

$$\begin{array}{rl}\frac{h(s)}{R_c}& \approx \frac{R_c-R_p\left(1-s^2/(2R_p^2)\right)+Cp(s)}{R_c}\\ & \approx 1-{\pi }_1+\frac{s^2}{2R^2}+{\pi }_2/2\\ & \approx 1-{\pi }_1+\frac{s^2}{2R^2},\end{array}$$

$${\theta }_{[a,b)}={\int }_a^b\varphi (z)\mathrm{d}z.$$

$$\theta ={\int }_{-\infty }^{\infty }\varphi (z)\mathrm{d}z,$$

$$F=n_0{\int }_{-\infty }^{\infty }f(z)\varphi (z)\mathrm{d}z,$$

$$\alpha (1-\theta )=\beta {\int }_{z_0}^A\varphi (z)\mathrm{d}z-U\varphi (z_0),$$

$$U\mathrm{d}\varphi /\mathrm{d}z=-\beta \varphi$$

$$\varphi (z)=\frac{\alpha \beta }{-U(\alpha +\beta )}\exp \left(-\beta (z-A)/U\right).$$

$$F=-\frac{\alpha n_0{p}_1}{\alpha +\beta }\frac{\left(e^{\gamma A}-1\right)+\left(\gamma U/\beta \right)}{1-\left(\gamma U/\beta \right)},$$

$$U=\frac{b(F_0-F)}{F+a},$$

$$a$$

$$b$$

$$F_0$$

$$\alpha (1-\theta )=\beta {\int }_A^{z_0}\varphi (z)\mathrm{d}z+U\varphi (z_0).$$

$$\varphi (z)=\frac{\alpha \beta }{U(\alpha +\beta )}\exp \left(-\beta (z-A)/U\right).$$

$$F=-\frac{\alpha n_0{p}_1}{\alpha +\beta }\frac{\left(e^{\gamma A}-1\right)+\left(\gamma U/\beta \right)}{1-\left(\gamma U/\beta \right)},$$

$$\theta ={\int }_A^B\varphi (z)\mathrm{d}z=\frac{U}{\beta }\varphi (A)\left(\exp \left(\beta (B-A)/U\right)-1\right),$$

$$\varphi (z)=\frac{\alpha \beta c}{-U(\alpha +\beta c)}\exp \left(-\beta (z-A)/U\right),$$

$$F=n_0{\int }_A^{B}f(z)\varphi (z)$$

Full text

1 Introduction

Intracellular transportation of vesicles frequently involves translocation through channels and into narrow pockets. An example of such a process that has received considerable experimental attention with the advent of live-cell imaging is the motor-driven transportation of protein-rich vesicular endosomes through the necks of dendritic spines in neurons [40, 47, 11]. The micron-sized vesicles deform strongly during transport as they are squeezed through the sub-micron sized necks by molecular motors. These vesicles play a critical role in transport of chemicals responsible for the normal functioning of neuronal synapses, and understanding the biophysical basis for and limits to their transport is thus an important problem. Indeed, this situation occurs not only in intracellular trafficking, but also in numerous cellular and microfluidic settings in which elastic bodies, such as manufactured elastic capsules [12, 13], hydrogels [32] or living cells [16, 8, 3] squeeze through narrow constrictions. While driving forces such as pressure gradients, fluid flow, molecular motors, or external fields may promote passage, the changes in shape required for passage result in energetic barriers to transportation, altering the passage time and eventually leading to non-passage of the vesicle. Such transport processes are thus determined by a subtle balance between forces originating from various physical sources and the geometry of the constriction.

In recent years, considerable progress has been made in studying this process computationally in the context of the translocation of elastic capsules through cylindrical [43, 41, 30, 44] and square [29, 36] micro-channels. These studies typically consider transit through open channels (i.e. infinite channels or those with periodic boundary conditions at each end) in contrast with transportation into a closed constriction, i.e. one with a blind end, as occurs in dendritic spines (Fig. 1(a) and (b)). Dendritic spines are bulbous structures protruding from dendrites. The heads of dendritic spines are the locations of postsynaptic densities, which are membrane receptor-enriched regions responsible for detecting neurotransmitters released into the synaptic cleft. Therefore a flux of membrane receptors to the spine head is required to maintain synaptic function, and these receptors are actively transported to dendritic spines on membrane-bound vesicles known as recyling endosomes. The nearly occlusive movements of these elastic vesicles through the narrow spine necks leads to large deformations of the vesicle. Moreover, as a consequence of the spine geometry, when the vesicle enters the constriction either passively due to an ambient pressure or due to the activity of motors, fluid enclosed by the spine must flow out of the same constriction in a narrow layer surrounding the vesicle. This leads to long-range interactions mediated by fluid incompressibility, a problem that is the particular focus of this study.

Here, we quantify how the translocation time through open and closed channels with a constriction depends on the constriction size, vesicle elasticity, and applied force. In Section 2 we outline the abstraction of the experimental observations of endosomal transport. In Section 2.1 we present a lattice Boltzmann method that we use to study the translocation dynamics of a soft vesicle, focusing on how the applied force and channel geometry impact the transit speed of the vesicle (Fig. 1(c)). In Section 3 we describe a simplified model based on lubrication theory that captures the essential features of translocation determined using lattice Boltzmann simulations. In Section 3.2, we apply the lubrication model to the phenomenon of intracellular cargo transportation by replacing the assumption of prescribed forcing by a realistic description of molecular motors, and in Section 4 we compare the results of simulations to experimental observations on vesicle transportation into dendritic spines. We conclude with a discussion of how our model captures the diverse behaviors observed in vesicle translocation through spines, including vesicle rejection, corking, and translocation, and we show how different behaviors can be obtained by varying model parameters such as the constriction geometry and the dynamics of the molecular motors.

2 Computational formulation of 3D endosomal transport

Translocation of elastic vesicles through narrow channels is critical in the context of neuronal regulation associated with synaptic activity. The vesicles, also called endosomes, serve to transport receptors from the cell body along dendritic spines and finally undergo exocytosis onto the spine membrane. Typical endosomes have diameters around 1-2 µm whereas spine necks have diameters of several hundred nm, and the forces driving translocation are generated by families of molecular motors that move bidirectionally [41]. The diversity of molecular constituents helps account for the rich set of behaviors observed experimentally in fluorescently labeled vesicles transiting into and out of the spine. These behaviors include positive and negative translocation velocities, and step-like transitions between these velocities [11] suggestive of a discrete set of possible values. Additionally, the small aspect ratio of the long narrow dendritic spines causes the vesicles to deform as they squeeze through the neck, as can be seen in the timelapse images and kymograph in Fig. 1(a) and (b). Although much progress has been made on cataloguing the types of cytoskeletal molecules present [25, 11], such as molecular motors, the interplay of geometry, hydrodynamics and motor activity that governs how vesicles are forced through the spine neck is still poorly understood.

To begin analyzing this system, we consider a simple model of translocation, wherein an elastic vesicle of radius $R_{p}$ is driven by an applied force $F$ through a rigid constriction in a channel filled with a Newtonian fluid. We define the constriction by a single scale, e.g. a radially symmetric bump of radius $R_{c}$ at its narrowest point (Fig. 1(c)). We define the confinement ratios $\pi_{1}$ by $\pi_{1}:=R_{p}/R_{c}$ and the dimensionless forcing $\pi_{2}$ by $\pi_{2}:=CF/(\pi R_{p}^{3})$, where the compliance $C$ represents the elasticity of the vesicle, so that the two dimensionless quantities $\pi_{1}$ and $\pi_{2}$ govern the translocation dynamics. In this work, we restrict attention to confinement ratios $\pi_{1}<1$, since this case already involves significant computational challenges in resolving the lubrication layer and since nearly spherical vesicles allow for several simplifying theoretical assumptions. In the Discussion section, we comment on why lessons learned from studying this case are likely to apply more generally to the case of extreme confinement ratios $\pi_{1}\gg 1$ in which the vesicle must undergo large deformations simply to fit inside the channel.

2.1 Lattice Boltzmann model

Before moving to a simplified mathematical description of this problem, it is useful to perform full-scale 3D simulations to get a sense of the elastohydrodynamical effects associated with the motion of a soft vesicle through a narrow constriction. We use the lattice Boltzmann approach as our fluid solver and a finite element method to compute the elastic response of the vesicular membrane. These are coupled via an explicit immersed boundary method for the interaction of the fluid and the elastic membrane. The 3D lattice Boltzmann solver considers the fluid as a cluster of pseudo-particles that move on a lattice under the action of external forces. The fluid is represented by a distribution function $\vec{f}_{i}$ that represents the probability of finding a pseudo-particle at position $\vec{r}$ with velocity in direction $\vec{e}_{i}$. The position and velocity spaces are both discretized on a lattice with spacing $\delta x$, and we use the so-called D3Q19 lattice, which has nineteen velocity directions [4, 42]. The time evolution of the force distribution ${f}_{i}$ is governed by

[TABLE]

where $\delta t$ is the discrete time step and $\Omega_{i}$ is the collision rate between the fluid pseudo-particles, which is approximated by the Bhatnagar-Gross-Krook operator, $\Omega_{i}=-\delta t(f_{i}-f_{i}^{eq})/\tau_{r}$. Here $\tau_{r}$ is the relaxation time, related to the dynamic viscosity by $\mu=\rho c_{s}^{2}\delta t(\tau_{r}/\delta t-1/2)$, where $c_{s}=1/\sqrt{3}\delta x/\delta t$ is the lattice speed of sound and $\rho$ is the density of the fluid. One can convert between lattice units and SI units using these relations; for our calculations, we set the lattice constant, timestep, mass, and relaxation time to unity (see [23] for more information on the conversion between physical and lattice units).

The equilibrium distribution is given by the truncated Maxwell-Boltzmann distribution,

[TABLE]

where $\vec{u}$ is the velocity vector and the $\omega_{i}$ are weight factors that result from the velocity space discretization. Using this relation, we can calculate macroscopic hydrodynamic quantities such as the local pressure and velocity.

2.2 Membrane model and geometry

Deviations from the equilibrium shape of the vesicle increase its total elastic energy, which consists of in-plane shear and area dilation terms, and an out-of-plane bending term. The total energy is thus given by the sum of these contributions:

[TABLE]

where $\epsilon_{s}=\kappa_{s}\left(I_{1}^{2}+2I_{1}-2I_{2}\right)/12+\kappa_{\alpha}I_{2}^{2}/12$ is the surface strain energy density, which depends on the surface elastic shear modulus $\kappa_{s}$ and the area dilation modulus $\kappa_{\alpha}$ as in [26, 27], while $\kappa_{B}$ is (the out-of-plane) bending modulus of the membrane and $H$ is the mean curvature. Here $(\lambda_{1},\lambda_{2})$ are the eigenvalues of the displacement gradient tensor $D$, the invariants $I_{1}$ and $I_{2}$ are defined by $I_{1}=\lambda^{2}_{1}+\lambda^{2}_{2}-2$ and $I_{2}=\lambda^{2}_{1}\lambda^{2}_{2}-1$, and we have used the nonlinear strain energy density proposed by Skalak [45] for biological membranes, which is valid for both small and large strains. We assume that the vesicle is slightly stretchable, with a maximum of $5\%$ change in total area. For a more in-depth overview of both the LB method and the membrane model, including its relation to the microscopic structure, the numerical evolution of the deformation gradient, and the resulting membrane forces, we refer the reader to [28, 26, 27].

2.3 Simulation results

To address the question of how the applied force and relative vesicle size affect translocation, we focus on a vesicle with fixed mechanical properties as given in Table 1. The parameters were chosen to assure numerical stability, but also to yield dimensionless values governing translocation that are comparable to typical values for bilayer vesicles. Note that the dimensional values of the moduli and applied force are higher than the ones of typical bilayer vesicles [5]. However, previous research [30] has shown that although the energetics of highly stretchable containers are very different than in the case of bending energy-dominated containers, the resulting phase diagrams are similar in terms of translocation dynamics, suggesting some form of universality between the two limits. In our simulations we assume that the fluids inside and outside the capsule are Newtonian and have the same material properties, i.e. the same viscosity and density. Table 1 shows the parameters used in our lattice Boltzmann simulations, both in dimensionless and SI units. Note that these values may be manipulated by varying the viscosity and the speed of sound of the medium; this is discussed in more detail in Narváez et al. [38]. The Reynolds number of these simulations is on the order of 0.1 – 1. Moreover, since the viscosity-dominated regime of lubrication theory is applicable whenever $h/R_{p}\ll 1$ and $(h/R_{p})^{2}\text{Re}\ll 1$ [1], and in our case $h/R_{p}\approx 0.1$ and $(h/R_{p})^{2}\text{Re}=10^{-3}$ – $10^{-2}$, we are indeed in the non-inertial lubrication regime.

The surface of the vesicle is triangulated to allow for efficient calculations of the deformations, and the number of faces has been fixed to 3380, which is sufficient to capture the deformations studied (see the convergence study in Appendix C). We present our results in conventional lattice units ($\delta x=\delta t=1$) in terms of the two dimensionless parameters mentioned above; the dimensionless confinement ratio $\pi_{1}:=R_{p}/R_{c}$ and the dimensionless force $\pi_{2}:=CF/(\pi R_{p}^{3})$ where we choose area dilation as the dominant energetic contribution, hence $\pi_{2}=5F/(\pi\kappa_{\alpha}R_{p})$.

We first assume that there is a constant body force on the vesicle. The total body force $F$ is implemented by applying a constant force uniformly to each node of the vesicle. To allow for the study of the lubrication layer between the vesicle and the channel wall, the fluid grid is made to be sufficiently fine, with size $196\times 196\times 384$, and we consider two different minimal neck radii with $R_{c}=30,40$. Figs. 2(a) and (b) show a representative sequence of vesicle deformations and the associated velocity field during translocation for both closed and open constrictions with $\pi_{1}=0.75$ and $\pi_{2}=0.13$.

2.4 Reverse flow due to a blind channel hinders translocation

To showcase the difference between transportation through open and closed constrictions, we compare the time-evolution of the position of the center of mass of a vesicle in these two situations (Fig. 2(c)). The dimensionless transit time through the constriction $\tau/\tau_{0}$ (where $\tau=R_{p}/U_{\text{min}}$ is measured relative to the transit time $\tau_{0}$ in free space) is up to an order of magnitude slower for the closed constriction than for the open constriction. As we will show further at the end of this section, the transit time $\tau$ sensitively depends on $\pi_{1}$, the vesicle radius relative to that of the constriction, and on the dimensionless applied force $\pi_{2}$. To explain this dramatic difference in translocation dynamics between open and closed constrictions, we now turn our attention to the boundary layer associated with the back-flow of fluid that occurs during translocation.

To visualize the fluid flow in the narrow region between the channel wall and the vesicle, in Figs. 2(a) and (b) we plot the axial velocity (in the $z$-direction) of the fluid and vesicle for both the open and closed geometries. For the open constriction, the surrounding fluid is dragged along with the vesicle and no back-flow of fluid occurs. For the closed constriction, however, fluid incompressibility demands that a narrow layer with reverse flow emerges to allow fluid to escape the pocket as the vesicle enters the constriction. This reverse flow hinders the passage of the vesicle and increases the amount of deformation the vesicle must undergo in order to squeeze through the constriction. In Fig. 2(d) we take a closer look at the velocity profile within the constriction and find that, by increasing the applied force, the amount of fluid that leaves the closed constriction per unit time (and hence the fluid velocity in the lubrication layer) increases as well. We note that rescaling the velocity of the surrounding fluid by that of the vesicle allows us to the collapse the velocity profiles onto a universal profile.

We remark that the emergent lubrication layer is cylindrically symmetric, as one would expect from the presence of a lift force that drives the deformable vesicle from the surrounding walls [10]. In Appendix C we also demonstrate this symmetry and measure the stability of the axisymmetric configuration by showing that the vesicle recovers its centered position in the constriction after small perturbations toward the wall.

2.5 Scaling of transit time

Our simulations show that the translocation of vesicles into closed constrictions is highly dependent on the driving force $F$ and the geometry of the vesicle and constriction via the radii $R_{p}$ and $R_{c}$, as illustrated in the inset Fig. 3(a). In Fig. 3(a) we measure the minimal lubrication gap $h_{0}$ between the vesicle and the channel as a function of $R_{p}$. For the case $R_{p}/R_{c}\ll 1$ this is a linearly decreasing function, but for $R_{p}/R_{c}\approx 1$, the size of the lubrication layer $h_{0}/R_{c}$ no longer decreases linearly with $R_{p}/R_{c}$. Rather, it is set by a balance between the fluid pressure exerted on the membrane and the elastic properties of the vesicle.

To quantify how this lubrication layer impacts the transit time $\tau$, we have measured $\tau$ for various values of the body force and confinement ratio $R_{p}/R_{c}$. As expected, increasing the radius of the vesicle lowers the minimal velocity in the constriction and hence increases $\tau$. If we rescale $\tau$ by the time $\tau_{0}$ it takes for a vesicle to move its own length (i.e. $\tau_{0}=6\pi\mu R_{p}^{2}/F$ according to Stokes law), we find that all the curves collapse for moderate values of $R_{p}/R_{c}$ where $R_{p}/R_{c}<1$ (Fig. 3(b)). For $1-R_{p}/R_{c}\rightarrow 0$ we find that $\tau/\tau_{0}$ depends strongly on $R_{p}/R_{c}$, i.e. a slight increase in vesicle size drastically increases the passage time. In the inset of Fig. 3(b) we indicate that $\tau/\tau_{0}\sim(h_{0}/R_{c})^{-2}$ for gap sizes between $0.05<h_{0}/R_{c}<0.5$. Although we cannot definitively establish a power law behavior from the full 3D simulations since the gap sizes range over only a single decade, we note that this power of 2 is slightly smaller than the power of 2.5 for the rigid vesicle case that will be derived using the lubrication model in Section 3.1.

We are prevented from studying this dependency for $R_{p}/R_{c}>0.95$ since doing so would demand resolving very narrow fluid layers. Note that in the regime $R_{p}/R_{c}>1$ in which undeformed vesicles do not fit within the constriction, it has been shown previously that for $R_{p}/R_{c}$ sufficiently large the vesicles get stuck in the constriction, i.e. the scaled transit time $\tau/\tau_{0}$ diverges [30].

3 Lubrication Model

As outlined in the previous section, the transportation of vesicles through narrow, closed constrictions involves the emergence of a thin lubricating layer. Such layers are challenging to resolve in fixed-grid simulations such as the lattice Boltzmann method described above. However, we can take advantage of results from lubrication theory valid precisely in this limit, building on previous formulations for the motion of pellets through fluid-filled tubes [33, 46].

Defining $h(z)$ to be the height of the vesicle above the channel wall as a function of the axial coordinate $z$, we assume axisymmetric flows and deformations, as justified by the results of 3D simulations presented in Section 2.1. The coordinate $z$ is defined so that $z\in(Z-R_{p},Z+R_{p})$, where $Z$ is the center of mass of the vesicle. Since $h\ll R_{c}$ in the case of interest, we may invoke the standard result from lubrication theory [1] that

[TABLE]

where $U=\mathrm{d}Z/\mathrm{d}t$ is the vesicle velocity in the $z$-direction. By incompressibility, the flux $Q$ through the gap must be equal through each cross section [1] so

[TABLE]

Rewriting (5) in terms of $\partial p/\partial z$ and integrating,

[TABLE]

Setting $z=Z+R_{p}$ in the above equation (6) determines the flow rate $Q$ in terms of the pressure drop $p(Z+R_{p})-p_{0}=\Delta p$, which is a function of the applied force $F$ via $\Delta p=F/(\pi R_{p}^{2})$, neglecting the viscous drag term of order $\mu UR(R/h)$ from the balance of forces since it is dominated by the pressure scale $\mu U/\rho(R/h)^{2}$ in the lubrication limit. (For the purpose of this scaling argument, $R\sim O(R_{c})\sim O(R_{p})$ as we are interested in the limit $R_{c}\approx R_{p}$.) This results in the equation

[TABLE]

By conservation of mass, we know that the backflow $Q$ is balanced by the fluid dragged forward by the vesicle, i.e.

[TABLE]

where we have assumed that the gap is sufficiently small that $R_{p}\approx R_{c}$. To close the systems of equations for $p(z)$, $h(z)$, $Q$, and $U$, we need a constitutive law relating the height to the pressure. We assume an approximately spherical vesicle and make the ansatz that

[TABLE]

where the channel radius $\widetilde{R}_{c}(z)$ is a function of position (related to the minimum channel radius via $R_{c}=\min_{z}\widetilde{R}_{c}(z)$, $C$ is the elastic compliance of the vesicle, and $p_{0}$ is the far-field pressure which we assume vanishes. Note that (9) is a drastically simplified form of the vesicle elasticity since it does not allow for axial extension and neglects constraints such as conserved vesicle volume. A derivation of the compliance $C$ based on the vesicle material properties is given in Appendix A.

Equations (5)-(9) give a system of four equations for the the unknown flux $Q$, and velocity $U$, the pressure $p(z)$, and the height $h(z)$, and fully describe the dynamics of the forced vesicle in the constriction. (The equivalent non-dimensionalized system is presented in Appendix B.) Rather than solving (5)-(9)simultaneously, we start with the vesicle position $Z=Z(t)$ and initial guesses $p^{0}(z)$, $h^{0}(z)$, $Q^{0}$, and $U^{0}$ and use an iterative method as follows:

[TABLE]

iterating until the system converges to fixed points $p^{\infty}(z)$, $h^{\infty}(z)$, $Q^{\infty}$, and $U^{\infty}$. The vesicle position is then updated using $\mathrm{d}Z/\mathrm{d}t=U^{\infty}$ and the process is repeated using the new position $Z(t+\Delta t)$. The converged values $p^{\infty}(z)$, $h^{\infty}(z)$, $Q^{\infty}$, and $U^{\infty}$ at the previous step are used as initial guesses for the next iteration. This ensures that the iterates converge reliably to a solution provided that the timestep $\Delta t$ is sufficiently small. We have found that the rate of convergence can be accelerated using Steffensen’s method [7].

3.1 Scaling of transit time at constant applied force

As shown in Fig. 4(a), there is good agreement between the trajectories computed by the lubrication model (LM) and 3D lattice Boltzmann (LB) simulations. Here, we have used the channel geometry illustrated in Fig. 4(b). Because of its analytical simplicity in comparison to directly solving the full 3D fluid-structure interaction problem, the lubrication model allows for a thorough exploration of the phase space associated with $\pi_{1}=R_{p}/R_{c}$ and $\pi_{2}=CF/(\pi R_{p}^{3})$ (Fig. 4(c)). (See Table 2 for the parameters used in the lubrication model.) For an elastic vesicle with compliance $C=5\times 10^{-9}\,\text{m}/\text{Pa}$, simulation of the lubrication model reveals a plateau in the transit times as $\pi_{1}\to 1$ (Fig. 4(d)). The reason for this plateau is that, when $\pi_{1}\sim 1$, the deflection term $Cp(z)$ dominates in the equation (9) for the height. It is precisely $\pi_{2}$ that controls the magnitude of this deflection term: we find the dimensionless minimum height $h_{0}/R_{c}\approx\pi_{2}/2$.

In the case $\pi_{2}/2\ll 1-\pi_{1}$ in which elastic deformations are negligible, the absence of deformations yields much narrower lubrication layers. We find a scaling relation $\tau/\tau_{0}\sim(h_{0}/R_{c})^{-5/2}$ in the small $h_{0}$ limit. (Recall that $\tau_{0}$ is the time required for the forced vesicle to move a distance equal to its own length in free space and that $h_{0}$ is the minimal gap spacing, which in the inelastic case is related to $\pi_{1}$ by $h_{0}/R_{c}\approx 1-\pi_{1}$.) To understand this scaling, we start from (7) and (8), which yield

[TABLE]

To estimate the integral $\int_{Z-R_{p}}^{Z+R_{p}}1/h^{n}(z)\mathrm{d}s$ for arbitrary $n$, we use the formula (9) for the height function, neglecting the elastic term, so that

[TABLE]

where making the approximation $\sqrt{1+s^{2}}\approx 1+s^{2}/2$ is justified since the integrand becomes singular near $s=0$ as $h_{0}\to 0$. It follows that the timescale $\tau\sim R/U\sim\mu R^{7/2}(h_{0}^{-3/2}+Rh_{0}^{-5/2})/F\sim(\mu R^{9/2}/F)h_{0}^{-5/2}$, since when $h_{0}\ll 1$, $h_{0}^{-3/2}\ll h_{0}^{-5/2}$. Since $R_{c}$ is constant and $\tau_{0}$ is constant in the limit $R_{p}\to R_{c}$, we have $\tau/\tau_{0}\sim(h_{0}/R_{c})^{-5/2}$.

Neglecting the elastic term is appropriate since for $s\approx 0$ and $R_{c}\approx R_{p}=:R$,

[TABLE]

where we have used $p(s)\approx F/(2\pi R^{2})$ at $s\approx 0$ in the second line and the assumption $\pi_{2}/2\ll 1-\pi_{1}$ in the final line.

In contrast to the scaling of the transit time with the minimal gap size $\tau\sim h_{0}^{-5/2}$ derived in the inelastic case, the scaling observed in the lattice Boltzmann simulations of Section 2.5 is $\tau\sim h_{0}^{-2}$. This is because those simulations take place in a regime that is neither inelastic nor dominated by membrane compliance, but rather in a transition regime, corresponding to the bend in the graph of Fig. 4(d), that gives rise to an apparent $h_{0}^{-2}$ scaling over a narrow range of $h_{0}$.

3.2 Active translocation by molecular motors

Having seen how the simplicity of the lubrication model allows us to understand the pressure-driven passage of deformable vesicles through narrow intracellular channels relatively easily, we now turn to complement our analysis by accounting for endosomal transport driven by the dynamics of molecular motors that have their own kinetics of binding and unbinding associated with force production.

Although the molecular motor dynein walking along microtubules plays an important role in transporting vesicles through dendritic branches, there is strong experimental evidence that the actin-myosin cytoskeleton dominates transport into the spines [11]. Therefore, to model motor activity, we use a force-velocity relation for myosin binding to actin filaments in muscle based on the classic work of Huxley [22] and more recent work by Lacker and Peskin [31, 21]. Although the details of molecular motor function undoubtedly differ between bidirectional intracellular transportation and muscle, the key model assumptions do not depend sensitively on the details of the individual motors, in that we assume (i) there are two independent motor species that are identical aside from the direction in which they push, and (ii) motors are fixed within the cortex and do not diffuse or travel with the vesicle, as allowed in some other models such as [50]. We also refer the interested reader to the related descriptions of molecular motor assemblies in [24, 9]. Defining $\phi(z)$ to be the population density of motors attached to a vesicle and $\theta_{[a,b)}$ as the fraction of attached motors with displacements satisfying $a\leq z<b$, it follows that

[TABLE]

The total fraction $\theta$ of attached motors is therefore

[TABLE]

and the total force exerted on the vesicle is

[TABLE]

where $n_{0}$ is the total number of motors and $f(z)$ is the functional form of the force exerted by an individual motor, which will be specified later on. Assume rate constants $\alpha$ and $\beta$ of motor attachment and detachment, respectively, and that connections are always formed at a displacement $z=A$, so that the motors tend to push the vesicle in the $-z$ direction with velocity $U<0$. (We will later discuss the case of motion in the $+z$ direction.) Now, considering only those attachments with displacement $z_{0}<z<A$ for fixed $z_{0}$, at steady-state it must follow that

[TABLE]

representing a balance between formation of new attachments at a rate $\alpha(1-\theta)$, detachment at a rate $\beta\int_{z_{0}}^{A}\phi(z)\mathrm{d}z$, and motion by advective flux at a rate $-U\phi(z_{0})$. From (20), one can derive the equation

[TABLE]

as well as the boundary condition $\alpha(1-\theta)=-U\phi(A)$, from which it follows that

[TABLE]

Using the form $f(z)=-p_{1}\left(\exp\left(\gamma z\right)-1\right)$ for the force exerted by an individual motor, we substitute into (19) the above result and carry out the integration in (19) to obtain the force-velocity relation for the ensemble of motors

[TABLE]

which can be rewritten in the familiar Hill-like form [21, 20] as

[TABLE]

where

[TABLE]

The case of motion in the $+z$ direction is posed as an exercise in [21], and we work it out here in detail. For $U>0$, the equation corresponding to (20) is

[TABLE]

If we assume that the rate $\beta$ of detachment is constant for all $z>A$, this results in the population density

[TABLE]

Using this population density and performing the integration in (19) results in the force

[TABLE]

which happens to be identical to (23), so that there is a single smooth function representing motion in both the $+z$ and $-z$ directions. However, as motors cannot be stretched to arbitrarily large displacement without breaking, we introduce a maximum displacement $B>A$ at which motors must detach. The fraction of attached motors then becomes

[TABLE]

and solving for the population density from (21) yields

[TABLE]

where the normalization constant $c=1-\exp\left(\beta(B-A)/U\right)$. Performing the integral

[TABLE]

with $f(z)$ as before results in the force-velocity curve

[TABLE]

Note that (30) is recovered in the limit $B\to\infty$. When $B$ is finite, the force-velocity curves for motion in the $+z$ and $-z$ directions are no longer described by a single equation, but it can be shown that the resulting piecewise formula is smooth at the transition $U=0$ [21]. All together then, for motors forming connections at $z=A$ and attempting to move the vesicle in the minus $z$ direction,

[TABLE]

where we have assumed that motors connect at a fixed displacement $z=A$. All motors attach at a fixed rate $\alpha$, detach at fixed rate $\beta$, and are stretched or contracted according to the vesicle velocity. All motors detach upon extension to $z=B$. See Fig. 5(b) for the form of the resulting force velocity curve using a representative value of the nondimensional quantity $\gamma(B-A)$.

We have introduced the notation $F_{A}$ to denote the force generated by the motor species that forms connections at $z=A$ considered thus far. A competing species that forms connections at $z=-A$ and disconnects at $z=-B$ behaves in the opposite way, with motion in the $+z$ direction corresponding to shortening, and motion in the $-z$ direction corresponding to lengthening. The equations in this case are

[TABLE]

We assume that the two competing species of motors are coupled only through the resulting velocity of motion. This assumption allows us to compute the net force on a vesicle simply by adding the forces from each species. Define $\phi_{1}$ to be the fraction of motors that form connections at $z=-A$ out of the $n_{0}$ total motors, so that the motor species that forms connections at $z=A$ makes up the fraction $\phi_{2}=1-\phi_{1}$ of the total. This results in a force $F=\phi_{1}F_{-A}+\phi_{2}F_{A}$, with the resultant force-velocity curve shown in Fig. 5(c).

Because of the kinetics of attachment and detachment, including two opposing species of molecular motors at equal concentrations leads to multistable dynamics in which vesicles can be trafficked in either direction or remain stationary. For example, setting $\phi_{1}=\phi_{2}=0.5$ results in a multistable behavior in which the motors can generate forward, reverse, or net zero motion (Fig. 5). The steady-state velocities are precisely those at which the force production curve intersects the viscous drag curve. To simplify the bifurcation analysis, we approximate the viscous drag curve by a line through the origin. The slope of the resulting viscous drag line depends on the degree of confinement; in free space, this is the familiar Stokes law $F=6\pi\mu R_{p}U$. The existence of multiple steady-states is illustrated in Fig. 5(b) by plotting the difference between the motor force-velocity curve and the viscous drag $F=\zeta U$ for a representative friction coefficient $\zeta$. Visual inspection reveals the existence of 5 steady states, 3 of which are stable. The stability of the rightmost steady state can be assessed in terms of the sign of the derivative at steady-state. For negative derivatives, if the velocity is increased slightly, the drag force exceeds the motor force and the vesicle slows down; if the velocity is decreased slightly, the motor force exceeds the drag and the vesicle speeds up. Therefore, the rightmost steady-state is stable. Likewise, positive derivatives indicate unstable steady-states.

As the fraction $\phi_{1}$ is increased or decreased, there is a sequence of bifurcations at which the multiple fixed points coalesce onto the single stable solution observed in the case of having only one motor species (Fig. 5(c)). For small $\phi_{1}$, there is a single stable velocity and it is in the $-z$ direction. As the fraction increases, there is a saddle-node bifurcation and a stable velocity in the $+z$ direction emerges in addition to the stable negative solution. When $\phi_{1}\sim 0.5$, a second saddle-node bifurcation gives rise to another stable solution that crosses from small negative to small positive velocities. In Fig. 5(d) we plot a complementary bifurcation diagram showing the dependence of the steady-state velocities on the height-dependent friction coefficient $\zeta$, under the assumption of equal numbers of both species of motors. For small values of $\zeta$, there are five steady-state velocities, of which three are stable. As the drag increases, two saddle-node bifurcations annihilate four steady states, leaving only the stationary solution.

Note that the approximation of linear viscous drag is used only to study the qualitative behavior; in our actual simulations, we always solve for the fluid drag that emerges from the coupling between fluid, elastic, and motor dynamics.

Multistability in mixtures of motors of opposing polarities has been previously observed in both experimental and theoretical studies [18]. Gilboa et al. [17] performed motility assays on an experimental actomyosin system with actin tracks of alternating polarities and observed bidirectional motion. Muller et al. [37] studied theoretically the effect of antagonistic motor species exerting forces in opposite directions, which they called tug-of-war, and observed multistability upon simulating a master equation for the distribution of bound states. Hexner and Kafri [19] also considered a similar tug-of-war scenario and analyzed the possible phase space, including bifurcation diagrams consistent with those of Fig. 5. More recent work includes the effect of thermal fluctuations on switching between states [35] and tug-of-war in the context of searches for a random target [39]. The present study builds on this previous work by investigating quantitatively the subtle interplay between the multistable dynamics of mixed motor systems and the effects of hydrodynamics associated with the confined geometry in which the cargo moves. Our results demonstrate that the steady state velocities depend not only on the fraction of motor species present, but also on the effective Stokes law. Since the viscous drag depends on the geometry of the constriction through which the cargo passes, this sets the force that must be generated by the aggregate of mixed-polarity motors. As the channel height varies, the velocity can jump abruptly from one branch of solutions to another, as we now demonstrate in the next section.

4 Vesicle trafficking into dendritic spines

As mentioned in the introduction, the dendritic spine is an example of a heterogeneous environment in which various biopolymers interact with different species of molecular motors that traverse filaments in opposite directions and species that are able to reverse their direction of motion. This diversity helps explain the richness of behaviors seen in vesicle transit through spines, including (i) one-directional transit through the spine (ii) bidirectional transit up to the spine neck and back, and (iii) corking in the spine neck.

To model the effect of a crowded intracellular environment, we use an elevated fluid viscosity $\mu=0.12\,\text{Pa}\cdot\text{s}$, which is 100x the viscosity of water. We also include a proxy for thermal fluctuations in the following manner: instead of taking the converged velocity $U^{\infty}_{j}$ at timestep $j$ as the initial guess for fixed point iteration at the next timestep $U^{0}_{j+1}$, we add noise by setting $U^{0}_{j+1}=U^{\infty}_{j}+\eta\xi$, where $\xi$ is a random number drawn from a zero-mean unit-variance normal distribution. This increases the likelihood that the solutions of fixed point iteration at subsequent time steps will lie on different solution branches.

Upon including the effects of competing molecular motor species our model is able to recapitulate at a qualitative level many of the phenomena observed in vesicular trafficking into dendritic spines. In these simulations $(\pi_{2}/2)/(1-\pi_{1})$ ranges from $1/6$ – $1/4$, so that according to the criterion of Section 3.1 vesicle elasticity plays only a moderate role and the qualitative analysis of Section 3.2 is applicable. Our simulations lead to bidirectional motion and transit times on the order of $10$ – $100$ seconds, in line with the experimental observations [11]. When the fraction $\phi_{1}=0.57$ of upward-directed motors is used, $R_{p}=0.96$ µm and $R_{c}=1.22$ µm, with a total motor force of 50 pN and noise magnitude $\eta=1.35\times 10^{-1}$ µm/s, the vesicle moves in a processive fashion up the spine (Fig. 6(a)), consistent with the single steady solution in the limit $\phi_{1}\to 1$. Given myosin’s stall force of around 3 pN [34], the motor forces of 50 – 200 pN used here correspond to having a few dozen motor molecules. When equal fractions $\phi_{1}=\phi_{2}=0.5$ of upward and downward-directed motors are used with the parameters otherwise fixed, the vesicle becomes trapped in the spine neck (Fig. 6(b)). This is reminiscent of experimental observations of vesicles trapped in spines and hypotheses that such corking of vesicles in spine necks may serve as a chemical insulator [40, 47], as has been studied theoretically in terms of anomalous diffusion along spiny dendrites (see [6] and references therein). When the fraction $\phi_{1}=0.46$ of upward-directed motors is used, either positive or negative velocities can occur at the same position in the channel depending on the random initial guess. In one representative simulation, the vesicle makes it through the spine neck on its first attempt but not on its second attempt, during which it turns around and must try again (Fig. 6(c)). This last simulation is done using a different set of parameters, with the minimum channel radius $R_{c}=2.15$ µm, vesicle radius $R_{p}=1.5$ µm, forcing $F_{0}=200$ pN, noise magnitude $\eta=1.35\times 10^{-2}$ µm/s, and $\pi_{5}=0.02$. We find analogous behavior using equal but opposite fractions with $\phi_{1}=0.54$, but with the vesicle spending most of its time sequestered in the spine (Fig. 6(d)).

Our simulations are subject to the boundary conditions that only non-negative velocities are allowed at the base of the channel, and only non-positive velocities are allowed at the end of the channel. Therefore the vesicle waits at the base until it acquires a large enough random kick and begins to move in the $+z$ direction into the spine. Similarly, the vesicle stops at the end of the spine and waits until it receives a sufficiently large random kick to be brought in the opposite direction by the same assembly of motors.

The observation of bidirectional transport in our simulations (Fig. 6(c)) is supported by experimental observations in [11] and previous studies of molecular motors present in dendritic spines. Restricting attention to actomyosin-based vesicle trafficking, which is thought to be the dominant mechanism for vesicle entry into spines [11], we are aware of two plausible mechanisms for bidirectional transport. First, myosin types II, V, and VI are reported to be present in spines [25], and whereas all other known myosin species walk toward the barbed (+) end of actin, myosin VI walks toward the pointed (-) end [48]. Second, actin of mixed polarities has been observed in spine necks [25], so that a single-directional myosin motor could move in either direction based on the actin bundle to which it is attached.

5 Discussion

By performing fully resolved 3D simulations using the lattice Boltzmann method, we have been able to study the translocation of an impermeable elastic vesicle into a closed tube filled with incompressible fluid. If the vesicle radius is nearly equal to the radius of the channel neck, the vesicle deforms into a teardrop-like shape and is separated from the channel walls by a thin lubrication layer. Because the fluid trapped in the pocket must escape, this fluid layer includes both forward and reverse flow, unlike the case of an open channel in which only forward flow is observed. This reverse flow decreases the size of the effective lubrication layer for the closed channel and increases the amount of force required for translocation at a given speed.

To capture the essential behavior, we have developed a lubrication model that involves only two nondimensional parameters involving the gap size and applied forcing and allows us to efficiently explore a large parameter range. We have found that the dependence of the transit time $\tau$ on minimum lubrication layer width $h_{0}$ scales as $\tau\sim h_{0}^{-5/2}$ in the inelastic regime, and these results have been validated by a combination of simulation and scaling analysis. We now make a brief digression to compare these results to those obtained in previous work on the forced transport of elastic containers through constrictions [30]. In that work, which focused on the case of large vesicles and identified a finite force under which the vesicle could not pass, the elastic properties of the vesicle dominated over fluid effects such as the formation of a lubrication layer. Here, the limit of interest is the one in which the vesicle size is nearly equal to the constriction size and lubrication effects play a crucial role, necessitating a more realistic model of fluid drag such as the the lubrication model we have used. Moreover, whereas in [30] the vesicle was assumed to pass through a series of known shapes, allowing for a computation of the energy and therefore the motor force required as functions of position through the channel, here we use a simple model for the vesicle shape changes (i.e. linear compliance). The lubrication assumption is invoked at all points through the vesicle transit, and the coupling of the shape to the pressure via the linear compliance results in a vesicle that slows down according to a power law but never stops moving through the channel. Combining this lubrication model with a more comprehensive series of allowed shapes would enable capturing the relevant fluid dynamics as well as the elastic barriers to translocation.

Developing more realistic models of the vesicle’s elastic response at large deformations that capture nonlinear stiffening effects would improve on the simple linear compliance model. Removing the assumption of small deformations would also allow us to explore the case $\pi_{1}>1$, which is of significant biological interest. As noted above, we believe the effects investigated here in the small deformations limit are applicable to the case of large deformations, and we expect the impeding effect of a closed end on the vesicle motion to become even more dramatic because of the narrower lubrication layer involved. We have not yet been able to investigate the case of large deformations since, in addition to the modeling challenges of accurately capturing vesicle shape and elasticity, it is difficult to resolve the smaller lubrication layers involved because of the maximum grid size allowed in the lattice Boltzmann simulations. Efforts are underway to overcome these challenges. In particular, some of us have begun developing a numerical method for fluid-structure interaction that combines a direct Navier-Stokes solver for the bulk flow with lubrication theory in regions with thin fluid layers [15].

We have added an additional level of realism to our model by including a mechanistic force-velocity relation for the forces exerted by molecular motors. The motor model we use includes two species with opposite polarities and exhibits multistable dynamics, consistent with previous studies in the literature of mixtures of molecular motors. Our work demonstrates that this multistability can interact with the channel geometry in fascinating ways, and our model is able to reproduce several of the behaviors observed experimentally in spines. The number and direction of steady state velocities can be controlled by varying the channel radius and the fraction of motors pushing in either direction, and these control mechanisms may be used by the cell to control dendritic spine maintenance, growth, and atrophy. As mentioned above, vesicle elasticity does not play a critical role in Section 3.2 and therefore is not an essential ingredient for the multistable dynamics. Performing an in-depth parameter space investigation in the future would make it possible to see how the behavior changes in regions where fluid flow and molecular motors are strongly coupled to changes in vesicle shape.

The present model makes several simplifying assumptions that must be refined for more specific quantitative models, e.g. biophysically detailed studies of vesicle trafficking into dendritic spines. For instance, we have left out the effect of fluctuations caused by Brownian motion and motor kinetics. The formulation of the lubrication model could be revisited to include thermal fluctuations, and the continuum theory of molecular motors used here could be replaced by a kinetic model to explicitly include the stochastic effects of motors coming on and off. Given the sub-micron scale of dendritic spine necks, we expect stochasticity to be an important ingredient in more quantitative studies of vesicle trafficking. Intracellular environments such as dendritic spines are crowded environments filled with various biomolecules, and the assumption of a surrounding Newtonian fluid could be improved by accounting for these non-Newtonian effects. Finally, the presence of a strong reverse flow through the lubricating layer depended crucially on the assumption of impermeability to water of both the spine wall and vesicle membrane. This is consistent with experimental evidence that mammalian neurons do not contain water-permeable membrane channels called aquaporins [2], but experimental investigation of the permeability of recycling endosomes is necessary to validate this assumption.

6 Acknowledgments

We thank the anonymous reviewers for their constructive feedback and assistance in improving this manuscript. The authors wish to acknowledge financial support under National Science Foundation grant DMS-1502851 (TGF), the Netherlands Organization for Scientific Research (NWO-FOM) within the program ”Barriers in the Brain: the Molecular Physics of Learning and Memory” (No. FOM-E1012M), and the NSF grant DMR 14-20570 through the Harvard MRSEC (CR, LM), and the MacArthur Foundation (LM).

Appendix A Determination of Compliance $C$

The elasticity of the vesicle in our lubrication model is captured by a single-parameter, the compliance, which determines the deflection of the vesicle in response to the pressure. To be able to use the model to simulate vesicles with known material properties, we need to determine the appropriate compliance. To do so, we consider a spheroidal neo-Hookean surface parameterized by coordinates $\mathbf{q}=(q_{1},q_{2})$ with the deformed configuration $\mathbf{X}=\mathbf{X}(\mathbf{q})$ and reference configuration $\mathbf{Z}=\mathbf{Z}(\mathbf{q})$ in $\mathbb{R}^{3}$. Following the formulation described in [14, 49], the elastic energy $E$ of such a surface with bulk modulus $\kappa_{\alpha}$ and shear modulus $\kappa_{s}$ is given by

[TABLE]

where $\lambda_{1}$ and $\lambda_{2}$ are the principal stretch ratios, i.e. $\lambda_{1}^{2}$ and $\lambda_{2}^{2}$ are the eigenvalues of the $2\times 2$ matrix $G_{0}^{-1}G$, where $G_{ij}=(\partial X_{k}/\partial q_{i})(\partial X_{k}/\partial q_{j})$ and $(G_{0})_{ij}=(\partial Z_{k}/\partial q_{i})(\partial Z_{k}/\partial q_{j})$, using the Einstein convention for summation over repeated indices. Let the reference configuration $\mathbf{Z}=(z_{1},z_{2},z_{3})$ be given by the sphere with radius $R$ so that

[TABLE]

where $q_{1}\in[0,\pi)$ and $q_{2}\in[0,2\pi)$ are the usual polar and azimuthal coordinates, respectively. Now, consider oblate spheroidal deformations $\mathbf{X}=(x_{1},x_{2},x_{3})$ parametrized by

[TABLE]

with $\sigma>1$. (Taking $\sigma<1$ would give a prolate spheroid.) Note that the axes are scaled so that the total volume enclosed by the surface is independent of $\sigma$. Direct computation gives

[TABLE]

Neglecting the poles which are a set of measure zero, we have therefore

[TABLE]

which implies that $\lambda_{1}^{2}\lambda_{2}^{2}=\det(G_{0}^{-1}G)=\sigma^{2}\cos^{2}q_{1}+\sigma^{-1}\sin^{2}q_{1}=\sigma^{-1}\left(1+(\sigma^{3}-1)\cos^{2}q_{1}\right)$. It follows that the elastic energy is

[TABLE]

where we have neglected the shear energy and used the fact that $\det{G_{0}}^{1/2}$ is the area element in reference coordinates.

Letting $\sigma=1+\epsilon$ for $\epsilon\ll 1$, $\sigma^{-1/2}\approx 1-\epsilon/2$ and $\sigma^{3}-1\approx 3\epsilon$. Therefore

[TABLE]

where we have used $\det{G_{0}}^{1/2}=a^{2}\sin q_{1}$ in the final equality. Evaluating the remaining integral is straightforward and results in a value of $16\pi/5$. Thus,

[TABLE]

so that the effective spring constant is $4\pi\kappa_{\alpha}/5$. Since the force $F$ on the oblate ellipsoid satisfies $F=(4\pi\kappa_{\alpha}/5)(a\epsilon)$ whereas the compliance is related to the pressure by $Cp=a\epsilon$, with $p\approx F/(4\pi R_{p}^{2})$,

[TABLE]

Taking representative values of $R_{p}=1$ µm and $\kappa_{\alpha}=1\times 10^{-3}\,\text{N}/\text{m}$ yields an estimate for the compliance $C\approx 5\times 10^{-9}\,\text{m}/\text{Pa}$. Note that we have kept only terms involving $\kappa_{\alpha}$ and have neglected the additional elastic moduli $\kappa_{B}$ and $\kappa_{s}$. The good agreement between the compliance derived in this appendix and the compliance estimated from the full 3D simulation justifies our assumption that the membrane’s elastic properties are essentially determined by the single modulus $\kappa_{\alpha}$ in the regime of interest.

Appendix B Nondimensionalization

Here, we present a nondimensionalized system of equations for the lubrication model. As stated in Section 2, there are two nondimensional groups $\pi_{1}$ and $\pi_{2}$ that govern the fluid-structure interaction dynamics: $\pi_{1}=R_{p}/R_{c}$, where $R_{c}$ is the radius at the narrowest point in the channel, and $\pi_{2}=CF/(\pi R_{p}^{3})$, where $C$ is the compliance and $F$ is the strength of forcing. Note that in the context of the 3D membrane model we may write $\pi_{2}=5F/(\pi\kappa_{\alpha}R_{p})$ according to (41). Of course, a complete investigation of the 3D problem would require additional dimensionless groups, e.g. those involving the length of the channel [30] and the additional elastic moduli $\kappa_{B}$ and $\kappa_{s}$. However, as mentioned previously, in the lubrication limit the essential behavior is captured by representing the channel geometry with the single parameter $R_{c}$ and the vesicle elasticity with a single parameter $C$ based on the area dilation modulus $\kappa_{\alpha}$.

We introduce the nondimensional variables $\widetilde{z}=z/R_{c}$, $\widetilde{h}=h/R_{c}$, $\widetilde{U}=6\pi R_{p}\mu U/F$, $\widetilde{Q}=6\pi R_{p}\mu U/(F(\pi R_{c})(R_{c}-R_{p}))$, and $\widetilde{p}=\pi R_{p}R_{c}p/F$. Rewriting the formulation (6)–(9) for the model under a fixed external force in nondimensional form, we have:

[TABLE]

Several interesting features of this model become readily apparent when written in the nondimensionalized form above. On the one hand, in the limit $\pi_{1}\to 1$, the product $\widetilde{Q}(1-\pi_{1})$ tends to zero so that one might erroneously assume it can be neglected. In fact, it is precisely in this limit that $\widetilde{h}\to 0$, so that the integrand $\widetilde{h}^{-3}$ become singular and the term $\widetilde{Q}(1-\pi_{1})\int_{0}^{\widetilde{z}}\widetilde{h}^{-3}(s)\mathrm{d}s$ must be retained. Another property of this model is that the dimensionless variable $\widetilde{U}$ is not constant, as would be the case in free space according to Stokes’ law. In confined geometries such as the channel considered here, its value depends on the nondimensional spacing $\widetilde{h}$.

Next, we consider how incorporating molecular motors into the lubrication model as in Section 3.2 affects the above nondimensionalized system. As described in [21], the motor model adds four nondimensional groups $\pi_{i},\,i=3,\,4,\,5,\,6$, with $\pi_{3}=\alpha/\beta$ the ratio of attachment and detachment rates, $\pi_{4}=\gamma A$ the nondimensional attachment position, and $\pi_{5}=\gamma(B-A)$ reflecting the maximum negative displacement of a motor. There is a fourth additional nondimensional group $\pi_{6}=(F/6\pi R_{p}\mu)/(\beta/\gamma)$ that gives the ratio of velocity scales between translocation and motor adhesion dynamics. Note that in the context of the motor model, the forcing constant to be used in $\pi_{2}$ is simply the stall force $F_{0}=\alpha p_{1}n_{0}/(\alpha+\beta)(\exp(\gamma A)-1)$.

In terms of these nondimensional units, the motor model may be written

[TABLE]

where $\widetilde{F}_{A}=F_{A}/F_{0}$. The form of $F_{-A}$ is analogous:

[TABLE]

Setting $\widetilde{F}=\phi_{1}\widetilde{F}_{A}+\phi_{2}\widetilde{F}_{-A}$, the motors are then incorporated in the lubrication model by replacing (43) by

[TABLE]

For our simulations in Section 4, we have used $\pi_{3}=1$, $\pi_{4}=4.7$, $\pi_{5}$ in the range $\pi_{5}=0.02$ – $0.1$, and $\pi_{6}=10$ – $18$. In terms of the key dimensional parameters, using a stall force of $F_{0}=50$–$200$ pN results in transit times on the order of $10$ – $100$ seconds for vesicles of radius $R_{p}=1$ – $2$ µm. This is consistent with the experimentally observed transit times of approximately $40$ – $60$ seconds [11].

Appendix C Lattice Bolzmann simulations

Here we discuss the convergence study that we have performed to validate the 3D simulation results obtained using the lattice Boltzmann method and show that the axisymmetric assumption of the lubrication model is supported by the 3D simulations. All the simulations in this paper have been performed on a $128\times 128\times 384$ grid. To ensure this resolution is sufficient to capture the narrow lubrication layer, we perform simulations using identical parameters on $64\times 64\times 192$ and $128\times 128\times 384$ grids ($\pi_{1}=0.62$ and $\pi_{2}=0.16$, and mechanical properties as in Table 1). In Fig. 7(a) we plot the velocity of the vesicle as a function of its position and find that the velocity profiles are nearly identical for both grids. This ensures that a system size of $128\times 128\times 384$ is sufficient to capture the translocation dynamics, at least for these system parameters. Similarly we have performed a convergence study on the mesh-size of the vesicle in Fig. 7(b) by calculating the velocity as a function of position for vesicles discretized using 8000, 9680 and 11520 mesh points. In this figure we show that all three mesh resolutions yield essentially the same trajectory. In all the simulations presented in this paper we have chosen the number of mesh points to be on the order of $480R_{c}$, which has been shown to give converged results in the parameter regime used [26, 28].

The flow-field around the vesicle is axisymmetric during translocation, as illustrated by the radially-symmetric velocity profile of Fig. 7(c), which shows a color-map of the $z$-velocity in the $x$-$y$ plane halfway through the constriction at $z=180$ (using parameters as in Fig. 2(a)). We have observed that this axisymmetric configuration is stable in the sense that axisymmetry is recovered after slight perturbations in the radial position of the vesicle.

In Fig. 7(d) we compare the time evolution of the vesicle for identical dimensionless groups $\pi_{1}=0.73$ and $\pi_{2}=0.27$, where the black curve corresponds to the moduli as in Table 1 and the red curve has an applied force and mechanical moduli that are one order of magnitude smaller ($F=0.5$, $\kappa_{\alpha}=0.1$, $\kappa_{s}=0.0015$ and $\kappa_{B}=0.0018$). Note that the rescaled velocity within the constriction is identical for both cases, ensuring that these dimensionless groups are the relevant groups to describe the translocation of the elastic vesicle through the constriction.