Remarkable structural transformations of actin bundles are driven by their initial polarity, motor activity, crosslinking, and filament treadmilling

Bundled actin structures play a key role in maintaining cellular shape, in aiding force transmission to and from extracellular substrates, and in affecting cellular motility. Recent studies have also brought to light new details on stress generation, force transmission and contractility of actin bundles. In this work, we are primarily interested in the question of what determines the stability of actin bundles and what network geometries do unstable bundles eventually transition to. To address this problem, we used the MEDYAN mechano-chemical force field, modeling several micron-long actin bundles in 3D, while accounting for a comprehensive set of chemical, mechanical and transport processes. We developed a hierarchical clustering algorithm for classification of the different long time scale morphologies in our study. Our main finding is that initially unipolar bundles are significantly more stable compared with an apolar initial configuration. Filaments within the latter bundles slide easily with respect to each other due to myosin activity, producing a loose network that can be subsequently severely distorted. At high myosin concentrations, a morphological transition to aster-like geometries was observed. We also investigated how actin treadmilling rates influence bundle dynamics, and found that enhanced treadmilling leads to network fragmentation and disintegration, while this process is opposed by myosin and crosslinking activities. Interestingly, treadmilling bundles with an initial apolar geometry eventually evolve to a whole gamut of network morphologies based on relative positions of filament ends, such as sarcomere-like organization. We found that apolar bundles show a remarkable sensitivity to environmental conditions, which may be important in enabling rapid cytoskeletal structural reorganization and adaptation in response to intracellular and extracellular cues.


