Active Curved Polymers form Vortex Patterns on Membranes

Recent in vitro experiments with FtsZ polymers show self-organization into different dynamic patterns, including structures reminiscent of the bacterial Z-ring. We model FtsZ polymers as active particles moving along chiral, circular paths by Brownian dynamics simulations and a Boltzmann approach. Our two conceptually different methods point to a generic phase behavior. At intermediate particle densities, we find self-organization into vortex structures including closed rings. Moreover, we show that the dynamics at the onset of pattern formation is described by a generalized complex Ginzburg-Landau equation.

Intracellular structuring is often facilitated by the active dynamics of cytoskeletal constituents. The origin of these driven dynamics and their impact on pattern formation has been extensively studied using artificial motility assays of cytoskeletal filaments [1][2][3][4]. Another intriguing example of self-organization due to driven filaments was reported recently by Loose and Mitchison [5] In vitro, the bacterial protein FtsZ forms membranebound, intrinsically curved polymers. These seem to exhibit treadmilling dynamics (consuming GTP) and, as a result, move clockwise on the membrane. Depending on the protein density, polymers cluster into dynamic structures such as rotating rings or jammed bundles, despite the absence of attractive interactions [6]. These ring structures are of particular interest, since in vivo, FtsZ builds the contractile Z-ring which drives cell division in a yet unknown way [7][8][9]. But also in the in vitro experiments, the pattern forming mechanism remain unclear even on a qualitative level.
Motivated by these experimental findings, we have studied pattern formation in a class of active systems, where particles move on circular tracks and interact only via steric repulsion. To assess the dynamics of this class, we consider two conceptually different models: First, we emulate active particles as elastic polymers with fixed intrinsic curvature that move with a constant tangential velocity [ Fig. 1(a)] and perform Brownian dynamics simulations. Second, we employ a kinetic Boltzmann approach, where point-like particles move on circular paths and undergo diffusion and binary collisions (with polar symmetry) according to a simplified collision rule [ Fig. 1(b)]. As a result, we identify different phases of collective behavior as a function of density and noise level. With both approaches, we find flocking into vortex patterns in the regime of intermediate density and noise strength. Our simulations for extended particles predict the formation of closed ring structures reminiscent of those found in Ref. [5], even in the absence of any attractive interactions. In the mesoscopic limit, our analysis yields that, close to the onset of vortex formation, the dynamics at onset of ordering is characterized by a novel generaliza- tion of the complex Ginzburg-Landau equation.
In our Brownian dynamics simulations we consider a system of M curved polymers of the same chirality embedded in a two-dimensional membrane of area A with periodic boundary conditions. Each polymer is described as an inextensible worm-like chain [10,11] of length L, persistence length p , and intrinsic curvature κ 0 . For a given polymer conformation r(s), parameterized in terms of arc length s, the overall bending energy is given by E bend = 1 2 p k B T L 0 ds [κ(s)−κ 0 ] 2 , where κ(s)=|∂ 2 s r(s)| denotes the local curvature. Excluded volume interaction is implemented by a repulsive truncated Lennard-Jones potential (for details see the Supplemental Material [12]). To assure motion of the filament contour on a circular track (apart from noise), polymers are propelled with a tangential velocity v 0 (s)=v 0 ∂ s r(s). This accounts for the effective motion of treadmilling in a simplified way [12]. Note that for this choice, the area explored by a circling polymer is minimal. In the free draining limit, the dynamics of the polymer system is then determined by a set of coupled Langevin equations for the contours r (m) (t, s) of each polymer m=1, 2..., M : balancing viscous friction with elastic and repulsive forces generated by the total energy E and Langevin noise η with zero mean and η(t, s) · η(t , s ) =4k B T ζδ(t − t )δ(s − s ).
To numerically solve the polymer dynamics we employ a bead-spring representation of the polymers [13,14]. For most simulations, we adapted length scales close to those observed in Refs. [5,8]: κ −1 0 =0.5 µm, L=0.9 µm, p =10 µm. The relevant dimensionless parameters that characterize the system are the reduced noise σ and density ρ. Here, σ:=k B T p /(ζv 0 L 2 ) relates thermal forces at length scale p with friction forces, and ρ:=(R 0 /b) 2 denotes the squared ratio of the radius of curvature R 0 =κ −1 0 to the mean polymer distance b= A/M . For dilute systems, ρ 1, our simulations show that each polymer is propelled on a circular path and collisions between polymers are infrequent; see Fig. 2(a) and Movie 1 in the Supplemental Material [12]. The positions of the polymers' centers of curvature r (m) cc are uncorrelated as in a gas, and we refer to this state as a disordered state. On increasing ρ, we observe that a significant fraction of filaments begin to collide and collect into localized vortex structures (vortex state). These ring-like structures are highly dynamic. They assemble and persist for several rotations, during which their centers of mass remain relatively static; see Fig. 2(b) and Movie 2 [12]. Despite our simplified kinetic assumption, the overall phenomenology resembles the FtsZ patterns observed by Loose and Mitchison [5], including vortex assembly, disassembly and localization. In the dense regime, ρ 1, where each polymer is likely to collide, these vortices are unstable. Instead, the polymers cluster and form jammed 'trains' that travel through the system in an irregular fashion; see Fig. 2(c) and Movie 3 [12].
In order to quantitatively distinguish between the various observed patterns and organize them into a 'phase diagram' we consider the pair correlation function g(d cc ) [15,16] [12]. The ensuing 'phase diagram' is shown in Fig. 2(d). As in other active systems [17][18][19][20][21][22][23][24], pattern formation is favored by increasing density and decreasing noise strength. Jammed states prevail only when density is high and noise level low. Note also that the structure of the phase diagram depends on the ratio of filament length L to radius of curvature R 0 . Polymers with an arc angle close to κ 0 L=2π (closed circles) retain a singlecircle structure and do not form any collective structures upon increasing ρ [Movie 4 [12]]. Conversely, reducing κ 0 L suppresses the formation of closed ring structures, due to inefficient alignment of short polymers. Instead, these polymers cluster into flocks which move on approx- imately circular paths [Movie 5 [12]]. Hence, we conclude that the range of arc angles of FtsZ polymers, κ 0 L≈0.6π, observed in vitro [5], facilitates the formation of closed polymer rings particularly well [ Fig. 2(b)]. In summary, closed polymer rings require explicit curvature and filament lengths larger than a certain threshold value. For other interactions than local, steric repulsion ring structures may also emerge [1,25,26]; straight, rotating rods may form vortex arrays but not closed rings [27]. We complement the Brownian dynamics simulations of active particles that are propelled on circular tracks by considering the mesoscopic limit of vanishing particle extension. To this end, we have employed a kinetic Boltzmann approach [20,[28][29][30][31][32][33] to determine the collective behavior and the corresponding phase transitions in this limit, irrespective of the microscopic details of the constituent particles. In detail, we simplified the active system to one consisting of spherical particles (of diameter d) moving clockwise with constant speed v 0 on circular orbits of radius R 0 . This accounts for both selfpropulsion and spontaneous curvature but neglects the finite extension of the polymers as compared to our Brownian dynamics simulations.
We further assume that a particle's orientation is altered by 'self-diffusion' as well as by local binary collisions. In self-diffusion, a particle's instantaneous orientation θ changes at rate λ into θ+η, where we assume η to be Gaussian-distributed with standard deviation σ. As in other particle-based active systems [29,31,34], binary collisions are modeled by a polar alignment rule where the orientations of the collision partners align along their average angle plus a Gaussian-distributed fluctuation; for simplicity, we take the same width σ as for self-diffusion.
The kinetic Boltzmann equation [20,[28][29][30][31][32][33] for the one-particle distribution function f (r, θ, t) then reads It describes the dynamics of the density of particles in phase-space element drdθ which is being convected due to particle self-propulsion and which undergoes rotational diffusion and binary particle collisions, as given by the collision integrals I d [f ] and I c [f, f ], respectively; for explicit expressions please see the Supplemental Material [12]. Note here the critical difference to field theories for straight-moving particles [29,[35][36][37]; there is an additional angular derivative in the convection term, which reflects the fact that the particles are moving on circular orbits. In the following we rescale time, space and density such that v 0 =λ=d=1. Then, the only remaining free parameters are the noise amplitude σ, κ 0 , and the mean particle densityρ=A −1 A dr π −π dθ f (r, θ, t) measured in units of λ/(dv 0 ), i. e. the number of particles found within the area traversed by a particle between successive self-diffusion events.
To identify possible solutions of the Boltzmann equation and analyze their stability, we performed a spectral analysis. Upon expanding the one-particle distribution function in terms of Fourier modes of the angular variable, f k (r, t)= π −π dθ e iθk f (r, θ, t), one obtains where explicit expressions for the collision kernels I n,k (σ) are given in the Supplemental Material [12]. For k=0, Eq. (2), yields the continuity equation ∂ t ρ=−∇·j for the local density ρ(r, t):=f 0 (r, t) with the particle current given by j(r, t)=v 0 (Ref 1 , Imf 1 ) T . In general, Eq. (2) constitutes an infinite hierarchy of equations coupling lower with higher order Fourier modes. A linear stability analysis of Eq. (2) enables further progress. Since I n,0 =0 for all n, a state with spatially homogeneous densityρ=f 0 and all higher Fourier modes vanishing is a stationary solution to Eq. (2) (disordered state). To linear order, the dynamics of small perturbations δf k with respect to this uniform state is given . For a polar collision rule, as considered here, only µ 1 can become positive, defining a critical density ρ c (σ) at µ 1 (ρ c , σ):=0 [Fig 3(a)]. Above threshold (ρ>ρ c ), the spatially homogeneous state is unstable, the particle current grows exponentially, and collective motion may emerge.
In close proximity to the critical density ρ c (σ) a weakly non-linear analysis yields further insights into the dynamics of the system and the ensuing steady states. Here we follow Ref. [28] and assume small currents f 1 1 at onset. Then, balancing of the terms in the continuity equation, the equation for f 1 , and terms involving f 1 in the equation for f 2 implies the scaling ρ−ρ∼f 1 , f 2 ∼f 2 1 as well as weak spatial and temporal variations ∂ x/y ∼f 1 , ∂ t ∼f 1 . To include the lowest-order damping term in f 1 , we retain terms up to cubic order in f 1 . This yields the following hydrodynamic equation for the complex particle where ∇:=∂ x +i∂ y . While this equation shows the same functional dependencies on local density and current as found in systems with straight propulsion [29], the coefficients α, ξ, ν, γ and β are now complex-valued (for explicit expressions please see the Supplemental Material [12]). This can be traced back to the angular convection term in Eq. (1), or equivalently to the corresponding phase-shift term in Eq. (2). As a consequence, the field theory of active systems with particles moving on circular orbits with defined chirality is generically given by a complex Ginzburg-Landau (GL) equation with convective spatial coupling as well as densitycurrent coupling. This constitutes a highly interesting generalization of the standard (diffusive) complex GL equations [38,39], and is qualitatively different to real GL-type equations that were previously applied in the context of self-propelled particles [28]. Above threshold, ρ>ρ c (σ), the active chiral hydrodynamics described by the generalized GL equation Eq. (3) exhibits a uniform oscillatory solution with f 1 =F 1 e iΩ0t , i. e. a state in which particles move on a circular (chiral) path with an an- 1/2 gives the particle density. However, a linear stability analysis of Eq. (3) shows that for densities slightly larger than ρ c this oscillatory solution is linearly unstable against finite wavelength perturbations in the current and density fields. Preliminary numerical solutions of the generalized GL equation, Eq. (3), take the form of rotating spots of high density that appear to show turbulent dynamics [12,40]. This is qualitatively distinct from the high-density bands found for straight-moving particles [19,41] and the vortex field of a fluid coupled to torque dipoles [42,43].
Far above threshold, closure relations such as those discussed above [28] may become invalid and with them the ensuing hydrodynamic equations. Therefore, we proceed with the full spectral analysis of the Boltzmann equation, Eq. (2), as detailed in the Supplemental Material [12]. First, we numerically calculate the spatially homogeneous solutions for all angular Fourier modes f k below some cutoff wave vector k max . For given values of ρ and σ and a desired accuracy ε of this mode truncation scheme, the cutoff is chosen such that |f kmax+1 |<ε. We find that forρ<ρ c (σ) a spatially homogeneous state where all modes but f 0 vanish is the only stable state. In contrast, above threshold (ρ>ρ c (σ)) there is a second solution for which |f 1 |>0. It corresponds to a polar ordered state whose orientation is changing periodically in time with frequency v 0 κ 0 . For moderateρ−ρ c , the amplitude quantitatively agrees with the result from the generalized GL equation; see Supplemental Material [12]. In a second step, we consider wave-like perturbations, δf k (q) with wave-vector q, of the spatially homogeneous oscillatory solution in a co-rotating frame. The largest real part of all eigenvalues of the corresponding linearized system for δf k then yields the linear growth rate S(q) [ Fig. 3(b)]. In accordance with the linear stability analysis of Eq.(3), we find that for densities slightly larger than ρ c a spatially homogeneous solution is unstable against finite wavelength perturbations. The dispersion relation S(q) exhibits a band of unstable modes, with the maximal growth rate S max decreasing as one moves away from the threshold ρ c [ Fig. 3(a,b)]. Actually, there is lobe-like regime in parameter space where S(q)<0 [ Fig. 3(a)], and hence a homogeneously polar ordered state with rotating direction is stable. We emphasize here that our stability portrait [ Fig. 3(a)] is independent of κ 0 and hence equally valid for systems of straight-moving particles. For our two approaches ( Fig. 3(a) and Fig. 2(d)), the onset to order is governed by a similar trend [12], common for active systems [24,44]: disorder prevails for low density or high noise, order is promoted for high density or low noise.
To determine the spatiotemporal dynamics in the regime where neither a spatially homogeneous state nor a homogeneously polar ordered state are stable, we resort to a modified version of the SNAKE algorithm [31] to numerically solve Eq. (1). It accurately reproduces the threshold value ρ c (σ) at which the spatially homogeneous state becomes unstable [ Fig. 3(c)]. Above threshold (ρ>ρ c ) we find that local density fluctuations quickly grow and evolve into stable swirls, i. e. disc-like flocks of high density and polar order moving on circular paths; see Fig. 3(d), and Movie 6 in Supplemental Material [12]. The radius of such a path is approximately given by R 0 . These swirl patterns closely resemble the swirling flocks observed in the Brownian dynamics simulations for short polymer arc angles [Movie 5 [12]], as well as our preliminary numerical solutions of the generalized GL equation, Eq. (3) [12,40]. Moreover, in accordance with the spectral analysis, we find a second threshold density, above which the system settles into a homogeneously polar ordered state with a periodically changing orientation [Movie 7 [12]]. The amplitude and frequency of the polar order agree with the numerical results of the spectral analysis to high accuracy [12], while the numer- ically determined phase boundaries differ. The SNAKE algorithm produces stable swirl patterns only in a parameter regime where our linear stability analysis yields significant growth rates. This is mainly due to spurious noise caused by the discretization of the angular variable, which tends to suppress inhomogeneities in the regime of small growth rates. Furthermore, the finite system size constricts the band of possible modes and allows only for patterns of sufficiently short length scales.
For active systems of circling particles that interact via steric repulsion, our microscopic and mesoscopic treatments strongly suggests that a phase of collective vortex structures is a generic feature. Within this class, our work shows that extended polymers which as a whole follow circular tracks can form closed rings. Concerning our motivation of circling FtsZ, further research is needed to elucidate the dynamics of treadmilling; yet our minimal kinetic assumption suggests that varying particle density alone suffices to regulate the patterns as observed by Loose and Mitchison [5]. Compared to systems of straight moving particles we find qualitatively new phenomena [12,40]. For those systems, it was already reported that (globally achiral) vortices can occur due to collisions of particles of asymmetric shape [45] or due to memory in orientation [25,46]. Some of our findings, like the polymer length dependence of patterns and the possible emergence of active turbulence [47,48], pose interesting questions for future work. Our analysis yields a mapping of the emergent dynamics onto a generalized Ginzburg-Landau equation, providing a connection between active matter and nonlinear oscillators [40].
We thank F. Thüroff  In their experiments [S1], Loose and Mitchison observe that FtsZ polymers undergo depolymerization and polymerization processes leading to an effective translation in the direction of the polymers' backbones. However, the underlying molecular details are unclear, as they involve many qualitatively and quantitatively unknown reactions and a yet unstudied interplay of different auxiliary proteins (e.g. FtsA, ZipA). Here, we neglect these details and focus on the collective effects of many FtsZ polymers retaining only their effective movement along circular tracks. To realize this kind of motion we assume an intrinsic particle velocity.

NUMERICAL IMPLEMENTATION OF BROWNIAN DYNAMICS
In the following, we discuss the details of the implementation of the Brownian dynamics simulations. We use a bead-spring model [S2, S3] that comprises the following discretization scheme: a polymer of length L is subdivided into N beads at positions r i = (x i , y i ) T (i = 1, 2, ..., N ), with N − 1 bonds of length a; the (normalized) bond vectors are given by ∂ s r ≈ ri+1−ri a =:t i ; the bending angle between two adjacent bonds is given by θ i = arccos(t i+1 ·t i ). The corresponding bending energy reads where θ 0 ≈ aκ 0 is the spontaneous bending angle. In the bead-spring model, neighboring beads are connected by stiff harmonic springs. The corresponding stretching energy is given by In the simulations, the spring constant k is chosen larger than all other force constants to account for the fact that biopolymers are nearly inextensible; as a consequence, stretching modes relax fast compared to other dynamic processes. At the same time, k cannot be chosen arbitrarily large as this would strongly limit the maximal simulation time T max (see below for values).
In the two-dimensional system of M polymers, we assume steric repulsion between adjacent polymer segments r (m) i (m = 1, 2, ..., M ). As an interaction potential we use a truncated Lennard-Jones potential [S4-S6] with r j |, the potential strength, and Θ(r) the Heaviside step function. At distances smaller than the bond length a, the potential is strongly repulsive.
In the Langevin description, the equation of motion is given by a force balance between elastic, active, thermal and dissipative terms. For the i-th bead of a polymer, the equation of motion reads . We employ the following implementation of the tangential propulsive force F prop = ζv 0 ∂ s r: For the integration of Eq. (S4) we use an Euler-Maruyama iteration scheme [S7] with sufficiently small time steps ∆ = 0.0001τ with the unit time τ = ζa 2 /(k B T ). In our simulations, we used the following set of parameters: L = 9a, p = 100a, k = 500k B T /a 2 , = 1k B T, θ 0 = 0.2, ζ = 1 and a periodic system of area A = 60a×60a (such that it can contain many consecutive polymer lengths). In the main text, the unit of length is set to a = 100 nm, such that L = 0.9 µm, p = 10 µm are roughly similar to FtsZ filaments. The noise strength σ = k B T p /(ζv 0 L 2 ) was varied as follows: we changed the temperature scale in the interval k B T ∈ [0, 1] for v 0 = 5, and for k B T = 1 varied v 0 in the range v 0 ∈ [1,5]. The maximal simulation times T max for all simulations in the main text were chosen such that the single polymer rotation time τ R = 2π/(κ 0 v 0 ) is much smaller. We took T max > 400τ R and T max > 700τ for our data to provide a sufficiently large sampling interval for both convective and diffusive motion. To consolidate the results, data were recorded for 10 independent simulation for each given set of parameters.

ANALYSIS OF THE PAIR CORRELATION FUNCTION
To analyze the patterns observed in the Brownian dynamics simulations, we consider the pair correlation function g(d cc ) [S8, S9] of center distances d cc = |r cc are the curvature centers of each polymer, generated by averaging over the local curvature and all local reference positions on a contour (see Fig. S1(a)). In contrast to the positions r (m) , the curvature centers do not oscillate due to self-propulsion and hence represent a more stable measure of particle position.

DERIVATION OF THE HYDRODYNAMIC EQUATIONS
To assess the dynamics at larger scales, we employed a kinetic Boltzmann approach. The corresponding generalized Boltzmann equation for f (θ, r, t) is given by Eq. (1). The self-diffusion and collision integrals I d and I c , respectively, are given by where S(ψ) = 4dv 0 | sin( ψ 2 )| is the scattering cross section for spherical particles of diameter d and velocity v 0 in two dimensions as detailed in Ref. [S10]. The collision integral represents ferromagnetic alignment of two particles with orientation φ 1 and φ 2 along their average angle θ = 1 2 (φ 1 + φ 2 ). The brackets denote an average over a Gaussian-distributed noise variable η. To obtain a dimensionless form we used the rescaling with ρ 0 = λ/(dv 0 ). Measuring time, space and density in units of λ −1 , v 0 λ −1 , and ρ 0 , respectively, allows to set d = λ = v 0 = 1. Then, the only remaining free parameters are the noise amplitude σ, κ 0 , and the mean particle densityρ = A −1 A dr π −π dθ f (r, θ, t). To proceed, we performed a Fourier transformation of the angular variable: f k (r, t) = π −π dθ e iθk f (r, θ, t). This leads to the Boltzmann equation in Fourier space, Eq. (2), where the Fourier transforms I n,k are given by S(|Φ|) P k cos(Φ(n − k/2)) − cos(Φn) .
The equation for f 1 then couples to the nematic order field f 2 via a term ∼ f * 1 f 2 of order f 3 1 , where the star denotes complex conjugate. Writing down contributions from Eq. (2) for k = 2 of quadratic order in f 1 yields an expression for f 2 as a function of f 1 . The expression for f 2 can then be substituted into Eq. (2) for k = 1 to obtain a closed equation for f 1 . Together with the continuity equation, the hydrodynamic equations for the density and the particle current read where ∇ := ∂ x + i∂ y . The coefficients are given by α := (I 0,1 + I 1,1 ), We note that the employed truncation scheme implies fast relaxation of the nematic order field f 2 such that ∂ t f 2 is assumed to be negligible on time scales of the dynamics of f 1 . f 2 is then slaved to f 1 via f 2 = −2ν∇f 1 + γf 2 1 .

Homogeneous isotropic state
To study the stability of the homogeneous isotropic state we substitute ρ =ρ + δρ and f 1 = δf 1 with the wave-like perturbations of the form δρ(r, t) ∼ δρ q e iq·r , where δρ q and δf 1,q are in general complex amplitudes that are assumed to be small. Periodic boundary conditions in our numeric solution impose |q| = n 2π L , n Z, where L = √ A and A is the area of the (quadratic) system. The linearized set of equations of motion for the perturbations δρ q (t), δf 1,q (t) and δf * 1,q (t) has the characteristic polynomial where S is the eigenvalue of the linearized set of equations for δρ q (t), δf 1,q (t) and δf * 1,q (t). We note that [ν] is positive for all densities. Forρ < ρ c , all coefficients in (S12), including the S-independent terms are positive, such that (S12) only yields S with negative real part. Thus, forρ < ρ c the homogeneous isotropic state is linearly stable against inhomogeneous wave-like perturbations. Forρ − ρ c > 0, the real part of S becomes positive where the fastest growing mode is always at q = 0.

Homogeneous oscillatory state
To study the stability of the homogeneous oscillatory solution we substitute small perturbations in the basis of the homogeneous oscillating solution: where the amplitudes δρ (0) , δρ (1) , δf (0) , δf (1) and δf (2) are again of the form (S11). Truncating at the lowest order of (ρ − ρ c ), which is α(ρ − ρ c ), yields a closed set of linear equations for the amplitudes. The eigenvalue with the largest real part of this linear system determines the growth rate S(q) of wave-like perturbations. We find that the dispersion relation yields positive S(q) for finite q (see Fig S2).

NUMERICAL LINEAR STABILITY ANALYSIS IN THE FULL PHASE SPACE
In the derivation and the stability analysis of Eqs. (S9) we rely on the assumption of small particle currents which might be justified at onset. However, this assumptions is in general questionable and not well justified for densities much larger than ρ c . To obtain a stability map for the full phase space (Fig. 3), we first calculated the homogeneous solution of Eq. (2) retaining only modes up to k max . Given some values ofρ and σ and a desired accuracy of this mode truncation scheme the cutoff is chosen such that |f kmax+1 | < . As a next step, we linearized Eq. (2) with respect to this solution and calculated the maximal growth rate S(q) of wave-like perturbations with wave vector q. If S(q) > 0 for some |q|, the homogeneous solution is unstable whereas if S(q) < 0 for all |q|, the corresponding homogeneous solution is stable.
Note that the homogeneous version of Eq. (2) (neglecting all gradient terms) is invariant under a phase shift f k → f k e ikv0κ0t . Choosing the orientation of the polar order at t = 0 to be aligned along the x-axis, Eq. (2) is solved by f k = |f k |e ikv0κ0t with the time and space independent amplitude |f k |. |f k | is then determined by the stationary homogeneous version of Eq. (2): This equation is identical to the stationary homogeneous Boltzmann equation for straight moving particles; i.e. where κ 0 = 0. Hence, the solutions for the amplitudes |f k | are identical to the solutions for the Fourier modes in systems of straight moving particles [S10]. To proceed, we truncate the infinite sum in Eq. (S14) at k max and calculate the solution of all |f k | with |k| ≤ k max . Fig. S3 depicts the solution for the amplitude |f 1 | as compared to the solution of the generalized Ginzburg-Landau equation as well as the SNAKE algorithm. The explicit solution for |f 1 | and higher modes justifies the scaling scheme used to derive Eqs. (S9) in the vicinity of ρ c [Fig. S3, inset]. For decreasing noise σ or increasing densityρ an increasing number of Fourier modes starts to grow [Fig. S3, inset]. In our numerical calculations we typically included 30 − 50 Fourier modes. The dashed region in Fig. 3(a) indicates the regime where we cannot find a nontrivial solution to Eq. (S14) by neglecting Fourier modes above the chosen k max = 50 and where we would have to choose a larger k max . With the substitution f k = (|f k | + δf k )e ikv0κ0t the linear system for δf k then reads Here, we performed a coordinate transformation to a frame rotating with angular frequency κ 0 such that ∇ → e ikv0κ0 ∇. Assuming wave-like perturbations as in Eq. (S11), we solved Eq. (S15) for the maximal eigenvalue and get the growth rate as a function of the wavenumber in the rotating frame (see Fig. 3(b)). The maximum taken over all wavenumbers |q| > 0 then defines the maximal growth rate S max of wave-like perturbations.
In agreement to previous results [S10], we found that the growth rate is maximal for q parallel to the particle current. The contour plot of S max as a function of ρ and σ yields the phase diagram Fig. 3(a). Note again, that our stability analysis and the resulting phase diagram Fig. 3(a) is independent of curvature and also valid for the well-studied system of propelled particles without curvature [S10, S12, S13]. Hence, Fig. 3(a) shows that the Boltzmann approach is capable of reproducing phases of all states observed in [S12, S14] including a transition from travelling wave patterns to global homogeneous order.

NUMERICAL SOLUTION OF THE BOLTZMANN EQUATION WITH SNAKE
In order to study the resulting steady states in the regime where our linear stability analysis predicts inhomogeneities, we numerically solved the generalized Boltzmann equation, Eq. (1). To this end we employed the SNAKE algorithm as introduced in Ref. [S15]. As tesselations we used a quadratic periodic regular lattice with equally sized angular slices. Circling propulsion was included by rotating the angular distribution of each lattice site with a frequency v 0 κ 0 in addition to the straight convection steps. The system was initialized with a disordered state with small random density fluctuations around the mean densityρ = A −1 A ρ(r, t). Changing κ 0 did not change the observed patterns qualitatively. In the limiting case of very small κ, we observed traveling wave patterns as reported in Refs. [S12, S14, S15]. For Fig. 3(c), Movie 6, and Movie 7 we used a lattice of of 200 × 200 grid points with lattice field size 2 and angular disretization of 24 angular slices; hence, A = 400 × 400 = 160000. In the swirl phase the swirl size grows for growingρ − ρ c whereas the radius of a swirl's motion stays at approximately κ −1 0 . Fig. S4 shows the parameter values ofρ and σ where the SNAKE algorithm exhibits steady swirl patterns together with the phase diagram obtained from the adapted mode truncation scheme.

REMARK ON THE SHAPE OF THE PHASE CURVES
When comparing the transition to order in the phase diagrams 2 and S4 it should be noted that our particlebased and continuum approaches are distinct in the following features: polymer fluctuations vs. effective diffusion, multi-particle collisions vs. binary alignment, extended polymers vs. point particles. The functional form of ρ c (σ) (S10) depends on the choice of diffusion and collision noise (e.g. equally Gaussian distributed). In contrast, the form of the transition line in our Brownian dynamics simulations depends on the choice of the phenomenological criteria (disordered states for d min cc ≈ 2R 0 , vortex states for d min cc > 2R 0 and train states without d min cc ). These differences result in different shapes of the phase boundaries. In addition, the observed patterns in the vortex phase are distinct. While for our particlebased model we find closed, rotating rings, dense, rotat- ing swirls are observed in the continuum model ( Fig. 2(b) and Fig. 3(d)). These differences are interesting and should be considered as part of the results we obtained. For example, these differences will guide future model building for specific models, e.g. the dynamics of FtsZ, as they emphasise what molecular details need to be accounted for. For the discussion of this work, however, our emphasis was on the topology of the phase diagram (similar trend of the onset to order) and the fact that in both models one finds a vortex phase.
Movie5.mp4: Brownian dynamics simulation with parameters as in Movie 3, except for a changed contour length L = 6, resulting in an polymer arc angle Lκ 0 = 1.2.
Movie6.mp4: SNAKE solution forρ = 0.2 and σ = 0.45 with κ 0 = 0.1. The colour code denotes the local density ρ/ρ. The orientation and length of the arrows indicates the orientation and amplitude of the local particle current. Movie7.mp4: SNAKE solution forρ = 0.75 and σ = 0.2 with κ 0 = 0.1. The colour code denotes the local density ρ/ρ. The orientation and length of the arrows indicates the orientation and amplitude of the local particle current.