On the force--velocity relationship of a bundle of rigid living filaments
Alessia Perilli, Carlo Pierleoni, Giovanni Ciccotti, Jean-Paul Ryckaert
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A ug On the force–velocity relationship of a bundle of rigid living filaments.
Alessia Perilli ∗ Department of Physics, Sapienza University of Rome, P.le Aldo Moro 5, I-00185 Rome, Italy;and, Department of Chemistry, ´Ecole Normale Superi´eure, rue Lhomond 24, 75005 Paris
Carlo Pierleoni † Department of Physical and Chemical Sciences, University of L’Aquila, Via Vetoio 10, 67100 L’Aquila, Italy;and, Maison de la Simulation, CEA–Saclay 91120 Gif sur Yvette, France
Giovanni Ciccotti ‡ Instituto per le Applicazioni del Calcolo “Mauro Picone” (IAC), CNR, Via dei Taurini 19, I-00185 Rome, Italy;and, Sapienza University of Rome, P.le Aldo Moro 5, I-00185 Rome,Italy; and, University College Dublin (UCD), Belfield, Dublin 4, Ireland
Jean-Paul Ryckaert § Department of Physics, Universit´e Libre de Brussels (ULB),Campus Plaine, CP 223, B-1050 Brussels, Belgium (Dated: February 22, 2018)In various cellular processes, biofilaments like F–actin and F–tubulin are able to exploitchemical energy associated to polymerization to perform mechanical work against an externalload. The force–velocity relationship quantitatively summarizes the nature of this process. By astochastic dynamical model, we give, together with the evolution of a staggered bundle of N f rigidliving filaments facing a loaded wall, the corresponding force–velocity relationship. We computesystematically the simplified evolution of the model in supercritical conditions ˆ ρ = U /W > ǫ = d W /D = 0, where d is the monomer size, D is the obstacle diffusion coefficient, U and W are the polymerization and depolymerization rates. Moreover, we see that the solution at ǫ = 0 is valid for a good range of small non–zero ǫ values. We consider two classical protocols:the bundle is opposed either to a constant load or to an optical trap set–up, characterized by aharmonic restoring force. The constant force case leads, for each F value, to a stationary velocity V stat ( F ; N f , ρ ) after a relaxation with characteristic time τ micro ( F ). When the bundle (initiallytaken as an assembly of filament seeds) is subjected to a harmonic restoring force (optical trapload), the bundle elongates and the load increases up to stalling (equilibrium) over a characteristictime τ OT . Extracted from this single experiment, the force–velocity V OT ( F ; N f , ρ ) curve is foundto coincide with V stat ( F ; N f , ρ ), except at low loads. We show that this result follows from theadiabatic separation between τ micro and τ OT , i.e. τ micro ≪ τ OT . Submitted to the Journal of Chemical Physics on August 22nd 2017 ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] I. INTRODUCTION
Cell motility in vivo is a large scale manifestation of the living character of the cytoskeleton bio–filaments network [1].In particular, F–actin filaments produce growing lamellipodium or filopodium structures where G–actin monomerspolymerize at the barbed end of filaments, directly in contact with the cytoplasmic membrane. The speed of themembrane deformation/displacement at the leading edge of the cell adjusts itself so that the force generated by thegrowing filaments compensates the resisting load coming from the membrane tension and to the crowded environmentaround the cell. For living filaments opposing a loaded mobile obstacle, the macroscopic force–velocity relationship, V ( F ), linking the obstacle velocity, V , only to the instantaneous applied load, F , quantitatively summarizes thecombined action of the elementary self–assembling processes. In such adiabatic conditions, implying a time scaleseparation between the self–assembling process and the response of the obstacle, the V ( F ) dependence could beprobed equivalently by different protocols like, to cite the two most frequently used, the constant force load (e.g.clamped force set–up), where one directly observes the steady state velocity, and the harmonic load (the samplegrows against an AFM cantilever or an optical trap), where the obstacle velocity can be followed as the load increasescontinuously up to stalling.Abiabatic conditions cannot be in general guaranteed and indeed, careful investigations on a branched actin networkgrowing against an AFM tip have shown that the recorded V ( F ) relationship can be function of the load history [2].The direct force–velocity relationship, V ( F ), is in any event widely used as a characteristics of network dynamics tocompare experimental measurements and modeling approaches for in–vitro [3] and in–vivo systems [4].To make progress on the rationalization of the conditions of validity of the widely used concept of force–velocityrelationship, V ( F ), we will restrict our considerations to a simple network where a bundle of parallel (proto)filaments(actin or tubulin) grows normally against a loaded obstacle. The general mechanism, linking work production and(de)polymerization kinetics of living bio–filaments, has been originally formulated theoretically by Hill for an incom-pressible bundle of N f parallel filaments pressing against a mobile obstacle [5]. Successively, when the filaments of thebundle are treated as independent and equivalent and when it is assumed that the depolymerization rate is unaffectedby the external load, the wall velocity V MF ( M F indicates the mean field character of this treatment) has beenwritten as [6, 7] V MF ( F ; ρ , N f ) = d (cid:20) U exp (cid:18) − F dN f k B T (cid:19) − W (cid:21) (1)where U = k on ρ and W = k off are the single filament bulk rate constants, related to bulk chemical rate constants k on and k off , for single monomer polymerization and depolymerization steps, ρ is the free monomer density, d is thesingle filament increment of contour length per incorporated monomer and F is the external force exerted on the wall.Supercritical conditions, where filament polymerization dominates over depolymerization, require ˆ ρ = ρ ρ c = U W > ρ c = k off k on is the critical value of the monomer density at which the bundle has no tendency to grow nor toshrink in absence of load.Eq.(1) predicts, for F = 0, a growth velocity of the free bundle V MF = d ( U − W ) > F s ,at which the velocity vanishes, is given by F Hs = N f k B Td ln ˆ ρ . (2)The notation F Hs reminds that this expression was originally established by Hill using thermodynamic arguments [8].Eq.(2) has been recently derived, in a special limit, by equilibrium Statistical Mechanics for a bundle of rigid filaments[9]. Indeed, it has been found that the statistical mechanics average of the wall position, taken over the equilibriumoptical trap ensemble, multiplied by κ T , converges exponentially fast for κ T → in–vitro the grafted bundle seed needed to follow its subsequent loadedgrowth. However, it is interesting to note the diversity in these few approaches. The growth of single grafted tubulinfilaments, which are bundles of 13 proto–filaments, was followed by imaging techniques [7, 10]. Regrouping ( F, V )data for different observation times and for different samples, a master force–velocity relation could be established. Inanother experiment using an acrosome bead complex of N f = 8 ÷
10 F–actin filaments held in an optical trap device,the growth of a bundle was followed in time against an harmonic load [11]. A rising signal finishing with a plateauwas observed but the final stationary force was surprisingly much lower than the expected stalling force, Eq.(2), itsvalue being close to the stalling force predicted for a single filament. The analysis in this experiment considers manyrelaxation curves, but in many cases data had to be eliminated due to interferences during the relaxation process withthe onset of escaping filaments. This happens because growing filaments undergo a large bending fluctuation whichallows them to start growing freely along the obstacle. The transient behavior, which can be converted into a V ( F )law by the derivation in time of the wall position, was not exploited. Finally, in a recent study [12], recording therate of radial distance between two colloidal particles separated by a growing grafted actin bundle, the force–velocityrelationship of actin bundles was established in constant load conditions.The outcome of the earliest experiment [7, 10], confirmed by the more recent experimental work [12], is that thevelocity, and hence the power of transduction of multi filament bundles, is much lower than predicted by Eq.(1).The discrepancy highlights the non–independence of elementary chemical steps at the tip of different filaments inthe bundle, with the effect of reducing the additivity of the action of each filament. The bundle model needs to bespecified and the dependence between chemical events and wall position for a given longitudinal seeds disposition hasto be quantitatively taken into account. This aspect is present in the multi–filament Brownian Ratchet (BR) models[7, 12–14] which generalize the single filament brownian ratchet model introduced by Peskin et al. [15], for which onefinds that the velocity vanishes for a load equal to Hill’s expression, Eq.(2) [7]. For these bundle models, the importantcharacteristics which distinguish the dynamical behavior of the bundle are the number of rigid living filaments, thelongitudinal disposition of the seeds of the filaments and the wall diffusion coefficient D which introduces a secondcharacteristic time τ D = d /D next to the chemical events time scale τ chem = W − . This fact suggests to introducethe parameter ǫ = τ D /τ chem to be able to discuss the condition of this second adiabatic separation (not to be confusedwith the one associated to the existence of V ( F )). For both experiments having probed the V ( F ) relationship, itwas found that data could be interpreted successfully with the model of a staggered bundle (= staggered longitudinalseed disposition [9]) of rigid filaments in very fast wall diffusion conditions ( ǫ = 0), a model we will denote as SRBR(Staggered Rigid Brownian Ratchet). On the contrary, for a similar model with an in registry bundle (unstaggeredlongitudinal seed disposition) [14], the predicted velocity was much too low with respect to the experimental data[7, 10, 12].In the stochastic dynamical models here considered, the force–velocity relationship depends parametrically, for agiven seed arrangement, on the number of filaments, the reduced free monomer concentration and the time scale ratio ǫ . In the case of constant load, the explicit form for the asymptotic force–velocity relationship, V stat ( F ; N f , ˆ ρ , ǫ ), forour models has been established by stochastic dynamics studies at finite ǫ [16] and at ǫ = 0 [7, 12]. In the latter casea simplified algorithm, exploiting the time scales separation, has been used for the staggered bundle case. Indeed, at ǫ = 0, the wall position distribution at given filaments configuration, is found to be time–independent and equal tothe equilibrium distribution of the wall position resulting from the 1D Brownian motion of a wall in the external loadfield, with the wall positions restricted to be greater than the position of the most advanced filament tip.Interestingly, we add that for the SRBR model ( ǫ = 0), successive theoretical developments [7, 12, 13] havegiven, with a very good approximation [12], two coupled closed expressions for the velocity V stat ( F ; N f , ˆ ρ ) and thedistribution of filament relative sizes (see Section III) g ( k ; F, N f , ˆ ρ ), for the stationary state.In this work, we consider the stochastic staggered bundle model of rigid filaments in supercritical conditions andperform a series of dynamical runs for different load conditions. We first look at the constant force case, treatingboth the stationary state itself and the asymptotic transient evolution to reach it. We next envisage the bundle, insimilar thermodynamic conditions, initially taken with very short filaments, subject to a harmonic load − κ T L , where L is the wall position and where κ T is trap strength (optical trap set–up). Mimicking the optical trap experiment[11], the bundle and the average wall position grow and reach stalling. We compare our computed longest relaxationtime with a theoretical approximate expression derived along the lines of the D´emoulin et al. theory. We derive andcompare the force–velocity relationship extracted from this optical trap relaxation with the one obtained in stationaryconditions. As expected, we found that the two coincide in adiabatic conditions, i.e. when the characteristic time ofthe optical trap relaxation is much larger then the characteristic time of the relaxation in the constant force case.Our algorithms follow the same lines of those used in previous studies. However, in our study we deal with anoptical trap load, while most studies (with an exception restricted to the ǫ = 0 case [17]) assume a constant load.Moreover, while algorithms for finite ǫ or ǫ = 0 are usually just assumed, we establish an explicit link showing howthe ǫ = 0 model is derived from the general finite ǫ case.In Section II we present the general Fokker–Planck model for a bundle of rigid filaments with an arbitrary seeddisposition, facing either a constant or a harmonic load and we derive the explicit wall algorithm (EWA) giving thesampling rules to generate stochastic trajectories for any finite ǫ case. We then use a perturbation expansion to derivethe ǫ = 0 model, still for the constant load or the harmonic load, and we derive the simplified implicit wall algorithm(IWA) giving the sampling rules in the ǫ = 0 case. Section III reports and discusses our results for constant forceand optical trap loads for the same bundle system generally using the ǫ = 0 approach, since the wall diffusion takesplace very quickly with respect to the mean time between (de)polymerization events. However, we also verify thatthe simplified algorithm is robust, since we find identical results in a reasonable range of ǫ non–zero values. SectionIV concludes with a summary of the main results and with some perspectives. II. MODEL AND IMPLEMENTATION
We consider a bundle of N f > x axis) to a fixed planar substratewall (along y and z directions). The filaments are modeled as discrete rigid linear chains with monomer size d andlength related to the number of attached monomers, i ≥
2, as L ci = ( i − d . Let h n be the location along x axisof the seed (first monomer) of the filament n close to the grafting plane ( − d/ < h n < d/ in–registry (or unstaggered ), where h n = 0, n = 1 , N f , and homogeneous (or staggered ), where seeds are regularly spaced as h n = (cid:20) n − . N f − . (cid:21) d n = 1 , N f . (3)A moving obstacle, a hard wall located at distance L from the parallel substrate wall, is loaded with a compressionalexternal force F bringing it into contact with the living filaments. We will consider two types of load, the constantforce F , and the optical trap setting, with F = − κ T L , where κ T is the trap stiffness and L the distance between thewalls.The bundle force for rigid filaments is impulsive. Its effect is taken into account by imposing a confining boundaryto the wall motion at the tip location of the longest filament.Filaments either grow by a single monomer polymerization step with bulk rate U , proportional to the free monomerdensity ρ , or shrink by a single monomer depolymerization step with bulk rate W . The ratio U /W = ˆ ρ is thefree monomer density divided by its critical value, i.e. the value at which the two bulk rates are equal. We will beinterested to supercritical conditions only (ˆ ρ > d to the wall, the polymerization rate becomes zero, while the depolymerization one isassumed to remain unchanged.The dynamics of the bundle of growing filaments against the loaded mobile wall presents two main time scales: thechemical one, τ chem = 1 /W ∼ /U , and the diffusive one, related to the diffusive motion of the wall and estimatedby τ D = d /D where D is the wall diffusion coefficient. The ratio of time scales, ǫ = τ D /τ chem , in typical in–vitro experiments is ǫ ≪
1, but might sometimes go close to 1 for a very large colloidal particle in a crowded environment.In the following, we establish a Fokker–Planck equation to describe the dynamics of an arbitrary bundle of inde-pendent rigid filaments subjected to a constant or harmonic load, for arbitrary value of ǫ . A. General Fokker–Plank equation for a bundle of rigid filaments against a constant or harmonic load
We describe the time evolution of N f filaments against a load in terms of the filament sizes and the wall position, { j , . . . , j N f , L } . The wall position must always lie beyond the tip of any filament – and so beyond the tip of themost advanced one, n ∗ with size j n ∗ . Defining X n ( j n ), the position of the tip of filament n and X ∗ that of the mostadvanced one, X n ( j n ) = ( j n − d + h n (4) X ∗ ≡ X ∗ ( j , . . . , j N f ) = max n =1 ,N f { X n ( j n ) } = X n ∗ ( j n ∗ ) , (5) L > X ∗ . (6)We assume that the joint probability distribution function P j ,...,j Nf ( L, t ) satisfies a Fokker–Planck equation in timemixing a continuous process in space for the wall position with a discrete process for filament sizes. For the modeldescribed above, we have ∂P j ,...,j Nf ( L, t ) ∂t + ∂∂L J j ,...,j Nf ( L, t ) = U N f X n =1 (1 − δ ,j n )Θ ( L − X n ( j n )) P j ,...,j n − ,...,j Nf ( L, t ) − N f X n =1 Θ ( L − X n ( j n + 1)) P j ,...,j n ,...,j Nf ( L, t ) + W N f X n =1 P j ,...,j n +1 ,...,j Nf ( L, t ) − N f X n =1 (1 − δ ,j n ) P j ,...,j n ,...,j Nf ( L, t ) (7)where Θ( x ) is the Heaviside step function and the probability current density is J j ,...,j Nf ( L, t ) = − D " ∂P j ,...,j Nf ( L, t ) ∂L − F ( L ) k B T P j ,...,j Nf ( L, t ) (8)In Eq.(8) the compressive force can be either a constant F < F ( L ) = − κ T L modeling the opticaltrap. The right–hand side of Eq.(7) represents the sink and source terms affecting the dynamics due to polymerizationand depolymerization events. Their explicit expression indicates that, in one step at fixed L , transitions are onlypossible between adjacent microscopic states, where ( N f −
1) filaments have the same size while the size of theremaining filament differs by ± L > X ∗ , and that the filament size cannotbe smaller than two.The general normalization condition for the distribution P j ,...,j Nf ( L, t ) is ∞ X j =2 · · · ∞ X j Nf =2 Z ∞ X ∗ dL P j ,...,j Nf ( L, t ) = 1 (9)while the boundary conditions on the probabilities are P j ,...,j Nf ( L, t ) (cid:12)(cid:12) L 1. We then substitute to the wall position L the discretevariable k = int (cid:20) Lδ (cid:21) ≡ int [ l ] . (12)In this way Eq.(7) will become a finite difference equation in k representing a discrete Markov chain in continuoustime d P dt = P Q (13)with P ( t ) = {P j ,...,j Nf ,k ( t ) } j n ∈ [2 , ∞ ) n =1 ,N f ,k ∈ [ int [ ( d + h Nf ) /δ ] , ∞ ) a vector field and Q the generator matrix of theMarkov chain. The elements of the matrix Q contain the (de)polymerization rates for the filaments, U j n ( L ) = U Θ ( L − X n ( j n + 1)) (14) W j n ( L ) = W (15)and the forward/backward jump rates for the wall; the expressions of these matrix elements are given in Appendix A.To circumvent the difficulty of solving analytically Eq.(13), one can produce a number of realizations of the discreteMarkov chain using any appropriate algorithm, in our case the Gillespie algorithm [19, 20]: given an initial conditionat time t , the state of the system is estimated in terms of the set of random variables { j , . . . , j N f , k } at time t producing statistically correct trajectories, from which the probability distribution function P j ,...,j Nf ,k ( t ) can beinferred by histograms. Starting from the initial state, the system is allowed to evolve by random steps involvingonly one reaction per time: one filament depolymerization or polymerization, or the wall forward or backward jump.Denoting by i the current state of the system, the reachable states i m are those differing from i for only onevariable by ± 1, namely { j , . . . , j n ± , . . . , j N f , k } or { j , . . . , j n , . . . , j N f , k ± } . It is straightforward to see that thenumber of these possible final states is 2 N f + 2. The transitions i → i m , m ∈ [1 , N f + 2], are described in Eq.(13)by the generator matrix elements Q i m i , the rates of going from i to i m . The corresponding diagonal element is Q i i = − P i m = i Q i m i [21]. The evolution of the system is determined by two random variables: the time to thenext reaction, τ , and the final state i m , or equivalently the index of the reaction, m ∈ [1 , N f + 2]. From generalMarkov chain theory, τ is known to be an exponentially distributed random variable: given the current state i , theparameter of the exponential distribution is given by − Q i i . Instead, the probability for the jump m linking states i and i m to take place is given by the ratio between Q i i m and | Q i i | [21]. The main loop of the algorithm followsthis scheme:0. The initial state i is specified in terms of the state vector { j , . . . , j N f , k } . We take for the initial value of k asmall fixed value and, for the filament, compatible initial sizes;1. The matrix elements Q i i m are calculated for any state i m reachable from i ;2. The time to the next move is determined using the so–called direct method, which follows from the standardinversion method of the Monte Carlo theory [22]: a random number r ∈ [0 , 1] is generated from the uniformdistribution and the time τ is taken as τ = 1 | Q i i | ln 1 r ; (16)3. The index of the next move is determined using the same method: a second random number r ∈ [0 , 1] isgenerated and the index m is taken as the smallest integer satisfying m − X n =1 Q i i n | Q i i | < r m X n =1 Q i i n | Q i i | ; (17)4. The sampled move is taken by updating the state vector i → i m and the time is incremented by τ ;5. Go back to 1, until a maximum time t max is reached;6. End the simulation.The state vector { j , . . . , j N f , k } is stored for the calculation of histograms and averages.This algorithm [18, 19], solving the Fokker–Planck Eqs.(7, 8), works for any seed disposition (staggered and un-staggered), for any finite value of the dimensionless parameter ǫ ≡ d W D and both for the two cases of constant forceand optical trap load. We will call it the explicit wall algorithm (EWA).In the next subsection we treat the specific, important, reference case of loaded bundles of rigid filaments in thelimit ǫ → 0. In this limit the wall re-equilibrates instantaneously after any change of the position of the most advancedtip of the bundle. The interest of this limit is justified since in in − vitro experiment with actin bundles / colloidalparticles (e.g. the optical trap experiment [11]) the typical value of the ratio of time scales is ǫ ≪ 1. We will see thatthe dynamics of the bundle simplifies for two reasons: the elimination of the fast motion of the wall permits to goto longer times and the dimensionality of the problem is reduced. The new algorithm, called implicit wall algorithm(IWA), becomes then decidedly more efficient. B. Treatment of the Fokker–Planck equation in the ǫ = 0 limit [23] Given the separation of time scales between the chemical events and the wall diffusion, it is convenient to rewriteEq.(7) in terms of dimensionless variables in order to put in evidence the ratio ǫ = τ D τ chem = W d D . Defining ˜ t = W t , x = Ld and f = F dk B T , multiplying Eq.(7) by d D and redefining the probability distribution functions, we get: ǫ ∂ e P j ,...,j Nf ( x, ˜ t ) ∂ ˜ t + ∂∂x e J j ,...,j Nf ( x, ˜ t ) = ǫ ( ˆ ρ N f X n =1 (1 − δ ,j n )Θ ( x − X n ( j n ) /d ) e P j ,...,j n − ,...,j Nf ( x, ˜ t ) − N f X n =1 Θ ( x − X n ( j n + 1) /d ) e P j ,...,j n ,...,j Nf ( x, ˜ t ) + N f X n =1 e P j ,...,j n +1 ,...,j Nf ( x, ˜ t ) − N f X n =1 (1 − δ ,j n ) e P j ,...,j n ,...,j Nf ( x, ˜ t ) ) (18)with e J j ,...,j Nf ( x, ˜ t ) = − ∂∂x e P j ,...,j Nf ( x, ˜ t ) − f ( x ) e P j ,...,j Nf ( x, ˜ t ) (19)the probability current density in the reduced units. In the ǫ → ǫ zero–th order approximation: ∂ e P (0) j ,...,j Nf ( x, ˜ t ) ∂x − ∂∂x h − f ( x ) e P (0) j ,...,j Nf ( x, ˜ t ) i = 0 . (20)By integrating in dx form X ∗ to ∞ and using the boundary conditions for the probability, one gets ∂ e P (0) ∂x = − f ( x ) e P (0) so that the general solution is: e P (0) j ,...,j Nf ( x, ˜ t ) = a ( j , . . . , j N f , ˜ t ) exp (cid:18) − Z ∞ x dxf ( x ) (cid:19) (21)On the other side, it is always possible to write the joint probability as the product of the marginal distribution forthe subset { j , . . . , j N f } times the conditional probability distribution for x : e P (0) j ,...,j Nf ( x, ˜ t ) = e P ( j , . . . , j N f , ˜ t ) e P ( x | j , . . . , j N f , ˜ t ) (22)Therefore, given the general solution (21), we can write it as e P (0) j ,...,j Nf ( x, ˜ t ) = e P ( j , . . . , j N f , ˜ t ) e P EQ ( x | j , . . . , j N f ) (23)as the x dependence, Eq.(21), is explicit and time–independent. The wall distribution e P EQ ( x | j , . . . , j N f ) is anexplicit, time independent, normalized, distribution for the wall position conditional to the set of filaments sizes.Explicit expression for the two normalized cases of constant load and optical trap are: e P EQ ( x | j , . . . , j N f ) = f exp( − f x )exp( − f X ∗ /d ) constant load r κ T π exp (cid:0) − ˜ κ T x (cid:1) erfc (cid:16)q ˜ κ T X ∗ /d (cid:17) optical trap (24)with ˜ κ T = κ T d k B T . From Eq.(24) we get the average wall position conditional to the bundle sizes { j , . . . , j N f } as: E ( x | j , . . . , j N f ) = Z ∞ X ∗ x e P EQ ( x | j , . . . , j N f ) dx = X ∗ d + 1 f constant load r κ T π exp h − ˜ κ T ( X ∗ /d ) i erfc (cid:16)q ˜ κ T X ∗ /d (cid:17) optical trap (25)Note that the full distribution at ǫ = 0, given by Eq.(23), is still a time–dependent function, since filamentsizes change by single monomer polymerization/depolymerization events; the infinite separation of the time scales( ǫ = 0) implies that after any chemical event, the wall immediately re–equilibrates according to the time–independentdistribution, Eq.(24), given the new set of filament sizes. To get the full distribution, we write e P j ,...,j Nf ( x, ˜ t ) as anasymptotic expansion in terms of the small parameter ǫ : e P j ,...,j Nf ( x, ˜ t ) = e P (0) j ,...,j Nf ( x, ˜ t ) + ǫ e P (1) j ,...,j Nf ( x, ˜ t ) + . . . (26)where e P (0) j ,...,j Nf ( x, ˜ t ) is given by Eq.(23). If we substitute this expansion, truncated to the first order, into Eq.(18),to the order ǫ we find the following equation: ∂ e P (0) j ,...,j Nf ( x, ˜ t ) ∂ ˜ t + ∂∂x e J (1) j ,...,j Nf ( x, ˜ t ) =ˆ ρ N f X n =1 (1 − δ ,j n )Θ ( x − X n ( j n ) /d ) e P (0) j ,...,j n − ,...,j Nf ( x, ˜ t ) − N f X n =1 Θ ( x − X n ( j n + 1) /d ) e P (0) j ,...,j n ,...,j Nf ( x, ˜ t ) + N f X n =1 e P (0) j ,...,j n +1 ,...,j Nf ( x, ˜ t ) − N f X n =1 (1 − δ ,j n ) e P (0) j ,...,j n ,...,j Nf ( x, ˜ t ) (27)Integrating both sides of this equation from x = X ∗ /d to ∞ , applying the boundary conditions Eq.(11) on e J (1) andthe normalization of e P EQ and using Eq.(23), we get: ∂ e P ( j , . . . , j N f , ˜ t ) ∂ ˜ t =ˆ ρ N f X n =1 (1 − δ ,j n ) Z ∞ X ∗ /d dx Θ ( x − X n ( j n ) /d ) e P EQ ( x | j , . . . , j n − , . . . , j N f ) e P ( j , . . . , j n − , . . . , j N f , ˜ t ) − ˆ ρ N f X n =1 Z ∞ X ∗ /d dx Θ ( x − X n ( j n + 1) /d ) e P EQ ( x | j , . . . , j n , . . . , j N f ) e P ( j , . . . , j n , . . . , j N f , ˜ t )+ N f X n =1 e P ( j , . . . , j n + 1 , . . . , j N f , ˜ t ) − N f X n =1 (1 − δ ,j n ) e P ( j , . . . , j n , . . . , j N f , ˜ t ) (28)where we couldn’t use the normalization condition for the terms where we have left the integration explicitly written.This equation describes a discrete process for the filament sizes in continuous time, which can be rewritten in avectorial form, similar to Eq.(13): dP dt = P Q (0) (29)with Q (0) generator matrix of the process, whose elements are given in Appendix B.The numerical solution of the Markov chain equation described by Eq.(29) follows exactly the same scheme describedabove for the general Fokker–Planck equation for ǫ > h L i t . Similar model and procedures have been used:i. for constant load option and in-registry [14] or staggered [7, 12, 13] bundles; ii. for optical trap only for staggeredbundles [17]. III. SIMULATIONS AND RESULTSA. Units, parameters and stochastic runs In our simulations, length, time and energy units are taken as d , W − , and k B T respectively. All quantities willbe mentioned in reduced units based on the above three fundamental units. For actin d = 2 . nm ; experimentalinformation for W gives W = 1 . s − ; and, at room temperature k B T = 4 . × − J . We choose to performour studies on a bundle of N f = 32 rigid filaments with a staggered disposition of seeds at a reduced densityˆ ρ = U W = 2 . 5. With reference to a wall constituted by a bead of micron size in water opposing the actin bundles[11, 12], experimental information gives for the adimensional parameter introduced in the previous section, the value ǫ = 5 . × − . Given the small value of ǫ , we performed the major part of our simulations in the ǫ = 0 limit withthe IWA algorithm. However, we have considered interesting to compare the results of the IWA algorithm with thoseof the EWA corresponding to a finite but small value of ǫ . With the very small experimental value of ǫ , EWA wouldbe highly inefficient, since the computer time would be essentially spent to study the wall diffusion next to a bundlewith quasi-fixed filament sizes. Since for the load-velocity relationship we need to sample both wall and filamentsizes, we decided to adopt a value of epsilon thousand times bigger, ǫ = 5 × − . This value, in fact, still permits togive a sufficient representation of the wall dynamics. Our EWA approach requires to discretize the space variable L with elementary steps δ = d/M . For M , we have adopted M = 100. To compute the solution of our Fokker–Planckequation, both for ǫ = 0 (IWA) or ǫ > L (i.e. k = int( L /δ )) and to sample the initial filament sizes for each trajectory of thestochastic dynamics according to the filament size equilibrium probability P eq ( j , j , ....j n , ....j N f ; L ), conditional tothe chosen wall location [9]. For initiating IWA runs, the initial filament sizes must be arbitrarily chosen and theinitial wall location then follows from its conditional distribution. B. Observables of interest 1) Wall positionThe wall position L is the quantity directly followed in time in real experiments and corresponds to the expectedvalue of the random variable ˆ L over the solution of the FP equation, h ˆ L i t . The calculation of this quantity is directin the EWA case, while it has to be determined in the IWA case through the instantaneous size distribution of j n , n = 1 , 32, implying ˆ L values ahead of the tip of the most advanced filament at X ∗ , given by Eq.(5), using Eq.(25).2) Relative size (in number of monomers) of filaments with respect to the leading one.In terms of the tip positions X n and X ∗ defined by Eqs.(4, 5), let us define the relative subset index m = 1 , N f − m ( n ) = mod (cid:18) X ∗ − X n d/N f , N f (cid:19) = mod (( j n ∗ − j n ) N f + n ∗ − n, N f ) n = 1 , N f ; n = n ∗ (30)This index represents in successive order the filament of order n nearest neighbor of n ∗ , second neighbor of n ∗ , etc.Therefore it gives an intrinsic order to the vector representing the relative size of each filament. Note that the dividendin the function mod in the given condition is always positive. Then we can define, for each filament n , the quantity k m = int (cid:20) X ∗ − X n ( m ) d (cid:21) = int (cid:20) j n ∗ − j n ( m ) + n ∗ − n ( m ) N f (cid:21) m = 1 , N f − d the relative distance from the mostadvanced tip of the first, second, etc. neighboring index.This vector of relative sizes is interesting because its time-dependent probability distribution reaches a stationaryvalue in the case of the wall subjected to a constant load.3) Density of relative size of N f − g ( k ) = N f − P N f − m =1 δ k,k m (32)At time t , the microscopic distribution will be g ( k, t ) = h ˆ g ( k ) i t . Specifically, we will characterize the internalstructure of the bundle either by g (0 , t ), the average probability densitythat the tip lies at a distance smaller than d from the tip of the most advanced filament, or the average relative size h k i t = P ∞ k =0 kg ( k, t ). We will denote by g ( k )and k av the time–asymptotic values of these quantities for constant load force dynamics. C. Constant force load We have computed the relaxation towards the stationary state for a homogeneous bundle of N f = 32 rigid livingfilaments at ˆ ρ = 2 . F . We have chosen various values of F in the range 0 . 25 with F s the stalling force, Eq.(2). For each load value, we have generally used the IWA algorithm toproduce 10 independent trajectories, starting at time 0 with all filament sizes set to the same value ( j n (0) = 500 , n =1 , j n = 2. That could happen0 F (k B T/d) V ( d W ) FIG. 1: Force–velocity relationship for a homogeneous bundle of N f = 32 rigid filaments at ˆ ρ = 2 . ǫ = 0). The V stat stationary velocity data points (red filled circles) are obtained as the asymptotic slope of h ˆ L i t for constant force runs at eachshown load value. Error bars are less than symbol sizes. The dashed green line is the D´emoulin et al. theoretical estimate of V stat ( F ) based on Eqs.(C1, C2). The blue continuous curve is the force–velocity relationship obtained by the optical traprelaxation at κ T = 0 . F = 29 . when F > F s with negative average velocities. In two cases, starting with L = 5 d , we have used the EWA algorithm,averaging over 10 independent trajectories. To determine the microscopic relaxation time of the bundle, we havefitted the asymptotic time evolution of the average wall position as h L i t = C + V stat t + C ′ exp ( − t/τ micro ). To getthe diffusion coefficient of the bundle Γ, we have also fitted the asymptotic behavior of the mean square elongation σ ( t ) = h ˆ L i t − h ˆ L i t ∼ t →∞ t [24].In Figure 1, we report V stat ( F ) together with the D´emoulin et al. prediction (Eqs.(C1, C2)) for a similar staggeredbundle of rigid filaments at ǫ = 0 in the same conditions [12]. This comparison shows that the theoretical predictionof V stat ( F ) represents quite accurately (the difference never exceeding 2%) the exact results obtained between zeroload and stalling conditions.In Figure 2, we collect transient times τ micro and the diffusion coefficient of the bundle, Γ. Note the consistencywithin Γ values obtained from IWA or EWA runs. τ micro results of the order of W − except at small loads where itdiverges: in the discussion section we will come back to this important point.Figures 3 and 4 show respectively, for the stationary state, the load dependent averages g (0; F ) and k av ( F )Eq.(31, 32). D´emoulin’s predictions for the same quantities are also shown in these two Figures, confirming theirquantitative accuracy. D. Optical trap Let us start this section with an important remark: for our model, the choice of κ T appears to be completelyarbitrary, although, of course, it should satisfy at least the condition that the final equilibrium value of the length ofthe bundle is much greater than d , h ˆ L i EQ /d ≫ 1, in order to avoid boundary effects. However, we will see below thatthis choice will guarantee the equivalence of the results of the optical trap set–up against the constant load, at leastfor non diverging τ micro .Figure 5 shows time–dependent averages, F t = κ T h ˆ L i t , for optical trap relaxations computed by EWA and IWA for κ T = 0 . 25 and only by IWA for κ T = 0 . L = 5 d , while in the IWA case the filament sizes all start at j n = 6. Theresults, obtained by the two algorithms for κ T = 0 . 25, are indistinguishable, confirming the validity of the simplifiedalgorithm. The EWA algorithm has been used with a value for ǫ of 0 . 05 that clearly indicates the validity of the1 F/F s τ m i c r o ( W - ) F/F s Γ ( d W ) (a) (b) FIG. 2: (a) Load dependence of the relaxation time τ micro for a homogeneous bundle of N f = 32 rigid filaments counteringa constant load F at ˆ ρ = 2 . 5. (b) Diffusion coefficient Γ of the bundle. Blue symbols (IWA) and red symbols (EWA) refer tothe stationary part of the constant load stochastic dynamics experiment mentioned in (a). simplified computation done in the ǫ = 0 limit. Note that the plateau values give the stalling force predicted by Hill,Eq.(2), within statistical error bars. Fluctuations of ˆ L at equilibrium is given, as expected [9], by σ eqL = p ( k B T /κ T ) .The vertical bars reported in the figure represent the standard deviation, κ T σ L ( t ), associated to the fluctuation of theforce. They remain bounded along the entire curve by the equilibrium value, indicating a limited fluctuation betweenindividual trajectories ˆ L ( t ). That’s relevant because experiments performed by an optical trap set–up usually referto single trajectory measurement [11], whose validity is guaranteed by the smallness of fluctuations.As Eq.(13) refers to a Markov process, one expects an asymptotic relaxation of h L i t as h L i t = h L i EQ + A exp ( λ t ) + · · · = h L i EQ + A exp (cid:18) − tτ OT (cid:19) + . . . (33)where A is the amplitude (dependent on initial conditions) of the slowest, non–zero, mode with eigenvalue λ = − τ OT of the generator matrix governing the dynamics of the system. In the same long time limit, one has h F i t = κ T h L i t = F s + A κ T exp ( λ t ) + . . . (34) h V i t = A λ exp ( λ t ) + . . . (35)and thus, formally one can express the longest relaxation time of the optical trap relaxation as τ OT = − λ − = κ T − lim t →∞ F s − h F i t h V i t (36)From the data in Figure 5 one gets for κ T = 0 . 25 and from EWA trajectories τ OT = 1185 ± 50, while from IWA τ OT = 1164 ± 10. For the only IWA case at κ T = 0 . τ OT = 654 ± 10. By numerical differentiation we havecalculated the slopes of h L i t ≡ l ( t ; κ T ), h V i t = d h L i t dt ≡ v ( t ; κ T ). Eliminating t from the pair of parametric equations[ h F i t = κ T l ( t ; κ T ) , v ( t ; κ T )], we can get the velocity as a function of the force, still a function of κ T . The force–velocity relationship for κ T = 0 . V stat ( F )previously established for the constant force load stationary state. Identical results are obtained from the relaxationwith κ T = 0 . 25 (not shown), indicating a weak dependence, if any, of the force–velocity relationship on κ T .2 F (k B T/d) g ( ) FIG. 3: Values at k = 0 of the relative size distribution, g(0), as a function of the external load with N f = 32, ˆ ρ = 2 . ǫ = 0. The red filled circles are obtained in the stationary regime of constant load runs at the shown values of F . Error barsare less than symbol sizes. The dashed green line is the D´emoulin et al. estimate of g (0; F ) based on Eq.(C2). The bluecontinuous curve is obtained for the optical trap by eliminating from g (0 , t ) = h ˆ g (0) i t and h F i t = κ T h L i t at κ T = 0 . t . E. Adiabaticity Figure 1 shows that, except at low forces (short time part), V stat ( F ) superposes well the velocity force relationshipextracted from the optical trap relaxation, independently of κ T . The identity between the stationary force–velocityrelationship with the one obtained by the relaxation process in the optical trap set–up is a clear indication of the factthat the optical trap set–up is working in adiabatic conditions, i.e. that we have a relaxation process happening inbetween stationary states. We can derive from this apparent adiabaticity, especially valid at long times when the loadchanges slowly in time, that τ OT ≡ κ T − lim t →∞ F s − h F i t h V i t = − " κ T (cid:18) ∂V stat ∂F (cid:19) F s − ≡ γκ T (37)where V stat ( F ) is the constant load force–velocity relationship and where γ , defined as minus the inverse of theslope of V stat ( F ) at stalling in Eq.(37), γ = − (cid:20)(cid:16) ∂V stat ∂F (cid:17) F s (cid:21) − , is a friction coefficient having a chemical (and nothydrodynamic) origin. The structure of the relaxation time expression Eq.(37) resembles that of an overdampedbrownian oscillator.Eq.(37) can be tested with our data. Using τ OT estimates mentioned earlier for the two values of κ T , we get threecompatible γ estimates (291 ± κ T = 0 . 25, 296 ± 12 for EWA at κ T = 0 . 25, and 295 ± κ T = 0 . γ = 293 ± V stat ( F ) at stalling. The numerical derivative estimated with our too spread data gives γ = 272;unfortunately this value is not sufficiently precise to be completely reliable. Certainly, the uncertainty provided bycomputing the left and right incremental ratios giving respectively γ = 202 and γ = 414 tells us that we are withinthe numerical uncertainty. As for the D´emoulin result, its approximate estimate of the slope leads to γ Dem = 281 . γ to the value characteristic of the mean field force–velocity relationship3 F (k B T/d) k a v FIG. 4: Average filament relative size for N f = 32 at ˆ ρ = 2 . ǫ = 0. The red filled circles, denoting k av , are obtainedin the stationary regime of constant load runs at each shown load value F . Error bars are less than symbol sizes. The dashedgreen line is the D´emoulin et al. theoretical estimate of k av ( F ) based on Eqs.(C4). The blue continuous curve is obtained forthe optical trap set–up by eliminating from h k i t = h P ∞ k =0 k ˆ g ( k ) i t and h F i t = κ T h L i t at κ T = 0 . t . t (1/W ) < F > t ( k B T / d ) ε=0.05 κ T =0.25, τ O.T. =1184 ε =0 κ T =0.25, τ O.T. =1164 ε =0 κ T =0.4511, τ O.T. =654 FIG. 5: Non–equilibrium relaxations of staggered bundles of N f = 32 rigid filaments growing in an optical trap at reduceddensity ˆ ρ = 2 . 5, all of them starting from initial conditions with the wall set to a value L ≈ d ≪ h L i EQ . The wall position h L i t and the associated root mean square deviation σ L ( t ) are found as a function of time in the figure where what iseffectively shown is the load evolution κ T h L i t and corresponding RMSD κ T σ L ( t ). The final plateau value of the relaxations iscompatible with the value F s = F Hs given by Eq.(2) indicated by an horizontal thin black line. The dashed lines represent thebest fit of an exponential asymptotic behavior Eq.(34), providing estimates of τ OT and hence of the chemical friction γ defined by Eq.(37). We find γ = 291 ± κ T = 0 . γ = 295 ± κ T = 0 . τ OT values obtained. γ MF = N f k B Td W we can define a new adimensional coefficient, C ( N f , ˆ ρ ), as C ( N f , ˆ ρ ) = γγ MF = d W N f k B T γ = 9 . ± . C = 1 not only in the MF case,but also in the single filament brownian ratchet in the ǫ = 0 limit because the force–velocity relationship is identicalto the MF expression for N f = 1. The D´emoulin estimate of C gives in our case ( N f = 32 , ˆ ρ = 2 . C Dem = 8 . τ micro . The characteristic time of the optical trap equilibration is τ OT . We have seen, in Figure 2, that, for F/F s > . 15, the typical microscopic relaxation time τ micro lies in the range ≈ (1 ÷ W − . Now we can explain what wehave anticipated at the beginning of this section: with the values we have chosen for κ T , corresponding to equilibriumsizes of the bundle well satisfying the condition L EQ /d ≫ τ OT result automatically to be two to three order of magnitude larger (see the valuesgiven in Fig.5). It is important to stress, however, that τ micro values diverge as F/F s → 0, a property paralleled bythe divergence in the same limit of k av . IV. CONCLUDING REMARKS In this work, we have considered, in a Markovian approximation, a stochastic dynamical model to compute theevolution and the statistical properties of a staggered bundle of N f rigid living filaments growing against a loadedwall. In the Fokker–Planck equations we have written down to give an explicit dynamics to our system, a parameter, ǫ = τ D /τ chem , plays a special role. Generally, the model has to be solved for values of ǫ relatively small. It is foundthat if we take the ǫ = 0 limit, the dynamics simplifies and the overall computations become much lighter. We haveshown numerically that the results obtained in a reasonable range of non–zero values of ǫ in the neighborhood ofzero, coincide with the results obtained using the limiting model and the simplified algorithm. This indicates therobustness of the ǫ = 0 limit. As a consequence, the major part of the computational work of the present paperhas been performed in this limit. As we have told before, for the loading of the wall, we considered two classicalprotocols: a constant load or an optical trap set–up, characterized by a harmonic restoring force. By a series ofcomputer experiments in the case of a constant load and by only one suitable relaxation calculation in the opticaltrap, we have obtained for the two protocols the classical force–velocity relationship. With the exception of the regionof very weak loads, we have found perfect coincidence of the results. We have been able to explain this universalityof the response of the system as a result of the time scale separation between the relaxation time needed by thewall to adjust to a change of the external force and the characteristic time needed by the chemistry to change theconformation of the bundle. This condition is violated when the load is very small and that is why the optical trapand constant load results differ, even dramatically, in that region. Our results suggest that experiments measuring theforce–velocity relationship with a harmonic load offer in principle, many advantages over the approach where constantforce set–ups are used. Indeed, only a single sample is needed to get a V ( F ) estimate over a large F window in thefirst case while a separate experiment and in general a specific sample is needed for each steady state at constantload F investigated. Alternative protocols are possible, like imaging techniques used in reference [10], but the rulesneeded to get adiabaticity are easily transposed. We have been also able to confirm the large scale validity of theapproximate theory developed by D´emoulin et al. [12] to compute the properties of our system.In this work, we have only considered rigid filament. Interpreting experimental data with rigid models impliesthat the semi–flexible character of living biofilaments has limited influence on the results. How the bundle dynamicsis affected by the flexibility is a delicate point, which is largely unknown and this, to some extent, hampers theconfidence in interpreting data with rigid filament models. Work is in progress to clarify the influence of flexibility onthe force–velocity relation. ACKNOWLEDGMENTS We thank G. Kozyreff for his help with the perturbation expansion of the F.P. equation in section 2B. Two of us(CP and JPR) are grateful to J. Baudry, J.F. Joanny et D. Lacoste for useful discussions. We thank G. Destr´ee fortechnical help. JPR thanks the University of L’Aquila for hospitality during a three months visit. CP is supported bythe Agence Nationale de la Recherche (ANR) under the project “HyLightExtreme”. AP is supported by a Mobility5Grant for PhD students from Sapienza University of Rome, and thanks the ENS for hospitality during a six monthsvisit. Appendix A: Discretized Fokker-Planck Equation for the wall–bundle system in an optical trap or constantforce set–up In this appendix we derive a proper discretization of the Fokker–Planck equation together with the elements of thegenerator matrix Q of the Markov process given by Eq.(13).To get the matrix elements which account for the discretization of the variable L , following the procedure introducedin [18], we concentrate only on the diffusive part of Eq.(7) for the wall position probability at given chemical state, P j ,...,j Nf ( L, t ) ≡ P j ( L, t ): ∂P j ( L, t ) ∂t = − ∂∂L J j ( L, t ) (A1)where J j ( L, t ) = − D (cid:16) ∂P j ( L,t ) ∂L + k B T d Φ dL P j ( L, t ) (cid:17) is the probability current, with d Φ dL = κ T L or − F for respectivelythe optical trap or the constant force set–up. We define the probabilities for the wall to be in the intervals ( l = L/δ ) k − / l < k + 1 / k + 1 / l < k + 3 / p k ( t ) = Z k +1 / k − / P j ( l, t ) dl (A2) p k +1 ( t ) = Z k +3 / k +1 / P j ( l, t ) dl. (A3)By defining the wall forward rate F k +1 / of going from k to k + 1 ( F k − / from k − k ) and the wall backward rate B k +1 / of going from k + 1 to k ( B k − / from k to k − p k ( t ) can be writtenas [18]: dp k ( t ) dt = F k − / p k − − ( F k +1 / + B k − / ) p k + B k +1 / p k +1 = − ( F k +1 / p k − B k +1 / p k +1 ) + ( F k − / p k − − B k − / p k )= − ( J k +1 / − J k − / ) (A4)where the rates F k ± / and B k ± / have to be derived by discretizing Eq.(A1). J k +1 / is the net probability fluxbetween sites k and k + 1 ( J k − / is between k − k ).If we now discretize Eq.(A1) using e.g. the central difference method ( f ′ k +1 / = ( f k +1 − f k ) /δ ) and comparethe resulting discrete equation with Eq.(A4), the forward and backward rates obtained will not respect the detailedbalance, a sufficient condition to reach equilibrium, while we expect the evolution of the Markov chain to lead to it,with each process balanced by its reverse.To overcome this difficulty, following [18], we can look for the stationary solution of Eq.(A1) and see if, by integrationover a proper interval of lengths, we can obtain an identification of the rates, bringing us to coefficients satisfying thedetailed balance.Looking at the definitions Eqs.(A2, A3), we see that to get p k and p k +1 from a solution of the stationary equation(A1) we need to solve it in the interval ( k − / , k + 3 / P j ( L, t )in terms of the stationary solution P EQ ( l ) of: D ddl (cid:18) dP EQ ( l ) dl + ∆Φ k +1 / k B T P EQ ( l ) (cid:19) = 0 (A5)in l ∈ ( k − / , k + 3 / d Φ /dl by the constant approximation ∆Φ k +1 / , with∆Φ k +1 / = Φ( k + 1) − Φ( k ) (A6)The general solution of Eq.(A5) is P EQ ( l ) = η exp (cid:16) − ∆Φ k +1 / k B T l (cid:17) + θ with η and θ constants. Plugging this expressioninto Eqs.(A2, A3), one can easily find η and θ in terms of p k and p k +1 . Then the (approximate) stationary solution6of the Fokker–Planck equation for the wall in the interval ( k − / , k + 3 / 2) is: P EQ ( l ) = ∆Φ k +1 / ( p k − p k +1 ) k B T (cid:16) exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − (cid:17) exp (cid:18) ∆Φ k +1 / k B T ( k − / (cid:19) exp (cid:18) − ∆Φ k +1 / k B T l (cid:19) + p k exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − p k +1 (cid:16) exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − (cid:17) l ∈ ( k − / , k + 3 / 2) (A7)From this equation we get the probability flux in the same interval: J EQ ( l ) = − e D dP EQ ( l ) dl − e D ∆Φ k +1 / k B T P EQ ( l ) = − e D ∆Φ k +1 / k B T p k exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − p k +1 (cid:16) exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − (cid:17) (A8)with e D = D/δ the diffusion constant in δ units. Comparing this current with the probability flux defined in Eq.(A4),we get the following forward and backward rates: F k +1 / = e D ∆Φ k +1 / /k B T exp (cid:16) ∆Φ k +1 / k B T (cid:17) − B k +1 / = e D − ∆Φ k +1 / /k B T exp (cid:16) − ∆Φ k +1 / k B T (cid:17) − . (A10)The same approach for the interval ( k − / , k + 1 / 2) can be used to get F k − / and B k − / .By direct substitution, we see that Eqs.(A9, A10) respect the detailed balance, F k +1 / P EQ ( k ) = B k +1 / P EQ ( k + 1).Substituting the appropriate expression for Φ( L k ), we have:∆Φ k +1 / = ( F δ constant load κ T δ (cid:0) ( k + 1) − k (cid:1) optical trap (A11)The non–zero elements of the generator matrix Q can now be written as follows: Q { j ,...,j n ,...,j Nf ,k }{ j ,...,j n ,...,j Nf ,k +1 } = F k +1 / (A12) Q { j ,...,j n ,...,j Nf ,k }{ j ,...,j n ,...,j Nf ,k − } = C k − / = ( B k − / if k − > X ∗ /d Q { j ,...,j n ,...,j Nf ,k }{ j ,...,j n +1 ,...,j Nf ,k } = U j n = ( U if k > X n ( j n + 1) /d Q { j ,...,j n ,...,j Nf ,k }{ j ,...,j n − ,...,j Nf ,k } = W j n = W (A15) Q { j ,...,j n ,...,j Nf ,k }{ j ,...,j n ,...,j Nf ,k } = − F k +1 / − C k − / − N f X n =1 ( U j n + W j n ) (A16)with X n ( j n ) and X ∗ given by Eqs.(4, 5). The row sums of this matrix are zero, as required for a generator matrix ofa Markov chain: X { j ′ ,...,j ′ n ,...,j ′ Nf ,k ′ } Q { j ,...,j n ,...,j Nf ,k }{ j ′ ,...,j ′ n ,...,j ′ Nf ,k ′ } = 0 . (A17)Eq.(13) represents hence a continuous time Markov process with discrete states; as for the variable L , the discretestates are approximations (exact in the δ → Appendix B: Elements of the ǫ = 0 generator matrix In this appendix we write explicitely the matrix elements of Q (0) , generator of the Markov process in the ǫ = 0limit Eq.(29). Since in this limit the integration in L allowed us to get rid of the continuous wall diffusion process,these elements can be written immediately: Q (0) { j ,...,j n ,...,j Nf }{ j ,...,j n +1 ,...,j Nf } = U j n = U A ( n ) ( j , . . . , j n , . . . , j N f ) (B1) Q (0) { j ,...,j n ,...,j Nf }{ j ,...,j n − ,...,j Nf } = W j n = W (B2) Q (0) { j ,...,j n ,...,j Nf }{ j ,...,j n ,...,j Nf } = − N f X n =1 ( U j n + W j n ) (B3)where A ( n ) ( j , . . . , j n , . . . , j N f ) is given in by: A ( n ) ( j , . . . , j n , . . . , j N f ) = Z ∞ X ∗ /d dx Θ ( x − X n ( j n + 1) /d ) e P EQ ( x | j , . . . , j n , . . . , j N f )= exp h − f (cid:16) X ∗ ′ − X ∗ (cid:17) /d i constant loaderfc h ( e κ T / / X ∗ ′ /d i erfc h ( e κ T / / X ∗ /d i optical trap (B4)where X ∗ ′ is the most advanced filament’s tip for the set of filament sizes { j , . . . , j n + 1 , . . . , j N f } . Eq.(B4) has beenderived previously for constant load [7, 12] and for optical trap load [17]. Appendix C: D´emoulin et al. prediction for V ( F ) and k av D´emoulin et al. [12] have proposed an approximate solution for the force–velocity relationship of staggered rigidfilaments subjected to a constant load F in the ǫ = 0 limit. They found: V ( F ) = dU N f N f exp (cid:18) − F dk B T (cid:19) + N f − X m =1 g (0) ( N f − m ) exp (cid:18) − F d ( N f − m ) N f k B T (cid:19) − dW N f g (0) N f − X m =1 m (1 − g (0)) m − + N f (1 − g (0)) N f − (C1)with the relative size distribution g ( k ) given by: g ( k ) = d ( U − W ) − VdU (cid:18) V + dW dU (cid:19) k k = 0 , ∞ . (C2)It can be verified that at stalling, F = F s = N f k B Td ln ˆ ρ , one gets V = 0 and g (0) = 1 − ˆ ρ − .For the comparison in the text, we need to compute k av and V ( F ) explicitely:1. k av : defining ξ = V + dW dU to simplify expressions, one gets from Eq.(C2): g ( k ) = (1 − ξ ) ξ k (C3) k av = ∞ X k =1 kg ( k ) = ξ − ξ (C4)2. V ( F ): Inserting g (0) = d ( U − W ) − VdU into Eq.(C1), we find for V ( F ) a polynomial equation in V to solve. Writing8 v = V /dW , we find: φ ( v ) = ˆ ρ N f N f exp (cid:18) − F dk B T (cid:19) + N f − X m =1 (cid:18) − v ˆ ρ (cid:19) ( N f − m ) exp (cid:18) − F d ( N f − m ) N f k B T (cid:19) − N f (cid:18) − v ˆ ρ (cid:19) N f − X m =1 m (cid:18) v ˆ ρ (cid:19) m − + N f (cid:18) v ˆ ρ (cid:19) N f − − v = 0 (C5)Eq.(C5) can be solved numerically using the Newton–Raphson method, for which the derivative of φ ( v ) withrespect to v is needed: φ ′ ( v ) = − ˆ ρ N f N f − X m =1 ( N f − m ) ( N f − m ) exp (cid:18) − F d ( N f − m ) N f k B T (cid:19) (C6)+ 1 N f ˆ ρ N f − X m =1 m (cid:18) v ˆ ρ (cid:19) m − (cid:18) v ˆ ρ − (cid:18) − v ˆ ρ (cid:19) ( m − (cid:19) − N f ( N f − (cid:18) − v ˆ ρ (cid:19) N f − − v , which can be taken as e.g. the Hill’s value.The solution of D´emoulin et al. equation for a bundle of N f = 32 filaments and supercritical density ˆ ρ = 2 . V ( F ) in Eqs.(C3, C4), for g (0) and k av in Fig.3 and 4 respectively.To predict within the present theory the value of τ OT , we need to compute the derivative of V ( F ) with respectto F at stalling, obtaining γ Dem . From Eq.(C1), we get dV /dF as an implicit function of V ( F ) and F . At stalling F = F s , V ( F s ) = 0, we obtain: ∂V∂F (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) F = F s = − d W k B T ˆ ρ − N f h P N f − m =1 ˆ ρ m (1 − ˆ ρ − 11 ( N f − m ) N f i N − f P N f − m =1 m ˆ ρ − m + ( N f − ρ − ( N f − − N − f P N f − m =1 m ˆ ρ − m [1 − (ˆ ρ − m − γ Dem = − ( dV /dF ) − s = 281 . 6, in agreement with our results (see main text). [1] T. Risler Cytoskeleton and Cell Motility in Encyclopedia of Complexity and Systems Science , Springer, 1738-1774, (2009).[2] S.H.Parekh, O. Chaudhuri, J.A. Theriot and D.A.Fletcher Nature. Cell biology , 1219 (2005).[3] D.B. Smith and J. Liu, Phys. Biol. , 016004 (2013)[4] C.H. Schreiber, M. Steward and T. Duke, Proc. Natl. Acad. Sci. USA (20), 9141 (2010).[5] T. L. Hill, Proc. Natl. Acad. Sci. USA (9), 5613 (1981).[6] T.E. Schaus and G.G. Borisy, Biophys. J. , 1393 (2008).[7] G. Sander van Doorn, C. Tanase, B.M. Mulder and M. Dogterom, Eur. Biophys. J. , 2 (2000).[8] T.L. Hill and M. W. Kirschner, Int. Review of Cyt. , (1982), 1-125.[9] A.Perilli, C. Pierleoni, G. Ciccotti and J.P. Ryckaert, J. Chem. Phys. , 245102 (2016).[10] M. Dogterom and B. Yurke Science , 856 (1997).[11] M. J. Footer, J.W.J. Kerssemakers, J.A. Theriot and M. Dogterom, Proc. Natl. Acad. Sci. USA , 2181 (2007).[12] D. D´emoulin, M-F. Carlier, J. Bibette, and J. Baudry, Proc. Natl. Acad. Sci. , 17845 (2014).[13] A. Mogilner and G. Oster, Eur. Biophys. J. , 235 (1999).[14] K. Tsekouras, D. Lacoste, K. Mallick, and J.-F. Joanny, New J. Phys. , 103032 (2011).[15] C. S. Peskin, G. M. Odell and G. F. Oster, Biophys. J. , 316 (1993).[16] R. Wang and A. E. Carlsson, New Journal of Physics , 113047 (2014).[17] A. E. Carlsson, Phys. Biology , 036002 (2008).[18] H. Wang, C. S. Peskin and T. C. Elston, J. theor. Biol. , 491 (2003).[19] D. T. Gillespie, The Journal of Physical Chemistry (25), 2340 (1977).[20] P. Gaspard and E. Gerritsma, J. Theor. Biol , , 672 (2007).[21] J. R. Norris, Markov chains , Cambridge University Press (1997).[22] D. T. Gillespie, Markov Processes: An Introduction for Physical Scientists. San Diego, Academic Press, INC. (1992). [23] G. Kozyreff, “Brownian ratchet with a linear restoring force” , unpublished notes (2016).[24] P. Ranjith, D. Lacoste, K. Mallik and J-F Joanny Bioph. Journal96