Pattern formation and the mechanics of a motor-driven filamentous system confined by rigid membranes

Pattern formation and the mechanics of a mixture of actin filaments and myosin motors that is confined by a rigid membrane is investigated. By using a coarse-grained molecular dynamics model, we demonstrate that the competition between the depletion force and the active force of the motors gives rise to actin accumulation in the membrane vicinity. The resulting actomyosin structure exerts pressure on the membrane, that, due to nematic alignment of the filaments, converges to a constant for large motor active force. The results are independent of filament length and membrane curvature, indicating the universality of this phenomenon. Thus, this study proposes a novel mechanism by which the compounds of the cytoskeleton can self-organize into a higher-order structure.


I. INTRODUCTION
The field of active matter has seen significant development in recent decades. Seminal studies include the finding of the ordered state in the Vicsek model [1], the theoretical study on Purcel's swimmer [2], and the experimental realization of active colloids [3][4][5]. Following these, a number of studies have found that the active system can self-organize into various patterns [6,7]. In addition, many studies, such as the ones on suspensions of microorganisms [8,9], examined the rheological properties, that are affected by coordinated activity. These studies suggest that the pattern formation observed in active matter systems may build a basis for complex emergent functions [10].
Among them, the actomyosin cortex is one of the most basic structures, as it maintains the cellular shape. It is a network underneath the cellular membrane, consisting of actin filaments and myosin motors as well as other as- * tarama.mitsusuke@phys.kyushu-u.ac.jp sociated proteins. Thus, to form the actomyosin cortex, these molecules first need to accumulate in the membrane vicinity. Interestingly, such accumulation of active matter under confinement has been reported in recent studies using very simple models of active colloids [49,50] and active rods [51,52]. The phenomenon has been compared with motility-induced phase separation [53,54]. However, while actomyosin and active colloids show similar active dynamics, there are significant differences between these systems. Firstly, unlike active colloids, actin filaments do not self-propel unless they treadmill (an effect not considered in our study) [55]. Secondly, actin filaments are very long compared to the typical active colloidal rods. They are filamentous molecules of about 10 nm in width and a few tens of nm to several µm in length [56]. In addition, their persistence length is about 10 µm, comparable to the size of a cell, which also makes them quite stiff. Generally, such a long stiff object experiences a strong depletion force when confined due to the resulting restriction of its rotational degree of freedom. Therefore, in the case of these long stiff objects that only becomes active when bound by the motors that produce force on them, accumulation at a confining membrane is both non-trivial yet highly relevant for biological systems.
The purpose of this study is to investigate pattern formation of an actomyosin network that interacts with a confining membrane and to analyze the mechanical function of the resulting self-organized structure. In particular, we focus on the competition between the depletion force and the motor active force. Previously, the dynamics of actomyosin has been studied by using a macroscopic continuum model referred to as active gel model [57][58][59]. However, in order to bridge the microscopic molecular information and macroscopic structures, we develop a mesoscopic molecular dynamics model, in a similar spirit of those in Refs [60,61], in which the molecular origin of the macroscopic mechanics were analyzed. This paper is organised as follows. We first define the coarse-grained molecular dynamics model of cytoskeleton in the next section, which is followed by the explanation of the simulation results in section III. Finally, Sec. IV is devoted to the summary and discussion.

