Formation of helical membrane tubes around microtubules by single-headed kinesin KIF1A

The kinesin-3 motor KIF1A is in charge of vesicular transport in neuronal axons. Its single-headed form is known to be very inefficient due to the presence of a diffusive state in the mechanochemical cycle. However, recent theoretical studies have suggested that these motors could largely enhance force generation by working in teams. Here we test this prediction by challenging single-headed KIF1A to extract membrane tubes from giant vesicles along microtubule filaments in a minimal in vitro system. Remarkably, not only KIF1A motors are able to extract tubes but they feature a novel phenomenon: tubes are wound around microtubules forming tubular helices. This finding reveals an unforeseen combination of cooperative force generation and self-organized manoeuvreing capability, suggesting that the diffusive state may be a key ingredient for collective motor performance under demanding traffic conditions. Hence, we conclude that KIF1A is a genuinely cooperative motor, possibly explaining its specificity to axonal trafficking.

In region A, motors are either strongly bound (state 1, black circles) or weakly bound (state 2, grey circles) to the MT. They are excited and decay with average rates ω , ω respectively. In region B, motors are detached from the MT (state 0, white circles) and diffuse freely outside the curved region of size r in the tip. Overlapping is allowed to account for the twodimensional diffusive motion of motors on the membrane. Motors detach from state 2 and attach from state 0 with mean rates ω 2 , ω a respectively. Finally, in region C, detached motors feel a soft repulsive potential V which prevents them to enter the curved region. b) Tube extracted from a GUV with surface density of motors ρ ∞ . An influx J + of motors enters the tube region at x (1) and an outflux J − enters back to the GUV at x (0) .  Bravais lattice representing the MT surface with directions r 1 , r 2 , forming an angle θ. The lattice unit is centered in each node. b) Lattice unit with lengths l 1 and l 2 and ratchet asymmetries a 1 , a 2 . c) Schematic description of a two-state noise-driven ratchet for the i-th direction and its equivalent description on a lattice. The transition p 12 is not depicted in this case. Motors can diffuse with rates u i , v i and advance with rates p i along the i-th direction. d) Two-dimensional diffusion on the lattice in the weakly bound state. e) Possible transitions between the strongly bound state (black circles) and the weakly bound state (grey circles) and its corresponding transition rates.

Supplementary Note 1: Geometry and energetics of a tubular helix
We consider a tubular helix of radius r (tube radius) winding on an imaginary cylinder of radius R (MT radius plus some extra space occupied by the motors, depicted as a dashed circle in Supplementary  Fig. 1a, left), such that r + R is the distance between the center of the MT and the center of the tube. The helix moves along the z-axis with an angular pitch p which is the distance the helix advances per radian along the z-axis. Each point on the surface of the helix is determined by the vector X. We can parametrize every point on the helix surface with two variables {s, φ}, where s is the arclength following the path of the helix through its center and φ is an angle which determines at a given arclength s the position of a point in the circle of radius r. If we call a the angle the helix moves per arclength unit, the distance the helix advances along the z-axis after moving a distance s along the helix will be pθ, where θ ≡ as. If we move an arclength unit, we are moving a distance ap in the z-axis and a distance aR 0 in the axis perpendicular to the z-axis, where R 0 ≡ R + r. In this way we find 1 = a 2 (R 2 0 + p 2 ). On the other hand, the angle ζ the center of the helix forms respect to the z-axis fulfills (see Supplementary Fig. 1b): The position of the center of the left-handed helix is r(s) = (R 0 sin θ, R 0 cos θ, pθ). We can parametrize the surface of the tubular helix X(s, φ) as: with φ ∈ [0, 2π), using the orthonormal triad {t, n, b} such that t(s) = ∂sr |∂sr| , n(s) = ∂st |∂st| and b(s) = t∧n |t∧n| (see Supplementary Fig. 1a, right). The mean curvature of the surface H(φ) can be determined through the first and second fundamental forms g µν (s, φ) = ∂ µ X · ∂ ν X and Π µν (s, φ) = ∂ µν X · n, µ, ν = s, φ; which in a matrix form read: The mean curvature can be calculated through the expression H(φ) = EN −2F M +GL 2 det(g) [1,2] which leads to: where K(φ) = − C cos φ r(1−rC cos φ) and C = R0 R 2 0 +p 2 is the curvature of the helical spine curve of the tube. We notice that for p → ∞ the curvature tends to the one of a cylinder of radius r i.e. 1/2r. The total surface of a tubular helix of lenght L and radius r is 2πrL, the same as for a cylinder. Consequently, the surface energy is 2πrLγ, and it does not depend on the pitch. Let us define the small quantity ≡ rC. This quantity is bounded below 1 in the experiments. Considering the first correction up to fourth order in , the free energy of the system can be approximated as: In the limit p → ∞ (or → 0), we recover the free energy of a cylindrical tube with extraction force F = 2π √ 2κγ and diameter 2r = κ/γ. Therefore, we conclude that the last two expressions for the cylindrical case are accurate to second order in and yield good approximations to the actual values for the helical tubes.

