The asymmetric simple exclusion process on chains with a shortcut revisited
aa r X i v : . [ phy s i c s . b i o - ph ] N ov The asymmetric simple exclusion process on chains with ashortcut revisited
Nadezhda Bunzarova , , Nina Pesheva , and Jordan Brankov , Bogoliubov Laboratory of Theoretical Physics,Joint Institute for Nuclear Research, 141980 Dubna, Russian Federation Institute of Mechanics, Bulgarian Academy of Sciences, 1113 Sofia, Bulgaria
We consider the asymmetric simple exclusion process (TASEP) on open networkconsisting of three consecutively coupled macroscopic chain segments with a shortcutbetween the tail of the first segment and the head of the third one. The model wasintroduced by Y.-M. Yuan et al [J. Phys. A 40, 12351 (2007)] to describe directedmotion of molecular motors along filaments. We report here unexpected results inthe case of maximum current through the network which revise the previous findings.Our theoretical analysis, based on the effective rates approximation, shows that thesecond (shunted) segment can exist in both low-density and high-density phases,as well as in the coexistence (shock) phase. Numerical simulations demonstratethat it is the last option that takes place - the local density distribution and thenearest-neighbor correlations in the middle chain correspond to a shock phase withcompletely delocalized domain wall. Surprisingly, the main quantitative parametersof that shock phase are governed by a real root of a cubic equation the coefficients ofwhich simply depend on the probability of choosing the shortcut. The unexpectedconclusion is that a shortcut in the bulk of a single lane always creates traffic jams.
Pacs:
Keywords : TASEP, traffic flow models, non-equilibrium phase transitions, traffic oncomplex networks, biological transport
Introduction. — The asymmetric simple exclusion process (TASEP) is one of the paradig-matic models for understanding the rich world of non-equilibrium phenomena. Devised tomodel kinetics of protein synthesis [1], it has found a number of applications to vehiculartraffic flow [2], biological transport [3], one-dimensional surface growth [4], forced motion ofcolloids in narrow channels [5], spintronics [7], transport of ’data packets’ in the Internet [6],current through chains of quantum dots [8], to mention some.Novel features of the TASEP have been found on networks consisting of coupled linearchains with nontrivial geometry. In the approach advanced in our work [9] each macroscopicsegment s of the network is considered in a stationary phase determined by its effectiveinjection α ∗ s and ejection β ∗ s rates. At that, exact in the thermodynamic limit results for thedensity profile are incorporated. The only molecular field type approximation used consistsin the neglect of correlations between different chain segments. This allows one to treatthe coupling between each two connected segments as coupling to reservoirs with certaineffective rates. The possible phase structures of the whole network are obtained as solutionsof the resulting set of equations for the effective rates that follow from continuity of current.The importance of our approach for modeling complex biological transport phenomena waspointed out by Pronina and Kolomeisky [10]. This method became very popular and wasused in a number of studies of TASEP and its generalizations on networks with differentgeometries, e.g., with junctions, bifurcations, intersections, interacting lanes [11]. Finite-sizeeffects on the density profile due to shifting the position of the double-chain section fromthe middle of the linear network were studied too [12].Here we consider the TASEP on open chain with a shortcut in the bulk, introducedas ’model A’ in [13]. The current through the shortcut is proportional to a probability q . It is convenient to consider the system as composed of three consecutively connectedmacroscopic chain segments and a shortcut between the tail of the first segment and thehead of the third one. In principle, the effect of a shortcut can easily be understood: thedecrease in the current through the shunted part (second segment) of the original chainleads to a sharp change of the particle density in the latter. If the chain without a shortcut( q = 0) is in the low-density (LD) phase, its bulk density ρ LDbulk < / J = ρ LDbulk (1 − ρ LDbulk ) < /
4. The shortcut takes a part J sc > J (2) = J − J sc has to be supported by still less bulkdensity ρ (2)bulk < ρ LDbulk in that segment. Similarly, when the initial chain is in the high-density(HD) phase with ρ HDbulk > /
2, the drop in the current through the second segment, causedby the shortcut, leads to a still higher bulk density in that segment, ρ (2)bulk > ρ HDbulk . In thesecases all the three segments remain in the same phase, though with different density of themiddle one.Not so clear, however, is the situation when the initial chain is in the maximum current(MC) phase with ρ MCbulk = 1 /
2. Now the drop in the current through the shunted segmentof the network can be compensated equally well by decrease or increase in its bulk density.Then, the middle segment is forced either in low-density, or in high-density phase, whichmay lead also to coexistence of LD phase on the left-hand side with HD phase on the right-hand side. This phase structure is additionally favored by the downward (upward) bend inthe density profile of the first (third) segment in the maximum current phase. In the caseof open system with variable total number of particles the coexisting phases are likely to beseparated by a completely delocalized domain wall. Such was the situation observed in eachof the equivalent segments in a double-chain section incorporated in the middle of a longlinear chain [9]. It seems plausible that the above mechanism of influence of the shortcuton the phase state of the shunted segment should be invariant with respect to the explicitstructure of the shortcut. In particular, one may consider a shortcut in the form of anadditional (shorter) chain connecting the last site of the first segment to the first site of thethird one. Since the length of the shortcut is irrelevant, we can include the case of parallelsegments with equal length considered in our work [9]. However, the authors of [13] haveclaimed that in the case of their ’model A’, the shunted middle segment can exist only in thehigh-density phase. This contrast in the conclusions motivated us to renew the study, bothanalytically and numerically, of the model. The results may have important implications forvehicular traffic flow control, as well as for biological transport in living cells.
Microscopic model. — Here we consider model A of a shortcut, suggested in [13], whenboth the injection α and ejection β rates at the open ends of the system are larger than1/2, so that the first and third segments are in the maximum current phase. The shortcut isbetween the last site of the first segment, with occupation number τ (1) L , and the first site of thethird segment, with occupation number τ (3)1 . According the rules of the random-sequentialalgorithm, when a particle at the last site of the first segment (with τ (1) L = 1) attempts tomove, the particle may jump along the main track to the first site of the second segmentwith rate (1 − q )(1 − τ (2)1 )(1 − τ (3)1 ) + (1 − τ (2)1 ) τ (3)1 , or take the shortcut to the first site of thethird segment with rate q (1 − τ (3)1 ), or stay immobile with rate (1 − q ) τ (2)1 (1 − τ (3)1 ) + τ (2)1 τ (3)1 .These rules lead to exact expressions for the stationary current through segment 2, J (2) = (1 − q ) h τ (1) L (1 − τ (2)1 )(1 − τ (3)1 ) i + h τ (1) L (1 − τ (2)1 ) τ (3)1 i = h τ (2) L (1 − τ (3)1 ) i , (1)and through the shortcut J sc = q h τ (1) L (1 − τ (3)1 ) i , 0 ≤ q ≤ Theoretical analysis. — In the effective rates analysis [9] one neglects the correlationsbetween sites belonging to different segments, so that the above expressions simplify to J (2) = ρ (1) L (1 − ρ (2)1 )[(1 − q )(1 − ρ (3)1 ) + ρ (3)1 ] = ρ (2) L (1 − ρ (3)1 ) , (2)and J sc = qρ (1) L (1 − ρ (3)1 ), where ρ (s) i = h τ (s) i i , s = 1 , ,
3, is the average value of the occupationnumber τ (s) i in a given stationary state. Within the above approximation effective injection, α ∗ s , and ejection, β ∗ s , rates for segment s = 1 , , J (s) = β ∗ s ρ (s) L = α ∗ s (1 − ρ (s)1 ), with α ∗ = α and β ∗ = β . Thus, taking into account that J (1) = J (3) = J (2) + J sc , one obtains α ∗ = α, β ∗ = 1 − ρ (2)1 + qρ (2)1 (1 − ρ (3)1 ) , (3) α ∗ = ρ (1) L [1 − q (1 − ρ (3)1 )] , β ∗ = 1 − ρ (3)1 , (4) α ∗ = ρ (2) L + qρ (1) L , β ∗ = β. (5)Here we have taken into account that α ∗ comes from both the expression for the current J sc and the last one in (2). Expressions (3)-(5) for the effective rates coincide exactly withequations (4) obtained in [13]. However, the results of our analysis in the case when the firstand third segments are in the maximum current phase are essentially different from thoseclaimed in [13].We confine ourselves to the study of possible phase structures of the type ( M, X, M ),when the first and third segments are in the maximum current phase M , and the secondsegment is in a low-density phase ( X = L ), high-density one ( X = H ), or on the coexistenceline ( X = C ). Note that the case ( X = M ) is excluded, since the presence of a shortcut( J sc >
0) implies J (2) < J (1) = J (3) = 1 /
4. To check the consistency of a given structure(
M, X, M ) with the corresponding conditions on the effective rates (3)-(5), we make useof the known, exact in the thermodynamic limit, values of the bulk density ρ (s)bulk and localdensities ρ (s)1 , ρ (s) L , in dependance on the thermodynamic phase of each segment s = 1 , , L ≫ ρ (1)bulk = 1 / , ρ (1)1 = 1 − / (4 α ) , ρ (1) L = 1 / (4 β ∗ ) (6) ρ (3)bulk = 1 / , ρ (3)1 = 1 − / (4 α ∗ ) , ρ (3) L = 1 / (4 β ) . (7)By inserting the expressions for ρ (1) L and ρ (3)1 into Eq. (4), we obtain α ∗ = 1 / (4 β ∗ ) − q/ (16 α ∗ β ∗ ) , β ∗ = 1 / (4 α ∗ ) , (8)and, from Eq. (2), J sc = 1 / − J (2) = q/ (16 β ∗ α ∗ ). Now we pass to the separate considerationof each of the possibilities X = L, H, C . Middle segment in the low-density phase. — In this case ρ (2)bulk = ρ (2)1 = α ∗ , J (2) = α ∗ (1 − α ∗ ) , ρ (2) L = α ∗ (1 − α ∗ ) /β ∗ . (9)Substituting the expressions for ρ (2)1 and ρ (2) L into Eqs. (3) and (5), we find β ∗ = 1 − α ∗ [1 − q/ (4 α ∗ )] , α ∗ = α ∗ (1 − α ∗ ) /β ∗ + q/ (4 β ∗ ) . (10)We have obtained a set of four nonlinear equations, see (8) and (10), for the four effectiverates β ∗ , α ∗ , β ∗ , and α ∗ . The solution depends on one free parameter, because one of theseequations is a consequence of the other three. From the first equation in (8), we obtainthe equation α ∗ = 1 / (4 β ∗ ) − J sc = 1 / (4 β ∗ ) − / α ∗ (1 − α ∗ ), where we have used therelationship J sc = 1 / − J (2) in combination with the expression for J (2) given in (9). Hence,we find β ∗ = 1 / [1 + 4( α ∗ ) ]. After substitution of this expression into the first equation (10),we solve the latter for α ∗ and obtain α ∗ = q [1 + 4( α ∗ ) ] / [4(1 − α ∗ ) ]. Finally, the secondequation in (8) yields β ∗ = (1 − α ∗ ) / { q [1 + 4( α ∗ ) ] } . Now one can readily verify that theabove expressions for β ∗ , α ∗ and β ∗ identically satisfy the second equation in (10).It remains to check the consistence of the results obtained with the conditions for( M, L, M ) phase structure of the network. The free parameter α ∗ has to satisfy the in-equality α ∗ < / β ∗ > / α > / α ∗ < β ∗ for the secondsegment to be in the low-density phase leads to the cubic inequality4 q ( α ∗ ) − α ∗ ) + (4 + q ) α ∗ − < . (11)This inequality has to be fulfilled simultaneously with the condition α ∗ > / β > / − α ∗ ) < ( q/ α ∗ ) ] . Therefore, the free parameter α ∗ has to obey the constraints qα ∗ [1 + 4( α ∗ ) ] < (1 − α ∗ ) < ( q/ α ∗ ) ] (12)which define a nonempty interval when α ∗ < /
2. As a simple consequence, in the case ofvanishing probability of the shortcut, q → − , the free parameter α ∗ → / − , which agreeswith the result α ∗ = ρ bulk = 1 / Middle segment in the high-density phase. — In this case the exact thermodynamicparameters of the second segment are ρ (2)bulk = ρ (2) L = 1 − β ∗ , J (2) = β ∗ (1 − β ∗ ) , ρ (2)1 = 1 − β ∗ (1 − β ∗ ) /α ∗ . (13)Substituting the above expressions for ρ (2)1 and ρ (2) L into Eqs. (3) and (5), we find β ∗ = q/ (4 α ∗ ) + [1 − q/ (4 α ∗ )] β ∗ (1 − β ∗ ) /α ∗ , α ∗ = 1 − β ∗ + q/ (4 β ∗ ) . (14)Taking into account Eqs. (8), we have again a set of four nonlinear equations for the foureffective rates. We shall solve these equations and show that one of them is a consequenceof the other three. As in the previous case, one of the effective rates will appear as a freeparameter in the solution. The second equation in (8) yields α ∗ = 1 / (4 β ∗ ). Combining thisresult with the second equation in (14), we express β ∗ as a function of β ∗ : β ∗ = qβ ∗ / [1 − β ∗ (1 − β ∗ )]. Taking into account that the current trough the shortcut is J sc = 1 / − J (2) ,in view of the present expression for the current J (2) , see Eq. (13), we can rewrite the firstequation in (8) as α ∗ = 1 / (4 β ∗ ) − / β ∗ (1 − β ∗ ). The substitution here of β ∗ yields α ∗ = [1 / ( qβ ∗ ) −
1] [1 / − β ∗ (1 − β ∗ )]. One can readily check that the above expressions for α ∗ , α ∗ , and β ∗ satisfy the first equation in (14) identically with respect to β ∗ = 1 − ρ (2) L .Finally, we check the conditions on the effective rates which imply the phase structure( M, H, M ). The condition α ∗ > β ∗ leads to the cubic inequality4 q ( β ∗ ) − β ∗ ) + (4 + q ) β ∗ − > , (15)which, together with β ∗ < / α > / β ∗ > /
2. The secondcondition leads to the inequality qβ ∗ > (1 / − β ∗ (1 − β ∗ )]. The right-hand side being -0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,10,200,250,300,350,400,450,500,55-0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,10,200,250,300,350,400,450,500,55 α * 2 = ρ (2)1 q FIG. 1. (Color online) Comparison of the numerically evaluated a ∗ = ρ (2)1 , shown by red stars,with the values of the appropriate root of the cubic equation (18), shown by blue circles, for tendifferent values of the rate q . non-negative, q → β ∗ → /
2, hence J (2) → / J sc →
0. The analysis of thisquadratic inequality is equivalent to β ∗ < / q/ − (1 / p q (4 + q ), the right-hand side ofwhich is less than 1/2. The condition α ∗ > / β > /
2) is satisfied whenever β ∗ < / Middle segment on the coexistence line. — The second segment can exist in a low- orhigh-density phase, depending on whether α ∗ > β ∗ or α ∗ < β ∗ , respectively. Naturally, weexpect the coexistence phase (shock phase) to take place at a common point in the closureof the above open sets, i.e., when the rates α ∗ = β ∗ coincide with an appropriate root ofthe cubic equation given by an equality sign in expressions (15) and (11). To prove this, weset α ∗ = β ∗ and assume the exact in the thermodynamic limit values of the current and thelocal densities at the endpoints of the second segment in the coexistence phase, J (2) = α ∗ (1 − α ∗ ) = β ∗ (1 − β ∗ ) , ρ (2)1 = α ∗ = ρ − ( J (2) ) , ρ (2) L = 1 − α ∗ = ρ + ( J (2) ) . (16)Here ρ ± ( J ) = (1 ± √ − J ) / J . Substituting the above expressions for ρ (2)1 and ρ (2) L into Eqs. (3) and (5), we obtain β ∗ = ρ + ( J (2) ) + qρ + ( J (2) ) / (4 α ∗ ) , α ∗ = ρ + ( J (2) ) + q/ (4 β ∗ ) . (17) x(x) FIG. 2. (Color online) Local density profile at α = β = 0 .
75 and q = 0 . x = i/L , where i = ( s − L + 1 , ( s − L + 2 , . . . , sL labels the sites in the segment s , s = 1 , ,
3. The shape of the density profile in the first and thirdsegments is typical for the MC phase, while that in the second segment closely resembles the lineardependence with the distance characteristic of the coexistence phase with completely delocalizeddomain wall. The predictions of the domain wall theory are shown by a blue line.
Inserting in the first equation 1 / (4 α ∗ ) expressed from the second equation in (8) with β ∗ = α ∗ , and replacing ρ + ( J (2) ) and ρ − ( J (2) ) by 1 − α ∗ and α ∗ , respectively, we obtain β ∗ = 1 − α ∗ + q ( α ∗ ) . On the other hand, dividing both sides of the first equation in (17) by4 β ∗ , and using J sc = 1 / − J (2) , we arrive at 1 / ρ + ( J (2) ) / (4 β ∗ ) + ρ − ( J (2) ) (cid:0) / − J (2) (cid:1) .Solving the above for β ∗ , and using that J (2) = ρ + ( J (2) ) ρ − ( J (2) ) = ρ + ( J (2) ) α ∗ , we obtain β ∗ = [1 + 4( α ∗ ) ] − . Clearly, β ∗ > / α ∗ < /
2. The equality the right-hand sides ofthe two above derived expressions for β ∗ leads to the equation4 q ( α ∗ ) − α ∗ ) + (4 + q ) α ∗ − . (18)Unexpectedly, the value of α ∗ = β ∗ is determined by a root of cubic equation, which issingular at q = 0. More precisely, these effective rates are given as function of q , 0 ≤ q ≤ q → + .A comparison of the values of α ∗ given by the appropriate root of (18) and ρ (2)1 evaluatedby computer simulations are shown in Fig. 1. From the second equation in (8) at β ∗ = α ∗ q=0.9 q=1.0 q=0,3 x F corr (x) q=0.7q=0.1 q=0.5 FIG. 3. (Color online) Position dependence of the nearest-neighbor correlations along the networkat different values of q . The normalized coordinate x = i/L is the same as in Fig. 2. we have α ∗ = 1 / (4 α ∗ ), so that α ∗ < / α ∗ > / Predictions of the domain wall theory. — An open chain with stationary current
J < / ρ − ( J ) and highdensity ρ + ( J ). According to the domain wall theory [15], on the coexistence line the twophases are separated by a completely delocalized domain wall. As a result, the averageddensity profile is linear, changing its value from ρ − ( J ) at the left end of the chain to ρ + ( J ) atits right end. This prediction is compared to numerical simulation data in Fig. 2 for externalrates α = β = 0 .
75, length of each segment L = 400, and q = 0 .
5. The data was averagedover 100 runs of length 2 attempted moves each. One sees a very good agreement betweenthe theoretical prediction ρ (2)1 = ρ − ( J (2) ) ≃ .
282 and the simulation result ρ (2)1 ≃ . ρ (2)400 = ρ + ( J (2) ) ≃ .
718 and the simulation result ρ (2)400 ≃ . F cor ( x ) = h τ (2) i τ (2) i +1 i − h τ (2) i ih τ (2) i +1 i as a function of the normal-ized distance x = i/L . The simulation results for all q show almost vanishing correlations inthe bulk of the first and third segments and a parabolic-like shape in the second segment,see Fig. 3. In the latter case the noticeable tilt of the ’parabolas’ to the right when q ≥ . G ( s,s +1) between the segments s and0 -0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1-0,010,000,010,020,030,040,050,060,070,080,09-0,1 0,0 0,1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 1,0 1,1-0,010,000,010,020,030,040,050,060,070,080,09 J sc , qmax F cor FIG. 4. (Color online) Comparison between the numerically estimated maximum value of thenearest-neighbor correlations in the second segment max x F cor ( x ), shown by red stars connectedwith a red line, and the current J sc through the shortcut, shown by blue circles connected with ablue line, at different values of the parameter q . s + 1, s = 1 ,
2. For example, we have numerically evaluated G (1 , ≃ − . q = 0 . q = 1 .
0, while G (2 , ≃ − . q = 0 . G (2 , ≃ . q = 1 .
0. Theoretically,the maximum value of F cor ( x ) is reached at the midpoint of the chain and equalsmax x F cor ( x ) = [ ρ + ( J (2) ) − ρ − ( J (2) )] / / − J (2) = J sc . The validity of this prediction of the domain wall theory is illustrated in Fig. 4.
Conclusions. — The model predicts that a shortcut in the bulk of a road carryingmaximum stationary current inevitably causes traffic jams characteristic of a shock phasewith completely delocalized domain wall. The main parameters of the average densityprofile and the nearest-neighbor correlations in the shunted segment are governed by a cubicequation with coefficients simply depending on the probability of choosing the shortcut.
Acknowledgements. — N.B. gratefully acknowledges support from a grant of the Repre-sentative Plenipotentiary of Bulgaria to the Joint Institute for Nuclear Research in Dubna.N.P. thanks Roumen Anguelov and Jean Lubuma for their hospitality at the University of1Pretoria, where a part of this work was also carried out. [1] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Biopolymers , 1 (1968).[2] K. Nagel, M. Schreckenberg, J. Physique I , 2221 (1992); D. Chowdhury, L. Santen, and A.Schadschneider, Phys. Rep. , 199 (2000); D. Helbing, Rev. Mod. Phys. , 1067 (2001).[3] A. Parmeggiani, T. Franosch, amd E. Frey, Phys. Rev. Lett. , 086601 (2003); A. Roux, G.Capello, J. Cartaud, J. Prost, B. Goud, and P. Bassereau, Proc. Natl. Acad. Sci. USA , 5394(2002); G. Koster, M. VanDuijn, B. Hofs, andM. Dogterom, Proc. Natl. Acad. Sci. USA ,15583 (2003); T. M. Nieuwenhuizen, S. Klumpp, and R. Lipowsky, Phys. Rev. E , 061911(2004); C. Leduc, O. Campas, K. B. Zeldovich, A. Roux, P. Jolimaitre, L. Bourel-Bonnet,B. Goud, J.-F. Joanny, P. Bassereau, and J. Prost, Proc. Natl. Acad. Sci. USA , 17096(2004); C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez, and J. Howard, PNAS , 6100 (2013); I. Neri, N. Kern, and A. Parmeggiani, Phys. Rev. Lett. , 098102 (2013).[4] J. Krug and H. Spohn, Phys. Rev. A , 4271 (1988); J. Krug, P. Meakin, and T. Halpin-Healy,Phys. Rev. A , 638 (1992); T. Sasamoto, J. Phys. A , L549 (2005).[5] T. Chou and D. Loshe, Phys. Rev. Lett. , 3552 (1999); Q.-H. Wei, C. Bechinger, and P.Leiderer, Science , 625 (2000); A. B. Kolomeisky, Phys. Rev. Lett. , 048105 (2007).[6] T. Huisinga, R. Barlovic, W Knopse, A. Schadschneider, M. Schreckenberg, Physica A ,249 (2001).[7] T. Reichenbach, E. Frey, and T. Franosch, New J. Phys. , 159 (2007).[8] T. Karzig and F. von Oppen, Phys. Rev. B , 045317 (2010).[9] J. Brankov, N. Pesheva, and N. Bunzarova, Phys. Rev. E , 066128 (2004).[10] E. Pronina, A. B. Kolomeisky, J. Stat. Mech.: Theory Exp., P07010 (2005).[11] P. Pierobon, M. Mobilia, R. Kouyos, and E. Frey, Phys. Rev. E , 031906 (2006); Z.-P.Cai, Y.-M. Yang, R. Jiang, M.-B. Hu, Q.-S. Wu, and Y.-H. Wu J. Stat. Mech.: Theor. Exp.,P07016 (2008); B. Embley, A. Parmeggiani, N. Kern, J. Phys.: Condens. Matter , 295213(2008); Z.-P. Cai, Y.-M. Yuan, R. Jiang, K. Nishinary and Q.-S. Wu - J. Stat. Mech.: Theor.Exp., P02050 (2009); X. Wang, R. Jiang, K. Nishinary, M.-B. Hu, and Q.-S. Wu - Int. J. Mod.Phys. C , 967 (2009); X. Wang, R. Jiang, M.-B. Hu, K. Nishinary, and Q.-S. Wu - Int. J.Mod. Phys. C , 1999 (2009); B. Embley, A. Parmeggiani, N. Kern, Phys. Rev. E , 041128 (2009); H.-F. Du, Y.-M. Yuan, M.-B. Hu, R. Wang, R. Jiang and Q.-S. Wu, J. Stat. Mech.:Theor. Exp., P03014 (2010); I. Neri, N. Kern, and A. Parmeggiani, Phys. Rev. Lett. ,068702 (2011); X. Song, L. Ming-Zhe, W. Jian-Jun, and W. Hua, Chin. Phys. B , 060509(2011); S. Xia, L. Tang, H. Wang, Cent. Eur. J. Phys., , 1077 (2011); C. Appert-Rolland, J.Cividini and H. J. Hilhorst - J. Stat. Mech., P10014 (2011); R. Chatterjee, A. K. Chandra,and A. Basu, Phys. Rev. E , 032157 (2013); A. Raguin,A. Parmeggiani,and N. Kern, Phys.Rev. E , 042104 (2013).[12] N. C. Pesheva and J. G. Brankov, Phys. Rev. E , 062116 (2013).[13] Y.-M. Yuan, R. Jiang, R. Wang, M.-B. Hu, and Q.-S. Wu, J. Phys. A , 12351 (2007).[14] B. Derrida, E. Domany, and D. Mukamel, J. Stat. Phys. , 667 (1992); B. Derrida, M. R.Evans, V. Hakim and V. Pasquier, J. Phys. A: Math. Gen. , 1493 (1993).[15] A. B. Kolomeisky, G. M. Sch¨uutz, E. B. Kolomeisky and J. P. Straley, J. Phys. A , 6911(1998); L. Santen and C. Appert, J. Stat. Phys.106