Simulated actin reorganization mediated by motor proteins

Cortical actin networks are highly dynamic and play critical roles in shaping the mechanical properties of cells. The actin cytoskeleton undergoes significant reorganization in many different contexts, including during directed cell migration and over the course of the cell cycle, when cortical actin can transition between different configurations such as open patched meshworks, homogeneous distributions, and aligned bundles. Several types of myosin motor proteins, characterized by different kinetic parameters, have been involved in this reorganization of actin filaments. Given the limitations in studying the interactions of actin with myosin in vivo, we propose stochastic agent-based models and develop a set of data analysis measures to assess how myosin motor proteins mediate various actin organizations. In particular, we identify individual motor parameters, such as motor binding rate and step size, that generate actin networks with different levels of contractility and different patterns of myosin motor localization, which have previously been observed experimentally. In simulations where two motor populations with distinct kinetic parameters interact with the same actin network, we find that motors may act in a complementary way, by tuning the actin network organization, or in an antagonistic way, where one motor emerges as dominant. This modeling and data analysis framework also uncovers parameter regimes where spatial segregation between motor populations is achieved. By allowing for changes in kinetic rates during the actin-myosin dynamic simulations, our work suggests that certain actin-myosin organizations may require additional regulation beyond mediation by motor proteins in order to reconfigure the cytoskeleton network on experimentally-observed timescales.


Introduction
Virtually all cells contain a cytoskeleton, a collection of structural filaments that are required for critical processes including division and migration [1]. The cytoskeleton consists of three major classes of filaments: actin filaments, microtubules, and intermediate filaments. The actin cortex, a thin meshwork of actin filaments just below the cell membrane, is a major constituent of the cytoskeleton. Indeed, actin accounts for 10% or more of a cell's total protein, making it one of the most abundant proteins [2]. Actin filaments are highly dynamic, growing and shrinking through the gain and loss of individual actin monomers. Actin filaments are also polar, with distinct polymerization kinetics at the two ends resulting in directionally biased filament growth. The end that favours actin monomer addition is called the barbed (or plus) end, while the end that is less favourable for polymerization is called the pointed (or minus) end. Cells rely on the dynamic nature of actin filaments to respond quickly to internal and external cues by reorganizing the actin cortex. These actin reorganizations can result in shape changes and variations in the mechanical properties of the cell [3][4][5][6]. Despite our understanding of the complex filament-level dynamics of actin, conditions that favor formation of specific actin network architectures are still poorly understood, and thus the formation of cortical actin networks is the focus of this study.
In addition to polymerization dynamics, actin filaments are transported by the activity of motor proteins, particularly those from the myosin superfamily. Members of the myosin superfamily bind to actin filaments and hydrolyze ATP during their power stroke to generate force, resulting in movement of individual filaments. In many cells, including the early C. elegans embryo, Type II myosins, also called conventional myosins, are implicated in reorganization of cortical actin filaments. The family of Type II myosins, which consists of the proteins NMY-1 and NMY-2 in C. elegans [7], assemble into mini-filaments, with multiple heads containing actin binding domains at either end of the mini-filament. This bipolar mini-filament structure allows Type II myosins to simultaneously bind two actin filaments, moving the bound filaments relative to each other. Type II myosins are plus (barbed) end directed motor proteins, meaning that they take a "step" towards the plus end of the actin filament during their power stroke, resulting in movement of the actin filament in the direction of its plus end. These myosins are also non-processive, meaning they release from the actin filament after a single power stroke and do not continue to "walk" along the actin filament. The released myosin diffuses until it finds another pair of actin filaments available for binding. In cells, myosin motors often assemble into myosin mini-filaments, which consist of small ensembles of myosin heads and can have increased processivity [8,9]. Myosin can also be prevented from performing a power stroke if the force applied to a bound myosin mini-filament is greater than its stall force.
Eukaryotic cells contain approximately 40 different myosin genes [10], and many of these other myosin motor proteins are thought to be involved in actin cytoskeleton organization. For instance, myosin V is an unconventional myosin which transports cargo as it moves along actin filaments [11], and plays a critical role in fission yeast cytokinesis [12,13]. Myosin VI, which moves towards the minus end of an actin filament [10], segregates to distinct spatial locations throughout the cell cycle [14], and is critical for cell proliferation in certain cancers [15]. While the dynamics of these other myosin motors are less well understood compared to Type II myosins, it is clear that kinetic parameters associated with these different myosins can vary widely. For instance, non-muscle myosin IIA motors are thought to have an individual head step size of 6 nm during the power stroke [16,17], whereas the myosin V motor has been found to have a mean step size of 36 nm [18], and the processive myosin VI motor has a broader distribution of step sizes, with mean forward steps of � 30 nm [19]. Similarly, the characteristic unbinding force of the non-muscle myosin IIA motor is assumed to be 12.6 pN in [17,20], while muscle myosin II has been found to have a measured average unbinding force of 9.2 pN [21], and the myosin-V unbinding force has been estimated as 3-5 pN [22]. Given the wide variability in kinetic parameters associated with myosin motor proteins despite their highly conserved genomic sequences [23,24], and a lack of comprehensive information about the dynamics and role of different motor proteins in regulating actin cytoskeleton organization, we focus here on a detailed study of the effect of different parameters associated with motor protein activity on actin networks. Other regulating factors such as actin-binding proteins, nucleators, cross-linkers, or changes in turnover dynamics have also been implicated in actin organization [25][26][27]. Here, we aim to understand and quantify how motor regulation (through variation of individual motor kinetic rates) influences changes in actin cytoskeleton organization. This allows us to characterize changes in the actin cytoskeleton observed in vivo without constraining the simulations to a specific motor protein and its properties.
One specific example of striking cortical actin reorganization occurs over the course of the cell cycle. For instance, in the early embryo of the nematode worm Caenorhabditis elegans (C. elegans), as reported in [28,29], and more recently in [30], cortical actin filaments are initially organized in an open meshwork, characterized by patches with few filaments. This open meshwork is reconfigured into a homogeneous, isotropic mesh. As the early embryo prepares for first division, cortical actin filaments are aligned at the middle of the embryo to form the cytokinetic ring, with filaments outside the cytokinetic ring orienting towards the division plane. These actin reorganizations occur over the course of approximately 15 minutes in the early embryo. This substantial reorganization of cortical actin through the cell cycle is a common theme in many cell types, from plants [31][32][33] to mammals [34,35]. Despite variation in the specific filament organization (as shown schematically in Fig 1), the in vivo observations suggest that reorganization of cortical actin is a common feature of many different cell types in many different organisms.
A number of mathematical models have been proposed to investigate the formation of higher order actin structures due to the activity of myosin. Continuum models consisting of PDEs have shown the ability of motor proteins with different characteristics to reproduce experimentally observed structures [36][37][38]. This approach allows for analysis of the corresponding model, but does not take into account the noisy interactions of individual filaments or motors. In addition, these continuum mathematical models are unable to capture the structural evolution of the interacting proteins at the molecular level. Stochastic models that explicitly simulate individual filament and motor dynamics have been used to yield insights into the dynamics of actin and myosin structures. Reviews and comparisons of existing agentbased cytoskeletal models are provided in [17,39,40]. While the different models vary in their implementation, these frameworks consistently show that changes in motor protein activity can induce different actin organizations [40][41][42][43][44]. These modeling approaches have yielded insights into the range of actin-based structures that can be formed in the presence of motor proteins such as Type II myosins, but is it not yet understood how motor proteins can efficiently and robustly transition between different actin organizations or how they might coordinate to establish observed actomyosin structures.
In this investigation, we use the stochastic simulation platform MEDYAN to simulate the organization and transition between different actin organizations. This modeling framework accounts for complex chemical dynamics of various chemical species, a mechanical polymer model incorporating filament turnover dynamics, as well as mechanochemical coupling. Additional features of this agent-based model framework and details on how it improves on prior coarse-grained models of cytoskeletal dynamics are provided in [17, Table 1]. Recent models have focused on addressing questions such as the impact of crosslinker density and filament rigidity on local shape deformation [27,40], the contribution of the actin network on cortex tension regulation [45], or the impact of motor and actin concentration on structure formation [44]. MEDYAN is more comprehensive in that it includes features such as polymerization and depolymerization of semi-flexible actin filaments, steric interactions, a realistic model of myosin mini-filament ensembles, and allowing for a three-dimensional simulation domain. This framework allows us to predict realistic cytoskeletal behavior for a large set of myosin motor parameters, as well as for interactions of two motor types with different characteristics. We therefore adapt data analysis measures to characterize a variety of simulated actin structures and quantify the time course of their formation, focusing on understanding how one or two motor populations with different kinetic characteristics may regulate the dynamic behavior. We find that a single motor protein is capable of producing a range of actin structures, with variations in motor protein step size, binding rates, stall force, and number of motor heads resulting in the greatest changes in actomyosin organization. In particular, changes in these parameters can produce actin structures ranging from tightly clustered foci to loose meshworks of filaments, which have been observed in different experimental settings, including during different stages of the cell cycle. When two motor protein populations with different kinetic parameters interact with the same actin meshwork, additional properties emerge: actin structures may adopt an intermediate organization, between the two extremes of the motor proteins acting alone; one motor protein may dominate, entirely dictating the structure of the actin meshwork with the second motor protein acting as passive cargo; and we also observe some motor protein segregation, with motor proteins occupying distinct spatial regions. Additionally, we find that transitioning between between actin structures can be achieved in a more timely manner when two motor proteins act together. Together, these results demonstrate the importance of cooperation between motor proteins to efficiently construct and reorganize actomyosin systems.