Superhelical effect on the tube pitch
Let us consider a helical tube of pitch P growing around a MT which has superhelicity of pitch P sh . The angle that the tube will advance per arclength unit will be: where p is the angular pitch, p is the angular pitch for a 13pf MT and ξ is the relative increase in angle which reads: where cot ζ sh = p sh /R 0 .

Estimation of the off-axis force exerted by the motors
Next we estimate the off-axis force exerted by the motors by using free energy arguments. From Eq. S5 we find that a lower bound of the excess free energy ∆F associated to the winding of the tube reads: This excess free energy is due to the work KIF1A motors perform in the off-axis direction W off along the tube. This work can be estimated to be: where F off is the total off-axis force exerted by the motors and N w is the winding number. Actually, the above expression underestimates the work by taking the radius of the displacement as that of the microtubule, and not that of the point at the membrane where the force is exerted. The winding number can be expressed in terms of the angle the helix moves per arclength unit a as N w = La/2π. Equating the last two expressions we obtain the total off-axis force: In our experiments 0.04 − 0.5. Which gives a lower bound of the total off-axis force in the range F off 0.04 − 2 pN.

Supplementary Note 2: In silico model for longitudinal tube pulling
We illustrate the problem of tube-pulling for a tube of radius r and extraction force F , considering N motors in the vesicle reservoir and assuming that motors use a single pf 1 extending the model presented in Ref. [3]. We distinguish two main regions in the system: the tube region and the vesicle region ( Supplementary Fig. 2b).

