Particle-Based Modeling of Living Actin Filaments in an Optical Trap

We report a coarse-grained molecular dynamics simulation study of a bundle of parallel actin filaments under supercritical conditions pressing against a loaded mobile wall using a particle-based approach where each particle represents an actin unit. The filaments are grafted to a fixed wall at one end and are reactive at the other end, where they can perform single monomer (de)polymerization steps and push on a mobile obstacle. We simulate a reactive grand canonical ensemble in a box of fixed transverse area A, with a fixed number of grafted filaments Nf, at temperature T and monomer chemical potential μ1. For a single filament case (Nf=1) and for a bundle of Nf=8 filaments, we analyze the structural and dynamical properties at equilibrium where the external load compensates the average force exerted by the bundle. The dynamics of the bundle-moving-wall unit are characteristic of an over-damped Brownian oscillator in agreement with recent in vitro experiments by an optical trap setup. We analyze the influence of the pressing wall on the kinetic rates of (de)polymerization events for the filaments. Both static and dynamic results compare reasonably well with recent theoretical treatments of the same system. Thus, we consider the proposed model as a good tool to investigate the properties of a bundle of living filaments.


Introduction
When a bundle of parallel actin filaments in supercritical conditions hits a moving wall subject to an opposing constant load force F L , the balance between the polymerization force, the load and the friction force from the solvent (usually negligible) leads to a stationary velocity v of the obstacle. Here, v indicates a coarse-grained velocity averaged over rapid fluctuations due to microscopic events, like the addition/removal of single monomers occurring at the bundle tip and to the usual wall Brownian motion. A stationary non-equilibrium state is only possible if the living filaments are rigid so that the distribution of the tip-wall distance(s) (for single or many filaments) becomes independent on the length of the bundle [1][2][3][4] and also independent of time in a time window much longer than the obstacle diffusion time and the chemistry time (inverse polymerization rate). When filament flexibility is taken into account, both the filament length and the wall position (with respect to the origin of the filaments) become relevant variables, and the dynamical state can no longer be stationary, not even on a coarse-grained time window [5][6][7].
For rigid living polymers pressing against a mobile wall, the Brownian ratchet (BR) model, in the fast wall-diffusion limit, provides the simplest theoretical rationalization of the conversion of chemical energy into mechanical work against the load [1]. A single rigid filament originating from a grafted seed, subjected to polymerization and depolymerization steps with size increment d at its unique active end, is perturbed by a fluctuating moving wall. Any polymerization attempt (with attempt rate U 0 ) is rejected if it leads to overlap with the wall, but is accepted otherwise. The depolymerization steps occur at a uniform rate W 0 , independently of the tip-wall distance. The wall is undergoing a 1D drifted Brownian motion limited on one side by the rigid filament tip and otherwise characterized by a friction coefficient ζ , a temperature T and a systematic force term F L directed towards the filament's seed. If the friction is low, βζ d 2 min[ U −1 0 , W −1 0 ] (where β = 1/k B T), the fast Brownian motion of the wall leads rapidly to a stationary distribution P(x) ∝ exp (−βF L x) in terms of the tip-wall distance x, between any pair of successive chemical events. The resulting stationary velocity is [1]: where the exponential term gives the probability for an attempted polymerisation step to be successful or equivalently the probability for the wall to create a gap of size d. No closed expression of the velocity v(F L ; ζ ) exists when the wall diffusion gets slower (as ζ gets larger), but it generally leads to a slight decrease of the velocity [8][9][10]. Generalization of this single rigid living filament model towards many filament bundles has been developed over the years [2,3,11], and the exact form of the stationary velocity v(F L ; N f ) for a homogeneous bundle of N f filaments with staggered seeds in the high diffusion limit has been established recently [2,4]. Bundles growing against membranes have also been simulated on the basis of ad hoc stochastic models [12,13] where the obstacle to polymerization is a fluctuating membrane. Filaments are treated as living rigid cylinders, which, in the latter case, have an excluded volume repulsion for the diffusing monomers, which are treated explicitly. Therefore, in the vast majority of applications where an actin bundle is pushing on a mobile obstacle, the filaments are assumed to be rigid.
When the filament's flexibility is taken into account, the situation is much more complex as both the fluctuations of the wall position and the tip bending fluctuations must couple to the chemical steps [3]. Flexible filaments pushing against a mobile wall subject to constant load were investigated recently by dynamic Monte Carlo [8,9]. In this simulation study, the filaments have bending degrees of freedom and can grow or shrink with respective probabilities U 0 and W 0 , the acceptance of a polymerization step being conditioned by a criterion of no overlap with the wall. The wall is undergoing a biased random walk as in the BR model. Emphasis was put on the effects of flexibility and wall diffusivity on the pseudo-stationary velocity v(F L ; L, ζ ) of a single filament in supercritical conditions where polymerization steps dominate. The L dependence of v(F L ; L, ζ ) for cases with finite persistence length L p was not explicitly investigated. The authors stress however the conditions on L and F L to observe the lateral escape of the polymerizing filament: the spectacular fluctuation, which they called pushing catastrophe, becomes possible when a bending fluctuation allows the filament to start growing parallel to the obstacle, preventing further conversion of chemical energy into mechanical work against the load.
To analyze the specific features of flexibility in living biofilaments, it is advantageous to replace the constant load setup by a space-dependent load leading to the establishment of a true equilibrium state. In a series of two papers [5,6], we have established the statistical mechanics framework of a bundle of living filaments in chemical equilibrium with free monomers, i.e., a reactive grand canonical ensemble (GC), within a two-chambered box mimicking an optical trap apparatus (see Figure 1).
The bundle with a fixed number N f of living filaments grafted to the fixed wall of the first chamber presses against a mobile wall subject to an harmonic restoring force −κ T L where the distance L is measured from the grafting wall. The whole box of fixed length L R and fixed transverse area A is coupled to a thermostat at temperature T, and both chambers, separated by the mobile partitioning wall (0 < L < L R ), contain free monomers at the same imposed chemical potential µ 1 > µ 1c where µ 1c is the value of the chemical potential at the critical state defined as the state where the filament in the bulk has no tendency either to grow or to shrink. By tuning the trap strength κ T at fixed (N f , L R , A, µ 1 , T), one can study the equilibrium properties of bundles of filaments for different (self-adjusted) average positions of the wall. The average optical trap length L for bundles of flexible filaments in ideal conditions can be roughly estimated by Hill's polymerization force value F H [14][15][16], the exact average force for rigid filaments [6], so that we expect [5,6]: where:ρ is the free monomer density reduced by its critical value.
In this paper, we exploit a 3D coarse-grained particle-based model (adapted from [17]), and we follow explicitly the 1D dynamics of the wall, of the grafted bundle filaments, consisting of linear assemblies of connected monomers, and of the diffusing free monomers, while the system is subjected to explicit reactive (de)polymerizing events modifying continuously the size of the filaments. To illustrate our simulation approach at fixed temperature T and at fixed chemical potential µ 1 , we treat two cases: a single grafted semi-rigid filament and a N f = 8 filaments bundle, pressing against a mobile wall subject to a harmonic restoring potential. The dynamics of the wall and of all monomers (both free monomers or monomers integrated into filaments) includes the solvent by a Langevin approach with monomer and wall friction coefficients ζ and ζ , respectively. The free monomer chemical potential is maintained by homogeneously adding/removing particles using a Widom insertion scheme. A Monte Carlo local move, with the attempt rate using Poisson reaction time statistics and the acceptance rate based on energetic criteria, deals with reactive monomer assembly/disassembly. The reaction scheme satisfies the detailed balance and leads to chemical equilibrium between the free monomers' concentration and the filament size populations, strongly influenced by the wall pressure on the filament. The average force and the average number of hitting filaments and their fluctuations are computed and analyzed in reference to theoretical predictions proving that our present model reproduces the equilibrium properties of the system well [6]. More importantly, the fluctuations of the wall position and of the filament size show that the system behaves dynamically as a Brownian damped oscillator where the friction term contains both the hydrodynamic component from the solvent and a chemical component related to the filament size changes, with an amplitude [6] of the order of: which corresponds typically to a few monomer sizes d in the actin case. We provide evidence that the chemical friction is originated by the BR mechanism invoked in the theoretical treatment of the system. In this work, we limit ourselves explicitly to rather rigid filaments, therefore limiting the effects of flexibility, with the specific aim to validate the various aspects of the model against known results. The properties of longer and more flexible bundles will be investigated in a future work.
The reactive grand canonical ensemble considers the box of size A × L R limited by two solid walls at x = 0 and at x = L R represented by planes in yellow, with fixed temperature T and fixed free monomer chemical potential µ 1 . The system, periodic in the y and z directions, is divided into two chambers by a mobile wall (depicted as the plane in dark blue) oriented parallel to the fixed walls and located at the variable x = L position (0 ≤ L ≤ L R ). The filament seeds (two first beads in dark blue) are grafted in the wall at x = 0 (on a centered square lattice for the illustrated case N f = 8), and filaments (variable part represented by light blue beads) grow in supercritical conditions within the left chamber at the expense of free monomers (red beads). In the second chamber, only free monomers are present at the imposed chemical potential µ 1 . The mobile wall is loaded by a restoring trap force κ T L indicated by the compressed spring symbol drawn in the right chamber, having a stored mechanical energy V OT = 1 2 κ T L 2 represented in the graph below. Equilibrium results from the cancellation of the sum of the L dependent restoring force and the bundle polymerization force. The molecular dynamics simulation treats explicitly the left chamber system only, the action on the moving wall by the free monomers in the right chamber being replaced by the average force p(µ 1 , T)A parallel to the load force where p is the pressure of a free monomers gas at (µ 1 , T). The wall dynamics and the individual monomer dynamics (both the free ones and those integrated into filaments) are described by Langevin dynamics assuming a solvent producing free draining with friction coefficients ζ for the wall and ζ for monomers. Instantaneous (de)polymerizing reactions follow a local Monte Carlo algorithm integrating or rejecting a free monomer at one filament tip. Free monomer chemical potential µ 1 is maintained in the simulation box by the addition/deletion of free monomers.
The paper is organized as follows. In Section 2, we describe the microscopic model for our reactive mixture of grafted filaments and free monomers. In Section 3, we analyze the various time scales for an actin bundle in the optical trap experiment. We also explain how we set up our model and discuss the choice of parameters that are adopted in the simulation. In Section 4, we present static and dynamic results for a single filament (N f = 1) and for a homogeneous bundle with eight filaments (N f = 8). Our conclusions and perspectives are gathered in the last Section 5. Details of the modeling of the chemical steps and the imposing of the free monomer chemical potential are provided in the Appendixes A and B, respectively.