Stochastic simulation framework for actin-myosin interactions
We carry out mechanochemical simulations of actin-myosin interactions using the MEDYAN (Mechanochemical Dynamics of Active Networks) modeling framework developed in [17]. This simulation package uses a coarse-grained representation of interacting semi-flexible polymers (actin filaments) in three dimensions. The cytoskeletal network mechanics are integrated with stochastic reaction-diffusion processes, whose dynamics are calculated using the next reaction method [46]. The simulation space is divided into compartments and diffusing molecules are assumed to be uniformly mixed within each compartment. Stochastic movement between compartments is used to model the diffusion and molecular transport of various chemical species.
Here we use MEDYAN to model actin filament polymerization phenomena, in addition to essential crosslinker (α-actinin) processes such as binding and unbinding. Additional active processes involving motor protein (such as myosin II minifilament) binding, unbinding, and walking are incorporated in the model. The force fields employed to model the actin filaments, as well as their interaction potentials with linkers and motors that characterize filament deformations, are detailed in [17]. To understand actin-myosin organization in the simulation domain, we take advantage of the mechanical modeling of actin filaments, which consist of cylinder units with equilibrium spacing [43]. Further details about the MEDYAN model framework and implementation can be found in [9,17,43,47].
We are especially interested in using this stochastic simulation framework to understand motor regulation and how it impacts cytoskeleton organization. While changes in other regulatory proteins or in turnover dynamics may also occur, here we focus on whether variation in parameters associated with motor protein activity can account for the diversity and dynamics of actin-based structures observed in vivo. In addition, we investigate the emergent actin-myosin organization in simulations where two populations of motors with distinct properties interact with actin filaments. This is motivated by studies into the maintenance of ring channels, circular openings in developing C. elegans oocytes, which suggest that the Type II myosins NMY-1 and NMY-2 act antagonistically to maintain a stable ring channel opening. Further, it was shown that NMY-1 and NMY-2 occupy spatially distinct regions near the ring channel opening [48]. Thus, we wish to investigate the conditions and kinetic parameters under which motor proteins may segregate into spatially distinct regions.

