Structural organization and energy storage in crosslinked actin assemblies

During clathrin-mediated endocytosis in yeast cells, short actin filaments (< 200nm) and crosslinking protein fimbrin assemble to drive the internalization of the plasma membrane. However, the organization of the actin meshwork during endocytosis remains largely unknown. In addition, only a small fraction of the force necessary to elongate and pinch off vesicles can be accounted for by actin polymerization alone. In this paper, we used mathematical modeling to study the self-organization of rigid actin filaments in the presence of elastic crosslinkers in conditions relevant to endocytosis. We found that actin filaments condense into either a disordered meshwork or an ordered bundle depending on filament length and the mechanical and kinetic properties of the crosslinkers. Our simulations also demonstrated that these nanometer-scale actin structures can store a large amount of elastic energy within the crosslinkers (up to 10kBT per crosslinker). This conversion of binding energy into elastic energy is the consequence of geometric constraints created by the helical pitch of the actin filaments, which results in frustrated configurations of crosslinkers attached to filaments. We propose that this stored elastic energy can be used at a later time in the endocytic process. As a proof of principle, we presented a simple mechanism for sustained torque production by ordered detachment of crosslinkers from a pair of parallel filaments.

In many cellular processes that involve the deformation of membranes or the movement of vesicles and organelles, the energy from biochemical reactions is converted into forces. The biological filaments called actin are one of the major force producing machineries of the cell. It is commonly believed that the elongation of these filaments at their tip is the only way actin filaments can exert force. However, the amount of force produced by this mechanism can only account for a small fraction of the force in key cellular processes, such as clathrin-mediated endocytosis. In this paper, we demonstrate that connecting actin filaments with each other with flexible proteins called crosslinkers is a new way to transform biochemical energy into mechanical energy, and that this stored mechanical a1111111111 a1111111111 a1111111111 a1111111111 a1111111111

Introduction
The cytoskeleton protein actin assembles into three major structures in yeast cells, including endocytic actin patches, actin cables, and the contractile ring [1,2]. In actin cables and the contractile ring, formin-nucleated actin filaments are crosslinked into long bundles with a length on the order of microns [3][4][5]. Computational models of these actin structures typically treat actin filaments as semi-flexible polymers that are connected by rigid segments [6][7][8][9]. In contrast, the organization of the actin network in actin patches formed during clathrin-mediated endocytosis is drastically different from that in actin cables or the contractile ring. The length of filaments in actin patches is strongly limited by capping and severing proteins [10], and mathematical modeling predicted that the average length of filaments is less than 200 nm [11]. Filaments of this length scale can be considered as straight rods, because the persistence length of actin filaments is on the order of 10μm [12][13][14], which allows them to sustain forces larger than 10pN without buckling [15].
At the endocytic actin patch, a small area of the flat plasma membrane invaginates towards the cytoplasm upon assembly of actin. In budding yeast, the invagination elongates up to 140nm in depth, and then is pinched off, releasing a tear-shaped vesicle [16]. Actin is essential for many of these steps, from the initiation of invagination to vesicle scission [17][18][19]. Despite extensive experimental work that characterized the overall dynamics of assembly, disassembly and ensemble movements of proteins of the actin meshwork [19][20][21][22][23][24][25][26], the precise structural organization of actin filaments within the endocytic patch remains unknown. Indeed, individual filaments are not resolvable even in electron micrographs, in which the actin network appears as a ribosome-exclusion zone, which is about 200nm in depth and 100nm in width [16].
Actin crosslinking proteins play a crucial role in determining the mechanical responses of the actin network to force perturbation [27][28][29]. Fimbrin (Fim1p) is the second most abundant protein recruited to the endocytic patch during clathrin-mediated endocytosis in fission yeast, after actin [20]. It has two actin binding domains that enable it to crosslink adjacent filaments. Deletion of fimbrin results in significant defects in endocytic internalization [24,30,31]. In vitro experiments have shown that fimbrin efficiently bundles long actin filaments, but bundling efficiency is reduced in the presence of capping protein as a result of decreased filament length [31]. It remains unclear how this length-dependent bundling activity arises and how this activity is related to the role of fimbrin during clathrin-mediated endocytosis.
Internalization of the endocytic membrane is hindered by the high turgor pressure (P * 0.8 × 10 6 Pa [32,33]) inside yeast cells [34]. Under such high pressure, theoretical studies suggest that the force needed to initiate membrane invagination is on the order of 3000pN [35,36] and actin polymerization is thought to provide the driving force. However, assuming no more than 150 filaments are simultaneously generating the force [11,20,21], each of these filaments has to generate a force of at least 20pN, an order of magnitude larger than the maximum polymerization force of *1pN of actin filaments measured in vitro [37]. This number of 20pN is likely an underestimate since the calculation here uses a very generous estimate for the number of growing filaments, up to 20-fold of what mathematical modeling predicts [11]. Therefore, actin polymerization alone is not enough to provide the force necessary to elongate a clathrin-coated pit. Even though type-I myosins participate in endocytosis, their low power output over a narrow range of forces suggest that they are more likely force-sensing tethers rather than force generators [38][39][40].
In this paper, we present a computational model for dynamic crosslinking of rigid actin filaments in conditions relevant to clathrin-mediated endocytosis. We show that kinetic and mechanical properties of the crosslinkers finely tune the structural transition of actin network between bundles and meshworks. In addition, we show that the chemical binding energy is converted into elastic energy upon binding of crosslinkers. The elastic energy stored in individual crosslinkers is significantly higher than their thermal energy. This surprising property is a consequence of the helical pitch of actin filaments, which leads to torsional strains between crosslinkers attached to a common pair of filaments. We discuss the mechanical implications of these torsionally stressed crosslinkers and propose a possible mechanism to generate directed rotation of filaments by orderly detaching the crosslinkers.