Tube region
A subset of motors N t ⊂ N are found in the tube region at time t. A motor i from this subset is found in state k i (t) at time t, where k i is a discrete stochastic variable. Motors can be detached from the MT (k i = 0), strongly bound to the MT (k i = 1) or weakly bound to the MT (k i = 2). In all cases, motors are also bound to the tube considered as a soft cargo. The dynamics of each motor is different depending on the region where it is found (A, B or C, see Supplementary Fig. 2a). Next we describe the dynamics in each region: Region A corresponds to the region in between the tube and the MT, were motors can be either weakly or strongly bound to the MT (i.e. k i = 0). In this region the dynamics reads [3]: where λ is the friction coefficient, such that the diffusion coefficient follows the Einstein relation is the potential motors feel depending on the state k i . In the strongly bound state (k i = 1), motors feel a periodic saw-tooth potential of asymmetry a and periodicity l. The results are not significantly affected by the particular form of the ratchet potential as shown in previous studies where similar results are obtained in a lattice model [3,4]. On the other hand, in the weakly bound state (k i = 2) they feel a constant potential U (x i , 2) = U 2 . W accounts for the motor-motor interaction potential and is taken as a truncated Lennard-Jones potential of hard-core size σ and energy minimum , where the former corresponds to the motor size and the latter is taken large enough to ensure that the interaction is effectively hardcore for motor-motor distances smaller than σ. The i-th motor interacts only with a subset of motors S(t) at time t which are bound to the MT, i.e.
. Finally x max corresponds to the position of the foremost motor in the system.
In region B, motors are detached from the MT (k i = 0) and undergo free diffusion on the membrane tube outside the curved region in the tip of size r. Defining ξ i = x i − x tip as the relative distance between the position of the i-th motor and the tip position x tip , the condition for a detached motor to be in region B is |ξ i | > r. The dynamics reads: where λ t , D t are the friction and diffusion coefficients on the membrane tube respectively. We notice there is no interaction potential for the motors in this region since we allow overlapping to account for the two-dimensional diffusive motion of motors on the tube. In this way, motors are no longer ordered respect to their label i. Region C corresponds to the curved region of the tip where |ξ i | < r. In this case we neglect noise and the dynamics simply read: (S13) 1 The case of Np pf is done by simply scaling the extraction force F and the surface density of motors ρ∞ by Np.
with a repulsive potential V in the form of a truncated Morse potential which is valid for |ξ i | < r and is zero otherwise. This is an ad hoc choice to simply prevent motors to enter in region C. We adjust the Morse parameters to ensure a soft short-ranged repulsive potential preventing the detached motors to enter the curved region. We associate the dynamics of the tube tip with the dynamics of the foremost motor i.e.ẋ tip =ẋ max . If the foremost motor detaches, the tube retracts with a retraction velocityẋ tip = −F/λ m until a new bound motor is found. λ m is an effective friction parameter which is inferred from the retraction velocity of motors v r 100 µm s −1 observed experimentally for F 20 pN [5]. Finally we estimate the number of bound motors in the tip over time by counting the number of consecutive bound motors pulling on the tube, where we define two motors as consecutive if they are at a distance less than δ D = 4D/ω. Also, we estimate the density of bound motors from the set of positions of the bound motors at each time step of the simulation using a smoothing technique. We define the density of motors at each point as the number of bound motors in a characteristic bandwidth divided by its length. The bandwidth size is taken 100 nm, which is similar to the experimental spatial resolution.
Next we describe the state dynamics. All state transitions in the system are stochastic with dwell times which are exponentially distributed. Motors are excited from state 1 to 2 from localized regions of size δ near the minima of the ratchet potential with average rate ω . On the other hand, decays from state 2 to 1 are delocalized with average rate ω. Attachment events occur with average rate ω a . However, they are not always possible due to excluded volume interactions in region A. Thus, we say that a motor i will only attach if it can find a free site i.e. j∈S W (x i − x j ) = 0. We include exponentially dependent detachment kinetics on the force. In vitro experiments using single-headed kinesin have shown that the detachment rate at zero load is much larger in the weakly bound state (∼ 1 s −1 ) than in the strongly bound state (∼ 0.01 s −1 ) [6]. We choose to neglect detachment from state 1 for simplicity. The addition of detachment in state 1 leads to similar dynamics in the system. The average detachment rate from state 2 of the i-th motor at time t, will depend on the passive forces the i-th motor feels over time: where ω 0 2 is the detachment rate at zero load from state 2,F i (t) is local time-averaging of the noisy signal of the passive forces F i (t) = − j∈S W (x i − x j ) − F δ xi,xmax and d is a characteristic length which is typically 2 to 4 nm for kinesin [6,7,8]. We used a simple smoothing technique by choosing F i (t) as the time average of F i in the region [t − τ, t], where τ is the window size. This is taken as τ = 100 ms, which is big enough to average the passive forces a motor feels during a hydrolysis cycle (∼ 10 ms), and smaller than the timescale of the tube motion (∼ 1 s). Variations of τ around this value did not affect significantly the dynamics of the system.