Actin dynamics and accessory protein settings
In our simulations, we consider a 2μm × 2μm × 0.2μm domain, with cubical compartments with side length of 0.2μm. This domain geometry reflects a region of the cell that is close to the membrane and thus essentially models a patch of the cell cortex where myosins are likely to have a significant impact on actin reorganization. A similar thin rectangular domain was used to simulate F-actin and myosin-II minifilament interactions at the cortex in [44] and was found necessary to reproduce some features of in vitro patterning.
As in our previous work [49], we carry out standard implementations of the model in [17], which is parameterized for actin filament polymerization and depolymerization, α-actinin cross-linking proteins, and non-muscle myosin IIa motor filaments; however, we use larger numbers of myosin motors, consistent with the myosin concentrations used in computational studies of actin bundles [43]. As in [17,20], myosin motors are modeled so that the mechanochemical effect of increased pulling is a catch bond; this is reflected by an exponential decrease in the motor unbinding rate with increase in the total stretching force experienced by the motor ( [17], Supporting Information S3 Text). A concentration of α-actinin crosslinking proteins is also included, and we refer to ( [17], Supporting Information S3 Text) for the mechanochemical model used for these linkers. Actin filaments are initialized as short polymers with random positions and orientations in the simulation domain. Actin turnover is incorporated in all simulations, with polymerization and depolymerization reactions at both ends of the filament and reaction constants chosen as in [17] and informed from experimental results in [50]. Excluded volume interactions between cylindrical units in neighboring polymers are included in all MEDYAN simulations [17].

Myosin motor parameters influence the emerging actin organization
We begin by considering the dynamic organization of actomyosin networks in the presence of one myosin motor population. In this manuscript, we refer to bipolar aggregates of myosin molecules as myosin motors. Such bipolar minifilaments stay bound longer and can work more efficiently compared to individual myosin molecules. To accurately model minifilaments based on the implicit properties of individual myosins, we use the parallel cluster model, which offers a rigorous statistical mechanics-based paradigm to understand the emergent behaviors of a group of myosins [20]. In this model, the kinetic parameters of binding, walking, stalling, and unbinding of myosins can be predicted based on individual myosin properties such as binding rate, unbinding rate, stall force, and unbinding force. In addition, we can also account for the variability in the number of myosins in a population of minifilaments.
In [17], motor parameters are chosen to model the behavior of non-muscle myosin IIA minifilaments. Since various myosin motors have been hypothesized to exert force on the actin cytoskeleton, we investigate the impact of different motor properties on cytoskeleton organization. We build on the simulation framework and baseline parameter values for nonmuscle myosin IIA motors in [17] (see Table 1) to uncover such differences in actomyosin organization. By changing one motor parameter at a time and characterizing the resulting organization using the data analysis methods described in § 5, we suggest potential mechanisms of motor regulation that may be responsible for changes in actin assembly. Since MED-YAN simulations are stochastic, the dynamic actomyosin organization may vary across simulations; unless otherwise noted, we consider ten independent stochastic runs for each parameter setting.