II. MODEL
We start by defining our coarse-grained molecular dynamics model of actin filaments and myosin motors. We model a filament by a linear series of discrete particles connected by elastic springs and a motor by a particle with two heads (Fig. 1). The two motor heads can bind to two different filaments and walk actively along them. This reproduces the motion of a bipolar motor such as the non-muscle myosin II minifilaments that possess multiple actin binding sites on both ends [62]. The equations of motion for the filament and motor particles are given by Here r f f,i represents the position of the ith particle of filament f . The index f is omitted hereafter for ease of notation. r mot m is the position of the motor m, and the coordinate 0 ≤ α mot m,h ≤ 1 along the filament segment r f ir f i+1 gives the position of the hth head of motor m as (Fig. 1). Then, the equation of motion for the motor head reads where r f i(i+1) = |r f i+1 −r f i | is the filament segment length. See Appendix A for the derivation. Inertia terms are neglected because of the small size and velocity of the molecules; The typical sizes of the cytoskeletal filament and the motor are less than a few µm and submicrometers, respectively. The sliding velocity of the motor head is submicrometers per second. Note that the filament polarity is taken into account by the order of the filament particle indices, which gives the cue for the motor  heads to walk actively towards the filament plus end, i.e., in increasing index i. The filament and motor particles experience friction and fluctuation from the surrounding cytosolic fluid, which satisfy the fluctuation-dissipation theorem; ξ X i = 0 and ξ X i,a (0)ξ X j,b (t) = 2γ X δ ij δ ab δ(t) with X = {f, mot}. The cytosolic friction is defined using a cylindrical approximation [63] as γ f = 3πη cyto (3d f + 2 f 0 /5) and γ mot = 6πη cyto (3d mot + mot 0 /5), where η cyto is the cytosolic viscousity and d f (d mot ) and f 0 ( mot 0 ) are the width and length of the filament segment (motor), respectively. On the other hand, since the motor heads are bound to and walk along filaments due to the motor walk force (MWF) f walk , they experience a sliding friction ζr f i(i+1) dα mot m,h /dt. The immobile membrane is also discretized by small particles [64], which confine the filament through the repulsive interaction with U rep (r; , σ, r * ) = 1−σ/r * (σ/r − σ/r * ) H(r * − r). rep , σ rep , and r rep * are the intensity, decay length, and cutoff length of the repulsive interaction. r memb k is the position of the membrane particle closest to the filament particle r f i . H(x) is the Heaviside step function that takes the value 1 for x > 0 and 0 otherwise.
The length and straightness of the filament are maintained by the stretching and bending elasticity acting between the filament particles: where U harm (r; κ, 0 ) is a harmonic potential with elastic modulus κ and rest length 0 , and U bend (r − , r + ; k) = −kr − ·r + is bending energy with bending rigidity k f,bend .
Here, we use the abbreviations r = |r|, r f ij = r f j − r f i , and Similarly, the stretching and bending forces act between the motor particles and motor heads: , is assigned to the filament particles with the geometric weight α mot m,h , which is necessary because the binding position may be in between the filament particles (Fig. 1b). Therefore, the force on the filament particle i from the motors reads Here the first and second summations are calculated over the motor heads bound to the filament segments r f i -r f i+1 and r f i−1 -r f i , respectively. This ensures force and torque conservation of the filament segments. Note that the counter force of the motor sliding friction and MWF, which act on the filaments, vanishes because of the force balance equation (3). Therefore, Eqs. (1)-(3) statistically satisfy the force-and torque-free conditions, which are required for active systems including migrating cells [12,13].
In addition to the equations of motion (1)-(3), we consider the following stochastic processes. The first one is the process of motor binding and unbinding. Unbinding of a motor head occurs in three cases. Firstly, a motor head that is bound to a filament unbinds stochastically at the rate of ω mot TO . Secondly, when a motor head reaches the end of a filament by actively walking along it, it unbinds instantaneously with the probability unity. Thirdly, a motor head unbinds with the probability unity when the other head of the motor unbinds. Here we assume that a motor takes either the bound state where both of the heads are bound to filaments or the free state where neither of the two heads are bound to filaments for simplicity. The free motors are assumed to diffuse sufficiently fast, so that they distribute uniformly. Then, a free motor binds to a randomly-selected pair of filaments whenever they find a position on each filament that are separated by the length of the motor mot 0 . Thus, this motor binding process alone causes no force on the system. The active force is generated in the cytoskeleton through the motor stretching and bending energies that are stored when motor heads move along filaments. The second stochastic process is filament turnover that takes place at a rate of ω f TO . The filament undergoing turnover is placed back into the system immediately at random position with random orientation, and all the motors previously bound to it become free. Although this model can be applied to both actin filaments and microtubules, in this paper, we focus on the system of actin filaments and myosin motors with an appropriate choice of the relevant parameters as summarized in Appendix B.
We solve the set of time-evolution equations in the following manner; First, we calculate the turnover of the filaments and motors, and the position of the motor head is updated by solving Eq. (3) with the Euler method. Then, Eqs. (1) and (2) are solved by using the fourthorder Runge-Kutta method. In the following, energy, length, and time are rescaled by using thermal energy U 0 = k B T , motor length l 0 = mot 0 , and motor cytosolic friction coefficient t 0 = γ mot ( mot 0 ) 2 /k B T . Finally we comment on the difference of this approach from existing models. There are several open source packages based on similar molecular dynamics models, including AFINES [65], aLENS [66], CyLaKS [67], Cytosim [68], and Medyan [69]. In these models, a motor is reduced to a harmonic potential with zero rest length or it is represented as a single segment with finite length since myosin motors often form minifilaments. In constrast to these pre-packaged simulation approaches, a motor in our model is modeled by two segments, which connect the motor particle with the two heads. This treatment is very similar to those in Ref. [63,70] and allows to include both the finite length of the motor and its cytosolic friction, which is mechanically more relevant to the actual situation since the cytosolic friction depends on the length of the motor.

