Dynamic Mechanisms of Cell Rigidity Sensing: Insights from a Computational Model of Actomyosin Networks

Cells modulate themselves in response to the surrounding environment like substrate elasticity, exhibiting structural reorganization driven by the contractility of cytoskeleton. The cytoskeleton is the scaffolding structure of eukaryotic cells, playing a central role in many mechanical and biological functions. It is composed of a network of actins, actin cross-linking proteins (ACPs), and molecular motors. The motors generate contractile forces by sliding couples of actin filaments in a polar fashion, and the contractile response of the cytoskeleton network is known to be modulated also by external stimuli, such as substrate stiffness. This implies an important role of actomyosin contractility in the cell mechano-sensing. However, how cells sense matrix stiffness via the contractility remains an open question. Here, we present a 3-D Brownian dynamics computational model of a cross-linked actin network including the dynamics of molecular motors and ACPs. The mechano-sensing properties of this active network are investigated by evaluating contraction and stress in response to different substrate stiffness. Results demonstrate two mechanisms that act to limit internal stress: (i) In stiff substrates, motors walk until they exert their maximum force, leading to a plateau stress that is independent of substrate stiffness, whereas (ii) in soft substrates, motors walk until they become blocked by other motors or ACPs, leading to submaximal stress levels. Therefore, this study provides new insights into the role of molecular motors in the contraction and rigidity sensing of cells.


Introduction
Cells modulate their properties and activities in response to the surrounding environment, via morphological rearrangements driven by cytoskeletal contractility and reorganization. Various quantitative studies using gels with tuned elasticity have provided insights into the understanding of how cells respond to matrix stiffness [1,2,3,4]. On soft substrates, cells generate low forces with randomly aligned actin filaments, leading to a weak response with wrinkles or strains of the substrates. By contrast, stiff substrates result in extensive cell spreading and enhance contractility with numerous stress fibers. Other experimental results collectively suggested that on stiff substrates, cells tend to deform intracellular structures rather than the substrate as seen in myosin/actin striations [5,6,7,8].
Several mechanisms governing such mechano-sensing of cells have been proposed in experimental studies, and multiple mechanisms likely exist involving different intracellular structures. For example, a large number of mechano-sensing molecular motifs that vary conformation over a range of mechanical forces transduce mechanical signals into biochemical ones [9,10,11]. It has also been believed that actomyosin contractility contributes to cell mechano-sensing [12,13,14]. For example, non-muscle myosins were shown to be crucial for stem cells to sense matrix elasticity [15]. Local forces acting on both integrin-mediated [16,17] and cadherin-mediated adhesions exhibit a similar relationship with stiffness [18,19].
Different phenomenological laws have been proposed to explain the substrate-dependent mechano-sensing. For example, a simple ''two-spring model'' predicted that stiffer environments lead to stronger traction forces [20]. A ''three-spring model'' was proposed later to explain the stiffness-dependent orientation of stress fibers in adherent cells [21]. To elucidate interactions of molecular motors with adhesion complexes in the mechanosensing process, a different theoretical model based on active matter theory was proposed [22]. It demonstrated that for short timescales (t%100 s), mechano-chemical transduction from the motors plays a dominant role since the adhesion complexes are unlikely to have enough time to recruit associated proteins. Concurrently, numerous computational models have been developed to elucidate the mechanisms of mechano-sensing. For instance, they showed that actin networks can adjust to mechanical environments by modulating cross-links within the networks [23], and also suggested a mechanism for stiffness-sensing of cells adhered to a compliant surface mediated by actin filament alignment in the direction of force application [24].
Taken together, these recent experimental, theoretical, and computational efforts have led to new insights about the structural reorganization of the cytoskeleton as well as the effects of extracellular stiffness on cell behaviors. However, little is known about the roles of actomyosin contractility in mechano-sensing on timescales of hundreds of seconds, which are biologically relevant. In this work, using a Brownian dynamics computational model [25,26], we investigate the large-scale contractile responses of an actomyosin network on timescales of hundreds of seconds, during which protein recruitment and responses from molecular motifs can occur, to elucidate one actomyosin-driven rigidity-sensing mechanism that functions under diverse conditions. Specifically, the effects of external elasticity on cytoskeletal contractility and network morphology are evaluated by systematically varying model parameters, e.g. the concentration and kinetics of motors. Our simulations successfully reproduce some of the large-scale mechano-sensing responses of cells such as active contractility and force generation, in good agreement with recent experimental observations [12,27,28]. Merely by modeling actin and myosin activity in the absence of proteins related to adhesion complexes, we predict both equilibrium and dynamic behaviors, indicating that actomyosin machinery can function as a stand-alone mechanism for the mechano-sensing of cells.

