Two-species active transport along cylindrical biofilaments is limited by emergent topological hindrance

Active motion of molecules along filamentous structures is a crucial feature of cell biology and is often modeled with the paradigmatic asymmetric simple exclusion process. Motivated by recent experimental studies that have addressed the stepping behavior of kinesins on microtubules, we investigate a lattice gas model for simultaneous transport of two species of active particles on a cylinder. The species are distinguished by their different gaits: While the first species moves straight ahead, the second follows a helical path. We show that the collective properties of such systems critically differ from those of one-species transport in a way that cannot be accounted for by standard models. This is most evident in a jamming transition far below full occupation, as well as in non-equilibrium pattern formation. The altered behavior arises because - unlike the case in single-species transport - any given position may be targeted by two particles from different directions at the same time. However, a particle can leave a given position only in one direction. This simple change in connectivity significantly amplifies the impact of steric interactions and thus becomes a key determinant of mixed species transport. We computationally characterize this type of hindrance and develop a comprehensive theory for collective two-species transport along a cylinder. Our observations show high robustness against model extensions that account for additional biomolecular features and demonstrate that even small fractions of a second species can significantly alter transport. This suggests that our analysis is also relevant in a biological context.


I. INTRODUCTION AND MOTIVATION
Efficient collective molecular transport is a vital prerequisite for a multitude of processes in cell biology on many different levels. Examples range from messenger RNA (mRNA) translation to organelle transport. Typically, highly functional molecular motors that transform chemical energy into stepwise mechanical translation perform this complex task by moving along filamentous structures such as the cytoskeleton or mRNA [1,2]. Of particular importance for intracellular organization are the cylindrically shaped and polarized microtubules. Kinesins-the molecular motors associated with microtubules-exhibit distinct efficient motility, as they are capable of "walking" processively over micrometer distances towards the microtubule end [1].
Many experimental studies have focused on elucidating the microscopic working mechanisms of molecular motors [3][4][5], but understanding their collective behavior [6] remains a challenging task. For this reason, concepts drawn from statistical physics and modeling have proven to be of much relevance, as they offer a means of linking the microscopic to macroscopic behavior. In this context, the totally asymmetric exclusion process (TASEP) [7,8] and extensions thereof have proven particularly fruitful. In its simplest formulation, the TASEP accounts for two central aspects of transport: active motion and steric interactions. It describes the motion of point particles along a one-dimensional lattice and therefore serves as an ideal basis to study collective transport of one-particle species along a single track. Despite its simplicity, this model shows a surprisingly rich phenomenology and captures many essential features of transport processes. Indeed, by now it has acquired the status of a paradigmatic model, not only for transport [8,9] but also for nonequilibrium physics, in general, comparable to that of the Ising model for equilibrium physics.
Motion of kinesins was initially thought to occur mostly along the so-called protofilaments-separate lanes oriented parallel to the axis of the cylindrical microtubule (see Fig. 1, green kinesin). While this motion along a single protofilament is a well-studied feature of one of the most prominent kinesins involved in intracellular transport, kinesin-1 [10], several in vitro studies have revealed notable exceptions to this behavior. Members of all three superfamilies of molecular motors-kinesin [11][12][13][14], dynein [15], and myosin [16]-have been shown to produce torsional force during their axial translation. While the precise basis for this observation has remained unclear, recent studies of kinesin-2 [17] and the microtubule depolymerizing kinesin-8 [18,19] strongly suggest that these molecular motors regularly switch protofilaments and show a bias towards one side. This effectively results in a helical motion along the microtubule [17,20,21] (see Fig. 1, orange kinesin). On a theoretical level, systems with a single species of kinesins that stochastically switch protofilaments have been studied recently by employing extended TASEP-like models with multiple lanes [22]. In that case, the qualitative behavior of collective transport along the cylindrical microtubule was reported to be widely conserved as compared to models for kinesins that track protofilaments. However, many different types of molecular motors are present in a single cell and, as lane switching is likely to be the rule rather than the exception, the general scenario is that several different molecular motor species should interact on a single cytoskeletal filament. This raises the question of how the interplay of distinct molecular motor species that show different gaits alters collective transport along a cylindrical structure.
Inspired by molecular transport on microtubules, we employ lattice gas models to find generic principles of collective transport by two species of particles that are distinguished by different gaits on a cylinder. Here, we demonstrate that the emerging behavior of such systems critically differs from collective behavior in the presence of a single species only. The simultaneous presence of molecular motors that follow a straight and a helical course, respectively, inevitably leads to crossings between their trajectories: A certain lattice site may be targeted from two different directions but, once occupied, can only be vacated in a single direction. This modification of connectivity in the network topology of potential particle movements amplifies the impact of steric interactions globally and thus hinders particle motion-an effect we call topological hindrance. Topological hindrance produces highly nontrivial correlations between the dynamics of particles of the different species. Specifically, the particle current and distribution on the filament are now dependent not only on the total number of particles but also on the fraction of the respective species. The impact of topological hindrance is most evident in the jamming of particle flows at densities far below full occupation. We present an analytical framework that quantifies topological hindrance and provides a theoretical basis for understanding twospecies transport along cylinders, much as the TASEP does for transport by a single motor species. Moreover, our model predicts nonequilibrium patterns in the particle distribution that have not been observed in classical models for single-species transport. To specifically target the robustness and biological relevance of topological hindrance, we further investigate extended models that account for specific biomolecular features. We find that topological hindrance still plays a key role for collective transport properties in these cases. While the extended models are too complicated for an analytic investigation, we can understand their behavior on the basis of our idealized model and our theory for topological hindrance. This paper is organized as follows: We begin with a review of the collective transport by a single species in Sec. II. Our model for collective transport by two species on a cylindrical structure is presented in Sec. III. The phenomenology of our model and the key differences to transport with a single species are discussed in Sec. IV, which also provides a qualitative explanation of topological hindrance. Furthermore, we address how the system dimensions influence topological hindrance. We then develop a theory to quantify the strength of topological FIG. 1. Active transport by two species along a cylindrical structure. (a) Kinesins (e.g., kinesin-1, green) may track a single lane of a microtubule (protofilament), but many molecular motor species have also been reported to regularly switch protofilaments in a biased fashion and thereby effectively undergo a helical motion (e.g., kinesin-2, orange). (b) Model implementation to study collective two-species transport along cylindrical structures. We consider a lattice with length L and width W. Periodic boundary conditions are employed along the transversal direction to account for a cylindrical geometry (e.g., microtubule). Particles of species T (lane tracking, green) hop to the neighboring lattice site on the right at rate ν T , and particles of species S (lane switching, orange) hop to the neighboring lattice site on the upper-right at rate ν S . All particles exclude each other. Particles enter the system with rates α T and α S at the left boundary and exit the system at the right boundary with rates β T and β S , respectively. Gray areas denote system boundaries.
hindrance for arbitrary particle densities in Sec. V, which yields the current-density relation for an arbitrary number of lanes and species fractions. This allows us to compute the complete phase diagram of our model, i.e., the particle density and current emerging in the system as a function of the control parameters. In Sec. VI, we discuss the robustness of our results against model modifications and their biological relevance, and provide a guideline for potential experimental verification. Finally, in Sec. VII, we relate our work to existing mathematical theories of driven systems and discuss its applicability and possible relevance in the biological context.