III.1. Accumulation in the membrane vicinity
First, we consider actomyosin dynamics inside a circular membrane. Without motors, the filaments are depleted near the membrane. The depletion zone is also observed when motors are introduced but no MWF is applied (Fig. 2a). For finite MWF, however, the filaments accumulate in the vicinity of the membrane, forming a structure that resembles the actomyosin cortex (Fig. 2b). Here, the membrane curvature is set to κ = 0.05, corresponding to a radius R = 20 (4µm), which is twice the filament length L f = 10 (2µm).
To investigate whether the observed accumulation is universal or a simple consequence of the confinement curvature, we varied the curvature of the membrane. Interestingly, the filament accumulation is found for planar membranes with zero curvature and even for an oppositely curved circular membrane with positive curvature, in which case the filaments are present outside the membrane (Figs. 2c and 2d).
To quantify this accumulation, we calculated the density distribution of filaments (DDF) ρ f defined as a function of the distance from the membrane d memb by where the function ] from the membrane and equals 0 otherwise: is the area of the region within the distance from the membrane, and d d gives the width of the average region, which we set to d d = f 0 /2. x t represents the time average. Since the filaments can bend, we calculate the DDF using the distance of all filament particles from the membrane where n k is the normal direction to the membrane.
The DDF reaches zero at the membrane and increases with d memb (Fig. 2e). When MWF is absent, ρ f almost vanishes at small d memb showing a peak at larger distance, which represents the fact that the filaments are depleted from the membrane vicinity. When MWF is switched on, however, the peak in ρ f appears closer to the membrane, as marked by the arrow heads. This indicates the filament accumulation in the membrane vicinity.
At equilibrium, filaments are depleted from the membrane vicinity where their rotational degrees of freedom are restricted (See Appendix C). This is also true when motors cross-link filaments but lack the active force generation (Fig. 2a). In order for filaments to accumulate in the membrane vicinity, MWF needs to overcome the depletion force. To quantify this point, we plotted the DDF peak distance,d memb , as a function of MWF. Althoughd memb also depends on the number of filaments and motors, we found that all data collapse on top of each other when the MWF is rescaled to the average MWF F walk = f walk N m /N f (Fig. 3a). This indicates the universality of the observed phenomena. As F walk increases, d memb shifts towards the membrane and becomes smaller than L f /2, at which the depletion force sets in, and thus, d memb < L f /2 corresponds to the filament accumulation in the membrane vicinity. We define the transition point by F walk = F walk * that givesd memb = L f /2. Note that the critical average MWF is smaller for larger membrane curvature κ (Fig. 3c). Presumably this effect occurs because the rotational degree of freedom of the filaments increases with the membrane curvature, which weakens the depletion force (See Fig. 7c in Appendix C). In addition, we note that,d memb in the limit of large F walk converges to a constant valued ∞ memb , which is smaller than L f /2 (Fig. 3d). In conclusion, MWF exerted on the cytoskeleton is able to overcome the depletion force, leading to a filament accumulation in the membrane vicinity that resembles the actomyosin cortex.