Introduction
Understanding emergent behaviors in the actin cytoskeleton is important, since key biological functions such as cellular growth, division, and motility depend on cytoskeletal dynamics. Actin networks are transient and malleable within a single cell, constantly forming and remodeling various micro-architectures at different cellular locations. One salient class of such actin structures are actin bundles [1,2], which appear both in the cell body as stress fibers [3][4][5], or in specialized cellular processes such as filopodia [6][7][8], stereocilia [9][10][11] and microvilli [11][12][13]. The formation and functionality of bundles is spatially and temporally regulated by various proteins that interact with actin, including nucleation factors [5,14], crosslinkers [15,16] and molecular motors [17]. Depending on the cellular context and specific influence of actin regulatory proteins, different bundle structures emerge [18][19][20], with distinct functional roles as elaborated next.
First, owing to the polarity of individual actin filaments, an actin bundle may be unipolar, apolar, sarcomeric or of graded polarity (S1 Fig in S1 Text) [1], being tethered at either one or both ends via a cross-linker or a molecular motor. Cellular structures of unipolar bundles are ubiquitous across living cells: they aid in cytoplasmic streaming in pollen tubes [21][22][23] as well as in fungal hypha formation [24]. Unipolar bundles are found in animal cells in filopodia [25] and proximal dorsal stress fibers [5,26]. These structures are usually crosslinked but do not contain myosin as active crosslinkers in vivo [27]. However, there is evidence for myosin minifilaments, which are polymers of 10-30 myosin units [28], walking through tracks of unipolar bundles [29] under in vitro conditions. Non-muscle myosin II-A isoform minifilaments (NMII-A) [30] incorporate into dorsal stress fibers when these stress fibers are connected to transverse arcs [27]. Recently, NMI has also been implicated in filopodial force generation in neural growth cones [31]. Other actin bundle organizations, such as apolar bundles, underlie the portions of stress fibers near the cell center in fibroblasts [32], in sections of contractile rings [33,34] and in sections of ventral and transverse stress fibers. Composite apolar bundles are mostly found as two unipolar precursor bundles interacting with one another due to filament sliding in response to NMII-A activity [29]. Finally, more finely organized polarity arrangements are found in sarcomeric ordering, which has great significance in generating contractility of stress fibers [1,2].
The large variety of actin fibers and their distinct functional roles have instigated a growing interest in how various factors, such as crosslinking activity and minifilament concentrations, affect the dynamics and stability of bundles having unipolar and apolar organization. Another salient property of many actin bundles namely, non-sarcomeric contractility [35][36][37], has also been a focal topic of experimental and computational studies. An agent based model was employed by Zemel et al. to study one dimensional unipolar and apolar bundles of 10μm length under different concentrations of motors moving according to prescribed force-velocity relations [38]. They found that apolar bundles readily undergo sliding motions resulting in internal sorting. However, the lack of spatial details in their model did not allow predictions of the stability of bundle morphology under the influence of myosin. Dasanayake et al. [39] simulated actin bundles of 10μm length, modeled in 2D, considering explicit filament stretching and bending in addition to minifilament stretching forces in a confined volume. They found that apolar bundles exert higher wall stresses than parallel filaments.
Nevertheless, it is still unclear how actin bundles behave at long time scales under a wide range of crosslinker and myosin concentrations, especially if they were to be modeled in three spatial dimensions. The latter point should be emphasized because studying actin bundles in 3D is crucial for reconciling with the phenomenology of bundles observed under in vivo cellular conditions. Recent studies in three dimensions by Kim and coworkers indicated that filament polarity plays an important role in bundle formation from disordered actin networks [40] and also in tension generation [40,41]. To shed further light on the mechanochemistry of actin bundles, in this work we set out to establish the fundamental principles that govern the stability of two fundamental bundle organizations, namely, purely unipolar and apolar bundles.
Understanding conditions that either stabilize or destabilize various bundle geometries will shed light on the basic principles that govern cytoskeletal organization, bringing insights into complex in vivo processes, such as cytoplasmic disassembly [42] and actin network turnover [43][44][45]. For example, stability of untethered bundles is likely to depend on their internal polarity structure as well as crosslinker and myosin motor conditions. Furthermore, cells modulate actin filament turnover through a host of mechanisms such as filament severing [46][47][48], branching [1,49] and capping [50,51]. As a result, in vivo turnover rates are orders of magnitude faster than those commonly observed in in vitro experiments [47,[52][53][54]. In general, actin turnover plays a critical role in stress relaxation of entangled actin networks at longer timescales. For example, skeletal myosin is known to fluidize actin networks at very low mole ratios (myosin head to total actin mole ratio, M:A 0.0039) [55]. On the other hand, myosin's catch bond behavior can lead to long residence times resulting in actin networks maintaining large internal stresses at high myosin concentrations. In addition, transient passive crosslinkers arrest network configurations and prevent relaxation, however, crosslinker unbinding events allow for slow reconfiguration dynamics [56,57]. The crosstalk between filament turnover and crosslinker mechanokinetics has not been comprehensively studied under a wide range of conditions. Hence, it is also necessary to explore stability of bundles under a broad set of treadmilling rates.
To address the above outlined problems, we have used MEDYAN (MEchanochemical Dynamics of Active Networks) mechanochemical force field [58] to study the stability and dynamics of untethered actin bundles under a diverse set of polarity arrangements and levels of α-actinin crosslinker, myosin minifilament and treadmilling conditions. MEDYAN is further elaborated in the Methods section below. Our main finding is that unipolar bundles preserve bundle morphology at a wider range of crosslinker and, myosin concentrations than apolar bundles, because the latter experience a morphological instability due to being more susceptible to myosin induced intra-bundle filament sliding and shearing. We also show that three salient microarchitectures eventually emerge when simulating various untethered bundles under different conditions-bundles, asters and bundle-aster hybrids. Asters are filamentous networks exhibiting radial polarity sorting, with barbed ends clustered towards the aster center. We also investigated how the resulting phase diagrams depend on the speed of treadmilling. We found that network catastrophes, characterized by poorly crosslinked low density networks with obscure morphologies, occur when crosslinking cannot keep up with filament extension. Overall, our studies demonstrate that while a stand-alone apolar bundle is stable under a significantly narrower set of conditions compared to unipolar bundles, this geometric arrangement can serve as an important precursor to rich network remodeling phenomena such as global polarity sorting and sarcomeric organization.

Methods
MEDYAN is a mechanochemical force field for simulating active matter, including cytoskeletal networks. It deeply integrates chemical and transport dynamics with network mechanics, treating these phenomena on equal footing. MEDYAN has emerged from earlier efforts to model actin bundle growth in filopodia [59][60][61][62][63], where both active and passive transport were shown to critically influence growth dynamics [60,62,63]. MEDYAN's time evolution is based on alternating reaction-diffusion and mechanical equilibration steps, where the former events are stochastically generated according to the next reaction method [64], while the conjugate gradient approach is used to achieve mechanical equilibration [58]. This propagation scheme takes advantage of the wide timescale separation between slow chemical processes and fast mechanical equilibration speeds within sub-micron length scale volumes containing a portion of an actin network [65]. Furthermore, MEDYAN can effectively model actin filament polymerization processes as well as explicit α-actinin (crosslinker), myosin minifilament (motor) binding and unbinding events, in addition to myosin walking and mechanochemical feedbacks such as catch and slip bond behaviors.
In contrast to MEDYAN, other cytoskeletal modeling strategies, such as Cytosim (Cytoskeletal Simulation) [66,67], AFINES (Active Filament Network Simulation) [68], and the model by Kim and coworkers [41] rely on the Langevin dynamics of cytoskeletal components, explicitly simulating thermal undulations at the expense of significant diminution of computational efficiency. Among these models, MEDYAN's treatment of general reaction-diffusion processes is most comprehensive (in particular, with regard to G-actin's diffusion and reactions). In addition, MEDYAN's mechanical potentials are the most elaborate, for example, having an analytical excluded-volume potential representing steric repulsions and also complex dihedral angle potentials at the dendritic actin network branch points. A detailed comparison of different cytoskeletal modeling approaches can be found in Popov et al. (in particular, see Table S1 of [58]).
In MEDYAN, the reaction volume is divided into voxels based on the Kuramoto length [58,69]. Diffusing molecules are assumed to be uniformly mixed within each voxel and can diffuse from one voxel to another. Actin filaments can polymerize and depolymerize from both plus and minus ends based on experimentally determined rate constants [70]. Filaments that polymerize towards the boundary experience a reduction in polymerization rate based on the Brownian ratchet model [71]. In MEDYAN, the growth propensity for an actin filamentous tip is based on the local, instantaneous concentration of diffusing G-actin in the tip's neighborhood [58], while, in comparison, actin filaments grow/shrink at a constant rate in Cytosim, [72]. Furthermore, MEDYAN can explicitly account for ATP, ADP.Pi, and ADP bound actin monomeric states, enabling more elaborate simulations of F-actin polymerization dynamics [73,74]. In summary, filament length fluctuations and filament treadmilling can be studied using MEDYAN at the resolution of a single actin monomer. As elaborated below, we found that these fluctuations play an important role in determining whether the bundle stays coherent or undergoes a morphological transformation.
In MEDYAN, mechanical modeling of actin filaments is based on cylinder units with equilibrium spacing l m 0 � l p ðl m 0 = 108 nm in this study, lp-persistence length) connected with neighboring cylinders at flexible hinges. Persistence length of actin was reported from experiments as 17μm [75]. The MEDYAN force-field prevents filaments passing through each other via a novel cylinder-cylinder repulsion potential that is analytical, in contradistinction to the more widely employed technique of other comparable force fields, which relies on finding the closest distance between two cylinders that is used to compute their mutual repulsion [58]. Various MEDYAN mechanical potentials, such as intra-cylinder stretching and inter-cylinder bending are shown in Fig 1. α-actinin and myosin molecules that are bound to actin filaments are modeled as springs connecting two actin filament sites within their respective binding distances (α-actinin binding distance is in the range 30-40 nm, and minifilament binding distance is in the range 175-225 nm).
Those reactions that are mechanosensitive, for example, unbinding kinetics of actin binding proteins such as crosslinkers and motors, are influenced by the local instantaneous stresses of the actin network, via corresponding modifications of the reaction rate constants. Motor/ crosslinker binding and unbinding are modeled as a single step chemical reaction in MED-YAN, while in other force fields, for example, Cytosim, these processes occur via two elementary steps, comprising of separate reactions for each end of the motor or crosslinker [76]. In MEDYAN, crosslinker unbinding kinetics is modeled as a slip bond while myosin unbinding kinetics is modeled as a catch bond based on the parallel cluster model [77]. The motor walking rate is given by a linear force-velocity relationship for motor walking events. Further aspects of mechanical equilibration and mechanochemical feedback loops in MEDYAN as well as the implementation details of the chemical model are provided in Supporting Information (Section 2 in S1 Text).
Initial structures for all our simulations were based on 2 μm long actin bundles, comprising 30 actin filaments that correspond to in vivo stress fibers [32], where the internal arrangements of filament polarities were in either unipolar or apolar geometries. Filaments were initially placed on a hexagonal lattice with a spacing of 35 nm as found in experiments [78,79]. Bundles were modeled with α-actinin and NMIIA minifilaments that can bind and unbind from actin filaments. Bundles were simulated at 7 different concentrations of α-actinin (α-actinin to total actin mole ratio referred to henceforth as α:A 0.01, 0.05, 0.1, 0.25, 0.4, 0.6, 0.8) and 6 different concentrations of myosin (myosin head to total actin mole ratio referred to henceforth as M:A 0.0225, 0.045, 0.09, 0.18, 0.225, 0.675). 8 trajectories, each 2000 seconds long, were generated for each of the 6x7 = 42 mole ratio pairs (α:A, M:A) studied.
In addition, we also studied how the speed of filament treadmilling influences bundle stability. Myosin mole ratios of 0.0225 0.045, 0.225 and 0.675 were considered at α:A 0.01, 0.1 and 0.4 to investigate how all observed non-treadmilling morphologies behave under different treadmilling conditions. Treadmilling speed was varied by simultaneously altering polymerization and depolymerization rates at both filament ends by the same factor χ. As a reference, χ = 1.0 corresponds to the in-vitro treadmilling rate. We chose the following χ values, (0.1, 0.3, 0.6, 1.0, 3.0, 6.0, and 10.0), hence mimicking treadmilling speeds that are both slower and faster than the in vitro rate. As a technical detail, in this work we have developed a flexible volume protocol that permits expanding and contracting the reaction volume along the X-axis, which allows avoiding artificial boundary effects on the bundle major axis in a computationally efficient manner. This technique is further elaborated in Supporting Information (Section 2.5 in S1 Text). 7 trajectories were generated for each combination of the myosin mole ratio, α-actinin mole ratio and χ factor.
In order to understand the underlying morphologies sampled, we devised a hierarchical clustering scheme. To carry out structural clustering of obtained bundle configurations, we first computed the distributions of plus end-plus end (Dis ++ ), minus end-minus end (Dis -), and plus end-minus end (Dis +-) Euclidean distances. Jensen Shannon divergences [80] between each of the 42 mole ratios taken pair-wise were used to construct initial -condition -specific dissimilarity matrices (S12 Fig in S1 Text). The complete linkage method [81] results in a hierarchical cluster, which we visualized as dendrograms (shown below). Supporting Methods Section 2.4.1 in S1 Text provides further details on the clustering algorithm.

