Pattern formation in active cytoskeletal networks

We present a one-dimensional model combining two of the main features of active biopolymer solutions that have only recently been investigated experimentally, namely the molecular motor driven active transport of filaments and the (visco-)elastic properties of filament networks held together by crosslinkers or entanglement effects. It is shown that the pattern forming mechanisms, associated with the motor-mediated transport of filaments, are substantially altered when coupled to a filament network: in the case of a permanent network, the long-range clustering of filaments changes either to stationary periodic filament density patterning or to propagating pulses. If the network is however viscoelastic, molecular motor activity can lead to traveling or standing filament density waves.


Introduction
The cytoskeleton of eukaryotic cells is constituted of filaments, actin and microtubules, and relies vastly on its ability to actively reorganize itself by polymerization processes and motor protein complexes [1,2]. The latter active, nonequilibrium process has become a rapidly growing field of research: in reconstituted filament-motor bundles, contracted states and spontaneous oscillations have been found [3]- [5]. In extended systems and in the dilute limit, the intermittent motor mediated filament transport leads to spontaneous self-organization [6]- [11].
A second important feature of biopolymers is that upon increasing filament densities, elastic networks and gels may be formed either by entanglement effects or by means of permanently crosslinking proteins. During recent years continuous research efforts, both on the experimental and theoretical side, have been devoted to the elastic and dynamic properties of such viscoelastic networks without motor activity [12]- [17]. Lately the interplay between active motor transport and an elastic network has started to be experimentally addressed [18,19], as well as the influence of a trace amount of crosslinks (not sufficient to build up an interconnected network) on the pattern formation processes [20]- [22]. Whereas the motor-mediated transport of filaments as well as passive networks have already been addressed separately to quite some extent, theoretical efforts to join both vitally important features have been started only recently [9,23]. Understanding such active crosslinked gels would not only be of obvious biological interest, it also comprises a hitherto unexplored material class, that could be called an 'active elastomer': liquid crystalline filament ordering is known to appear in biopolymer solutions at sufficiently high densities [24]. Crosslinking a dense ordered filament solution to a gel results in a passive elastomer that has already been recently discussed in the case of actin in [25]. Even more tantalizing is that the temporary action of molecular motors in