Model of crosslinked rigid actin filaments
We model actin filaments as rigid cylindrical rods with subunits that carry a helical pitch ( Fig 1A). The length of filaments is restricted to the range of actin filament size during endocytosis (81nm to 216nm), around 2 orders of magnitude shorter than their persistence length *10μm. Filaments of this length scale remain virtually straight under *10pN of force and untwisted under *100pN Á nm of torque (see the Methods section). Therefore we can neglect bending and twisting, and describe the motion of a filament by its translational velocity V c of Actin filaments are modeled as rigid rods made of subunits carrying an orientational vector O i that represents the normal vector to the binding surface of the actin crosslinker. For visualization, each subunit is depicted as a sphere but is actually a disk-shape with a diameter of b = 6nm and a height of δ = 2.7nm. The orientation of a filament is described by the unit vector N, pointing from the pointed end (-) to the barbed end (+), and the first subunit's orientational vector M = O 1 . Two consecutive subunits have an angle of π14/13 in their orientations. (B) Crosslinker turnover is described by stochastic formation and breakage of bonds between two actin subunits in different filaments, with rate constants k f and k b , respectively. (C) Each crosslinker is modeled as a combination of three springs, one extensional spring with stiffness κ ext , which characterizes the stretchiness l c between the two actin binding domains, and two torsional springs with stiffness κ tor , which characterizes the flexibility of the angles θ i and θ j between the axis of the crosslinker (which links the centers of each subunits it is attached to) and the vector normal to the binding surface of each actin subunit in each filament ( center of mass, and angular velocity O relative to the center of mass. Details of the model can be found in the Methods section. We model actin crosslinkers as elastic springs that connect two actin subunits in different filaments. Each spring has three elastic components: one in extension, which represents how much the spring is stretched, and two in torsion, each one representing how much the orientation of the axis of the crosslinker linking both actin subunits deviates from the vector normal to the binding interface of each actin subunit (Fig 1C). This elasticity is the simplest model to account for: (i) the intra-molecular flexibility between both actin binding domains, and (ii) the flexibility in the binding interface between each actin binding domains and the actin subunits they are bound to. Specifically, the elastic energy E of a crosslinker is composed of an extensional part E ext and a torsional part E tor , E = E ext + E tor . The extensional energy accounts for (i), where κ ext denotes the extensional stiffness, l c denotes the length, and l 0 denotes the rest length of the crosslinker (Fig 1C). The torsional energy accounts for (ii), where κ tor denotes the torsional stiffness, θ i and θ j denote the angles between the actin subunits and the axis of the crosslinker (Fig 1C). Note that our model does not take into account the stiffness of the rotation of filaments around the axis of the crosslinker. Doing so would force us to consider a preferred relative orientation of the filaments with each other, and strongly favor either bundles or meshworks. No experimental values for the extensional and torsional stiffnesses of fimbrin are available in the literature. For our simulations, we used values within a few orders of magnitude of stiffnesses measured for fascin and antibodies [41,42]. Crosslinker turnover is modeled as Poisson processes with a crosslinker formation rate constant k f and a breakage rate constant k b , which increases exponentially with the total elastic energy ( Fig  1B). A more detailed description is presented in the Methods section. The length of filaments is kept constant in a given simulation and the total number of actin subunits in the filaments is fixed to N actin = 7000 for all simulations [20]. The maximum occupancy of crosslinkers on filaments is constrained to remain under 25% (or 1 crosslinker for 4 subunits), which is equivalent to a maximum of 875 attached crosslinkers, close to the peak value *900 measured experimentally [20]. We initiate each simulation with uncrosslinked filaments that are randomly positioned and oriented. Reflecting boundary conditions are imposed to ensure filaments stay in a cubic box of 500nm in size. Simulations are performed using the reference values listed in Table 1 unless otherwise mentioned.
To quantify the organisation of actin filaments, we introduced the global and local nematic order parameters S global and S local . These quantities characterize the degree of alignment between filaments in the entire simulation space, or in a local neighborhood, respectively. Their values range from 0 to 1, and a larger value indicates that filaments are more aligned with each other globally, for S global , or locally, for S local . Note that, in practice, the minimum reachable value for S local is usually close to 0.4 (S6 Fig). The detailed mathematical definition of these parameters can be found in the Methods section.
Actin assembly and disassembly takes tens of seconds during endocytosis. We set the total simulation time to be 50s. Within this period, in most of our simulations, the metrics reach steady state. However, there are cases where the system is in a transient state, mostly due to slow convergence of the global nematic order parameter S global (S8 Fig). Therefore, we chose the local nematic order parameter S local averaged from 40s to 50s as the major metric to characterize the organisation of actin networks. Conclusions based on S local are more robust than based on S global .

