{'fraction_non_alphanumeric': 0.03662590525940876, 'fraction_numerical': 0.014602873085598956, 'mean_word_length': 4.500544188071397, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Driven Transport on open filaments with inter-filament switching processes

We study a two filament driven lattice gas model with oppositely directed species of particles moving on two parallel filaments with filament switching processes and particle inflow and outflow at filament ends. The filament switching process is {\it correlated} such that particles switch filaments with finite probability only when oppositely directed particles meet on the same filament. This model mimics some of the coarse grained features observed in context of microtubule (MT) based intracellular transport, wherein cellular cargo loaded and off-loaded at filament ends are transported on multiple parallel microtubule (MT) filaments and can switch between the parallel microtubule filaments. We focus on a regime where the filaments are weakly coupled, such that filament switching rates scale inversely as the length of the filament. We find that the interplay (off)loading processes at the boundaries and the filament switching process leads to some distinctive features of the system. These features includes occurrence of variety of phases in the system with inhomogeneous density profiles including localized density shocks, density difference across the filaments and bidirectional current flows in the system. We analyze the system by developing a mean field (MF) theory and comparing the results obtained from the MF theory with the Monte Carlo (MC) simulations of the dynamics of the system. We find that the steady state density and current profiles of particles and the phase diagram obtained within the MF picture matches quite well with MC simulation results. These findings maybe useful for studying multi-filament intracellular transport.

We study a two filament driven lattice gas model with oppositely directed species of particles moving on two parallel filaments with filament switching processes and particle inflow and outflow at filament ends. The filament switching process is correlated such that particles switch filaments with finite probability only when oppositely directed particles meet on the same filament. This model mimics some of the coarse grained features observed in context of microtubule (MT) based intracellular transport, wherein cellular cargo loaded and off-loaded at filament ends are transported on multiple parallel microtubule (MT) filaments and can switch between the parallel microtubule filaments. We focus on a regime where the filaments are weakly coupled, such that filament switching rates scale inversely as the length of the filament. We find that the interplay (off)loading processes at the boundaries and the filament switching process leads to some distinctive features of the system. These features includes occurrence of variety of phases in the system with inhomogeneous density profiles including localized density shocks, density difference across the filaments and bidirectional current flows in the system. We analyze the system by developing a mean field (MF) theory and comparing the results obtained from the MF theory with the Monte Carlo (MC) simulations of the dynamics of the system. We find that the steady state density and current profiles of particles and the phase diagram obtained within the MF picture matches quite well with MC simulation results. These findings maybe useful for studying multi-filament intracellular transport.

I. INTRODUCTION
One-dimensional driven diffusive systems, unlike their equilibrium counterparts are known to exhibit boundary induced phase transitions [1][2][3]. Such systems have also served the purpose of providing a framework for studying wide class of driven biological phenomenon ranging from transport across biomembranes [4], to transport on individual cellular filament [5][6][7][8][9] and cytoskeletal filament network [10].
Filament based intracellular transport involves oppositely directed motors, which use multiple arrays of cytoskeletal filaments, to actively transport cellular cargoes such as mitochondria, endosomes, and pigment granules [11,12]. It has been observed that long-distance cellular cargo transport on microtubule (MT) filaments is achieved by sets of oppositely directed motor proteins, e.g; dynein and kinesin, which attach to the cellular cargoes and transport them actively along these filaments. The transport itself is determined by different processes at play at the molecular level, e.g; the motor processivity, directional switching dynamics of the cargo carried by the motors, the underlying filament organization, the (un)binding characteristics of the motors to the filament and the boundary input(output) rate of cargoes at filament ends [13]. It has also been observed that both long distance regulated transport [13] and phenomenon of jamming arise out of the collective action of the these motor proteins [14,15].
One of the approaches to study transport in such systems has been to describe it in terms of coarse-grained driven lattice gas models wherein the MT filament is considered as a 1-d lattice, the interactions between the transported cargoes are included via excluded volume effect and the underlying driven stochastic dynamics due to various processes incorporated in the description [8]. Some of the previous theoretical attempts have focused on the interplay of stochastic directional switching mechanisms, directional hopping of individual cargoes and the effect of the input and output of the cargoes on the boundaries on a single filament [8,9,16]. However for a variety of biological situations such as axonal transport in neurons, cargo transport takes place on a multiple parallel array of MT filaments [15,17]. For example, in vitro studies on cultured neurons have revealed that cargo switching between neighbouring filaments occurs in axonal transport of mitochondria on neurons [18]. On theoretical grounds, it has been argued that even without considering the effects of the boundaries, the interplay of the translation process on filaments and the filament switching processes can manifest in form of a phase transition between a an inhomogeneous jammed phase of the transported cargoes and a freely flowing phase with homogeneous cargo density in each filament [19]. Thus studying the role of multiple filaments in determining the transport properties of such systems is of considerable importance.
Driven transport on parallel lattices have been studied theoretically in different contexts [19][20][21][22][23][24][25][26][27][28][29][30][31][32], and the particle switching dynamics between adjacent lanes have also been taken into account explicitly in some cases [19-22, 27-30, 32]. In this paper we will focus on how the transport along two parallel filaments is affected by the interplay of boundary inflow and outflow of particles at the filament ends and filament switching dynamics of particles. Before we proceed describing the mode in detail, we wish to highlight a few aspects of transport that have been observed in the context of intracellular transport : (a) Experimental studies, such as the one on endosomal transport on MT reveal that cellular cargoes can switch between neighbouring filaments [18,33]. (b) Experiments suggest that cellular cargoes traveling in opposite direction on the same MT can also cross each other and continue with their translational motion along the same filament [33]. (c) Cargo transport is also dependent on loading and offloading of cellular cargoes at the filament ends [13,34]. Motivated with these experimental observations, in the model that we study, we focus on the interplay of the boundary driven processes of particle input and output at filament ends with the active motor driven cargo translocation on the filaments and filament switching processes. Accordingly, we will consider two parallel filaments with oppositely directed species of particles which translate on the lattices with a specified hopping rate. The oppositely directed particles are also allowed to pass through each other with a certain specified rate on the same lattice. The particles are also allowed to switch between the lattices with certain finite probability only when oppositely directed species meet each other on the same lattice, so that the switching between the lane is a correlated process [19,20]. Thus implicitly we will consider that individual motors carrying the cargo have a propensity to switch between different filaments when they experience a force, when hindered by an oppositely directed particle moving on the same filament. This in turn can decrease the binding affinity of the motor to the filament and induce it to switch and bind to the neighbouring filament. Finally we will allow particles to enter and exit the filament ends with prescribed rates.
In Section II we specify the model, the dynamical rules on the lattice and the corresponding equations of motion for the system. In Section III we set up the Mean Field (MF) continuum equations for the system and analyze the boundary conditions at the ends of the lattice. In Section IV we analyze the different phases possible for the system, obtain the corresponding MF steady states density profiles for the particles and construct the MF phase diagram for the system and compare these results with the Monte Carlo (MC) simulation results. Finally in Section V we discuss these results in context of multifilament intracellular transport. We consider two parallel cellular filaments, represented by two finite and parallel one dimensional lattices of length L with N sites, with lattice spacing ǫ = L/N . The cellular cargoes transported along these filaments will be referred to as particles, and will be characterized by two different species. Specifically, along each filament, the transported cargo can either be a (+) particle which moves from left to right on the filament or a (−) particle that moves from right to left. Without loss of generality, the two filaments are labeled as 1 and 2. The instantaneous state of the system is described in terms of the occupation numbers, which indicate the spatial localization of the two species of particles on the two parallel filaments. Specifically, n + i,1 corresponds to a occupation number of a particle at site i moving to the right on filament 1. The maximum allowed occupancy at any lattice site is 1 so that each lattice site is occupied either by a (+) particle, (−) particle or is vacant, (0).

II. THE MODEL
The dynamics of this system can be expressed in terms of the movements allowed for the particles. For the sites in the bulk, in each individual filament, a (+) particle can hop from site i to site i + 1 with a rate α if that site is vacant. Similarly a (−) particle can hop from site i to site i − 1 with an identical rate α if it is vacant. If the i-th site on filament 1 is occupied by a (+) and if the neighboring site to the right, i.e; site i + 1 is occupied by an oppositely directed (−), then two different processes can occur: with rate β the two particles can swap their positions, and the (+) particle moves to the site i + 1 and the oppositely moving (−) particle moves to site i; while with rate π o 12 the (+) particle can switch the to filament-2, at the same corresponding site, identified by its index i, if the site is vacant. Similarly, a (−) particle from filament 1 at site i + 1 can switch to the site of filament 2 with the same site index i+1 with a rate µ o 12 if that site on the other filament is vacant. Identical processes that we have described for filament 1 also happens for filament 2. The rates of filament switching processes from filament 2 to filament 1 for the (+) and (−) particles are π o 21 and µ o 21 for (+) and (−) particles, respectively. Particle switching between the two filaments can be understood as arising from the stronger loading force experienced by the motor proteins which carry the cargo when they push against an oppositely directed cargo. This leads to an increase in the rate of motor detachment (along with the cargo) from the filament and offers the possibility of a subsequent reattachment of the molecular motor to the neighbouring site on the other filament [35]. Since we are interested in the regime where the filaments are weakly coupled, we probe the regime where the filament switching processes compete with loading and offloading processes at filament boundaries. We systematically implement it by choosing the filament switching rates at individual lattice site such that they scale inversely with system size so that we have, Although the bulk processes are analogous to those introduced in Ref. [19,20], which focused on the collective behavior of cargoes moving along filaments with a closed ring morphology and with overall particle number conservation, we will concentrate on the behavior of such particles for open filaments. In this configuration the overall particle number is not conserved, and the motion on incoming and outgoing particles at the filament ends must be accounted for.
For the boundary sites at the filament ends, a (+) particle can enter the left end of the filament (with site label i = 0) of filament 1 with a rate α + 1 if it is vacant, and it can leave from the right end of the filament (with site label N − 1) of filament 1 with a rate β + 1 . Similarly a (−) particle can enter the right end of the filament 1 with a rate α − 1 if it is vacant and it can leave from the left end of the filament with a rate β − 1 . Similar processes also occur in filament 2 with the corresponding rates being α + 2 , β + 2 , α − 2 and β − 2 respectively. All the the dynamic processes that characterize the model are schematically depicted in Fig. 1.

III. MEAN FIELD EVOLUTION EQUATIONS
The time evolution for the average occupation number for the two oppositely directed species along each individual filament can be expressed in terms of gain and loss terms arising from translation and filament switching processes. As described in Ref. [20], these terms involve the averages of the local occupation numbers for the particles at each site as well as various combination of averages of two-point correlators of site occupation numbers, which account for the role of particle correlations in the particle collective dynamics.
We set α = β = 1 and choose π 12 = π 21 = µ 12 = µ 21 = π, which correspond to a symmetric scenario where the propensity to switch filaments is the same for both (+) and (−) species and it is symmetric about filament 1 and filament 2. The Mean Field (MF) evolution equations are obtained factorizing the two point correlators of the occupation numbers. The continuum MF evolution equations are derived by rescaling the total length L to 1 and letting N → ∞ so that ǫ → 0 [6,20].
Correspondingly, p 1 (x), p 2 (x), n 1 (x) and n 2 (x) are then the average densities as a function of the relative position in the filament, x. The MF continuum equations in the bulk can be expressed as, where we have displayed terms up to first order in ǫ. The corresponding expression for the currents of each species read, A. Steady State profiles From Eqs.(1)-(4) , the steady state profiles corresponding to the continuum MF evolution of the molecular motors satisfy, (12) which govern the bulk profiles of (+) and (−) particles in the two filaments. One can rewrite them in terms of the fluxes of the total number of particles and its difference introducing Since in the bulk the total current for the combined system comprising of the two filaments is isolated, J p and J n have to be spatially constant. Subtracting Eq. (10) from Eq. (9) we identify an additional conservation law, which corresponds to the sum of the absolute values of the currents of opposite species along one filament. As shown in Appendix B, it is useful to reorganize the three independent conserved quantities, J p , J n and J 1 and which remain spatially uniform in the bulk, in terms of three new parameters, Using these new conserved quantities, the equations for the density profiles in the bulk can be decoupled, and express them in terms of one single density variable. Specifically, we can write where η ± p1 , µ ± p1 and ν ± p1 are functions of p 1 alone. Their explicit functional dependence is provided in Appendix B. Similarly, we can get decoupled differential equations for n 1 , p 2 and n 2 , The explicit form of the decoupled differential equations which govern the density profiles and the relevant solutions are discussed in Appendix B. In Appendix B, Eq.(B1-B4), provides the mathematical expression for the quantities present in Eq. (16)(17)(18)(19). Appendix B also describes the different sets of density profiles that can be obtained from the previous set of equations. We have found that there are 16 different branch solutions corresponding to the decoupled differential equation. The choice of a particular boundary condition corresponding to a particular phase, restricts the possible choices to 8. Finally as discussed in Appendix B, physical considerations such as bounds on the physical value of density selects an unique solution to these differential equations for each set of prescribed boundary conditions. Therefore, the continuous MF equations determine the density profiles of the particles in the two filaments, once the boundary conditions are prescribed. Although Eqs. (16)(17)(18)(19) for the different species densities decouple in the bulk, they are coupled through the boundary condition of Eqs. (13)(14)(15). In the next subsection we discuss the set of boundary conditions satisfied by these differential equation corresponding to a particular phase.

B. Boundary conditions
The allowed phases that characterize the state of transport on the two filaments is controlled by the particle input and output at the boundaries. The steady state density profiles are determined by first order differential equations. This is due to the fact that the diffusive contribution is of higher order in the lattice spacing, ǫ, and their contribution drops in the continuum limit, ǫ → 0. As a result, the system cannot fulfill, generically, the input and output boundary conditions and boundary layers are expected [6]. One then must determine which of the fluxes at the filaments' ends determine the bulk density profiles and under which conditions coexistence between density phases can develop along the filaments.
For the system composed of two filaments, we have 8 different particle entrance or exit rates at the boundaries, that we express as α ± 1/2 ,β ± 1/2 . However there are only 4 boundary conditions to be specified for the complete solution of the steady state differential equations. In order to figure out the possible physically relevant boundary conditions, it is useful to build upon the boundary conditions that are satisfied by a closely related 1-D lattice gas model [1,2], which has the similar translocation dynamics along the filament as our model, but which does not allow for inter-filament exchange processes. In fact the model discussed in Ref. [1,2] is exactly the same as ours for the particular case of π = 0, which corresponds to a situation where the inter-filament switching dynamics of the particles is turned off. Since for our case, the filaments are weakly coupled thus it is expected that the boundary conditions satisfied for our two-filament system are the various possible combination of the boundary conditions that are satisfied for individual lattice for the case studied in Ref. [1,2]. However we would like to stress that although the boundary conditions are obtained as simple combination of boundary conditions of the individual lattices, the resultant density and current profile obtained by spatially integrating the steady state differential equation would be qualitatively different due to the lattice switching term in the bulk.
We enumerate the possible combination of the boundary conditions and the resultant phases for each of those particular combinations.
Filament 1 in HL phase and filament 2 in HL phase ((HL) 1 − (HL) 2 ): When the bulk current of (+) matches with the output current of (+) at the right boundary and the input current of (−) matches with the bulk current of (−) at the right boundary both for filament 1 and 2, the resultant phase corresponds to a situation where the (+) particles are in high density(H) phase and (−) particles are in low density(L) phase (refered to as HL phase), in both the filaments. The boundary conditions that are satisfied in the continuum limit are, where J + 1R , J − 1R , J + 2R and J + 2R refer to the currents for (+) and (−) particles in filaments 1 and 2 at the right(R) boundary respectively, while p 1R , n 1R , p 2R and n 2R refer to the densities of (+) and (−) particles in filaments 1 and 2 at the right boundary. Using Eq.(20), we can find the expression of the boundary densities at the right end of both filaments in terms of the entry and exit particle rate By symmetry there can be another phase where both in filament 1 and 2, (−) are in high density phase while the (+) are in low density phase and the current conditions are satisfied at the left boundary. Further one can find a situation where for filament 1 the current at the left boundary for (+) and (−) matches with the bulk current, while for filament 2, the current at the right boundary for (+) and (-) matches with the bulk current. Similarly, there exists a phase for where for filament 1, the current at the right boundary for (+) and (−) matches with the bulk current, while for filament 2, the current at the right boundary for (+) and (−) matches with the bulk current. For all these 4 different phases, the structure of the boundary condition is exactly similar.
Filament 1 in LL phase and filament 2 in LL phase : : When the bulk current of (+) particle matches with the input current of (+) at the left boundary and input current of (−) matches with the bulk current of (−) at the right boundary both for filament 1 and 2, the resultant phases corresponds to LL phase in both the filaments. The boundary conditions satisfied by the currents are, while the expression for the currents at the other boundaries read Using Eq. (22) and Eq. (23), we can express the boundary densities as a function of entry and exit rates and the currents [2], As opposed to the continuity of the overall particle fluxes at filament's ends due to particle conservation, J pL = J pR and J nL = J nR , the current on the left and right end of an individual filament track will in general differ due to particle filament switching. In order to determine the densities for the two types of particles in this phase, we will assume that the (+) current in left boundary of filament 1 equals the (+) current in right boundary of the same filament, Analogously, we equate the currents at the left and the right boundary in filament 2, J ± 1L = J ± 1R and J ± 2L = J ± 2R . This is a reasonable assumption because tracks are weakly coupled, as has been checked with MC simulation.
This fact allows us to obtain an expression for the boundary densities in this phase. Explicitly, from Eq. (22) and Eq.(24), the boundary densities for both the lanes read , These coupled algebraic equations for each filament track can be numerically solved to get the corresponding boundary densities from which the density profiles can be numerically derived.
Filament 1 in LL phase and filament 2 in HL phase : ((LL) 1 − (HL) 2 ): When for filament 1, the bulk current of (+) matches with the input current of (+) at the left boundary and input current of (−) matches with the bulk current of (−) at the right boundary and the bulk current of (+) matches the output current of (+) at the right boundary and the input current of (−) matches with the bulk current of (−) at the right boundary in filament 2, the resultant phase for the system corresponds to LL phase in filament 1 and HL phase in filament 2.
In this case we again assume the continuity of the fluxes separately for the two particle types along filament 1, Accordingly, the equations for the boundary densities are similar in form to those expressed in Eq. (25) and the boundary densities for filament 1 can be obtained numerically as discussed earlier.
For filament 2, the densities are determined by Eq. (21). Due to the symmetry in the swapping rates between the two filament, a second phase with analogous structure is feasible,where filament 1 in HL phase and filament 2 in LL. There is still a further symmetry associated with the HL phase in any of the two filaments, i.e; if the current of (−) and (+) at the left boundary matches with the currents in the bulk, the structure of the boundary conditions remains unaltered. In order to illustrate this, consider HL phase in a particular filament, then the corresponding densities are determined by Eq.(21) and the boundary condition at x = 1 is satisfied, so that (+) is in high density phase and (−) is in low density phase. But analogously we could have a situation where the boundary condition at x = 0 is satisfied with (+) particles in Low density phase and (−) in high density phase. This situation would correspond to a different overall phase, but structure and the form of boundary density would be the same as Eq . (21) with the only difference that the indexes of (+) and (−) in the expression for the boundary density in Eq .(21) is interchanged. Thus this structure of boundary conditions would correspond to 4 distinct phases.

IV. PHASES AND PHASE DIAGRAM
Since the boundary fluxes control the particle fluxes in the bulk, once we have determined the expression for the densities at the filament boundaries in terms of the input and output rates at the filament boundaries, the MF density and current profiles in the bulk can be determined from Eq. (16)(17)(18)(19). However, unlike the case of model in Ref. [1,2] where the steady state density profiles and the corresponding phases in the bulk are solely determined by the boundary fluxes, the density profiles and phases now emerge from the interplay of the boundary processes and particle filament switching. This distinctly alters the nature of the density profiles and the topology of the phase diagram for the system under study.
Specifically, a first major consequence of particle exchange between filaments is that the density and current profiles in the bulk are no longer spatially homogeneous. In fact for certain range of input and output particle rates, the competition between the bulk and the boundary processes can result in density shocks along the filaments which are localized in the bulk.
Moreover, particle change between filaments also allows for phase coexistence in the bulk apart from the pure phases which satisfies only one set of boundary conditions. This feature of phase coexistence occurs only because the current profiles in the bulk for a set of boundary conditions are not homogeneous so that different MF solutions for the current, intersect each other at specific spatial location in the bulk of the system. The phase coexistence in the bulk happens when the currents of the different MF solutions arising out of different boundary conditions match at a point in the bulk of the twofilament system. In that case part of the system obeys one set of the boundary conditions while the other half obeys another set of boundary condition and these set of solutions are joined by the condition of continuity of current at a particular location in the bulk of the lattice. In general the system selects the set of MF steady state solution for which the corresponding current is minimum. This holds true as long as the MF current profiles do not attain the maximal current value in the bulk. In this paper we have restricted our analysis to the region of parameter space of entry and exit rates of particles for which the condition for maximal current is not reached.
As a result of these new features, the topology of the resulting phase diagram changes qualitatively with respect to the collective behavior in the absence of inter-filament particle exchange. In the following subsection we first describe the procedure to find the MF density and current profiles and determine the resultant phase. Subsequently we discuss the topology of the resultant phase diagram for this system and obtain the equations for the phase boundaries.

A. Density profiles and emerging phases
In order to find the MF density and current profiles using Eq. (16-19), we have to first determine the three independent conserved currents in the system e.g; J p , J n and J 1 . Subsequently, we have to provide the appropriate values of the densities at the boundaries to completely specify the solutions for the individual species.
(HL) 1 −(HL) 2 phase: Here, first we determine the values of the boundary densities at the right end of both filaments, e.g; p 1R , n 1R , p 2R and n 2R using Eq. (21). Thus J p , J n and J 1 can be determined at the right boundary. The entire density profile can now be found out by evolving the MF solution from the left end of both filaments using Eq. (16)(17)(18)(19). In Fig.2, we show the comparison of the MF profile with the MC simulations for this phase. It illustrates that both the density and the current profiles in the two filaments are not spatially homogeneous in contrast to similar phases in Ref. [1,2]. A similar procedure can be used to find out the profiles for the corresponding (LH) 1 − (LH) 2 phase.
(LL) 1 −(HL) 2 phase: The values of the boundary densities, p 1L and n 1R , can be determined by numerically solving Eq.(25) for the first filament. For the second filament we use Eq.(21) to determine p 2R and n 2R . Thus, both J n and J 2 are identified at the right boundary. To determine J p we use the method of successive iterations. In the first iteration we set p 1R = p 1L and obtain J p to get the density profiles. Again these profiles are not uniform due to inter-filament switching processes. Consequently, the density value obtained at x = 0 by evolving Eq. (16) from x = 1 is not same as p 1L . Hence, we take the difference between these two values at x = 0, and this difference is added to p 1R and evolve it to get the new density profile. We repeat this process until convergence is reached [36]. This procedure allows then to derive the entire density profile by evolving the MF solution from the right end of both filaments using Eq. (16)(17)(18)(19). An analogous procedure is used to identify the profiles for the case of (LL) 1 − (HL) 2 phase, (LL) 2 − (LH) 1 phase and (LL) 1 − (LH) 2 phases.
(LL) 1 − (LL) 2 phase: For this phase coexistence, we need to determine the values of the boundary densities p 1L , n 1R , p 2L and n 2R by solving Eq. (25) and use these values to extract the associated boundary fluxes, J p and J n . To this end, we assume n 1L = n 1R , a symmetry that holds in the absence of particle filament exchange [2]. This relation allows us to determine J 1 and consequently fix the value of n 2L . Again, we use the process of successive iterations for determining both n 1L and n 2L . The entire profile can be found out by evolving densities from the left end of both filaments using Eqs. (16)(17)(18)(19)) (See Fig.3). (HL) 1 − (SL) 2 : In this phase while filament 1 is in HL phase, there is phase coexistence in the second filament, such that at the right end of the filament, the boundary condition corresponding to HL phase is satisfied while at the left end of filament 2, the boundary condition corresponding to LL phase is satisfied and this phase is characterized by a density discontinuity of (+) species in filament 2 which results in a shock profile in that filament. For phase coexistence, the profile has to be such that the current for (+) corresponding to HL solution equals the current for (+) for the LL solution at a particular spatial location between the filaments ends. This current continuity condition follows from the fact that in the bulk the total added current of the two filaments has to be conserved as there is no particle exchange with the environment. The phase coexistence can be thought of as a mixture of the (HL) 1 − (HL) 2 phase with (HL) 1 − (LL) 2 . Accordingly, as one changes the parameters corresponding to input and output rate of particles, starting from a pure (HL) 1 − (LL) 2 phase, the system can evolve into (HL) 1 − (HL) 2 through an intermediate phase coexistence region corresponding to the (HL) 1 − (SL) 2 phase in parameter space, where an incipient shock of (+) particles originating at the right end of the filament eventually reaches the left end of the filament on change of parameters in the phase diagram. In order to determine the density profiles for this case, we make us of the fact that we know p 1R , n 1R , p 2R and n 2R because the right end of both filaments are in the HL phase. This property allows us to identify the three independent conserved currents e.g; J 1 , J p and J n for the entire filaments and the density profile from the right end of the filament can be plotted using Eqs. (16)(17)(18)(19). For the left end of filament, having determined p 1L from Eq. 25, the remaining densities at the left boundary can determined using the three conserved currents. Now the LL profile can be simply determined starting from the left end of the filament, using Eqs. (16)(17)(18)(19). The spatial location in the bulk for which the current for this solution matches with the current for the other solution obtained for HL−HL phase defines the position of shock, as shown in Fig.4. (LL) 1 − (SL) 2 : In this phase while filament 1 is in LL phase, there is phase coexistence in the other filament. For filament 2, at the right end the boundary condition corresponding to HL phase is satisfied while at the left end, the boundary condition corresponding to LL phase is satisfied. In order to determine the density profile for this case, we use Eq .(25) to identify n 1R at the right end of filament 1. Analogously, we also know p 2R and n 2R because filament 2 is in the HL phase. At the left end of both the filaments are in LL phase so that p 1L and p 2L are known using Eq . (25). Therefore it follows that J 2 , J p and J n are known for the entire filaments and all the densities at both filament boundaries are thus determined. Now the density profiles originating from both left and right end of the filament can be plotted using Eqs. (16)(17)(18)(19) separately. The point at bulk for which the current for both these solution match defines the shock position. (SL) 1 − (SL) 2 : This arrangement, composed by three coexisting phases in the bulk, is characterized by shocks for (+)-particles in both filaments. In both filaments the right end is in (HL) 1 − (HL) 2 phase, while the re-gion close to the left ends is a (LL) 1 − (LL) 2 phase. In between a (HL) 1 − (LL) 2 phase develops. The two density shocks in the bulk separate these three regions. For the phase region adjoining the right end of the filament, the boundary condition corresponding to (HL) 1 − (HL) 2 phase is satisfied and the boundary densities are determined by Eq. (21) and the MF density profiles are obtained by evolving the MF solution from the right end of the filament using these boundary densities. Accordingly, J 1 , J p and J n are known for the entire filaments. This analysis has shown how inter-filament switching process leads to a wealth of new inhomogeneous phases, allowing for phase coexistence, and the possibility of shocks in both the filaments as opposed to the collective dynamics of transport in the absence of such filament interactions.