III.2. Pressure
One major function of the cytoskeleton is force generation in the cell. Thus, we are interested in the pressure that the actomyosin produces on the membrane as a function of MWF and calculate it by where A memb is the length (area in 3d) of the membrane. The pressure depends on the number of motors and it scales with the number of filaments N f (Fig. 3b). In fact, all the data of the scaled pressure P = P/N f collapse except for the intermediate F walk . Peculiarly, in the intermediate regime of F walk , the scaled pressure P takes a larger value for a smaller number of motors N m /N f , although the data of P collapses for each ratio N m /N f (Fig. 3e). Another interesting point is that the scaled pressure converges to a constant value P ∞ for large F walk . P ∞ decreases as the membrane curvature increases (Fig. 3f). This convergence means that the active force that the motors generate on the actomyosin accumulation in the membrane vicinity is not directly converted to the pressure on the membrane.

III.3. Filament structure in the accumulation
To understand the reason why the pressure on the membrane become a constant for large MWF, we investigate the structural order of filaments inside the accumulation. To this end, we measure the polar and nematic order of the filaments with respect to the membrane normal direction in the membrane vicinity. The polar and nematic order parameters of the filaments in the mem- where, by considering the relevant neighborhood, the cutoff distance is set tod memb ≤ 1.25, which is about half ofd ∞ memb . The magnitude and direction of the polar and nematic order parameters are obtained by using the relationship P = p(cos φ, sin φ), N = n cos 2θ sin 2θ sin 2θ − cos 2θ .
Firstly, the magnitude of the polar order p decreases as F walk increases (Fig. 3g), with its direction perpendicular to the membrane φ ≈ 0 as shown in Fig. 3h. Secondly, the nematic order n increases for F walk F walk * (Fig. 3i), with its direction parallel to the membrane θ ≈ π/2 as depicted in Fig. 3j. These results indicate that the filaments approach the membrane as a result of MWF for F walk F walk * , whereas they tend to align nematically parallel to the membrane for F walk F walk * . This high nematic order of filaments aligning parallel to the membrane is the reason whyd memb and P converge to constant values for large F walk since the parallel filaments can slide along the membrane without pushing it. This also explains why P ∞ decreases as κ increases, since for larger membrane curvature the filaments nematically aligned parallel to the membrane have more chance to point away from the membrane (See Fig. 7 in Appendix C). Note that the high values of n for

III.4. Effect of filament length
In this section, we study the impact of the filament length. We performed simulations with shorter (L f = 8) and longer (L f = 13) filaments than those used so far (L f = 10). The results are summarized in Fig. 4. The data collapse of both the DDF peak distanced memb and the scaled pressure P is unaltered by the filament length. In particular, the odd increase in P for smaller motor number at intermediate F walk is also observed. These results are evidence of the universality of the observed phenomena in the cytoskeleton active dynamics.
For even shorter filaments, the collapse of data became worse, although the filament accumulation in the membrane vicinity is still observed for large MWF as shown in Fig. 5. One possible reason is that the cytosolic friction decreases with filament length, and thus, the cytosolic fluctuation increases. In all cases, however,d memb as well as P still converge to constant values (d ∞ memb and P ∞ ) for large F walk .

III.5. Filament motion
Finally, to understand how the filaments move and accumulate in the membrane vicinity, we plot the time evolution of the polar and nematic order parameters, and the ratio of the tangential and normal speed and velocity (i.e., signed speed) with respect to the membrane in Fig. 6.
Starting from a uniform distribution, filaments first form polar order in the bulk, i.e., far from the membrane (purple region around t 20 in Fig. 6a). This drives the filaments towards the membrane with a velocity given by the perpendicular component of the speed |v ⊥ | . At the same time, the filaments also move in the parallel direction to the membrane as indicated by the high value of the tangential speed |v | . In fact, these two components are comparable in strength (with the tangential component slightly higher than the normal one), causing a diagonal motion. In contrast, the tangential component of the signed velocity almost vanishes indicating a tangential motion of the filaments in two opposing directions, while the perpendicular component of the signed velocity again indicates the net filament motion towards the membrane (blue region in Fig. 6h corresponding to the purple region in Fig. 6a). This motion leads to the accumulation of the filaments in the membrane vicinity.
When the filaments accumulate in the membrane vicinity at the later stage t 20 , the polar order decreases and the nematic order appears. The filaments in the accumulation slide along the membrane as indicated by a finite value of their tangential speed. Note the similar motion was reported for self-propelled rods [51] although the propulsion mechanism is different from actomyosin. Since the filaments form nematic order, however, the tangential motion vanishes on average ( v ≈ 0) due to the existence of the counter-moving filaments.
Finally, we estimate the width of the high nematic order region in the membrane vicinity. To this end, we measured the maximum nematic order parameter at each time after the system has reached a steady state (for t ≥ 40). The time average of the maximum nematic order parameter is n max ≈ 0.75. We define the high nematic order region in the membrane vicinity as the region with n ≥ n max /2. The boundary with n = n max /2 is displayed by the gray solid line in Fig. 6c. We estimate the width of the high nematic order region as the time-average of the boundary distance d memb . The obtained length scale is d memb ≈ 5.12, which is plotted by the black dotted line in Fig. 6c. This width is approximately given by half the filament length b memb = L f /2, at which the filaments start to experience the depletion force. This agrees with the picture that the filaments are driven by the bulk polar order to accumulate in the membrane vicinity, and then exhibit a nematic order due to their interaction with the membrane.