Long filaments form bundles and short filaments form meshworks
We first studied how filament length influences the structure of actin clusters. For crosslinking rates k f below 0.1s −1 , the number of attached crosslinkers remained small, and actin filaments did not organize into higher order assemblies (S1 Fig). When the crosslinking rate k f was high enough, initially disconnected short filaments (81nm) quickly formed small clusters that eventually coalesced into three to four larger clusters (Fig 2A). The number of attached crosslinkers rapidly saturated to a dynamic steady-state (Fig 2C), as the crosslinkers underwent constant turnover. Filaments within the cluster are organized into a disordered meshwork, which is characterized by a small local nematic order parameter (S local * 0.4, Fig 2B). Long filaments (216nm) rapidly aligned with their neighbors into small bundles (Fig 2D), as indicated by the fast convergence of the local nematic order parameter S local to its steady state value around 1 (Fig 2E). Throughout the simulation, these locally aligned filaments remained connected with each other, and the number of clusters remained small ( Fig 2F). These connected bundles slowly adjusted their orientations to eventually coalesce into a few large bundles.
When filament length was increased from 135nm to 162nm, the structure of the actin network transitioned from meshwork to bundle, indicated by the sharp increase of the local nematic order parameter S local from *0.5 to above 0.9 ( Fig 2G). The global nematic order parameter S global showed similar trend as S local , but with a smaller magnitude, since filaments formed 2 to 4 independent clusters ( Fig 2H). Altogether these results show that crosslinked actin filaments with a size and crosslinker density comparable to what is measured during

Crosslinkers with high extensional stiffness drive bundle formation
Next, we explored the influence of the mechanical properties of actin crosslinking proteins on the organization of actin filaments. The organization of medium length filaments (135nm) varied dramatically for different combinations of κ ext and κ tor (Fig 3A-3C). For low extensional stiffness κ ext , actin filaments organized into a meshwork (Fig 3A), while they formed bundles for higher κ ext values ( Fig 3B). The transition between meshwork and bundle was tightly controlled, since the local nematic order parameter had a sharp increase around κ ext = 0.1pN/nm ( Fig 3E, red) at given κ tor = 10pN Á nm Á rad −1 . This sharp transition was even more pronounced for longer filaments (189nm), as reported in both local and global nematic order parameters ( Fig 3E and S2A Fig, orange). In contrast, when filaments were short (81nm), the transition was relatively smooth (Fig 3E, blue). From the above results, we conclude that crosslinkers with large extensional stiffness favor bundle formation. This result can be intuitively explained by the following simplified but heuristic example involving only two filaments. If two filaments are initially aligned with each other, a slight change in orientation between both filaments results in the stretching of crosslinkers bound at different positions, which leads to a restoring torque to realign the filaments. The torque is proportional to the extensional stiffness and the distance between the positions of attached crosslinkers, thus stiffer crosslinkers create larger realignment torque than softer crosslinkers. Longer filaments not only have more crosslinkers, but also crosslinkers that are more distantly positioned, therefore more extended and more inclined to create a larger torque that will restore the parallel alignment of filaments.

Crosslinkers with high torsional stiffness disfavor bundle formation
We next investigated the impact of torsional stiffness κ tor on the organization of actin networks, keeping the extensional stiffness at a relatively small value (κ ext = 0.1pN/nm). Torsional stiffness had virtually no influence on the organization of short filaments (81nm) (Fig 3F, blue). However, medium length filaments (135nm) had a sharp transition from bundle to meshwork at κ tor = 10pN Á nm Á rad −1 (Fig 3F, red). A similar trend was also observed for long filaments (189nm) (Fig 3F, orange). At highest torsional stiffness tested (κ tor = 100pN Á nm Á rad −1 ), crosslinker attachment lifetime was extremely short because their torsional energy often became much larger than the critical energy E c that modulates the detachment rate (see Eq 16), therefore the number of attached linkers was significantly reduced (S2H Fig). The uncrosslinked actin network formed at this regime was different from the connected meshwork formed at κ tor = 10pN Á nm Á rad −1 , though their nematic order parameters S local were similarly low.
The above results show that crosslinkers with high torsional stiffness disfavor bundle formation. This result can be explained by the fact that, when torsional stiffness is high, formation of several crosslinks between two aligned filaments results in very high torsional energies, due to the frustrated interactions between crosslinkers. This will become clear later in the paper. Therefore, it is energetically more favorable to form a few but over-stretched crosslinkers with many distant filaments than to form many but under-stretched crosslinkers with a few function of filament length. (H) Number of attached crosslinkers (magenta, left axis) and number of clusters (black, right axis) as a function of filament length. In (G) and (H), for each simulation, the means of the metrics were calculated from the data between 40s to 50s and the error bars indicate standard deviation over 10 simulations.

Turnover of crosslinkers is necessary for large bundle formation
We have shown that high extensional stiffness of crosslinkers favors bundle formation. When the breakage rate k 0 b was reduced to 0.01s −1 , even for large extensional stiffness (1pN/nm), filaments formed a structure where small bundles were interconnected but did not align with each other (Fig 3D), as indicated by the relatively low local nematic order parameter at k 0 b ¼ 0:01s À 1 compared with S local at k 0 b ¼ 10s À 1 (Fig 3G). The difference in the global nematic (E) Local nematic order parameter S local as a function of the extensional stiffness κ ext , with κ tor = 10pN Á nm Á rad −1 and k 0 b ¼ 10s À 1 . (F) Local nematic order parameter S local as a function of the torsional stiffness κ tor , with κ ext = 0.1pN/nm and k 0 b ¼ 10s À 1 . (G) Local nematic order parameter S local as a function of the strain-free breakage rate k 0 b , with κ tor = 10pN Á nm and κ ext = 1pN/nm. In (E-G), simulations were performed for filaments of various lengths: 81nm (blue), 135nm (red), 189nm (orange). For each simulation, the means of S local were calculated from the data between t = 40s to 50s and the error bars indicate standard deviation over 10 simulations. S local corresponding to the networks in panels (A-D) are identified by the red symbols with corresponding shapes.
https://doi.org/10.1371/journal.pcbi.1006150.g003 order parameters S global was even more pronounced for long filaments (189nm) (S2C Fig,  orange). At very high breakage rate, the lifetime of bonds between filaments was so short that filaments formed an essentially random, uncrosslinked network. Altogether, we conclude that crosslinker turnover is essential for bundle formation, as alignment of bundles requires the breaking of the bonds that disfavor filament alignment.