B. Phase boundaries
In the previous subsection we have illustrated that any pair of pure phases are mediated by a phase coexistence region in the phase diagram. Whenever the current solutions corresponding to two different phases cross each other at a particular spatial location along the filament, the system exhibits phase coexistence and the spatial location of the shock coincides with the location on the filament where the two different current solutions intersect. The system selects the combination of those set of steady state density profiles for which the corresponding current is minimum at any given spatial location in the bulk. Exploiting this insight, we can now determine the entire phase diagram and the corresponding phase boundary by using the condition when the minimum value of the current along the filaments is allowed at one of their ends. Accordingly, the phase boundary separating any two regions in the phase space corresponds to the parameter values for which the location of the intersection of the different current solutions occurs at either of the filament ends. We obtain the MF phase boundary using the conditions for allowed phase in particular parameter range of entry and exit rate of particles and compare these results with MC simulations. We will concentrate on the phase diagram when varying (α + 1 − β + 1 ) while holding the other parameters constant. This choice clearly illustrates the qualitative new scenarios that filament switching brings into collective transport. Moreover, the scheme described can be straightforwardly extended to analyze the phase diagram when varying other sets of control parameters.. Boundary between (HL) 1 − (LL) 2 phase with (HL) 1 − (SL) 2 : Here the phase boundary is determined by the condition, the arguments as in the previous case. However, the condition of current matching of (−) species on filament 2 is done at x = 0 and this determines the phase boundary between the two phases.
Boundary between (LL) 1 − (LL) 2 phase with (SL) 1 − (LL) 2 : The condition to determine the phase boundary is given by Here, J 1+ (LL)1−(LL)2 (1) is the current of (+) species in filament 1 at the right boundary when the system is in (LL) 1 −(LL) 2 phase and J 1+ (HL)1−(LL)2 (1) is the current of (+) species in filament 1 at the right boundary when the system is in (HL) 1 − (LL) 2 phase. For the (LL) 1 − (LL) 2 phase the densities -p 1L and p 2L is known from Eq. (25), which identifies J p . For the (HL) 1 − (LL) 2 phase p 1R , n 1R , and n 2R are determined by Eq. (21) and this is used to find J 1 and J n . Having obtained these fluxes, the entire density and current profile for the two phases can be found out and the current matching condition for the (+) species in filament 1 determines the phase boundary in this case.
Boundary between (HL) 1 − (LL) 2 phase with (SL) 1 − (LL) 2 : Condition for this boundary in terms of boundary currents is given by Here, J 1+ (LL)1−(LL)2 (0) is the current of (+) species in filament 1 at the left boundary when the system is in (LL) 1 − (LL) 2 phase and J 1+ (HL)1−(LL)2 (0) is the current of (+) species in filament 1 at left boundary when the system is in (HL) 1 − (LL) 2 phase. For the (LL) 1 − (LL) 2 phase the densities -p 1L , n 1R , p 2L and n 2R are known from Eq. (25), which identifies J p and J n . For the (HL) 1 − (HL) 2 phase, p 1R and n 1R are determined from Eq. (21) and this is used to find J 1 . Having obtained these fluxes, the entire density and current profile for the two phases can be found out and the current matching condition for the (+) species in filament 1 at the left boundary determines the phase boundary.
Boundary between (SL) 1 − (LL) 2 phase with (SL) 1 − (SL) 2 : The phase boundary derives from the parameters for which, is satisfied, where J 2+ (HL)1−(LL)2 (1) is the current of (+) species in filament 2 at the right boundary when the system is in the (HL) 1 − (LL) 2 phase and J 2+ (HL)1−(HL)2 (1) is the current of (+) species in filament 2 at the right boundary when the system is in the (HL) 1 − (HL) 2 phase. At the right boundary the system is in the (HL) 1 −(HL) 2 phase, and therefore p 1R , p 2R , n 1R and n 1R can be determined using Eq. (21). Consequently ,J p , J n and J 1 are known and J 2+ (HL)1−(HL)2 (1) can be found out. Since the left end of the filament 2 is in the LL phase, p 1L is known, and the entire density and current profile can be determined. From these profiles we identify the remaining flux, J 2+ (HL)1−(LL)2 (1). Matching the current for the two profiles at the right boundary determines the phase boundary.
where J 1+ (LL)1−(LL)2 (0) is the current of (+) species in filament 1 at the left boundary when the system is in the (LL) 1 − (LL) 2 phase and J 1+ (HL)1−(LL)2 (1) is the current of (+) species in filament 2 at the left boundary when the system is in the (HL) 1 − (LL) 2 phase. At the right boundary the system is in the (HL) 1 −(HL) 2 phase. Therefore p 1R , p 2R , n 1R and n 2R can be determined using Eq. (21), while n 2R can be determined using Eq. (25), which identifies J p , J n and J 1 . Using this information the entire density and current profile All the other possible combination of phases boundaries can be obtained by simply interchanging the labels of filament 1 and 2 and using exactly the same conditions for phase boundaries that have been described in this subsection.
From all these conditions, we can now build a complete phase diagram. Fig.6 and Fig.7 show the comparison of MF and MC phase diagram in different phase plane cuts, as a function of (α + 1 − β + 1 ) for fixed values of the rest of the control parameters. The MF phase boundaries that have been obtained exhibit fairly good agreement with the phase diagrams obtained by MC simulations, which shows that MF is quantitatively accurate to describe the different phases that characterize transport intros system. The contrast between Fig.6 and Fig.7 also shows that the topology of the phase plane can drastically be altered by tuning the parameters corresponding to the particle entry and exit rates although the filaments themselves are weakly coupled through the filament switching processes. Further, we also see that changing the particle entry and exit rates in one filament can affect the phases in the neighbouring filament. For instance in Fig.7 as one increases the value of entry rate of (+) particles in filament 1, keeping the exit rate of (+) particles fixed, the resultant phase in the other filament passes over from an LL phase to and SL phase with a shock developing on this other filament. The phase diagram also shows the possibility to have shock reentrant phases. For a given entry rate, the increase in the exit rate naturally favours a transition from HL to LL phases, but these transitions are modulated by the developments of shocks. As a result, in the transition from HL-LL to LL-LL is mediated by an intermediate region of SL phases, and for large enough entry rates, the LL phase is destabilized by the development of an SL phase as the exit rate increases. ] are without filament switching (π = 0) as discussed in [1,2]. [(b) and (d)] are with filament switching (π = 1). In (b); α − 1 = 0.4, β − 1 = 0.3, α + 2 = 0.8, α − 2 = 0.2, β + 2 = 0.25, β − 2 = 0.7. In (d); α − 1 = 0.2, β − 1 = 0.7, α + 2 = 0.2, α − 2 = 0.2, β + 2 = 0.2, β − 2 = 0.6. Parameter regime of (a) is same as that of (c) and parameter regime of (b) is same as that of (d).
The resulting phase diagram for this system can topologically be very distinct from the phase diagrams for two species transport in the absence of particle switching between filaments. To emphasize this fact, in Fig. 8 shows the MF phase diagram for a system in the absence of filament switching (π = 0), Figs.8a and 8c, with corresponding predictions when particle filament exchange is allowed Figs.8b and 8d keep'the same weak exchange rates between the two filaments and modify the correspond to a choice of entry and exit rates. We have chosen representative parameters to show the differences associated to particle filament exchange. For Figs.8a and 8b one illustrates that one of the new features introduced by particle exchange between filaments is the appearance of phase coexistence regions sandwiched between pure phases; a scenario forbidden when opposite particles displace along a unique filament [1,2]. Figs. 8(c) and Fig.8(d), for a different set of parameter values, show that particle filament exchange can have a deeper impact on the phase diagram topology. In this case we see that, while in the absence of filament switching dynamics the system always in a pure phase, the effect of filament switching of particles manifests as enriching the phase behaviour for the system which allows for phase coexistence and presence of shocks in the two filaments.
As illustrated in these two examples, generically we find that two pure phases are always connected by a phase region where the system exhibits phase coexistence and bulk localized shocks.