Microscopic Model of the System in an Optical Trap Setup
Our simulations deal with actin proteins, modeled as spherical units, which can be either free monomers or part of a grafted F-actin filament, a semi-flexible linear assembly of monomers. The system, illustrated in Figure 1, is enclosed in a rectangular box of length L R along the x laboratory axis, having a square transverse area A = L 2 (in the yz−plane). Periodic boundary conditions apply in transverse directions only. The box is limited in the x direction by two fixed walls located at x w = 0 and x w = L R . This box is divided into two chambers of volumes LA for Chamber 1 and (L R − L)A for Chamber 2, separated by an additional mobile transverse partition wall of area A, localized by a variable position L along x (0 < L < L R ), denoted as the mobile wall.
This two-chamber structure allows us to mimic the experimental setup of an actin bundle in an optical trap [18]. We consider an actin bundle with a fixed number N f of filaments, grafted in the fixed wall of Chamber 1 and in contact with a bath of reacting free monomers, pressing against the moving wall (the colloidal particle trapped within an optical trap). In Chamber 2, a bulk solution of free monomers in contact with the same free monomer bath is pressing on the opposite side of the moving wall. All interacting units and the mobile wall (bead) are embedded in a viscous solvent, which will be included in the particle and wall's equations of motion according to the classic Langevin approach. Importantly, the mobile wall is further subject to an external load taken as a harmonic restoring force −κ T L where κ T is the trap force constant. All forces applied to the moving wall appear in a 1D Langevin equation dealing with the wall dynamics, directly coupled to the monomer dynamics. The two chambers are directly connected to a common thermostat (Langevin bath at fixed temperature) and to a common reservoir of free monomers at fixed chemical potential µ 1 , which can exchange free monomers with both Chamber 1 and Chamber 2.
The dynamics of monomers is coupled to instantaneous single monomer (de)polymerization reactions taking place locally at the filament tips, thus only in Chamber 1. In supercritical conditions, the filaments tend to grow and, thus, finally press on the mobile wall. Thanks to the opposite load force increasing with L, a stable equilibrium (including chemical equilibrium) is attained when the wall position oscillates around a stable value: we are then sampling the optical trap ensemble for a bundle of living filaments pressing against a mobile wall, within the framework of the reactive grand canonical ensemble [6].
Although the simulation could include the two coupled subsystems occupying the two chambers, we have restricted our microscopic simulations to the primary box only. The pressure exerted by the free monomers in the secondary chamber is taken into account by adding an average pressure force pA on the mobile wall, where p = p(µ 1 , T) is the pressure of a pure free monomer solution at the imposed chemical potential µ 1 and temperature T. Pressure force fluctuations on the moving wall due to free monomers in the secondary box are thus neglected.
In the following subsections and in the two Appendices, we detail the model of the reactive mixture embedded in Chamber 1.

