Spontaneous pulsing states in an active particle system

Spontaneous pulsing states in an active particle system

Sarah Klein, Cécile Appert-Rolland and Martin R. Evans Laboratoire de Physique Théorique, Université Paris-Sud and CNRS UMR 8627 - bâtiment 210, 91405 Orsay Cedex, France Fachrichtung Theoretische Physik, Universität des Saarlandes, D-66123 Saarbrücken, Germany SUPA, School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9 3FD, Scotland mevans@staffmail.ed.ac.uk

We study a two-lane two-species exclusion process inspired by Lin et al. (C. Lin et al. J. Stat. Mech., 2011), that exhibits a non-equilibrium pulsing phase. Particles move on two parallel one-dimensional tracks, with one open and one reflecting boundary. The particle type defines the hopping direction. When only particles hopping towards the open end are allowed to change lane, the system exhibits a phase transition from a low density phase to a pulsing phase depending on the ratio between particle injection and type-changing rate. This phase transition can be observed in the stochastic model as well as in a mean-field description. In the low density phase, the density profile can be predicted analytically. The pulsing phase is characterised by a fast filling of the system and - once filled - by a slowly backwards moving front separating a decreasing dense region and an expanding low density region. The hopping of the front on the discrete lattice is found to create density oscillations, both, in time and space. By means of a stability analysis we can predict the structure of the dense region during the emptying process, characterised by exponentially damped perturbations, both at the open end and near the moving front.


July 31, 2019

1 Introduction

Driven diffusive systems are studied within physics as models manifesting non-equilibrium states. They exhibit a richer variety of behaviour and structure of phases than the traditional equilibrium states [1] such as boundary-induced phase transitions, even in one dimension [2]. Amongst other applications, they have been used to model biological processes. For example, the first role of the totally asymmetric exclusion process with open boundaries (TASEP) was as a model of the dynamics of ribosomes along mRNA [3]. Since then the TASEP and its variations have become the generic starting point for the modelling of intracellular transport effected by molecular complexes exhibiting steric effects [4, 5, 6]. Furthermore, at a larger scale, active media, such as self-phoretic colloids or motile bacteria may exhibit motility induced phase separated states [7, 8, 9]. At still larger scales biological or artificial birds can exhibit cohesive flocking states driven by long wavelength sound modes [10, 11].

Many of these non-equilibrium models have an oscillatory or alternating nature not exhibited by equilibrium states. For example the collective direction of a flock may switch to and fro [12, 13, 14], fish may exhibit intermittency between schooling and milling [15], proteins may coherently oscillate from end to end of bacteria [16], and molecular motors may exhibit collective pulsing states [17]. An understanding of how non-equilibrium conditions allow such alternating states is only just emerging and forms an important challenge of non-equilibrium statistical physics.

Returning to the classical case of the TASEP, this model has now been generalised to multilane and multispecies versions to incorporate various biological features [4, 6]. For example, in the case of molecular motors, the switching of motors between lanes [18, 19], finite processivity [20], exchange of motors with the cytoplasm [21], and the change of motor states from active to inactive have been considered [22, 23, 24, 25, 26]. Such systems are often only tractable via mean-field (MF) theory, the simplest implementation of which results in a set of partial differential equations for the densities of different species on different lanes. It has been shown in many systems that the MF equations and furthermore a reduction to an effective single lane description gave an accurate description of steady state behaviour in a wide range of cases [19, 27]. Thus an important task is to determine if and when one can go further with a MF description and predict non-equilibrium oscillatory behaviours.

Systems in the half-open tube geometry have been of particular interest in the modelling of protrusions such as filopodia, lamellapodia, or microvilli [28, 29, 25], tubulation [30] and fungal hyphae growth [31]. The boundary condition at the closed end of the tube controls the length and may be either growing, treadmilling or no flux. The model of [23, 24, 25] which incorporates transitions between active and inactive motor states and treadmilling boundary conditions (retrograde flow) of a fixed length tube, revealed pulsing behaviour in which bunches of inactive motors move to the open end of the tube. In [32] a TASEP based model comprising 13 lanes and inspired by experimental observations of dynein dynamics on a hyphal tip [33, 34] was considered (see also [35]). For certain parameter regions the model displays so-called pulsing states in which the system alternatively fills and empties, as a result of the motion of some fronts separating low- and high-density regions.

We will study here a simplified111Actually the model considered here is a special case of the very general model introduced at the beginning of [32]. two-lane version of the latter model [32] that also exhibits these pulsing states but is more amenable to analysis. Our aim is to identify the minimal ingredients required to give pulsing states and to gain some theoretical understanding through analysis of a MF description.

The paper is organised as follows. In section 2 we introduce a two-lane two species stochastic model and present numerical simulations which illustrate that two phases of the model—low density and pulsing—exist. In section 3 we study the MF approximation to the model and show that the same two phases also exist. In the remainder of the paper we choose to concentrate on the MF version of the model. We then proceed in section 4 to analyse first the low density phase. within the MF approximation. In section 5 and 6 we turn to the pulsing phase and analyse the filling and emptying phases of the oscillatory behaviour respectively. The analysis reveals some subtle aspects of the oscillatory state and interestingly we will show that the discreteness of the lattice plays an unexpected role in the dynamics of the system. We then carry out a stability analysis of flat profiles in the high density region in section 7 and we examine the transition between the low density phase and the pulsing phase in section 8. Finally we conclude with a discussion in section 9.