V. CONCLUSIONS
To summarize, we have studied a multi-filament driven system with oppositely directed species of particles when the filaments are weakly coupled. Particle filament switching processes constitute correlated events because particles can only swap filaments, with a prescribed finite probability, when oppositely directed particles meet on the same filament. This aspect of filament switching mimics cellular cargo switching between neighbouring filaments during intracellular transport. We find that the interplay of the entry and exit processes of particles at the filament boundaries has a profound impact in the collective organization of the two species of displacing particles, leading to a variety of new scenarios. Specifically, we have identified the development of phase coexistence on the filaments, inhomogeneous density profiles, density shocks localized in the bulk and bidirectional current flows in the system. We have developed a mean field theory (MF) to characterize these phenomena, and have shown that the steady state density and current profiles of particles and the phase diagram obtained using a MF formulation match reasonably well with the Monte Carlo (MC) simulation results.
While in this paper we have focused on the implications that weak coupling between filaments would have on transport when there are particle input and output in both filaments, it would be interesting to explore the regime of strong filament coupling, where we expect a weaker spatial inhomogeneity in the particle density profiles. Further, for many biological situations such as transport in axons, it is a priori not clear whether boundary loading and off-loading of cellular cargo happens for all the parallel filaments or for specific filaments. In such situation one needs to determine the steady state distribution of cargoes and the resultant phases when one of the filament has much higher particle entry rates than the other.
As an extreme case, we have considered a situation where one filament has closed boundaries. Starting from a random distribution of particles in both filaments, we have observed the development of phase segregation between (+) and (−) particles in the closed filament. This phenomenon happens only due to the correlated lane switching process. The resulting phase segregated state in the blocked filament does not have any flux. Following the time evolution of such a system shows that starting from a random configuration, all the vacancies are expelled and eventually the filament stops exchanging particles with the other filament as the blocked filament attains a jammed configuration, with the (+) particles piling up from the left end of the filament while the (−) particles pile up on the right end . Understanding the transition of this phase segregated jammed steady state to the steady states discussed in this paper as one slowly increases particle input and output for the filament which is initially closed at the boundaries, remains as an interesting open challenge. 5000 time swaps.
In order to determine the phase boundaries by MC simulations, we use the fact that all phase transitions between pure phases are mediated by phase coexistence regions with shocks in the density profile. The phase boundaries can be determined numerically by tracking when the location of shock reaches the filament boundary. However, due to the finite size effects of the system, the shock has a certain finite width. We determine the shock width and track the position of the midpoint of the shock to identify the phase boundaries numerically and decide when they have reached a filament end. We have used system size of N = 5000 for determining the location of phase boundaries. For a fixed β + 1 we have varied α + 1 in steps of 10 −3 and this sets the accuracy of the phase boundaries obtained numerically.
Appendix B: Choice of branches for the MF solution for density To choose a proper branch uniquely from the potential eight solutions that can be derived from the Mean Field solutions for a particular variable we have to look carefully at the density profiles of various species in a particular phase. In Eqs. (16)(17)(18)(19), notice that various combinations of η ± 's, µ ± 's and ν ± 's appear in this set of differential equations which govern the spatial density profile for each of the individual species. Each of the individual combinations of η ± 's, µ ± 's and ν ± appearing in these set of differential equations are exclusively functions of a single variable. For example if we consider η ± p1 , then it appears in Eq. (16), where the explicit form reads as η + p1 = 1 2 + 1 4 + p 1 (1 − p 1 ) − J 1 and η − p1 = 1 2 − 1 4 + p 1 (1 − p 1 ) − J 1 . Each η + p1 and η − p1 could separately be combined with each of the two different values of µ ± p1 and ν ± p1 that appear in the expressions of Eq.(B1). As there are eight possible combinations of η ± p1 , µ ± p1 and ν ± p1 , therefore there would be eight potential solutions corresponding to the choice of a particular boundary condition. Similarly there would be eight possible solutions for Eqs. (17)(18)(19).
To do the classification we have to look at the various expressions for η ± 's, µ ± 's and ν ± 's and they read as, Among the eight different branches of Eq.(16) which arise due to eight possible combinations of η ± p1 , µ ± p1 and ν ± p1 , we take that particular branch of the equation which is a combination of η − p1 , µ − p1 and ν − p1 . This branch would be classified as Sol−1. The entire nomenclature is classified in Table−I for different solutions of Eq. (16). The same nomenclature is true for the different solutions of Eqs. (17-(19).
Let us assume that filament 1 is in LL phase whereas filament 2 is in HL phase. Therefore the density of (+) and (−) particles in both the filaments and hence the variables will have the following values.