Bistability, oscillations and bidirectional motion of ensemble of hydrodynamically-coupled molecular motors

We analyze the collective behavior of hydrodynamically coupled molecular motors. We show that the local fluxes induced by motors displacement can induce the experimentally observed bidirectional motion of cargoes and vesicles. By means of a mean--field approach we show that sustained oscillations as well as bistable collective motor motion arise even for very large collection of motors, when thermal noise is irrelevant. The analysis clarifies the physical mechanisms responsible for such dynamics by identifying the relevant coupling parameter and its dependence on the geometry of the hydrodynamic coupling as well as on system size. We quantify the phase diagram for the different phases that characterize the collective motion of hydrodynamically coupled motors and show that sustained oscillations can be reached for biologically relevant parameters, hence demonstrating the relevance of hydrodynamic interactions in intracellular transport.

The transport of macromolecules, vesicles and organelles inside cells relies on the active motion of molecular motors along biofilaments [1].When motors pull on organelles or vesicles, the fluid flow that they induce provides motor-motor hydrodynamic interactions (HI) hence providing an additional mechanism for molecular motors coupling [2] that, possibly, is responsible of the experimentally observed enhanced velocity of motors pulling on fluid vesicles [3].In particular, the motion of these cargoes has been observed to be mono-and bi-directional, the latter relying on the presence of motors pulling the cargoes in opposite directions [4][5][6][7][8][9][10].While previous studies have focused preferentially on rigidly coupled motors [11][12][13][14], in this letter we show that the bidirectional motion observed in experiments can be induced and controlled via the HI induced by motors active displacement.Exploiting a mean-field approach we identify the key parameters controlling the onset of the bidirectional motion and we show the relevance of HI for biologically relevant scenarios.
While, a few motors at the cargo tips can pull against the cargo drag, the rest of the motors can move along the cargo.Their net motion is responsible for the onset of HI.We model the molecular motor dynamics exploiting the two-state-model [11] that regards molecular motors as particles with two internal states.In order to account for motors pulling in opposite directions, we model the two families of motors as a single family of effective motors that in the "bound" state experience a symmetric peri-odic force f (x) = −∂ x V = f 0 cos(2πx/L) = f 0 f (x) of period L [14].With rate ω of f (x) motors jump to a "weakly bound" state, in which they diffuse freely.Motors bind with rate ω on (x).We describe the system in terms of densities of bound (ρ(x)) and weakly bound (σ(x)) motors.The total density, ρ(x) + σ(x), cannot exceed the maximal filament occupation prescribed by excluded volume, and in order to keep analytical insight, we assume that the motors remain in a dilute regime.The motion of molecular motors occurs in the low Reynolds regime and it generates a fluid flow whose magnitude reads γv(x) = f (x) + f (y)ρ(y)W (x, y)dy (1) where W (x, y) is the dimensionless Oseen tensor accounting for HI [23] and γ is the single motor drag coefficient.According to Eq. ( 1), the induced fluid velocity is characterized by the local flow due to the force that the filament exerts on the molecular motor (first term in rhs of Eq. 1) and the collective hydrodynamic flow induced by the rest of the molecular motors at the position of the reference motor.This separation between local and collective induced flows is characteristic of softly interacting motors, and it is absent for rigidly coupled motors [14].Since the HI extends over distances large compared to the spatial variation of the motor-filament interaction, for a large system size Λ, Λ ≫ L, we assume a mean field HI which holds when f (x) = f 0 and ρ(x) = ρ 0 , and becomes exact for an infinite system, Λ → ∞, if the motor density follows the spatial dependence imposed by the motor/filament periodic force.In this regime, reasonable when boundary effects can be disregarded, the spatial dependence of the hydrodynamic long-range coupling, W (x, y), experienced by the motors becomes negligible and it can be described by an effective dimensionless parameter, k.For example, for motors of linear size R moving in a three-dimensional environment, which grows logarithmically with system size.Alternatively, when motors are pulling on a membrane-coated cargo, such as organelles or vesicles, their tails are linked to molecules embedded in the fluid-like membrane (characterized by a 2D viscosity η 2D,mem ) quite more viscous than the cytoplasm (characterized by a 3D viscosity η 3D,cyt ).When the intrinsic dynamics of motor tails does not affect motors dynamics, we can identify the dynamics of the motors with that of the tracers and therefore we can calculate the 2D HI.In this regime, we expect the 2D HI to dominate the three-dimensional flows induced in the cytoplasm for motor-motor separations smaller than ) where 2 ln l r is the dimensionless Oseen tensor in 2D and M = min l e , Λ 2 .Eqs. (3),(4) capture the diverging nature of HI for increasing system sizes Λ.Therefore, for large system sizes, Λ/L ≫ 1, we can neglect the O(1) terms in Eq. (2) and Eq.(1) reduces to where we have introduced the dimensionless density ρ(x) = Lρ(x) (and similarly we introduce σ(x) = Lσ(x)).
Accordingly, we can write, in dimensionless units, the evolution of the mean bound, ρ(x) = ρ(x) ρ,σ , and weakly bound, σ(x) = σ(x) ρ,σ , motor densities [24] ρ where ψ(x) ρ = ψP [ρ(x)]dx stands for the average of ψ over the probability distribution of density profiles ρ(x) while ψ(x) x = (1/L) L 0 ψ(x)dx corresponds to the spatial average of ψ over the filament period (see Suppl.Mat).The first term in the rhs of Eqs. ( 6) describes motor advection according to the local velocity the motors are exposed to, the second term describes the motor binding kinetics while the third term in the evolution equation for weakly bound motors describes the fact that motors can detach from the filament with rate k off and bind with rate k on , proportional to the available space c ∞ − ρ(x) − σ(x) where c ∞ is the motor concentration in bulk.[25].
The motor flux in Eqs. ( 6) is proportional to λ = f 0 /ωLγ, the ratio of the typical time a motor needs to slide down the potential, Lγ/f 0 , and the characteristic inverse hopping time ω = ω on (x) + ω of f (x) x /2.The bound and weakly bound rate densities depend on the filament structure.We consider the simple, periodic form ωon/off = max(∆ω on/of f ∓ δω on/of f sin(2πx/L), 0) (7) which interpolates between regimes where the hopping rates are essentially homogeneous along the filament, ∆ω on/of f ≫ δω on/of f , or only take place at localized regions on the filament, ∆ω on/of f ≪ δω on/of f .The latter, together with an opposite sign between the bound and weakly bound rates accounts for the fact that binding and unbinding processes are localized at different regions along the filament.Normally, unbinding rates are more spatially localized than bounding rates, implying δω of f /∆ω of f > δω on /∆ω on .
For fast bulk kinetics k on,off → ∞, the concentration of weakly bound motors, σ is homogeneous along the filament.In this regime, Eqs.(6) decouple and the collective behavior of the motors is controlled by the evolution of ρ(x), in terms of rescaled binding rates that keep the same functional dependence as in Eq. ( 7).We have numerically solved Eq. ( 8) using a Lax-Wendroff scheme with periodic boundary conditions applied at the ends of a period of the ratchet potential.Fig. 1.(a) shows that the configuration where motors do not have a net velocity, stable for weak HI, becomes unstable at a critical coupling, k o , above which a net motor current, that breaks left/right symmetry, develops; a similar scenario has been described for rigidly coupled motors [14].Despite the apparent similarity between the emergence of net motion for soft and rigidly coupled motors, the different underlying physical mechanisms responsible for symmetry braking lead to significant differences in collective motor dynamics.While hydrodynamically coupled motors are characterized by a non-monotonous dependence of k o on λ, see as the sufficient condition [27].for symmetry breaking and onset of net motor currents.Since k ≥ 0, Eq. ( 9) identifies an interval λǫ[0; λ max ] for which symmetry breaking occurs, with λ max = 1 π δω ∆ω .According to Eq. ( 9), k diverges for both λ = 0, λ = λ max leading to the existence of a minimum value k = k c at λ = λ c = 1 2 λ max .For the parameters used in Fig. 1.(a) the stability analysis predicts λ c ∼ 8 • 10 −2 and k c 70 in good agreement with numerical results.For rigidly bound motors, the second term in the denominator of Eq. ( 9) disappears, leading to an inverse proportionality between the coupling constant k and the dimensionless forcing λ, consistent with the results shown in Fig. 1.(b).This difference has significant implications.For example,while for λ > λ c a reduction in λ will favor the onset of net fluxes, when λ < λ c decreasing λ will hinder, or even prevent, the development of a net motor flux.
In the opposite regime, when the exchange of molecular motors with the bulk is negligible, k on,off = 0, the total number of motors moving along a filament is conserved.When the evolution of bound, ρ(x), and weakly bound, σ(x), motors are coupled, Eqs. ( 6) must be solved consistently.In this regime motors tend to accumulate spatially, leading, in some cases, to large local motor densities.It is known that conservation of the overall number of motors moving along a filament promotes cluster [2] and shock wave formation [15].The morphological details of these structures are sensitive to excluded volume and short range interactions.However, for binding and unbinding rates sharply peaked at the potential extrema, ∆ω on,of f < δω on,of f , the development of regions of high molecular motor density only occurs for large coupling parameters.Hence, we can address the instability of the homogeneous, quiescent molecular motor profile avoiding the development of shock waves.We observe that this quiescent configuration destabilizes above a threshold coupling parameter, k 1 , characterized by a Hopf bifurcation, as shown in Fig 2 .(a)(for k 1 = 55).Above k 1 the stable state is characterized by a non-zero mean velocity and an oscillation of frequency, ω v , as shown in Fig. 2.(a) [28].Motor velocity oscillations emerge as a result of the periodic change in the density of the bound motors.While moving under the action of the driving potential, the fluid flow generated by bound motors advects weakly bound motors along.After reaching the bottom of the potential, bound motors cease to move and jump to the diffusive state with rate ω of f .This leads to an increase of diffusive motors, hence inducing an overall decrease in the average motor velocity, which relies in the small fraction of bound motors still displacing.Once diffusing motors reach the hopping region, they bind strongly to the filament at a rate ω on , starting a new cycle.As a result of this alternate, correlated motor exchange, both strongly, ρ(x), and weakly bound, σ(x), motor densities develop traveling waves.
As we increase k/λ a second bifurcation is observed, k 2 , above which molecular motors exhibit bistability.In this regime the motor density increases gradually where the motor states can be switched, leading eventually to an overall motor density that exceeds the maximum occupancy.In order to explore this regime we then redistribute uniformly the motor excess [29].
For the system shown in Fig. 2, above the threshold value k 2 = 95 the motor velocity still oscillates with frequency ω v around a non vanishing mean velocity, v.However, in contrast to the regime k 1 < k < k 2 , where the sign of v is fixed, for k > k 2 the sign of v changes with frequency Ω v .For time scales larger than Ω −1 v , the average motor velocity vanishes and a bistable behavior emerges (see Fig. 1 in Suppl.Mat.), analogous to the one experimentally observed [4][5][6][7][8][9][10]16].This second transition is captured by the dimensionless time, τ (k) = 2v(k)/ω 0 ǫ, defined as the ratio between the characteristic hopping time 1  2 ω −1 0 [30], and the time a particle spends in the region in which the hopping rate is non vanishing, ǫ, being pushed at a speed v(k) [31].As shown in Fig. 2.(b), when τ (k) 0.2 the time needed to jump between the two states is much smaller than the time that a particle spends in crossing the hopping region, therefore the majority of the motors rebind and the system undergoes sustained oscillations.On the contrary, when τ (k) 0.2 part of the motors in the weakly bound state cannot jump back to the bound state and do not contribute to the next cycle.The loss of active motors affects subsequent oscillations. Thee effects sum up until the system switches the direction of the average velocity on time scales of the order of Ω v .The number of oscillations between two subsequent switching events, N ω = 2ω v /Ω v , decreases for increasing k, as shown in the inset of Fig. 2.(a).For increasing k the two time scales approach, N ω → 1, and the bistability disappears; subsequently motors remain in a quiescent state.
The bistable state observed is typical of soft, hydrodynamically coupled motors for which the local density of both bound and weakly bound motors can adjust dynamically.This feature is absent both if we disregard weakly bound motor dynamics (see Fig. 1), or for rigidly coupled motors [14], when motor density rearrangements are suppressed.Bistability can be recovered for rigidly coupled motors [14], also when weakly bound motors are in contact with a reservoir (data not shown), only if the hopping dynamics is noisy.However, this mechanism vanishes for large system sizes, for which the noise becomes negligible.On the contrary, for hydrodynamically coupled motors bistability arises by increasing the system size, encoded in k, and persists for a finite range of values of k as shown in Fig. 2.(a), regardless of system size.
The collective phases identified for hydrodynamically  we obtain l = η 2D,mem /η 3D,cyt ∈ [10 −1 , 10] µm that lies within the typical range of biological situations for which the typical velocity is v ∼ 0.1 µm [19].The hopping rate can be assumed ω 0 ≃ 10 2 α s −1 , α being the inverse of the efficiency [32].For these values of parameters we get τ ≃ 0.2 that fits in the range of values identified by Fig. 2.(b).By inverting Eq. ( 4) we can calculate the dependence of Λ upon k when the 2D contribution dominates over the 3D.Fig. 3.(a) shows that for l ∼ 1µm, systems as small as Λ ∼ 0.5 µm can undergo hydrodynamically-induced symmetry breaking, while slightly larger systems Λ ∼ 0.7 µm develop bistability.For larger values of l the 3D contribution dominates and we expect hydrodynamically-induced symmetry breaking for systems of order of 10 − 100 µm emphasizing the relevance of the 3D HI for larger systems such as neurons, or in technological applications as in microfluidic devices.The onset of symmetry breaking and bistability depends on the parameters governing motors dynamics, namely, the force motor can provide, f 0 and their hopping rate ω.The ratio of f 0 and ω is encoded in the dimensionless parameter, λ.As shown in Fig. 3.(b) the values of both k 1 and k 2 decrease upon increasing the strength of the motors or decreasing their hopping rate the latter being easily controlled in experiments by tuning the ATP concentration.
In conclusion, the HI between bidirectional molecular motors strongly affects their dynamics.For intermediate system sizes, HI triggers the onset of net motor currents [33].For these regimes the motor velocity has been observed to be oscillatory, about a non vanishing aver-age, or bistable, when motors switch their direction over larger time scales, leading to an overall vanishing currents, as observed experimentally [5,[8][9][10].Such features rely on the local variation of the motor density, typical of HI, absent when this degree of freedom is neglected (rigid coupling) or reduced, as happens when the molecular motors are not conserved, e.g. through bath exchange.The typical system sizes over which the soft HI leads to symmetry breaking, or to bistability, are compatible with typical biologically relevant sizes, Λ ∈ [0.1, 10] µm, typical for Golgi apparatus displacement [8] or bistable cargo transport [5,9,10].P.M. acknowledges MINECO for financial support under projects FIS 2015-67837-P and Institute Curie for financial support.I.P. acknowledges MINECO and DURSI for financial support under projects FIS 2015-67837-P and 2014SGR-922, respectively, and from Generalitat de Catalunya under Program Icrea Acadèmia.ferent states of the same motor, their sum cannot exceed a maximum density, ζmax, corresponding to complete local occupancy.Accordingly, for non-dilute regimes, a global factor on the rhs of Eq. ( 6) must be included to ensure such a constraint.
[26] The rigid coupling regime is obtained by solving Eq. ( 8) disregarding the local contribution, f (x).
[27] Eq. ( 9) constitutes only a sufficient and not necessary condition since it has been obtained imposing that the dynamic system of linearized perturbations around the quiescent state has at least three real solutions without constraining their signs (see Suppl.Mat.) [28] For peaked hopping rates, the relevant time scale determining the oscillations is controlled by the hopping rates, ω0 rather than by their average values, ω [29] When exceeding the maximum value, ρmax, the excess of motors is equally redistributed along the period.
[30] The prefactor 1  2 accounts for the hopping rate variable shape.
[33] Such a coupling is amplified for non processive motors that are more prone to show collective effects due to their low duty ratio [20].

DISCUSSION ABOUT THE USE OF THE OSEEN TENSOR FOR DESCRIBING MOLECULAR MOTOR HYDRODYNAMIC COUPLING
In order to describe the hydrodynamic coupling between molecular motors pulling on a common cargo covered by a fluid-like membrane we have exploited the Oseen tensor (see Eqs.( 1),(2) in the main text).Indeed the Oseen tensor captures the far-field velocity profile in an unbound homogeneous fluid.However, in the present case molecular motors are running on a microtubule and the fluid flow they generate will be distorted by the presence of both the microtubule and the cargo motors are pulling on.Accordingly, in order to properly capture the presence of these boundaries, we should account for the appropriate image set.In some cases (e.g. for an infinite, solid planar wall) the presence of boundaries makes the leading contribution in the Green's function of the Stokes equation decay faster than the Stokeslet (captured by the Oseen tensor).However, even if the molecular motor has a diameter smaller than the cross section of a microtubule (even if not much smaller) the screening induced by the microtubule will be weaker than the one corresponding to a planar, infinite wall.Accordingly, we expect that the details of the solid cylinder will modify the induced hydrodynamic field only quantitatively (and in hydrodynamics such changes are usually gradual except very close to contact) with respect to the flow induced in an unbound system.Therefore, building the general framework on the Oseen tensor makes the derivations of the central equations and the analysis of their properties simpler and more transparent.Indeed, our approach allows us to identify the relevant coupling parameter, k.Even when the Oseen tensor is no longer relevant, the coupling parameter, k, remains the relevant parameter and only its functional dependence on the system geometry will be modified.Accordingly, our model is robust to the details of the hydrodynamic coupling and its results can be also meaningful in the context of coordinated molecular motor motion controlled by short range, unsteady hydrodynamic coupling (described in Ref. [7]).

DERIVATION OF EQS.(6) IN THE MAIN TEXT
The dynamics of molecular motors is an intrinsic stochastic process.Therefore in order to gain analytical insight it is useful to study the "typical" behavior, i.e. averaged over the realization of the noise that affects the active stepping of the molecular motors.Following the approach developed in Ref. [14,21] we write down the functional Smoluchowsky equation governing the time evolution of the probability distribution P [ρ, σ] of the dimensionless fields ρ(x) and σ(x): where v ρ (x) and v σ (x) are the local effective velocity field experienced by bound and weakly bound motors respectively and in the mean-field approach are defined as: where with ... ρ,σ we mean the average over all possible density distribution ρ(x), σ(x) each of which weighted with probability P and with ... x the average over all positions x, x ∈ [0, L] weighted with constant probability p(x) = dx/L.From Eq. 10 we can write down the evolution equation for the average values of the fields ρ(x) and σ(x): By performing a mean-field approximation, namely assuming and using the definition we can rewrite the previous equations as:

DERIVATION OF EQS.(8) IN THE MAIN TEXT
For fast bulk kinetics k on,off → ∞ Eq. ( 20) reads: substituting the last expression in the first of Eqs.( 6) of the main text we obtain Eq.( 8) of the main text.

DERIVATION OF EQS.(9) IN THE MAIN TEXT
The bifurcation portrait shown in Fig. 1 of the main text has been obtained by numerically solving Eq.( 8) (via a Lax-Wendroff method [22]) with periodic boundary conditions applied at the end of the extrema of the ratchet potential.In the following we derive an approximated scheme that allow us to discuss, from an analytical perspective, the bifurcations obtained numerically.We start expressing ρ by its Fourier series: ρ(x) = ρ 0 + Σ n ρ n cos(2πnx/L) + ρn sin(2πnx/L) (22) provides insight into the spontaneous symmetry breaking shown in Fig. 1 of the main text.By substituting Eq. ( 22) into Eq.(8) of the main text and using Eq.( 7) of the main text we obtain the following cascade of equations: The set of Eqs.(23a)-(23e) is an infinite cascade that, in order to be treated analytically should be truncated.As a closure of the cascade we assume that the Fourier modes higher than the first ones (namely ρ 1 and ρ1 ) are decaying fast enough so that their amplitudes can be regarded as constants in the time scale in which ρ 1 and ρ1 are evolving.Accordingly, we regard ρ 2 and ρ2 as parameters whose values are obtained from the numerical solution of Eq.( 8) and we get a closed system for ρ 1 , ρ1 .Such a closure leads to good a match with the full numerical solutions.In particular, such a system has a zero-velocity solution: and, for ρ 1 = 0 a, possibly, moving solution characterized by where ρ 1 is obtained by solving: The necessary condition in order to observe a spontaneous symmetry breaking is that Eq. ( 27) must have three real solutions.For the case ρ2 = 0 (that coincides with the outcome of our numerical solutions) we can write down a necessary condition for having three real solutions of Eq. ( 27), namely: For rigidly coupled motors the first term in Eq. ( 28) vanishes and Eq. ( 28) is a necessary and sufficient condition [14] for the onset of a moving solution.In this case, by rearranging eq.28, the necessary and sufficient condition for symmetry breaking reads: According to Eq. ( 29) by increasing the dimensionless forcing λ we can diminish the minimum value of the coupling above which symmetry breaking occurs.For soft-hydrodynamic coupled motors the necessary, but not sufficient, condition for symmetry breaking is For small values of ρ 2 , as obtained from the numerical solutions of Eq.( 8), we can assume ρ 2 ≃ 0 and, using Eq.(23a), the last expression simplifies to i.e.Eq.( 9) of the main text.Since k ≥ 0, Eq. (31) identifies an interval λǫ[0; λ max ] for which symmetry breaking occurs, with λ max = 1 π δω ∆ω .According to Eq. (31), k diverges for both λ = 0, λ = λ max leading to the existence of a minimum value k = k c at λ = λ c = 1 2 λ max .Therefore Eq. (31) provides a prediction of the critical value, λ c above which we can observe spontaneous symmetry breaking.Substituting the values of δω/∆ω used in the numerical solutions of Eq.( 8) that is, the full set of Eqs.(23), into Eq.(31) we obtain λ c ≃ 8 • 10 −2 that is in good agreement with the value obtained from the full solution of Eq.( 8) shown in Fig. 1.Such an agreement between the numerical solutions and the approximated analytical stability analysis underlines the relevance of the slowest modes in the control of the instability.We remark that the linear stability analysis performed here is general and it is valid for all models that share the same functional form once linearized about the non-motile state.Therefore, since such a functional will not be affected by accounting for more detailed models for the dynamics of the molecular motors or by solving the Stokes equation beyond the Oseen approximation, the agreement between the linear stability analysis and the numerical solution shows that the phenomenology captured by our model is robust and will persist also for more detailed models that go beyond the mean-field approximations approach.

BISTABLE REGIME
As explained in the main text, Fig. 1 shows that velocity oscillates with frequency ω about a non vanishing value whose sign switches with frequency Ω.

FIG. 3 :
FIG.3:(a): minimum system size, Λ, calculated from Eq. (4) for symmetry breaking (solid lines) and for bistable onset (dashed lines) as a function of l = η2D,mem/η3D,cyt for motor pulling on membranes.The membrane-embedded tracer size, R, is of the order of the membrane thickness ∼ 4nm leading to R ∼ 1/2L.Thicker lines stands for larger values of k: k = 50, 70, 90 (solid lines) and k = 95, 110 (dashed lines) respectively.(b): values of the coupling parameters, k1, k2 at which the two bifurcations occur as a function of the motor properties encoded in λ = f0/(ωLγ).

FIG. 4 :
FIG.4: Time dependence of the velocity of bound motors in the bistable regime.