2 Stochastic Particle Model

2.1 Definition of the model

In this section we define a stochastic particle model, inspired by Lin et al. [32], which exhibits pulsing behaviour in a certain parameter range. We consider a lattice composed of two parallel lanes of length with an open boundary on the left and a closed boundary on the right (see Figure 1). Two types of point-like particles hop stochastically on the lattice. Plus particles are injected at the left end of each lane with rate , and hop from site to site to the right with rate . Minus particles hop only towards the left, also with rate , and leave the system with rate at the left boundary. Additionally, particles change type (and consequently reverse their hopping direction) with rate . We refer to the type-changing process as transmutation. The particle position is exclusive so that a single lattice site can only be occupied by exactly one particle or else be vacant.

In addition to these dynamical processes, a minus particle can change lane if the next site (to the left) on the same lane is occupied, and the next site on the other lane is empty. This transition occurs with the same rate as the hopping within one lane. The plus particles on the other hand do not change lanes. Thus, the lane changes of minus particles couple the dynamics of the two lanes.

Figure 1: Schematic drawing of the possible transitions in the two lane model. Red (plus) particles hop with rate from left to right. Green (minus) particles hop with rate from right to left. Minus particles change lane with rate if the next site on their lane is occupied while it is empty on the other lane. Particles transmute with rate . Injection of plus particles occurs with rate on both entrance sites on the left, if they are empty. Minus particles located on the leftmost sites leave the system with rate .

If the ratio of hopping rate to transmutation rate is too small, plus particles will transmute and turn back towards the left end before having travelled through the whole system. In this paper, we consider only the case where is at least of order , in order to allow particles to reach the reflecting boundary after they have entered at the left.

If the ratio is too large, minus particles are impeded from leaving the left-hand end and the dynamics slows down. To keep simulations within reasonable time, we restricted ourselves to the case , i.e. at most of order .

2.2 Phases of the system

For sufficiently small the system is in a low density phase. This phase is characterised by a low density of particles in the bulk of the system, possibly with a boundary layer of accumulated plus particles at the reflecting boundary at the right-hand end (see Figure 2 (a)).

With increasing a phase transition from a stationary state with zero flux to a pulsing state occurs. This state is characterised by a cyclic behaviour. First the system rapidly fills up with plus particles. Once the total density throughout the system is close to , an emptying process begins and a front separating a high density and a low density region moves from the reflecting boundary towards the open one (see the spatio-temporal plot in Figure 2 b). When the system is almost empty, the filling process starts again. These cyclic density changes result in a change of the sign of the current, as the system empties and fills and so particles flow in and out. Just as in the low density phase, a small accumulation of plus particles persists near the wall on the right. The filling and emptying process continues indefinitely with a regular period which depends on .

Note that the pulsing phase can exist only if there is an asymmetry in the dynamics of plus and minus particles [32]. The role of this asymmetry is to induce in dense regions a net particle flow towards the exit that will be responsible for the emptying of the system. This net flow can be obtained either by allowing only minus particles to change lane — a feature that will make them more motile than plus particles in dense regions — or by producing more minus particles than plus particles through asymmetric transmutation rates, as was considered in [32]. In the latter case pulsing states can be observed even when both types of particles can change lanes. In the present paper, we only consider the case of asymmetric lane changing rules.

A third trivial phase is observed for : starting from an empty lattice plus particles fill the whole system and can never exit the system.

(a) (b)
Figure 2: (a-b) Spatio-temporal plot of one lane in the two-lane Particle Model, in the (a) low density phase, and in the (b) pulsing phase. Red (plus) particles move to the right, green (minus) particles to the left. Parameters: , for both, (a) and , (b) and .

3 The Mean-Field model

In this section we introduce a mean-field approximation to the stochastic particle model of the previous section and we will show that it exhibits the same phases and phenomenology, in particular, the same pulsing behaviour. Importantly, the MF equations are amenable to analysis and allow us a deeper insight into the mechanisms involved. Consequently, the MF model in its own right becomes the subject of our study.

3.1 Definition

We now write down the equations for the time evolution of the densities in the MF approximation where all correlations between densities are ignored. To this end we replace the boolean occupation numbers at site by the densities and , for plus and minus particles respectively. In the MF model we assume not only that density correlations are factorised, but also that the densities are the same on both lanes for a given location (this latter assumption is implicit in the fact that and depend only on but not on the lane). In this way we arrive at the following system of bulk and boundary equations.

Bulk equations:


Left boundary:


Right boundary:


In equations (1, 1b), one recognises the bulk fluxes of plus/minus particles through the link , which we denote by


Note that, as fluxes are measured to the right, the flux of plus particles gives a positive value for whereas the flux of minus particles gives a negative one for so that and are always of opposite sign.

The boundary conditions coming from (23b) can be rewritten


Equations (1-3b) can be considered as defining a dynamical system in its own right, which we refer to as the MF model in the remainder of this paper. This set of equations can directly be solved numerically.