Model
We use a previous agent-based model [26] based on Brownian Dynamics to simulate active cross-linked actin networks as systems that generate force as well as sense surrounding mechanical conditions. In this approach, we explicitly take into account actin filaments, ACPs, and molecular motors and their local interactions. To facilitate understanding of the results predicted in this study, we briefly present their main features.

Formation of an Active Actin Network
Active actin networks with motors are generated in a similar fashion to previous studies [26]. Monomers of actin (G-actins), passive ACPs, and motors are assembled into a network via reversible reactions in a 3-D cubical domain with periodic boundary conditions in all directions. Actin can exist in either monomeric or filamentous form while ACPs and motors can exist in three states: monomeric (free), inactive (partially bound), and active (bound to two filaments) states. Note that following the initial formation of the network, monomeric ACPs and motors are implicitly considered via their local concentration and secondorder reaction equations. After concentrations of G-actin, ACPs, and motors reach a dynamic steady state, residual G-actins are deleted with actin assembly/disassembly deactivated for simplicity. A geometrically identical network is used in all simulations to isolate the effects of the stiffness of the surrounding medium and other parameters ( Figure 1A). To vary the concentration of motors, they are removed from networks or added as monomers at the beginning. The average filament length (SL f T) is ,2 mm, actin concentration, C A , is 12 mM, density of ACPs, R ACP ( = C ACP /C A ), is 0.01, and the initial width of the cubical domain is 5.0 mm. Density of motors, R M ( = C M /C A ), is 0.02 unless specified.

Mechanics of Actin Filaments, ACPs, and Motors
Actin filaments comprise cylindrical segments of length 140 nm (r 0,A ), and both the ACPs and motors are represented by two arms parallel to each other spanning between cross-linked actin filaments a distance of 70 nm (26r 0,ACP ) and 140 nm (26r 0,M ), respectively. Motions of the network components are governed by the Langevin equation: where m is the mass of each element (actin, ACP, or motor), r is the element's location, f is the friction coefficient, t is time, F B is a thermal force satisfying the fluctuation-dissipation theorem, and F is a net deterministic force including extension, bending, and repulsive forces. Since inertia of all elements is negligible on the length and time scales of interest, positions of the elements are updated using the Euler integration scheme: where Dt is a time step. Extension and bending of the cylindrical segments constituting actin filaments, ACPs, and motors are computed using simple quadratic potentials, denoted by subscripts ''s'' and ''b'' respectively: where r is bond length, k s is extensional stiffness, h is bending angle, k b is bending stiffness, and the subscript 0 denotes an equilibrium (zero-force) value. As in our previous studies [26], bending stiffnesses are introduced to restrict actin filament bending (k b,A ), keep the two arms of ACP (k b,ACP,1 ) or motor (k b,M,1 ) parallel, and maintain the right angle between the axis of a filament and the arm of ACP (k b,ACP,2 ) or motor (k b,M,2 ). Specific values of the geometrical and mechanical parameters are listed in Table S1. In addition, the repulsive force is responsible for volume-exclusion effects by which actin filaments cannot pass through each other, which is calculated by the following harmonic potential, U r , depending on the minimum distance, r 12 , between two cylindrical segments [26]: where k r is the strength of repulsive effects, and r c is the diameter of cylindrical segments. Then, the repulsive force is distributed to the two ends of the actin segment based on the relative location on the segment where r 12 is measured.

