On the force–velocity relationship of a bundle of rigid living filaments.
Abstract
In various cellular processes, biofilaments like F–actin and F–tubulin are able to exploit chemical energy associated to polymerization to perform mechanical work against an external load.
The force–velocity relationship quantitatively summarizes the nature of this process.
By a stochastic dynamical model, we give, together with the evolution of a staggered bundle of rigid living filaments facing a loaded wall, the corresponding force–velocity relationship.
We compute systematically the simplified evolution of the model in supercritical conditions at , where is the monomer size, is the obstacle diffusion coefficient, and are the polymerization and depolymerization rates.
Moreover, we see that the solution at 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 a harmonic restoring force.
The constant force case leads, for each value, to a stationary velocity after a relaxation with characteristic time .
When the bundle (initially taken as an assembly of filament seeds) is subjected to a harmonic restoring force (optical trap load), the bundle elongates and the load increases up to stalling (equilibrium) over a characteristic time .
Extracted from this single experiment, the force–velocity curve is found to coincide with , except at low loads.
We show that this result follows from the adiabatic separation between and , i.e. .
Submitted to the Journal of Chemical Physics on August 22nd 2017
m\quinn_mref:n #1 \seq_new:N ł_quinn_mref_seq \cs_new:Npn \quinn_mref:n #1 \seq_set_split:Nnn ł_quinn_mref_seq , #1 \seq_pop_right:NN ł_quinn_mref_seq ł_tmpa_tl ( \seq_map_inline:Nn ł_quinn_mref_seq LABEL:##1, exp_args:NV LABEL:ł_tmpa_tl ) \ExplSyntaxOff
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 monomers polymerize at the barbed end of filaments, directly in contact with the cytoplasmic membrane. The speed of the membrane deformation/displacement at the leading edge of the cell adjusts itself so that the force generated by the growing filaments compensates the resisting load coming from the membrane tension and to the crowded environment around the cell. For living filaments opposing a loaded mobile obstacle, the macroscopic force–velocity relationship, , linking the obstacle velocity, , only to the instantaneous applied load, , quantitatively summarizes the combined action of the elementary self–assembling processes. In such adiabatic conditions, implying a time scale separation between the self–assembling process and the response of the obstacle, the dependence could be probed 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 sample grows against an AFM cantilever or an optical trap), where the obstacle velocity can be followed as the load increases continuously up to stalling.
Abiabatic conditions cannot be in general guaranteed and indeed, careful investigations on a branched actin network growing against an AFM tip have shown that the recorded relationship can be function of the load history [[2]]. The direct force–velocity relationship, , is in any event widely used as a characteristics of network dynamics to compare 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–velocity relationship, , 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 incompressible bundle of parallel filaments pressing against a mobile obstacle [[5]]. Successively, when the filaments of the bundle are treated as independent and equivalent and when it is assumed that the depolymerization rate is unaffected by the external load, the wall velocity ( indicates the mean field character of this treatment) has been written as [[6], [7]]
(1) 
where and are the single filament bulk rate constants, related to bulk chemical rate constants and , for single monomer polymerization and depolymerization steps, is the free monomer density, is the single filament increment of contour length per incorporated monomer and is the external force exerted on the wall. Supercritical conditions, where filament polymerization dominates over depolymerization, require where is the critical value of the monomer density at which the bundle has no tendency to grow nor to shrink in absence of load.
Eq.(1) predicts, for , a growth velocity of the free bundle while the stalling force , at which the velocity vanishes, is given by
(2) 
The notation 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 equilibrium optical trap ensemble, multiplied by , converges exponentially fast for to the Hill’s prediction.
Experimental measurements of the force–velocity relationship for multi–filament bundles (tubulin or actin) [[10], [7], [11], [12]] are not many, reflecting the difficulty to prepare in–vitro the grafted bundle seed needed to follow its subsequent loaded growth. However, it is interesting to note the diversity in these few approaches. The growth of single grafted tubulin filaments, which are bundles of proto–filaments, was followed by imaging techniques [[10], [7]]. Regrouping data for different observation times and for different samples, a master force–velocity relation could be established. In another experiment using an acrosome bead complex of 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 plateau was observed but the final stationary force was surprisingly much lower than the expected stalling force, Eq.(2), its value being close to the stalling force predicted for a single filament. The analysis in this experiment considers many relaxation curves, but in many cases data had to be eliminated due to interferences during the relaxation process with the onset of escaping filaments. This happens because growing filaments undergo a large bending fluctuation which allows them to start growing freely along the obstacle. The transient behavior, which can be converted into a law by the derivation in time of the wall position, was not exploited. Finally, in a recent study [[12]], recording the rate of radial distance between two colloidal particles separated by a growing grafted actin bundle, the force–velocity relationship of actin bundles was established in constant load conditions.
The outcome of the earliest experiment [[10], [7]], confirmed by the more recent experimental work [[12]], is that the velocity, 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 in the bundle, with the effect of reducing the additivity of the action of each filament. The bundle model needs to be specified and the dependence between chemical events and wall position for a given longitudinal seeds disposition has to be quantitatively taken into account. This aspect is present in the multi–filament Brownian Ratchet (BR) models [[7], [12], [13], [14]] which generalize the single filament brownian ratchet model introduced by Peskin et al. [[15]], for which one finds that the velocity vanishes for a load equal to Hill’s expression, Eq.(2) [[7]]. For these bundle models, the important characteristics which distinguish the dynamical behavior of the bundle are the number of rigid living filaments, the longitudinal disposition of the seeds of the filaments and the wall diffusion coefficient which introduces a second characteristic time next to the chemical events time scale . This fact suggests to introduce the parameter to be able to discuss the condition of this second adiabatic separation (not to be confused with the one associated to the existence of ). For both experiments having probed the relationship, it was found that data could be interpreted successfully with the model of a staggered bundle (= staggered longitudinal seed disposition [[9]]) of rigid filaments in very fast wall diffusion conditions (), a model we will denote as SRBR (Staggered Rigid Brownian Ratchet). On the contrary, for a similar model with an in registry bundle (unstaggered longitudinal seed disposition) [[14]], the predicted velocity was much too low with respect to the experimental data [[10], [7], [12]].
In the stochastic dynamical models here considered, the force–velocity relationship depends parametrically, for a given 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, , for our models has been established by stochastic dynamics studies at finite [[16]] and at [[7], [12]]. In the latter case a simplified algorithm, exploiting the time scales separation, has been used for the staggered bundle case. Indeed, at , the wall position distribution at given filaments configuration, is found to be time–independent and equal to the equilibrium distribution of the wall position resulting from the 1D Brownian motion of a wall in the external load field, 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 (), successive theoretical developments [[13], [7], [12]] have given, with a very good approximation [[12]], two coupled closed expressions for the velocity and the distribution of filament relative sizes (see Section III) , for the stationary state.
In this work, we consider the stochastic staggered bundle model of rigid filaments in supercritical conditions and perform a series of dynamical runs for different load conditions. We first look at the constant force case, treating both the stationary state itself and the asymptotic transient evolution to reach it. We next envisage the bundle, in similar thermodynamic conditions, initially taken with very short filaments, subject to a harmonic load , where is the wall position and where 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 relaxation time with a theoretical approximate expression derived along the lines of the Démoulin et al. theory. We derive and compare the force–velocity relationship extracted from this optical trap relaxation with the one obtained in stationary conditions. As expected, we found that the two coincide in adiabatic conditions, i.e. when the characteristic time of the 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 an optical trap load, while most studies (with an exception restricted to the case [[17]]) assume a constant load. Moreover, while algorithms for finite or are usually just assumed, we establish an explicit link showing how the 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 seed disposition, facing either a constant or a harmonic load and we derive the explicit wall algorithm (EWA) giving the sampling rules to generate stochastic trajectories for any finite case. We then use a perturbation expansion to derive the 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 case. Section III reports and discusses our results for constant force and optical trap loads for the same bundle system generally using the approach, since the wall diffusion takes place very quickly with respect to the mean time between (de)polymerization events. However, we also verify that the simplified algorithm is robust, since we find identical results in a reasonable range of non–zero values. Section IV concludes with a summary of the main results and with some perspectives.
Ii Model and Implementation
We consider a bundle of living filaments, grafted normally (say along the axis) to a fixed planar substrate wall (along and directions). The filaments are modeled as discrete rigid linear chains with monomer size and length related to the number of attached monomers, , as . Let be the location along axis of the seed (first monomer) of the filament close to the grafting plane (). For a bundle of many filaments, two seed dispositions are usually considered: in–registry (or unstaggered), where , , and homogeneous (or staggered), where seeds are regularly spaced as
(3) 
A moving obstacle, a hard wall located at distance from the parallel substrate wall, is loaded with a compressional external force bringing it into contact with the living filaments. We will consider two types of load, the constant force , and the optical trap setting, with , where is the trap stiffness and the distance between the walls.
The bundle force for rigid filaments is impulsive. Its effect is taken into account by imposing a confining boundary to the wall motion at the tip location of the longest filament.
Filaments either grow by a single monomer polymerization step with bulk rate , proportional to the free monomer density , or shrink by a single monomer depolymerization step with bulk rate . The ratio is the free monomer density divided by its critical value, i.e. the value at which the two bulk rates are equal. We will be interested to supercritical conditions only (), where the filaments tend to grow against the loaded wall. When a filament tip gets closer than to the wall, the polymerization rate becomes zero, while the depolymerization one is assumed to remain unchanged.
The dynamics of the bundle of growing filaments against the loaded mobile wall presents two main time scales: the chemical one, , and the diffusive one, related to the diffusive motion of the wall and estimated by where is the wall diffusion coefficient. The ratio of time scales, , in typical in–vitro experiments is , 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 independent rigid filaments subjected to a constant or harmonic load, for arbitrary value of .
ii.1 General Fokker–Plank equation for a bundle of rigid filaments against a constant or harmonic load
We describe the time evolution of filaments against a load in terms of the filament sizes and the wall position, . The wall position must always lie beyond the tip of any filament – and so beyond the tip of the most advanced one, with size . Defining , the position of the tip of filament and that of the most advanced one,
(4)  
(5)  
(6) 
We assume that the joint probability distribution function satisfies a Fokker–Planck equation in time mixing a continuous process in space for the wall position with a discrete process for filament sizes. For the model described above, we have
(7)  
where is the Heaviside step function and the probability current density is
(8) 
In Eq.(8) the compressive force can be either a constant or an elastic force modeling the optical trap. The right–hand side of Eq.(7) represents the sink and source terms affecting the dynamics due to polymerization and depolymerization events. Their explicit expression indicates that, in one step at fixed , transitions are only possible between adjacent microscopic states, where filaments have the same size while the size of the remaining filament differs by unit, taking into account the restriction , and that the filament size cannot be smaller than two.
The general normalization condition for the distribution is
(9) 
while the boundary conditions on the probabilities are
(10)  
(11) 
To simplify the treatment of the continuous–discrete structure of Eq.(7), we discretize, following reference [[18]], the wall position with a grid step , with , integer, . We then substitute to the wall position the discrete variable
(12) 
In this way Eq.(7) will become a finite difference equation in representing a discrete Markov chain in continuous time
(13) 
with a vector field and the generator matrix of the Markov chain. The elements of the matrix contain the (de)polymerization rates for the filaments,
(14)  
(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 discrete Markov chain using any appropriate algorithm, in our case the Gillespie algorithm [[19], [20]]: given an initial condition at time , the state of the system is estimated in terms of the set of random variables at time producing statistically correct trajectories, from which the probability distribution function can be inferred by histograms. Starting from the initial state, the system is allowed to evolve by random steps involving only one reaction per time: one filament depolymerization or polymerization, or the wall forward or backward jump. Denoting by the current state of the system, the reachable states are those differing from for only one variable by , namely or . It is straightforward to see that the number of these possible final states is . The transitions , , are described in Eq.(13) by the generator matrix elements , the rates of going from to . The corresponding diagonal element is [[21]]. The evolution of the system is determined by two random variables: the time to the next reaction, , and the final state , or equivalently the index of the reaction, . From general Markov chain theory, is known to be an exponentially distributed random variable: given the current state , the parameter of the exponential distribution is given by . Instead, the probability for the jump linking states and to take place is given by the ratio between and [[21]]. The main loop of the algorithm follows this scheme:

The initial state is specified in terms of the state vector . We take for the initial value of a small fixed value and, for the filament, compatible initial sizes;

The matrix elements are calculated for any state reachable from ;

The time to the next move is determined using the so–called direct method, which follows from the standard inversion method of the Monte Carlo theory [[22]]: a random number is generated from the uniform distribution and the time is taken as
(16) 
The index of the next move is determined using the same method: a second random number is generated and the index is taken as the smallest integer satisfying
(17) 
The sampled move is taken by updating the state vector and the time is incremented by ;

Go back to 1, until a maximum time is reached;

End the simulation.
The state vector is stored for the calculation of histograms and averages.
This algorithm [[19], [18]], solving the Fokker–Planck Eqs.\mrefFPeq,eq:pcdF, works for any seed disposition (staggered and unstaggered), for any finite value of the dimensionless parameter and both for the two cases of constant force and 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 the limit . In this limit the wall reequilibrates instantaneously after any change of the position of the most advanced tip of the bundle. The interest of this limit is justified since in experiment with actin bundles / colloidal particles (e.g. the optical trap experiment [[11]]) the typical value of the ratio of time scales is . We will see that the dynamics of the bundle simplifies for two reasons: the elimination of the fast motion of the wall permits to go to longer times and the dimensionality of the problem is reduced. The new algorithm, called implicit wall algorithm (IWA), becomes then decidedly more efficient.
ii.2 Treatment of the Fokker–Planck equation in the limit [[23]]
Given the separation of time scales between the chemical events and the wall diffusion, it is convenient to rewrite Eq.(7) in terms of dimensionless variables in order to put in evidence the ratio . Defining , and , multiplying Eq.(7) by and redefining the probability distribution functions, we get:
(18)  
with
(19) 
the probability current density in the reduced units. In the limit, it is legitimate to replace Eq.(18) by its simpler zero–th order approximation:
(20) 
By integrating in form to and using the boundary conditions for the probability, one gets so that the general solution is:
(21) 
On the other side, it is always possible to write the joint probability as the product of the marginal distribution for the subset times the conditional probability distribution for :
(22) 
Therefore, given the general solution (21), we can write it as
(23) 
as the dependence, Eq.(21), is explicit and time–independent. The wall distribution is an explicit, 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:
(24) 
with . From Eq.(24) we get the average wall position conditional to the bundle sizes as:
(25) 
Note that the full distribution at , given by Eq.(23), is still a time–dependent function, since filament sizes change by single monomer polymerization/depolymerization events; the infinite separation of the time scales () implies that after any chemical event, the wall immediately re–equilibrates according to the time–independent distribution, Eq.(24), given the new set of filament sizes. To get the full distribution, we write as an asymptotic expansion in terms of the small parameter :
(26) 
where 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:
(27)  
Integrating both sides of this equation from to , applying the boundary conditions Eq.(11) on and the normalization of and using Eq.(23), we get:
(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 a vectorial form, similar to Eq.(13):
(29) 
with 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 described above for the general Fokker–Planck equation for .
As already mentioned, in this case the algorithm is more efficient since it spans longer times (we have integrated out the fast variable) and it has to treat a reduced number of variables.
The solution of Eq.(29) and the conditional probability for the wall position Eq.(24), give the necessary information needed to compute all timedependent ensemble averages, as e.g. . Similar model and procedures have been used: i. for constant load option and inregistry [[14]] or staggered [[13], [7], [12]] bundles; ii. for optical trap only for staggered bundles [[17]].
Iii Simulations and Results
iii.1 Units, parameters and stochastic runs
In our simulations, length, time and energy units are taken as , , and respectively. All quantities will be mentioned in reduced units based on the above three fundamental units. For actin ; experimental information for gives ; and, at room temperature . We choose to perform our studies on a bundle of rigid filaments with a staggered disposition of seeds at a reduced density . 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 . Given the small value of , we performed the major part of our simulations in the limit with the IWA algorithm. However, we have considered interesting to compare the results of the IWA algorithm with those of the EWA corresponding to a finite but small value of . With the very small experimental value of , EWA would be highly inefficient, since the computer time would be essentially spent to study the wall diffusion next to a bundle with quasifixed filament sizes. Since for the loadvelocity relationship we need to sample both wall and filament sizes, we decided to adopt a value of epsilon thousand times bigger, . This value, in fact, still permits to give a sufficient representation of the wall dynamics. Our EWA approach requires to discretize the space variable with elementary steps . For , we have adopted . To compute the solution of our Fokker–Planck equation, both for (IWA) or (EWA), we need to fix initial configurations. Our choice for EWA has been to fix the wall location (i.e. ) and to sample the initial filament sizes for each trajectory of the stochastic dynamics according to the filament size equilibrium probability , conditional to the chosen wall location [[9]]. For initiating IWA runs, the initial filament sizes must be arbitrarily chosen and the initial wall location then follows from its conditional distribution.
iii.2 Observables of interest

Wall position
The wall position is the quantity directly followed in time in real experiments and corresponds to the expected value of the random variable over the solution of the FP equation, . The calculation of this quantity is direct in the EWA case, while it has to be determined in the IWA case through the instantaneous size distribution of , , implying values ahead of the tip of the most advanced filament at , given by Eq.(5), using Eq.(25).

Relative size (in number of monomers) of filaments with respect to the leading one.
In terms of the tip positions and defined by Eqs.\mrefXn,Xstar, let us define the relative subset index given by
(30) 
This index represents in successive order the filament of order nearest neighbor of , second neighbor of , etc. Therefore it gives an intrinsic order to the vector representing the relative size of each filament. Note that the dividend in the function in the given condition is always positive. Then we can define, for each filament , the quantity
(31) 
Each component of this vector represents in discrete units of monomer size the relative distance from the most advanced tip of the first, second, etc. neighboring index.
This vector of relative sizes is interesting because its timedependent probability distribution reaches a stationary value in the case of the wall subjected to a constant load.

Density of relative size of filaments with respect to the leading one
This quantity is defined by the microscopic observable
(32) 
At time , the microscopic distribution will be . Specifically, we will characterize the internal structure of the bundle either by , the average probability densitythat the tip lies at a distance smaller than from the tip of the most advanced filament, or the average relative size . We will denote by and the time–asymptotic values of these quantities for constant load force dynamics.
iii.3 Constant force load
We have computed the relaxation towards the stationary state for a homogeneous bundle of rigid living filaments at pressing against a constant load . We have chosen various values of in the range with the stalling force, Eq.(2). For each load value, we have generally used the IWA algorithm to produce independent trajectories, starting at time with all filament sizes set to the same value . We have chosen this value to avoid to fall at later times at the lower boundary . That could happen when with negative average velocities. In two cases, starting with , we have used the EWA algorithm, averaging over independent trajectories. To determine the microscopic relaxation time of the bundle, we have fitted the asymptotic time evolution of the average wall position as . To get the diffusion coefficient of the bundle , we have also fitted the asymptotic behavior of the mean square elongation [[24]].
In Figure 1, we report together with the Démoulin et al. prediction (Eqs.\mrefvf_dem,g_k) for a similar staggered bundle of rigid filaments at in the same conditions [[12]]. This comparison shows that the theoretical prediction of represents quite accurately (the difference never exceeding ) the exact results obtained between zero load and stalling conditions.
In Figure 2, we collect transient times and the diffusion coefficient of the bundle, . Note the consistency within values obtained from IWA or EWA runs. results of the order of except at small loads where it diverges: in the discussion section we will come back to this important point.
iii.4 Optical trap
Let us start this section with an important remark: for our model, the choice of appears to be completely arbitrary, although, of course, it should satisfy at least the condition that the final equilibrium value of the length of the bundle is much greater than , , in order to avoid boundary effects. However, we will see below that this choice will guarantee the equivalence of the results of the optical trap set–up against the constant load, at least for non diverging .
Figure 5 shows time–dependent averages, , for optical trap relaxations computed by EWA and IWA for and only by IWA for . In the EWA case, the relaxations start from a bundle size, short with respect to the final equilibrium value, i.e. , while in the IWA case the filament sizes all start at . The results, obtained by the two algorithms for , are indistinguishable, confirming the validity of the simplified algorithm. The EWA algorithm has been used with a value for of that clearly indicates the validity of the simplified computation done in the limit.
Note that the plateau values give the stalling force predicted by Hill, Eq.(2), within statistical error bars. Fluctuations of at equilibrium is given, as expected [[9]], by . The vertical bars reported in the figure represent the standard deviation, , associated to the fluctuation of the force. They remain bounded along the entire curve by the equilibrium value, indicating a limited fluctuation between individual trajectories . That’s relevant because experiments performed by an optical trap set–up usually refer to 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 as
(33) 
where is the amplitude (dependent on initial conditions) of the slowest, non–zero, mode with eigenvalue of the generator matrix governing the dynamics of the system. In the same long time limit, one has
(34)  
(35) 
and thus, formally one can express the longest relaxation time of the optical trap relaxation as
(36) 
From the data in Figure 5 one gets for and from EWA trajectories , while from IWA . For the only IWA case at . By numerical differentiation we have calculated the slopes of , . Eliminating from the pair of parametric equations , we can get the velocity as a function of the force, still a function of . The force–velocity relationship for , shown in Figure 1, turns out to be equivalent, except at small loads, to previously established for the constant force load stationary state. Identical results are obtained from the relaxation with