The pulsing states which are found in the Particle Model of section 2 are also observed in the MF model, as shown in Figure 3. The Particle Model and the MF model are not completely equivalent, as can be seen for example from the slight differences in the time-averaged density profiles in Figure 3 (a),(b). The comparison of Figure 2(b) and 3(c) shows also the role of fluctuations at the end of the emptying process in the Particle Model.

First we give an overview of the various phases exhibited by the system in the MF description, before studying them in more detail in later sections. All figures or analytical calculations presented in the remainder of this paper will refer to the MF model.

3.2 Phase diagram

(a) (b)


Figure 3: Time-averaged density profiles for (a) or (b) obtained from Monte Carlo simulations of the two lane particle model (PM-red/-green) and from the MF equations (MF-black). In (c) a spatio-temporal plot of the MF model is shown.
Parameters: , , and .
Figure 4: Phase diagram for the MF model with and . The system is for either in a pulsing or a low density phase. For the system fills up with plus particles and remains blocked in a state of density one (red line). Purple circles (resp. green square) belong to the low (resp. pulsing) phase. The straight line is a guide to the eyes that separates the different types of symbols and has slope . We checked for , and that this slope does not depend on . For , we found no dependence of the slope for and .
Figure 5: Average density per site , averaged over the whole system and over time, in the MF model. Here , , .

We consider the MF model only in the regime where is at least of order (as for the particle model). In the MF model we find the same three phases identified in Section 2 : a low density phase, a pulsing phase, and a singular line for , where the system fills with plus particles and remains blocked in a state with density one. In the remainder of this paper we only consider the non-trivial case where . In Figure 4 we see that the phase boundary between the low density phase at high /low and the pulsing phase at low /high appears to be a straight line in the plane. For the explored parameter values (see Figure 4), we did not find any dependence of the slope on or . Thus the system’s phase is determined by the ratio between and , at least when they are small compared to .

As an order parameter for the system we take the total density per site,


averaged over the whole system and over time, which we denote . Then as illustrated in Fig 5, the order parameter exhibits a discontinuity at the transition.

Figure 6: Density profile in the low density phase, for , , and with . The inset shows the good agreement between the prediction (17a,17b) and the numerical results.

3.3 Total Current

It is useful to consider the net flux of particles or total current through the link which is given by


As a consequence of mass conservation, we have the continuity equation


In a stationary state, is constant with respect to and equal to zero because of the reflecting boundary of the system:


where the superscript indicates that the variable is taken in the stationary state. Therefore, in a stationary state, at the entrance the net flux must also be zero


This yields the relation


As we always consider that , one finds from (11) that in a stationary state


Note that this is true only for , while can be either small or large (large meaning here close to ), depending on the value of .

4 Analysis of the low density phase

In this section we assume not only that , but also that we are in the low density phase, so that in the stationary state, not only but also is small. To fix ideas, we assume in the remainder of this section that


Assuming smooth variations of the densities in the system, one can derive the partial differential equations describing the behaviour of the system in the continuous space limit. This description is indeed expected to be relevant in the low density limit for which there are no sharp shocks (except near the reflecting boundary).

We introduce a continuous space coordinate and assume that is large. Then, setting and , the Taylor expansion of Eqs (1, 1b) gives


It should be noted that corresponding equations linearised in the densities have been considered in [32] appendix A.

We are interested in the stationary state of (14a, 14b), so that the left-hand side of the equations is zero. For simplicity, we will omit the superscript , indicating the stationary state, in the following.

We also need the Taylor expansion of the zero-flux condition (9) :


We know from (11) and (13) that the density of minus particles on the first site is of order . As in the low density phase there are no steep gradients (apart from near the right boundary), density derivatives will also be at most of this order.

We thus write the densities as


where we assume that and equivalently for . We will then solve equations (14a, 14b 15) order by order (see appendix) resulting in


Figure 6 shows a comparison between direct simulation of the MF equations and the prediction (17a, 17b), with excellent agreement. The lowest order of this calculation recovers that of [32]. However, when considering higher order terms in of the densities, the non-linear terms in (14a, 14b) play a role.

5 Analysis of the pulsing phase: the filling stage

We now turn to the analysis of the behaviour in the pulsing phase which involves filling and emptying stages as part of its cycle. First we show some numerical results and then compare those with analytical expressions. For the analytical expressions we assume that


This scaling ensures that we are in the pulsing phase (see Figure 4).

As can be seen from Figure 3 (c) the filling stage begins with an empty system (, ). Then, as the system rapidly fills, it is divided into two regions: a low density region to the left with , and a high density region to the right with , .

(a) (b)
Figure 7: Density profiles (a) and fluxes (b) for and in the filling stage. Note that actually we plot to allow a direct comparison with . The inset in (b) shows that the fluxes converge towards in the high density region in the filling stage.

In Figure 7 (a) the densities and are shown at a given time during the filling stage of the system. The two regions are separated by a front structure which moves to the left and occupies a few lattice sites over which peaks to one and dips to zero. The peak in corresponds to plus particles that cannot enter the high density region to the right.