Dynamic Behaviors of ACPs and Motors
We assume that each motor in our simulation corresponds to a single myosin minifilament consisting of multiple myosin II molecules. As described in our previous studies [29], motors in the active state walk along actin filaments toward a barbed end at a rate, k w , depending on the extensional force acting on the arm, F F s~+ U s : where d w 's and l w 's are time constants and mechanical sensitivities for walking of motors (Table S1), respectively, andt t is a unit vector locally tangent to an actin segment in the direction of a pointed end. Although motors in this study mimic a myosin minifilament consisting of numerous myosin II molecules, Eq. 6 and the values of d w 's and l w 's are adopted from a single-molecule experiment examining myosin V under 1 mM ATP. Our intention was to model generalized motor activity, and we chose myosin V because it has been extensively characterized. Nevertheless, the loaddependent walking rate of the minifilament is still qualitatively similar to that of myosin V, justifying the use of Eq. 6 for roughly mimicking myosin minifilament behavior. As seen in Eq. 6, only tension (r §r 0,M ) directed to a pointed end (F F s : t t §0) affects k w , resulting in a stall force, ,4 pN, beyond which motors cease walking.
In addition, as in [30], ACPs and motors are able to unbind in a force-dependent manner following Bell's equation: where k 0 u is the zero-force unbinding rate coefficient for ACPs (k 0 u,ACP ) or motors (k 0 u,M ), and l u is the mechanical sensitivity for unbinding of ACPs (l u,ACP ) or motors (l u,M ) ( Table S1). Note that although unbinding of motors is also one of the phases of walking, these two events are considered separable for systematic analysis. If the arm of motors reaches the barbed end of a filament by walking, it remains there until it unbinds.

Boundary Conditions of the 3-D Computational Domain
After obtaining the network, actin filaments crossing the domain boundaries are severed and permanently clamped with periodic boundary conditions deactivated in all directions. During the measurement of strain and stress, the boundaries also act as sticky surfaces to take the binding between actin filaments and membrane into account; if either end of an actin filament is located within 30 nm of a boundary, the end is irreversibly clamped.
Normal stress (s) on each boundary is the sum of normal forces exerted by actin filaments clamped on the boundary, divided by area. s is used to compute movement of the boundaries in simulation; we assume that the domain is surrounded by an elastic medium with identical Young's modulus, E, on all boundaries. Each boundary (assumed planar) experiencing s is displaced a distance corresponding to a strain s/E.
Simulations begin with zero stress on all boundaries and proceed over time for 200 s. At this point, the network reaches a steady state stress in most cases and we define this to be the plateau stress. However, in a few cases, stress continues to slowly rise even after 200 s.

Measurement of network stiffness
In order to measure the stiffness of networks, we applied differential sinusoidal normal displacement of amplitude 280 nm to the networks and calculated the responding stress. For the purpose of this calculation, all the motor and actin cross-linking dynamics was deactivated to probe the instantaneous network stiffness, avoiding any progressive time-dependent changes in the network. Under these conditions, the networks exhibit a predominantly elastic response as indicated by the small phase delay between the applied strain and the responding stress ( Figure S1). We then calculated network stiffness by dividing the amplitude of stress by that of strain.

Results
Here, we investigate the role of molecular motors as rigidity sensors and predict the contractile (normal) stress and strain of actomyosin networks tethered to 3-D cubical domains. We examined these as a function of the various kinetic parameters and concentrations of motors as well as different elasticity of the surrounding medium.

Network morphology and stress evolution depend on substrate stiffness
The initial network ( Figure 1A) starts from a zero stress condition. Due to motor activity, and depending on substrate stiffness (E), the network shrinks to different extents at different rates. Lower E leads to shrunk and concentrated networks with highly bent actin filaments ( Figure 1B). On the other hand, higher E prevents the domain contraction, forming heterogeneous networks with highly stretched filaments ( Figure 1B). These differences in network morphology have been reported in experiments where they found, for different cell types, that Factin networks tend to be denser and less organized on more compliant substrates [31,32]. Stress (s) in all cases rapidly increases at the beginning although the rate of increase gradually falls, rising at a much slower rate by ,200 s in most cases ( Figure 2A). Recognizing that stress continues to rise after this time, but constrained by computational resources from extending the calculations further, we use the value of stress at 200 s as a reference, and denote it as the ''plateau stress'', s p . For E,3 kPa, s p is proportional to E but becomes relatively constant for E$3 kPa, which corresponds well to literature [3,12,16,27] ( Figure 2B). The maximum of s p is ,420 Pa. The initial slope of stress, _ s s 0 , measured at t,10 s increases swiftly for E,3 kPa and slower for E.3 kPa ( Figure 2D). The initial strain rate, _ e e 0 ( = _ s s 0 / E), decreases with greater E ( Figure 2C); since contraction is associated with energy expenditure to overcome the internal friction and the rupture of cross-links, cells contracting against softer substrates will experience larger energy dissipation, leading to the slower rise in stress. The ''plateau strain'', e p (strain at plateau stress), decreases with greater E falling below 0.05 for E.10 kPa ( Figure 2E). s p and e p with various E show a first zone where s p rapidly changes, followed by a period of slower increase, which agrees well with [33]. We also measured mechanical power, P = s 0 _ e e 0 V where s 0 is initial stress corresponding to _ e e 0 , and V is the instantaneous volume of the domain. P exhibits a bimodal dependence on s 0 , having a peak at E,0.6 kPa ( Figure 2F). At this peak, the network exerts 40% of the maximum s 0 with intermediate _ e e 0 (,0.015 s 21 ), compared to cases with other E.

