{'fraction_non_alphanumeric': 0.03505662462020072, 'fraction_numerical': 0.014892735475554737, 'mean_word_length': 4.498671054296924, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 8, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Dynamics of Vesicles Driven Into Closed Constrictions by Molecular Motors

We study the dynamics of a model of membrane vesicle transport into dendritic spines, which are bulbous intracellular compartments in neurons driven by molecular motors. We reduce the lubrication model proposed in [Fai et al, Active elastohydrodynamics of vesicles in narrow, blind constrictions. Phys. Rev. Fluids, 2 (2017), 113601] to a fast-slow system, yielding an analytically and numerically tractable equation that is equivalent to the original model in the overdamped limit. The key parameters of the model are the ratio of motors that prefer to push toward the head of the dendritic spine to the ratio of motors that prefer to push in the opposite direction. We perform a numerical bifurcation analysis in these parameters and find that steady-state vesicle velocities appear and disappear through a number of saddle-node bifurcations. This allows us to identify the region of parameter space in which multiple stable velocities exist. We show by direct calculations that for sufficiently close vesicle-to-spine diameter ratios, there can only be unidirectional motion. The critical vesicle-to-spine diameter ratio predicted by our analysis, at which there is a transition from unidirectional to bidirectional motion, is consistent with experimental observations of vesicle trajectories in the literature.


Introduction
Intracellular transport is required by cells for many essential purposes, e.g. to deliver proteins to their functional locations and to maintain concentration gradients. The relatively long lengths of axons and dendrites means that diffusion is too slow to serve as a practical transport mechanism in neurons, making directed transport particularly important. Not only size, but also geometry, may play an important role. Bulbs known as dendritic spines cover dendritic trees of pyramidal neurons and are ubiquitous in the cortex: up to 85% of the cortex consists of pyramidal neurons [14], and up to 85% of excitatory inputs terminate on dendritic spines [34]. Spines and their distinctive shapes-with narrow necks of significant curvature separating the spines from the dendritic branch-are likely to have a significant influence on transport.
Dendritic spines protrude 0.5-1µm [38] from the dendritic shaft and often appear mushroom-like in shape. Spines exhibit a thin 50-100nm [21] neck that opens into a head that is up to several hundred nanometers greater in diameter [29]. A key aspect of dendritic spine remodeling and maintenance includes the transport of surface proteins on vesicles that are squeezed through the narrow neck into the spine head [10,35,40,16]. Recent experiments have shown that movement is not always unidirectional (translocation), but includes no movement (corking), and rejection [35,40]. The mechanisms underlying these changes in direction are not well understood.
There are generally two approaches taken by theoretical studies of molecular motordriven intracellular transport. In the first, the motor complex is assumed to have switching rates determined by some underlying Markov process, and consequences of varying switching rates are explored, such as the distribution of motor complex velocities [30,26], and mean first passage times to transport targets on dendritic morphologies [6,32,31,33,7]. In the second, individual motors are modeled as part of a greater complex or population, which generatively produces bidirectional motion despite the assumption of symmetry [23,20,2,36]. However, the effect of constrictions on motor dynamics has typically been held constant or neglected.
Constricted motion is of general interest in engineering and biology with related problems appearing in manufactured elastic capsules [12,13], hydrogels [28], the movement of living cells [4,9,18], and axonal transport [43,25,39]. In terms of theory, the effect of constrictions on molecular dynamics is a very recent topic of interest. In Fai et al. 2017, it was found that constrictions enrich the dynamics of motor populations, allowing motors to switch between multidirectional and unidirectional motion [16]. The goal of this paper is to thoroughly explore the dynamics of the model derived in [16] by classifying the bifurcations of the underlying ODE model. The paper is organized as follows. In Section 2, we briefly review the derivation of the lubrication model and its nondimensionalization, and we briefly discuss two equivalent versions, one of which (the "slow" subsystem) is the model considered in [16], and other of which (the "fast" subsystem) is the model we explore in depth in this paper. In Section 3, we numerically establish the existence and robustness of multistability through a bifurcation analysis of the "fast" subsystem. Some of these numerical results are corroborated in Section 4, where we analytically establish the existence and stability of particular velocities as a function of key parameters. We conclude the paper with a discussion that includes estimates of realistic parameter regimes of this model, and the resulting behaviors predicted for different types of dendritic spines. All code and data and documentation for reproducing figures in this paper are available on our github repository at https://github.com/youngmp/park_fai_2020

Lubrication Model
The lubrication model we consider in this paper was originally derived in [16], where the authors consider an idealized model of vesicle translation through the neck of a dendritic spine. In contrast to previous studies of transport through constrictions in periodic or unbounded tubes, the authors close off one end to model transport into spine-like geometries. Objects entering so-called closed constrictions are surrounded by fluid backflow, which introduces a large velocity gradient, allowing for a simplification of the Navier-Stokes equations through lubrication theory [1]. Variables and parameters are defined as follows. The center of mass of the vesicle is Z and its radius is R p . The function p(z) is the pressure exerted on the vesicle at position z ∈ [Z − R p , Z + R p ]. The flux (Q) equations follow from the observation that the amount of fluid carried by the vesicle must match the amount of fluid escaping the closed constriction, i.e., at position Z, the cross-sectional area of the vesicle is balanced by the fluid backflow determined by the distance h between the vesicle and spine wall. This results in the following reduced axisymmetric model of vesicle trafficking: Here, F (U ) is the net force from the molecular motors and is a key feature that contributes to the tug-of-war dynamics of the model. While the viscocity of water is on the order of µ = 0.69mPa S at a body temperature of 37 degrees Celcius, the intracellular environment is densely packed with proteins, filaments, and organelles [35,41,19], which may increase the effective viscosity by several orders of magnitude. To reflect this assumption, we take µ = 1.2Pa S [16]. The vesicles we consider in this paper are recycling endosomes, which serve to replenish surface proteins and vary from 0.5-4µm in diameter [10,16]. Recycling endosomes are greater in size relative to smaller vesicles on the order of ten to hundreds of nanometers, which may serve other purposes and have been observed in dendritic spines (based on data from VAST Lite [5]). For simplicity, and because we do not consider non-recycling endosomes, we will refer to recycling endosomes as vesicles throughout this paper. Our model assumes there are two species of myosin motors that are identical except for one difference: one species prefers to push the vesicle up towards the spine head and the other prefers to push the vesicle down away from the spine head (Equivalently, one could assume a single species of myosin motors but allow for actin filaments that are aligned parallel or antiparallel to the spine neck). Forces exerted by each species depend on the vesicle velocity. Following the notation of [16], the net motor force is written and F A (U ) are the force-velocity curves of motors that push towards and away from the spine head, respectively. The parameter φ represents the ratio of motor populations: φ = 0 corresponds to only downwards-pushing motors, φ = 1 corresponds to only upwards-pushing motors, and φ = 0.5 corresponds to equal numbers of motors pushing up and down.
For a given species, when the vesicle moves in the preferred direction, the motors attach and detach with intrinsic rates α and β, respectively. In the non-preferred direction, the motors not only detach due to the rate β but are subject to yield effects: the motors extend up to a finite extension, beyond which the motors yield and no longer exert a force. The force p(z) exerted by each motor depends on their individual position z and is generally a monotonically increasing function with p(0) = 0. In the present study, we use where p 1 and γ are the motor force parameters (note that the position z in Equation (5) represents the relative position of an individual motor, which is distinct from the z used in the height function Equation (4). We will no longer reference p(z) to eliminate confusion, and any further reference to position z will refer to the absolute position used in Equation (4)). With this choice of force-extension, in the limit of large motor number, the forces in the preferred and non-preferred directions are functions of velocity as characterized by the force-velocity curves. For upwards-pushing motors, the force velocity curve, F −A (U ), is given by where c(U ) = 1 − exp[β(B − A)/U ]. Because the downwards-pushing motors follow precisely the same rules but for opposite signs in force and velocity, it follows that F A (U ) = −F −A (−U ). We refer the reader to [16,22] for details on the derivation of the force-velocity functions F −A (U ) and F A (U ).

Nondimensionalized Lubrication Model
To enable an analysis of the dynamics of Eqs. (1)-(4), we first reduce the equations: we plug in Equation (1) into Equation (4) and Equation (3) into Equation (2), yielding a system of two equations for the velocity U and the height between the vesicle and the constriction wall h(z): Next, we nondimensionalize Equations (7) and (8) and take z = z/R p , h = h/R p , and is the stall force. Note that here we use tildes to denote dimensionless quantities (tildes will be dropped later on).

Fast-Slow Lubrication Model
From this point on we work exclusively with the nondimensional system unless explicitly stated. Therefore, we write Z = Z, U = U , h = h, ζ =ζ and F = F . Note that Equation (10) includes a term to account for vesicle compliance with a prefactor of π 2 . Representative parameters of vesicle compliance reveal the nondimensional compliance π 2 to be relatively small, on the order of π 2 ≈ 0.09 [16]. To a first approximation, we take π 2 ≈ 0, which greatly simplifies Equations (9) and (10). This low compliance limit yields the fast-slow system, where F is the dimensionless net motor force, U is the dimensionless vesicle velocity, Z is the dimensionless vesicle position, ζ is the dimensionless drag that captures information about the constriction geometry (Equation (11)), and ε is a dimensionless mass term that may be zero, and equals zero in the overdamped limit. Example dynamics of Equation (13) are shown in Figure 1. In Figure 1A, we show a representative idealized dendritic spine. The base of the spine is marked by a vertical gray dashed line positioned at the dimensional position of −5µm, with a base diameter of 6µm. The spine transitions linearly into the constriction, which has a radius of R c = 1.22µm, and the length of the constriction is taken to be 5µm. The vesicle, shown as a black circle, has a radius of R p = 0.96µm. The first initial condition we consider starts the vesicle at the base of the spine with positive initial velocity (U 0 , Z 0 ) = (0.43µ/s, −5µm). Solutions to this initial condition are shown by black curves. As the vesicle moves into the constriction, confinement effects at the neck greatly reduce the velocity of translocation ( Figure 1B, black). However, the vesicle position increases until it reaches the end of the channel ( Figure 1C, black). We show another initial condition which starts at the base with negative initial velocity, (U 0 , Z 0 ) = (−0.3µ/s, −5µm), and denote solutions of this initial condition in red. The velocity of the vesicle remains negative until it hits the nopenetration boundary condition, which we impose at the two ends of the channel. When the vesicle hits the base of the spine at −5µm, the velocity is instantaneously reset to zero and the position to −5µm. This zero velocity happens to be in the basin of attraction of a negative velocity, so the vesicle remains at the base. Different initial conditions reveal different long-time dynamics, suggesting multistability.
We show another representative spine in Figure 1D, where the vesicle and constriction values are 1.5µm and 2.15µm, respectively (all other geometry parameters are the same as in Panel A). The two initial conditions are (U 0 , Z 0 ) = (0.17µ/s, −5µm) (black) and (−0.1µ/s, 0µm) (red). Note that while the first initial condition (black) successfully translocates, the second initial condition (red) does not ( Figure 1F). Again, the differing dynamics as a function of initial conditions suggests the system is multistable. We remark that while we use piecewise linear channels throughout this paper, any constriction geometry is allowed as long as the vesicle radius is close to the channel radius, i.e., R p ≈ R c .
As a starting point for our analysis, we explore system (13) under two equivalent limits that reveal different aspects of the dynamics. Viewed in the "slow" time t, taking the limit ε → 0 yields the slow subsystem, The dynamics of Equation (14) exist on the critical manifold S 0 defined by The goal of this paper is to classify the bifurcations of this manifold in order to understand the set of dynamics of this system. Fenichel theory guarantees that for ε sufficiently small, the dynamics of Equation (13) closely follow the dynamics on the slow manifold [17,8].
While in the present study we operate largely within the overdamped limit, where ε = 0, we sometimes take ε > 0 to integrate the equations numerically using standard methods as in Figure 1. Details of this approach are included in Appendix A.1.  Figure 1A,D. The critical manifolds were computed by applying a range of U and Z values to the critical manifold equation in Equation (15) and extracting points that most closely satisfy the condition. The dynamics on the manifold are determined by a hybrid system: for a given set of initial conditions (Z 0 , U 0 ), the fast dynamics instantaneously carry the solution to the nearest stable manifold. Along the stable manifold, the slow dynamics evolve according to (14), until the solution reaches a fold, at which point the fast dynamics instantaneously carry the solution to the next stable manifold. These hybrid dynamics are shown by dashed gray arrows for the fast dynamics and black arrows for the slow dynamics in Figure 2A,D.
Let s = t/ε, and call s the "fast" time. The fast subsystem (often referred to as the layer problem) yields a substantially more tractable version of the lubrication model. In particular, the fast subsystem is much easier to analyze using bifurcation theory. A straightforward application of the chain rule yields Thus, from the perspective of the fast time s, Z is a slow variable. Letting ε → 0 yields the fast subsystem, These equations are much simpler than the slow subsystem and is equivalent to the quasisteady-state approximation. Because Z is constant in this limit, we only need to analyze the following one-dimensional ODE: where ζ can be treated as a parameter. The bifurcations of Equation (18) are related to the slow subsystem (Equation (14)) through the viscous drag term. Because ζ is a function of position Z, there exists a mapping from the critical manifold to the bifurcation curve. For example, as the vesicle center of mass approaches the center of the constriction, viscous drag increases monotonically ( Figure 2C). At this stage, the bifurcation curve and critical manifold closely resemble scaled versions of each other (Figures 2A,B and Figure 2D,E when Z ∈ [−5µm, 0µm]). Beyond the center of the constriction, the viscous drag term decreases monotonically and the critical manifold resembles a reflected version of the bifurcation curve (Figure 2A,B and Figures 2D,E when Z ∈ [0µm, 5µm]). Thus understanding the bifurcation curves and the viscous drag terms are sufficient to understand the critical manifold of the overdamped system. With this mapping in mind, we turn to a thorough numerical analysis of the bifurcations of this system. In this section, we perform a numerical bifurcation analysis of the fast subsystem by following roots of the right-hand side of Equation (18):

Bifurcations of the Force-Velocity Curve
Details of the numerics are given in Appendix A.2.

Bifurcations in φ-ζ
We begin with single-parameter bifurcation diagrams in φ and ζ by fixing one parameter and varying the other. In Figures 3A-C, we fix ζ at three different values and follow equilibria as a function of φ 1 . The symmetry of these curves about φ = 0.5 comes from our choice of force-velocity curves for competing motors: identical force-velocity curves are used for both species, so the existence of fixed points for φ > 0.5 are the same but with opposite sign when φ < 0.5. As ζ decreases from Panels A to C, the saddle-node denoted by the orange star occurs at progressively smaller values of φ 1 , and the saddle-node denoted by the yellow + occurs at progressively greater values of φ 1 . The change results in the creation of multistable velocities (stable velocities are black and unstable velocities are red). In Panel A, there exist values of φ with only 1 or 3 fixed points. In Panels B and C, the change in the position of the folds gives way to the existence of 5 fixed points. In addition, we find that φ above and below 0.5 tend to yield positive and negative velocities. There are some exceptions, such as Panel C, where φ > 0.5 can result in negative velocity. This observation is in fact the well-established tug-of-war effect. These one-parameter bifurcation diagrams give us a good starting point of how stable velocities change as a function of two parameters as well as some insight into the shape of the bifurcation surface. Note that in Panel A, the parameter range in φ for which there is only 1 solution is much greater relative to the range of φ for where there are 3 solutions. We therefore call the single-velocity solutions more "robust" relative to the multistable solutions. We can also conclude from Panel C that the tug-of-war effect is more likely to be seen for lower values of viscous drag. We will later discuss robustness in a similar way, where relatively greater parameter ranges correspond to increased robustness.
In Figures 3E, F, we fix φ 1 at 0.48 and 0.49, respectively, and vary ζ. These bifurcation diagrams are less intuitive, but nevertheless can be understood as slices through the bifurcation surface described above. More importantly, these bifurcation diagrams reveal some information about the underlying critical manifold. In 3E, we see that if the vesicle has a sufficiently positive initial velocity at low drag, it maintains a positive velocity as the vesicle moves into the spine neck, and the drag grows due to the constriction. At a critical drag denoted by the orange star, the vesicle instantaneously switches velocity in the opposite direction by jumping down to the lower stable branch and eventually exits the constriction through the base of spine at −5µm. Similar discontinuous behavior is seen in 3F: with sufficiently positive initial velocity, the vesicle initially moves towards the constriction before jumping down to the middle stable branch where the velocity is near zero. In this scenario, the vesicle remains stuck for long times. Note that by symmetry of these functions, it follows that the bifurcation diagrams look identical with a change of sign for φ = 0.52, 0.51. Therefore, the model predicts that motor-driven transport through constrictions will generally push the vesicle towards the spine head as long as the initial condition is sufficiently far into the constriction and that upwards-pushing motors are dominant.
While one-parameter diagrams are useful, we wish to understand how multistability changes in the entire φ-ζ parameter space. We address the question of multistability by noting how each one-parameter bifurcation changes as a function of an additional parameter. This process naturally partitions the φ-ζ parameter space into multistable regions. Our one-dimensional system only produces saddle-node (SN) bifurcations, which produce or destroy pairs of fixed points. We can therefore track the number of fixed points throughout the parameter space by tracking these saddle-node bifurcations. This process yields Figure 3D, in which we suppress fixed points and only display bifurcation points and the number of fixed points in each region. Each of the four colored curves correspond to a saddle-node bifurcation, and are labeled with a unique color and symbol. As expected, the number of fixed points changes across the various saddle-node curves. For example, panel A shows that as φ 1 increases, the total number of fixed points changes in the order 1 → 3 → 1 → 3 → 1. The same can be observed in panel A by tracing the slice at ζ = 1.75 × 10 −4 Kg/s. The two-parameter diagram completely characterizes the total number of fixed points for each region of the parameter space and shows that we can expect multistability in much of the displayed parameter space. If viscous drag is sufficiently small (in the bottom region of panel D), multistability exists for a wide range of motor ratios φ 1 . As viscous drag increases, the range of multistability becomes much smaller as fixed points disappear through saddle-node bifurcations. For sufficiently large viscous drag ζ, multistability ceases to exist as the saddle-nodes disappear through cusps, and there exists only one stable velocity for all motor ratios. The critical vesicle-to-spine diameter ratio at the cusps is roughly 2µm/3µm.
In terms of the bifurcation surface, the parameter slices in panels A-C show that for several values of fixed ζ, the bifurcation surface contains four folds. By choosing greater values of ζ, we find that the folds of the surface eventually flatten out through a pair of cusp bifurcations. Beyond this cusp bifurcation, for greater values of ζ, we expect that the velocities U are negative when φ 1 < 0.5, positive when φ 1 > 0.5, and zero when φ 1 = 0.5. From this geometric intuition, it follows that when confinement effects are sufficiently large, the velocity of the vesicle is determined purely by the ratio of upwards-and downwardspushing motors. We make this observation more rigorous in Section 4.

Robustness in π 4 , π 5
The two-parameter figure in φ-ζ provides a complete description of how multistability changes as a function of motor ratio and constriction geometry when motor parameters π 3 − π 6 are fixed. However, there may be variations in motor parameters due to the existence of multiple motor types such as myosin V and VI [10] and variations in ATP and ADP concentration, which is known to differentially modulate myosin motor dynamics [42].
As explained in detail in Appendix A.3, the cusp bifurcation separates the parameter space between multistable velocities and globally stable or unstable solutions. Indeed, cusps generally serve as a sufficient condition for the existence of hysteresis. Therefore, understanding the cusp bifurcation may provide important insights into controllability. For a given ζ, it is possible to track the location of these cusps as a function of π 4 and π 5 . The result of this process for various π 4 and ζ is shown in Figure 4A. Each curve represents a cusp bifurcation as a function of π 4 , π 5 for a given ζ (we briefly describe how we determined the location of these cusps in Appendix A.3).
Given parameters π 4 , π 5 chosen somewhere in the π 4 -π 5 parameter region below the ζ = 0 level curve in Figure 4, the ζ level curves suggest that there exists a cusp bifurcation at some ζ * (π 4 , π 5 ) > 0. Because cusps are a sufficient condition for the existence of multistability, it follows that multistable states exist for appropriate choices of φ 1 and ζ ≤ ζ * . The ζ = 0 level curve represents a sufficient condition for the loss of multistability: noting that for each fixed point, ∂ζ/∂π 5 < 0 in a neighborhood about the fixed point, i.e., ζ decreases as π 5 increases, it follows that for any π 4 , π 5 in at least a neighborhood above this curve, there can exist no cusp bifurcation with a positive ζ.
Various points in Figure 4A are marked B-E, with corresponding two-parameter bifurcation diagrams shown in the remaining subplots, Figures 4B-E. The point marked by in Figure 4A represents the π 4 and π 5 values for Figure 3. These diagrams provide evidence that for smaller π 5 and greater π 4 , the region of multistability tends to increase. Recalling that π 5 is the maximum displacement of a motor in its non-preferred direction, smaller π 5 implies that bidirectional motion due to noise can be made more likely by allowing the motor to detach earlier. Next, π 4 is the initial motor attachment position in either the preferred or non-preferred direction, and we find that the area of multistability increases as motors have greater initial extension. Together, we predict that strong initial attachment forces combined with greater yield effects can result in more frequent directional switching.

Existence
The existence of a stationary solution U = 0 is straightforward to prove by inspection for φ = 0.5. In this section, we expand about this solution to determine the existence of solutions when φ is near the equal motor ratio φ = 0.5 and when ζ is large. We let φ = 0.5 +φ andζ = 1/ζ and explore the cases whereφ andζ are small. These limits inform us of the local linear behavior of the velocity function U = C(φ,ζ) by writing First we derive an equation for the small deviation φ = 0.5 +φ, where 0 < |φ| 1. Constant velocity solutions U = C must satisfy Solving forφ yields,φ Our nondimensional parameters are always greater than zero, so 0 < 1 − e π 4 and 0 < π 6 , and ζ > 0. It follows that the derivative ofφ with respect to C is nonzero: e π 4 π 6 + ζ e π 4 − 1 = 0.
We then obtain a local equation for C as a function ofφ by invoking the inverse function theorem: We have seen that in some parameter regimes, viscous drag grows to large values during translocation ( Figure 2C,F), but nonzero velocities persist due to unequal motor ratios ( Figure 1C,F). To establish this inverse relationship between velocity and drag, we Taylor expand about infinity withζ = 1/ζ as defined above. Then solving forζ in terms of C yieldsζ To derive a local equation for C as a function ofζ, we examine the derivative ofζ with respect to C: dζ dC .
So long as φ = 1/2 and π 4 > 0, this derivative is well-defined, and we can invoke the inverse function theorem to write the local equation for C as a function ofζ.
After a trivial substitution we arrive at the desired equation, Combining the local velocity estimates, we arrive at the velocity equation as a function of φ and ζ: This existence equation is valid for large ζ and anyφ ∈ [−1/2, 1/2], or for smallφ and any ζ. Equation (21) provides justification for the effects observed from simulations earlier in the paper. For ζ large, the velocity must be small and proportional to ζ −1 , and the sign of the velocity is determined by the sign ofφ, that is, it is the ratio of motors that determines the direction of motion for large drag. In the case of small drag, this equation is only valid for small φ, but the equation yields the same intuition that is consistent with our one-parameter bifurcation diagrams, i.e., that the sign of the velocity is determined by the dominant motor species. Finally, we can make additional predictions. In the case of large drag, or small drag with near-equal (but non-equal) motor ratios, increasing motor attachment position and increasing the velocity scale ratio between translocation and motor adhesion dynamics (π 4 and π 6 ) will increase the velocity of translocation. In addition, velocity does not locally depend on the ratio of motor attachment and detachment parameters and the maximum displacement of each motor (π 3 or π 5 ).

Linear Stability
For our one-dimensional problem, a linear stability analysis is sufficient to understand the stability of fixed points, which follows from the slope of the total force function: Note that if the derivative of F is bounded, it follows that the eigenvalue is negative (λ < 0) for sufficiently large drag forces. Therefore, all velocities satisfying force-balance are stable. We can rule out multistability for large drag because, if there is more than one fixed point, there must also be more than one unstable fixed point. This is impossible by the negativity of λ shown above. By using continuity of the force-velocity curves (Appendix B.1), the proof follows by contradiction: if there are only two fixed points that are both stable, the slope of the net force function must be negative at each point. By the intermediate value theorem, there must exist a third point between them, which contradicts the original claim of two stable points. So either one point must be unstable or there is only one fixed point. This argument can be extended to eliminate any number of stable fixed points in the case of large drag.
Recalling that π 4 , π 6 > 0, we generally expect constant velocity solutions to be stable for φ small.

Discussion
In this paper, we fully characterize the dynamics predicted by a model of vesicles driven into closed constrictions. We cast the system into a fast-slow system and perform a twoparameter bifurcation analysis on the fast subsystem, and then determine how the cusps in the resulting bifurcation surface vary with the motor parameters. The model predicts multistability (i.e. bidirectional motion) for smaller values of viscous drag and unidirectional motion for greater values of the viscous drag corresponding to tight constrictions. We remark that while dendritic spines exhibit a great diversity of morphologies such as thin spines that are often less than 2µm in length with neck diameters ranging from 0.06-0.2µm [3] and mushroom-shaped spines that are often less than 1µm in width [37], most can be captured by our reduced axisymmetric lubrication model. Moreover, while vesicles of recycling endosomes range in diameter from small to large (roughly in the range from 0.5µm to 2µm [10]) and deform strongly in the spine neck, the most important requirement for us to apply the theory of this paper is that the diameter of the vesicle and spine wall are similar. While detailed electron microscopy images exist of dendritic spines [24], it is a significant challenge to efficiently search and classify recycling endosomes. Therefore, for now we rely on published endosome images to approximate the physiologically relevant parameter ranges.
In Figure 5A,B we show two representative experimental images from [10] of spines containing recycling endosomes. The spine is outlined by a thin yellow line and the endosome is shown traveling through the spine neck in a series of four time-lapse images with an associated kymograph on the right. The scale bar is 2µm in all panels. Using these images, we determine approximate regions where physiological parameters may lie in theφ − ζ parameter space (Figures 3D and 5C). Figure 5C is the same as Figure 3D but with ζ plotted on a log scale. The time-lapse images provide virtually no information about the ratio of motors, so we make no restrictions on φ for now.
We estimate the viscous values from these images by estimating the height between the vesicle and spine wall in Figure 5A and assume a constant constriction for simplicity. Through a crude manual approximation, we estimate the height in panel A to be at most 0.1µm. This height is substantially smaller than the heights considered in this paper, and therefore yields a relatively greater viscous drag value of ζ ≈ 2 × 10 −2 Kg/s (we assume  Figure 5: Time-lapse images of recycling endosomes adapted from [10] and available under the CC BY NC ND license. A: a recycling endosome translocates through a thin spine in a series of four time-lapse images (left). A kymograph is shown to the right. The vesicle is roughly 1µm in diameter, and the distance between the vesicle and neck wall is at most 0.1µm. B: a recycling endosome translocates into a stubby spine in a series of time-lapse images, with the associated kymograph on the right. The vesicle is roughly 0.5µm in diameter, the distance between the vesicle and neck wall is roughly 0.15µm, and the distance between the vesicle and head wall is roughly 0.5µm. All scale bars 2µm. C: approximate ranges of drag (red transparent) superimposed on the two-parameter diagram from Figure 3 with the drag ζ plotted on a log scale. Labels A and B correspond to Panels A and B.
the same microscopic motor parameters as in Figure 3). This value corresponds to a point far above the cusps, where only a single-velocity solution exists ( Figure 5C, top red band labeled "A"). This range lies well within the unidirectional regime. Because the time-lapse images in Figure 5A show the endosome traveling towards the spine head, we infer that the upwards-pushing motors are dominant. The only region where the vesicle could switch direction is at the base of the spine, where there some distance between the diameter of the cell wall and the recycling endosome is possible. Thin spines generally are much smaller in diameter relative to recycling endosomes, so that the distance between the vesicle and spine wall is small. We conclude that once the vesicle has entered a thin spine, we generally expect unidirectional movement. Figure 5B shows an example of a stubby spine. Here, the endosome is smaller, with a diameter range from 0.5µm to 0.6µm. The lower bound of the height range is difficult to establish, but we estimate that the height between the endosome and neck is roughly 0.05-0.1µm, which yields viscous drag values in the range ζ ≈ 1 × 10 −4 to 6 × 10 −4 Kg/s. Recalling that the cusps in Figure 3 are on the order of 2 × 10 −4 Kg/s (where the critical vesicle-to-spine diameter ratio is roughly 2µm/3µm), we find that both unidirectional and multidirectional motion are possible ( Figure 5C, middle red band labeled "B, neck"), where unidirectional solutions take a much greater portion of the parameter range in a linear scale. Another property to consider is that stubby spines have shallow constrictions that lead into a large head, where the vesicle may spend a substantial amount of time relative to the neck. It is therefore worth considering whether multistable solutions exist past the neck. In the head, the height between the endosome and head ranges from roughly 0.1µm to 0.5µm, which yields a drag range of 3 × 10 −6 Kg/s to 1 × 10 −4 Kg/s, which places the vesicle squarely in a multistable regime ( Figure 5C, lower red band labeled "B, head"). Strikingly, the kymograph shows unidirectional movement through the neck, and bidirectional motion in the spine head.
The viscous drag in these representative examples varies over four orders of magnitude and our model predicts dramatically different qualitative behaviors over this range. According to our model, thin spines with large vesicles will exhibit unidirectional movement due to the large drag experienced by the vesicle. In contrast, stubby spines will exhibit multidirectional motion, especially if the recycling endosome is of smaller size. Interestingly, published images of vesicle movement in spines appear to confirm these claims.
A drawback of the present study is the use of mean-field models. In particular, the forcevelocity functions used in this paper rely on the limit of large numbers of myosin motors [22], and therefore questions about noise cannot be addressed in this framework. One possible approach to overcome this limitation is to replace the current mean-field model of molecular motors by a discrete model; such a discrete model is used to compute mean passage times in Allard et al [2]. Finally, numerical experiments are another feasible approach to address the question of mean first passage times of a translocating vesicle. Similar approaches have been performed in [11] to pursue questions on motor and intermediate filament parameter effects on intermediate filament transport.

Acknowledgments
The authors acknowledge support under National Institute of Health grant T32 NS007292 (YP) and National Science Foundation grant DMS-1913093 (TGF).

A Numerics
In this section we detail the numerical methods and parameters of the fast-slow lubrication model.

A.1 Integration
For convenience, we restate the original fast-slow lubrication model here: We solve Equation (23) numerically by taking ε nonzero and integrating forwards in time (as we have done using the forward Euler method). This approach is beneficial because velocities satisfying instantaneous force-balance, i.e., U such that F (U ) − ζ(Z)U = 0, are fixed points of a one-dimensional ODE (assuming Z is constant). The question of convergence is then very simple: if an initial condition is within the basin of attraction of a stable velocity, the initial condition will converge to this velocity.
Another benefit of taking this approach is that for ε sufficiently small, system (23) is sufficiently similar to the overdamped system when ε = 0, as guaranteed by Fenichel theory [17,8]). When ε 1, there is a separation of timescales so that very small times steps are required for numerical stability. We find that choosing ε = 1 works well in practice, as the dynamics follow the slow manifold but do not require impractically small time steps. This approximation works even at moderate values of ε (e.g., at ε = 1) because, whereas we have nondimensionalized velocity in terms of the free space prediction of Stokes' law, in confined geometries there is a small parameter which may be obtained by rescaling ε by the effective drag. For a given choice of ε, we typically take the time step dt to be smaller by one or two orders of magnitude.