Long timescale morphologies of actomyosin bundles are primarily determined by their initial polarity and the myosin concentration
A wide array of steady state network morphologies emerge when non-treadmilling bundles with unipolar initial organization (non-treadmilling-BUInit) and apolar initial organization (non-treadmilling-BAInit) were simulated under a broad set of α-actinin and myosin concentrations (Figs 2 and 3 and S1 and S2 Movies). A trajectory is considered to have reached steady state when the network's radius of gyration reaches a stationary value (see S10 Fig in S1 Text). We found that for any combination of α-actinin and myosin ratios, marked differences are observed between steady state network morphologies of apolar and unipolar bundles. More specifically, antiparallel filaments show a strong tendency to mutually slide in response to myosin activity as a consequence of the latter being unidirectional walkers (S2 Movie). To quantitatively characterize the resulting network morphologies, we applied a novel structurebased clustering analyses that we have developed in this work (see the Methods section and the Supporting Information, Section 2.4.1 in S1 Text), which revealed dominant network morphologies preferred under various conditions, as elaborated below.
The resulting dendrograms reveal three broad clusters closest to the root, colored in green, red and blue, pointing to three major network morphologies (Figs 2-4 3 in S1 Text). Both order parameters delineated well the bundled morphologies from the aster-like morphologies. However, the hierarchical classification strategy introduced in this work was also able to identify the intermediate bundle-aster morphologies.
Some reflection on the internal structure of these dendrograms shows that within each initial polarity arrangement, NMIIA concentration is the main driver of the resulting network morphology (see Fig 4). Specifically, the same highest level clusters are formed from configurations with similar M:A ratios, while finer-grained additional clustering is determined by other factors, such as the crosslinker (α-actinin) concentration. The effect of myosin is primary because at high motor concentrations inter-filament distances become significantly widened outside of the α-actinin binding compatibility zone of 30-40 nm. This, in turn, crucially depletes the network of crosslinker binding sites (S4 Fig in S1 Text), hence, greatly diminishing parallel alignment among actin filaments and, subsequently, destabilizing the bundle phase.
It is interesting to compare our finding of the depletion of crosslinkers with increasing myosin concentration with the recent simulations that studied sorting of two different crosslinkers along an actin bundle. Freedman et al. have established that a significant length difference between two actin binding proteins, when combined with a specific range of the filament bending moduli, can result in spatial sorting of crosslinkers along a bundle [83]. Here, we found that active myosin walking may significantly increase overall inter-filament separation, thereby reducing the number of sites available for crosslinker binding in the bundles studied   When simulations were started with unipolar initial conditions, bundle-like states were stable when M:A � 0.09 (21/42 cases = 50% of all M:A/α:A mole ratios studied). On the other hand, apolar initial arrangements result in stable bundle-like states only in~20% (7/42) of the cases (i.e. when M:A = 0.0225). Thus, unipolar bundles are stable under a significantly wider range of conditions than their apolar counterparts. We tracked this difference primarily to myosin activity in apolar bundles giving rise to two thin polarity sorted sub-bundles that are together much longer than the initial bundle length and, furthermore, mutually interact via their barbed ends (see Figs 3 and 4). The resulting thinner bundles are more susceptible to bending deformations than the corresponding unipolar bundle. Consequently, at the time scale of about 30 minutes probed in this work, unipolar and apolar bundles arrive at different metastable morphologies despite being under the same crosslinker and myosin conditions. Presumably, if these systems were to be ergodic, then bundles are expected to eventually evolve to identical steady state configurations regardless of the initial polarity arrangement. However, as shown in our previous works [84,85], cytoskeletal dynamics maybe sufficiently glassy such that only metastable states are reachable over laboratory timescale.
Under very high myosin activity (M:A = 0.675), both unipolar and apolar bundles undergo a morphological collapse, preferring radially symmetric aster-like structures (the last column in Figs 2 and 3). We note in passing that in cells, asters are primarily found in microtubule networks as radially symmetric structures with filament plus ends spatially clustered together [86,87]. Actin networks are also expected to form radially polarity sorted asters [88]. Such structures have been found in in vitro treadmilling actin networks subject to skeletal muscle myosin (M:A 0.1), fascin and myosin [89] or just skeletal muscle myosin (M:A 0.02) alone [90,91], with the filament plus ends oriented towards the center of the aster.
The difference in threshold myosin ratios between skeletal muscle myosin and NMIIA for the onset of aster-like structures is partly explained by the increased processivity of muscle myosin. Actin asters were observed in-vitro when cytoskeletal structures were destabilized using Cytochalasin D [42]. Recently, they were shown to be essential in fission yeast cells during fusion [92,93]. Fission cell actin asters are considered to be formed due to Fus1 nucleators and multimerization of Myo51 and Myo52 [93]. Overall, we found that experimentally observed salient network morphologies of treadmilling networks can also be sampled in our simulations under non-treadmilling conditions.