Baseline parameter simulations.
We begin by describing the actomyosin organization under the baseline parameter values in Table 1 using the data analysis methods in § 5. In these baseline simulations, the myosin motors dynamically organize the actin cytoskeleton into 1-2 clusters with some stray filaments; Fig 2A illustrates several time snapshots of a sample actin-myosin simulation (see S1 Video for a complete sample simulation). Actin filaments organize both into parallel bundles as well as into aster-like structures under baseline conditions. The actin organization can be described using the radial distribution function (see § 5.3), which is shown in Fig 2B for the same simulation times, averaged over ten stochastic MEDYAN runs. This measure illustrates that the inter-monomer distance distribution changes from a wide peak at medium distance (radius) values to a large peak at small values corresponding to the filaments that are clustering together, as well as a flatter peak at larger distances. This reflects that, through time, the actin cylinders making up the filaments can belong to different actin structures or clusters.
Additional methods of characterizing the actin cytoskeleton organization are provided in Fig 3A. The actomyosin network radius of gyration (described in § 5.1 and introduced for this system in [17]) shows the effect of stochasticity in the model simulations, as different runs exhibit a range of behaviors, with some simulations leading to a small increase in the radius of gyration (decreased network contractility) and others leading to a small decrease in the radius of gyration (increased network contractility). Since there is no global alignment of filaments in the simulation domain, the orientational order parameter (described in § 5.2 and introduced for this system in [17]) does not show significant changes through time. Finally, the actin organization shows an overall clustered distribution in the spatial statistics measure described in § 5.7 (as opposed to a regular or spatially random distribution of actin cylinders).
The myosin motor organization is visualized in Fig 3B using a three-dimensional motor localization plot as a function of time and of distance from each motor to the center of the domain (further described in § 5.4); this average radial motor localization with respect to the middle of the domain does not change significantly through time. However, the measure defined in § 5.5, which calculates the area of the boundary polygon around the myosin motors,  (Table 1). Baseline simulation snapshots at times (A1) t = 100, (A2) t = 1000, and (A3) t = 2000s, with actin filaments depicted as long red polymers, plus ends of filaments indicated in dark red, myosin motors represented as medium-length dashed blue lines, and cross-linkers shown as short black lines. The parameter values used are d step = 6 nm, k head.bind = 0.2/s, F s = 100 pN, range heads = 15-30 (see Table 1). shows a steady decline through time, indicating that the motors are overall localizing in space as they cluster actin filaments into tighter actomyosin structures. This is further confirmed using the spatial statistics-based measure described in § 5.7, which increases through time and therefore suggests that the distribution of the motor protein pattern becomes increasingly more clustered through time. This aggregation of myosin motors at the core of clustered actin structures is consistent with simulated [27,39,40,43,44], in vitro [27,44,[52][53][54][55], and in vivo [56] experimental observations on actin reorganization by myosin motors.
In the following sections, we present variations in motor parameters that lead to significant changes in cytoskeleton organization as compared to the baseline conditions. In Table 2, we summarize how these parameters affect microscale aspects of myosin minifilament behavior  Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 independent stochastic runs. See Fig 2 for  in the MEDYAN model as well as how they impact network-level filament and motor organization.

3.1.2
Step size. We refer to the physiological binding distance of a single myosin motor head d step as the motor step size. In the MEDYAN model, this parameter affects the base walking rate of the motors: k 0 1À r r k head;bind , where ρ is the motor duty ratio, d total is the distance between binding sites on the model actin cylinders, and k head,bind is the single head binding rate [17]. Fig 4A shows the actin-myosin organization at the final time of sample simulations with small (3 nm) and large (12 and 36 nm) myosin step sizes relative to the baseline value. The small step size leads to considerably more spread out filament organization, with some filament alignment at the domain boundaries, consistent with the motor localization in Fig 4B. On the other hand, the larger step sizes lead to more compact contractile actin-myosin clusters, with motors localized inside these clusters. This is due to the fact that increasing the step size leads to an increase in the base motor walking rate, so that myosin motors have better access to filaments and therefore lead to enhanced contraction of the network. We find that this behavior of the actin-myosin organization is similar across additional stochastic runs, as illustrated by the time-series measures in Fig 5. The actomyosin network radius of gyration and the area of the motor boundary both increase at small step sizes, reflecting the relaxing of the filaments into a more homogeneous distribution for small step sizes. For larger step sizes, the radius of gyration and the motor area decrease through time, showing faster establishment of contractile clusters.
The large stepsize simulations thus result in the formation of aster-like structures with radial polarity sorting, as the barbed (plus) ends of the actin filaments point primarily towards the aster centers. This is consistent with previous MEDYAN modeling results for myosin IIdriven actin organization [43], as well as with other computational models for myosin II and myosin VI which find that, in disordered bundles, filaments slide with each other due to motor and cross-linker forces and establish local polarity sorting [39,41,44]. In addition, aster-like filament assemblies where the barbed ends point towards the core of the structure have also been observed in simulations [39,43], in vitro [27,39,52,53,55] and in vivo [56] model systems.  Table 1). https://doi.org/10.1371/journal.pcbi.1010026.g004