One also observes a peak in near the wall, as particles cannot go further to the right. The structure of this boundary layer is actually more complex because and oscillate out of phase near the wall while converging to a value of approximately one half in the bulk. Indeed, in this high density region, fluxes are almost zero. Note, however, that in absolute value, is double , a feature that will be explained in the next section.

Looking at the profile of the flux per bond (Figure 7 (b)) in the low density region, it becomes obvious that the positive flux is (in absolute value) greater than the flux out of the system and so the overall density increases with time until the density is almost one at every site (cf. Figure 3 (b)).

As a first step in understanding the filling stage, we wish to predict the velocity of the moving shock front. To do so, we notice that, in the low density region, the density of minus particles is much smaller than the density of plus particles. For this reason we neglect the outflux of minus particles to predict the filling velocity. If we assume that the density profiles in the low density region are stationary, so that is constant in this region, then the inflow of particles feeds the high density region of density almost , which grows with velocity . Thus we obtain


For the scaling given in (18) we have good agreement between the predicted filling velocity and the numerical results (Figure 8). If increases one observes deviations from this value. A measurable outflux of minus particles (Figure 7 (b)) reduces the filling velocity. However, a numerical study (not presented here) of the density profile in the low density region of the filling phase reveals that it is not stationary in time. Some further effort would be required to take into account this non-stationarity.

Figure 8: dependence of the filling velocity when the system is in a pulsing state (here for , , ). The dashed blue line shows the analytic prediction of eq. (19), valid for small .

6 Analysis of the pulsing phase: the emptying stage

As can be seen from Figure 3 (c) the emptying stage begins with a full lattice (, ). Then, as the lattice empties the system is divided into two regions: a high density region to the left with , and a low density region to the right with , . The two regions are separated by a front structure which occupies a few lattice sites over which peaks to one. The front structure moves to the left with a speed very much less than the speed of the front in the filling stage (see previous section).

In Figure 9 we examine in more detail the density profiles and particle fluxes in the two regions. Figure 9 (a) shows the densities and at a single moment in time. From the values of the densities given in the caption, we already see that the deviations from in the dense region on the left-hand side of the shock and from in the dilute region on the right-hand side are very small.

In Figure 9 (b) as well as , the modulus of the net mass current, are shown. From the numerical data we find two interesting facts: First, in the dense region the flux of minus particles is double the flux of plus particles (in absolute value). We shall explain this fact in the next section. Furthermore, is constant in the dense region whereas it decreases slightly (in absolute value) with increasing in the dilute region.

The first observation implies that is negative, confirming that there is a net mass flow out of the system at the left end. The net mass flow is still negative but small in the low-density right-hand region although the net individual currents are larger than in the dense region. The linear increase in in the low density region corresponds to a linear decrease in in this low density region.

(a) (b)
Figure 9: Density profiles (a) and fluxes (b) in the emptying stage for , and . Note that actually we plot and to allow a direct comparison with . The densities and for this parameter set are typically in the plateau of the dense region, and in the low density phase (see Figs 11 and 13).

6.1 Shock dynamics

We idealise the moving shock seen in Figure 9 (a) as a discontinuity (possibly with a complex shape, as long as it is of finite spatial extent) separating two domains of constant densities. Our aim is to predict the shock velocity through Rankine-Hugoniot equations, which give a relation between the densities and the fluxes on both sides of the shock. In doing so, we implicitly assume that mass is conserved for each species, so that the source terms corresponding to transmutations will not be taken into account.

We denote the fluxes of plus and minus particles on both sides of the shock by


where the index refers to the left- and right-hand side of the shock, respectively. Furthermore, we define the net current on each side of the shock as


In a uniform domain (flat density profiles), the equations (14a,14b) become


Thanks to the source term on the right -hand side of (22), the system relaxes towards a state in which , which is indeed what is observed at lowest order in the numerics. We shall thus assume that in each domain on each side of the shock. We know from the numerics that these values are close to and , respectively. More precisely, using the knowledge gained from the numerical results for the densities (Figure 9(a)) we make the following assumption


with (in Figure 9(a), we find ).

Inserting this in the flux definitions (20) we find


which explains why the flux of minus particles is double that of plus particles to lowest order in .

The Rankine-Hugoniot relations are given by


where we assume that the front for both species moves at the same speed . These equations simply express mass conservation for each species at the shock. Using the expressions (23-24) and observing, furthermore, that and are small compared to , so that we can neglect second order terms in , we obtain


We can now eliminate from the equations (26) by assuming that so as to obtain, to leading order, the relation

(a) (b)
Figure 10: Dependencies of the emptying front velocity . (a) against the deviation from density 1/2 in the dense region and (b) against the transmutation rate . The points on both figures correspond to the same simulations so that the value of for the parameter values used throughout this paper () is . In (a) the line corresponds to the Rankine-Hugoniot expression given in eq. (28) and shows very good agreement. Furthermore, we find that the velocity is also proportional to (the line in (b) is a fit with slope ).
Parameters: , .

and an expression for the front velocity


In Figure 10 we test the predicted relation (28) numerically and find that indeed and are proportional to each other, with a proportionality coefficient well predicted by the Rankine-Hugoniot prediction (28). Furthermore, we checked that does not depend on the system size (see for example Figure 11 (c)). However, we do not have a prediction for the value of . In Figure 10 (b), we observe that —and thus also —is proportional to .