A.2 Continuation
General continuation strategies in one and two parameters can be found in Chapter 10.3 of Kuznetsov [27]. We use XPPAUTO [15] (version 8) to generate our bifurcation diagrams unless stated otherwise. For the generation of the one-and two-parameter diagrams, we use the following numerical values: . The displayed Dsmax value may make AUTO run multiple passes over some branches, in which case we reduce Dsmax to 0.001 preventing this issue while incrementing at a reasonable pace. We include many more details in a mini tutorial on our github repository at https://github.com/youngmp/park_fai_2020.

A.3 Cusps
Saddle-nodes may be computed by finding tangencies in the right-hand side of Eq. (18) given a solution U = C with the constraint that U satisfies (??). The conditions for a saddle-node bifurcation are [27]: We denote any U that simultaneously satisfy these equations by U * . We can simplify the search for saddle-nodes by writing the system We can just as easily useζ and look for fixed points in (Ū ,ζ) as a function of φ. It is simply a matter of preference. We useŪ andφ to emphasize the difference between (25) and the original fast subsystem (18). This new system eliminates one parameter when it comes to computing bifurcation diagrams. Fixed points of (25) correspond to saddlenode bifurcations in the fast subsystem (18). TheŪ nullcline of (25) shows the stable fixed points in (18), and theŪ nullcline intersections with theφ nullcline correspond to saddle-nodes. Thus, the phase space of (25) corresponds to one-parameter bifurcations in (18), and the one-parameter bifurcations in (25) correspond to two-parameter bifurcations in (18). In particular, saddle-node bifurcations of (25) correspond to cusp bifurcations in (18). We exploit the latter fact to generate the diagram of cusp bifurcations ( Figure 4). In practice, we compute intersections in the nullclines of Equation (25) and track where the intersections of these nullclines change in number. Suppose that we are interested in finding cusp bifurcations as a function of π 4 and π 5 for a given ζ. On one side of the cusp bifurcation, Equation (25) shows two fixed points and on the other side shows zero. By defining a function that returns +1 on one side of the bifurcation and −1 on the other, we can use a root-finding method such as Brentq to determine π 5 where there is a transition in the fixed point number to high numerical precision.