Binding rate.
We denote the per-head motor binding rate by the on-rate k head,bind . Increasing this parameter leads to an increase in the base walking rate of the motors, but also affects the base filament unbinding rate in a nonlinear way according to the parallel cluster model for non-processive motors [20] used in MEDYAN [17]. The small binding rate sample simulation in Fig 6A shows a slightly more spread out cytoskeleton organization, whereas the large binding rate simulation (0.4/s) displays compact clusters with fewer free filaments than in the baseline case. In general, increasing the on-rate leads to an increase in the motor's duty ratio (the proportion of time that a head spends in the bound state) and therefore yields a larger number of bound heads, so that myosin motors reside on the filaments longer and contract them. However, further increasing this rate to 0.8/s leads to considerably looser and more spread out actomyosin organization, with no noticeable clustering (see S2 Video). The reason for this observation is that the motors reside on the filaments much longer and are generally less mobile, so that they play more of a cross-linking role in the dynamic actin organization.
These observed behaviors are consistent across simulations, as illustrated by the time-series measures in Fig 6B. Although there is more variability across model runs for the 0.8/s binding rate, the actin organization relaxes into a more homogeneous distribution in this parameter setting, while the myosin motors spread out across a larger portion of the domain. We note that, unlike the linear change in actomyosin behavior as a result of varying the step size d step in § 3.1.2, the system behaves nonlinearly as the head binding rate k head,bind increases; this is due to the fact that the latter parameter can impact multiple mechanisms in the model (base walking rate, base filament unbinding rate), thus providing additional and more nuanced insights on the impact of myosin motor parameters on the cytoskeleton organization. It is also worth noting that the actomyosin structures in Fig 6A are characterized by significant fractions of the actin filament barbed ends pointing outward towards the domain boundaries. This is similar to in vivo experimental observations in cultured fibroblasts [56], where actin filaments located at the periphery of control cells were found to be oriented with their barbed ends outwards. The radius of gyration decreases as motor step size is increased, reflecting the formation of actin clusters (Fig 4A). (A2) The area of the polygon enclosing the motor proteins also decreases as motor step size increases, due to motor proteins binding to clustered actin filaments. The parameter values used are k head,bind = 0.2/s, F s = 100 pN, range heads = 15-30 (see Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 stochastic runs. See Fig 3 for fil;walk is the base walking rate of the motors, F ext is the external force or tension experienced by the myosin filament, and α is a parameter that tunes the strength of the dependence on the external force [17]. Fig 7A1 shows that increasing the stall force is associated with an increase in clustering and contractility, since the walking rate stays larger for higher external forces experienced by the motor. While there is more variability in the stochastic runs associated with this parameter, Fig 7B1 is consistent with this observation  Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 stochastic runs indicating high variability between different runs. See Fig 3 for Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 stochastic runs. See Fig 3 for  that larger stall forces lead to a more contractile actin network and to tighter spatial segregation (on average) of the motors. The actin structures emerging from changes in the stall force parameter resemble patches and aster-like assemblies observed in previous computational [39,40,43,44] and experimental [26,27,54,57] studies investigating the actin reorganization by myosin II motors.
3.1.5 Number of motor heads. We let N t represent the number of heads in the myosin filament, and vary the range of this parameter. N t influences the base filament unbinding rate according to the parallel cluster model for non-processive motors [20] used in MEDYAN [17], and is directly proportional to the zero force residence time. Fig 7A2 shows that allowing for exactly 2 heads of the myosin minifilament (a dimer configuration) leads to a more spread out actin organization, with alignment at the domain boundaries, while larger numbers of heads (minimum 30 and maximum 45) yield more compact cluster organization. This observation is also summarized using the measures characterizing the dynamic actin and myosin organization in Fig 7B2. This observation is consistent with the fact that increasing the number of heads N t leads to a decrease in the base filament unbinding rate, so that myosin motors are less likely to unbind from actin and therefore more consistently organize actomyosin clusters. The myosin heads activity was previously found to be essential for actin organization, since inactivated myosin II motors were not able to produce patterns in vitro when studying interactions of actin with skeletal muscle myosin and the fascin bundling protein [53]. Our results on increased actin network contractility with an increase in myosin motor heads also agrees with in vitro studies in [54], which found that larger clusters consisting of myosin II molecules lead to a more dynamic and contractile network, while smaller motor clusters underwent less rearrangement and were more mechanically stable. As in [54], we similarly find that myosin motors are embedded within the network of actin filaments and contribute to the formation of actin structures, as illustrated in Fig 7A2. 3.1.6 Open meshwork organization. While some of the parameter settings investigated above illustrate an opening in the actin organization and some filament alignment at the boundaries, another means of generating an open actin-myosin meshwork is to reduce the number of myosin motor minifilaments in the MEDYAN model simulations. S3 Video shows the progression to a more open meshwork in the simulation domain with a decrease in the motor number as well as with a reduction in the motor stepsize. This suggests that regulation of the availability of the active motors may be one way to generate the open cytoskeletal meshworks observed in vivo, such as in early stages of the cell cycle. Lower motor concentrations have also previously been found to give rise to less reorganization and contractile structure formation than higher motor concentrations in simulations in [17,39], and to generate bundles of mechanically stable networks in in vitro studies of actomyosin networks in [52,54].

Two-motor populations contribute to tuning of cytoskeletal organization and reveal dominant motors
Experimental observations have shown that several types of myosin motors [48] or multiple populations of the same myosin motor with characteristics that depend on the local cellular environment [58] may regulate actomyosin organization. Motivated by these observations, we study the impact of two-motor populations with different motor parameters on the organization of actin cytoskeletal networks in MEDYAN stochastic simulations. We consider the same total number of motors as in § 3.1, divided into equal numbers of motors for each of the two motor populations of interest.

Tuning behavior of motor populations with different parameters.
In many of the model simulations performed, we find that the actomyosin organization is tuned so that various measures of cytoskeletal network behavior lie in-between those corresponding to the behaviors of interactions with a single motor population (i.e., characterized by a single motor parameter set). Two sample examples for interactions with motors with 3 vs. 6 nm step sizes as well as with dimers (exactly 2 heads at each end of the minifilament) vs. motors with 30 − 45 heads are shown in Fig 8; as in the case of simulations with one motor population, we represent actin filaments as red polymers and cross-linkers as short black lines, while here myosin motors with the first parameter value are shown in dashed blue lines and with the second parameter value, in green dashed lines. In both examples, the distribution of actin inter-monomer distance distribution at the final simulation time is balanced between the distributions resulting from simulations with each of the single motor populations. Similarly, the measure that quantifies myosin motor localization with time (calculated for all motors in the simulation) reflects the same observation that actomyosin behavior is tuned compared to the singlemotor population settings. As before, the resulting actin structures show previously-observed characteristics, such as the outward orientation of actin barbed ends at the periphery (Fig 8A2) or the formation of bundle-aster hybrids with both types of myosin motors embedded in the actin network [43,54,56].

Dominant behavior of certain motors.
In certain two-motor population simulations, we find that one of the motors dominates the dynamics and is able to re-position the other motor population. In these parameter settings, the dominant motors dictate the overall actin organization. To illustrate this, we build on the network contractility measure in § 5.1 to determine the first time in each simulation when the radius of gyration increases or decreases by a certain threshold amount (determined based on the relative increase or decrease in contractility observed for that parameter). The box plots in Fig 9A show Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 stochastic runs.

Motor segregation in certain parameter regimes.
In few of the parameter regimes investigated, we found evidence of some spatial segregation of the two motor populations interacting with actin filaments. In Fig 10, we use the measures described in § 5.6 and 5.7 to analyze the interactions of actin filaments with motors with 3 vs. 36 nm step sizes as well as with motors with on-rates of 0.1/s vs. 0.8/s. The measures described in § 5.6 rely on finding the two-dimensional boundary polygons around the point clouds consisting of each motor's center locations; we then compute the intersection area between these boundary polygons for the two motor species as well as the distance between the centroids corresponding to the two polygons. Fig 10A shows that the intersection area of the boundaries for motors with the two different step sizes decreases in time compared to baseline simulations, thus suggesting that the motors might segregate in space; however, the distance between the centroids of the motor boundaries does not change significantly. This is because the actomyosin network is consistently organized in a tight cluster for this motor combination, with motors in both categories localizing closer together through time. To further clarify the distribution pattern of the two motors, we apply the spatial statistics measure described in § 5.7 to each motor population. The right panel of Fig 10A shows that the 36 nm step size motor has an even tighter cluster distribution within this actomyosin network. Fig 10B suggests that there is distinct spatial segregation of the motors with 0.1 vs. 0.8/s binding rates, given that the normalized intersection  Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 independent stochastic runs.
https://doi.org/10.1371/journal.pcbi.1010026.g009 PLOS COMPUTATIONAL BIOLOGY area of their boundaries decreases and the centroid distance between these motor boundaries increases through time. We further confirm this by visualizing the spread measure for each motor population in the right panel: the small on-rate motor forms a cluster through time, whereas the large on-rate motor is distributed throughout the simulation domain given its less mobile behavior (as also observed in § 3.1.3).

Transitions in motor parameters reflect the remodeling ability of the cytoskeleton
To understand the capacity of the actomyosin network to re-organize under myosin motor regulation in response to stimuli, such as during the cell cycle or in cells where local ATP abundance is altered, we implement a MEDYAN framework where the myosin motor binding rate can change at a specified time point during the simulation. In particular, we consider the setting where the myosin binding rate switches between 0.8/s and 0.4/s. We chose these parameter values for our study since, as shown in Fig 6 in  § 3.1.3, the 0.8/s binding rate leads to loose actomyosin organization whereas the 0.4/s binding rate generates organizations with compact clusters.
In Fig 11, we explore changes in this motor parameter 4000 s into simulations that last a total of 9000 s, in order to allow the system to equilibrate. In § 3.1.3, we found that a relatively  Table 1). Solid and dashed lines indicate the average and shaded areas indicate the standard deviation over 10 independent stochastic runs.
https://doi.org/10.1371/journal.pcbi.1010026.g010 large binding rate (0.8/s) results in a spread out actin organization, with less mobile motors given their long residence time on the filaments. Fig 11A and 11B shows that the average actin contractility measure undergoes a very slow decrease following the switch to motor binding rate 0.4/s and that the myosin motors slowly become more localized in the simulation domain (see S4 Video for a sample simulation). Similarly, this parameter change leads to more localized actin organization, as illustrated by the better contouring of one peak (corresponding to a more compact cluster) in the radial distribution function in panel C1 of Fig 11. This is reminiscent of in vitro models of the actomyosin cortex where the network contracts into foci that then further coalesce with each other [59]. Switching from binding rate 0.4/s to 0.8/s shows a slow relaxation from the network's contractile behavior (Fig 11A) but no considerable change  Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 independent stochastic runs.
https://doi.org/10.1371/journal.pcbi.1010026.g011 in the myosin organization (Fig 11B) or in the pairwise distances between actin cylindrical segments (Fig 11C). Regulation of the binding rate earlier in the process has a similar impact on actin contractility and myosin motor localization (see S1 Fig). The switch from 0.8/s to 0.4/s binding rate, which in this case occurs in a more transient state of the network (2000 s into a 7000 s simulation), shows a more distinct progression from two peaks to one peak in the radial distribution function, given the occasional transient two-cluster actin organization driven by the 0.4/s myosin binding rate. These observations suggest that a significant change from a contractile actin-myosin network organization may require control from additional regulatory processes.

