# Modeling Cytoskeletal Traffic: An Interplay between Passive Diffusion and Active Transport

Izaak Neri, Norbert Kern, Andrea Parmeggiani

► To cite this version:

Izaak Neri, Norbert Kern, Andrea Parmeggiani. Modeling Cytoskeletal Traffic: An Interplay between Passive Diffusion and Active Transport. Physical Review Letters, 2013, 110 (9), pp.098102. 10.1103/PhysRevLett.110.098102. [hal-00805162](https://hal.science/hal-00805162)

HAL Id: hal-00805162

https://hal.science/hal-00805162

Submitted on 4 Jun 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Modelling cytoskeletal traffic: an interplay between passive diffusion and active transport

Izaak Neri¹,², Norbert Kern¹,², and Andrea Parmeggiani¹,²,³,⁴

¹Université Montpellier 2, Laboratoire Charles Coulomb UMR 5221, F-34095, Montpellier, France

²CNRS, Laboratoire Charles Coulomb UMR 5221, F-34095, Montpellier, France

³Université Montpellier 2, Laboratoire DIMNP UMR 5235, F-34095, Montpellier, France

⁴CNRS, Laboratoire DIMNP UMR 5235, F-34095, Montpellier, France

(Dated: October 25, 2018)

We introduce the totally asymmetric exclusion process with Langmuir kinetics (TASEP-LK) on a network as a microscopic model for active motor protein transport on the cytoskeleton, immersed in the diffusive cytoplasm. We discuss how the interplay between active transport along a network and infinite diffusion in a bulk reservoir leads to a heterogeneous matter distribution on various scales. We find three regimes for steady state transport, corresponding to the scale of the network, of individual segments or local to sites. At low exchange rates strong density heterogeneities develop between different segments in the network. In this regime one has to consider the topological complexity of the whole network to describe transport. In contrast, at moderate exchange rates the transport through the network decouples, and the physics is determined by single segments and the local topology. At last, for very high exchange rates the homogeneous Langmuir process dominates the stationary state. We introduce effective rate diagrams for the network to identify these different regimes. Based on this method we develop an intuitive but generic picture of how the stationary state of excluded volume processes on complex networks can be understood in terms of the single-segment phase diagram.

PACS numbers: 87.16.A-, 87.16.Uv, 05.60.Cd, 89.75.Hc

## Introduction
Statistical physics has been very successful in deducing macroscopic properties of materials from the interactions between their microscopic components. Active matter systems on the other hand, such as active colloids, animal flocks or cytoskeletal assemblies, are prone to develop out-of-equilibrium patterns. These spatial heterogeneities are in fact essential to life processes [1].

<figure><img src="image_2.png" /><figcaption>FIG. 1. (color online) Statistical physics model of cytoskeletal transport, capturing the competition between active transport of particles through a network (cytoskeleton) with diffusion in a bulk reservoir (cytoplasm).</figcaption></figure>

Here we would like to initiate a microscopic statistical physics approach to describe collective organization of molecular motors in cells. Motor proteins navigate actively throughout the cell [2, 3] along the cytoskeleton, a network of filamentous assemblies spanning the cytoplasm. These proteins can exert forces, depolymerize filaments and transport biological cargos along the cytoskeleton [4], and thus play an important role in the assembly, self-organization and functioning of cells [5]. Single-molecule properties of the motors are well-studied and can now be measured accurately [6]. But understanding how motors collectively self-organize remains a very impor-

tant step in developing a microscopic vision of intracellular organization.

Macroscopic approaches to study intracellular motor protein transport have been developed [7], including efforts to introduce microscopic aspects [8, 9], but a generic microscopic picture of cytoskeletal transport has yet to emerge. Here we present a tentative approach using a minimalistic model for cytoskeletal active transport, which consists of directed motion of particles (motor proteins) along the network (cytoskeleton) and diffusion in the bulk (cytoplasm), see Fig. 1. We model the directed motion along the cytoskeleton using the totally asymmetric simple exclusion process (TASEP) [10] along a disordered directed network [11]. The binding and unbinding of particles between the network and the bulk is represented via Langmuir kinetics (LK) [12]. We consider the particles in the bulk to be infinitely diffusive which, as we will show, is a relevant limiting case.

TASEP is a fundamental model in non-equilibrium physics [13], but is also used in more applied topics, such as modelling macromolecules moving through capillary vessels [14], electrons hopping through a quantum-dot chain [15] and vehicular traffic [16]. LK on the other hand is a well known fundamental equilibrium process in chemical physics [12].

Our model constitutes a generalization of transport through closed networks [11] to open systems, as they are relevant for cytoskeletal transport. It may also be seen as generalization to a large scale network of the totally asymmetric exclusion process with Langmuir kinetics (TASEP-LK) on a single segment [17], which has been shown to quantitatively describe in-vitro experiments of motor proteins [18]. Our study differs from previous work [8] in that we consider exclusion interactions and disordered networks, both essential points to the physical picture we develop.

The fundamental question we address is how spatial hetero-
<figure><img src="image_3.png" /><figcaption>FIG. 2. (color online) Sketch of the method to study transport through complex networks (a), here for TASEP-LK. Each segment of the network (b) is considered to connect to two reservoirs with effective rates $\alpha_{\text{eff}}$ and $\beta_{\text{eff}}$ (c), which are set by the junction densities.</figcaption></figure>

geneities emerge in such open active systems, here due to the competition between diffusion in a reservoir (which spreads particles homogeneously) and active transport along a network (which generates heterogeneities [11]).

We will show how this competition is governed by the hopping rate for particles on the network, the exchange rates with the reservoir and the filament length. One major result of our work is that there exist three regimes of transport through complex networks which are linked to the scale at which spatial density heterogeneities arise: the network, the segments or the sites.

**Microscopic model for cytoskeletal transport** We represent the cytoskeleton as a network of directed segments of L sites each, connected by junction sites, see Fig. 2 (a). We use random networks to reflect the topological complexity of the cytoskeleton, a standard approach in modelling networks [19]. Specifically we use Erdös-Rényi graphs of mean connectivity c, as in [11]. Whereas the specific topology is not important for our qualitative results, the fact that networks are irregular, i.e. that the number of incoming $c_v^i$ and outgoing $c_v^o$ segments of the junction v differ, is relevant.

In each directed segment particles move according to TASEP-LK rules [17]: particles hop uni-directionally at rate p, subject to exclusion interactions. Furthermore, particles obey binding/unbinding kinetics with attachment rate $\omega_A$ and detachment rate $\omega_D$ at every site along the network. Particles in the reservoir are assumed to diffuse infinitely fast.

The phase diagram of TASEP-LK has been determined [17] for a single segment connecting a reservoir with entry rate $\alpha$ to a reservoir with exit rate $\beta$ (see supplemental materials [20]). TASEP-LK is best characterized in terms of the dimensionless parameters $\Omega_A = \omega_A L/p$ and $\Omega_D = \omega_D L/p$, where $\Omega_D$ relates the distance an isolated particle typically moves before detaching to the filament length [17]. It is convenient to consider the parameters $\Omega = (\Omega_A + \Omega_D)/2$, characterizing the total exchange between reservoir and segment, and the ratio $K = \Omega_A/\Omega_D$, which sets the equilibrium Langmuir density $\rho_L = K/(K+1)$.

The $(\alpha, \beta)$-phase diagram of TASEP-LK is fully characterized in terms of the parameters $\Omega$ and $K$ (see [17, 20] and Fig. 4). In the following we refer to low density (LD), high density (HD) and maximum current (MC, or M in [17])

phases. All of these reduce to TASEP phases with a flat density profile for $\Omega = 0$ [10]. Furthermore LD and HD zones can coexist on the same segment separated by a domain wall leading to a coexistence phase (LD-HD).

**Mean field method for TASEP-LK on networks** We analyze TASEP-LK by extending the mean field arguments for TASEP presented in [11]. In this approach every segment $(v, v')$ connecting two vertices v and v' is considered to be governed by effective entry and exit rates $\alpha_{\text{eff}}^{\text{eff}}$ and $\beta_{\text{eff}}^{\text{eff}}$, see Fig. 2 (c). These rates are in turn determined by the average densities $\rho_v$ and $\rho_{v'}$ at the junction sites [11, 24]: $\alpha_{v,v'}^{\text{eff}} = p \rho_v/c_v^o$ and $\beta_{v,v'}^{\text{eff}} = p (1 - \rho_{v'})$. Balancing the currents at the junctions leads to the following closed set of equations in $\rho_v$:

$$ \frac{\partial \rho_v}{\partial t} = \sum_{v' \to v} J^{-} \left[ \frac{\rho_{v'}}{c_{v'}^o}, 1 - \rho_v \right] - \sum_{v' \leftarrow v} J^{+} \left[ \frac{\rho_v}{c_v^o}, 1 - \rho_{v'} \right] \quad (1) $$

where the sums are over incoming (outgoing) segments, and $J^{\pm}$ [$\alpha/p, \beta/p$] are the currents entering (leaving) a segment with rates $\alpha$ ($\beta$). The expressions for $J^{\pm}$ are readily available [17, 20]. Due to the Langmuir process the current is not constant along the segment such that $J^{-} \neq J^{+}$. Solving Eqs. (1) yields the complete stationary state of all segments in the network.

The overall particle density on a network immersed in a reservoir is equal to the Langmuir density $\rho_L$, set by the ratio $K$ [23]. In the following we discuss the physical phenomena at fixed K as the exchange parameter $\Omega$ smoothly interpolates from purely active transport in standard TASEP ($\Omega = 0$) to the diffusion dominated limit ($\Omega \to \infty$). The data shown here has been obtained using $K = 1.5$, which is a reasonable value for motor proteins [22] and theoretically not specific in any way: our analysis is general, and only very high and very low values of K require an additional discussion [23].

**Decoupling due to particle exchange** In principle the continuity equations (1) couples the densities $\rho_v$ and $\rho_{v'}$ of those junctions which are linked by a segment $(v, v')$: it is this coupling which makes the transport problem global, requiring to analyze the whole network simultaneously.

Here we exploit one feature of TASEP-LK, which we state by saying that the binding-unbinding process can 'decouple' the currents at the segment boundaries. Indeed, for the composite LD-HD phase it is known [17] that the in-current $J^{-}$ depends on the in-rate $\alpha^{\text{eff}}$ only, whereas the out-current $J^{+}$ depends on the out-rate $\beta^{\text{eff}}$ only. Any LD-HD segment therefore lifts one coupling constraint in Eqs. (1), since the in/out currents $J^{\pm}$ are determined locally by the junction densities $\rho_v$ and the local connectivity.

When complete decoupling is achieved, as is expected at high values of $\Omega$, one can directly deduce the junction densities for an arbitrary network. As presented in the supplement [21], this leads to an exact solution of the mean field Eqs. (1).
<figure><img src="image_4.png" /><figcaption>FIG. 3. (color online) Average current as a function of the exchange parameter Ω, for three different (average) connectivities ($c = 2, c = 10$ and $c \to \infty$) at $K = 1.5$. Agreement between simulation (symbols, for $L = 400$) and mean-field results (solid lines, Eq. (1)) is excellent. The (red) dashed line is a simplified mean-field result (Eq. (2)) while the dotted line denotes the Langmuir current $J/p = \rho_{\ell}(1 - \rho_{\ell})$. Results are for a single graph instance of $O(10^2)$ number of junctions.</figcaption></figure>

For $K > 1$ we have

$$ \rho_v = \begin{cases} c_v^0 \left( \frac{c_v^i - 1}{c_v^i c_v^i - 1} \right) & \frac{\rho_v}{c_v^i} \le \frac{1}{K+1}; c_v^i \ne 1 \\ \frac{c_v^0}{2} \left( 1 - \sqrt{1 - \frac{c_v^i}{c_v^i}} \right) & \frac{\rho_v}{c_v^i} \le \frac{1}{K+1} \\ \frac{1}{2} + \sqrt{\frac{1}{4} - \frac{c_v^0}{c_v^i} \rho_{\ell}(1 - \rho_{\ell})} & \rho_v \ge \frac{c_v^0}{K+1}; \rho_v \ge \frac{1}{2} \end{cases} \quad (2) $$

In Fig. 3 we compare these analytical predictions based on complete decoupling to both a (numerical) full mean-field solution to Eqs. (1) and to simulations. For low values of Ω the fully coupled description Eqs. (1) is necessary (especially for low mean connectivity $c$). Surprisingly, the decoupled description Eq. (2) is excellent even down to relatively low values of Ω. This is an important result, as it shows that single segment TASEP-LK [17] suffices to describe transport through any complex network for a wide Ω range. For comparison we have also indicated the average mean-field current in the $c \to \infty$ limit, i.e. the TASEP-LK current which is maintained even when all junctions are blocked ($\alpha^{\text{eff}} = \beta^{\text{eff}} = 0$). This constitutes a lower bound to any current in any network.

**Effective rate plots** Here we introduce effective rate plots as a way to understand intuitively the physics of active transport through networks by allowing to visualize the whole stationary transport state of the network. In Figs. 4 we map the effective rates ($\alpha_s, \beta_s$) of each segment $s$, obtained by numerically solving of the full mean field Eqs. (1), onto the single segment phase diagram. Note that the scattering of effective rates is due to the irregularity of the networks considered.

For sufficiently high $c$ Figs. 4 reveal that the effective rates cluster close to the origin, in the LD-HD phase; this explains why the simplified Eqs. (2) work well for high $c$ in Fig. 3. When increasing Ω the TASEP-LK phase diagram changes. In particular the LD phase reduces in favour of the LD-HD

phase. Moreover, at high Ω one notices a specific alignment of the effective rates as given by the decoupled Eqs. (2).

In the following section we show how effective rate plots allow to rationalize the scale at which density heterogeneities appear in the network.

<figure><img src="image_5.png" /><figcaption>FIG. 4. (color online) Effective rate diagrams for irregular graphs of mean connectivity $c$ (at $K = 1.5$, for given values of Ω). Effective rates follow from solving Eqs. (1). Single graph instances consisted of $O(10^2)$ junctions.</figcaption></figure>

**Heterogeneities for TASEP-LK on irregular networks** The parameter Ω regulates the way particles distribute along the network, at overall density $\rho_{\ell}$, and thus determines how heterogeneities develop. We characterize the stationary state from the effective rate plots by determining the fraction of segments occupying the corresponding phases, see Figs. 5. A complementary point of view is given in Fig. 6 with the distribution $W(\rho_s)$ of the mean segment densities $\rho_s$ in the network. From these figures we conclude that heterogeneities develop throughout the network in three successive regimes (we discuss the case $K > 1$ and refer to [23] for $K \le 1$):

<figure><img src="image_6.png" /><figcaption>FIG. 5. (color online) Fraction of segments in LD, HD and LD-HD phases, in the effective rate diagram of mean connectivity (a) $c = 2$ and (b) $c = 10$, for $K = 1.5$. The transition $\Omega_c$ between the network and the segment regime (denoted by the vertical dotted line) is determined by the condition that $n_{\text{LD}}$ vanishes. The graph instances used are those of Fig. 4.</figcaption></figure>
(i) the network regime, for low exchange rates $\Omega$, is characterized by the presence of LD and HD segments. The distribution $W(\rho_s)$ is marked by the LD and HD peaks (whereas LD-HD coexistence segments are distributed evenly over the intermediate density range). This bimodality implies a strongly heterogeneous density at the network scale.

(ii) in the segment regime, for intermediate exchange rates all LD segments have disappeared in favor of LD-HD segments. The distribution $W(\rho_s)$ is dominated by the LD-HD peak. Although all segments have similar average densities, the presence of domain walls implies strong inhomogeneities on the segment scale.

(iii) in the Langmuir regime, for large exchange rates $\Omega$, the Langmuir phase dominates. All segments have homogeneous densities except for small regions near the boundaries, and no heterogeneities arise beyond the scale of a few sites.

<figure><img src="image_7.png" /><figcaption>FIG. 6. (color online). The MF distribution $W(\rho_s)$ of segment densities for an irregular graph instance with mean connectivity (a) $c = 2$ and (b) $c = 10$, of $O(10^4)$ junctions, for $K = 1.5$. At low values the bimodal distribution of [11] is identified. When increasing $\Omega$ the center peak gradually grows while the edge peaks shrink. Network heterogeneities eventually disappear at intermediate $\Omega = \Omega_c$, leading to a unimodal density distribution.</figcaption></figure>

The transition between the network and segment regimes is sharp (identified as the point $\Omega_c$ where all LD segments disappear, $n_{LD} = 0$). Moreover this transition has an upper bound $\Omega^* = 1/2 + f(K) \ln[f(K)/(1/2 + f(K))]$ (with $f(K) = |K - 1|/(2(K + 1))$), for which the LD phase disappears from the one-dimensional phase diagram [23]. In contrast, the crossover from the segment regime to the Langmuir regime is progressive.

**Conclusions** We have analyzed active transport on a disordered network which is immersed in an infinitely diffusive bulk reservoir, as may be considered a simple model for cytoskeletal transport by molecular motors. The dimensionless parameter $\Omega$ characterizes the competition between active and passive transport, and how the associated collective effects lead to strong spatial inhomogeneities.

Three regimes arise, according to the scale at which heterogeneities appear in the network: a *network*, a *segment* and a *Langmuir* (site-dominated) regime, see Fig. 7. Interestingly, these scales also set the complexity that characterizes the theoretical analysis. In the network regime TASEP-LK transport is coupled throughout the network, whereas in the opposite

<figure><img src="image_8.png" /><figcaption>FIG. 7. Three regimes with heterogeneities on a network, segment and site scale, respectively, according to the exchange parameter $\Omega$ (see main text for $\Omega_c$ and $\Omega^*$). Particle densities are coded in grey scale. Pictorials on top indicate the scale at which continuity equations are solved.</figcaption></figure>

Langmuir regime the physics is essentially determined by the attachment/detachment process. In the intermediate segment regime decoupling implies that the transport characteristics of each segment follow readily from those of a single segment.

Effective rate plots allow to intuitively understand transport processes through networks from the single-segment transport characteristics; from the scattering of rates over both the LD and HD zones we can directly deduce the role of strong heterogeneities, see Fig. 4. This approach yields valuable a priori insight into yet more complex excluded volume transport such as TASEP with extended particles [25], TASEP with multiple species [26] and, as we show in [23], bidirectional motion [27]: as in TASEP-LK, the single-segment phase diagram serves as a basis for deducing the behaviour on the network. From the effective rate diagram approach it becomes also clear that our results extend to types of disorder other than topological which are relevant to biological systems, e.g. disorder in the actions of particles at the junctions [24, 28]. An interesting open question is how finite diffusion [29] could be handled within our approach.

Several conclusions may be relevant for modelling cytoskeletal transport. First, we have shown that strong inhomogeneities in the spatial distribution of motor proteins form for a wide range of parameters. Even in the case of infinitely fast diffusion considered here they resist the equalizing effect of bulk diffusion. Inhomogeneities would therefore be even more relevant for finite diffusion. Second, the presence of some exchange of motors between the cytoskeleton and the cytoplasm may in fact simplify a theoretical description, since the approximation of decoupling Eqs. (2) yields excellent results for large enough $\Omega > \Omega_c$ (see Fig. 3 and 7). Third, our analysis hints at a way to regulate the spatial distribution of motors, and therefore their cargos in the cell, by way of modifying the exchange parameter $\Omega$. It can be controled in independent ways, via $\omega_D$ or $\omega_A$ (through the bare biochemical rate or through the motor concentration), or via the length dependence in $\Omega$: regulating the cytoskeleton mesh size, for example by crosslinker proteins, would allow to control the
length scale of heterogeneities. Values reported in the litera-
ture [18, 22] show that the values used here ($\Omega, K \sim O(1)$)
is a reasonable order of magnitude, and it is therefore tempt-
ing to speculate that a moderate regulation of $\Omega$ might indeed
allow to provoke a crossover between the various regimes in
living cells.

We acknowledge support from ANR-09-BLAN-0395-02
and from the Scientific Council of the University of Mont-
pellier 2. We thank C. Leduc for discussing and useful refer-
ences.

[1] S. Ramaswamy, Annu. Rev. Condens. Matt. Phys. **1**, 323-345 (2010); M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, R. A. Simha, http://arxiv.org/pdf/1207.2929v1.pdf.; T. Vicsek, A. Czirok, E. Ben-Jacob, I. Cohen, O. Shochet, Phys. Rev. Lett. **75**, 1226-1229 (1995)

[2] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, P. Walter, *Molecular biology of the cell*, edition 5, Garland science, Taylor & Francis Group, New York USA (2008).

[3] J. L. Ross, M. Y. Ali, and D. M. Warshaw, Curr. Opin. Cell Biol., **20**, 4147 (2008).

[4] J. Howard, *Mechanics of Motor Proteins and the Cytoskeleton*, Sinauer Associates, Inc. (2001).

[5] D. A. Fletcher, R. D. Mullins, Nature **463**, 485-492 (2010).

[6] S. Manley, J. M. Gillette, G. H. Patterson, H. Shroff, H. F. Hess, E. Betzig, J. Lippincott-Schwartz, Nat. Methods **5**, 155-157 (2008); J. Yoo, T. Kambara, K. Gonda, H. Higuchi, Experimental Cell research **314**, 3563 (2008); P. Pierobon, S. Achouri, S. Courty, A. R. Dunn, J. A. Spudich, M. Dahan, G. Capello, Bio-physical Journal, **96**, 4268 (2009)

[7] F. Nédélec, T. Surrey, A. C. Maggs, Phys. Rev. Lett. **86**, 3192-3195 (2001); P. K. Trong, J. Guck, R. E. Goldstein, Phys. Rev. Lett. **109**, 028104 (2012)

[8] S. Klumpp, R. Lipowsky, Phys. Rev. Lett. **95**, 268102 (2005); S. Klumpp, T. M. Nieuwenhuizen, R. Lipowsky, Biophys. J. **88**, 3118-3132 (2005)

[9] P. Greulich, L. Santen, The European Physical Journal E: Soft Matter and Biological Physics **32**, 2, 191-208 (2010).

[10] J. T. MacDonald, J. H. Gibbs, A. C. Pipkin, Biopolymers **6**, 1 (1968); B. Derrida, E. Domany, D. Mukamel, J. Stat. Phys. **69**, 667 (1992); T. Chou, K. Mallick, R. K. P. Zia, Reports on Progress in Physics **74**, 116601 (2011).

[11] I. Neri, N. Kern, A. Parmeggiani, Phys. Rev. Lett. **107**, 068702 (2011).

[12] R. Fowler, "Statistical Mechanics", Cambridge University Press (1936).

[13] K. Mallick, J. Stat. Mech. P01024 (2011).

[14] T. Chou, D. Lohse, Phys. Rev. Lett. **82**, 3552 (1999).

[15] T. Karzig and F. von Oppen, Phys. Rev. B **81**, 045317 (2010).

[16] A. Schadschneider, Physica A **285**, 101 (2001); D. Helbing, Rev. Mod. Phys. **73**, 1067 (2001).

[17] A. Parmeggiani, T. Franosch, E. Frey, Phys. Rev. Lett. **90**, 086601 (2003); A. Parmeggiani, T. Franosch, E. Frey, Phys. Rev. E **70**, 046101 (2004).

[18] C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez, J. Howard, Proc Natl Acad Sci U S A. **109**, 6100 (2012).

[19] A. Barrat, M. Barthélemy, A. Vespignani, *Dynamical Processes on Complex Networks*, Cambridge University Press (2008).

[20] Supplementary Material 1: The full mean field equations

[21] Supplementary Material 2: Decoupling at high exchange rates: the simplified mean field equations

[22] V.C. Varga, C. Leduc, V. Bormuth, S. Diez, J. Howard, Cell **138**, 1174-1183 (2009); L. Reese, A. Melbinger, E. Frey, Biophys. J. **101**, 2190-2200 (2011)

[23] I. Neri, N. Kern, A. Parmeggiani, *to be submitted* (2012).

[24] B. Embley, A. Parmeggiani, N. Kern, Phys. Rev. E **80**, 041128 (2009).

[25] L. B. Shaw, R. K. P. Zia, K. H. Lee, Phys. Rev. E **68**, 021910 (2003).

[26] S. Muhuri, I. Pagonabarraga, EPL **84**, 58009 (2008)

[27] S. Sandow, Phys. Rev. E **50**, 2660-2667 (1994).

[28] A. Raguin, A. Parmeggiani, N. Kern to be submitted

[29] M. Ebbinghaus, L. Santen, J. Stat. Mech., P03030 (2009); H. Hinsch, E. Frey, Phys. Rev. Lett. **97**, 095701 (2006); K. Tsekouras, A B Kolomeisky, J. Phys. A: Math. Theor. **41**, 465001 (2008).