II. REVISITING ONE-SPECIES TRANSPORT
We start with a summary of the TASEP, which is one of the most fundamental models used to describe active transport of sterically interacting agents along defined pathways. While exact results relating to its properties have had a major impact on the field of nonequilibrium physics, in general [23][24][25][26][27], its applications cover a vast variety of transport processes, such as ion channels [28,29], spin transport [30,31], traffic flow [32], mRNA translation [33], and intracellular transport [8]. Here, the TASEP will serve as the starting point to treat molecular transport along cytoskeletal structures in the presence of a single motor species.
The model is defined as follows: Point particles populate a one-dimensional lattice along which they can hop stochastically at a rate ν to the neighboring lattice site on their right. To account for steric interactions, the particles exclude each other, such that a lattice site can only be occupied by a single particle. The particles enter the lattice from the left at a rate α and leave the lattice on the right at a rate β. This can be interpreted as connecting the lattice with two particle reservoirs with fixed densities ρ L ¼ α on the left and ρ R ¼ 1 − β on the right.
With these definitions, the TASEP allows one to study macroscopic properties that emerge in transport processes. Two central observables are the average particle density ρðx; tÞ (average particle distribution) and the average particle flux Jðx; tÞ, where x denotes the position in the system and t the time. Particles cannot be created or annihilated within the lattice. Therefore, a spatial difference in the particle flux must lead to a temporal change in the particle distribution in the ensemble average. This is reflected by a continuity equation that describes the system's temporal evolution on a macroscopic scale [34], Both the TASEP and the model considered in this paper are ergodic, and the Perron-Fobenius theorem holds true [35]. This means that they will evolve into a unique nonequilibrium steady state, on which we focus from now on.
For the TASEP and many other transport models, there is a unique connection between the current and the density, the current-density relation JðρÞ. The existence of this unique function means that the local current J is completely determined by the density ρ; J depends on the in rate α and out rate β only implicitly via the density. This phenomenological approach based on a unique current-density relation goes back to the work of Lighthill and Whitham [36] and is at the heart of many theories for various transport models [22,[37][38][39]. For the TASEP, the currentdensity relation is given by [23,24] This equation reflects the fact that particles may only move to empty lattice sites: Neglecting correlations, the probability of finding a particle (given by ρ) must be multiplied by the probability that a lattice site is empty (given by 1 − ρ) to obtain the current. Although this is just a heuristic mode of argumentation, it can be proven that the corresponding Eq. (2) is an exact relation for the TASEP on an infinite-dimensional lattice [23][24][25]. This relation also implies that the particle current vanishes if the density is either zero (no particles) or one (full occupation). In the latter case, the motion of every particle except for the last one is blocked, and consequently, no hopping within the lattice can occur. We refer to such a system as jammed.
Based on the current-density relation, it is possible to describe the complete macroscopic behavior in terms of the control parameters of the TASEP in the stationary state. Specifically, it uniquely determines the average density and particle current of the stationary state associated with a specific choice of control parameters α and β. A very successful theoretical concept in this context is the extremal current principle [37,38,40]. It states that, given a set of possible densities, the system will always realize the one corresponding to the extremal current. As explained above, the lattice boundaries are effectively connected to particle reservoirs of density ρ L ¼ α and ρ R ¼ 1 − β on the left and right ends, respectively. According to the extremal current principle, the steady-state density and current are then given by The current-density relation of the TASEP, Eq. (2), exhibits only a single maximum. Consequently, the extremal current principle predicts boundary-induced transitions between three different phases: Eq. (3) implies that either the left reservoir density ρ L , the right reservoir density ρ R , or the density corresponding to the sole maximum in the currentdensity relation, ρ MC ¼ 0.5, is a valid solution and therefore realized in the system. On an intuitive level, when the left TWO-SPECIES ACTIVE TRANSPORT ALONG … PHYS. REV. X 8, 031063 (2018) 031063-3 reservoir density is realized, a lack of particles determines the system's behavior, as the in rate is too small to create large particle jams. The corresponding phase is called the low-density (LD) phase. When the out rate is so small that it is the factor that limits transport, particles start to jam at the lattice end. Ultimately, the corresponding outflux will determine the overall particle current. The system then realizes the right reservoir density. This phase is called the high-density (HD) phase. Third, if neither a lack of particles nor the out rate is limiting transport, the particles' motion itself and thus the maximal possible current and the corresponding density constitute a constraint. This phase is the maximal-current (MC) phase. It is worth noting that the above derivation of phase transitions is independent of the microscopic rules of a system and has also been applied successfully to drivenlattice-gas models other than the TASEP [40,41]. Therefore, the extremal current principle suggests that the phase behavior of a driven lattice gas is qualitatively identical to that of the TASEP as long as the curve characteristics of the current-density relation are conserved: a region where the current shows a monotonic increase with increasing density, a region where the current monotonically decreases with increasing density, and a single maximum.
The TASEP has been generalized in various ways. Several studies investigated systems with a more involved geometry such as two lanes [30,31,[42][43][44], junctions [45][46][47][48][49][50][51], or networks [52][53][54][55][56]. In particular, Curatolo et al. [22] treated TASEPs with an arbitrary number of parallel lanes on a cylinder, where particles are also allowed to switch lanes. The authors find that, in this case, the system reduces in many ways to the single-lane TASEP, which justifies its application for the case of a single species of lane-switching molecular motors on a microtubule. Besides these studies addressing more complicated geometries, several authors have treated multiple species of particles, typically in opposite directions [39,[57][58][59][60][61][62]. Yet, how a mixture of different species of molecular motors-that naturally may move with different gaits-behave collectively on the cylindrical microtubule structure remains elusive.

III. MODELING TWO-SPECIES TRANSPORT
As discussed in the previous section, the TASEP adequately describes molecular transport along a cylinder in the presence of a single species of molecular motors. We now turn to a mixture of two species of molecular motors that are distinguished by different gaits. Specifically, we address the question of how to describe collective transport in the presence of molecular motors that move parallel to the cylinder axis and molecular motors that follow a helical path as suggested by experiments [17][18][19]. To this end, we study the stochastic model with Markovian dynamics illustrated in Fig. 1. We consider a two-dimensional lattice composed of W parallel lanes, each with a length of L lattice sites. The system is populated by two different particle species: a lane-tracking species (T species) and a lane-switching species (S species). Particles of the lanetracking species stochastically hop at rate ν T to the neighboring lattice site on the right while staying on the same lane. In detail, using Latin letters to denote the lane index and Greek letters to denote the site index, hopping of T particles is described by i → i and μ → μ þ 1. They represent, for example, molecular motors that track a single protofilament. Particles of the S species change lanes with every hopping event and stochastically move at rate ν S to the neighboring lattice site on the upper-right, i.e., i → i þ 1 and μ → μ þ 1. To implement a cylindrical structure, periodic boundary conditions in the transversal direction are imposed. Thus, members of the S species hop from the uppermost to the lowermost lane. In this way, the laneswitching particles represent molecular motors that move in spirals around the microtubule. Furthermore, particles are subject to steric interactions. In the model, they exclude each other, and hopping events can only occur if the corresponding lattice site is empty. At the left boundary, an empty site is filled with a particle of the respective species at rates α T and α S . Conversely, at the right boundary, particles of the T and S species leave the lattice at rates β T and β S .
Note that the assumptions of lane switching in each step and the absence of random particle attachment and detachment (Langmuir kinetics) are simplifications from a biological point of view. Their aim is to isolate the basic principles of two-species transport along a cylinder, which, in turn, allows us to develop an analytic theory for topological hindrance. To bridge back to biological systems and prove the relevance of our concepts, we discuss extended models in Sec. VI.
To describe the state of the system, we use occupation numbers n T i;μ , n S i;μ ∈ f0; 1g for lattice site μ on lane i. Here, 1 indicates that the lattice site is occupied by a T or S particle, whereas 0 stands for the absence of the respective species. We focus our analysis on the average particle distribution (density) ρ and the average particle current J that emerge in the system with the above stochastic rules. More specifically, we define ρ X i;μ ≔ hn X i;μ i as the ensemble averaged occupation of site μ on lane i by a particle of type X ∈ fT; Sg. The current of S particles at site μ on lane i is defined as the average number of S particles hopping onto this site per unit time: J S i;μ ≔ ν S hn S i−1;μ−1 ð1 − n T i;μ − n S i;μ Þi. Equivalently, the current for the T species J T is defined as Because of periodic boundary conditions along the transversal direction, we implicitly use the identification i ¼ W þ 1 ≡ 1 and i ¼ 0 ≡ W in these relations and in the following. For later convenience, we also define the occupation number irrespective of the particle species, n i;μ ≔ n T i;μ þ n S i;μ . The temporal evolution of average occupations in the bulk of the system (μ ≠ 1; L) is then given by the master equations At the boundary sites μ ¼ 1, L, the equations read By rescaling time, it is possible to set one hopping rate equal to ν X ¼ 1 without loss of generality. For simplicity, we focus on identical hopping rates ν S ¼ ν T ¼ 1 throughout this paper.

IV. PHENOMENOLOGY OF TOPOLOGICAL HINDRANCE
A. Two-species transport cannot be reduced to one-dimensional or single-species transport In this work, we focus on steady-state properties, i.e., dhn X i;μ i=dt ¼ 0 and likewise for other moments. Equations (4) show that the average occupation of a lattice site depends on higher moments, which ultimately leads to an unclosed hierarchy of equations. This typically precludes an exact solution. The simplest analytic approach to treat Eqs. (4) and (5) is the mean-field approximation [8], where correlations are neglected by factorizing (second) moments, hn X i;μ n Y j;ν i ¼ hn X i;μ ihn Y j;ν i, which closes the hierarchy. The mean-field approximation works successfully for the TASEP and various similar models, and it therefore has acquired the status of a standard method for the analytical treatment of driven lattice gases [8]. We employ this factorization scheme for the total current J i;μ ≔ J S i;μ þ J T i;μ . Furthermore, our system is irreducible [63]. For a continuoustime Markov process, this suffices to show that the stationary state is unique [64]. Since there is only one stationary state, it has to adapt the symmetry of the system. Thus, all macroscopic quantities have to be independent of the lane number, and we therefore omit the lane index i in the following. This reasoning is further validated by stochastic simulations as shown in Appendix D 1. The mean-field current-density relation then reads Here, we also define the total particle density ρ μ ≔ ρ T μ þ ρ S μ at site μ. Equation (6) is identical to the current-density relation of the TASEP, Eq. (2). In particular, it predicts that the current is independent of the fraction of spiraling molecular motors and depends only on the total density. To relate the current and density to the system's control parameters, it is useful to introduce new quantities. First, the total in rate is given by α ≔ α T þ α S . Second, since particles enter the system independently, the fraction of the current contributed by lane-switching particles is δ ≔ α S =ðα T þ α S Þ.
Because of current conservation, this current fraction is spatially constant. In the mean-field analysis, J S =J ¼ ρ S =ρ holds true, such that the current fraction δ of the S species also equals the density fraction of S particles. Therefore, we can consider ρ, J instead of ρ S , ρ T , J S , and J T , and use δ to compute the respective fraction in the mean-field analysis. Third, we can identify an effective out rate β of particles irrespective of their species. The average dwell time of a particle at the last site is With these definitions of α and β, we then obtain the full phase-diagram predicted by a naïve mean-field approximation, which recovers all TASEP relations.
To test this mean-field analysis, we perform stochastic simulations based on Gillespie's algorithm [65] for a system with two lanes. Irrespective of the initial conditions, the dynamics converge to a unique stationary state that is characterized by a particle density ρ averaged over the whole system and a particle flux J that we numerically computed for various values of α, β, and δ. The result for a system composed of two lanes is shown in Figs. 2(a) and 2(b). In clear contradiction to the mean-field analysis, Eq. (6), we observe a strong dependence of the average current and density of particles on the fraction of spiraling molecular motors in the system. These findings falsify the mean-field approximation and show that the current is not uniquely determined by the density but carries an explicit dependence on the species fraction δ. For fixed δ, however, a unique current-density relation J(ρðα; β; δÞ; δ) can still be found, as shown in Fig. 2(a). Unlike the species fraction δ, the in rate α and out rate β of particles are parameters that change the density only locally and therefore the current only implicitly. The resulting current-density relation J(ρðα; β; δÞ; δ), as shown for W ¼ 2 lanes in Fig. 2(a), exhibits a symmetry upon exchanging the species. This is explained by an invariance of particle dynamics when lattice sites are relabeled [66]. For δ ¼ 0 and δ ¼ 1, i.e., in the presence of a single species only, we recover the TASEP current-density relation. If both particle species are present, the relation changes and the current is always less than in a single-species setup. Most interestingly, in a mixed system, the current vanishes already at densities below full occupation [ Fig. 2(a), white area]; densities above a maximal value of ρ max ðδÞ (bold gray line) are not realized in the steady state. Moreover, a very small fraction of lane-switching particles already cause this effect: Approximately 2%-5% of a second species (see Fig. 2) is sufficient to cause significant deviations from the single-species and mean-field results. Besides this critical difference in the maximum density of a mixed and a one-species system, the qualitative shape of the current-density relation remains unchanged when the species fraction is varied. As shown in Fig. 2 for the special case of W ¼ 2 lanes, we observe the same curve characteristics as for the TASEP. These are a single maximum corresponding to the MC phase, a LD region where the current monotonically increases with increasing density, and a HD region where the current monotonically decreases with increasing density. The extremal current principle, Eq. (3), then suggests that the phase diagram of our model and the TASEP are topologically equivalent, although the positions of the phase boundaries might differ.
As the mean-field approximation fails to describe twospecies transport, the species current fraction δ cannot be equated with the species density fraction. This complicates any further discussion as ρ S and ρ T have to be treated separately. However, as shown in Appendix D 2, in the context of this paper, no significant deviations of the density species fraction from the current species fraction occur. Therefore, for the remainder of this work, we mostly use ρ and implicitly assume a conversion to ρ S and ρ T via δ. Furthermore, we refer to δ simply as a species fraction.
In summary, our simulation data show that, first, the dynamics of the two species are correlated, as manifested by the failure of the mean-field approach and the dependence of the current and density on the species fraction δ. Second, correlations of the species dynamics always lead to a reduction of the average particle current; hence, the presence of a second species always hinders transport. We conclude that our model cannot be reduced to the TASEP and that naïve mean-field approaches are incapable of capturing the emerging macroscopic transport behavior.

B. Topological hindrance
To obtain insight into the interactions that reduce the particle flux in the presence of both particle species, we take a qualitative look at particle configurations that might arise in our model. First, we consider jammed configurations where all particles (except for the very last one) are blocked and where, consequently, the average particle current vanishes. As discussed in Sec. II, for one species transport, a jammed system can only arise trivially if every lattice site is occupied. In contrast, our model exhibits jamming at densities below full occupation. A possible jammed configuration of a mixed system is shown in Fig. 3(a). Although the system is jammed, there exist lattice sites that are not occupied. These sites are inaccessible to any particle of this configuration [ Fig. 3(a), crosses] and ultimately lead to jamming below the full occupation. In our model, two particles from different species can simultaneously target a single site. Consequently, if that site is occupied, a single particle can block the motion of two other particles [ Fig. 3(b), configuration 1]. This, in turn, can create sites that are empty but nevertheless inaccessible. Moreover, if particles of both species try to access a single site that is empty and accessible [ Fig. 3(b), configuration 2], their motion is restricted insofar as either one of them, but not both, 3. Illustration of topological hindrance for a system with W ¼ 2 lanes at the right lattice end. (a) State where the motion of each particle in the bulk is blocked (jammed state) and only the rightmost particle is allowed to move (arrow). Empty lattice sites (crosses) are present, which lead to average densities below the full occupation for jammed systems. (b) The second species, which changes the connectivity of the network of possible particle motions, amplifying the impact of steric hindrance in our model. Bracketed area 1 shows how a single particle can now block up to two other particles. Area 2 shows how particles might align such that motion has to occur sequentially. The same symbols as in Fig. 1 are used to illustrate the system and particles. can hop at a time. These configurations hence act as an intrinsic bottleneck and will further reduce the particle current on average. To assess these considerations also quantitatively, we specifically compute the number of inaccessible lattice sites and, with it, the maximal density ρ max as a proof of principle for a system with δ ¼ 1=2 and W ¼ 2 lanes. The detailed computations are shown in Appendix A and lead to a value of ρ theory max ¼ 0.74. This is in very good agreement with simulations that yield ρ sim max ¼ 0.73, which validates the qualitative arguments given above.
In summary, when the lane-switching species is added to the system, the network topology of possible particle motions is changed. Any given site is now potentially connected to two other sites, although leaving a site is only possible in one direction at a time. This creates a two-toone connectivity. The single-species model, on the other hand, shows one-to-one connectivity or, for the general case of particles that may stochastically switch lanes, n-ton connectivity. The two-to-one connectivity of our model amplifies the impact of steric interactions relative to singlespecies transport and therefore hinders motion. First, steric interactions may now occur at more points in space as compared to single-species transport, and thus give rise to inaccessible lattice sites [see Fig. 3(b), group 1]. Second, at intrinsic bottlenecks, steric interactions can now act at more points in time, as particles may have to hop one after the other [see Fig. 3(b), group 2]. Throughout this paper, we refer to these phenomena as topological hindrance.

C. Influence of the number of lanes
The first important question to address is how the density of a jammed system ρ max is influenced by the number of lanes W. Our analytic approach for W ¼ 2 lanes, however, is not feasible for large numbers of lanes as the complexity of the underlying mathematical problem increases rapidly. Figure 4 shows numeric results for the maximal system density obtained from stochastic simulations of systems with up to W ¼ 50 lanes. Strikingly, the maximal density ρ max decreases rapidly with an increasing number of lanes. Therefore, the overall impact of topological hindrance grows accordingly and becomes the major determinant of the system's dynamics. This is of importance as large numbers of lanes are often encountered in biological contexts. For example, microtubules are typically composed of 13 lanes, where our model already exhibits a maximal density of ρ max ≈ 0.2. Let us also emphasize that, in analogy to the two-lane system, current reduction and, likewise, jamming due to topological hindrance practically saturate at low fractions of the second species (see Appendix D 3).
To improve our conceptional understanding of the decrease of the maximal density with increasing numbers of lanes, we can estimate the jamming density ρ max ðWÞ. To this end, we consider groups of lattice sites with equal site index μ but different lane index i, which corresponds to a column of our two-dimensional lattice. Those columns can be occupied by one up to W particles. The respective probabilities in a jammed system are unknown. A rough estimate is given by a maximum entropy approach, which assumes that all numbers (up to W) are equally likely to occur. Then, the maximal density for a system with W lanes can be estimated as where HmðWÞ denotes the Wth harmonic number. This scaling argument, despite its simple nature, turns out to be very accurate for the case of symmetric particle mixtures δ ¼ 1=2. While simulations show that the assumption of uniformly distributed numbers of particles per column is per se incorrect, deviations from the uniform distribution for highly and sparsely occupied columns seem to cancel out, leading to a correct mean value. An intuitive justification for this is that, for example, a completely filled column is compatible with all configurations to its left. Nonetheless, it is very unlikely to occur because only a very restricted number of states may appear to its right. On the other hand, a column with a single particle in it can occur to the right of any other state, but it is only compatible with a few configurations to its left. Hence, this state is very rare as well. On average, high and low occupancies are both suppressed. The effects balance out, making the mean a good approximation. Note here that the jamming of particles below full occupation shows strong similarities with the jamming transition observed in the Biham-Middleton-Levine model [67]. The latter is a (often deterministically described) model of two perpendicularly crossing pedestrian flows that exhibits a sharp transition from a phase with finite flow to jamming. A rigorous mathematical treatment of this transition is, however, lacking for the Biham-Middleton-Levine model [68]. In summary, our analysis suggests that topological hindrance becomes much more significant for large numbers of lanes, where we expect it to dominate the complete dynamics. In the Sec. V, we use a novel analytic method to construct the complete current-density relation and, thus, relate the impact of topological hindrance to arbitrary densities. For this theory, the maximal density is also a key parameter as it incorporates the full impact of the number of lanes into the relation.

V. A THEORY FOR TOPOLOGICAL HINDRANCE
In the following, we present an approximation method that is designed to quantify the impact of topological hindrance for arbitrary densities. This will enable us to derive the current-density relation of our model, which successfully takes those correlations into account that lead to the failure of mean-field arguments. To do so, we show that the current-density relation can be split into a meanfield contribution, as derived in Sec. IV, and a correction term that depends on the local density. This correction term can be associated with the inaccessibility of an empty lattice site and therefore quantifies topological hindrance at a certain density. We present a construction method for this hindrance function that, in turn, allows us to compute the current-density relation.

A. The hindrance function
First, we split the particle currents into a mean-field contribution and a correction term. The particle density is independent of the lane index i, such that we find Here, covðn X i;μ ; n Y l;ν Þ ¼ hn X i;μ n Y l;ν i − hn X i;μ ihn Y l;ν i denotes the covariance of two occupation numbers. According to the master equations, Eqs. (4) and (5), the particle currents J S and J T must be conserved on each lane in the stationary state. In addition, the current must be independent of the lane index i due to the rotational invariance imposed by the system's cylindrical symmetry. Hence, the total current J ¼ J S þ J T is conserved within the whole system: for arbitrary lanes i ∈ f1; …; Wg. The first summand is identical to the mean-field current derived in Sec. IV, whereas the covariances account for correlations. We define the hindrance function H at site μ as Then, using definition (9), the current can be rewritten as Note that no approximations have been made so far. Only symmetries and conservation laws have been employed.
Considering the current-density relation of the TASEP, J ¼ ρð1 − ρÞ, we find that strict accessibility of empty lattice sites is reflected by the 1 in the second factor. In contrast, Eq. (10) reveals that, in our model, an empty lattice site can be accessed only with probability 1 − H μ . Thus, the hindrance function H can be viewed as a correction to the accessibility of empty lattice sites and reflects topological hindrance in the system. As the current-density relation, Eq. (10), defines an unclosed set of difference equations, we have to employ approximation methods to proceed with our analysis. Often, moment closure techniques are applied to get an approximation for higher moments and thus the covariances. This would, in our case, also fix the hindrance function that depends on them. However, based on our discussion of topological hindrance and the maximal density in Sec. IV, we expect correlations to be long-ranged and not constrained to certain subsegments of the lattice. Therefore, such methods do not seem promising, especially for large numbers of lanes. Instead, we derive a theory similar to the considerations of Lighthill and Whitham [36,39] for our model based on Eq. (10). Specifically, we consider the limit of large system lengths L ≫ 1. Then, we can replace the discrete site index μ by a continuous spatial variable x ≔ μ=L ∈ ð0; 1. We perform a Taylor expansion in x, which we truncate at first order in the lattice spacing ϵ ≔ 1=L ≪ 1. This is justified for small spatial variations of the density profile [69]. Indeed, our stochastic simulations confirm this assumption for large aspect ratios L=W [see Fig. 8(b) and Sec. V D]. Taking everything together, the result is a current-density relation that depends on the local density but not on its gradient, The current is spatially conserved (and therefore does not explicitly depend on x), and hence the dependence of the hindrance function H on the lattice position must be implicit via the local density ρ. Thus, our current-density relation can be written as At this point, the explicit dependence of H on ρ is still unknown. To make progress, we derive an approximating function for the hindrance function H that captures its physical properties. This function can then be used to predict the current [70].

B. Constructing the current-density relation
Since the current-density relations of the TASEP and our model show similar curve characteristics, we expect their qualitative phase behavior to be comparable because of the extremal current principle. We split the approximating function for topological hindrance H into two corresponding regimes: H LD is associated with the monotonically increasing part of the current-density relation, and H HD with its monotonically decreasing part. These two functions can be expanded independently around extremal density values for which the behavior of H can be derived using mean-field arguments. Each extremal case isolates one central aspect of transport limitation: first, the lack of particles for ρ → 0, and second, jamming in the case of ρ → ρ max . In the intermediate regime, both transport limiting factors come into play.
For vanishing densities ρ → 0, particles do not interact and are thus uncorrelated. In this limit, we do not expect any topological hindrance, and therefore, the hindrance function H vanishes. Furthermore, we can compute the derivative of H by calculating the change in current due to the addition of particles in the mean-field approximation (see Appendix B). These two conditions allow us to determine the coefficients of an expansion of the function H LD up to linear order in ρ: Interestingly, the result is independent of the number of lanes W and depends on the species fraction δ in the simplest nontrivial way that fulfills the species exchange symmetry δ → 1 − δ.
In contrast to the low-density limit, particles are totally correlated for ρ → ρ max . A particle can only be located at a certain lattice site if the respective target site is occupied. For this limit, we can again find the corresponding derivative of the hindrance function using mean-field arguments (see Appendix B). Note that Hðρ max Þ is known from the definition of the maximal density ρ max itself. Again, we can determine the coefficients of an expansion up to linear order: Since the maximal density depends on the number of lanes, as discussed in Sec. IV, the same holds true for the highdensity approximation of the hindrance function. This is a major difference compared to the low-density limit. However, the low-and high-density components both show the simple dependence on δð1 − δÞ.
A closer look reveals that linear approximations for the low-and high-density regimes do not, in general, intersect on the interval ½0; ρ max . This would result in a discontinuous current-density relation that is unreasonable. We conclude that an expansion up to at least second order in ρ is necessary. Therefore, we have to look for further physical conditions to uniquely determine the additional coefficients of the expansion. This can be achieved by imposing a differentiable transition between the highdensity and low-density regimes. At this intersection point, neither a lack of particles nor jamming is the limiting factor for transport, and hence, the transition should take place at the density ρ MC , which corresponds to the maximal current. Furthermore, both derivatives must vanish at ρ MC . This can be translated into three conditions: Combining all properties derived for the hindrance function, we can uniquely approximate H LD and H HD up to second order in the bulk density [71]. The resulting hindrance function is shown in Fig. 6(a) for a system with two lanes and species fraction δ ¼ 1=2; the corresponding current-density relations for multiple numbers of lanes W are shown in Fig. 5. As can be seen, the results of our approximations are in very good agreement with the stochastic simulations and capture the main characteristics of the current-density relation. At this point, let us also emphasize that the theory can be applied for general δ and W. It does not involve any explicit dependence on the number of lanes W but depends on ρ max , which, in turn, depends on W and also on δ. Therefore, ρ max is a key quantity in our theory, which determines the currentdensity relation and hindrance function and thus quantifies topological hindrance in general.
Since the current is very well approximated with this method, the exact relation given by Eq. (9) allows us to infer the underlying correlations. This leads to a prediction for the covariances, covðn S i−1;μ−1 þ n T i;μ−1 ; n i;μ Þ ≈ ρHðδ; W; ρÞ; TWO-SPECIES ACTIVE TRANSPORT ALONG … PHYS. REV. X 8, 031063 (2018) where we have accounted explicitly for the dependence of the hindrance function on the control parameters W and δ. This prediction of the covariances is confirmed by the stochastic simulation shown in Fig. 6(b), as expected from the method's accurate prediction of the current.
In summary, our analysis shows that modifications of the current-density relation in our system are caused by particle correlations. These correlations are accounted for by a hindrance function H, which quantifies the degree of inaccessibility of empty lattice sites and therefore quantitatively characterizes topological hindrance. We construct the hindrance function in terms of a gradient expansion. This expansion is determined by physical constraints with respect to the high-and low-density limits, as well as the transition to the maximal-current phase. Numerical simulations validate our method.

C. Phase behavior of collective two-species transport
So far, we have computed the relation between the average density and particle current but have not addressed how these observables are connected to the system's control parameters. In this section, we determine the response of the current and density to a change in the control parameters α and β and thereby derive the complete phase diagram. As mentioned above, the extremal current principle implies the existence of the three phases as in the TASEP (low-density, high-density, and maximal-current phase), due to the similarity in the form of the currentdensity relation. Then, for various driven lattice gas systems, current conservation at the respective boundary is used to relate the system's current and density to in and out rates. Because of nontrivial correlations in our system, it is, however, not possible to simply treat the boundaries as reservoirs with density ρ L ¼ α and ρ R ¼ 1 − β as is done for the TASEP. In Appendix C, we present ways to account for correlations at the edges of the system. This involves determining the response of the current and density to the in rate α in the low-density phase and to the out rate β in the high-density phase.
The LD-MC and HD-MC transitions take place when the respective current matches the maximal bulk current. The latter can be computed with the methods outlined in the previous section. Analogously, the LD-HD transition takes place when both currents are identical. The resulting phase diagram for the special case of two lanes and a species fraction of δ ¼ 1=2 is shown in Fig. 7. The theoretical boundaries are calculated based on the analytically determined value for the maximal density for this choice of parameters (see Appendix A) and hence are derived without any free parameter. The results agree very well with the data obtained from stochastic simulations. Changing δ does not affect the existence of phases but interpolates between the boundaries known for the TASEP and those presented in Fig. 7. Increasing the number of lanes W lowers the maximal density (see Fig. 4) and increases the parameter range of the maximal-current The theory for the currentdensity relation (dashed lines) agrees very well with simulation data (symbols). Although deviations from the TASEP currentdensity relation (solid line) increase with a growing number of lanes, the three regimes of the current-density relation are preserved: the low-density regime (positive derivative with respect to ρ), the high-density regime (negative derivative), and a unique maximum. Note that the theory depends on ρ max , which in turn depends on W and δ. For the theoretical lines, a fit-free approach based on the analytic result for ρ max (see Appendix A) was employed for two lanes. For three and four lanes, results were taken from the simulations to ensure that deviations originate from the theory itself not the estimate of ρ max . The particle current refers to an arbitrary lane of the system in order to allow for a comparison to the TASEP. phase, while the corresponding current is reduced. Hence, as predicted by the extremal current principle, given the number of lanes and species ratio, the phase diagram is topologically equivalent to the one obtained for the normal TASEP [72]. Note that our theory holds for general values of the system width W and the species fraction δ.

D. Nonequilibrium pattern formation
The construction of the hindrance function shown above relies on an expansion of the density ρ in terms of the spatial variable x ¼ μ=L. Specifically, we assume that gradients in the density are negligible, which is justified for slowly varying density profiles. While this assumption as illustrated in Fig. 8(b) leads to accurate results in most cases, we observe surprising exceptions to this behavior for systems with a small aspect ratio L=W: If we sufficiently decrease the system length L or increase the system width W, our stochastic simulations reveal spatial oscillations of the stationary average particle distribution. As shown in Fig. 8(a), the density ρ i;μ ¼ hn i;μ i oscillates along the longitudinal direction (Greek index μ) with a wavelength equal to the width of the system W. Because of rotational symmetry, these wavelike patterns of the density profile are equally present for each lane i. Furthermore, oscillations are sustained for thousands of lattice sites almost without any decay; therefore, they show remarkable robustness. This behavior is in stark contrast to the TASEP, which (except for boundary layers) exhibits a constant density profile. In general, pattern formation is rarely observed in lattice gas models and has so far been found predominantly in the form of segregation or localization effects [62,[73][74][75].
As stated above, the continuous approximation Eq. (12) is, by virtue of its construction, incapable of describing these varying density profiles. Also, the species current fraction δ may, in this case, significantly differ from the species density fraction (see Appendix D 2). It is, however, worth noting that our theory still provides a good approximation for the current-density relation in terms of the total system density (i.e., when averaging over an oscillating density profile) for symmetric species fraction δ ¼ 1=2. To describe the density profile itself, it is necessary to include a dependence not only of the local density but also its spatial change in the hindrance function H. The problem is comparable to the boundary layers of the TASEP, which cannot be captured by a first-order continuous theory but are predicted by exact or higher order solutions [8,44,76]. In future work, it would be interesting to address if such a construction or some alternative analytic method is capable of reproducing these intriguing patterns.

VI. ROBUSTNESS AND BIOLOGICAL RELEVANCE
To establish a generic theory for two-species transport and to enable an analytic study of topological hindrance, we previously made several assumptions that do not necessarily hold in a physical or biological context. To probe the robustness of our results and the relevance of topological hindrance further, we now turn to extended versions of our model that account for biomolecular features of molecular motors.
First, kinesins are not likely to strictly follow a unique pathway; i.e., lane switching does not occur at every step but stochastically. Adding this modification to our model mitigates effects such as a strict maximal density and phase FIG. 7. Phase diagram for the special case W ¼ 2, L ¼ 4096, and δ ¼ 1=2. Colors denote density (left panel, blue) and current (right panel, orange) obtained from stochastic simulations. Three different phases can be identified: the LD phase where current and density only depend on the in rate, the HD phase where current and density only depend on the out rate, and the MC phase that is independent of both boundary parameters. The theory (bold lines) accurately predicts the phase boundaries and is valid for general δ and W as soon as the maximal density is known. Here, ρ max was taken from the analytic approach (see Appendix A). transitions since any inaccessible site of the original model can now be accessed, in principle, on long timescales. Nonetheless, topological hindrance also has a significant impact in this case. To demonstrate this, we perform simulations with a four-lane system in which laneswitching particles switch lanes with a probability of θ ¼ 0.2 and track their respective lane otherwise. The results are shown in Fig. 9. While this system does not strictly follow the current-density relation as predicted by our theory, it demonstrates that topological hindrance has a strong influence on collective transport as the particle current is heavily suppressed: For example, the average particle current is reduced by more than 50% at densities around the maximal density of the original model. For higher lane numbers, current reduction is even more pronounced. This shows that topological hindrance can also have substantial influence on systems where particles rarely switch lanes but that are composed of many lanes.
To challenge our results further, we perform stochastic simulations of an extended model that accounts for other important biomolecular features of molecular motors: In addition to the stochastic lane switching described above, particles may attach (rate ω X a ) and detach (rate ω X d ) randomly along the lattice (Langmuir kinetics) [77,78] and occupy two lattice sites to account for a dimeric structure [79] of most molecular motors (see illustration in Fig. 10). In analogy to our analysis of the minimal model shown in Fig. 2, we increase the particle attachment rate ω S a of the spiraling motor species while keeping the attachment rate of the tracking species ω T a constant. An example for the corresponding average particle current per lane is depicted in Fig. 10. Even in this greatly expanded model, topological hindrance causes a significant decrease in the particle current already at low fractions of the second species. In particular, the impact of topological hindrance follows the same principles as discussed in Sec. IV: It increases rapidly with an increasing number of lanes, and it sets at already small fractions of a second species. Taking all results together, we conclude that topological hindrance is most likely also of relevance in biological contexts, and the theory described here allows one to estimate its impact. Another interesting question in this context is whether a phenomenon similar to the density patterns that we observe in the main model exists in a biological context. As this strongly depends on the system considered, at this stage, we cannot provide a complete analysis. Nevertheless, to give an outlook on the impact of model extensions on pattern formation, a brief discussion is given in Appendix D 4.
Moreover, we would like to add that random particle attachment and detachment may, in general, modify our results. In particular, we expect that adding Langmuir kinetics to our model will lead to the emergence of additional phases [77,78] such that our results-while providing a mathematical foundation-cannot be transferred directly [80]. Whereas a complete analysis of the jamming behavior is out of scope of this work, the analysis of our fully extended model (Fig. 10) strongly suggests that whenever transport occurs over significant length scales, topological hindrance can be expected to significantly impact the particle current, in general.
To provide a further link from our theoretical work to biological systems, we finally address how our results could be tested experimentally. Two central predictions that our model makes in this context are (1) current reduction and (2) periodic particle jams. In case (1), the average particle current in the system with two species is expected to be significantly smaller compared to a single-species system (see Fig. 10). This could be measured experimentally using time-resolved single-molecule techniques. In case (2), topological hindrance is related to spatial correlations in the particle arrangement and movement, which may lead to patterns in the particle density. These could be resolved experimentally using techniques described previously, for example, by Maurer et al. [81]. One hurdle to overcome is to construct suitable total internal reflection fluorescence (TIRF) microscopy setups that do not obstruct the helical pathways of the spiraling motor species, similar to those described in Refs. [20,82]. For models similar to the one presented in this work but with more realistic particle dynamics (e.g., random particle attachment or detachment, nonstrict lane hopping, dimeric particles), we expect current suppression due to topological hindrance to be stronger for parameters that lead to bigger differences in both curves (shaded area). Thus, our theory provides qualitative insights in more complex models also. Note that the simulations suggest a noninjective current-density relation at the onset of jamming of the original model. Other parameter values are δ ¼ 1=2, L ¼ 16 384.

VII. SUMMARY AND CONCLUSION
In this work, we have studied collective transport of two particle species with distinct gaits and therefore different directions of motion on a cylinder. As a key result, we found that the presence of different stepping modes enhances steric interactions and compromises particle dynamics in ways not seen in single-species systems. This additional hindrance changes the macroscopic transport behavior in a way that-to the best of our knowledge-cannot be accounted for by previous models or analyses. The combination of a lane-tracking and a lane-switching species gives rise to unexpected phenomenology in the following fashion: (1) Systems with a mixed population of both species always jam at densities far below full occupation. (2) Closely related to that, a system with a mixed population is always characterized by a lower average particle current at a given total system density. (3) The proportion of the second species present has a crucial influence on this behavior: Current reduction sharply sets in at very small fractions of the second species and practically saturates for fractions larger than δ ≈ 5%. (4) In contrast to single-species transport, the system width W is an important determinant of collective dynamics. Increasing W rapidly decreases the density at which particles jam and also rapidly increases the degree of current reduction. (5) The average particle distribution (density profile) shows wavelike patterns for small aspect ratios L=W of the system. This again contrasts with single-species transport and many transport models, as those typically exhibit a spatially constant or linear density profile. The above observations can be traced back to the following microscopic origin: While for transport with a single species each position can be equally accessed and vacated from a certain number of directions, our model shows a more intricate connectivity. A single lattice site may now be accessed from two different directions but can be vacated in one direction only. This simple change in the network topology has substantial influence on collective behavior. It creates inaccessible lattice sites and slows down particle motion due to intrinsic bottlenecks at points where trajectories intersect. This concept, that steric interactions are amplified by the network topology, forms the basis of our theoretical approaches, which sheds light on the following issues.
(1) We provide a detailed analysis of the jamming process. The decrease of the jamming density with increasing system width W is well approximated by ρ max → lnðWÞ=W for large W and symmetric species mixtures δ ¼ 1=2. This explains the relevance of the system width W for topological hindrance and suggests convergence to a vanishing jamming density for large system widths. We specifically compute the jamming density in the case δ ¼ 1=2 and W ¼ 2, which validates our understanding of the jamming process. (2) Even at densities well below the jamming transition, many lattice sites are inaccessible. We show that inaccessibility is quantified by a hindrance function, which naturally arises from particle correlations. In doing so, we obtain the current-density relation for arbitrary species fractions and system widths. This method is not restricted to our specific system. It reproduces the correct current-density relation of the TASEP, and we expect it to be applicable to other TASEP-like systems that exhibit cylindrical symmetry. Dimers either enter the system at the left boundary at rates α S and α T or randomly attach to the lattice at respective rates ω S a and ω T a . Lane-switching particles may either attach with both legs in the same lane at rate ð1 − θÞω S a or, alternatively, with feet in neighboring lanes at rate θω S a . Particles exit the system at the right boundary at rates β S and β T or stochastically detach from the lattice at rates ω S d and ω T d . (b) Simulation results show the current for systems with different numbers of lanes W without any spiraling motion (empty symbols) and with rare stochastic lane switching (filled symbols) of θ ¼ 0.1. If no lane switching is possible, the current is not affected by the number of lanes, and it reduces to the single-species TASEP with Langmuir kinetics and dimeric particles [77][78][79], as identical hopping rates for both species are chosen. However, if the second species is allowed to switch lanes, a rapid current reduction similar to the one described above for the main model can be observed. Note that the current is no longer conserved within the system, such that a spatial average over a lane is taken. Other parameter values include (3) The current-density relation is used to compute the phase diagram of our model. In this way, we globally relate central macroscopic observables, namely, the average particle current and distribution, to the model's control parameters-the particle in and out rates. An important conclusion of our study is that the reduction to one dimension is a priori not possible for collective transport of differently moving agents along a cylinder. Although dynamics for each species alone are effectively one dimensional and are successfully described by mean-field methods, the joint system deviates from this behavior. Our analysis overcomes this limitation and establishes a framework for mixed populations comprising distinct particle species. It further identifies topologically amplified hindrance as a central determinant in this context. The derivation of the hindrance function not only quantifies the strength of topologically amplified hindrance but also implements an effective mapping to a one-dimensional system. This mapping is useful in two ways. First, it reduces the complexity of the mathematical problem by integrating out the transversal spatial dimension. Second, it allows one to infer possible microscopic interactions whenever only a one-dimensional projection is visible. For example, in a biological context, it is still technically challenging to resolve the motion of molecular motors below the diameter of their respective filament. Hence, the projection of particle positions on the filament contour is the standard information accessible in experiments. This holds particularly true in experimental setups with a multitude of moving particles, where super-resolution techniques are not always feasible. As our work shows, it is a highly nontrivial task to relate macroscopic information back to individual microscopic interactions, and our theory might yield valuable insights in this respect. But our considerations also raise the question of how much actual biological systems, such as molecular transport along microtubules, are affected by topological hindrance. While some of our specific results (e.g., a strict maximal density or phase transitions) are weakened by processes like random particle attachment or detachment and nonstrict lane hopping, we still expect topologically amplified interactions to play a critical role. Our analysis regarding model modifications in Sec. VI suggests that, whenever transport occurs over significant length scales and timescales, topological hindrance most likely changes collective behavior substantially. An important conclusion in the biological context is that particle jamming due to topological hindrance takes place at significantly lower particle densities as compared to single-species transport: Figure 10 suggests that even slight deviations from the straight stepping behavior (10% in this case) reduce the current to half of the value as compared to the single-species case at identical cytosolic particle concentrations. Thus, jamming transitions might also be more relevant in vivo than previously thought. Therefore, it might be revealing to probe this phenomenon experimentally in multispecies setups and to speculate about its biological implications.
In these ways, our work has various consequences for intracellular transport, and it yields a generic theory for collective motion of differently moving agents. On a theoretical level, our analysis is in line with results of the TASEP for single-species transport. But since the TASEP has relevance beyond the field of transport processes, our results also form a bridge to other intensively studied fields of statistical physics: Our model constitutes a simple system far from thermal equilibrium that can be investigated in great detail with the theory developed in this paper. It exhibits a wide range of phenomena that are insufficiently understood and have rarely been observed in analytically well-accessible models such as driven lattice gases. First, in contrast to many other driven lattice gas models, the average particle distribution exhibits complex pattern formation. While a detailed mathematical study of these patterns is beyond the scope of this work, our construction of the hindrance function and the corresponding reduction to one dimension can form a basis for a more rigorous investigation. Second, our model features a nontrivial jamming transition. While jamming also has severe and intricate implications for the TASEP and similar models, the origin of jams in our model is not only the consequence of mere overcrowding but also of the spatial arrangement of the constituents. This is of interest as jamming transitions arise in many different problems, such as traffic flow, granular media, or glassforming liquids [83]. Nonetheless, it is still a demanding and ongoing task to establish a concise theoretical framework for these phenomena [84,85], although it is an open question whether a general connection exists [86,87]. The challenge involved in finding a mathematical formulation of jamming processes also becomes evident in view of the Biham-Middleton-Levine model [67] for traffic flow: Although the model has been extensively studied for more than 15 years, a rigorous analytic approach to describe its jamming transition is still elusive and considered to be an unsolved mathematical problem [68,88,89]. In this light, our study opens a new perspective on jamming processes as well as pattern formation in terms of simple, analytically accessible models, and it further elucidates the intriguing principles of collective behavior emerging in nonequilibrium systems. One of the critical differences of our model in comparison to one-species transport is complete jamming of particles at densities far below full occupation. In this section, we focus on these completely jammed configurations and determine their number of inaccessible lattice sites and thus the maximal density ρ max .
Complete jamming occurs only for vanishing out rates. One might expect that it is possible to focus on the case β ¼ 0. This results in a much easier stochastic process, which we refer to as the filling process. Our stochastic simulations, however, show that a more careful analysis is required when β vanishes: The average system density exhibits different values for arbitrarily small out rates β → 0 than for out rates that strictly equal zero, β ¼ 0. For example, for W ¼ 2 and δ ¼ 1=2, our simulations show that the density converges to a constant value of ρ max ≈ 0.73 for small but nonzero values of β (tested for different values of β up to 10 −10 ). Setting β ¼ 0, the system realizes a different density of ρ filling ¼ 0.79. Our analysis of completely jammed configurations hence has to treat the cases β ¼ 0 and β → 0 separately.
To obtain a relation between the densities corresponding to β ¼ 0 and β → 0, it is of relevance to understand how such a discontinuous behavior of the density can arise. Ultimately, it can be traced back to an instability of subconfigurations that arise in the filling process. A possible configuration is shown in Fig. 11(a). The filling process creates configurations that contain particles that do not block other particles [crossed particles in Fig. 11(a)]. Consequently, it is possible to remove these loose particles from a completely jammed system without allowing further particle motion. In the stationary state, however, each particle that is removed from the system has to be compensated by a new particle, on average. Therefore, loose particles are suppressed in any stationary state that exhibits a nonvanishing, yet arbitrarily small, stationary current. The system undergoes a discontinuous phase transition when the out rate changes from zero to a finite value.
Based on this reasoning, it is possible to derive analytic expressions of the average system density created by the filling process ρ filling and the maximal (dynamic) system density ρ max corresponding to β → 0. For simplicity, we restrict our discussion to W ¼ 2 lanes and δ ¼ 1=2 throughout this section. Our analysis includes the following steps: We first focus on the filling process with β ¼ 0. For this stochastic process, it is possible to obtain an exact analytic solution for the probability of finding an arbitrary configuration and therefore for ρ filling . The result further allows us to compute the probability p loose of finding a loose particle in a configuration created by the filling process. When a nonvanishing current is established, the presence of loose particles is suppressed, and the density decreases correspondingly, which will be computed in the last step of our derivation.

The filling process
The filling process is defined as the stochastic process arising for β ¼ 0. Typically, in jammed configurations, the system properties are independent of the in rate α, which suggests that we should also consider the limit α → 0. Then, we can assume that particles arrive at the right lattice end independently, which further simplifies our FIG. 11. Illustrations for the filling process, β ¼ 0, α → 0. (a) During the filling process, loose particles (crossed particles) are created. Loose particles are particles that can be removed from a jammed configuration without allowing further particle motion. Thus, they are suppressed in any stationary state that exhibits an arbitrary nonvanishing current. (b) If a column is fully occupied (box in group 1), the next particle (shaded particles) cannot pass this column; therefore, a truncation of the state space with respect to particles to the right is possible. This is referred to as double closure. If a linear array of at least two particles of the same species is followed by a particle of the other species, either a double closure is created or a single occupied column can be identified (boxes in groups 2 and 3), beyond which none of the following particles (shaded particles) can pass. This also allows one to truncate configurations to the right of this column, and it is referred to as a single closure. (c) Transition matrix of the stochastic process for β ¼ 0, α → 0. The state space reduces to arrays of particles of the same species that follow a double (solid end) or single closure (shaded end). A single particle can only arise after a double closure, as single closures immediately create arrays of length two. In (a) and (b), a gray area denotes the system boundary.
discussion. Our mathematical formulation should describe the sequence of lattice occupations when particles are added sequentially. Specifically, a (discrete) time step is defined by a new particle getting stuck in the system. A naïve way to denote the state space would be to enumerate the 3 2×L possible lattice occupations. However, this leads to a very irregular and high-dimensional state space that precludes an exact solution. Figure 11(b) shows several configurations that are jammed at the lattice end but empty at the front. Light-colored particles denote the possible positions at which the next particle can end up. In all cases, the next particle cannot pass the column that is marked with a box. This means that its final position is independent of the configuration to the right of this box. We can identify such a column where particles cannot pass in two different cases: First, whenever a column is filled completely [box in group 1, Fig. 11(b)], no more changes can be made to the configuration to the right of this column. We refer to this column as double closure. Second, whenever a linear array of particles of the same species is followed by a particle of the other species, a column that cannot be passed can be identified. In this case, either a double closure may be created or a new kind of closure arises. Specifically, if no double closure is created, we can always identify a halfoccupied column [boxes in groups 2 and 3, Fig. 11(b)], beyond which no further change is possible. We refer to this column as a single closure. The above reasoning then suggests the following truncation scheme for the state space: Whenever a closure occurs, we can truncate the state space with respect to the occupation to the right of it. This significantly reduces the complexity of the state space. The remaining possible configurations are linear arrays of a single species to the left of either of the two different types of closures. Given this reduced state space, it is possible to write down the resulting (infinite-dimensional) transition matrix that characterizes the underlying stochastic process. The result reads where the states are enumerated as shown in Fig. 11(c). This matrix encodes the evolution of probabilities for occupations after a closure and the frequencies at which double as well as single closures occur. For an infinite system, the occupational probabilities of the stochastic process will converge to a steady state, which corresponds to the eigenvector of the stochastic matrix M with eigenvalue one. The result reads With these steady-state probabilities for lattice configurations, we can derive an exact result for the average particle density in an infinite system. We compute the frequencies of blocks of subconfigurations with a closure at the left being added to the system. These frequencies can be determined by the transition rates and the steady-state probabilities. For example, an array of two T particles following a double closure [second configuration in the enumeration scheme of Fig. 11(c)] can be "closed" by an arriving S particle such that a stable block dimension 2 × 2 occupied by three particles is created. As the probability of an S particle arriving at the respective site is 1=4, the corresponding frequency for this block to occur is 1=4P 2 . We weigh these frequencies with the block lengths and densities, which yields the system density in the limit L → ∞. Our analytic approach results in a value of ρ filling ¼ 0.79. Notably, we made no approximations in this computation. Consequently, it should coincide with the system densities for β ¼ 0 and α → 0. Indeed, simulations with respective parameter values confirm our result.

Maximal dynamic density ρ max
As illustrated above, the presence of unstable, loose particles will lead to a different density as soon as β goes from zero to arbitrarily small values. Therefore, computing the number of loose particles that arise in the filling process is key to determine the density decrease as soon as the system becomes dynamic. In the first step to compute the maximal density ρ max , we determine the probability of finding a loose particle. In the second step, we show how to correctly replace the corresponding occupational density and thereby estimate ρ max .
Loose particles can only arise whenever both lanes are occupied at a given position μ, i.e., within a doubleoccupied column. The probability of double-occupied columns created by the filling process is easily computed from the filling density. For further convenience, we introduce a new random variable c μ ∈ f0; 1g that equals 1 if the μth column is fully occupied (irrespective of the particle species, i.e., for n X 1;μ ¼ 1 ∧ n X 2;μ ¼ 1) and 0 otherwise. Then, where Pðc μ ¼ 1Þ is the probability of finding a fully occupied column. Then, a loose particle occurs in a double-occupied column in two different situations: first, when this column is followed by a column to its left that is occupied by only one particle. This will always result in a loose particle in the double-occupied column. The second situation is if the double-occupied column is followed by another double-occupied column to its left where a particle from each species is present. Then, the two particles of the left column are directed to a single site, which creates a loose particle next to this site. To compute the corresponding probabilities, we use the exact solution of the filling process. In detail, we can determine the probability of a double closure being followed by a fully occupied column. This scenario is identical to the conditional probability of a double-occupied column being followed by another double-occupied column and, using the above solution of the filling process, is given by pðc μ ¼ 1jc μ−1 ¼ 1Þ ≈ 0.64. In half of these cases, the fully occupied column on the left is populated by particles from different species. This corresponds to the second case described above and therefore creates a loose particle. If the double-occupied column is not immediately followed by another double-occupied column [corresponding to (1 − pðc μ ¼ 1jc μ−1 ¼ 1Þ)], the first case as described above occurs, and again a loose particle is created. Taking all these considerations together leads to Here, we use the definition of the conditional probability, Knowing the probability of a site being occupied by a loose particle, it seems tempting to subtract the corresponding occupations from the filling density ρ filling . But particle rearrangements in a dynamic system will continuously create new loose particles. Hence, we have to replace sites occupied by a loose particle with sites of a suitable effective density ρ eff .
To this end, consider a hole (i.e., an empty site) that propagates from right to left through the system. In each hopping step, the hole might either hit a column that contains a loose particle with probability p loose or one that contains only stable particles with probability 1 − p loose . As only loose particles contribute to a modification of the density, we are only interested in the first case. Let the average time between two consecutive holes arriving at any given position in the steady be denoted by T. A fully occupied column that contains a loose particle will, on average, remain in this state for a time T before a hole arrives. Then, the hole either hits the loose or the stable particle in this column. If the stable particle is removed, the density is not changed, and the column remains at full occupation until the next hole arrives, i.e., for another time T. On the other hand, the column remains at half occupation for a time T if the loose particle is removed. This leads to an effective (time-averaged) density of Hence, the filling density ρ filling and the maximal dynamic system density ρ max are related via The result is ρ max ≈ 0.74, which is in very good agreement with stochastic simulations that yield ρ max ≈ 0.73. This validates our understanding of topological hindrance for jammed configurations. Because of approximations in the above calculation, our computed value for ρ max is, opposed to the one for ρ filling , not an exact result. Specifically, we implicitly assume that for nonvanishing β → 0, occupations are created in the same way as they are in the filling process β ¼ 0. The good quality of our result therefore indicates that, as a hole propagates through a jammed system, empty lattice sites are filled by particles in a mostly uncorrelated fashion. Note that the filling density does, in principle, depend on the particle in rate α. If particles are injected with a high frequency, mutual interactions-while traveling to the jammed end-lead to correlations and sorting effects that ultimately cause varying values for ρ filling . This dependency of the filling density on α is, however, irrelevant for our computation of the maximal (dynamic) density ρ max : Correlations that build up during the particles' motion can also be neglected in the case β → 0 as the system is almost completely jammed.
The theoretical arguments of this section addressing a transition from the filling density to the maximal dynamic density are further illustrated with the following consideration. Figure 12 shows the temporal evolution of the average system density with a very small out rate, starting from an empty system. The out rate was chosen such that, on average, the system completely fills with particles before a first particle exits. Indeed, we observe that the density gradually increases until ρ filling is reached. In particular, the density overshoots the maximal dynamic density ρ max . As soon as particles leave the system, loose particles are removed, which causes the average density to decrease and fluctuate around ρ max . These fluctuations are also explained by our above argumentation: As single holes travel through the system, particles are reordered, and loose particles are transiently created, which leads to an (also transient) increase of the system density.

APPENDIX B: DERIVATIVES OF THE HINDRANCE FUNCTION
In the following, we calculate the derivatives of the hindrance function for the two extremal cases ρ → 0 and ρ → ρ max , as used in Sec. V. To do so, we determine the change in current dJ in response to a change in density dρ at a given lattice site and equate it with the derivative of the current-density relation, Eq. (12), In the mean-field approximation, a lattice site is occupied with a particle with probability ρ irrespective of the neighboring lattice sites. We find an S particle at any given lattice site with probability δρ and a T particle with probability ð1 − δÞρ, respectively. Increasing the density by a value dρ, we have two contributions to the change of the current: first, a decrease since particles are more likely to be blocked, and, second, an increase because the additional density itself may contribute to the overall current. Both contributions can again be split into two separate cases. We are either dealing with a focused state [see Fig. 3(b), group 1], where two particles are blocked simultaneously, or with the case where only one particle is affected. Two particles from neighboring lanes are focused on the same lattice site with probability ρ 2 δð1 − δÞ, resulting in a current reduction of −ρ 2 δð1 − δÞdρ. A single particle is blocked with probability ρδ(1 − ρð1 − δÞ) þ ð1 − δÞρð1 − δρÞ. Here, ρδ(1 − ρð1 − δÞ) is the probability of only one S particle being blocked and ð1 − δÞρð1 − δρÞ of only one T particle being blocked. The same contributions can be derived for the increase in current due to the removal of a single particle. Summing this all up, we obtain This expression reduces to the derivative of the TASEP current-density relation for δ ¼ 0 and δ ¼ 1.
The mean-field result is not correct, in general, for our system but becomes exact in the limit ρ → 0. Using Eq. (B1), our result for the derivative of the hindrance function reads As, by definition, HðρÞ=ρ → H 0 ð0Þ for ρ → 0, we obtain FIG. 12. Temporal evolution of the system density for α → 0, β → 0, β ≪ α. As the number of simulation steps increases (x axis), particles enter the system, "pile up" at the end, and create an average density close to ρ filling ≈ 0.79 (solid line). As soon as particles leave the system, the presence of loose particles is suppressed, which leads to a drop of the average density until the maximal dynamic density ρ max ≈ 0.73 is reached. The dashed line depicts the result of the analytic approach for ρ max described in this section. Here, the density is obtained from an ensemble and spatial average. A single realization (inset) shows that the average system density (i.e., spatial average) oscillates around ρ max , which can be explained by large reordering events if stable particles are removed after several loose particles are extracted. Simulation parameters include α ¼ 0.02, β ¼ 10 −4 , δ ¼ 1=2, We proceed in a similar way for the case ρ → ρ max . However, instead of an uncorrelated system, we face total correlations between particle occupations; a site is only occupied if the one in front is also. Filling the last hole in a jammed configuration will not increase the current, but it will prevent the last possible particle motion. This means that we are left with those parts of Eq. (B2) that have a negative contribution. Evaluating the corresponding equation at the maximal density ρ max , we find Again, using Eq. (B1) and the definition of the maximal density Hðρ max Þ ¼ 1 − ρ max , we obtain for the derivative of Hðρ max Þ. Note that, as opposed to the limit ρ → 0, this is not an exact result but should be interpreted as a refined mean-field approximation.

APPENDIX C: DERIVATION OF CURRENT AND DENSITY RESPONSE TO THE CONTROL PARAMETERS IN THE HIGH-AND LOW-DENSITY PHASES
In this section, we derive the dependence of the current and density on the in rate α and out rate β. As for the TASEP, collective behavior in the low-density phase only depends on the in rate α, whereas the high-density phase is determined by the out rate β.

Low-density phase
For many driven lattice gases, current conservation is used to relate the system's current and density to in and out rates. For the case of the TASEP, the low-density phase is dictated by an influx J ¼ αð1 − ρ 1 Þ. Furthermore, there is no boundary layer at the left lattice end in the low-density phase such that the density at the first site equals the bulk density. Thus, because of current conservation, ρ 1 ¼ α for the TASEP. For our model, however, we face a new problem. Simulations suggest that, even in the low-density phase, the occupation of a first lattice site is not equal to the bulk density and is also not identical to the reservoir density α. This is due to nontrivial correlations between occupations at the left lattice end (see Fig. 13), which we compute in the following. Indeed, the density at the first site is always slightly higher than α. Also, when, for example, the in rate of T particles, α T , is increased, the density of S particles at a first lattice site ρ S i;1 increases and vice versa. The latter observation, in particular, contradicts equating ρ S i;1 ¼ α S . This behavior can be understood in the following way. Consider, for example, a T particle at the first lattice site. This T particle can never prevent an S particle to hop in front of it; the T particle can be blocked by the S particle, but it cannot block the S particle. Therefore, this T particle will, on average, stay longer at the first lattice site as compared to the TASEP, and the corresponding density ρ i;1 will be higher than α. The exact flux-balance condition for T particles, Eq. (5), reads Again, correlators of the form hn X 1 i;1 n X 2 i;2 i lead to an unclosed set of equations. To overcome this problem, we approximate second moments. Note that all particles enter the system from the left reservoir in an uncorrelated fashion. The way T particles can block and are blocked by other T particles follows the same scheme as in the TASEP. Therefore, it is reasonable to assume that correlations balance (as for the TASEP). Thus, we factorize hn T i;1 n T i;2 i ≈ hn T i;1 ihn T i;2 i for the same species. This is different for hn T i;1 n S i;2 i, where a TASEP-like balance of the blockage and getting blocked no longer holds true as argued above.
Explicitly writing down the master equation for hn T i;1 n S i;2 i in the steady state leads to As particles that enter the system are uncorrelated, we can set hn T i;1 n S j;1 i ¼ ρ T i;1 ρ S j;1 for i ≠ j. Furthermore, we are FIG. 13. Second moment hn T i;1 n S i;2 i of occupations of different species at the first two lattice sites as predicted by our theory (lines) and stochastic simulations (symbols) for the special case δ ¼ 1=2 and W ¼ 2. Deviations between theory and simulation only occur close to the phase transition at α ≈ 0. 22. interested in a solution that is valid in the low-density phase, where we expect α to be small and ρ i;1 to be of order (but not necessarily equal to) α, i.e., ρ i;1 ∝ OðαÞ. If correlations are weak, we can also assume that hn i;μ n j;ν i ∝ Oðα 2 Þ and hn i;μ n j;ν n k;ξ i ∝ Oðα 3 Þ. The transitions from the low-density phase to the maximal-current and high-density phases take place at lower values of α as compared to the TASEP (cf. Fig. 5), which ensures very low values for the density in the low-density phase in our model. Then, upon truncating at third order in α, we arrive at an equation for hn T i;1 n S i;2 i that is valid for low densities: Here, we also assume that ρ S i;2 ≈ ρ S i;1 [90]. An analogous equation can be derived for hn S i;1 n T i;2 i. Inserting the results in Eq. (C1) and its analogue for the influx of S particles, we find Equations (C4) can be interpreted as a two-species TASEP, where particles from different species cannot block each other when entering the system. This means that, for low densities, the possibility of overtaking particles from the other species when entering the system effectively acts as if no exclusion was present for different species; particles from the same species still exclude each other. Equations (C4) can be solved for ρ T i;1 and ρ S i;1 , which fixes the current and bulk density in the system.
For the case of identical in rates, the solution of Eqs. (C4) is of a particularly simple form, and the total current is given by Even though this result is independent of the number of lanes, the actual current in the system can never become higher than the respective maximal current of the bulk, which, in turn, depends on the number of lanes. As soon as the bulk limits transport, the transition to the maximalcurrent phase is triggered. For higher in rates, the conservation of the maximal current dictates the density at the first lattice sites. A comparison with simulation data for different numbers of lanes is shown in Fig. 14. The lowdensity theory agrees with the data up to the point of the phase transition. The exact value of the in rate α at which the transition into the MC phase takes place depends on the number of lanes.
In summary, Eq. (C4) correctly describes the behavior of our system in the low-density phase that is found to be independent of the number of lanes W of the system. Dependencies on the number of lanes are only relevant for phase transitions that can be derived by applying the extremal current principle.

High-density phase
In the high-density phase, the system's current is dictated by the right boundary and determined by an outflux of particles J ¼ βρ L . Here, ρ L denotes the density at the last lattice site for an arbitrary lane i. Because of current conservation, this outflux must equal the bulk current, where ρ denotes the density in the bulk. As for the lowdensity phase, we lack knowledge about the value of the density at the dominating boundary, in this case ρ L . Clearly, the mean-field result ρ L ¼ 1 − β cannot be used. This becomes evident when we consider jammed systems that exhibit a maximal density ρ max that, in general, does not equal 1 − β. However, we expect the density at the right lattice end to be of the order of the bulk density, ρ L ≈ ρ. Using this assumption in Eq. (C6) we obtain a unique solution for the particle current. The result is shown in Fig. 15 and is in agreement with stochastic simulations.
(a) Current plotted against the in rate α in the LD phase for different numbers of lanes W. For low in rates, the current increases monotonically until it becomes constant at a point that depends on the number of lanes; the respective current corresponds to the maximal current of the bulk. Low-density theory agrees with simulation results for the LD phase. The low-density theory is independent of the number of lanes. (b) Density at the first lattice site ρ i;1 plotted against the in rate α for different numbers of lanes. In the LD phase, the density is dictated by the left boundary and agrees with the low-density theory. After the transition to the maximal-current phase, the density is determined by conservation of the maximal current. Again, the corresponding value for the maximal current depends on the respective number of lanes, which leads to a distinct behavior for different numbers of lanes in the MC phase.

Rotational symmetry of the density profile
In the main text, we discussed how density profiles have to be identical for every lane because any stationary state has to be unique. However, the time needed to reach this stationary state can be very long since the system may get stuck in metastable configurations that break symmetry. Figure 16 shows that indeed all lanes admit identical density profiles in our simulations. This demonstrates that simulation times are chosen sufficiently long to ensure convergence to the (rotationally invariant) stationary state.

Relation between the species fraction of the current and the density
In Sec. IV of the main text, we have shown that, in the mean-field approximation, the current species fraction δ is identical to the density species fraction ρ S =ρ. However, this result does not hold true in general. To assess the difference between the species fraction of the density and that of the current, we perform stochastic simulations for a wide range of parameter values. As illustrated in Fig. 17, our simulations show that deviations between both quantities are relatively small (below 10%). This justifies using δ as the general species fraction in this paper.
There are, however, important exceptions where this equivalence cannot be employed. If we regard systems with a small aspect ratio L=W and significantly asymmetric species mixtures, the species fractions of current and density strongly differ. Figure 18 shows differences of up to 60% between these quantities in stochastic simulations, with corresponding parameter values. In all of these cases, the naïve assumption ρ S ≈ δρ clearly fails, and the individual densities have to be considered separately.

Current-density relation for a system with a large number of lanes
The current density relation for a large number of lanes shows very drastic modifications as compared to the TASEP (see Fig. 19). Most of the possible density values are not realized at all in the stationary state. Furthermore, the maximal current is far below the one known from the TASEP. As for low numbers of lanes, a very small fraction of the second species (2%-5%) is already sufficient to drastically reduce transport efficiency.
FIG. 15. Average system current as a function of the out rate β in the high-density phase for different numbers of lanes for the special case δ ¼ 1=2. Comparison between stochastic simulations (symbols) and theoretical predictions (lines). The current increases monotonically with increasing out rate. For sufficiently high out rates, the system enters the maximal-current phase, and the current becomes constant. The assumption of ρ L ¼ 1 − β with a corresponding out flux of J ¼ βð1 − βÞ (uppermost dashed line) fails to describe our model. FIG. 17. Relative deviation of the density species fraction from δ for the special case of W ¼ 2 lanes and L ¼ 4096 sites. Color denotes simulation results for the relative deviation and thus the error that is made when assuming a density species fraction of value δ. Lines indicate the phase boundaries as predicted by our theory.

Robustness of the density patterns with respect to model extensions
One of the most intriguing features of our main model is the formation of regular density patterns. However, while this is clearly an interesting phenomenon to study in the context of nonequilibrium physics, it is an open question whether such patterns may occur in a biological systems. While a rigorous and complete analysis of pattern formation in this context is out of the scope of this work, we would like to provide an outlook on potential changes due to different model extensions. Figure 20 shows how the density patterns change as compared to Fig. 8(a)  Since no differences between monomers and dimers are observed, we restrict our discussion to the dimeric case for the rest of this section.
Including nonstrict lane hopping leads to a spatial decay of oscillatory patterns [ Fig. 20(b)]. On the other hand, random particle attachment smoothens the density profiles globally [ Fig. 20(c)]. Finally, random particle detachment (even for comparably high detachment rates) still results in oscillatory patterns [ Fig. 20(d)]. This can be explained in the following way: To establish density patterns, particles of different species need to jam periodically because of particles they have encountered before. For the main model, the typical length scale needed for this to happen is the width of the lattice. In the case of nonstrict lane hopping, the period of circulating once around the system is stochastic, which gradually destroys the "phase" of our density oscillations. If we add random particle detachment to the system, length scales are not changed: Some particles may leave the system, which makes jams less likely but does not alter the length scale on which those jams appear. Random particle attachment critically changes this phenomenology: A certain lattice position no longer corresponds to a fixed run length of a specific particle. Therefore, the stationary density profile becomes translationally invariant, which, in turn, has to suppress density patterns when averaged over large time windows. FIG. 18. Relative deviation of the density species fraction from δ for the special case W ¼ 13 lanes and L ¼ 512 sites. Color denotes simulation results for the relative deviation and thus the error that is made when assuming a density species fraction of value δ. Lines indicate the phase boundaries as predicted by our theory.
FIG. 19. Current-density relation obtained from stochastic simulations with W ¼ 13 lanes and L ¼ 16 384 sites. The emerging current of an arbitrary lane (blue) shows a dependence not only on the average system density but also on the fraction of lane-switching particles, δ. The system jams at a maximal density ρ max ðδÞ (bold gray line); densities between ρ max and full occupation are not realized in the stationary state. (c) As soon as dimers are allowed to attach randomly on the lattice (ω S a ¼ ω T a ¼ 5 × 10 −5 ), the system loses its characteristic length scale and the density profile becomes smooth. (d) Random dimer detachment (ω S d ¼ ω T d ¼ 10 −3 ) still allows for the formation of density patterns. While the wavelength stays the same, the amplitude of the oscillations is reduced. Other parameters include W ¼ 13, L ¼ 512, α ¼ 0.6, β ¼ 0.2, δ ¼ 1=2.