The Actin Bundle-Free Monomer Reactive Mixture
The N f actin filaments in the reacting chamber are grafted normally to the fixed wall at x w = 0. Free monomers are present in the same volume, and their number fluctuates in time because of chemical reactions with the filaments and because of the exchange of particles with the reservoir at fixed free monomer chemical potential µ 1 . Filaments are thus allowed to grow (shrink), their contour length L c changing by d (monomer size in F-actin) each time a monomer is incorporated into (removed from) the linear structure. They can exert a direct force on the mobile wall when their contour lengths approach the distance ≈ L between the two walls.
Let N t (t) be the total number of monomers in the simulation box at some time t. We will need to trace not only the monomers' positions, but also their insertion/deletion for grand canonical statistics and their chemical character, since they can continuously interchange between free and assembled monomers. A pair of dynamical topological indices (n, k) are assigned to each monomer in the system. The first index n ∈ [0, N f ] indicates the specific filament, n = 0 meaning that the monomer belongs to the solute bath. The second index k ∈ [1, j n ] indicates the rank position of the monomer in the specific filament n. Monomers with indices k = 1 and k = 2 are the first two permanent and fixed monomers forming the seed and used to graft the filaments to the wall. Monomers with k = j n (t) represent the tip monomer of filament n at time t. For n = 0, j 0 (t) represents the number of free monomers in the system at that time; hence N t (t) = ∑ N f n=0 j n (t). Index j 0 (t), hence N t (t), fluctuates because of the coupling to an external solute bath at fixed free monomer chemical potential µ 1 [6]. Between single particle exchange events, N t is constant, but the chemical character and, therefore, the topological indices of the monomers can change because of the occurrence of chemical events. The first two monomers of each filament, (n, 1) and (n, 2), have fixed coordinates (h n , y n , z n ) and (h n + d, y n , z n ), respectively. Any particular filament length j n is restricted to j n ≥ 2. In the special case of a single filament (N f = 1), we set h n = 0 while the y n and z n seed coordinates are irrelevant given the translation symmetry in the transverse plane. For N f > 1, we locate the position of filament seeds in the grafting plane on a planar body-centered square lattice of unit cell size a = L /l where l is an integer. Choosing N f = 2l 2 , the surface density (transverse area per filament) is (L ) 2 /(2l 2 ). The longitudinal disposition of the filament is of primary importance [6]. Setting h n = 0 for all n defines a bundle in the registry (unstaggered). In our application to a multi-filaments bundle, we will adopt the homogeneous bundle case (staggered disposition of seeds) setting, for n even, h n = ( n−0.5 N f − 0.5)d.