Vesicle region
The vesicle region is described as a motor reservoir with surface density of motors ρ ∞ . Motors diffuse on the vesicle and eventually they enter the tube region through the boundary x = x (0) ( Supplementary  Fig. 2b, red line). Hence, in the boundary we have an influx of motors J + (t) which are bound to the MT and to the tube. Experimental evidences show that it is feasible to neglect the influx of motors only bound to the tube [9]. On the other hand, there is also a flux of motors leaving the tube J − (t) by diffusion. Since KIF1A is able to make large backward excursions in the weakly bound state, for practical reasons it is important to ensure motors will not fluctuate near the boundary x = x (0) . Hence, we let motors appear at x (1) = x (0) +δ D , (Supplementary Fig. 2b, orange line) and let x = x (0) act as an absorbing boundary condition for the motors that leave the tube region. The number of motors in this region N t (t) will depend on time through the flux balance: Experimental evidences indicate that Eq. (S15) reaches a quasi-steady state [9]. Far from the tip, in the quasi-steady state the density of motors bound to the tube in the mean-field limit reads [5]: where in our case ω 0 d = ω 0 2 /(1 + β) is the motor detachment rate at zero load and β = ω/ω . Therefore, at the boundary x (1) , the average influx of motors will be J + = ρ b V 0 . In our simulations, a new motor will be introduced in the system stochastically every certain time taken from an exponential distribution with mean dwell time J + −1 . Finally, motors crossing the boundary x (0) will be incorporated in the vesicle reservoir.

Discussion on realistic parameters for the in silico model
The most determinant characteristic of single-headed KIF1A is that its movement relies on diffusion.
In vitro experiments have reported diffusion coefficients in the range of 20 to 40 nm 2 ms −1 which involve motor excursions much larger than the ratchet periodicity of 8 nm [10] . We will consider 20 nm 2 ms −1 as a reasonable value. The asymmetry a of the ratchet is an adjustable parameter for the model which is difficult to grasp from experiments. The asymmetry reduces the overall velocity of the system and it can lead to non-trivial effects specially in the limit of weak noise [11]. For our purposes, we adjust this parameter to 20 % of the periodicity length. Finally, the motor size σ is carefully chosen to avoid possible commensurability effects [11,12]. The characteristic rates ω and ω are found in the literature within the range of hundreds of Hz. Whereas ω is a parameter coming from the affinity between the motor domain and the MT, ω depends on ATP concentration in the solvent. Experimental data suggest that ω ≤ 250 s −1 and ω 250 s −1 [10,13]. We choose β such that the velocity of a single KIF1A at zero load is similar to the experimental gliding velocities ∼ 80 nm s −1 (see Methods). The resulting value is β 7.5. The detachment rate of KIF1A has been found to be ω d ∼ 0.1 s −1 . Using our value of β we get ω 0 2 1 s −1 , which is in agreement with the results in Ref. [6].
The radius of the tube r and the threshold force to extract a tube F depend on the bending rigidity κ and the surface tension of the membrane γ through the expressions r = κ/(2γ) and F = 2π √ 2κγ [14]. κ is assumed to be roughly constant in experiments whereas γ can vary substantially. Although in principle the surface tension can be adjusted in vitro by changing the osmotic pressure difference in the system, the statistical dispersion of γ from vesicle to vesicle is large and makes it difficult to control this parameter. The typical range of γ implicitly obtained through our data analysis is 3 × 10 −4 − 10 −1 pN nm −1 . The density of motors in the vesicle can be obtained assuming that each lipid occupies a surface of approximately 0.4 nm 2 . In the experiments we used different fractions of biotinylated lipids, in the range 0.01-1 mol %. This range corresponds to 200 − 20000 µm −2 . The diffusion coefficient on the tube is much larger than the diffusion coefficient in the weakly bound state (D t D). Typically, D t 1 µm 2 s −1 [9]. Finally, the allowed range of attachment rates is reported to be 0.1 s −1 ≤ ω a ≤ 10 s −1 [13], and we will take an intermediate value.