B Properties of the force-velocity curve
The force-velocity curves are very well-behaved but it is not always obvious how. In this section we briefly show some properties used in the text.
In particular, there are two possible problem areas: when U = 0 and when U = 1/π 6 . The latter case is especially important to consider because numerical problems occasionally arise when evaluating U = 1/π 6 directly. In the first case, the left limit is straightforward: The right limit depends on the behavior of the exponentials. Noting that exp(−π 5 /π 6 U ) → 0 as U → 0 + , it is straightforward to check that the left and right limits agree independent of the parameters π i , and therefore F A is continuous at U = 0. Continuity of F −A follows by definition.
In the second case, when U = 1/π 6 ≥ 0, we take a close look at the second line of Equation (26). Potential problems arise in the term (1 − e π 5 e −π 5 /π 6 U )/(1 − π 6 U ). For convenience, let v := 1 − 1/(π 6 U ) so that v → 0 as U → 1/π 6 . Then the term becomes 1 − e π 5 e −π 5 /π 6 U ) We have used Taylor expansion of the exponential about zero. So the term converges to π 5 − 1 as v → 0 in either side of the limit. It follows that the function F A is continuous at U = 1/π 6 , as desired.

B.2 Limits
Consider π 4 → 0. This limit occurs when attachments occur at zero position (A = 0), or when the force-scaling parameter γ = 0. Both cases appear to be somewhat unrealistic: in the former case, newly attached crossbridges will apply no force (Equation (5)), and in the latter case, the force exerted by a single motor is always zero. Indeed, we find that the term e π 4 − 1 appears in the denominator of several terms and competes with U when U is small. For |U | 0, the term (e π 4 − 1) −1 dominates, and we expect F A (U ) → ∞ as π 4 → 0. When |U | 1, we expect F A (U ) → −1 as U → 0 for each π 4 small.
In contrast, the dimensional force-velocity curve converges for π 4 small and diverges for π 4 . A nondimensionalization using F 1 = F 0 /(e π 4 − 1) would remove this problem of divergence, but this rescaling is not necessary when considering bifurcations: for any fixed point F (U, π 4 ) = 0, where F is the nondimensional net force, multiplying both sides by the scaling factor (e π 4 − 1) yields (e π 4 − 1)F (U, π 4 ) = 0, so in fact scaling preserves fixed points as a function of π 4 .