IV. DISCUSSION
To summarize, we demonstrated that the competition between the depletion force and MWF leads to an accumulation of filaments in the membrane vicinity, which resembles the actomyosin cortex. The self-organized structure exerts pressure on the membrane, which converges to a constant value for large MWF because of nematic alignment of the filaments parallel to the membrane. Interestingly, for intermediate MWF, the pressure increases when the number of motors is decreased. We highlighted the universality of the phenomenon by showing the data collapse of the DDF peak position and the scaled pressure (P = P/N f ) as functions of the average MWF F walk = f walk N m /N f . Qualitatively the same results are obtained for different membrane curvatures. Moreover, the results are quantitatively unaltered by the filament length, except for very short filaments for which the pressure profile changes while the structure formation is still observed. Our results provide a novel insight into the self-organization of cytoskeleton active matter into higher-order structures and its emergent mechanical function.
Our model is relatively simple, such that our findings can in principle be tested against experiments on artificial cells [71,72]. However, to our knowledge, experiments on the self-organization of actomyosin cortex in artificial cells remain challenging [73], and oftentimes depletion agents are used instead to study the dynamics with the cortex. In comparison to the actomyosin cortex of real cells, our model omits many elements. For instance, actin nucleators and severing proteins are not included, which are thought to play a major role in the cortex formation [74]. Instead, we simplified the polymerization/depolymerization processes with a stochastic turnover. As a result, the self-organized structures tend to fade away for very large turnover rates, leading to a uniform distribution of actin filaments.
Recent studies reported that contraction occurs due to additional passive crosslinkers between filaments without force generating activity [75][76][77]. This effect also appeared in our simulation when crosslinkers were introduced, but was observed at positions both close to and far away from the membrane starting form a uniform distribution. This contraction occurs even without the membrane. Therefore, to reproduce a contractile actomyosin cortex in simulation, both motors and crosslinkers need to be introduced to a filamentous network that is enriched in the membrane vicinity due to membrane-associated nucleation and biased polymerization. To maintain the localization of the contractile actomyosin cortex in the membrane vicinity, additional linkers between the filaments and membranes are necessary.
Finally, the extension of our system to three dimensions is important. In 3d, we expect defects to appear in the accumulation, which move dynamically. Such defects were previously found in the experiment using microtubules [78]. Further, several studies using epithelial monolayers, another type of active system, focused on the dynamics of the defects, which causes cell exclusion and cell death due to the local stress [15]. Therefore, it would be interesting to see if the defects in the accumulation of actomyosin also leads to the local change in the pressure. The position of the motor head h (= 1, 2) of motor m bound to a filament segment connecting r f i and r f i+1 is given by where 0 ≤ α m,h ≤ 1 measures the relative position along the filament segment. This allows us to take into consideration the situation where the motor head is bound to a position in between the filament particles. Then, the time derivative of Eq. (A1) consists of two terms: The first term represents the sliding motion of the motor head along the filament segment, whereas the second term is the velocity of the filament at the motor head position. The latter corresponds to the transport of the motor head caused by the translation of the binding filament. Note that the left-hand side of Eq. (A2) is the velocity of the motor head measured in the Lab frame.
Since we assume that the heads of the bound motors are always on the filaments, we only need to consider the time evolution of α m,h to specify the position of the motor heads. The motor heads experience motor stretching and bending force as well as a motor walk force (MWF) f mot,walk m,h . We assign the former two to the particles of the filament segments to which the motor heads bind, to keep the motor heads on them. Then, the time-evolution equation of the motor head h is given by where ζ is the sliding friction coefficient between the motor head and filament. The molecular motor has an ability to walk actively along the filament by consuming chemical energy in the form of ATP. This effect is included in the model as the MWF The sign specifies the direction of the motion; the plus (minus) sign corresponds to active walking towards the plus (minus) end of the filament. The direction of the motor walk is characteristic to the type of motor. In the case of the non-muscle myosin II motors that we consider in this study, the walking takes place towards the plus end of the actin filaments: f mot,walk . Here, we make a simplification on the process of motor head stepping forward along a filament, which is modelled by a constant walk force, since we are interested in the long time scale dynamics. This simplification allows us to neglect the transient unbind of the motor head while stepping forward. Note that, since the motor heads are moving along the filament, MWF acts in the direction parallel to the filament segments.
The force from the motor head on the filament is given by The first and second terms represent the counterpart of the sliding friction force between the motor head and filament and the MWF, respectively, so that the law of action and reaction is satisfied. The last two terms are the motor stretching and bending force acting on the motor head h, which is assigned on the filament particles. From Eq. (A4), however, the first two terms on the righthand side of Eq. (A6) cancel and the remainder reads Since the filament is modeled by discrete particles, this force should be split to the two edge particles (r f i and r f i+1 ) depending on the geometric weight α m,h as See Fig. 1b. Note that α m,h measures the position of the motor head bound to in between the filament particles. This assignment of the force ensures the conservation of the force and torque acting on the filament segment due to f mot (r mot m,h ). If there is more than one motor head on the filament segment r f i -r f i+1 , then the right hand side of Eq. (A8) is summed up over them. As a result, the filament element i experiences the force from the motors given by where the first and second summations correspond to the contribution from the motor heads on the filament segments r f i−1 -r f i and r f i -r f i+1 , respectively. Here, M i represents the list of the motor heads bound to the filament segment r f i -r f i+1 .