Phase diagram of filament organization
Parameter dependence of the actin network structure is summarized in Fig 4 where we plotted the local nematic order parameter S local and the number of attached crosslinkers N attach as a function of crosslinking rate constant k f and filament length L. For a combination of low extensional and high torsional stiffnesses (κ ext = 0.1pN/nm, κ tor = 10pN Á nm Á rad −1 ), values of S local are concentrated either close to 1 or close to 0.5 (Fig 4A, yellow and blue regions, respectively). These two regions are separated by a narrow transition band around S local = 0.75 (Fig 4A, green). Values of N attach are clustered either close to the saturation number 875, or less than 100, divided by a transition band around N attach = 300 (Fig 4B, yellow, blue and green regions, respectively). Therefore we chose the lines S local = 0.75 and N attach = 300 as the boundaries to define the phase diagram ( Fig 4C). These two lines divide the parameter space into three regions: (1) Above the line S local = 0.75, actin filaments are locally aligned into a Structural organization and energy storage in crosslinked actin assemblies bundle; (2) Between the lines S local = 0.75 and N attach = 300, filaments form a crosslinked, disordered meshwork; (3) Below the line N attach = 300, filaments are essentially uncrosslinked over the entire simulation time. For larger extensional stiffness (κ ext = 1pN/nm), the relative positions of the regions are similar, but most of the phase diagram corresponds to bundles (Fig 4F), and only the shortest filaments (81nm) form meshworks. In all cases, the existence of three regions in these phase diagrams requires moderate to high breakage rates. When the breakage rate is low (k 0 b ¼ 0:01s À 1 ), meshworks occupy the entire parameter space (S3 Fig).

Crosslinked actin networks store elastic energy
Though crosslinkers in our model were not active elements, we found that crosslinkers rapidly became stretched in length, and twisted in angle. To quantify the crosslinkers' deformations, we introduced the extensional strain = (l c − l 0 )/l 0 , which measures the relative change of a crosslinker's length l c from its rest length l 0 , and the torsional strain θ, which measures the angle between the crosslinker and its bound actin subunits. For all stiffness values, the distributions of the extensional strain and the torsional strain θ significantly deviated from the corresponding Boltzmann distribution for a single independent free spring (Fig 5A-5C). Strikingly, the extensional strain peaked at *0.5 but not zero (Fig 5A-5C, top), indicating that the crosslinkers were stretched on average. The distribution had a narrower width for higher extensional stiffnesses κ ext . The peak of torsional strain θ decreased with increasing torsional stiffness, while the widths of the distributions were essentially the same as in the Boltzmann distribution (Fig 5A-5C, bottom).
The deviation from the Boltzmann distribution can be primarily accounted by the coupling between crosslinkers that are attached to the same pair of filaments. Indeed, a simpler 1D example of the Brownian motion of two particles each subject to a spring follows similar properties (S5 Fig). In this example, each spring generates a force of −κx if the particle is displaced from the equilibrium position 0 to x. If the movement of the two particles were independent (S5A Fig), the joint distribution of their positions x 1 and x 2 would simply be the product of the Boltzmann distribution of an individual particle. However, if the movements of the two particles are coupled, for instance, subject to the constraint x 1 − x 2 = 2x 0 (S5B Fig), the position distribution of particle 1 becomes p(x 1 ) = Cg(x 1 )g(x 1 − 2x 0 ), where C is the normalization constant. Therefore particle 1 is displaced by x 0 due to the coupling with particle 2. Similarly, if two rigid filaments are bound by several crosslinkers, the extensional and torsional strains of these crosslinkers are coupled, and this coupling gives rise to significant strains in the crosslinkers. We then determined the dependence of the elastic energy stored in the crosslinkers on the crosslinker stiffness. At fixed torsional stiffness κ tor = 10pN Á nm Á rad −1 , the average magnitude of the extensional strain || decreased with increasing extensional stiffness κ ext (Fig 5D). However, the average extensional energy per crosslinker E ext ¼ 1 2 k ext l 2 0 2 increased (Fig 5E), from *0.3k B T for soft extensional springs up to *5k B T for stiff ones, with only a weak dependence on filament length (Fig 5E). The average torsional energy per crosslinker reads E tor ¼ 1 2 k tor ðy 2 i þ y 2 j Þ, with θ i and θ j being the average torsional strains at the two actin subunits. They have the same magnitude θ. At fixed extensional stiffness κ ext = 0.1pN/nm, E tor increased with torsional stiffness κ tor from *0.2k B T to up to *10k B T, with again a weak dependence on filament length (Fig 5G). Both extensional and torsional energies vary relatively smoothly with the corresponding stiffnesses over three orders of magnitude, which is in stark contrast 2k B T , where C 1 and C 2 are the normalization constants for 2 (−1, 1) and θ 2 (0, π) respectively. Note that for the extensional strain, when we calculate the Boltzmann distribution, the energy contribution from with the sharp structural transition between meshwork and bundle when stiffnesses are varied over the same range (Fig 3E and 3F).
Compared with the increase of the extensional energy E ext with κ ext (*8 − 17 fold), the increase of torsional energy E tor with κ tor (*40 − 50 fold) was more pronounced. This was likely due to stronger coupling between torsional strains of crosslinkers than between extensional strains. Elastic energy plateaued and then slightly decreased for very high stiffnesses (κ ext = 10pN/nm or κ tor = 100pN Á nm Á rad −1 ), which was the consequence of higher detachment rate leading to a smaller number of attached linkers, thus reducing the frustrated interactions (S2 Fig).

