Dynamics of Vesicles Driven Into Closed Constrictions by Molecular Motors
TThe Dynamics of Vesicles Driven Into Closed Constrictions byMolecular Motors
Youngmin Park ∗ and Thomas G. Fai Department of Mathematics, Brandeis University, Waltham, MA 02453
Abstract
We study the dynamics of a model of membrane vesicle transport into dendriticspines, which are bulbous intracellular compartments in neurons driven by molecularmotors. We reduce the lubrication model proposed in [Fai et al, Active elastohydrody-namics of vesicles in narrow, blind constrictions. Phys. Rev. Fluids, 2 (2017), 113601]to a fast-slow system, yielding an analytically and numerically tractable equation thatis equivalent to the original model in the overdamped limit. The key parameters ofthe model are the ratio of motors that prefer to push toward the head of the dendriticspine to the ratio of motors that prefer to push in the opposite direction. We performa numerical bifurcation analysis in these parameters and find that steady-state vesiclevelocities appear and disappear through a number of saddle-node bifurcations. Thisallows us to identify the region of parameter space in which multiple stable velocitiesexist. We show by direct calculations that for sufficiently close vesicle-to-spine diameterratios, there can only be unidirectional motion. The critical vesicle-to-spine diameterratio predicted by our analysis, at which there is a transition from unidirectional tobidirectional motion, is consistent with experimental observations of vesicle trajectoriesin the literature.
Intracellular transport is required by cells for many essential purposes, e.g. to deliver pro-teins to their functional locations and to maintain concentration gradients. The relativelylong lengths of axons and dendrites means that diffusion is too slow to serve as a prac-tical transport mechanism in neurons, making directed transport particularly important.Not only size, but also geometry, may play an important role. Bulbs known as dendriticspines cover dendritic trees of pyramidal neurons and are ubiquitous in the cortex: up to85% of the cortex consists of pyramidal neurons [14], and up to 85% of excitatory inputs ∗ Corresponding author [email protected] a r X i v : . [ q - b i o . CB ] A p r erminate on dendritic spines [34]. Spines and their distinctive shapes–with narrow necksof significant curvature separating the spines from the dendritic branch–are likely to havea significant influence on transport.Dendritic spines protrude 0.5–1 µ m [38] from the dendritic shaft and often appearmushroom-like in shape. Spines exhibit a thin 50–100nm [21] neck that opens into ahead that is up to several hundred nanometers greater in diameter [29]. A key aspect ofdendritic spine remodeling and maintenance includes the transport of surface proteins onvesicles that are squeezed through the narrow neck into the spine head [10, 35, 40, 16]. Re-cent experiments have shown that movement is not always unidirectional (translocation),but includes no movement (corking), and rejection [35, 40]. The mechanisms underlyingthese changes in direction are not well understood.There are generally two approaches taken by theoretical studies of molecular motor-driven intracellular transport. In the first, the motor complex is assumed to have switchingrates determined by some underlying Markov process, and consequences of varying switch-ing rates are explored, such as the distribution of motor complex velocities [30, 26], andmean 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 popula-tion, which generatively produces bidirectional motion despite the assumption of symmetry[23, 20, 2, 36]. However, the effect of constrictions on motor dynamics has typically beenheld constant or neglected.Constricted motion is of general interest in engineering and biology with related prob-lems appearing in manufactured elastic capsules [12, 13], hydrogels [28], the movement ofliving cells [4, 9, 18], and axonal transport [43, 25, 39]. In terms of theory, the effect ofconstrictions on molecular dynamics is a very recent topic of interest. In Fai et al. 2017, itwas found that constrictions enrich the dynamics of motor populations, allowing motors toswitch between multidirectional and unidirectional motion [16]. The goal of this paper is tothoroughly explore the dynamics of the model derived in [16] by classifying the bifurcationsof the underlying ODE model.The paper is organized as follows. In Section 2, we briefly review the derivation ofthe lubrication model and its nondimensionalization, and we briefly discuss two equivalentversions, one of which (the “slow” subsystem) is the model considered in [16], and otherof which (the “fast” subsystem) is the model we explore in depth in this paper. In Sec-tion 3, we numerically establish the existence and robustness of multistability through abifurcation analysis of the “fast” subsystem. Some of these numerical results are corrob-orated in Section 4, where we analytically establish the existence and stability of partic-ular velocities as a function of key parameters. We conclude the paper with a discussionthat includes estimates of realistic parameter regimes of this model, and the resultingbehaviors predicted for different types of dendritic spines. All code and data and docu-mentation 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 theauthors consider an idealized model of vesicle translation through the neck of a dendriticspine. In contrast to previous studies of transport through constrictions in periodic orunbounded tubes, the authors close off one end to model transport into spine-like ge-ometries. Objects entering so-called closed constrictions are surrounded by fluid backflow,which introduces a large velocity gradient, allowing for a simplification of the Navier-Stokesequations 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 pressureexerted on the vesicle at position z ∈ [ Z − R p , Z + R p ]. The flux ( Q ) equations follow fromthe observation that the amount of fluid carried by the vesicle must match the amountof fluid escaping the closed constriction, i.e., at position Z , the cross-sectional area of thevesicle is balanced by the fluid backflow determined by the distance h between the vesi-cle and spine wall. This results in the following reduced axisymmetric model of vesicletrafficking:[ p ( z ) − p ]6 µ = U (cid:90) zZ − R p h ( s ) ds − Q πR c (cid:90) zZ − R p h ( s ) ds, (1) Q = 2 πR (cid:32) U (cid:90) Z + R p Z − R p h ( s ) ds − F ( U )6 πR µ (cid:33) (cid:44)(cid:32) (cid:90) Z + R p Z − R p h ( s ) ds (cid:33) , (2) Q = − πR c U, (3) h ( z ) = ˜ R c ( z ) − (cid:113) R p − ( z − Z ) + C [ p ( z ) − p ] . (4)Here, F ( U ) is the net force from the molecular motors and is a key feature that contributesto the tug-of-war dynamics of the model. While the viscocity of water is on the order of µ = 0 . µ = 1 . . µ m in diameter [10, 16]. Recyclingendosomes are greater in size relative to smaller vesicles on the order of ten to hundreds ofnanometers, 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 considernon-recycling endosomes, we will refer to recycling endosomes as vesicles throughout thispaper.Our model assumes there are two species of myosin motors that are identical exceptfor one difference: one species prefers to push the vesicle up towards the spine head andthe other prefers to push the vesicle down away from the spine head (Equivalently, onecould assume a single species of myosin motors but allow for actin filaments that are3ligned parallel or antiparallel to the spine neck). Forces exerted by each species dependon the vesicle velocity. Following the notation of [16], the net motor force is written F ( U ) = φF − A ( U ) + (1 − φ ) F A ( U ), where F − A ( U ) and F A ( U ) are the force-velocity curvesof motors that push towards and away from the spine head, respectively. The parameter φ represents the ratio of motor populations: φ = 0 corresponds to only downwards-pushingmotors, φ = 1 corresponds to only upwards-pushing motors, and φ = 0 . α and β , respectively. In the non-preferred direction, themotors not only detach due to the rate β but are subject to yield effects: the motors extendup to a finite extension, beyond which the motors yield and no longer exert a force. Theforce p ( z ) exerted by each motor depends on their individual position z and is generally amonotonically increasing function with p (0) = 0. In the present study, we use p ( z ) = p ( e γz − , (5)where p 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 inthe 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 forcesin the preferred and non-preferred directions are functions of velocity as characterized bythe force-velocity curves. For upwards-pushing motors, the force velocity curve, F − A ( U ),is given by F − A ( U ) = αn p αc ( U )+ β e γA ( − e − β ( B − A ) /U e γ ( B − A ) ) − ( − e − β ( B − A ) /U ) (1+ γU/β )1+ γU/β , U < αn p α + β e γA − − γU/β γU/β , U ≥ , (6)where c ( U ) = 1 − exp[ β ( B − A ) /U ]. Because the downwards-pushing motors follow pre-cisely 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-velocityfunctions F − A ( U ) and F A ( U ). To enable an analysis of the dynamics of Eqs. (1)–(4), we first reduce the equations: weplug in Equation (1) into Equation (4) and Equation (3) into Equation (2), yielding asystem of two equations for the velocity U and the height between the vesicle and the4onstriction wall h ( z ): U = F ( U )6 πRµ (cid:82) Z + R p Z − R p R c h ( s ) + Rh ( s ) ds , (7) h ( z ) = ˜ R c ( z ) − (cid:113) R p − ( z − Z ) + C µ (cid:34) U (cid:90) zZ − R p h ( s ) ds − Q πR c (cid:90) zZ − R p h ( s ) ds (cid:35) . (8)Next, we nondimensionalize Equations (7) and (8) and take (cid:101) z = z/R p , (cid:101) h = h/R p , and (cid:101) U = 6 πR p µU/F , where F = (exp( γA ) − αp n / ( α + β ) is the stall force. Note thathere we use tildes to denote dimensionless quantities (tildes will be dropped later on).Plugging the nondimensionalized terms into Equations (7) and (8) yields, (cid:101) U = (cid:101) F ( (cid:101) U ) / ˜ ζ ( (cid:101) Z ) , (9) (cid:101) h ( (cid:101) Z + (cid:101) z ) = π ( (cid:101) Z + (cid:101) z ) − (cid:112) − (cid:101) z + π (cid:101) U (cid:90) (cid:101) z − (cid:101) h − ( (cid:101) Z + s ) + (cid:101) h − ( (cid:101) Z + s ) ds, (10)where π = R c ( (cid:101) z ) /R p , π = CF / ( πR p ), and˜ ζ ( (cid:101) Z ) = (cid:90) − (cid:101) h − ( (cid:101) Z + s ) + (cid:101) h − ( (cid:101) Z + s ) ds. (11)The function ˜ ζ is the viscous drag coefficient produced by the constriction geometry. Weshow examples of this function and discuss its importance for the dynamics of this systemin the next section. Concrete examples of ˜ ζ are shown later on.The nondimensionalized net motor force (cid:101) F ( (cid:101) U ) = φ (cid:101) F − A ( (cid:101) U ) + (1 − φ ) (cid:101) F A ( (cid:101) U ) consists oftwo functions (cid:101) F − A ( (cid:101) U ) and (cid:101) F A ( (cid:101) U ), which are related by (cid:101) F A ( (cid:101) U ) = − (cid:101) F − A ( − U ) using thesame symmetry arguments in the dimensional equations. It is straightforward to show thatthe nondimensionalized force-velocity curve is (cid:101) F A ( (cid:101) U ) = − π (cid:101) U ( e π − − − π (cid:101) U , if (cid:101) U < − ( π +1) π (1 − e − π /π (cid:101) U )+1 [ e π (1 − e π e − π /π (cid:101) U )] − (1 − π (cid:101) U )(1 − e − π /π (cid:101) U )( e π − − π (cid:101) U ) , if (cid:101) U ≥ . (12)The function includes numerous parameters representing various microscopic motor prop-erties: π = α/β (ratio of attachment and detachment rates), π = γA (the nondimensionalattachment position), π = γ ( B − A ) (the maximum displacement of a motor in its non-preferred direction), and π = ( F / [6 πR p µ ]) / ( β/γ ) (the ratio of velocity scales betweentranslocation and motor adhesion dynamics). When conversions back to dimensional forcesare needed, we write F X = (cid:101) F X F for X = A, − A .5 R p R c µ m µ m A . . U ( µ m / s ) B Initial Condition 1Initial Condition 2 − Z ( µ m ) CD t (s) . . U ( µ m / s ) E t (s) − Z ( µ m ) F Figure 1: Constriction geometry and resulting dynamics for different parameter sets. A:The initial spine diameter (6 µ m) decreases to the neck radius R c = 1 . µ m. The vesicle(black circle, radius R p = 0 . µ m) begins at the base of the channel (dashed vertical grayline) and moves in the direction of the arrow for initial condition 1. B,C: Resulting velocity U ( µ /s) and position Z ( µ m) plotted over time (s) for two different initial conditions( U , Z ) = (0 . µ /s , − µ m) (black) and ( − . µ /s , − µ m) (red). φ = 0 . π = 1, π = 4 . π = 0 . π = 10, F = 50. D,E,F: same as A,B,C, but with R c = 2 . µ m, R p = 1 . µ m, φ = 0 . π = 1, π = 4 . π = 0 . π = 10, F = 200. Initial spinediameter is 6 µ m. The two initial conditions are ( U , Z ) = (0 . µ /s , − µ m) (black) and( − . µ /s , µ m) (red). Simulation parameters ε = 1, dt = 0 .
02, integrated using forwardEuler (see Appendix A.1).
From this point on we work exclusively with the nondimensional system unless explicitlystated. Therefore, we write Z = (cid:101) Z , U = (cid:101) U , h = (cid:101) h , ζ = ˜ ζ and F = (cid:101) F . Note thatEquation (10) includes a term to account for vesicle compliance with a prefactor of π .Representative parameters of vesicle compliance reveal the nondimensional compliance π to be relatively small, on the order of π ≈ .
09 [16]. To a first approximation, we take π ≈
0, which greatly simplifies Equations (9) and (10). This low compliance limit yieldsthe fast-slow system, dZdt = U,ε dUdt = F ( U ) − ζ ( Z ) U. (13)6here F is the dimensionless net motor force, U is the dimensionless vesicle velocity, Z isthe dimensionless vesicle position, ζ is the dimensionless drag that captures informationabout the constriction geometry (Equation (11)), and ε is a dimensionless mass term thatmay be zero, and equals zero in the overdamped limit.Example dynamics of Equation (13) are shown in Figure 1. In Figure 1A, we show arepresentative idealized dendritic spine. The base of the spine is marked by a vertical graydashed line positioned at the dimensional position of − µ m, with a base diameter of 6 µ m.The spine transitions linearly into the constriction, which has a radius of R c = 1 . µ m,and the length of the constriction is taken to be 5 µ m. The vesicle, shown as a blackcircle, has a radius of R p = 0 . µ m. The first initial condition we consider starts thevesicle at the base of the spine with positive initial velocity ( U , Z ) = (0 . µ /s , − µ m).Solutions to this initial condition are shown by black curves. As the vesicle moves into theconstriction, 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 thechannel (Figure 1C, black). We show another initial condition which starts at the basewith negative initial velocity, ( U , Z ) = ( − . µ /s , − µ m), and denote solutions of thisinitial condition in red. The velocity of the vesicle remains negative until it hits the no-penetration boundary condition, which we impose at the two ends of the channel. Whenthe vesicle hits the base of the spine at − µ m, the velocity is instantaneously reset to zeroand the position to − µ m. This zero velocity happens to be in the basin of attraction ofa negative velocity, so the vesicle remains at the base. Different initial conditions revealdifferent long-time dynamics, suggesting multistability.We show another representative spine in Figure 1D, where the vesicle and constric-tion values are 1.5 µ m and 2.15 µ m, respectively (all other geometry parameters are thesame as in Panel A). The two initial conditions are ( U , Z ) = (0 . µ /s , − µ m) (black)and ( − . µ /s , µ m) (red). Note that while the first initial condition (black) successfullytranslocates, the second initial condition (red) does not (Figure 1F). Again, the differingdynamics as a function of initial conditions suggests the system is multistable. We re-mark that while we use piecewise linear channels throughout this paper, any constrictiongeometry 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 limitsthat reveal different aspects of the dynamics. Viewed in the “slow” time t , taking the limit ε → dZdt = U, F ( U ) − ζ ( Z ) U. (14)The dynamics of Equation (14) exist on the critical manifold S defined by S := { ( Z, U ) ∈ R | F ( U ) − ζ ( Z ) U } . (15)7he goal of this paper is to classify the bifurcations of this manifold in order to understandthe 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 ε > − . − . . . . Z ( µ m) − . − . . . . U ( µ m / s ) A Critical Manifold
Stable ManifoldUnstable Manifold ζ (Kg/s) × − − . − . . . . U ( µ m / s ) B Bifurcation Curve − . − . . . . Z ( µ m) . . . ζ ( K g/ s ) × − C Viscous Drag − . − . . . . Z ( µ m) − . − . . . . U ( µ m / s ) D Critical Manifold ζ (Kg/s) × − − . − . . . . E Bifurcation Curve − . − . . . . Z ( µ m) ζ ( K g/ s ) × − F Viscous Drag
Figure 2: The mapping between the bifurcation diagram and critical manifold throughthe viscous drag function. A: An example of the critical manifold in the phase space of U and Z . Black arrow heads denote the direction of motion on the slow manifold. Graydashed arrows indicate the direction of motion in the fast system. B: An example of aone-parameter bifurcation diagram, where steady-states U are plotted as a function of ζ .C: The relationship between viscous drag ζ and position Z . The dimensional positions Z = − µ m through Z = 0 µ m, we expect the critical manifold to resemble a version ofthe bifurcation diagram given by the mapping between drag ζ and position Z . Beyond theconstriction, from dimensional positions Z = 0 µ m through Z = 5 µ m, the critical manifoldresembles a reflected version of the bifurcation curve. A–C: Parameters as in Figure 1A–C.D–F: Parameters as in Figure 1D–F.In Figure 2A,D we show examples of the critical manifold for the corresponding ge-ometries in 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 points8hat most closely satisfy the condition. The dynamics on the manifold are determined bya hybrid system: for a given set of initial conditions ( Z , U ), the fast dynamics instanta-neously carry the solution to the nearest stable manifold. Along the stable manifold, theslow dynamics evolve according to (14), until the solution reaches a fold, at which pointthe fast dynamics instantaneously carry the solution to the next stable manifold. Thesehybrid dynamics are shown by dashed gray arrows for the fast dynamics and black arrowsfor the slow dynamics in Figure 2A,D.Let s = t/ε , and call s the “fast” time. The fast subsystem (often referred to asthe 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. Astraightforward application of the chain rule yields dZds = εU (16) dUds = F ( U ) − ζ ( Z ) U. (17)Thus, from the perspective of the fast time s , Z is a slow variable. Letting ε → dZds = 0 dUds = F ( U ) − ζ ( Z ) U. These equations are much simpler than the slow subsystem and is equivalent to the quasi-steady-state approximation. Because Z is constant in this limit, we only need to analyzethe following one-dimensional ODE: dUds = F ( U ) − ζU, (18)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 mappingfrom the critical manifold to the bifurcation curve. For example, as the vesicle centerof mass approaches the center of the constriction, viscous drag increases monotonically(Figure 2C). At this stage, the bifurcation curve and critical manifold closely resemblescaled versions of each other (Figures 2A,B and Figure 2D,E when Z ∈ [ − µ m , µ m]).Beyond the center of the constriction, the viscous drag term decreases monotonically andthe critical manifold resembles a reflected version of the bifurcation curve (Figure 2A,Band Figures 2D,E when Z ∈ [0 µ m , µ m]). Thus understanding the bifurcation curves andthe viscous drag terms are sufficient to understand the critical manifold of the overdampedsystem. With this mapping in mind, we turn to a thorough numerical analysis of thebifurcations of this system. 9 Bifurcations of the Force-Velocity Curve
Figure 3: Two parameter bifurcation diagram in φ and ζ . Saddle-node (SN) bifurcationsare shown in (D) as colored branches with a unique color and symbol for each branch.Numbers in (D) indicate the total number of fixed points in the corresponding region ofparameter space. Subplots A, B, C, E, F, show one-parameter slices of the two-parameterdiagram. Saddle-nodes are labeled with the corresponding branch color and symbol. Thecritical vesicle-to-spine diameter ratio at the cusps is roughly 2 µ m/3 µ m.In this section, we perform a numerical bifurcation analysis of the fast subsystem byfollowing roots of the right-hand side of Equation (18): f ( U ) = φF − A ( U ) + (1 − φ ) F A ( U ) − ζU. (19)10etails of the numerics are given in Appendix A.2. φ - ζ We begin with single-parameter bifurcation diagrams in φ and ζ by fixing one parameterand varying the other. In Figures 3A–C, we fix ζ at three different values and followequilibria as a function of φ . The symmetry of these curves about φ = 0 . φ > . φ < . ζ decreases from Panels A to C, the saddle-node denoted by the orange star occursat progressively smaller values of φ , and the saddle-node denoted by the yellow + occursat progressively greater values of φ . The change results in the creation of multistablevelocities (stable velocities are black and unstable velocities are red). In Panel A, thereexist values of φ with only 1 or 3 fixed points. In Panels B and C, the change in the positionof the folds gives way to the existence of 5 fixed points. In addition, we find that φ aboveand below 0.5 tend to yield positive and negative velocities. There are some exceptions,such as Panel C, where φ > . φ for which there is only 1 solution is much greater relative to the rangeof φ 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 thetug-of-war effect is more likely to be seen for lower values of viscous drag. We will laterdiscuss robustness in a similar way, where relatively greater parameter ranges correspondto increased robustness.In Figures 3E, F, we fix φ at 0.48 and 0.49, respectively, and vary ζ . These bifurca-tion diagrams are less intuitive, but nevertheless can be understood as slices through thebifurcation surface described above. More importantly, these bifurcation diagrams revealsome information about the underlying critical manifold. In 3E, we see that if the vesiclehas a sufficiently positive initial velocity at low drag, it maintains a positive velocity asthe vesicle moves into the spine neck, and the drag grows due to the constriction. At acritical drag denoted by the orange star, the vesicle instantaneously switches velocity inthe opposite direction by jumping down to the lower stable branch and eventually exitsthe constriction through the base of spine at − µ m. Similar discontinuous behavior isseen in 3F: with sufficiently positive initial velocity, the vesicle initially moves towards theconstriction before jumping down to the middle stable branch where the velocity is nearzero. In this scenario, the vesicle remains stuck for long times. Note that by symmetryof these functions, it follows that the bifurcation diagrams look identical with a change ofsign for φ = 0 . , .
51. Therefore, the model predicts that motor-driven transport through11onstrictions will generally push the vesicle towards the spine head as long as the ini-tial condition is sufficiently far into the constriction and that upwards-pushing motors aredominant.While one-parameter diagrams are useful, we wish to understand how multistabilitychanges in the entire φ - ζ parameter space. We address the question of multistability bynoting how each one-parameter bifurcation changes as a function of an additional param-eter. This process naturally partitions the φ - ζ parameter space into multistable regions.Our one-dimensional system only produces saddle-node (SN) bifurcations, which produceor destroy pairs of fixed points. We can therefore track the number of fixed points through-out the parameter space by tracking these saddle-node bifurcations. This process yieldsFigure 3D, in which we suppress fixed points and only display bifurcation points and thenumber of fixed points in each region. Each of the four colored curves correspond to asaddle-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 φ increases, the total number of fixed points changes in the order1 → → → →
1. The same can be observed in panel A by tracing the slice at ζ = 1 . × − Kg/s.The two-parameter diagram completely characterizes the total number of fixed pointsfor each region of the parameter space and shows that we can expect multistability inmuch of the displayed parameter space. If viscous drag is sufficiently small (in the bottomregion of panel D), multistability exists for a wide range of motor ratios φ . As viscousdrag increases, the range of multistability becomes much smaller as fixed points disappearthrough saddle-node bifurcations. For sufficiently large viscous drag ζ , multistability ceasesto exist as the saddle-nodes disappear through cusps, and there exists only one stablevelocity for all motor ratios. The critical vesicle-to-spine diameter ratio at the cusps isroughly 2 µ m/3 µ m.In terms of the bifurcation surface, the parameter slices in panels A–C show that forseveral values of fixed ζ , the bifurcation surface contains four folds. By choosing greatervalues of ζ , we find that the folds of the surface eventually flatten out through a pair ofcusp bifurcations. Beyond this cusp bifurcation, for greater values of ζ , we expect that thevelocities U are negative when φ < .
5, positive when φ > .
5, and zero when φ = 0 . π , π The two-parameter figure in φ - ζ provides a complete description of how multistabilitychanges as a function of motor ratio and constriction geometry when motor parameters π − π are fixed. However, there may be variations in motor parameters due to theexistence of multiple motor types such as myosin V and VI [10] and variations in ATP and12igure 4: Cusp bifurcations as a function of π and π . A: For each ζ >
0, there existcusp bifurcations along a set of π , π . Example level curves are plotted for ζ = 0 , . × − , . × − , . × − . We take 4 representative pairs of π , π labeled B–E and showthe corresponding two-parameter bifurcation diagrams in B–E. The point labeled with (cid:63) corresponds to Figure 3.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 parameterspace 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. Fora given ζ , it is possible to track the location of these cusps as a function of π and π . Theresult of this process for various π and ζ is shown in Figure 4A. Each curve represents acusp bifurcation as a function of π , π for a given ζ (we briefly describe how we determinedthe location of these cusps in Appendix A.3).13iven parameters π , π chosen somewhere in the π - π parameter region below the ζ = 0 level curve in Figure 4, the ζ level curves suggest that there exists a cusp bifur-cation at some ζ ∗ ( π , π ) >
0. Because cusps are a sufficient condition for the existenceof multistability, it follows that multistable states exist for appropriate choices of φ and ζ ≤ ζ ∗ . The ζ = 0 level curve represents a sufficient condition for the loss of multistability:noting that for each fixed point, ∂ζ/∂π < ζ decreases as π increases, it follows that for any π , π in at least a neighborhood abovethis curve, there can exist no cusp bifurcation with a positive ζ .Various points in Figure 4A are marked B–E, with corresponding two-parameter bifur-cation diagrams shown in the remaining subplots, Figures 4B–E. The point marked by (cid:63) inFigure 4A represents the π and π values for Figure 3. These diagrams provide evidencethat for smaller π and greater π , the region of multistability tends to increase. Recallingthat π is the maximum displacement of a motor in its non-preferred direction, smaller π implies that bidirectional motion due to noise can be made more likely by allowing themotor to detach earlier. Next, π is the initial motor attachment position in either thepreferred or non-preferred direction, and we find that the area of multistability increases asmotors have greater initial extension. Together, we predict that strong initial attachmentforces combined with greater yield effects can result in more frequent directional switching. The existence of a stationary solution U = 0 is straightforward to prove by inspectionfor φ = 0 .
5. In this section, we expand about this solution to determine the existenceof solutions when φ is near the equal motor ratio φ = 0 . ζ is large. We let φ = 0 . φ and ˆ ζ = 1 /ζ and explore the cases where ˆ φ and ˆ ζ are small. These limitsinform us of the local linear behavior of the velocity function U = C ( ˆ φ, ˆ ζ ) by writing C ( ˆ φ, ˆ ζ ) = ∂C∂ ˆ φ (cid:12)(cid:12)(cid:12)(cid:12) (0 , ˆ φ + ∂C∂ ˆ ζ (cid:12)(cid:12)(cid:12)(cid:12) (0 , ˆ ζ + O ( ˆ φ ˆ ζ, ˆ φ , ˆ ζ ) . First we derive an equation for the small deviation φ = 0 . φ , where 0 < | ˆ φ | (cid:28) U = C must satisfy0 = φF − A ( C ) + (1 − φ ) F A ( C ) − ζC = 12 F − A ( C ) + 12 F A ( C ) + ˆ φ [ F − A ( C ) − F A ( C )] − ζC. Solving for ˆ φ yields, ˆ φ = 12 F − A ( C ) + F A ( C ) − ζCF A ( C ) − F − A ( C ) . < − e π and 0 < π ,and ζ >
0. It follows that the derivative of ˆ φ with respect to C is nonzero: d ˆ φdC (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) C =0 = 12 e π π + ζe π − (cid:54) = 0 . We then obtain a local equation for C as a function of ˆ φ by invoking the inverse functiontheorem: C ( ˆ φ ) = 2 ˆ φ e π − e π π + ζ + O ( ˆ φ ) . (20)We have seen that in some parameter regimes, viscous drag grows to large valuesduring translocation (Figure 2C,F), but nonzero velocities persist due to unequal motorratios (Figure 1C,F). To establish this inverse relationship between velocity and drag, weTaylor expand about infinity with ˆ ζ = 1 /ζ as defined above. Then solving for ˆ ζ in termsof C yields ˆ ζ = CφF − A ( C ) + (1 − φ ) F A ( C ) . To derive a local equation for C as a function of ˆ ζ , we examine the derivative of ˆ ζ withrespect to C : d ˆ ζdC (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) C =0 = 1( e π −
1) [2 φ − . So long as φ (cid:54) = 1 / π >
0, this derivative is well-defined, and we can invoke the inversefunction theorem to write the local equation for C as a function of ˆ ζ . C (ˆ ζ ) = ˆ ζ ( e π −
1) [2 φ −
1] + O (ˆ ζ ) . After a trivial substitution we arrive at the desired equation, C ( ζ ) = ζ − ( e π −
1) [2 φ −
1] + O (cid:0) ζ − (cid:1) . Combining the local velocity estimates, we arrive at the velocity equation as a function ofˆ φ and ζ : C ( ˆ φ, ζ ) = 2 ˆ φ ( e π − (cid:104) ( e π π + ζ ) − + ζ − (cid:105) + O (cid:0) ζ − (cid:1) + O ( ˆ φζ − , ˆ φ , ζ − ) . (21)This existence equation is valid for large ζ and any ˆ φ ∈ [ − / , / φ and any ζ . Equation (21) provides justification for the effects observed from simulations earlierin the paper. For ζ large, the velocity must be small and proportional to ζ − , and thesign of the velocity is determined by the sign of ˆ φ , that is, it is the ratio of motors thatdetermines the direction of motion for large drag. In the case of small drag, this equation15s only valid for small φ , but the equation yields the same intuition that is consistent withour one-parameter bifurcation diagrams, i.e., that the sign of the velocity is determinedby the dominant motor species. Finally, we can make additional predictions. In thecase of large drag, or small drag with near-equal (but non-equal) motor ratios, increasingmotor attachment position and increasing the velocity scale ratio between translocationand motor adhesion dynamics ( π and π ) will increase the velocity of translocation. Inaddition, velocity does not locally depend on the ratio of motor attachment and detachmentparameters and the maximum displacement of each motor ( π or π ). For our one-dimensional problem, a linear stability analysis is sufficient to understand thestability of fixed points, which follows from the slope of the total force function: λ = F (cid:48) ( C ) − ζ. (22)Note that if the derivative of F is bounded, it follows that the eigenvalue is negative ( λ < λ shown above. By using continuity of the force-velocity curves (Appendix B.1), theproof follows by contradiction: if there are only two fixed points that are both stable, theslope of the net force function must be negative at each point. By the intermediate valuetheorem, there must exist a third point between them, which contradicts the original claimof 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 caseof large drag.In the special case C = 0, the eigenvalue can be computed explicitly to yield λ = e π π − e π − ζ. In the physically relevant parameter regime, the conditions π , π > ζ ≥ λ < φ = 0 .
5, we use Equation (20) to rewrite the eigenvalue: λ = F (cid:48) (cid:18) φ e π − e π π + O ( ˆ φ ) (cid:19) − ζ = F (cid:48) (0) + 2 ˆ φ e π − e π π F (cid:48)(cid:48) (0) − ζ + O ( ˆ φ )= e π π − e π + 2 ˆ φ e π − e π π (cid:20) e π (2 φ − π e π − (cid:21) − ζ + O ( ˆ φ ) . φ = 0 . φ , the term 2 φ − φ , making the secondterm order O ( ˆ φ ). Therefore, the eigenvalue equation reduces to λ = e π π − e π − ζ + O ( ˆ φ ) . Recalling that π , π >
0, we generally expect constant velocity solutions to be stable forˆ φ small. In this paper, we fully characterize the dynamics predicted by a model of vesicles driveninto closed constrictions. We cast the system into a fast-slow system and perform a two-parameter bifurcation analysis on the fast subsystem, and then determine how the cusps inthe resulting bifurcation surface vary with the motor parameters. The model predicts mul-tistability (i.e. bidirectional motion) for smaller values of viscous drag and unidirectionalmotion for greater values of the viscous drag corresponding to tight constrictions.We remark that while dendritic spines exhibit a great diversity of morphologies such asthin 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], mostcan be captured by our reduced axisymmetric lubrication model. Moreover, while vesiclesof recycling endosomes range in diameter from small to large (roughly in the range from0.5 µ m to 2 µ m [10]) and deform strongly in the spine neck, the most important requirementfor us to apply the theory of this paper is that the diameter of the vesicle and spine wallare similar. While detailed electron microscopy images exist of dendritic spines [24], it isa significant challenge to efficiently search and classify recycling endosomes. Therefore,for now we rely on published endosome images to approximate the physiologically relevantparameter ranges.In Figure 5A,B we show two representative experimental images from [10] of spinescontaining recycling endosomes. The spine is outlined by a thin yellow line and the en-dosome is shown traveling through the spine neck in a series of four time-lapse imageswith an associated kymograph on the right. The scale bar is 2 µ m in all panels. Usingthese images, we determine approximate regions where physiological parameters may liein the φ − ζ parameter space (Figures 3D and 5C). Figure 5C is the same as Figure 3Dbut with ζ plotted on a log scale. The time-lapse images provide virtually no informationabout the ratio of motors, so we make no restrictions on φ for now.We estimate the viscous values from these images by estimating the height betweenthe 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 most0.1 µ m. This height is substantially smaller than the heights considered in this paper, andtherefore yields a relatively greater viscous drag value of ζ ≈ × − Kg/s (we assume17
B AB , Neck B , Head C Figure 5: Time-lapse images of recycling endosomes adapted from [10] and available underthe CC BY NC ND license. A: a recycling endosome translocates through a thin spinein a series of four time-lapse images (left). A kymograph is shown to the right. Thevesicle is roughly 1 µ m in diameter, and the distance between the vesicle and neck wallis at most 0.1 µ m. B: a recycling endosome translocates into a stubby spine in a seriesof time-lapse images, with the associated kymograph on the right. The vesicle is roughly0 . µ m in diameter, the distance between the vesicle and neck wall is roughly 0.15 µ m, andthe 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 diagramfrom Figure 3 with the drag ζ plotted on a log scale. Labels A and B correspond to PanelsA and B.the same microscopic motor parameters as in Figure 3). This value corresponds to a pointfar above the cusps, where only a single-velocity solution exists (Figure 5C, top red bandlabeled “A”). This range lies well within the unidirectional regime. Because the time-lapseimages in Figure 5A show the endosome traveling towards the spine head, we infer thatthe upwards-pushing motors are dominant. The only region where the vesicle could switchdirection is at the base of the spine, where there some distance between the diameter of thecell wall and the recycling endosome is possible. Thin spines generally are much smallerin diameter relative to recycling endosomes, so that the distance between the vesicle andspine wall is small. We conclude that once the vesicle has entered a thin spine, we generallyexpect unidirectional movement.Figure 5B shows an example of a stubby spine. Here, the endosome is smaller, witha diameter range from 0 . µ m to 0 . µ m. The lower bound of the height range is difficultto establish, but we estimate that the height between the endosome and neck is roughly0.05–0.1 µ m, which yields viscous drag values in the range ζ ≈ × − to 6 × − Kg/s.Recalling that the cusps in Figure 3 are on the order of 2 × − Kg/s (where the criticalvesicle-to-spine diameter ratio is roughly 2 µ m/3 µ m), we find that both unidirectional andmultidirectional motion are possible (Figure 5C, middle red band labeled “B, neck”), where18nidirectional solutions take a much greater portion of the parameter range in a linearscale. Another property to consider is that stubby spines have shallow constrictions thatlead into a large head, where the vesicle may spend a substantial amount of time relativeto the neck. It is therefore worth considering whether multistable solutions exist pastthe neck. In the head, the height between the endosome and head ranges from roughly0.1 µ m to 0.5 µ m, which yields a drag range of 3 × − Kg/s to 1 × − Kg/s, whichplaces the vesicle squarely in a multistable regime (Figure 5C, lower red band labeled “B,head”). Strikingly, the kymograph shows unidirectional movement through the neck, andbidirectional motion in the spine head.The viscous drag in these representative examples varies over four orders of magnitudeand our model predicts dramatically different qualitative behaviors over this range. Ac-cording to our model, thin spines with large vesicles will exhibit unidirectional movementdue to the large drag experienced by the vesicle. In contrast, stubby spines will exhibit mul-tidirectional 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 force-velocity 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 possibleapproach to overcome this limitation is to replace the current mean-field model of molecularmotors by a discrete model; such a discrete model is used to compute mean passage timesin Allard et al [2]. Finally, numerical experiments are another feasible approach to addressthe question of mean first passage times of a translocating vesicle. Similar approaches havebeen performed in [11] to pursue questions on motor and intermediate filament parametereffects on intermediate filament transport.
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 lubricationmodel. 19 .1 Integration
For convenience, we restate the original fast-slow lubrication model here: dZdt = U,ε dUdt = F ( U ) − ζ ( Z ) U. (23)We solve Equation (23) numerically by taking ε nonzero and integrating forwards intime (as we have done using the forward Euler method). This approach is beneficial becausevelocities 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 ofconvergence is then very simple: if an initial condition is within the basin of attraction ofa stable velocity, the initial condition will converge to this velocity.Another benefit of taking this approach is that for ε sufficiently small, system (23) issufficiently similar to the overdamped system when ε = 0, as guaranteed by Fenichel theory[17, 8]). When ε (cid:28)
1, there is a separation of timescales so that very small times stepsare required for numerical stability. We find that choosing ε = 1 works well in practice, asthe 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, whereaswe 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 ε bythe effective drag. For a given choice of ε , we typically take the time step dt to be smallerby one or two orders of magnitude. A.2 Continuation
General continuation strategies in one and two parameters can be found in Chapter 10.3 ofKuznetsov [27]. We use XPPAUTO [15] (version 8) to generate our bifurcation diagramsunless stated otherwise. For the generation of the one- and two-parameter diagrams, weuse the following numerical values:
Ntst = 15 (default)Nmax = 200 (default)Npr = 50 (default)Ds = 1e-10Dsmin = 1e-10Ncol = 4 (default)Dsmax = 0.1
All other numerical values remain at default values. We adjust
Ntst , Nmax , Npr , Par Min ,and
Par Max as needed (our mini tutorial on the repository explains what parameters are20ppropriate). The displayed
Dsmax value may make AUTO run multiple passes over somebranches, in which case we reduce
Dsmax to 0.001 preventing this issue while incrementingat a reasonable pace. We include many more details in a mini tutorial on our githubrepository 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 asaddle-node bifurcation are [27]: f ( U ) = F ( U, φ ) − ζU = 0 ,f (cid:48) ( U ) = F (cid:48) ( U, φ ) − ζ = 0 ,f (cid:48)(cid:48) ( U ) = F (cid:48)(cid:48) ( U ) (cid:54) = 0 . (24)We denote any U that simultaneously satisfy these equations by U ∗ . We can simplify thesearch for saddle-nodes by writing the system¯ U (cid:48) = F ( ¯ U , ¯ φ ) − ζ ¯ U . ¯ φ (cid:48) = F (cid:48) ( ¯ U , ¯ φ ) − ζ. (25)We can just as easily use ¯ ζ and look for fixed points in ( ¯ U , ¯ ζ ) as a function of φ . It issimply a matter of preference. We use ¯ U and ¯ φ to emphasize the difference between (25)and the original fast subsystem (18). This new system eliminates one parameter whenit comes to computing bifurcation diagrams. Fixed points of (25) correspond to saddle-node bifurcations in the fast subsystem (18). The ¯ U nullcline of (25) shows the stablefixed points in (18), and the ¯ U nullcline intersections with the ¯ φ nullcline correspond tosaddle-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 bifurcationsin (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 wherethe intersections of these nullclines change in number. Suppose that we are interested infinding cusp bifurcations as a function of π and π for a given ζ . On one side of the cuspbifurcation, Equation (25) shows two fixed points and on the other side shows zero. Bydefining a function that returns +1 on one side of the bifurcation and − π where there is a transitionin 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 thissection we briefly show some properties used in the text.21 .1 Continuity
The continuity of the force-velocity curve (Equation (12)) is not immediately obvious. Werewrite the force-velocity curve verbatim for convenience: (cid:101) F A = − π (cid:101) U ( e π − − − π (cid:101) U , if (cid:101) U < − ( π +1) π (1 − e − π /π (cid:101) U )+1 [ e π (1 − e π e − π /π (cid:101) U )] − (1 − π (cid:101) U )(1 − e − π /π (cid:101) U )( e π − − π (cid:101) U ) , if (cid:101) U ≥ . (26)In particular, there are two possible problem areas: when (cid:101) U = 0 and when (cid:101) U = 1 /π . Thelatter case is especially important to consider because numerical problems occasionallyarise when evaluating (cid:101) U = 1 /π directly. In the first case, the left limit is straightforward:lim (cid:101) U → − (cid:101) F A = − . (27)The right limit depends on the behavior of the exponentials. Noting that exp( − π /π (cid:101) U ) → (cid:101) U → + , it is straightforward to check that the left and right limits agree independentof the parameters π i , and therefore (cid:101) F A is continuous at (cid:101) U = 0. Continuity of (cid:101) F − A followsby definition.In the second case, when (cid:101) U = 1 /π ≥
0, we take a close look at the second line ofEquation (26). Potential problems arise in the term (1 − e π e − π /π (cid:101) U ) / (1 − π (cid:101) U ). Forconvenience, let v := 1 − / ( π (cid:101) U ) so that v → (cid:101) U → /π . Then the term becomes1 − e π e − π /π (cid:101) U )1 − π (cid:101) U = − (1 − e π v )(1 − v ) v = − [1 − (1 + π v + O ( v ))](1 − v ) v = − [1 − (1 + v ( π −
1) + O ( v ))] v = v ( π −
1) + O ( v )) v = π − O ( v ) . We have used Taylor expansion of the exponential about zero. So the term converges to π − v → (cid:101) F A is continuous at (cid:101) U = 1 /π , as desired. B.2 Limits
Consider π →
0. This limit occurs when attachments occur at zero position ( A = 0), orwhen 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 in22he latter case, the force exerted by a single motor is always zero. Indeed, we find thatthe term e π − U when U is small. For | U | (cid:29)
0, the term ( e π − − dominates, and we expect F A ( U ) → ∞ as π →
0. When | U | (cid:28)
1, we expect F A ( U ) → − U → π small. F A ( U, π →
0) = ∞ U < , − U = 0 , −∞ U ≥ . In contrast, the dimensional force-velocity curve converges for π small and diverges for π . A nondimensionalization using F = F / ( e π −
1) would remove this problem of di-vergence, but this rescaling is not necessary when considering bifurcations: for any fixedpoint F ( U, π ) = 0, where F is the nondimensional net force, multiplying both sides bythe scaling factor ( e π −
1) yields ( e π − F ( U, π ) = 0, so in fact scaling preserves fixedpoints as a function of π . References [1] David J Acheson. Elementary fluid dynamics, 1991.[2] Jun Allard, Marie Doumic, Alex Mogilner, and Dietmar Oelz. Bidirectional sliding oftwo parallel microtubules generated by multiple identical motors.
Journal of mathe-matical biology , pages 1–24, 2019.[3] Jon I Arellano, Ruth Benavides-Piccione, Javier DeFelipe, and Rafael Yuste. Ultra-structure of dendritic spines: correlation between synaptic and spine morphologies.
Frontiers in neuroscience , 1:10, 2007.[4] Josephine Shaw Bagnall, Sangwon Byun, Shahinoor Begum, David T Miyamoto, Vi-vian C Hecht, Shyamala Maheswaran, Shannon L Stott, Mehmet Toner, Richard OHynes, and Scott R Manalis. Deformability of tumor cells versus blood cells.
Scientificreports , 5:18542, 2015.[5] Daniel R Berger, H Sebastian Seung, and Jeff W Lichtman. Vast (volume annotationand segmentation tool): efficient manual and semi-automatic labeling of large 3d imagestacks.
Frontiers in neural circuits , 12:88, 2018.[6] Paul Bressloff and Jay Newby. Directed intermittent search for hidden targets.
NewJournal of Physics , 11(2):023033, 2009.[7] Paul C Bressloff and Jay M Newby. Metastability in a stochastic neural networkmodeled as a velocity jump markov process.
SIAM Journal on Applied DynamicalSystems , 12(3):1394–1435, 2013. 238] Henk W Broer, Tasso J Kaper, and Martin Krupa. Geometric desingularization of acusp singularity in slow–fast systems with applications to zeemans examples.
Journalof Dynamics and Differential Equations , 25(4):925–958, 2013.[9] Sangwon Byun, Sungmin Son, Dario Amodei, Nathan Cermak, Josephine Shaw,Joon Ho Kang, Vivian C. Hecht, Monte M. Winslow, Tyler Jacks, Parag Mallick,and Scott R. Manalis. Characterizing deformability and surface friction of cancercells.
Proceedings of the National Academy of Sciences , 110(19):7580–7585, 2013.[10] Marta Esteves da Silva, Max Adrian, Philipp Sch¨atzle, Joanna Lipka, Takuya Watan-abe, Sukhee Cho, Kensuke Futai, Corette J Wierenga, Lukas C Kapitein, and Casper CHoogenraad. Positioning of ampa receptor-containing endosomes regulates synapse ar-chitecture.
Cell reports , 13(5):933–943, 2015.[11] JC Dallon, C´ecile Leduc, Sandrine Etienne-Manneville, and St´ephanie Portet. Stochas-tic modeling reveals how motor protein and filament properties affect intermediatefilament transport.
Journal of theoretical biology , 464:132–148, 2019.[12] Geoffrey Dawson, Edgar H¨aner, and Anne Juel. Extreme deformation of capsules andbubbles flowing through a localised constriction.
Procedia IUTAM , 16:22–32, 2015.[13] Wynter J Duncanson, Thomas E Kodger, Sahab Babaee, Grant Gonzalez, David AWeitz, and Katia Bertoldi. Microfluidic fabrication and micromechanics of permeableand impermeable elastomeric microbubbles.
Langmuir , 31(11):3489–3493, 2015.[14] Guy N Elston. Specialization of the neocortical pyramidal cell during primate evolu-tion. 2007.[15] G. Bard Ermentrout.
Simulating, analyzing, and animating dynamical systems: aguide to XPPAUT for researchers and students , volume 14. SIAM, 2002.[16] Thomas G Fai, Remy Kusters, Jens Harting, Chris H Rycroft, and L Mahadevan.Active elastohydrodynamics of vesicles in narrow blind constrictions.
Physical ReviewFluids , 2(11):113601, 2017.[17] Neil Fenichel. Geometric singular perturbation theory for ordinary differential equa-tions.
Journal of differential equations , 31(1):53–98, 1979.[18] Sylvain Gabriele, Marie Versaevel, Pascal Preira, and Olivier Th´eodoly. A simplemicrofluidic method to select, isolate, and manipulate single-cells in mechanical andbiochemical assays.
Lab on a Chip , 10(11):1459–1467, 2010.[19] Edward G Gray. Axo-somatic and axo-dendritic synapses of the cerebral cortex: anelectron microscope study.
Journal of anatomy , 93(Pt 4):420, 1959.2420] Thomas Gu´erin, J Prost, and J-F Joanny. Motion reversal of molecular motor assem-blies due to weak noise.
Physical review letters , 106(6):068101, 2011.[21] Kristen M Harris and John K Stevens. Dendritic spines of ca 1 pyramidal cells inthe rat hippocampus: serial electron microscopy with reference to their biophysicalcharacteristics.
Journal of Neuroscience , 9(8):2982–2997, 1989.[22] Frank C Hoppensteadt and Charles S Peskin.
Modeling and simulation in medicineand the life sciences , volume 10. Springer Science & Business Media, 2012.[23] Frank J¨ulicher and Jacques Prost. Cooperative molecular motors.
Physical reviewletters , 75(13):2618, 1995.[24] Narayanan Kasthuri, Kenneth Jeffrey Hayworth, Daniel Raimund Berger, Richard LeeSchalek, Jos´e Angel Conchello, Seymour Knowles-Barley, Dongil Lee, Amelio V´azquez-Reina, Verena Kaynig, Thouis Raymond Jones, et al. Saturated reconstruction of avolume of neocortex.
Cell , 162(3):648–661, 2015.[25] Thomas J Koehnle and Anthony Brown. Slow axonal transport of neurofilamentprotein in cultured neurons.
The Journal of cell biology , 144(3):447–458, 1999.[26] Ambarish Kunwar, Suvranta K Tripathy, Jing Xu, Michelle K Mattson, PreethaAnand, Roby Sigua, Michael Vershinin, Richard J McKenney, C Yu Clare, Alexan-der Mogilner, et al. Mechanical stochastic tug-of-war models cannot explain bidi-rectional lipid-droplet transport.
Proceedings of the National Academy of Sciences ,108(47):18960–18965, 2011.[27] Yuri A Kuznetsov.
Elements of applied bifurcation theory , volume 112. SpringerScience & Business Media, 2013.[28] Yang Li, Ozan S Sarıyer, Arun Ramachandran, Sergey Panyukov, Michael Rubin-stein, and Eugenia Kumacheva. Universal behavior of hydrogels confined to narrowcapillaries.
Scientific reports , 5:17017, 2015.[29] CA Miermans, RPT Kusters, CC Hoogenraad, and C Storm. Biophysical model ofthe role of actin remodeling on dendritic spine morphology.
PloS one , 12(2):e0170113,2017.[30] Melanie JI M¨uller, Stefan Klumpp, and Reinhard Lipowsky. Tug-of-war as a coopera-tive mechanism for bidirectional cargo transport by molecular motors.
Proceedings ofthe National Academy of Sciences , 105(12):4609–4614, 2008.[31] Jay Newby and Paul C Bressloff. Random intermittent search and the tug-of-warmodel of motor-driven transport.
Journal of Statistical Mechanics: Theory and Ex-periment , 2010(04):P04014, 2010. 2532] Jay M Newby and Paul C Bressloff. Directed intermittent search for a hidden targeton a dendritic tree.
Physical Review E , 80(2):021913, 2009.[33] Jay M Newby and James P Keener. An asymptotic analysis of the spatially inho-mogeneous velocity-jump process.
Multiscale Modeling & Simulation , 9(2):735–765,2011.[34] Esther A Nimchinsky, Bernardo L Sabatini, and Karel Svoboda. Structure and func-tion of dendritic spines.
Annual review of physiology , 64(1):313–353, 2002.[35] Mikyoung Park, Jennifer M Salgado, Linnaea Ostroff, Thomas D Helton, Camenzind GRobinson, Kristen M Harris, and Michael D Ehlers. Plasticity-induced growth ofdendritic spines by exocytic trafficking from recycling endosomes.
Neuron , 52(5):817–830, 2006.[36] St´ephanie Portet, C´ecile Leduc, Sandrine Etienne-Manneville, and John Dallon. De-ciphering the transport of elastic filaments by antagonistic motor proteins.
PhysicalReview E , 99(4):042414, 2019.[37] W Christopher Risher, Tuna Ustunkaya, Jonnathan Singh Alvarado, and CaglaEroglu. Rapid golgi analysis method for efficient and unbiased classification of den-dritic spines.
PloS one , 9(9), 2014.[38] W. Christopher Risher, Tuna Ustunkaya, Jonnathan Singh Alvarado, and CaglaEroglu. Rapid golgi analysis method for efficient and unbiased classification of den-dritic spines.
PLOS ONE , 9(9):1–8, 09 2014.[39] Cynthia L Walker, Atsuko Uchida, Yinyun Li, Niraj Trivedi, J Daniel Fenn, Paula CMonsma, Roxanne C Larivi´ere, Jean-Pierre Julien, Peter Jung, and Anthony Brown.Local acceleration of neurofilament transport at nodes of ranvier.
Journal of Neuro-science , 39(4):663–677, 2019.[40] Zhiping Wang, Jeffrey G. Edwards, Nathan Riley, D. William Provance, Ryan Karcher,Xiang dong Li, Ian G. Davison, Mitsuo Ikebe, John A. Mercer, Julie A. Kauer, andMichael D. Ehlers. Myosin vb mobilizes recycling endosomes and ampa receptors forpostsynaptic plasticity.
Cell , 135(3):535 – 548, 2008.[41] Rafael Yuste.
Dendritic spines . MIT press, 2010.[42] Dennis Zimmermann, Alicja Santos, David R Kovar, and Ronald S Rock. Actin ageorchestrates myosin-5 and myosin-6 run lengths.
Current Biology , 25(15):2057–2062,2015.[43] Herbert Zimmermann. Accumulation of synaptic vesicle proteins and cytoskeletalspecializations at the peripheral node of ranvier.