In the Rankine-Hugoniot calculation, we have assumed constant domains on both sides of the shock. We shall see now that this is not fully accurate, and that a more careful study reveals some oscillations, both in the dense and in the dilute region.

6.2 Dense region

In the dense region, to leading order, the density of both species are equal and far away from the front are approximately equal to one half.

6.2.1 Oscillations in the dense region

Inspecting Figure 11(a), we observe that in the bulk of the dense region, both densities reach a plateau value . On both sides of the dense region, we observe deviations from these plateau values. The semi-logarithmic plots and in Figure 11(b) show that at the entrance, the deviations from the plateau densities decay exponentially when the particles penetrate the system. On the right -hand side of the dense region, we find some spatial oscillations in the density profiles that are also damped exponentially in the bulk.




Figure 11: Density profiles and in the dense region of the emptying stage. (a) Zoom around density on the density profiles and , for and . Only the sites near the boundaries of the dense region are shown (with a translation of the right sites in the case , in order to overlap with the data). One sees the decay of density perturbations to the bulk values at the entrance and the oscillatory behaviour on the right of the figure (left-hand side of the shock position ). (b) and for , the straight lines show the results obtained by the stability analysis of section 7. (c) for (black) and (red) (the data points are translated horizontally in order to have the same position as for ). Parameters: , , . For this parameter set, we use the same for the two different lattice lengths.

In the simulations (see Figure 11(a)), the oscillations in the dense region to the left-hand side of the shock have a wavelength of about seven sites both for and . Furthermore, we observe that and oscillate in opposite phase, with an almost perfect symmetry. Indeed, and overlap up to at least . Again, this holds both for and .

In Figure 11(c) we superimpose the profiles for two different lattice lengths and . We can see that the oscillations on the left-hand side of the shock appear to have the same exponential decay for the two lattice lengths, the same wavelength and even the same structure. A zoom in Fig 11(a) (not shown) shows also that there is no significant difference between the two systems sizes either for the wavelength or for the amplitude of the oscillations. Thus, the oscillatory behaviour appears to be independent of system size. The exponential decay at the entrance of the dense region is also independent of system size.

Measurements of the slopes of for and in Figure 11(c) give at the entrance a decay constant and close to the shock a decay constant .

We shall see in section 7 how a stability analysis can explain these observations.

Figure 12: Illustration of the shock moving by one site. Given are the density profiles for (empty symbols) and (filled symbols) at different times. We can observe how the density profiles are deformed when the shock moves one step to the left. The profile at is the same as at translated by one lattice unit to the left. Here, , , , and we find time units. The time interval between two successive profiles is .

6.3 Dilute region

In the emptying stage of the pulsing phase there is a dilute region to the right of the moving shock (see Figures 9). The analysis of Section 6.1 assumes equal low densities to the right of the shock. However, in this section we carry out a more detailed numerical study that reveals that, though they are on average of the same order, the densities and actually differ and oscillate in time.

Figure 13 illustrates the temporal variations of both densities at one given lattice site  located in the dilute region. We observe fast oscillations in and as a function of time. Neither nor change in a sinusoidal way, rather they show highly non-linear behaviour.

For the parameter set used in Figure 13, the typical value around which the densities oscillate is , a value which differs from the value predicted by the Rankine-Hugoniot relation (27), . This is not surprising, as obviously, the assumption that and are equal and uniform is too naive. Consideration of the source terms (which have been ignored in the analysis of Section 6.1) and the complex structure of the density profiles would be crucial to find a solution in this region. It seems out of reach to derive the profiles analytically, but several interesting observations can be made.

The period of the oscillations in Figure 13 is about 6 time units, which corresponds to the time period for the shock to move from one lattice site to the next one, as shown on Figure 12. Thus, this oscillatory behaviour seems directly related to the hopping of the shock on the discrete lattice. We checked that the behaviour was robust to changes in the time-step of the numerical iteration scheme used to solve the MF equations.

In Figure 14 we plot the density profiles in the dilute region at closely separated time intervals. One sees that the shock emits microscopic density pulses which correspond to the density peaks at a given site seen in Figure 13. The steepness of the peaks is to be related to the rather abrupt way in which the shock moves from one site to the next (see Figure 12).

Figure 13: Temporal oscillations of (red) and (green) at a fixed lattice site in the dilute region of the emptying state. Both densities exhibit non-sinusoidal oscillations.
Parameters: , , .
Figure 14: The two densities in the dilute region of the emptying stage, (a) and (b) at various equidistant times. The time interval shown time units is less than a period in Figure 13, or equivalently, less than the time interval shown in Figure 12. One observes, close to the shock, a microscopic pulse propagating to the right with a decaying amplitude as it propagates. One can also guess the presence of another microscopic pulse around previously emitted and already strongly damped. The parameters values are , , .

7 Stability analysis

In this section we consider perturbations about constant density solutions of the MF equations within a linear stability analysis. In many works [19, 27] it has been shown that a linear analysis of MF equations of the form of eq. (30) can give insight into the full solutions.

7.1 Perturbation about constant profiles

The discrete bulk equations (1-1b) accept as a stationary solution any constant profile provided that