Network stiffness tracks the generated stress
It has been recently found that cell stiffness tracks substrate stiffness over a range of stiffnesses before reaching a constant value [34]. We measured the steady-state stiffness of networks at each E. Network stiffness (E n ) was found to be proportional to (and nearly equal to) s p over the entire range of E, but proportional to E only up to a value of E,3 kPa ( Figure 2B). This tendency is consistent with the direct proportionality between prestress and G9 (or K9) of passive actin networks observed in experiments [35].

Effects of motor concentration
Motor density is varied by adjusting the initial concentration of motors in the network, R M . In all cases, s increases with E and then exhibits a much slower rate of increase at high E ( Figure 3A), but compared to the control case (R M = 0.02), the tendency is less clear in the other cases, especially for low R M where the dependence between s p and E weakens. For low R M (,0.02), s p tends to be higher with greater R M , in good agreement with literature [12,36,37]. With a maximum at R M = 0.02 for most E, s p drops for higher R M . High contractile activity of the network enhances the rate of stress generation. Thus, _ s s 0 and _ e e 0 increase for all values of E until R M reaches the optimal level explained above ( Figure 3B and C).

Effects of unbinding and walking behaviors of motors
Motor unbinding is explored by varying the zero-force unbinding rate (k 0 u,M ) and the processivity (l u,M ). On the other hand, motor walking is studied by varying the sensitivity (l w ) which is equivalent to variation of the stall force.
Although higher k 0 u,M results in more frequent unbinding, it does not necessarily lead to lower s p since unbinding can also help stalled motors due to blocking effects to bind to other binding sites so that they can keep walking. However, too frequent unbinding prevents motors from remaining attached to filaments for enough time to generate large stress. The early phase of s evolution is practically unaffected by changes in k 0 u,M while the later phase is  Figure 4D-F). Note that l w refers to both l w,1 and l w,2 , and their values are varied simultaneously. In response to variations of l u,M , the typical tendency of s(E) is conserved in most cases except that with l u,M = 76l Ã u,M where the dependence on E is noticeable only at low E since motors can bear very small forces ( Figure 4A). Therefore, _ s s 0 and _ e e 0 deviate from the control case only for high values of l u,M ( Figure 4B and C). Interestingly, P(s 0 ) and _ e e 0 (s 0 ) after normalization are collapsed into a unique curve ( Figure S3A and C). Regarding motor walking, l w determines the stall force at which the walking rate of motor approaches nearly zero (Eq. 6). Lower values of l w lead to higher stall forces, so the motors can overcome the applied forces and thus walk longer distances. Multiplying l w by 0.1 increases s p nearly three times compared to the control case ( Figure 4D). Besides, due to higher _ s s 0 , the domain shrinks more notably ( Figure 4E and F). Again, P(s 0 ) and _ e e 0 (s 0 ) collapse well into a single curve after normalization ( Figure S3B and D).