Initial polarity arrangement plays a key role in the evolution of the treadmilling bundle morphology
Next, we investigated the combined effect due to network turnover and mechanokinetics (crosslinker, and NMIIA) by including actin filament polymerization and depolymerization processes in bundle simulations analogous to the systems discussed above. In the treadmilling study, steady state was defined when the average filament length fluctuations reach their stationary values (see S11 Fig in S1 Text). To systematically modulate filament treadmilling, we varied polymerization and depolymerization rates at both barbed and pointed ends by a factor χ between 0.1 and 10.0, where lower χ values cause slower treadmilling. We modeled a reaction volume having dimensions of 4 μm x 1.5 μm x 1.5 μm, with an initial total actin concentration of 5 μM, where the simulation box can expand and contract along the major axis based on the instantaneous bundle length (more details are provided in the Supplementary Information, Section 2.5).
We investigated both unipolar and apolar bundles at M:A mole ratios of 0.0225, 0.09, 0.225 and 0.675, whereas α:A was sampled at 0.01, 0.1 and 0.4. These values were carefully chosen to capture the different salient network morphologies manifested under non-treadmilling conditions (non-treadmilling-BUInit and non-treadmilling-BAInit), namely, BL, AL and ABI states. 7 trajectories, each 2000 seconds long, were generated for each of the 84 triad conditions (α:A, M:A, χ), starting from either unipolar (χ-BUInit) or apolar (χ-BAInit) initial conditions. Similar to non-treadmilling bundles, the above-mentioned classes of network morphologies were determined by clustering trajectories from 84 triads using the same clustering protocol (S12 Fig in S1 text). We classify connected networks with morphologies that do not belong to BL, AL or ABI as either Type A catastrophes if the networks are fragmented into smaller clusters or Type B catastrophes, where filaments in the network are poorly connected (S5-S9 Figs in S1 Text).
We found that these catastrophes emerge from the interplay between inter-filament connectivity and treadmilling. Treadmilling is characterized by a net depolymerization at minus ends and filament growth at the plus ends. In particular, increasing χ leads to a faster rate of filament growth at plus ends (and faster depolymerization at minus ends). In these systems, network stability is assured as long as the newly formed filament segments are effectively crosslinked. Under conditions where the latter does not take place, filament treadmilling dominates the system's structural evolution, leading to poorly connected networks. When well-structured initial bundle configurations result in a highly fragmented network, we denote such transitions as catastrophes (see S5 Fig in S1 Text). On the other hand, at higher mole ratios of α-actinin or myosin, the rate of inter filament connections is enhanced, which prevents network catastrophes.
Having established how inter-filament connections influence network stability, we next investigated the effect of treadmilling speed. Our clustering analysis indicates that treadmilling networks starting from unipolar initial configurations attain BL and ABI morphologies (i.e. non-aster, non-fragmented morphologies) in 33/84 (~40%) cases, compared to 19/84 (~21%) cases for the apolar initial configurations. On the other hand, the apolar initial arrangements lead to aster-like structures in 39/84 (46%) cases as opposed to 20/84 (24%) cases for the unipolar cases. Taken together, these results demonstrate that treadmilling bundles that were evolved from unipolar initial configurations are less likely to undergo morphological collapse than those evolved from apolar initial configurations.
Based on the network morphologies observed at the sampled triad combinations, we suggest the following two phase diagrams, shown in Fig 5, indicating dominant network morphologies as functions of M:A and χ. To justify the choice of these two order parameters, we point out that α-actinin dynamics determine bundle stability in conjunction with myosin, however, affecting only very weakly the final network morphology of stable networks. Fig 5A suggests that treadmilling unipolar networks sample similar network morphologies to the non-treadmilling cases for χ �1.0. At larger χ values, network treadmilling dominates as crosslinkers and minifilaments cannot effectively connect filament segments that are formed, leading to network catastrophes (S6-S9 Figs in S1 Text, see snapshots from simulations with χ >1.0) On the other hand, treadmilling apolar bundles lead to aster-like states (under non-catastrophic triads). The BL, ABI morphologies sampled by treadmilling apolar bundles have rich diversity due to the interplay between myosin activity and treadmilling as explained later.
We finally discuss the weak influence of crosslinkers on the morphology of a treadmilling network. At α:A = 0.01 and under low myosin mole ratios (0.0225, and 0.09 shown in Figures  S6 and S7), we found that unipolar bundles treadmilling at χ = 1.0 result in network catastrophes. However, significantly increasing crosslinker concentration (e.g. between 10 and 40 folds) can rescue such networks to prefer bundle-like/intermediate morphologies (S3 and S4 Movies). These findings are in qualitative agreement with the observations by Bidone et al [40] that at high crosslinker concentrations bundles form robustly from networks obtained from a wide range of initial orientational biases (at M:A = 0.08).