Appendix B: Parameter values
The simulation parameters for the filaments and motors are set as summarized in Table I, unless otherwise stated. Some of the parameters are compared with existing experimental measurements. The bending rigidity of the motors is set to the value same as that of the filaments, since, to our knowledge, it has not been measured for myosin minifilaments yet. The parameters of the confining membrane and system size are set as follows. For the circular membrane with negative curvature, the radius is set as R = 20, which leads to the curvature κ = −0.05. For the plane membrane with zero curvature κ = 0, the distance between membranes is L x = 30 and the length of the membrane L y = 41.75 with periodic boundaries. For the circular membrane with positive curvature, in which case the filaments are present around the membrane, the radius is R = 20, and the width and height of the simulation box is set L x = L y = 70.9 with periodic boundaries. This leads to the curvature κ = 0.05.

Appendix C: Depletion force
To provide an intuitive idea of the depletion force, we show schematic sketches depicting a filament placed far from and close to the membrane in Fig. 7. A filament placed far from the membrane can rotate freely due to thermal noise (Fig. 7a), whereas the rotational degree of freedom of a filament placed close to the membrane is restricted (Fig. 7b). This entropic penalty tries to keep the filament away from the membrane, giving rise to a depletion force. The threshold distance at which the filament starts to experience the depletion force is approximately given by half of the filament length d memb = L f /2. The cutoff distance of the repulsive interaction between the membrane and filament particles lies within the threshold distance in this study.
In Fig. 7c, we show the effect of the membrane curvature on the depletion force. As the membrane curvature decreases, the filament rotation is more restricted. Therefore, the filaments experience a stronger depletion force for decreasing the membrane curvature. As a result, in the membrane vicinity, filaments originally placed parallel to the membrane at the same distance are more likely to point away from the membrane due to fluctuation as the membrane curvature increases.