Discussion
Cells are capable of adapting their properties through a variety of mechanisms. In this study, we investigated the role of molecular motors as rigidity sensors and also propose how these motors can induce such precise mechano-sensing. Despite the oversimplifications and assumptions made in the described system which only includes the dynamics of motors and ACPs, the simulated actin network exhibits macroscopic contractile behaviors remarkably similar to several experiments [12,27,28], further indicating that microscopic properties of individual constituents govern the network responses, and that motors play a central role in mechano-sensing. However, this does not negate the significance of other factors such as actin dynamics and structures, biochemical signaling, and adhesions dynamics in the cell's response and adaptation to mechanical cues. Rather, we demonstrate here that actomyosin machinery can be one of several possible mechanisms for cell rigidity-sensing phenomena. Nevertheless, we explored a wide range of parametric spaces in order to study how different parameters of the model influence cell adaptation, finding that those parameters affecting the kinetics of motors are the most critical for cellular adaptation to substrate stiffness.
We probed diverse contractile large-scale characteristics by evaluating network morphology, plateau stress (s p ), the initial increasing rate of stress ( _ s s 0 ), the initial rate of strain (_ e e 0 ), and mechanical power (P) over a wide range of parameters -the stiffness of the surrounding environment (E), the concentration (R M ) and the dynamics of motors (k 0 u,M , l u,M , and l w ). It was observed that softer substrates lead to shrunken and dense networks, whereas stiffer ones result in heterogeneous networks with minimal domain contraction ( Figure 1B). Overall, the qualitative pattern of stress evolution is quite consistent with multiple recent experimental works [12,14,27]; s increases with time rapidly at first but reaches s p in most cases before 200 s (Figure 2A). In addition, we found that s p is roughly proportional to E for E,3 kPa and becomes relatively constant for E.3 kPa ( Figure 2B). The existence of a transition to a slower rate of stress increase can be explained by the mechanisms that cause the motors to slow or stall: (i) all of the next binding sites in a barbedend direction are already occupied (blocking), (ii) reaching the motor stall force, or (iii) reaching the barbed end of an actin filament. Figure 5A demonstrates that only a small fraction of motors reach the barbed end of a filament for all E, so this would have little direct influence on s p . Blocking, on the other hand, is observed over the entire range of E, but is especially prevalent at lower E. This is due to the greater distance that motors need to walk before reaching their maximum force, combined with the tendency for all constituents (filaments, motors, and ACPs) to increase in density under large negative strains. For stiffer substrates, material strains are smaller and motors walk shorter distances before attaining the stall force. For E.3 kPa, s can reach s p determined by the stall force that motors can exert maximally while at lower E, s p is limited by the blocking effect which progressively decreases as E increases ( Figure 5B). This transition from blocking at low E to limitation due to motor stall force at high E constitutes a mechanism by which cells can sense substrate stiffness. Forces transmitted along the cytoskeleton and across adhesion complexes will vary according to the generated stress, leading to varying degrees of conformational change in these stress-bearing proteins. Since conformation determines biochemical activity, factors such as exposing cryptic binding sites, changes in binding affinity, or phosphorylation, for example, will modulate signaling activity and therefore, cell function. This mechanism is further confirmed by effects of C A on s p ( Figure S4). At E = 2560 Pa, we varied C A but maintained C M at 0.24 mM, meaning that R M decreases with higher C A . Note that C M = 0.24 mM is equivalent to R M = 0.02 in the control case with C A = 12 mM. s p is proportional to C A even with the same number of motors ( Figure S4A). More actins can provide the motors with greater space to walk, leading to fewer stalling events by blocking but more frequent stalling by applied forces ( Figure S4B). In addition, the effects of average filament length (SL f T) demonstrate the mechanism. We found that networks attain stress roughly proportional to the cube of SL f T ( Figure S5A). With shorter filaments, more motors reach the barbed ends while with longer filaments, motors are more likely to stall by attaining their maximum level of force generation ( Figure S5B). Interestingly, the percentage of blocking events remains relatively constant, independent of SL f T. The increase in the number of motors exerting maximum stall forces results in higher s p , consistent with our mechanism above.
Several recent studies suggested that the proportionality between s p and E could be correlated with _ s s 0 [12,27,28], which is consistent with our observation that _ s s 0 increases swiftly for E,3 kPa but slows down for E.3 kPa ( Figure 2D). From a biological point of view, cells are polarized and migrate in the direction of higher stiffness inducing the faster increase of traction force [3,38,39,40,41]. The rise of _ s s 0 with increasing E is in the same order of magnitude of the experimental findings [12,27,28]. Although the time required to reach s p is somewhat different between ,200 s in our simulations and ,600 s in experimental studies, it is common that the time needed to reach s p is relatively independent of E in both, as found in [12] (Figure 2A). The inverse proportionality between e p and E ( Figure 2E), consistent with experiments [33], supports that cells on stiff substrates tend to rearrange intracellular structures rather than deforming the substrate, as in [5,6,7,8]. Numerical results of _ e e 0 depending on E ( Figure 2C) and of P(s 0 ) ( Figure 2F), show that the emergent behavior of the network follows Hill's equation for muscle contraction [42]. We fit the _ e e 0 (s 0 ) curve ( Figure S6) with the relation s 0 za ð Þ _ e e 0 zb ð Þ~c where a = 10.0 Pa, b = 0.0103 s 21 , and c = 1.0023 Pa/s. In experiments using various types of muscle cells [43] and myoblasts cells [12], introducing a shape factor r~a=F max~b =V max &0:25 normalized the data. For our control case, we obtained r 1~a =s 0, max &0:12 and r 2~b =_ e e 0,max &0:4. We found that the values of these factors are regulated by parameters, such as k 0 u,M , l u,M , and l w , although normalizing the data leads to similar curves ( Figure S3A-D).
We provided further insights regarding the effects of motor concentration (R M ). It was observed s p attains a maximum at R M = 0.02 ( Figure 3A). This could be attributed to the limited binding sites for motors on actin filaments in our simulations. However, note that R M in our model corresponds to the concentration of multimerized myosin II structures in cells rather than that of individual molecules. Therefore, R M = 0.02 is actually a very high density of large motor structures, and thus cells are likely to have such an optimal R due to blocking effects between the large aggregates of myosins. A refined model including multiple myosin heads per motor and multiple binding sites per actin segment would help to clarify this issue.
Considering the variable extent of multimerization of myosin II molecules into a minifilament or thick filament, we evaluated the effects of the zero-force rate of motor unbinding (k 0 u,M ) ( Figure S2) and the mechanical sensitivity for unbinding (l u,M ) and walking (l w ) on network contraction. Motors with high l u,M more readily unbind, and those with high l w are more likely to be stalled at small forces.   (Table S1), and numbers in the legends indicate n. Note that A and B share a legend with C, and D and E share a legend with F. s p tends to be higher with lower l u,M and l w which correspond to more processive and stronger motors, respectively. doi:10.1371/journal.pone.0049174.g004 curves is conserved well in both cases. By contrast, with high k 0 u,M , s p is substantially reduced while _ s s 0 and _ e e 0 are not affected ( Figure  S2A-C).