Treadmilling apolar bundles attain a diverse set of network morphologies at long time scales
Under low myosin activity conditions, treadmilling bundles starting from an apolar arrangement generate a remarkably diverse set of final network morphologies, primarily based on how filament barbed ends are spatially localized (see Fig 6). These emergent network geometries arise from the tug-of-war between myosin activity and filament treadmilling. On the one hand, treadmilling apolar bundles are subject to the continuous process of plus end extension and minus end retraction. On the other hand, myosin activity drives mutual sliding of neighboring filaments. If the rate of plus end extension is slower than myosin sliding, myosin activity dominates, leading to two unipolar bundles connected at their plus ends (Fig 6A) or to networks with overlapping plus end segments (Fig 6B). Conversely, under vigorous plus end extension conditions compared with myosin sliding, the network transitions to polarity sorted bundles interacting at their minus ends via myosin (Fig 6C and S5 Movie).
Modulation of χ also controls the overall distribution of myosin minifilaments in such bundles. The analysis of the computed spatial distributions of myosin minifilaments and α-actinin under low myosin concentrations (Fig 6D), indicates that myosin spatially segregates close to pointed ends, characteristic of sarcomeric ordering [18, 36,94]. The latter arises when myosin minifilaments interact with minus ends flanked on either side by plus ends.