A model for an active filament network
A solution of cytoskeletal filaments forms networks in the presence of crosslinkers 3 , which are actively reorganized and set under stress in the presence of motors. To get a treatable model for such a complex system, several approximations have to be made. One assumption we use here, is that the average crosslinking time is large compared to the average runlength of the motors on filaments. In this limit, it is appropriate to divide the total filament density into a fluid part, made up of filaments not in permanent contact with the network and either moved by molecular motors or randomly (referred to as 'free' filaments hereafter) and into an elastic part supposed to be formed of crosslinked filaments, which form a (visco-)elastic network ('network' filaments). This approach is related to two-fluid models [27,28], well known for complex fluids. To overcome the limit of permanent crosslinkers, one can account for transition rates between free and network filaments, as will be done in the future. Such transition kinetics is known to favor spatially periodic patterns in related systems, e.g. self-assembling biofilaments [29]. For a one-to-one experimental realization of this model, one can imagine an in vitro system where crosslinking times much longer than in vivo are achievable, or even use artificial crosslinks.
A second approximation made-which however has no effect on the simplest case of a perfectly polar bundle solely investigated in this work-is that the network itself remains perfectly isotropic. The total filament system (i.e. free and network filaments together) may however become anisotropic upon the reorganization of the free filaments. We discuss this approximation in more detail below.
To highlight the new properties due to the (visco-)elastic background we will consider the simplest case of a 1D system, i.e. an active filament bundle structure. We use an approach commonly accepted in the literature for 1D motor-filament solutions without (visco-)elastic background [4,5,30] 4 . Thereby the free filament bundle axis is referred to as the x-direction. Since free filaments are polar with respect to the action of molecular motors, an active free filament bundle is distinguished by number densities c + (x) and c − (x) of filaments with their plus ends pointing either in the positive or negative x-direction. As the (visco-)elastic network is as well composed of filaments, it locally can be described by a second class of filament number densities, d ± (x). While motor-induced interactions among network filaments are neglected due to their large meshsize, motors may couple network filaments to free filaments (just as they couple free filaments). Therefore free filaments can be actively transported along the network in addition to active transport of free filament pairs. The d-filaments will be used to define the interactions between the free filaments and the network. On large scales, however, the network dynamics is best captured by a 1D displacement variable u(x). To get a closed set of equations the spatial modulations of the d-filament number density around its mean value d 0 will be expressed in terms of u(x), yielding coupled equations for the 'free' filament density profiles and the displacement variable.
Assuming filaments with a fixed length , the dynamics of the 'free' filaments is governed by conservation laws wherein D is an effective diffusion coefficient. Provided that two-filament interactions are prevalent (corresponding to a low motor density or duty ratio), the motor-induced currents incorporate both the interactions between free filaments and the interaction of a free and a network filament. This leads to the decomposition J ± = J ± 1 + J ± 2 + J ζ u , where J ± 1 stands for interactions between free filaments and J ± 2 for free filament/network interactions, respectively. As is derived later, J ζ u is caused by frictional effects the biopolymer network exerts on the free filaments. Both currents J ± 1,2 can be derived from micro-force balance, as has been done in [30] for the interaction of free filaments, which suggests the currents J ± 1 = J ±± cc + J ±∓ cc to be Here α and β are averaged cross-sliding velocities characterizing the interaction strength of parallel and antiparallel filaments. On larger scales, the α-contributions favor stationary contracted states of filaments, while the β-contributions lead to a phenomenon called polarity sorting, and also to propagating solitary density modes [4,5,31]. Since both the free filaments are exposed to the same local restoring force-namely the friction with the surrounding fluidthe center of mass of a free filament pair is invariant. This however does not hold if a network filament interacts with a free filament: different restoring forces, the aforementioned frictional force and the elastic contribution in the case of the network filament, lead to an effective barycentric transport, as can be seen from figure 1. The displacement of the center of mass is equivalent to an additional current J χ u , yielding are the natural extensions of equations (2) for motor-mediated interactions between d ± -and c ± -filaments. The current J χ u is given by It makes allowance for the fact that the viscoelastic force density [(λ + 2µ) + τ ∂ t ] ∂ 2 x u enters the microscopic force balance, when an active overlap between a network filament and a free filament is met. This can also be derived from microscopic force balance as is made plausible in the appendix. Another aspect which has to be taken into account is the frictional interaction between free filaments and the network as a whole if the network slides by at a velocity ∂ t u(x). This is done by the contribution with ζ = ζ /η. The dynamic equations for the free filaments then read Figure 1. (a) Schematic representation of the forces exerted by a motor complex between a free filament (light gray shaded) and a network filament (dark gray shaded) of the same orientation. (b) In contrast to interacting free filament pairs, the center of mass x CM of the free/network filament pair is displaced if a molecular motor has actively coupled the considered filaments. The microscopic forces depicted in (a) result in an effective barycentric transport of the free filament whose center of mass is shifted from x to x CM while the network filament's center of mass is hardly affected. The barycentric transport of the free filament that gives rise to the current J χ u strongly depends on the viscoelastic force density As aforementioned the dynamics of the (visco-)elastic network, locally constituted of the d ± -filaments, is best described by a displacement variable u(x). Force balance then implies ∂ x σ + f = 0 with σ an appropriate stress and f a force density yet to be defined. The force density f comprises the frictional force density −γ ∂ t u due to network frictions with the surrounding solvent, the drag force density f drag c accounting for the change in friction induced by the active free filament currents and the motor-mediated driving force density f u resulting from interactions between free and network filaments. This yields where γ denotes the viscous friction coefficient between the network and the surrounding solvent that in principle depends on the 3D morphology of the network but is considered 6 to be a constant in our 1D model 5 . The stress σ can be decomposed into σ = σ E + σ F . σ E stems from the elastic network deformations characterized by ∂ x σ E = (λ + 2µ)∂ 2 x u with Lamé coefficients λ and µ. τ being a friction coefficient, the frictional contribution σ F = τ ∂ x ∂ t u to the total stress σ accounts for frictions within the elastic network [32,33]. According to the Voigt-Kelvin body [34] the fraction (λ + 2µ)/τ then defines a time constant governing the formation/relaxation of strain within the network if an external force is applied/released. Network dynamics is governed by the following dynamical equation for the network's displacement variable u(x) In analogy to J ± , the force density f u originating from the motor-mediated interaction between free and network filaments reads Here η is the friction coefficient a free filament experiences if moved through the surrounding solvent. The free filament currents given by equations (2) acting on the network filaments alter the drag force −γ ∂ t u according to with ν denoting an adequate friction coefficient.
In order to get a closed set of equations, the spatial variation of the density d(x) needs to be interpreted in terms of the displacement field u(x). As the local density d(x) is increased in regions where the displacement field u(x) changes from high to low values and vice versa, spatial modulations of d(x) are proportional to ∂ x u(x) which allows the network density to be approximated (to leading order) by The parameter δ relates the spatial modulation of the network filament density d(x) around its mean value d 0 to the displacement field u(x).
In equation (11), we have made the approximation that the network itself remains perfectly isotropic (i.e. d + (x) = d − (x)). In the general 1D case, the local difference between the plus and minus filaments can be in principle nonzero, leading to an anisotropy. However, one has to keep in mind that we deal with a coarse-grained model, which is valid only on length scales above several filament lengths, and that the network is assumed to be tightly crosslinked, so that plus and minus filaments incorporated in the network cannot be moved independently as in the case considered originally in [5]. If the network is deformed and relaxes, due to the tight crosslinking of the filaments of the two possible directions, plus and minus, both locally undergo approximately the same displacements. In the 2D or 3D case, equation (11) has to be generalized. In 1D, anisotropy can be generated if-for a finite crosslinking time-transitions between network and free filaments are included, as will be discussed in future work.