where the last inequality comes from the fact that the total density cannot be larger than one.

We now wish to perturb this solution and to explore its linear stability. We denote the time-dependent perturbation about a constant profile with


On linearising in and , the bulk equations (1,1b) become


where it should be noticed that a factor could be factorised in the flux terms. These equations may be written in matrix form as


where the two-by-two matrix is defined by


The time dependence of the amplitudes , is thus determined by the eigenvalues of the matrix which are functions of and .

The equation for the eigenvalues of is quadratic in and is fourth order in , and it is in general difficult to write the solutions in a simple form. The solutions for provide information about the temporal stability while the dependence governs the behaviour of spatial perturbations. In the following we will consider two simple relevant cases separately, first by assuming a homogeneous perturbation in space () and second by assuming a stationary but spatially varying perturbation ().

7.2 Temporal stability of plateau solutions

First we analyse whether plateau (constant density) solutions are dynamically stable. For this reason we choose which gives us an easy equation for the eigenvalues namely


with the obvious roots


The corresponding eigenvectors are given by


The first eigenvector corresponds to a perturbation consisting of a uniform and equal shift of the two densities. The fact that confirms that any constant profile with is a stationary state.

The second eigenvector gives the same shift to both densities but with opposite sign. In this case the perturbation decreases exponentially with time, and both densities converge towards the constant value . The system is thus dynamically stable with respect to this type of perturbation.

For a combination of these shifts, densities will converge, but towards and not towards . From this we can conclude that the system is marginally stable with respect to uniform shifts in the densities.

7.3 Spatial stability

We now seek spatial perturbations around a constant profile that are stationary. They will be given by solutions for of the characteristic equation det M  for (stationary solutions). It is of no surprise that one solution is , as was already found in the previous subsection. But here we are not interested in this solution, which leads to uniform perturbations.

In order to find the other solutions, one needs to solve a cubic equation


The solutions obviously depend on the bulk plateau density . We find that for a low density , the three roots are real. Interestingly, above a critical density two of these real roots become a complex conjugate pair.

For the plateau density of Figure 11, namely we find that the real root has a modulus smaller than . This corresponds to perturbations that are exponentially damped for increasing lattice index. Such a perturbation is indeed observed on the numerical density profiles of Figure 11, near the entrance of the system. The damping rate predicted by the linear stability analysis (slope of the straight line in Figure 11 (b)) is in perfect agreement with the numerics.

The two other complex conjugate roots and correspond to oscillatory density perturbations. We find that above a density such that , the modulus of and is larger than . This explains that the oscillatory perturbations observed for on the right of the dense region (see Figure 11 (a)) are damped when the lattice index decreases. The modulus gives the damping of the oscillatory perturbation when it penetrates into the bulk of the dense region. Although not perfect, there is still a good agreement with the damping observed in the numerics (see Figure 11 (b)).

Besides, the wavelength of the oscillations predicted by the linear analysis for the parameters of Figure 11 is , in good agreement with the numerical result of Figure 11 (a).

In conclusion, although still has to be determined from the numerics, the linear analysis explains to a large extent the structure of the density profiles in the dense region. Again, however, the shock is moving so the assumption of stationarity of profiles is not fully accurate. Instead there is non-stationarity orginating from the motion of the shock between the discrete lattice sites. It would be of interest to understand better how the propagation of spatial oscillations in the profile and the dynamics of the shock front are coupled.

8 Transition between the low density phase and the pulsing phase.

In this section we try to understand how the transition between the pulsing and low density phase is characterised. We saw in Figure 5 that the transition is discontinuous with respect to an order parameter defined by the average density through the system. In Figure 15 we show the density profiles for values of close to the phase transition. One observes that with increasing the single peak in close to the reflecting boundary detaches from the boundary and at the same time a peak in arises. Once both densities approach the value one half at a certain site of the lattice, the transition from the low density phase to the pulsing phase occurs.

We can understand the phase transition better from the particle picture. In the low density phase a boundary layer of plus particles exists at the right reflecting boundary. After some time these plus particles transmute and can then leave the system again. If is low these minus particles have enough time to leave the dense region at the reflecting boundary. On increasing , more and more plus particles arrive before those in the boundary layer transmute into minus particles and leave. In that case those minus particles close to the boundary cannot leave the system anymore. The density in both, plus and minus particles, increases and approaches 1/2. This mutual blocking is necessary for the system to fill up and then trigger the pulsing behaviour.

Figure 15: The two densities at the reflecting boundary close to the phase transition. We show and for different values of the injection rate as the discontinuous transition is approached at .
Parameter: , , .

9 Conclusion

In this work we have presented two models which exhibit a spontaneous pulsing behaviour. At first we considered a particle model inspired by [32] where two types of particles hop along two lattice lanes of length . At the left boundary, plus particles are injected and minus particles leave while the right boundary is a reflecting one. The two types of particles differ in their hopping direction but also in their lane changing behaviour. Plus particles cannot change lanes whereas minus particles can if the next site in the same lane is already occupied by another particle. Along the whole system particles can transmute and convert their type. For large enough so that particles invade the system, the phase diagram is divided into two phases, a low density phase and a pulsing phase where the system shows spontaneous pulsing in the sense that it fills up until the density is almost one and then empties again.