Discussion
Actin bundles are important for cellular stability, growth and mechanosensing. While prior experimental and modeling research has primarily focused on bundle formation processes [40,72,95,96], in this work we have addressed the stability and temporal evolution of various bundle configurations. We used MEDYAN, a mechano-chemical forcefield based on molecular principles, to simulate bundle dynamics in 3D. The dimensionality of the model is crucial as filament deformations and mutual interactions are markedly dimension-dependent. In this study, we comprehensively analyzed how α-actinin and myosin influence the stability and morphological transformations of unipolar and apolar actin bundles. We discovered that at time scales of about 2000 seconds, non-treadmilling unipolar bundles are stable under a wider range of crosslinker and myosin mole ratios compared to apolar bundles. At high myosin mole ratios, we observed aster-like states characterized by interacting barbed ends grouped in the center of the cluster with radially emanating pointed ends.
We also investigated how actin turnover affects bundle morphology fates, by developing and applying a simulation protocol that allows moving system boundaries. Our results indicate that treadmilling bundles, both unipolar and apolar, undergo network catastrophes when the network's ability to form inter-filament connections is insufficient compared to the treadmilling speeds. In vivo cytoskeletal networks undergoing fast treadmilling may be able to avoid such undesired fragmentation using additional mechanisms such as filament capping and actin filament nucleators.
Interestingly, at high myosin concentrations, even quick treadmilling does not rescue the network from transitioning to aster-like configurations. On the other hand, at low myosin activity, initially apolar bundles explore a diverse spectrum of network organizations, which we attributed to the tug-of-war between minifilament activity and filament treadmilling. Under certain conditions, interesting, biologically relevant architectures emerge, such as sarcomeric-like organization. We are not aware of prior models that resulted in the spontaneous assembly of sarcomeric arrangements without imposing spatial restrictions on crosslinkers. In particular, previous attempts to reproduce sarcomeric distribution of treadmilling apolar networks relied on various assumptions, such as preferential binding of passive crosslinkers near plus ends [18] or considered systems with both plus and minus end directed motors [97].
Finally, we reflect on the biological consequences that follow from this work. We found that treadmilling bundles with apolar initial configuration are poised to undergo a remarkable morphological response to the perturbations of the environment, such as alterations of myosin activity or treadmilling rates. On the one hand, this level of sensitivity to parameters might be potentially detrimental to the overall stability of the cellular actin network. On the other hand, if only a small fraction of the cytoskeleton is organized as apolar bundles, the latter can sensitively respond to various signaling cues that affect the local concentrations of actin binding proteins. Thus, we propose that apolar bundles might be crucial to cytoskeletal robustness and adaptation in scenarios that demand drastic structural reorganization. The ability to rapidly change network morphology might be important in certain cellular functions where force production or rapid cellular reorganization are necessary. Overall, the optimal choice of bundle architecture should be determined by the specific cellular processes that it affects: for example, in the case of cargo transport [98] or protrusive growth [99], a structurally stable unipolar bundle would be the optimal choice, while contractile elements that require frequent reorganization would be more robust when apolar bundles are incorporated into their architectures [5,100].
In summary, we simulated actin bundles in 3D, while explicitly accounting for excluded volume interactions, diffusion of actin, crosslinker and NMIIA proteins, and numerous chemical and mechanical processes that enhance the model's realism. We systematically studied the influences of initial bundle polarity, concentrations of myosin and α-actinin and the network turnover rate, finding a remarkably rich palette of bundle evolution trajectories, from stable bundle states to asters and sarcomeric organizations. In future works, additional effects may be considered, such as actin filaments transiently tethering to the substrate and nucleation of filaments via formins or Arp2/3, which will bring us closer to achieving a more complete understanding of bundle dynamics under in vivo conditions.