The Potential Field
The potential field employed comprises the intrafilament part and the monomer-wall part, since we limit the present study to ideal gas conditions. However a term for monomer-monomer interactions could be included easily.
The intrafilament part for a filament of j n monomers is a sum of (j n − 2) stiff bonding potentials U b (r) forcing the distance r between adjacent monomers (with logical indices (n, k) and (n, k + 1)) to remain close to d (the first grafted bond is rigid with size d). The filament persistence length L p is imposed via a sum of (j n − 2) three body bending terms U bend (θ) implying the bending angle θ between adjacent bonds (i.e., the three monomers with topological indices (n, k), (n, k + 1), (n, k + 2)). Explicitly, we have taken: with κ = 5370k B T and k s = 4000k B T/d 2 . We first note that, when d is set to the monomer size (d = 2.7 nm), the chosen value of κ = L p k B T/d, for k B T = 1, leads to the typical F-actin experimental value L p = 5370 × 2.7 nm = 14.5 µm [19]. As for the bonding spring constant k s , we notice that an F-actin filament of one micron length has a stiffness to longitudinal compression of (44 ± 5) pN/nm (see Footnote b of As the force constant unit is k B T/d 2 = 0.57 pN/nm, one gets k s = (29,000 ± 1500)k B T/d 2 . Our choice to take a value of k s , which is quite a bit lower, is dictated by the need to avoid very fast intramolecular modes, which would require a small time step in the numerical integration of the equations of motion.
Like in a previous study [20], we adopted a purely repulsive 9-3 potential for the monomer-wall interaction: where s = |x − x w | and x w = 0, L corresponds to the position of grafting and the moving wall, respectively, and w = 0.1k B T. The wall potential is truncated beyond its minimum at s c = 3 1/6 d, where it vanishes. The constant 0 = 13.644143k B T in Equation (5) represents the energy gain as a new bond is created within the filament during a polymerization step (and conversely, an energy loss when a bond is removed during a depolymerization step) [5]. As shown in [21], this quantity enters in the size independent bulk equilibrium constant K 0 of the (de)polymerization chemical reaction, established by the chemical equilibrium condition µ i+1 = µ i + µ 1 , where µ i (i ≥ 2) indicates the chemical potential of a filament of i monomers. In ideal conditions, the equilibrium free monomer density ρ 1 and the distribution P i of filaments of size i are linked by [5,6,22]: where the q's are a single grafted filament or free monomer canonical partition functions and V is the system volume. For the intramolecular model defined by Equations (5) and (6), we established in [21] that: Given the adopted values of the parameters, Equation (10) gives K 0 = 39.07144d 3 . At the critical state P i+1 /P i = 1, ∀i ≥ 2, and Equation (9) provides: which is the value of the critical density of free monomers corresponding to a critical chemical potential: where Λ is the de Broglie thermal length of the free monomers. We define a convenient effective chemical potential µ * 1 as: Equation (11) provides: as a critical value of the effective chemical potential.

Monomer and Wall Equations of Motion.
The dynamics of each monomer i (here, i is a short notation for the pair of indices (n, k)) is described by a Langevin equation: where M is the monomer's mass and ζ the solvent friction coefficient with associated random force S i . F intra i is the sum of intrafilament forces imposing linear connectivity and bending flexibility, only present for monomers (n, k) with n > 0. F inter i is the excluded volume (EV) term. Our theoretical developments consider explicitly EV terms even if we have disregarded them in the applications presented in the Result section. Finally, F w i is the interaction of monomer i = (n, k) with the confining walls in thex direction.
The wall dynamics follows the 1D Langevin equation: where M w is the mass of the wall and ζ the solvent friction coefficient on the moving wall with associated random force R. F bun w is the total force on the wall due to the grafted filaments, and F m w is the total force exerted by free monomers. pA is the average pressure term applied to the moving wall due to the free monomers in Chamber 2 of Figure 1 with p = p(µ 1 , T).
The fluctuation dissipation theorem requires that the random forces S i (t) and R(t) satisfy: which sets up the temperature T of the thermal bath. Since we will consider the case of no excluded volume, we have:

Numerical Integration of the Equations of Motion
The monomer and wall equations of motion, Equations (15) and (16), are stochastic second order differential equations. To numerically integrate those equations, we exploited the algorithm proposed by Vanden-Eijnden and Ciccotti [23] (Equation (23) in [23]).
For a one-dimensional system with velocity v and position x, the integrator reads: Here, γ is the friction coefficient divided by the mass, φ n and η n are independent Gaussian variables with zero mean and unitary variance and σ = 2k B Tγ/m. We have adopted a time step of h = 5.33 × 10 −5 τ D required to cope with the vibrational frequencies of our model.

The Free Actin Chemical Potential
The system is surrounded by a reservoir of free monomers at fixed effective chemical potential µ * 1 or equivalently at fixed bulk reduced densityρ 1 : obtained from Equations (13) and (3). We control µ * 1 by adding/removing free monomers along the dynamical trajectory of the system. These addition and removal attempt steps, chosen with equal probability, are performed at randomly-selected Poisson times with adjustable rate ν GC , using a standard algorithm [24]. We used a version with the extra condition that if the added or removed free monomer turns out to be chemically reactive (susceptible to polymerize with a filament tip), the move is automatically refused.
Appendix B provides the justification and detailed procedure. In the present case of a confined ideal system, we get an effective rate of acceptance of the order of 90%.

The (De)Polymerization Steps
Our system is reactive with grafted filaments subject to micro-reversible (de)polymerization steps involving the release or the capture of one free monomer at the filament reactive end. A detailed account of the procedure to perform these (de)polymerization steps and its justification are given in Appendix A.

Numerical Experiments versus Real Experiments
We exploit the, so far unique, experimental results on an actin bundle in an optical trap [18] to estimate all experimental parameters and relevant length and time scales. The complexity of the system under consideration induces a number of widely-separated (many orders of magnitude) characteristic time scales. We need to contract this separation considerably in order to capture all relevant time scales in the same simulation. We argue that this does cause major problems as far as the order of the various characteristic times is respected, and their separation is large enough. We select the set of parameters gathered in Table 1 in order to reproduce, as best as possible, the experimental situation [18]. In the following, we justify our choices.
solvent friction on the obstacle We select the monomer size d, the thermal energy k B T and the wall diffusion time: the time needed for the colloidal bead with diffusion constant D to diffuse in pure water solvent over a distance d, as units of length, energy and time, respectively. The monomer size in the F-actin filament is d = 2.7 nm; the thermal energy at T = 300 K is k B T = 4.142 × 10 −21 J. The unit of time τ D for the experimental system can be obtained by estimating the diffusion through the friction experienced by a colloidal bead of radius R = 1µm in pure water (η 0.001 Pa s). By the Stokes law, the friction is Table 1.
The persistence length L p /d in the model is fixed to a typical experimental value L p = 14.5µm. Another relevant length scale is the typical absolute length of (un-crosslinked) actin filaments pressing on the wall. The equilibrium value L /d is roughly given by Equation (2) in terms of the external parameters N f ,ρ 1 (or equivalently µ * 1 ) and κ T [6]. To remain in the non-escaping regime where the obstacle can effectively stop polymerization, L /d must remain in the range 20-100 for the range ofρ 1 explored [6]. Hence, all external parameters are adjusted in the simulations to cope with the experimental situation. We thus note that all length scales probed in our mesoscopic simulations are representative of corresponding experimental quantities. K 0 d −3 is the bulk equilibrium constant of the (de)polymerization reaction, which, using Equation (11), fixes the free monomers' critical concentration ρ 1c . The physically relevant quantity is, however,ρ 1 = ρ 1 /ρ 1c , the reduced free monomer density, which lies in the range 1.5-4.0 in most in vitro experiments. This quantity is directly linked to µ * 1 using Equation (3), that is the chemical free energy, which can be converted into useful mechanical work by compressing the trap. The optical trap strength used in experiments is of the order of κ T ≈ 0.008 pN/nm [18]. As shown in [6], in the present range of values ofρ 1 , we need to use trap strengths per filament in the range κ T /N f ≥ 0.017 pN/nm to avoid the escaping filament regime, of the same order of magnitude as in the experiments. Energy scales, length scales and, thus, typical trap and polymerization forces are realistic in our simulations.
We now investigate the much more complex situation with time scales. All of the relevant ones are gathered in Table 2 and discussed below as relative quantities with respect to τ D defined earlier in Equation (24).

Experimental vs. Model Time Scales
In this mesoscopic system, there are a number of characteristic time scales that need to be considered for realistic modeling. From fast to slow modes, they are: the intrafilament modes, bond stretching, τ s , and bond bending, τ b ; the free monomer characteristic times, inertial and diffusion times, τ in fm and τ fm , respectively; the obstacle characteristic times: inertial, oscillatory (due to the trap) and diffusion times, τ in W , τ T and τ D , respectively, of the latex bead in the experiment (of the planar wall in our model); the characteristic time between two successive chemical events at the tip of the filaments τ chem in the bulk.
An estimate of the intramolecular times for F-actin can be obtained by assuming a bead-spring model similar to our present one. The characteristic times of the stretching and the bending modes are respectively τ s = 2π √ M/k s and τ b = 2π √ Md 2 /2κ where M is the mass of the bead (For the bending mode, we consider the single triatomic molecule with fixed central atom and fixed energy . For actin M = 42 kDa = 6.977 × 10 −23 kg, hence τ s = 12.9 ps and τ b = 21.2 ps using the values of k s = 29,000 k B T/d 2 and κ = 5370k B T (see Section 2.2). In terms of the obstacle diffusion time reported in Table 1 (our unit of time), we have τ s = 3.69 × 10 −7 τ D and τ b = 6.06 × 10 −7 τ D .
In our simulation, we fixed M = 0.003556(τ 2 D k B T/d 2 ), k s = 4000k B T/d 2 and κ = 5370k B T and obtain τ s = 5.9 × 10 −3 τ D and τ b = 3.6 × 10 −3 τ D . Our choice of parameters provides intramolecular times, relative to τ D , roughly four orders of magnitude larger than in the experimental system. However, we believe that the remaining gap of three orders of magnitude is sufficient to decouple intramolecular modes from the slower modes in the system, still providing an affordable working framework. Furthermore, our choice of k s leads to contour length fluctuations about three times too large. Indeed, for a filament of N 50 monomers, their amplitude is ≈ √ k B TN/k s = 0.1d (assuming N springs of stiffness k s in series and assuming that the average vibrational potential energy is (1/2)k B T). Such small fluctuations should affect only marginally the rate of chemical (de)polymerization controlled by the wall position and the filament bending fluctuations.
The next characteristic times concern the free monomers, namely G-actin proteins in water solvent. The typical radius of a G-actin protein, considered as a spherical particle, is r fm = 2.9 nm = 1.074d [25]. From the Stokes law, we get ζ = 6πηr fm 5.5 × 10 −11 J s/m 2 . The inertial time of free monomer is thus τ in fm = M/ζ 1.3 ps = 3.6 × 10 −8 τ D and the diffusion time τ fm = r 2 fm /D fm = r 2 fm ζ/k B T = 1.1 × 10 −7 s = 3.14 × 10 −3 τ D , five orders of magnitude larger than τ in fm . In our simulation r fm = d, and we have chosen ζ = 0.5τ D k B T/d 2 , so that τ in fm = M/ζ = 7.1 × 10 −3 τ D and τ fm = 0.58τ D . Here, the separation between the inertial and the diffusion time scales is only three orders of magnitude. More importantly, the relative diffusion time with respect to the obstacle diffusion time τ D is quite a bit larger than in reality, meaning that free monomers diffuse much too slowly in our model. This could influence the rate of chemical events, hence the rate of supercriticality of the solution, since a chemical reaction needs available free monomers in the reacting region to occur. We obviate this deficiency by operating a continuous removal/addition of free monomers in the system in order to simulate the grand canonical ensemble at fixed µ * 1 . We have adopted a rather high insertion/removal attempt frequency ν GC = 187.5τ −1 D (see Section 2.5) to ensure a free monomer density in the system as uniform as possible. Since the total number of free monomers in the system N 1 is of the order of 100 (see Table 3) and since the degree of acceptance of these moves is of the order of 90%, during the characteristic time τ D , each free monomer is inserted and removed once on average.
The next characteristic times concern the moving obstacle. In the absence of the F-actin bundle, the obstacle is subject to the harmonic potential from the optical trap and to the friction from the water solvent. For a harmonic oscillator with friction, the overdamped condition (real solutions of the characteristic equation) reads [26]: where M w is the mass of the moving object and τ T = 2π √ M w /κ T is the period of the obstacle oscillations in the optical trap. Assuming a mass M w = 4π 3 R 3 ρ ≈ 4.2 × 10 −15 kg for a latex bead of radius 1 µm (we have assumed the density of latex roughly equal to the density of water) and the typical trap strength κ T = 0.008 pN/nm [18], the experimental values are τ in w = 2.23 × 10 −7 s = 6.3 × 10 −3 τ D and τ T = 1.44 × 10 −4 s = 4.1τ D .
In our model, we have chosen M w = 0.010667(τ 2 D k B T/d 2 ) and ζ = τ D k B T/d 2 = 1; hence, τ in w 10 −2 τ D , τ T = 4.7τ D and τ T = 1.2τ D for the adopted κ T values, respectively 0.019375k B T/d 2 and 0.275k B T/d 2 for the N f = 1 and the N f = 8 cases.
This analysis shows that, both in the experimental case and in our modeling, the obstacle, even in absence of the bundle, undergoes an overdamped motion. This conclusion is in agreement with the transient behaviors observed in Reference [18]. Moreover, our model reproduces correctly the ratio between the free oscillation time of the obstacle in the trap and obstacle diffusion times, while the inertial time is somewhat slower than in the experiments, but still two orders of magnitude faster than the diffusion time.
The last important characteristic time is related to chemical reactions at the tips of the filaments. We introduce a chemistry characteristic time τ chem = W −1 0 for which the standard experimental value W 0 = 1.4 s −1 leads to τ chem /τ D = 2.04 × 10 4 . This indicates that the wall diffusion is much faster than the time scale over which an F-actin filament fluctuates in size (Note that the chemical reaction is a very fast event, usually considered instantaneous. Here, τ chem gives an estimate of the time interval between two successive reactions). It is useful to mention here that these conditions are precisely assumed to establish the rigid filament velocity-load relation in Equation (1) and its many filament generalization provided by the Mogilner-Oster model ( [2][3][4]).
As we will show in the next section, the obstacle dynamics in the presence of the F-actin bundle undergoes a motion at equilibrium with a characteristic time, which we denote τ L , roughly 1000 times larger than τ chem , the largest time discussed so far. This is the signature of the Brownian ratchet mechanism governing the interaction between the chemically-fluctuating bundle and the loaded obstacle. In order to get quantitative information on the obstacle motion, we need to follow the dynamics of the system for many τ L . With the adopted time step h, τ L corresponds to roughly half a billion time steps; therefore, any trajectory should contain several billions of time steps in order to compute τ L . This is the reason why we had to compress considerably the experimentally wide range of characteristic times as discussed in this section. Figure 2 illustrates the occurrence of several, widely-separated time scales in the system dynamics for the single filament case discussed in the Result section.

Illustrative Results
We report results from two simulations, one for a single filament (N f = 1) and one for a homogeneous bundle of N f = 8 filaments, the latter case being the typical system observed in the in vitro experiments of Reference [18]. Except for the values of N f and κ T , all experiments are done with the same external parameters, i.e., a transverse area A = 36d 2 , a temperature k B T = 1 and the free monomer chemical potential βµ * 1 = −3.6654 + ln 2.5 = −2.7492, corresponding toρ 1 = 2.5. The pressure force due to the free monomers in the second chamber of the optical trap ensemble in Equation (16) is fixed by Equation (19). All internal parameters, described in Section 3 are taken to be similar in the two experiments.
The choice of κ T requires some care, and we exploit for guidance the theoretical predictions for rigid filaments, Equations (2) and (4). The non-escaping regime for actin filaments atρ 1 = 2.5 is limited to an extreme value L max = 70.2d [6]. By choosing κ T = 0.019375kT/d 2 for N f = 1, we expect to be in the condition L = 47.3d < L max − 3 k B T κ T = 70.2d − 21.6d = 48.6d, which guarantees that the whole range of L values probed during the wall fluctuations lies within the non-escaping regime [6]. For N f = 8, we take κ T = 0.275k B T/d 2 , which provides L ≈ 26.6d with σ L = 1.90d for rigid filaments (see the Supplemental Material of Reference [6]), again in the good range of L values.
As discussed quantitatively in Reference [6], short flexible filaments behave like rigid filaments, while deviations are observed for longer filaments even before entering in the escaping regime. Further dynamical effects of filament flexibility will be discuss in a future publication [7]. We limit the present investigation to nearly rigid cases since our main aim here is to validate the simulation approach against known theories and/or existing experiments.
To generate the initial configurations for our optical-trap simulations, we have first fixed the obstacle wall at the expected average value given by Equation (2) and let the filaments, with initially only two monomers, grow during the dynamics. After a long enough simulation, the system reaches the statistical equilibrium at given (µ * 1 , T, L, A). This state is the initial configuration for a new run with the mobile wall according to the equations of motion (16). The wall position is monitored in time to let the wall equilibrate before starting the production run. This procedure is repeated in parallel for a number of equivalent initial configurations (with only two monomers per filament) with a different random number seed, which ensures that we obtain statistically independent, but equivalent realizations of the dynamics, for the purpose of statistics.
For N f = 1, we ran 32 independent trajectories for a time interval of 5.6 × 10 4 τ D per trajectory. The observed wall relaxation time, obtained from the time correlation function of the wall fluctuation, is τ L 950τ D ; hence, each trajectory covers more than 50 wall relaxation times.
For N f = 8, we ran 50 independent trajectories of 1.3 × 10 4 τ D . In this case, the observed wall relaxation time is τ L 1400τ D ; hence, each trajectory covers roughly 10 wall relaxation times. Figure 2 shows a portion of an equilibrium trajectory for the N f = 1 case, in a time window of τ L ≈ 900τ D , roughly one wall relaxation time. Chemical reactions are easily detected by the discontinuous variations of the number of bonds k(t) in the filament. The normal component of the filament end-to-end vector X(t) and the instantaneous wall position L(t) show that the wall executes a Brownian ratchet type dynamics, polymerization steps being clearly allowed by natural excursions of the wall towards larger L values. Note that in general, X(t)/d is slightly smaller than k due to the filament bending. Occasionally, the opposite is seen when the filament is relatively straight and subject to a positive fluctuation of its contour length due to bond stretching.

Static Properties
In Table 3, we report relevant equilibrium averages of the static properties for the bundles. We compute the average trap length L and its distribution, the average size of the bundle I = ∑ N f n=1 j n /N f (a collective property), the average size of the single filament in the bundle i = j n (an individual filament property), the average longitudinal component of the end-to-end vector of the single filaments X = x j n − x 1 and the average modulus of the transverse component of the end-to-end vector of the single filaments R ⊥ = (y j n − y 1 ) 2 + (z j n − z 1 ) 2 . Table 3. Equilibrium averages of static properties of the systems. N 1 indicates the number of free monomers in the system, L the average position of the wall, L H the wall position according to Equation (2), I the average size of the bundle and i the average size of a single filament and X and R ⊥ the longitudinal and transverse components of the end-to-end vector of the single filaments, respectively.   [6]. We recently demonstrated the statistical mechanics foundations of this Gaussian distribution, with maximum at L H = F H /κ T and width σ H = √ k B T/κ T . In particular, F H can be shown to be the average force exerted by the bundle of rigid living filaments on a hard mobile wall, within the optical trap grand canonical ensemble and in the continuous limit. The present filament model has bending flexibility adjusted to F-actin, small contour length compressibility (with amplitude of ≈ 0.1d) and a discretized contour length with step ≈ d. Moreover, the filament-wall interaction is a soft repulsive potential Equation (8) with a wall effective core at s * = 0.846d defined by U w (s * ) = k B T. In Reference [6], we considered a model of a hard obstacle and discussed the behaviors of a fully-rigid and a flexible (only bending flexibility adjusted to F-actin) model, with a discrete contour length step of d. For a single rigid filament, spectacular oscillations on the spacial scales d in P(L) were observed as a result of discreteness effects (see Figure 3 of Reference [6]). For the flexible model, they were rounded off to some extent by the flexibility. The features of the distribution in Figure 3 for N f = 1, for the present model, show a similar behavior, with an expected overall Gaussian shape and oscillations on the scale of d, rounded off by the filament flexibility, filament compressibility and softness of the wall repulsion. In Reference [6] the flexibility was shown to produce an increase of the average force by a few percent, with respect to the rigid case. Our present model provides an average wall position (average force) in agreement with the rigid model results within the present statistical uncertainty. This shows the adequacy of the present model, at least for static equilibrium properties.
The N f = 8 case in the right panel of Figure 3 shows a much smoother behavior, as expected for a homogeneous bundle. Results for the average trap length and its variance are reported in Table 3. We note that the present bundle model provides a slightly larger average trap length ( L > L H ) and fluctuation. A similar behavior was observed in Reference [6] for the model with flexible, but incompressible filaments in the presence of a hard wall (see Figure 5a in that reference for a similar case and its Supplemental Material for numerical values of the system with the same external constraints). This observation reinforces our previous finding that flexibility enhances the average bundle force. Furthermore, taking as a reference Hill's value for the average trap length L H , the deviation in the present flexible and compressible filament model in the presence of a soft wall ( L /L H − 1 = 0.015) is roughly two-times the deviation of the flexible, but incompressible model with the hard wall ( L /L H − 1 = 0.006; see the Supplemental Material of Reference [6]), suggesting that the compressibility and the soft character of the obstacle further enhances the effects of flexibility.
Direct calculations of the average force exerted by the filaments on the wall give F bun = 7.4413k B T/d, F 2 bun = 153.95(k B T/d) 2 , in very good agreement with the estimate from κ T L = 7.437(k B T/d), the equivalence being easily proven by statistical mechanics [6]. In Reference [6] for a discrete worm-like chain (WLC) with seeds at h n , hitting a hard wall located at L, we introduced a crossover index z n by the expression: The filament n with size z n (and contour length L c = (z n − 1)d) is the longest filament starting at h n , which does not hit the wall at L. The filament n of size i n can then also be specified by a relative index: where m > 0 indicates a filament hitting a wall and m ≤ 0 a filament that does not interact directly with the wall.
In the present study, we deal with a slightly different model, as we have a soft repulsive wall that starts at a distance s c ≈ 1.20d from the geometrical position of the wall at L, and in addition, the contour length is no longer strictly constant. We adopt here the same definitions of z n and m in terms of h n and the geometrical wall position L, keeping in mind that a filament will now hit the wall strongly when i n ≥ z n (m ≥ 0). For the case i n ≤ z n − 2 (m ≤ −2), the filament (with a fixed contour length of at most (z n − 2)d) does not interact with the wall, as its tip remains further than 0.8d from the plane at L − h n − y c where the wall repulsion would start to be felt. The small filament contour fluctuations (of the order of 0.1d) would not substantially change the situation. Finally, the case m = −1 is clearly a case of weak interferences with the wall. Figure 4 reports the relative size distribution Q m for the two experiments. We observe a exponential rise Q m ∝ρ m 1 (to a very good accuracy) for negative m values up to m = −1 inclusive. The decay of Q m for m ≥ 0 reflects a strong wall influence and, hence, the decrease towards zero of the wall factor α(i n , L − h n ) < 1 introduced in References [5,6], and giving Q m ∝ α(z n + m, L − h n )ρ m 1 . The vertical shift between the two exponential rises reflects the strong difference in the tails for positive m, which affects the normalization (for N f = 8, Q 1 is out of scale and not significantly different from zero).

Analysis of Chemical Kinetics
During the experiments, we monitor different counters in which we cumulate information for the N f filaments at each time step h. The m index of each individual filament is computed, and a counter aiming at the Q m calculation is increased by one. If a (de)polymerization step is attempted at that time step for a filament of relative size m, we add one to a counter of (de)polymerization attempts for a filament index m. If a (de)polymerization step is accepted and realized at that time step for a filament of relative size m, we add one to a counter of (de)polymerization successful steps for a filament of index m. At the end of the run, the chemical events' counters with index m are divided by the Q m counter and by the time step h to get the best estimates of the attempt (de)polymerization rates U att m and W att m and of the effective (de)polymerizing rates U m and W m . At equilibrium, the total rate of polymerization must be equal to the total rate of depolymerization, (28) Figure 5 shows the effective (de)polymerization rates W m , U m computed in both experiments (N f = 1, 8). Agreement is observed between the two experiments. Furthermore, the ratio of the two plateau values, corresponding to the bulk (de)polymerization rates, provides the value of 2.54, in tight agreement with the supercritical condition: imposed by fixing the free monomer chemical potential µ * 1 with the external reservoir. Attempt rates are given for the N f = 1 case in the Appendix (see the Figure in

Wall Relaxation and Related Equilibrium Time Correlation Functions
In Figure 6 for the bundle of N f = 8 filaments, we plot the time correlation function of the equilibrium fluctuations of several related properties. The autocorrelation function of wall fluctuations C LL (t) = δL(t)δL(0) exhibits a single exponential behavior well represented by C LL (t) = 3.785 exp(−t/τ L ) with τ L = 1375τ D . A similar single exponential behavior (with a slightly larger zero time value; see Table 3) is followed by the autocorrelation function of the fluctuations of the bundle size C I I (t) proving that the slow relaxation of the wall fluctuation is tightly related to the chemical events and to the Brownian ratchet mechanism induced collectively by the bundle. On the other hand, fluctuations of single filament properties, like the longitudinal component of the single filament end-to-end vector C XX (t) and the single filament size C ii (t), exhibit a relaxation with two characteristic times that can be well fitted by a linear combination of exponentials. The slow relaxation time is ≈ τ L , while the fast relaxation time is τ fast 30τ D .
For the single filament case, we observe a single exponential relaxation (not shown) with a characteristic time of τ L 950τ D .
Finally, in Figure 7, we report for both N f = 1 and N f = 8 the autocorrelation function of the fluctuations of the length of the transverse single filament end-to-end vector. Again, a linear combinations of two exponential decays represent well the relaxations, the long relaxation time being again the trap length relaxation time for the respective cases. Quite different values are observed for the two short relaxation times; for N f = 1, we have τ fast 25τ D , while for N f = 8, it is τ fast 2.7τ D . This short characteristic time is related to the relaxation of the intramolecular degrees of freedom, and the observed difference between the two cases can be qualitatively justified observing that the eight-filament bundle system has an average trap length quite shorter than the one-filament system, with shorter and more rigid filaments, hence with faster relaxation times. Note that a wide separation between filament equilibration times and the wall characteristic time is at the basis of the adiabatic approximation used in [7] for flexible filaments.  Figure 7. Time correlation function of the fluctuation δR ⊥ (t) = R ⊥ (t) − R ⊥ of the modulus R ⊥ = (y j n − y 1 ) 2 + (z j n − z 1 ) 2 of the transverse part of the end-to-end vector of the single filament of the bundle with N f = 8 (red circles) and N f = 1 (blue diamonds).

Conclusions
In this work, we have presented a particle-based model to simulate the dynamics of a grafted bundle of parallel living actin filaments pushing on a mobile wall subject to load. The underlying statistical mechanics framework is the reactive grand canonical ensemble dealing with a two-chambered system with a mobile and loaded partition wall separating a supercritical mixture of G-actin monomers and living flexible F-actin filaments from a pure bulk solution of G-actin. Our goal is to provide a working model to study, e.g., the growth of F-actin filopodia against a resisting membrane.
Here, we specialize the model to represent an optical trap setup, a growing bundle of semi-flexible actin filaments pressing on a wall subject to a restoring force increasing linearly with the wall position L. This apparatus allows in particular to measure the bundle polymerizing force [18] (see Figure 1). Our work provides a theoretical tool able to predict quantitatively the link between the wall motion and the bundle's filament kinetics, the connection empirically described by the velocity-load relationship. In addition, we can also analyze the critical conditions under which the occurrence of filaments escaping laterally, hence no longer participating in the conversion of chemical free energy into work on the loaded wall, can be safely avoided. Our simulation approach is then complementary to optical trap experiments on the actin bundle polymerization forces, experiments that have proven so far to be difficult to perform and to interpret, mainly as a result of the flexible character of the pushing filaments.
Our model actin proteins are diffusing (mesoscopic) point-particles able to assemble/disassemble into linear semi-flexible structures by explicit single monomer (de)polymerization steps, under conditions where the free monomer chemical potential is imposed uniformly in the system. The chemical steps at the filament tips are governed by stochastic rules which fix the bulk (de)polymerizing rates and automatically deal with a realistic influence of the obstacle wall on the chemical rates when the tip is approaching the wall. Although the method works if excluded volume (EV) interactions are considered, the present work is limited to ideal system conditions (no EV) so that monomers interact only with the walls and, when integrated within a linear filament, are subjected to bonding and bending interactions with first and second neighboring monomers.
We have shown that the present mesoscopic approach with explicit consideration of filament bending flexibility allows us to tackle all relevant length scales realistically, in particular the filament length and the filament persistence length. Conversely, given the wide spread of the realistic time scales from the very fast filament bending modes to the slow wall relaxation kinetics induced by the Brownian ratchet mechanism of the pushing bundle, one has to condense the times scales on a narrower range. We have accelerated artificially the chemical reactions by a factor of ≈ 1700 to make our calculations feasible within a single simulation approach. Despite such an acceleration of one particular process, we have preserved the hierarchy of characteristic times associated with the various dynamical processes taking place within the system.
We illustrated the approach considering a single filament system and a system with a homogeneous bundle of eight filaments. We studied static properties and fluctuation dynamics at equilibrium at the super critical conditionρ 1 = 2.5. The trap strength κ T was chosen in each experiment, such that important flexibility effects are avoided, the strategy being to first test the method on cases where a rigid filament approach provides a semi-quantitative, but robust reference.
We have shown that our results are coherent with those predicted theoretically both for a strictly rigid model and for a model of flexible worm-like chains hitting a hard wall in the grand canonical ensemble [6].
The present approach needs further tests in conditions where flexibility effects become more relevant (lower κ T , longer filaments) and in conditions where the acceleration of the chemistry is reduced to test the influence of the ratio of characteristic times governing chemical steps and wall diffusion in the absence of the bundle. More ambitious developments could envisage the replacement of the rigid obstacle wall by a flexible membrane, the consideration of additional proteins developing branched networks or capping active actin filaments to stop polymerization. Consideration of free monomer concentration inhomogeneities and the distinction between ATP and ADP actin complexes could also be envisaged, showing the rich potential of our model.
• Ω s (n, j n ) is a spherical layer centered at the location of the tip monomer (n, j n ) with volume: having interior and exterior radii R min = 0.8909d and R max = 2d, respectively (These bounds were kept unchanged with respect to [17], where they were motivated by the presence of EV interactions modeled as shifted Lennard-Jones potentials. In the present case where EV are not included, these parameters are merely an arbitrary choice for defining the tip's volume of capture for diffusing free monomers). Any free monomer located within Ω s (n, j n ) at time t is defined as a reactive free monomer with respect to filament n at time t. Note that several distinct free monomers can be reactive with the same filament's tip and that a single free monomer can be reactive with respect to the tips of several filaments.
• Ω c (n, j n ) is centered on the point where, for a given microscopic configuration of the j n − 1 first monomers, both the vibrational part of the bonding potential and the bending vibrational potential energies involving the tip monomer (n, j n ) are zero. This point is located at a distance r = d from monomer (n, j n − 1) and is collinear with the bond axis u connecting monomers (n, j n − 2) and (n, j n − 1). The portion of space Ω c (n, j n ) is delimited in space by a spherical shell centered on monomer (n, j n − 1) with extreme radii r min and r max and by a cone with apex at the same monomer (n, j n − 1) location, with symmetry axis u and with an angle opening θ max , where all of these r and θ extreme values correspond to the distances or angles for which the purely intramolecular vibrational term in Equation (5) or in Equation (6) is equal to 3k B T. This portion of space has a volume of: A tip monomer (n, j n ) will be defined as reactive at time t if it lies within Ω c (n, j n ), which will be the case most of the time along equilibrium trajectories.
We request that the reactive steps satisfy the micro-reversibility condition at equilibrium: where J and I are the two microscopic states separated by the reversible single step (de)polymerization reaction, π(J) ∝ e −βu(J) is the equilibrium Boltzmann weight of state J and P(J → I) represents the conditional probability to go from state J to state I. As in the usual Metropolis method, the conditional probability is split into the product of two probabilities P(J → I) = T(J → I)A(J → I), an a priori sampling probability T(J → I) and an acceptance probability A(J → I). The two distinct states are more precisely defined by the following two processes: • The polymerization step: Consider the filament n with j n monomers in a given microscopic configuration and a free monomer located at some specific position r s ∈ Ω s (n, j n ), and consider the process of sampling new coordinates for this free monomer r c ∈ Ω c (n, j n + 1) and adding the necessary intramolecular interactions in order to make this monomer the new tip of an extended filament of size j n + 1. We denote J the initial microscopic state and I the final state, which are illustrated respectively in Figures A1 and A2; • The depolymerization step: Consider the same filament n with j n + 1 monomer in the same microscopic configuration (state I), and consider the process of moving the tip monomer at its original position prior to the polymerization step r s ∈ Ω s (n, j n ) and removing the intrafilament interactions involving monomer j n + 1. The microscopic state has changed from I into J in this case.
R max min R r max rmin Figure A1.
State J: The depolymerized state for filament n with size j n where the reacting monomer of index k shown in blue is free, but is necessarily located within the spherical layer Ω s (n, j n ) limited by spheres of radii R min and R max centered on monomer (n, j n ) shown with dashed-blue lines (see Equations (A1) and the related text for details). The dashed sphere represents the position of the reacting monomer after a possible polymerization ending with configuration in state I (see Figure A2). Green filled circles represent other free monomers in the environment, while connected full circles represent the grafted filament end monomers in the depolymerized case. Reproduced Figure A2. State I: The polymerized state for filament n with size (j n + 1) where the reacting monomer with intrinsic index k shown in dark blue is the tip of the filament and is, necessarily, located in a low intramolecular energy portion of space Ω c (n, j n + 1) resulting from the intersection of a conic volume (shown in green) of angle θ max and a spherical layer (shown in red) limited by spheres of radii r min and r max (see Equation (A2) and the related text for details). The dashed sphere represents the position of the reacting monomer after of a possible depolymerization ending with configuration in state J (see Figure A1). Reproduced from J. Chem. Phys. 136, 114901 (2012) [17]; http://dx.doi.org/10.1063/1.3694672, with the permission of AIP Publishing.
Given the reactive state J (I), the a priori probability of state I (J) is the product of the probability to attempt the reaction and the probability to generate the new position in the final state. We have defined the attempt rates ν pol and ν dep as: where each rate is proportional to the volume of the region where the new coordinates will be located and ν is a rescaling free parameter with the dimension of inverse time (related to an adjustable energy barrier of the (de)polymerization reaction). In Equation (A5), 0 is the bonding energy term of Equation (5), which acts as a penalty for the depolymerization rate (it is related to the difference between the two energy minima at the two sides of the barrier).
For each microscopic configuration of the system generated during the dynamics at discrete time t l = l h, all reactive states are systematically detected. We assume that the times at which the possible reactions occur are independent of each other and distributed according to a Poisson law P α (t) = ν α e −ν α t (α = pol, dep). The probability to attempt a reaction during the next time step h is the cumulative function of P α (t), namely F α (h) = (1 − e −ν α h ). To decide whether to attempt a given reaction, we extract a uniformly-distributed random number and compare it with F α (h) of that specific reaction. Given the values of h and of the rates ν α , in practice, F α (h) ≈ hν α 1, so that in the majority of cases, no reaction is attempted for all reactive monomers, and the dynamical integration algorithm is applied once more to the unaffected microscopic configuration. From time to time, one event is selected, and that unique reaction is attempted by sampling the new location of the reacting monomer uniformly within the target relevant reacting portion of space, defined by the nature of the transition. Therefore, the a priori conditional probabilities are: As in the usual Metropolis algorithm, the detailed balance condition Equation (A3) is satisfied by [27]: which is the usual condition on the energy difference except for the term 0 . Note that it appears in both equations with the same sign as u(I). The state I differs from J because it has one more bond with associated bonding energy Equation (5), where a − 0 term cancels the analogous term in both Equations (A8) and (A9) and associated bending energy Equation (6). Moreover, the values of non-bonding interactions (interaction with the wall and excluded volume interactions when present) involving the displaced monomer are different between the two states. Note that a different, but equally correct algorithm could avoid the penalty term exp(−β 0 ) in the depolymerization attempt rate and let it in the final acceptance probability. In this case, however, depolymerization events would be sampled much more frequently, but also rejected more frequently with a loss of efficiency of the procedure. During a chemical step, we do not advance the time of the dynamics, no matter if the reaction is accepted or rejected. Therefore, we treat reactions as instantaneous processes occurring always at the beginning of the generic time step. This is a possible way of coupling normal dynamics with MC steps. A possible justification is that during a chemical event, at most one particle is displaced while the others are still.

Appendix A.2. Prediction of Macroscopic Rates for the Ideal System
It is instructive to consider the ideal mixture treatment of the (filaments + free monomers) reactive system with the aim to provide computable expressions for the rates U m and W m introduced in Section 4.2.1. For the sake of simplicity in the notation, we limit our derivation to the single filament case; the extension to the multi-filament system is straightforward for non-interacting filaments. Let us consider the filament in state J with j monomers corresponding to m = j − z (here, h n = 0) that polymerize to the state I with j + 1 monomers, hence to m + 1 (moving its tip closer to the wall). Because the filaments are rather rigid and the obstacle is a plane, these processes take place most probably in a specific slab parallel to the wall. The rate at which such a process progresses can be written as the product of three different terms: (1) the number of free monomers available for the reaction N r (m); (2) the rate of the single polymerization event ν pol given in Equation (A4); (3) a global acceptance for the proposed processes A p (m) defined below. The number of reactive monomers N r (m) in the slab of index m is the integral over Ω s (m) of the free monomer density profile ρ 1 (r). Near the wall (m ≥ −2), the density profile is modified by the repulsive interaction with the wall, and for an ideal system (no particle correlations), it takes the form ρ 1 (r) = ρ 1 exp[−βU w (r)] where U w is given in Equation (8). Therefore, we can write: A p (m) = Ω s (m) dr J ρ 1 (r J ) Ω c (m+1) dr I Min 1; e −β[ vib (r I )+U w (r I )−U w (r J )] V c Ω s (m) dr ρ 1 (r) where r J and r I represent the position of the reactive monomer in states J and I, respectively, and vib (r I ) represents the vibrational contribution to the energy from the additional bond in state I. A similar expression can be cast for the depolymerization rate W m , where configurations in state I with j monomers, corresponding to m = j − z, go to state J with j − 1 monomers, corresponding to m − 1. In this case, however, the number of reactive monomers must be replaced by the number of reactive tips N t (m), i.e., the fraction of tip states with large enough statistical weight according to our definition above: where V is the volume of the entire system, and the wall interaction term ensures that tip positions are confined in the correct region of space. In analogy with the polymerization process above, we can write: Those expressions becomes particularly illuminating when m ≤ −3, i.e., far from the confining wall, and the wall potential vanishes. In that case, both U m and W m become independent of m and equal to their corresponding bulk values U 0 and W 0 , respectively. Far from the wall, the free monomer density profiles become constant and equal to ρ 1 and N r (m) = V s ρ 1 . The integral over r J in the numerator of Equation (A12) cancels the one in the denominator. Moreover, the Min[. . . ] operator in Equation (A12) always selects the exponential since the vibrational energy, the only surviving term, is always positive. As for the depolymerization rate, N t (m) still presents the ratio of two integrals, while the expression of A d (m) Equation (A16) becomes trivial since now Min[. . . ] = 1, again, because of the sign of the vibrational energy. We have: where the identity K 0 = e β 0 V dr exp (−β vib (I)), following from the definition of K 0 in Equation (9), was derived in Equation (27) of [17]. We note that V s V c , and the two integrals in Equation (A18) are roughly equal, since Ω c is the portion of space contributing most to the integral in the denominator for the added bond; we can estimate the bulk (de)polymerization rates as: We have computed numerically U att m , W att m and U m , W m using Equations (A11), (A15), (A10) and (A14) during our calculations, averaging also over the specific position of the reactive tip in state J at fixed m and over the wall fluctuations in our optical trap experiment. Results for the single filament system are shown in Figure A3 for both the attempt rates and the effective rates (already shown in Figure 5). We observe as U att m is roughly constant and equal to 1.2 far from the wall (m ≤ −3), a value in agreement with V s V c /(V s + V c )ρ 1 ν as indicated by Equation (A17). Closer to the wall, U att m decreases roughly linearly between m = −2 and m = 0, where it reaches a new constant value around 0.75. This decrease is therefore the combined effect of reducing the available volume of the spherical shell of the tip when approaching the wall and also of reducing the free monomer density because of the repulsion with the wall. As better seen in Figure 5, the final polymerization rate for m ≥ 0 vanishes because no more polymerizations are accepted. Conversely, W m and W att m are always very close to each other.  Figure A3. Relative size polymerization and depolymerization rates for the N f = 1 case at ρ 1 = 2.5. The comparison between the attempt rates U att m (blue triangles) and W att m (gray squares) and the effective rates U m (red circles) and W m (superimposed to W att m ).