7
To study now the essential effects of the network on the free filament dynamics, we consider in the following a perfectly polar bundle by setting The free plus-filaments with density c + (x) are denoted by c(x) hereafter. In the absence of a network the system is then known to trigger a subcritical, long-wavelength bifurcation to contracted states [5]. In the remaining sections, we show that both an elastic and a viscoelastic network significantly alter the pattern forming processes in motor-filament systems: the modified model already allows for oscillatory instabilities, even though there are no antiparallel interactions present, as well as for pattern forming instabilities with finite wavelength.

Patterns in an active elastic bundle
In case of an elastic network, i.e. τ = 0, the equations of the previous section reduce to These equations are rescaled by introducing dimensionless coordinates x = x/ and t = t D/ 2 , normalized density c + = c + and displacement variable u = u/ , as well as for the elastic modulus. The primes are omitted for the sake of simplicity.
The nonlinear integro-differential equations (12) are analyzed for a system of length L with periodic boundary conditions. The stability of the spatially homogeneous basic state, corresponding to a filament density c(x) = c 0 and a vanishing (or constant) displacement field u(x) = 0, can be analyzed by linearizing equations (12) around this state, c(x) = c 0 + k c k (x) and u(x) = k u k (x), with small perturbations c k , u k ∝ e σ t+ikx and k = 2π n/L (n ∈ Z) being a discrete wave number. Moreover c −k = c k as well as u −k = u k holds because c(x) and u(x) are real fields. The resulting set of homogeneous linear equations yields a quadratic polynomial in σ as a solvability condition. If the real part of at least one of these solutions σ 1,2 = κ 1,2 ± iω becomes positive for some wave number k, the basic state (c 0 , 0) T is said to be unstable. In the case that the imaginary part ω vanishes, the instability is stationary, whereas for finite ω an oscillatory bifurcation is met [35]. The motor parameter α is identified as the control parameter governing transitions from the homogeneous basic state to spatially and/or temporally modulated states.
The neutral curve α S (k) characterizing the stationary bifurcation is obtained by solving the neutral stability condition κ(k) = σ (k) = 0 for the control parameter α. Defining = ηd 0 /2 + νc 0 and = χ d 0 sin(k) − ηk/2, the instability threshold α S± c is given by the minimum of α S (k), which is located either at a vanishing critical wave number k S c = 0, thus implying a longwavelength instability, or at a finite wave number k S c that yields spatially periodic patterns. At the onset of the oscillatory bifurcation σ 1,2 = ±iω and the critical value α O c of the control parameter is determined to be the minimum of the neutral curve For smaller values of E an oscillatory long-wavelength instability is met, with the typical growth rate κ(k) depicted in (d). E-values beneath the solid line yield a long-wavelength, stationary instability whose wave number dependent growth rates are depicted in (e), an instability type which is known from oriented filament bundles without network. Within the area between the solid and dotted lines lying beyond the dashed line, a long-wavelength, oscillatory instability competes with a finite wave number stationary one, with growth rates displayed in (c).
with an appendant critical wave number k O c = 0. Herein is the neutral curve yielding the onset α F c ≈ 1/c 0 of the long-wavelength stationary instability in the absence of a network [30]. Within the limit (d 0 → 0) and disintegrating networks (i.e. E → 0), the stationary thresholds α S± c disappear whereas the onset α O c of the oscillatory bifurcation converges to α F c with the frequency ω simultaneously vanishing. Therefore one regains for vanishing network densities the well-known long-wavelength, stationary bifurcation, for motor activities α > α F c , characteristic of oriented filament bundles [30]. A typical bifurcation scenario in the (α, E)-plane is sketched in figure 2(a). Altogether the elastic network adds to the long-wavelength bifurcation at α F c three new instability types as evinced by the diverse growth rates κ 1,2 (k) plotted as a function of the wave number k in figures 2(b)-(e). Provided that the rescaled barycentric transport coefficient χ is large enough, the elastic network increases the instability thresholds of the homogeneous filament distribution towards motor activities α higher than α F c , depending on E, either beyond the dashed or the dotted line in figure 2(a). If χ is too small the solid line abuts the threshold α F c and the only instability left is the long-wavelength stationary one, known already from oriented bundles without network elasticity [30], cf figure 2(e). The same allegation holds for small network elasticities. Figure 3 shows numerically obtained density profiles c(x, t) of the free filaments beyond the respective instability thresholds as a function of x for a fixed system length L. The diversification of the latter instability by the network is hence mostly driven by the barycentric transport imposed by the network upon the free filaments. If filament bundles with both orientations, i.e. comprising c + -as well as c − -filaments, are considered instead of oriented ones, the barycentric transport coefficient χ becomes less relevant as the possibility of polarity sorting introduces an additional microscopic transport mechanism captured by the motor activity β = 0. In the case of very stiff networks an increase in α leads to a shortwavelength, stationary instability whose wave number-dependent growth rates are shown in figure 2(b). Within the nonlinear regime, this instability leads to stripe patterns as shown in figure 3(a). Looser polymer networks, corresponding to lower E-values, favor a transition from the homogeneous filament distribution to a long-wavelength, oscillatory mode characterized by the growth rate κ(k) in figure 2(d) and rather irregular dynamics as depicted in figure 3(c). This oscillatory instability can compete with a short-wavelength stationary one, whose associated growth rates are displayed in figure 2(c). Within the nonlinear regime free filament density peaks, resulting from a coarsening process at an early stage, oscillate between extremal positions as shown in figure 3(b). This type of nonlinear behavior occurs in a large region of the (α, E)plane bounded by the solid, dashed and dotted lines. For small rescaled elasticities E the longwavelength, stationary instability characteristic of free oriented filament bundles prevails in the presence of a network. The appendant nonhomogeneous stationary numerical solution is pictured in figure 3(d).