Possible mechanism for torque production by release of the stored elastic energy
We have shown that crosslinking of filaments leads to elastic energy stored in the crosslinkers. How can this energy be transformed into mechanical work? Here, we propose a mechanism for torque generation through orchestrated detachment of crosslinkers by studying a simple model with only two filaments. Let us consider a pair of short filaments, where every other subunit of each filament is crosslinked (Fig 6). For simplicity, we assume the two filaments are parallel and consider only the rotation of filaments around their axes. We show that consecutive detachment of crosslinkers from the pointed end to the barbed end lets the filament rotate in the same direction by π/13 for each detachment (Methods).
Building on this simple proof of principle, we show that sustained directional rotation can be achieved with any filament length and crosslinker spacing and configurations such that (i) the angles between two consecutive crosslinkers along the pair of filaments have the same sign, and (ii) the sum of all these angles is smaller than 2π. Under these conditions, breakage of crosslinkers from one end to the other produces directional torque (see Methods).

Organization of actin filaments in diffraction limited assemblies
In this paper, we showed that highly crosslinked actin networks made of rigid filaments (< 200nm) can form either disordered meshworks or ordered bundles, depending on the filament length and the mechanical and kinetic properties of the crosslinkers. A recent in vitro study of short actin filaments (200nm) showed that with increasing density of filamin, the initially sparsely distributed actin filaments condensed into a spindle-shaped aggregate, in which the organization of actin filaments displayed nematic order [50]. This observation is consistent with the phase diagram in Fig 4C of our study for filaments of 200nm with increasing crosslinking rate.
We investigated the possible structures formed by actin filaments in the presence of elastic crosslinkers. However, our study does not take into account other actin regulating proteins involved in endocytosis, such as the Arp2/3 complex and capping protein, and thus does not completely resolve the question of the organization of filaments crosslinked with fimbrin at steric interactions, which lead to the empty region at highly negative strains in the histogram, is neglected. (D, E) Average absolute value of the extensional strain (D) and the corresponding extensional energy E ext (E) as a function of the extensional stiffness κ ext . (F, G) Average torsional strain θ (F) and the corresponding torsional energy E tor (G) as a function of the torsional stiffness κ tor . In (D-G), simulations are performed for filaments of various lengths: 81nm (blue), 135nm (red), and 189nm (orange). For each simulation, the means of the energy were calculated from the data between 40s to 50s and the error bars indicate standard deviation over 10 simulations. Energies corresponding to the networks in panels (A-C) are identified by the red symbols with corresponding shapes. https://doi.org/10.1371/journal.pcbi.1006150.g005 Structural organization and energy storage in crosslinked actin assemblies the site of endocytosis in yeast. Even with this limitation, these results provide valuable insights about possible actin filament architectures for endocytosis and other cellular processes that involve short actin filaments. Using the rate constants for fimbrin that have been measured in vitro [31,45], and the fimbrin concentration in fission yeast cytoplasm (3.7 μM) [20], these values correspond to rates in our model k f = 0.2s −1 and k 0 b ¼ 0:04s À 1 . Our simulations suggest that the slow off-rate of fimbrin should favor an assembly of actin filaments into a meshwork (S3 Fig). Further simulations with branched filaments, and with geometries and dynamics more representative of endocytosis will tell us which type of structure is present at endocytic Structural organization and energy storage in crosslinked actin assemblies sites. In addition, further experimental characterization of the mechanical properties of fimbrin and other crosslinkers will be key to understanding the self-organization of actin filaments in diffraction limited structures, and to test the predictions of our simulations.