Discussion
Cortical actin undergoes dynamic reorganization in many organisms, adopting a wide variety of configurations from homogeneous meshes to spatially localized clusters. Experimental observations show that many members of the myosin motor family (such as myosin II, V, and VI) may be involved in cytoskeleton remodeling. Even within one myosin family, small differences in motor activity between isoforms can result in different organization of the actin cortex [6]. Motivated by these observations and using an established stochastic agent-based modeling framework, here we investigate: a) the impact of varying kinetic parameters for motor protein activity on the organization of the actin cytoskeleton, b) the segregation and complex structure generation driven by two motor populations interacting with the same actin meshwork, and c) whether regulation of motor behavior at a certain time in development can lead to remodeling of the actin cytoskeleton. We propose data analysis measures that assess both the dynamic remodeling of the actin network as well as the clustering and segregation of the myosin motor proteins. While motivated by understanding the interactions of certain myosin motors that have been found to collectively organize the cytoskeleton, we do not constrain our model to specific motor proteins and their corresponding properties, but rather aim to understand how regulation of individual motor kinetic rates can lead to diverse actin network structures.
Overall, we find that cytoskeleton organization is highly sensitive to the kinetics of interacting motor proteins. Here we focus on the role of several key kinetic parameters (Table 1): binding rate, stall force, motor step size, and the number of heads per minifilament. We also studied the influence of motor parameters such as the stretching force constant, the per-head unbinding force, and the reaction range of the motor binding reaction on actin organization; we find that the ranges considered for these parameters (see Table 1) do not yield significantly different actomyosin organization from the baseline.
When acting individually, a single type of motor protein can produce a range of actin organizations when interacting with an ensemble of actin filaments. Tight clusters of actin filaments with a smaller radius of gyration are associated with large motor step sizes, higher stall force, and higher number of heads. Homogeneous networks, where the actin filaments are loosely spread on the domain, are associated with lower values of those parameters. Variations in the per-head motor binding rate suggest an optimal intermediate value is needed to organize the actin filaments into tight clusters, with low and high values of this parameter resulting in a loose network. We observe rich network remodeling by the myosin activity, including the emergence of bundles, asters, and homogeneous distributions that have been previously observed in computational and experimental studies. Features of the simulated cytoskeleton organization that are shared with previous work include filament polarity sorting in bundles and asters [27,41,43,44], polar asters with barbed ends oriented towards the center [43,52,53,55,56], and localization of myosin motors in the center of the asters [27,40,[52][53][54][55][56]. These similar results confirm the accuracy of our model simulations, which further provide a comprehensive study of how individual motor parameters regulate these distinct patterning regimes.
In addition, we also take of advantage of the advanced features of the MEDYAN stochastic modeling framework to give insight into cytoskeleton network dynamics driven by two motor populations. We find that when two distinct populations of motor proteins interact with the actin filaments, unexpected behaviours can emerge. Many motor combinations appear to compromise, for instance when the motors have different step sizes or motor heads, with network measures taking on values that lie between the values when the motors act alone. This suggests that some motor protein properties are complementary and can work together to produce novel organizations of actin filaments. In contrast, some combinations of motor protein populations appear to act antagonistically, with one motor protein dominating the organization of the meshwork. For instance, motor proteins with a large step size dominate motor proteins with a small step size, with network measures in the presence of both motor proteins almost indistinguishable from those when the motor with a large step size acts alone. In this case, the less dominant motor acts as a passive cargo, being transported to particular regions in the domain by the activity of the dominant motor; this also means that the passive motor is only able to find actin filaments to bind to based on the cytoskeleton organization imposed by the dominant motor. This relationship between the motors can be advantageous if the less dominant motor requires a particular filament configuration or spatial localization but cannot achieve the required dynamics when acting alone. Further, in the presence of two motor protein populations, spatial segregation can be achieved, where motors that are initially homogeneously distributed on the domain self organize into distinct regions. This type of motor protein sorting has been previously observed [48], and appears to rely on differences in step size and binding rates in our study. Spatial exclusion of filament-bound cytoskeletal components has been an active area of research. Computational studies on spatial segregation of crosslinker populations have shown that effective separation is dependent on actin polymerization rate and mechanical properties such as the size difference between two crosslinkers and the bending modulus of actin filaments [60]. Here, we show that kinetic parameter differences are sufficient to achieve spatial segregation in actin networks with heterogeneous motor populations.
Dynamic transitions in actin network organization can be realized when motor protein kinetics are changed as a result of differential gene expression or regulatory pathway dynamics. These kinetic rate changes could be due to inactivation or degradation of one type of motor protein with simultaneous activation or synthesis of another motor protein with different kinetics, or regulatory proteins hiding or exposing functional sites on the motor protein. We find that, when activated in the presence of a particular organization of actin filaments, some motor proteins are unable to remodel the existing actin meshwork, while other motor proteins are able to reconfigure the actin filaments into new structures. This could be due to the fact that, on top of potential regulation by motor proteins, dynamic actin remodeling is likely to also be influenced by other regulating factors. For instance, recent in vivo studies have found that the organization of the cytoskeleton into actin patches, rings, and asters may be primarily driven by the Arp2/3 complex in living HeLa cells and in immune cells [26,57]. However, myosin II is still found to play a significant role in the progression and disassembly of actin patterns [57].
By rigorously quantifying actin filament organization in response to motor protein kinetics in detailed stochastic simulations, we have demonstrated a range of possible actin-based structures. Acting individually, motor proteins can produce many of these structures, but efficient transitions between structures require motor proteins to act together. Changes in motor protein kinetics and cooperation between motor protein types may account for the large scale changes in actin filament organization observed in vivo.

Network contractility
To assess the contractile behavior of each actomyosin network generated, we calculate the network radius of gyration [17]. This measure has been shown to be effective in determining the filament contractility in MEDYAN simulations [17]. Each filament in MEDYAN is stored in terms of the locations of the coarse-grained cylinder segments (monomer units) that make up each actin filament. We let n be the total number of cylindrical segments from all actin filaments in one time frame of a MEDYAN simulaiton. Let r i = (x i , y i , z i ) be the location of the ith cylinder (with coordinates in 3-dimensional space) and determine the geometric center of the ensemble of all cylinders as r GC = (mean(x i ), mean(y i ), mean(z i )). Then the network radius of gyration for that time frame is defined as [17]: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi 1 n When evaluated at a time-series of MEDYAN simulation frames, the network radius of gyration has a decreasing pattern when the contractile behavior increases through time, and increases when the contractility of the network decreases through time. In visualizations of this contractility measure, we normalize the network radius of gyration R G by dividing by its value at the first time point in the simulation.

Global network alignment
To determine the alignment of the actin filaments in the actomyosin network, we calculate an orientational order parameter of the system of actin filaments [17]. This involves setting up the ordering tensor: and defining the orientational order parameter S as the largest eigenvalue of this tensor; we note that this measure has also been used to assess the alignment of molecules in liquid crystal models [61]. Here, N is the number of filaments in the actomyosin network and I 3 is the 3 by 3 identity matrix. The normalized direction vectorr i of each actin filament i is calculated based on the 3-dimensional locations of the two filament ends (the locations of the first and last segments rather than the locations of each cylindrical segment between the ends), which allows for this value to reflect the alignment of the network even when the filaments are bending [17]. If the locations of the ending cylinders are given by (x 1 , y 1 , z 1 ) and (x 2 , y 2 , z 2 ), then with l the length of r i . The value S = 0 of the largest eigenvalue of Q corresponds to random alignment of the filaments in the system, while the value S = 1 indicates a perfectly aligned filament network [17,61].

Radial distribution function
We calculate a variation of the radial distribution function to understand the distances between emerging structures in the simulated polymer network. This involves computing the distances between all pairs of actin cylindrical segments in the simulation and binning the distances into a distribution. Letting r i = (x i , y i , z i ) be the location of the ith cylinder as in § 5.1, we determine the matrix of pairwise distances Z, where Z ij = d(r i , r j ) is the Euclidean distance (L 2 norm) between actin monomer unit i and j. Noting that the distance between segments ranges from 0 nm to 2000 ffi ffi ffi 2 p nm (the maximum distance if the actin segments are at opposite corners of the domain), we divide this range into 50 bins and denote the centers of the bins by radius j . For each time frame t, we then define: Here 1 A ðxÞ is the indicator function with value 1 when x 2 A and 0 when x = 2 A. N s is the number of cylindrical segments at time t and therefore the normalization is done by dividing by the number of pairs of actin cylindrical segments.

Motor localization
To quantify the spatio-temporal localization of motors in the simulation domain, we divide the simulation domain into cylindrical annuli and determine the number of motors bound to filaments in each such volume and at each time. The thin z dimension of the simulation domain gives the height of each cylinder and the circles centered at x = y = 1000 nm with radii 0, 250, 500, 750, and 1000 nm provide bounds between the annuli. Note that the first volume is actually a cylinder with center at (1000, 1000), while the last volume extends outside the boundaries of the cubic simulation domain (this is no concern since motor proteins will simply not be found there). Using the locations of the centers of the minifilaments m x , m y , m z , we record the number of myosin motors that are bound to filaments at each time point and count how many are located in each cylindrical annulus.

Motor spread
We aim to quantify the spread of a myosin motor population in MEDYAN actin-myosin simulations. To simplify computation and due to the small height of the domain, we restrict our attention to the centers of the minifilaments in x-y space (m x , m y ). We apply the boundary and polyshape functions in MATLAB to the population of motors, thus generating a 2-dimensional boundary polygon around the motors. We use a default shrink factor of 0.5 for the boundary method, which means that the resulting polygon around the myosin motors is tighter than the convex hull of the points. Let the polygon around the motor population at time t be denoted by PðtÞ; then we introduce a measure for the myosin motor spread in the domain: where L = 2000 nm is the side of the square domain and thus the polygon area is normalized by the area of the 2-dimensional simulation domain considered.

Motor segregation measures
Building on the framework for the motor spread measure in § 5.5, we introduce two measures for determining the segregation of two motor populations in MEDYAN actin-myosin simulations. We similarly consider the centers of the motor minifilaments in two dimensions (m x , m y ). Using the boundary and polyshape functions in MATLAB for each population of motors as above, we generate two-dimensional boundary polygons around each motor population. We denote the polygon around the first motor population at time t by P 1 ðtÞ and the one around the second motor population by P 2 ðtÞ. Let (c x,1 , c y,1 ) and (c x,2 , c y,2 ) be the twodimensional positions of the centroids of polygons P 1 ðtÞ and P 2 ðtÞ. We then define the following measures for the normalized intersection area and the distance between the centroids of the two polygons: where L = 2000 nm is the side of the square domain. The measure A int is normalized against the area of the union of polygons P 1 ðtÞ and P 2 ðtÞ, so as to capture the intersection area relative to the space that both motor populations cover. The measure D cent is similar to the separation distance measure proposed in [62] for the distance between F-actin and myosin-II fluorescence areas.

Spatial statistics
Spatial statistical methods are useful in understanding the distribution patterns of proteins [63]. Motivated by the use of protein pattern analysis in microscopy images as described in [63], we use the K-Ripley function to understand how random, cluster, or regular distributions may form in the simulation domain for actin monomer units and myosin motor proteins. As in the previous method, we focus on the locations of proteins in the x-y space. For actin, we sample 30% of the monomer units along each filament (as done in [49]) to obtain the corresponding point process. For motor proteins, we directly use the locations of the centers of the myosin minifilaments. To analyze these point processes, we calculate the K-Ripley function (using the spatstat function in R), which measures the number of neighbors within a certain radius r to a point [63]: where λ is the density of points in the studied region, d ij is the distance between points i and j, and N is the number of points in the dataset. We record the normalized form H(r) of the K-Ripley function: HðrÞ ¼ ffi ffi ffi ffi ffi ffi ffi ffi ffi KðrÞ p r À r ð8Þ and note that H(r) = 0 for complete spatial randomness, H(r) > 0 for clustering, and H(r) < 0 for regularity in the distribution of the point process. We therefore record the signed area under the curve of H(r) at nine time points throughout the simulation (every 250 seconds); larger values for this measure correspond to more clustered patterns in the distribution of actin monomers or of myosin motor proteins.  Table 1). Solid lines indicate the average and shaded areas indicate the standard deviation over 10 independent stochastic runs.

Supporting information
(EPS) S1 Video. Evolution of the cytoskeleton network for baseline parameters. Sample actinmyosin cytoskeleton organization using MEDYAN for the baseline parameters in Table 1.