Patterns in an active viscoelastic bundle
Crosslinks among network filaments generally have a finite lifetime. If their lifetime is of the same timescale as the active transport processes, the static network becomes viscoelastic, behavior that is incorporated in the model by a time constant (λ + 2µ)/τ as done in a simple Voigt-Kelvin model. The associated equations read wherein the frictional contributions ν and ζ have been omitted, a choice that does not alter the overall scenario. Moreover the viscoelastic contribution τ ∂ t u(x + ξ ) to the force density [(λ + 2µ) + τ ∂ t ] ∂ 2 x u entering the barycentric flux (4) of the free filaments will be neglected subsequently as it would only quantitatively modify the existence domains of the encountered instability types. The equations can once again be rescaled using the dimensionless coordinates of the previous section with the supplement τ = τ/(γ 2 ).
A linear stability analysis analogous to the one performed in the last paragraph reveals that the viscoelasticity of the network additionally triggers beyond some critical motor activity α W c , lying at the minimum of the neutral curve a bifurcation to traveling (TW) or standing (SW) free filament density waves with a critical wave number k W c . The corresponding generic growth rate is displayed in figure 4(d) as a function of the wave number k. Due to the network's viscoelasticity long-wavelength perturbations become always damped, provided that the elasticity E is large enough, as evinced by the growth rate shapes in figures 4(b)-(d). Additionally, as the total filament density-dx[c(x) + d(x)] = const-is a conserved quantity for the investigated model all growth rates are diffusive in the long-wavelength limit. This mechanism proves to be utterly important since the coarsening process associated with a long-wavelength instability is suppressed and finite wave number TW and SW patterns are stabilized. Here, the viscoelasticity of the network thus permits the filament-motor system to produce stable patterns. It is however known that long-wavelength modes yielded by conservation laws can render TW or SW states unstable in the strongly nonlinear regime as shown, e.g. in [36] for a model of mobile, charged ion channels embedded into a biomembrane.
Otherwise the phase diagram, figure 4(a), looks rather similar to the one discussed previously for elastic networks. The only differences are: upon increasing the motor activity α, the homogeneous filament distribution now firstly is unstable against the aforesaid shortwavelength, oscillatory instability for small E-values, whereas large rescaled elasticities still favor a transition to a finite wave number stationary bifurcation with modes growing according to κ 1,2 (k) shown in figure 4(b). A further increase in α leads to a coupling of the latter instability is the instability threshold without network. For stiff networks with large E the homogeneous state becomes unstable, beyond the dashed line, via a finite wave number, stationary bifurcation whose growth rates κ 1,2 are plotted in (b) as a function of the wave number k. In the case of smaller E the latter state bifurcates through a finite wave number, oscillatory instability characterized by a typical growth rate shape depicted in (d). In between the dotted, dashed and long-dashed lines the short-wavelength, oscillatory instability competes with a stationary one, setting in at finite wave number, as evinced in (c) by the according growth rates κ 1,2 (k). The remaining areas within the (α, E)-plane exhibit exactly the same dynamics as in case of a non-viscoelastic network in figure 2(a).
to the finite wave number oscillatory one as may be deduced from the growth rate depicted in figure 4(c). Large rescaled elasticities E and finite τ -values generally favor the Hopf bifurcation to traveling or standing free filament density waves, which survive in the nonlinear regime as can be seen from figure 5. Within the limit τ → ∞ the onset α W c of the latter wave patterns, visualized by the long-dashed line in figure 4(a), abuts the stationary threshold α F c characteristic of free filament bundles. In fact, for τ → ∞, the force applied by motor mediated interactions to the network filaments is dissipated by prevalent, frictional effects within the network. Thus within this limit network strain cannot build up and the free filament dynamics is hardly influenced.