Mechanisms of energy storage by actin crosslinkers
Our simulations demonstrate that individual actin crosslinkers are able to store up to 10k B T of elastic energy, which is one order of magnitude higher than the elastic energy stored in an uncoupled spring in a thermal bath (1.5k B T), and about half of the energy released by ATP hydrolysis (*25k B T).
To get a better sense of the amount of energy stored in the crosslinkers, one can make a comparison with the energy necessary to deform the plasma membrane into an endocytic vesicle. A back of the envelope calculation estimates the work needed to create a cylindrical invagination of R t = 25nm in radius and D t = 140nm in depth [16] against the turgor pressure P * 0.8 × 10 6 Pa is PpR 2 The results of our model suggest that crosslinking filaments once with *900 crosslinkers, around 10 4 k B T energy could be stored, or about 1/6 of the total energy needed.
One may wonder where this large elastic energy comes from. We can show in a simplified model that the chemical binding energy of crosslinkers is indeed the source of the elastic energy. Let us consider two filaments with fixed positions and orientations, each having N subunits. In the following, we will consider the chemical balance between configurations where there is either n = 0 or 1 crosslinker between the two filaments. The rate at which a crosslinker is formed is k f Γ(1), where Γ(n) denotes the number of possible pairings to form n crosslinkers between subunits in the two filaments. The value of Γ(1) depends on the orientations and positions of the two filaments, and the mechanical properties of the crosslinkers, and could vary from 1 for very stiff crosslinkers and orthogonal filaments, to N 2 for infinitely soft crosslinkers and parallel filaments. Here we assume Γ(1) = N, which implies that for each subunit in one filament, there is a unique subunit in the other filament that is within the reaction distance to allow a crosslinker to be formed. When one crosslink is formed, its detachment rate is k 0 b e E=E c , assuming the crosslinker bears an elastic energy E. If the attachment rate is greater than the detachment rate, the system is more likely to be crosslinked, even though the elastic energy E stored in the crosslinker tends to drive down the crosslinker occupancy. Noting m eff E c ln ðNk f =k 0 b Þ the effective chemical binding energy that tends to drive up the occupancy of crosslinkers, this condition can be expressed as comparison between the effective chemical binding energy and the elastic energy μ eff > E. Using the parameters k 0 b ¼ 10s À 1 , k f = 1s −1 , E c = 10k B T and N = 50, we estimate μ eff = 16k B T for the first crosslink formation, which is larger than the 10k B T of elastic energy per crosslinker computed in our simulations. If we now consider the case where there are n crosslinkers formed between the filaments, the effective energy for binding an extra crosslinkers becomes m eff ðnÞ ¼ E c ½ ln ðk f =k 0 b Þ þ ln ðGðn þ 1Þ=GðnÞÞ. The second term in the bracket represents an entropic contribution that comes from the different ways of building n or n + 1 crosslinks between subunits of both filaments. This simplified two filament system illustrates how in our simulations with multiple filaments crosslinker occupancy is driven up by a similar entropic contribution in the chemical binding energy.
Our model considers "slip-bond" crosslinker detachements, i.e. crosslinkers are more likely to detach if force and torque are exerted on them (Eq 16). However, we could consider the case of "catch-bond" detachments, where crosslinkers are less likely to detach under force and torque, as it has been shown for some cytoskeleton proteins [40,51,52]. In this case, we expect that the elastic energy stored in crosslinkers would be larger than what we have observed in our simulations, and the conformational change required for the catch-bond behavior would increase the available energy limit.
The main reason energy storage is possible is that short actin filaments are rigid, which creates geometrical constraints on bound crosslinkers, forcing virtually all of them to fluctuate around average lengths and angles that are different from their rest lengths and angles. This implies that the crosslinkers are rigid enough to store elastic energy, but not as rigid as the filaments, so that filaments cannot be twisted or bent when crosslinked, or the distance over which filaments are twisted and bent in order to reduce the frustration, noted L tb , is much longer than the filament length L. In the opposite case, when filament length L ) L tb , filaments could form bundles in which individual filaments are twisted. Experiments conducted at this regime suggest that the frustrated interaction serves as a mechanism to control the size of the bundle [53] and cooperative binding of actin crosslinkers [54]. The two regimes have been theoretically studied by C. Heussinger and G. Grason [55].

How can the stored energy be used productively in cellular processes?
Only a small fraction of the force necessary to deform the plasma membrane during clathrinmediated endocytosis in yeast can be accounted for by actin polymerization alone. We predict that at least some of the missing force can come from the conversion of the elastic energy stored in the crosslinkers into force and/or torque.
In this paper, we proposed a specific mechanism for torque production by orchestrated detachment of crosslinkers. This mechanism is different from the Brownian ratchet mechanism of force production that is directly coupled to ATP hydrolysis [56,57]. However, the ordered detachment has to be coupled to a non-equilibrium process to provide the information necessary for the ordered detachment. Treadmilling of filaments coupled to ATP hydrolysis could play such a role. To estimate the order of magnitude of free energy necessary to provide this information, let us consider a pair of short actin filaments (e.g. 50-subunit long) that are crosslinked by 10 crosslinkers. The free energy cost of detaching the crosslinkers in a specific order among all the 10! possibilities is *k B T ln 10! = 15k B T, which is only a small fraction of the energy provided by the ATP hydrolysis of two actin filaments undergoing treadmilling (50 × 25 = 1250k B T).
We stress that this ordered detachment is only one possible mechanism to use the energy and more mechanisms need to be discovered. Future work with more realistic models for endocytosis or other actin-based processes will likely uncover new orchestrated mechanisms for force production. In our system, crosslinkers are driven up to a mechanically pre-stressed state by chemical binding energy. In principle, in order to release the elastic energy stored in crosslinkers, change in energetics of crosslinker binding/unbinding is necessary to induce collective detachment of crosslinkers. When actin polymerization is considered, ATP-bound actin is incorporated at the barbed end and undergoes hydrolysis after incorporation. The nucleotide content change could alter the crosslinkers' binding affinity, causing rapid detachment of crosslinkers, which is accelerated by the pre-stressed structure. Elastic energy released during this process could be converted into work by the reorganization of actin filaments.

Conclusion
We developed a computational model to study the dynamic assembly of actin filaments mediated by elastic crosslinkers. The organization of actin filaments were classified into either a meshwork or a bundle, characterized by their nematic order parameter and the number of attached crosslinkers. We showed that the elastic energy stored in crosslinkers increased with their stiffness due to coupling between crosslinkers bound to rigid filaments. As a proof of principle, we showed that the elastic energy could be converted into mechanical work by orchestrated detachment of crosslinkers between two parallel filaments. Our work provides a new perspective to study the mechanisms of force and torque production by actin filaments, in addition to the traditional end polymerization. It also provides an alternative energy source to account for the insufficient force production by actin polymerization during clathrin-mediated endocytosis.