Inspired by this explicit particle model we introduced a discrete MF model. It is defined by the master equations (13b) in the particle densities (for plus particles) and (for minus particles). We find that it also exhibits either a very regular pulsing behaviour or a low density phase (excluding the trivial case where particles do not change their type). We start our analysis in this low density phase, for which we can predict the stationary density profile.

The pulsing phase exhibits a more complex spatial structure and temporal behaviour, which makes the analysis more involved. In the filling stage we can predict analytically the filling velocity for a small transmutation rate. In the emptying stage the system’s dynamics is very slow. A front moves to the left, separating a dense region on the left and a dilute region on the right. In the dense region (between the entrance and the shock position), both densities deviate from only by very small values, typically or smaller. In a first approximation we neglect transmutation and predict the shock’s velocity as a function of this deviation, by a mass conservation argument. The agreement with the numerics is very good (Figure 10). In the dense region, the density for both types of particles is mostly constant in the bulk. However, some deviations from the plateau densities are observed at both ends of the dense region. Near the entrance, these deviations decay purely exponentially, while near the shock the perturbations additionally exhibit some oscillations in the dense region.

A linear analysis in the dense region allows us to predict the nature of the perturbations, the damping rate for each of them, and the wavelength of the oscillatory ones, with a good agreement with the numerics.

In the dilute region, one observes that some microscopic density pulses are released from the shock and propagate through the system to the right boundary. This microscopic pulse emission is connected to the hopping of the shock from one lattice site to the next one.

Eventually, the transition from the low density phase to the pulsing phase appears to occur through the detachment of an accumulation of particles from the reflecting boundary. The resulting transition is discontinuous.

For this system, several challenging questions are still open. Though we were able to characterise the low density phase, and to explain the structure of density profiles in the pulsing phase, in the emptying stage of the pulsing phase we are still lacking a full prediction of the shock velocity—and hence of the pulsing period. The fact that each particle type density is not conserved locally adds a further complication. It would be of interest to be able to take it into account in some generalised mass conservation relations. The presence of steep gradients separating constant regions makes it difficult to use perturbative calculations.

We have also seen how the motion of the shock between the discrete lattice sites gives rise to microscopic density pulse emissions which in turn imply non-stationary profiles. There are many systems, for example involving molecular motor dynamics on filaments, where discreteness of a spatial lattice is imposed. Therefore would be of interest to further explore the mechanisms of such microscopic pulses.

We thank Danielle Hilhorst, Frédéric Lagoutière, Christian Tenaud, Henk Hilhorst and Pascal Viot for fruitful discussions. This work was supported by “Investissements d’Avenir” of LabEx PALM (ANR-10-LABX-0039-PALM) – Funding agency : Investissement d’Avenir LabEx PALM – Grant number : ANR-10-LABX-0039 S.K. was supported by the Deutsche Forschungsgemeinschaft (DFG) within the collaborative research center SFB 1027 and the research training group GRK 1276. M.R.E. acknowledges partial support from EPSRC under grant number EP/J007404/1 and thanks LPT for hospitality.

Appendix A Order by order expansion in the low density phase

The lowest order of the expansion (15) yields


while the lowest order of Eqs (14a,14b) gives (in the stationary state)


As a result, one finds that and . Then we get the constant value of from the entrance boundary


We can show that the next order brings no new term, so that


We are thus left with


To next order, assuming that and , the no flux condition (15) gives


and the bulk equations yield


The sum of these last equations gives


where the last equality is obtained by using (45) into (47).

Now combining (11) and (45) taken for , we get