Supplementary Note 3: Lattice model for a single KIF1A motor
Here we present a phenomenological approach to describe the motion of a single-headed KIF1A motor moving in a two-dimensional lattice under an external load by taking into account the presence of two asymmetry parameters. In order to describe the motion of a single KIF1A motor on a 2D MT lattice, we extend the lattice description in Ref. [3] by considering a 2D oblique Bravais lattice (Supplementary Fig. 3a) with directions r 1 and r 2 forming an angle θ. The vector r 1 describes the on-axis motion of kinesin reflecting the polarity of the MT. On the other hand, the vector r 2 shows the off-axis motion determined by the bias due to the interaction of the motor domain and the MT lattice. The mechanism of motion can be understood by considering a properly generalized two-state model in two dimensions.
We define the state of the motor as k = 1, 2 depending on whether the motor is strongly bound (k = 1) or weakly bound (k = 2) to the MT. In the strongly bound state, we consider the motor feels a superposition of two asymmetric ratchet potentials forming an angle θ with asymmetries a i and periodicities l i (Supplementary Fig. 3b, c) where i = 1, 2 indicate the directions in the lattice. Each node in the lattice indicates the landscape minima. A motor in the strongly bound state can be excited with rate r = ω to the weakly bound state in which it undergoes two-dimensional diffusion on the lattice with rates u i , v i , i = 1, 2 ( Supplementary Fig. 3d). These rates depend linearly on an external i is the diffusion rate and D is the one-dimensional diffusion coefficient in the weakly bound state. When the motor decays, it can fall in one of the four possible regions depicted in Supplementary Fig. 3b. If the motor falls in the red region, it performs a q transition binding strongly to the MT in the same node ( Supplementary Figs. 3c and 3e.). However, if the motor falls in one of the three green regions, it will move to a new node in the lattice. The motor will make a p 1 transition if it falls in the upper green region, a p 2 transition if it falls in the left green region, and a p 12 transition if it falls in the dark green region. The latter probabilities are equal to the decay rate ω times the probability of falling in a given region, which can be directly obtained calculating the areas in Supplementary Fig. 3b. We define a given transition rate from state {k, R} to {k , R } as Γ(k , R |k, R). The different transition rates read: Γ(1, R + r 1 + r 2 |2, R) = p 12 = ω 4 (1 − 2ā 1 )(1 − 2ā 2 ) Γ(1, R + r 1 |2, R) = p 1 = ω 4 (1 − 2ā 1 )(1 + 2ā 2 ) Γ(1, R + r 2 |2, R) = p 2 = ω 4 (1 + 2ā 1 )(1 − 2ā 2 ) Γ(1, R|2, R) = q = ω 4 (1 + 2ā 1 )(1 + 2ā 2 ) Γ(2, R|1, R) = r = ω whereā i ≡ a i /l i . The mean velocities in the r 1 and r 2 directions will be: v 1 = l 1 σ ss 2 (p 12 + p 1 + u 1 − v 1 ) v 2 = l 2 σ ss 2 (p 12 + p 2 + u 2 − v 2 ) (S18) where σ ss 2 ≡ 1/(1 + β) is the steady state probability of the motor to be found in state k = 2 and β ≡ ω/ω . Substituting the form of the different rates, we can write the last expressions as: where v d i ≡ l i d i and f s i = ω 2di (1 − 2ā i ) is the dimensionless stall force in the direction r i . We notice we recover the results for the one-dimensional case in each direction [3]. The last expression can be rewritten in terms of the velocity at zero load v i (0) and stall force F s i as: where v i (0) = (l i ω/2)(1 − 2ā i )/(1 + β) and F s i = (k B T ω/l i d i )(1 − 2ā i ). Using the KIF1A parameters from the Supplementary Table 1