Model of actin filaments
Actin filaments are modeled as rigid cylindrical rods with diameter b and length L. The position of a filament is represented by its center of mass C. A unit vector N pointing from the filament's pointed end to the barbed end indicates the orientation of the filament. The i-th subunit (counting from the pointed end) carries a unit vector O i , which is normal to the binding surface with a crosslinker (Fig 1A). We assume all O i -s are perpendicular to the filament's orientation N. Based on the atomic structure of actin filaments [43], two consecutive subunits O i and O i+1 span an angle of 14π/13 calculated counter-clockwise from O i to O i+1 . This means two consecutive subunits on different strands have their binding interface in almost opposite directions, and two consecutive subunits on the same strand have their binding interface at an angle of 2π/13 (28˚). We arbitrarily choose the filament's rotational vector M as the normal vector of the first subunit M = O 1 . Thus the orientational degree of freedom of the filament is fully captured by three orthonormal vectors N, M and N × M.
The motion of a filament is described by its translational velocity V c and angular velocity O, which are defined by the following equations: The velocities V c and O are governed by the force-balance and torque-balance equations: Here the 3 × 3 matrices X t and X r denote the frictional matrix associated with translational and rotational motion of the filament, respectively. The vectors F e and T e denote the total deterministic force and torque generated by crosslinkers or induced by steric interactions between filaments. The vectors F s and T s denote the stochastic force and torque, which obey the fluctuation-dissipation relations: Here the subscript indicates the element of the vectors or matrices. The frictional matrices are anisotropic, and given by [58]: where x t k and x t ? are the frictional coefficients for translational movement parallel with and perpendicular to the filament's central axis, and x r k and x r ? are the corresponding frictional coefficients for rotation. The 3 × 3 identity matrix is denoted by I, and denotes the outer product of two vectors. The anisotropic frictional coefficients depend on filament length L and diameter b via the relations [43]: Here η denotes the viscosity of the medium.
To account for the steric interaction between filaments, if the shortest distance r min between two filaments is less than the diameter b of a filament, a constant repulsive force f st is applied along the lines connecting the two nearest points.
In each time step, we calculate all the forces and torques acting on a filament and determine the translational velocity V c and angular velocity O of the filament according to Eqs (6) and (7). The center of mass of a filament is then updated as: The updated orientations are: where Rot (OΔt) denotes the rotation matrix defined by the vector OΔt. The rotation vector M(t) is updated in the same way.

Model of actin crosslinking proteins
Each actin crosslinking protein is modeled as an elastic spring that bridges two actin subunits in two separate filaments. The crosslinking of two unoccupied subunits proceeds with a rate constant of k f , as long as the subunits are less than r c apart. The breakage of an established crosslink is assumed to follow a "slip-bond" mechanism and occurs with an energy-dependent rate constant: where k 0 b denotes the strain-free breakage rate constant, E denotes the total elastic energy, and E c denotes the critical energy that determines the sensitivity of the bond breakage on the forces and torques.
The elastic energy E of a crosslinker that bridges actin subunits in filaments α and β is a function of the positions, orientations and rotations of both filaments, as well as its positions in the filament, E = E(C α/β , N α/β , M α/β ). The force generated by the crosslinker on filament α reads: To determine the torque generated by the crosslinker on filament α, we choose three orthnormal vectors e 1 , e 2 , e 3 and virtually rotate filament α by an infinitesimal angle ϕ i around the axis e i . These operations are equivalent to applying the following infinitesimal changes to the orientational vectors of filament: The elastic energy correspondingly has an infinitesimal change E ! E + ΔE. The torque then reads: Forces and torques acting on filament β are derived in a similar way. The total elastic force and torque are obtained by summing (17) and (20) over all the crosslinkers bound to the filament.

Events performed in a single simulation step
At each time step Δt of the simulation, we perform the following operations: 1. Calculate all the deterministic and stochastic forces and torques on each filament. These forces and torques determine the translational and angular velocities of filaments according to Eqs (6) and (7). The positions and orientations of filaments are then updated according to Eqs (14) and (15).
2. For each pair of filament α and β, determine the number N a of crosslinkers to be attached between them by drawing a random number from the Poisson distribution, PoisðN a mon N b mon k f DtÞ, where N a mon and N b mon are the number of actin subunits in each filament. Then randomly pick N a pairs of subunits (i, j) with i in filament α and j in filament β. If they are not already occupied by a crosslinker and distance between them is shorter than the interacting distance r c and the total number of already occupied subunits in each filament is less than g max N a=b mon , a crosslinker is built between the subunits (i, j). 3. For each already existing crosslinker, determine whether it will be detached by comparing a uniformly distributed random number u on the interval [0, 1] with the energy-dependent breakage rate k b in Eq (16). If u < 1 − e −k b Δt , the crosslinker is detached.
In our simulation, we always set the time step Δt at least 100 times smaller than the relaxation time of the spring t ¼ minðx t k=? =k ext ; x r k=? =k tor Þ to ensure that we correctly capture the dynamics of the springs. For computational reasons, we used a high viscosity value η = 10Pa Á s, such that the relaxation time τ * 0.01s for κ ext = 0.1pN/nm and κ tor = 10pN Á nm Á rad −1 . We tested values of lower viscosity down to η = 0.1Pa Á s. There is no significant difference in the local nematic order parameter S local and the elastic energies E ext and E tor between η = 0.1Pa Á s and η = 10Pa Á s. However, the global nematic order parameter S global for long filaments is increased to 1 at lower viscosity (S7 Fig). This is because the enhanced diffusion increases the probability of filaments moving close to each other. As a result, the separated bundles observed at high viscosity merge into a single bundle when the viscosity is low, increasing the global nematic order parameter.