Conclusion
We elucidated one mechanism by which cells can modulate their properties and respond to the surrounding environment via cytoskeleton contractility, using an agent-based computational model. Although the model is based on molecular-level processes, macroscopic behaviors of the active cross-linked actin networks agree well with the response of cells probed in experimental quantitative studies [1,2,3,4]. We found that the biphasic relation between substrate stiffness and the level of generated forces [5,6,7,8] is attributable to a transition from stalling due to steric hindrance or ''blocking'' in soft substrates to that due to stall forces in stiff substrates. In addition, we showed that in response to increases in substrate stiffness, the contraction rate of cells increases while the corresponding contraction velocity decreases, also consistent with experiments [28]. All of these suggest that actomyosin contractility is one plausible stand-alone mechanism capable of contributing directly to cell mechano-sensing [12], consistent with various experimental findings that myosins are crucial for cells to sense surrounding matrix elasticity [14,15], and that cell responses to rigidity of the external matrix reflect adaptation of the actomyosin machinery to load following Hill's relation [12]. Fraction of motors stalled due to: (i) high applied forces (black), (ii) blocking (gray), or (iii) arrival at barbed ends of filaments (white) at steady state as a function of C A . At low C A , ,50% of motors are not stalled since many of them lie in the inactive state due to lack of network percolation. As C A increases, motors are more likely to be stalled due to high forces as opposed to blocking. (TIF) Figure S5 Effects of average actin filament length (SL f T). (A) s p increases dramatically as SL f T is increased. (B) Fraction of motors stalled due to: (i) high applied forces (black), (ii) blocking (gray), or (iii) arrival at barbed ends of filaments (white) at steady state as a function of SL f T. As SL f T increases, more motors are stalled due to attaining their maximum force while fewer motors are stalled due to arrival at barbed ends. The number of motors stalled due to blocking remains nearly constant regardless of SL f T. Table S1 List of model parameters. Numbers in parentheses are corresponding dimensionless values as defined in the text. ''*'' on symbols indicates reference values of parameters studied in the sensitivity analysis. Values marked by ''1'' are adopted from given literature with adjustment based on assumption that motors in this study consist of many myosin II molecules. Note that to increase length (N c ) and time scales (Dt), k s,A used in this study is 4 times smaller than that in our previous works, but it was confirmed that the results are virtually unaffected by the 4-fold decrease. (DOC)