Mechanosensitivity
Although the long-wavelength, stationary instability is already encountered in oriented fiber bundles lacking a (visco-)elastic network, where it is interpreted as a contracted state, it gains within the framework of the present (visco-)elastic model a striking new quality: its threshold, given by the minimum of equation (13), depends on both the mean network density d 0 and the network elasticity E. As does the stationary instability at finite wave number, the newly emerging oscillatory instabilities, that result from the coupling of the free filament dynamics to the (visco-)elastic network also depend on both parameters d 0 and E as can be seen from α O c . This stipulates that the aforesaid density instabilities can be controlled through compression and dilation-either externally by applying a force or internally by molecular motor activityof the elastic network or of the filament bundle in our 1D case. Compressing or dilating the network alters the mean filament density d 0 of the network, thereby changing the selforganization process of the free filaments and might be used to rapidly switch between the different instability scenarios discussed throughout this work. Addressing pattern formation by compression or dilation of a biopolymer network thus seems to be effectively equivalent to steering pattern formation by modifying motor activity through ATP concentration. It is known that cell movement needs adhesion to the substrate to transduce force, and it might be speculated that such external forces via the mechanism described above interplay with the internal cytoskeleton structure: periodic contraction waves in lamellipodia [37], as well as spontaneous oscillations in a muscle bundle have been reported [38], where the bundle's elasticity might be relevant too.