Metrics to describe local and global organization of filaments
We characterize the structure of actin clusters by introducing local and global nematic order parameter S local and S global . We map the connections between filaments into an undirected graph, with filaments being the nodes, and the number of crosslinkers being the value of the edges connecting two nodes. Filaments in a connected component of the graph are said to form a cluster if the number of filaments in the component is more than 10. The nematic order parameter S for a group of filaments is the maximum eigenvalue of the following matrix [59]: where G denotes the number of filaments in the group, N α denotes the orientational vector of filament α. For global nematic order parameter S global , the group in Eq (21) includes all the filaments. For a particular filament α Ã , S a Ã local is defined by grouping the filament α Ã and its connected nodes in Eq (21). The local nematic order parameter S local is the average of S a local over all the filaments that have at least 2 connected nodes.
Both S local and S global are in the range of [0, 1]. Values of S local close to 1 indicate that filaments are locally aligned with their connected neighbors. Values of S global close to 1 indicate that all the filaments are aligned. In general S local is greater than S global , and reaches steady state more rapidly, because filaments that are in close proximity can rapidly align, but it takes time for distant clusters of filaments to collide and reorient. S local is also more consistent over different simulations than S global , as reported by smaller error bars for S local than for S global (e.g. Fig  2G). If the viscosity of the medium η is reduced to 1Pa Á s, filaments form a single cluster and the error bars of S global become comparable with S local (S7B and S7C Fig).
Note that in a sparsely connected network with filaments in random orientation, S local * 0.5 (S1B Fig) is higher than one should expect (*0). This artifact is due to the fact that the sum in Eq (21) is done over a very small number of filaments (*3). We confirmed this property by numerically calculating the nematic order parameter for three unit vectors with random orientations. The resulting distribution of S local has a peak at 0. 45 (S6 Fig). This almost uncrosslinked network should be distinguished from the densely connected actin meshworks, which have local nematic order parameters S local in the same range (* 0.4 Fig 2B) but possess a large number of attached crosslinkers. Therefore, the number of crosslinkers in the meshwork is required to distinguish these two structures.

Validity of the rigidity assumption
In our model, we assumed that filaments are rigid. This rigidity assumption implies that (i) thermal fluctuations, and (ii) forces and torques exerted by crosslinkers do not significantly bend or twist the filaments. We can verify a posteriori that these conditions are actually fullfiled in our simulations. Indeed, the maximum force produced by a crosslinker in our simulations is *10pN when the extensional stiffness κ ext reaches 10pN/nm, and the maximum torque is *100pN Ánm when the torsional stiffness κ tor reaches 100pN Á nm Á rad −1 . Given the persistence length of actin filament for both bending and twisting is L p * 10μm [12,14,60,61], for a filament of length L = 200nm which consists of N = L/δ = 74 subunits, the angular change between two consecutive subunits due to thermal fluctuation is arccos ðe À L=L p Þ=N ¼ 0: 15 . The angular change due to bending caused by a force of f = 10pN in the middle of the filament when the two ends are fixed is arctan ð fL 2 48L p k B T Þ=N ¼ 0:15 [62].
The twisting angle by a torque of T = 100pN Á nm is ð TL L p k B T Þ=N ¼ 0:37 . Therefore, it is safe to consider filaments as stiff, and the energy stored in crosslinkers would not be dramatically different even if the finite stiffness of filaments was taken into account.

Direct rotation of filaments by consecutive detachment of crosslinkers
We consider the rotation of two parallel filaments around their axes by consecutive detachment of crosslinkers from the pointed end to the barbed end. We assume that every other subunit of each filament is crosslinked, such that the i-th crosslinker has an angle of with its attached actin subunit. The torque generated by the i-th crosslinker on the filament thus is −κ tor θ i . Here the crosslinker label i is ordered according to their distance to the pointed end of filaments. At torque balanced state, i.e., P n i¼1 y i ¼ 0, we have θ 1 = −(n − 1)π/13. Upon detachment of crosslinker 1, the total torque becomes imbalanced and the filament makes a rotation of angle Δϕ to reach a new torque balanced state, i.e., P n i¼2 ðy i þ D0Þ ¼ 0. This leads to Δϕ = −π/13. Similarly we can show that attachment of a new crosslinker at the (n + 1)-th position θ n+1 = θ 1 + n2π/13 will cause the filament to rotate the same angle in the same direction as caused by detachment of the first crosslinker.
Note that even though the rotation angle Δϕ is independent of the number of attached crosslinkers, it is required that crosslinkers are present in large enough number or are stiff enough to ensure that rotation will be significantly larger than thermal fluctuations, i.e. The above calculation can be easily extended to situations with more relaxed conditions than Eq (22). The angular rotation Δϕ (i) of the filament upon detachment of the i-th crosslinker satisfies the recursive relation: Directed rotation requires that Δϕ (i) have the same sign for all i. This condition is equivalent to considering that the angles between consecutive crosslinkers (θ i − θ i−1 ) have the same sign for all i.

Derivation of the effective binding energy μ eff
We consider a simplified model in which, when there are n crosslinkers formed, each crosslinker stores an elastic energy of E n . In fact, E n is varied among different crosslinkers and dependent on their positions and orientations in the filaments. Here we assume E n only depends on n, and the probability distribution P(n, t) for the number of crosslinkers is governed by where F(n|n − 1) denotes the number of ways to build the n-th crosslinker, given there are already n − 1 crosslinkers formed. The steady state distribution reads P ss ðnÞ ¼ Pð0Þ where GðnÞ ¼ 1 n! Q n i¼1 Fðiji À 1Þ is the number of possible ways to build n crosslinkers. By comparing the distribution of P ss (n + 1) with P ss (n), we have P ss ðn þ 1Þ P ss ðnÞ ¼ e À E nþ1 À m eff ðnÞ where μ eff (n) is the effective binding energy defined in the text. If μ eff (n) > E n+1 , the system is driven up to the (n + 1)-state, storing an elastic energy of E n+1 .