Collecting all terms, we have eventually



  • [1] D. Mukamel. Phase transitions in nonequilibrium systems. In Soft and Fragile Matter: Nonequilibrium Dynamics, Metastability and Flow. IoP publishing, Bristol, 2000.
  • [2] J. Krug. Boundary-induced phase transitions in driven diffusive systems. Phys. Rev. Lett., 67:1882, 1991.
  • [3] C. MacDonald, J. Gibbs, and A. Pipkin. Kinetics of biopolymerization on nucleic acid templates. Biopolymers, 6:1–25, 1968.
  • [4] T. Chou, K. Mallick, and R. K. P. Zia. Non-equilibrium statistical mechanics: from a paradigmatic model to biological transport. Reports on progress in physics, 74:116601, 2011.
  • [5] Y. Aghababaie, G.I. Menon, and M. Plischke. Universal properties of interacting Brownian motors. Phys. Rev. E, 59:2578 – 2586, 1999.
  • [6] C. Appert-Rolland, M. Ebbinghaus, and L. Santen. Intracellular transport driven by cytoskeletal motors : General mechanisms and defects. Phys. Rep., 593:1–59, 2015.
  • [7] F. Peruani, J. Starruss, V. Jakovljevic, L. Søgaard-Andersen, A. Deutsch, and M. Bär. Collective motion and nonequilibrium cluster formation in colonies of gliding bacteria. Phys. Rev. Lett., 108:098102, 2012.
  • [8] M.E. Cates. Diffusive transport without detailed balance in motile bacteria: does microbiology need statistical physics? Rep. Prog. Phys., 75:042601, 2012.
  • [9] J. Stenhammar, A. Tiribocchi, R.J. Allen, D. Marenduzzo, and M.E. Cates. Continuum theory of phase separation kinetics for active brownian particles. Phys. Rev. Lett., 111:145702, 2013.
  • [10] J. Toner and Y. Tu. Flocks, herds, and schools: A quantitative theory of flocking. Phys. Rev. E, 58:4828, 1998.
  • [11] S. Ramaswamy. The mechanics and statistics of active matter. Annu. Rev. Condens. Matter Phys., 1:323–345, 2010.
  • [12] A. Czirók, A.-L. Barabási, and T. Vicsek. Collective motion of self-propelled particles: Kinetic phase transition in one dimension. Phys. Rev. Lett., 82:209, 1999.
  • [13] O.J. O’Loan and M.R. Evans. Alternating steady state in one-dimensional flocking. J. Phys. A: Math. Gen., 32:L99, 1999.
  • [14] A.P. Solon and J. Tailleur. Revisiting the flocking transition using active spins. Phys. Rev. Lett., 111:078101, 2013.
  • [15] D.S. Calovi, U. Lopez, S. Ngo, C. Sire, H. Chaté, and G. Theraulaz. Swarming, schooling, milling: phase diagram of a data-driven fish school model. New Journal of Physics, 16:015026, 2014.
  • [16] M. Howard, A.D. Rutenberg, and S. de Vet. Dynamic compartmentalization of bacteria: Accurate division in E. Coli. Phys. Rev. Lett., 87:278102, 2001.
  • [17] J.S. Berg and R.E. Cheney. Myosin-X is an unconventional myosin that undergoes intrafilopodial motility. Nat. Cell Biol., 4:246–250, 2002.
  • [18] C. Schiffmann, C. Appert-Rolland, and L. Santen. Shock dynamics of two-lane driven lattice gases. J. Stat. Mech., page P06002, 2010.
  • [19] M.R. Evans, Y. Kafri, K.E.P. Sugden, and J. Tailleur. Phase diagrams of two-lane driven diffusive systems. J. Stat. Mech., page P06009, 2011.
  • [20] A. Parmeggiani, T. Franosch, and E. Frey. Totally asymmetric simple exclusion process with Langmuir kinetics. Phys. Rev. E, 70:046101, 2004.
  • [21] S. Klumpp and R. Lipowsky. Traffic of molecular motors through tube-like compartments. J. Stat. Phys., 113:233–268, 2003.
  • [22] K. Nishinari, Y. Okada, A. Schadschneider, and D. Chowdhury. Intracellular transport of single-headed molecular motors KIF1A. Phys. Rev. Lett., 95:118101, 2005.
  • [23] I. Pinkoviezky and N.S. Gov. Modelling interacting molecular motors with an internal degree of freedom. New J. Phys., 15:025009, 2013.
  • [24] I. Pinkoviezky and N.S. Gov. Transport dynamics of molecular motors that switch between an active and inactive state. Phys. Rev. E, 88:022714, 2013.
  • [25] I. Pinkoviezky and N. S. Gov. Traffic jams and shocks of molecular motors inside cellular protrusions. Phys. Rev. E, 89:052703, 2014.
  • [26] T. Reichenbach, T. Franosch, and E. Frey. Exclusion processes with internal states. Phys. Rev. Lett., 97:050603, 2006.
  • [27] A.I. Curatolo, M.R. Evans, Y. Kafri, and J. Tailleur. Multilane driven diffusive systems. J. Phys. A: Math. Theor., 49:095601, 2016.
  • [28] Y. Lan and G.A. Papoian. The stochastic dynamics of filopodial growth. Biophys J., 94:3839–3852, 2008.
  • [29] K. Wolff, C. Barrett-Freeman, M.R. Evans, A.B. Goryachev, and D. Marenduzzo. Modelling the effect of myosin X motors on filopodia growth. Phys. Biol., 11:016005, 2014.
  • [30] J. Tailleur, M.R. Evans, and Y. Kafri. Nonequilibrium phase transitions in the extraction of membrane tubes by molecular motors. Phys. Rev. Lett., 102:118109, 2009.
  • [31] K.E.P. Sugden, M.R. Evans, W.C.K. Poon, and N.D. Read. Model of hyphal tip growth involving microtubule-based transport. Phys. Rev. E, 75:031909, 2007.
  • [32] C. Lin, G. Steinberg, and P. Ashwin. Bidirectional transport and pulsing states in a multi-lane ASEP model. J. Stat. Mech., page P09027, 2011.
  • [33] P. Ashwin, C. Lin, and G. Steinberg. Queueing induced by bidirectional motor motion near the end of a microtubule. Phys. Rev. E, 82:051907, 2010.
  • [34] M. Schuster, S. Kilaru, P. Ashwin, C. Lin, N.J. Severs, and G. Steinberg. Controlled and stochastic retention concentrates dynein at microtubule ends to keep endosomes on track. EMBO Journal, 30:652–664, 2011.
  • [35] R. Juhász. Dynamics at barriers in bidirectional two-lane exclusion processes. J. Stat. Mech., page P03019, 2010.
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
Add comment
Loading ...
This is a comment super asjknd jkasnjk adsnkj
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test description