Conclusions
Numerous experimental and theoretical studies have been devoted either to the mechanical properties of filament networks or to the molecular motor triggered self-organization of filaments. In this work, we have investigated the influence of networks on the latter selforganization process by proposing a simple model that couples the active filament dynamics to either a purely elastic or a viscoelastic filament network.
The cytoskeletal structures that attracted most attention in recent years were the aster and vortex patterns found in dilute reconstituted solutions [6,7]. They have been successfully described as defect structures in the polar filament orientation field both by mesoscopic models in the dilute regime [11,39] and macroscopic hydrodynamics models [9,40]. To study the nonlinear dynamics of extended patterns and coarsening processes, the latter attempt proves difficult, since the generic approach is based on linear irreversible thermodynamics and yields only geometrical nonlinearities, while it is challenging to decide which additional nonlinearities are important [41]. Therefore, we start from a model for the dilute regime, based on [5], and extend it by accounting for the elastic background by means of motor-induced filament-filament interactions only, which naturally yields the leading order, quadratic nonlinearities.
While stationary and propagating contractile states were reported previously in suspended filament-motor systems [4,5,31], new instability types and nonlinear competitions between the latter ones emerge when introducing the (visco-)elastic network. Especially, in the viscoelastic case, TWs and SWs can be stabilized for specific elastic constants of the network. Moreover, all mentioned pattern forming processes are rendered mechanosensitive through the coupling to the network: the states as well as their onset crucially depend on the network density, elasticity and viscoelastic time constant, all these parameters being naturally relevant as a cell is subjected to external forces and continuously reorganizes its cytoskeleton.
Various generalizations of this model should be considered in the future: one relevant and interesting point is to account for a finite lifetime of the crosslinkers and describing the system at timescales much larger than the average crosslinking time. Firstly, the two-fluid separation has then to be generalized by transition rates between 'free' and 'network' filaments, naturally introducing the finite lifetime into the problem. As discussed above, this also allows the network to become anisotropic, even in 1D. Secondly, on these timescales also the network part of the model becomes a viscoelastic fluid rather than a viscoelastic solid and a Maxwell-like model, as suggested in the framework of [40], should be more appropriate.
Secondly, the current 1D model might be generalized to higher dimensions. The methods that can be used to achieve this are principally known [8,10,11]. In the case of quasi-permanent crosslinkers, this would also be the first step towards the active elastomer, since the broken symmetry variable of the polar filament orientation should nontrivially couple to the (visco-) elastic field [42]. In forthcoming work, we will also show that pattern formation can be conceived as an instrument to stiffen (visco-)elastic filament networks. with F act,a and F act,b being the active forces exerted by the motor on the two filaments, respectively (see figure A.1(a)). Since no external forces are involved, the center of mass of the two filaments is invariant, leading to F act = F act,a = −F act,b . Each force exerted on one filament by the motor is balanced by an opposing force acting on the other filament and both filaments cover the same distance x. Equation (A.1) then leads to J ++ cc = 1 η F act c + (x, t), with F act being the representation of F act in the mean field. Given that an overlap of the two filaments is a precondition for the active interaction and that due to force balance the forces F act have to obey distinct symmetry relations, the currents are given by equation (2a). For a more detailed and rigorous derivation we refer to [30].
However, the situation is thoroughly changed if a free filament interacts with a network filament, leading to as the equation of motion for the network filament, with F elast being the elastic restoring force. Consequently the distance x covered by the free filament is larger than in the case without network ( x), whereas the distance x NW the network filament covers during an active interaction is smaller resulting in an effective barycentric transport (see figure A.1(b)). Since the displacement of the free filament and the network filament are counter-directed, x − x NW = 2x (A.4) applies here. Taking the derivative with respect to t on both sides and inserting equation (A.2) and F act,a = −F act,b from above, gives an equation forẋ which can be used to calculate the associated current While in the mean field the current J ++ dc again takes the form of equation (3a), the current J χ ,++ u is due to the viscoelastic force density entering the force balance. It virtually transmits the viscoelastic force density to the free filament. Considering a network plus-filament at x + ξ and a free filament at x this gives rise to as the current associated with the barycentric transport.