Oct 7, 2006 - Theoretical modeling of these systems has been achieved by the use ..... the lattice is homogeneous and given as Langmuir isotherm ÌºL =...

0 downloads 10 Views 276KB Size

No documents

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models Hauke Hinsch1 , Roger Kouyos2, and Erwin Frey1 1

2

Arnold Sommerfeld Center and CeNS, Department for Physics, Ludwig-Maximilians-Universit¨ at M¨ unchen, Theresienstrasse 37, D-80333 M¨ unchen, Germany Theoretical Biology, Eidgen¨ ossische Technische Hochschule Z¨ urich, Universit¨ atsstrasse 16, CH-8092 Z¨ urich, Switzerland

Abstract. Motor proteins are key players in intracellular transport processes and biological motion. Theoretical modeling of these systems has been achieved by the use of step processes on one-dimensional lattices. After a comprehensive introduction to the total asymmetric exclusion process and some analytical tools, we will give a review on different lines of research attracted to the aspects of this systems. We will focus on the generic properties of a coupling between the exclusion process and Langmuir bulk kinetics that induce topological changes in the phase diagram and multi-phase coexistence.

1

Introduction

The identification of motion as a manifestation of biological life dates back to the earliest records of science itself. The Greek physician Erasistratos of Ceos studied biological motion on the length scale of muscles already in the 3rd century BC. He imagined muscles to function in the way of a piston contracting and relaxing from pneumatic origin. It was not until the invention of the microscope in the 17th century by van Leeuwenhoek that this theory could be devalidated with Swammerdams observation that muscles contract at constant volume. Concerning biological motion on a microscopic scale, scientists favored concepts of “living forces” for many centuries until this was finally ruled out by the observations of the Scottish botanist Robert Brown in 1828 who found all kind of matter to undergo erratic motion in suspensions. A satisfactory explanation was provided by Einstein in 1905 by the interaction with thermally fluctuating molecules in the surroundings. However, the molecular details remained unknown in the fog of low microscope resolution. Modern experimental techniques [1] have lately revealed the causes of sub-cellular motion and transport. Today we know that every use of our muscles is the collective effort of a class of proteins called myosin that “walk” on actin filaments. Generally spoken, we refer to all proteins that convert the chemical energy of ATP (adenosinetriphospate) in a hydrolysis reaction into mechanical work as molecular motors. These motors are highly specialized in their tasks and occur in a large variety: ribosomes move along mRNA strands while translating the codons into proteins, dynein is responsible for cilia motion and axonal transport, and kinesins play a

2

Hinsch et al.

key role in cytoskeletal traffic and spindle formation (for an overview see [2] and references therein). While the exact details of the molecular structure and function of motor proteins [3] remain a topic of ongoing research, on a different level attention was drawn to phenomena that arise out of the collective interaction of many motors. Early research along this line was motivated by mRNA translation that is managed by ribosomes. Ribosomes are bound to the mRNA strand with one subunit and step forward codon by codon. The codon information is translated into corresponding amino acids that are taken up from the cytoplasm and assembled into proteins. To increase the protein synthesis many ribosomes can be bound to the same mRNA strand simultaneously. This fact might induce collective properties as was first realized by MacDonald [4] who set up a theoretical model for the translation of highly expressed mRNA. The importance of effects caused by the concerted action of many motors can be deduced from a very simple example that has yet drastic consequences: the slow down of ribosomes due to steric hindrance caused by another ribosome in front – comparable to an intracellular traffic jam that might significantly slow down protein synthesis. A theoretical approach to collective phenomena in intracellular traffic will try to simplify the processes of molecular motion down to a single step rate rather than focus on the chemical or mechanical details on the molecular level of motor steps. Then it becomes possible to model and analyze the behavior of several motors with the tools of many-body and statistical physics. We will start this review with a short introduction on this single step model in Sec. 2 before we introduce the total asymmetric exclusion process (TASEP) as a theoretical model for intracellular transport. Sec. 3 describes the stationary states and density distributions and their phase diagram as a function of boundary conditions. After a review on several recent extensions in Sec. 4, we will focus on the competition of TASEP and bulk dynamics in Sec. 5. Before concluding, Sec. 6 contains further recent developments.

2

Model and Methods

In the quest for a theoretical model for the motion of molecular motors the first and simplest choice may be the use of a Poisson process. The “Poisson stepper” is assumed to be an extensionless object advancing stochastically in discrete steps along a one-dimensional periodic lattice. The process is uni-directional as the position of the stepper can be described by x(t) = a n(t) with the discrete step size a and the random variable n(t) being a sequence of growing integers. Step events occur stochastically with a rate r constant in both space and time. Consequently, the average time between two steps is then given by the dwell time τ = 1/r and the probability to find the “Poisson stepper” at a position n after time t by the Poisson distribution [5]. After we have defined a model for the translocation of a single motor, we proceed with our original task which aims at the understanding of collective properties of many motors. Of course, more elaborate models have been estab-

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

3

lished [6,7] that account for several rate limiting steps – examples are the ATP supply or the availability of amino acids for ribosomal mRNA translation. However, the very basic “Poisson stepper” is chosen for reasons of simplicity and in order to prevent unnecessary molecular details from masking collective effects. Still, the validity and limitations of this simplification have to be kept in mind. Being supplied with the dynamics of a single motor, a stage for the concerted action of many can now be set up. This was first done in a pioneering work by MacDonald [4] and is now widely know as the total asymmetric exclusion process (TASEP). It consists of an one-dimensional lattice (Fig. 1) with N sites labeled by i = 1, · · · , N and with a spacing of a = L/N , where L is the total length of the lattice. For convenience, L is often set to 1 and the lattice spacing then referred to as ε = 1/N .

Fig. 1. Schematic model of TASEP (particles are injected with rate α, move exclusively to the right being subject to hard-core exclusion, and are removed with rate β)

Particles have an extension of the size of the lattice spacing and are subjected to hard core exclusion due to steric hindrance. Therefore the occupation number ni of site i can only take the values 0 or 1. Particles on the lattice attempt jumps to their right neighboring site with a rate r, which will be set to unity in the following. Hereby, a reference time scale is set. The effective frequency of jumps can be much smaller than r when attempted jumps are rejected due to an already occupied target site. The attempted jump rate to the left is zero, since we deal with a total asymmetric exclusion process, in contrast to the asymmetric exclusion process or the symmetric exclusion process, where the jump rate to the left is non-zero or even equal to the jump rate to the right. Unless one uses periodic boundary conditions, specific dynamic rules have to be defined at the boundaries, which play a crucial role in the solution of the process. Among different other conditions (reflective, open with a blockage) the most common type are open boundaries, which we will use as well: at the left boundary (i = 1) particles attempt to attach with a rate α, while they detach at the right boundary (i = N ) with rate β. This is equivalent to two additional sites i = 0 and i = N + 1 at the boundaries, which are connected to the system by the bulk dynamics described above, and are constantly set to the density α and 1 − β respectively.

4

Hinsch et al.

In spite of its simplicity, TASEP shows a wide range of interesting properties. Since it was propelled into the scope of statistical physicists, it has become a paradigm for non-equilibrium physics. In contrast to equilibrium systems it lacks detailed balance but evolves into a steady state where a non-vanishing current is maintained between boundaries. Upon varying these boundary conditions, TASEP was found to exhibit phase transitions which – following general theorems [8] – are not even allowed for one-dimensional equilibrium systems in the absence of long-range interactions. However, the analysis of non-equilibrium systems is considerably complicated by the lack of universal concepts like the Boltzmann-Gibbs ensemble theory. Feasible methods exist nevertheless and will be explained in the next section.

3

Density and Current in Stationary States

In analyzing exclusion processes research can focus on a multitude of different properties. The probably most obvious to address is the density and current distribution in the stationary state. This is motivated by two reasons. On the one hand, one intuitively attributes a strong importance to density information with respect to the biological background as e.g. the ribosome density is connected to the rate of protein synthesis. On the other hand, promising experimental techniques can measure motor densities and may allow for validation of theoretical models. Of course, quite extensive research has also been attracted to a multitude of different properties like correlation functions [9,10], relaxation properties [11] or super-diffusive spreading of fluctuations [12] which will not be the topic of this review. We will focus on analytical methods (supported by numerical simulations) that are designed to investigate spatial density distributions in the stationary state of the system. To this end we will introduce some basic tools that have proven useful in the exploration of TASEP properties. These are based on mean-field approximations and reproduce many results that can also be derived exactly. We are well aware that this approach neglects correlations as included in the exact solutions that have been achieved for the TASEP density profile by applying either recursion relations [13] or a quantum Hamilton formalism with Bethe ansatz [14]. 3.1

Quantum Mechanics and Statistic Properties

As an introduction we will outline some general statistical properties. At any given moment, the system can be found in a certain configuration µ made up of the occupation numbers at each lattice site. The next occurring stochastic event (i.e. the jump of one particle to a neighboring site) will therefore change the system to another configuration µ′ . The transition probability pµ→µ′ is independent of the way the system had reached the initial configuration. Since there is no memory of the system’s history, but any transition probability solely depends on the preceding state, TASEP is a Markov process. In order to describe the

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

5

system’s evolution, we can then use a master equation for the probability to find the system in a certain state. X dP (µ) = [ωµ′ →µ Pµ′ (t) − ωµ→µ′ Pµ (t)] , dt ′

(1)

µ 6=µ

where the ωµ→µ′ are the transition rates from one configuration µ to another µ′ . How can we now translate this general property into a description of TASEP? To this end we will use a convenient notation, which applies methods from the quantum mechanics toolbox in order to formulate the master equation in terms of operators. It was introduced as “quantum Hamiltonian formalism” and allows for exact solutions [14]. We introduce operators n ˆ i (t), which act as occupation number operators, measuring the presence (ni = 1) or absence (ni = 0) of a particle at site i. This results in the Heisenberg equation (for an introduction see e.g. [14]) d n ˆ i (t) = n ˆ i−1 (t)[1 − n ˆ i (t)] − n ˆ i (t)[1 − n ˆ i+1 (t)] , dt

(2)

where the first term on the right hand side constitutes the jump of a particle from the left neighboring site to site i (and thus a particle gain) and the second term a jump from site i to the adjacent lattice site on the right (a particle loss). Note the intrinsic exclusion constraint in both terms that prevents jump events if the destination site is occupied, i.e. the expression in brackets equals zero. If one expresses these gains and losses in current, it becomes convenient to use the current operator ˆji (t) = n ˆ i (t)[1 − n ˆ i+1 (t)] . (3) This allows to rewrite (2) as a discretized form of a continuity equation with the discrete divergence ∇ˆji (t) = ˆji (t) − ˆji−1 (t): d n ˆ i + ∇ˆji (t) = 0 . dt

(4)

Similar equations for the boundaries are readily derived in the same way. Since we are interested in the average density on a certain lattice site, we need to compute the time (or ensemble) average of the operators. Equation (2) gives an equation of motion for the operator and allows to solve for the time evolution of hˆ ni (t)i. In executing the ensemble average of (2) two-point correlation functions like hˆ ni−1 (t)(1−ˆ ni (t))i appear on the right hand side. These correlation functions again are connected to higher order correlations via their equations of motion. The resulting infinite series of correlation functions can be solved exactly for special cases only. Generally, one is required to use mean-field approaches. 3.2

Mean-Field Solution

The rather blurry term “mean-field theory” is based on the concept of using time or space averages (e.g. by neglecting temporal or spatial fluctuations) and has

6

Hinsch et al.

found a wide range of applications in statistical physics (see [15]). In this chapter, we will explain its implementation for the TASEP and show a possible solution. To point out the use of mean-field theory in TASEP, we look again at the average of the operator n ˆ i (t). We are interested in the stationary state and therefore averaging signifies either a time or an ensemble average, since the system is ergodic. Then the average returns the density at the considered site i as ̺i = hˆ ni i. Performing the average over (2), leaves us with the difficulty of the infinite series of correlation functions mentioned earlier. The mean field approximation consists now in neglecting any correlations by setting e.g. hˆ ni n ˆ j i = hˆ ni ihˆ nj i (see [16]). In our case, we obtain for the current hˆ ni (t)(1 − n ˆ i+1 (t))i = hˆ ni (t)i(1 − hˆ ni+1 (t)i) .

(5)

This allows then, to rewrite (2) in the stationary state (d̺i (t)/dt = 0) as 0 = ̺i−1 (1 − ̺i ) − ̺i (1 − ̺i+1 ) .

(6)

Obviously, (6) could easily be solved numerically, since it forms a system of N difference equations. However, it is possible to reduce the set of equations to arrive at an explicit solution by making two assumptions. First, note that the stationary state condition d̺i (t)/dt = 0 implies a conservation of current throughout the bulk, as can be seen from (6). Ergo, we just have to solve one equation out of the set, to determine the stationary, homogeneous current (and density). For that purpose, we use the continuum approximation, which turns the spatial lattice variable quasi-continuous. This is achieved for a large number N of lattice sites on a lattice of normalized length L = 1. The distribution of sites then approaches a continuum as ε = L/N ≪ 1 and x = i/N is rescaled to the interval 0 ≤ x ≤ 1. Thereby an expansion in powers of ε is allowed. 1 ̺(x ± ε) = ̺(x) ± ε∂x ̺(x) + ε2 ∂x2 ̺(x) + O(ε3 ) 2

(7)

Using this expansion in (6) and the corresponding equations for the boundaries, results in the following first-order differential equation if we neglect all terms with higher orders in ε: (2̺ − 1)∂x ̺ = 0 . (8)

The corresponding boundary conditions are ̺(0) = α and ̺(1) = 1 − β. Because we have a first order differential equation that needs to satisfy two boundary conditions, we are evidently concerned with an over-determined boundary value problem. There are three solutions to (8): ̺bulk (x) = 1/2 does not satisfy either boundary condition (except for the special case α = β = 1/2), while ̺(x) = C can satisfy either the left or the right boundary condition, resulting in ̺α (x) = α and ̺β (x) = 1 − β, respectively. 3.3

Phase Diagram and Domain Wall Theory

To obtain a general solution satisfying the boundary conditions, both solutions need to be matched. Consider the case α, β < 1/2. To meet the boundary condition at both sides, the global density function ̺(x) has to be ̺α (̺β ) in an

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

7

environment close to the left (right) boundary. Since both ̺α and ̺β are uniform, the two solutions do not intersect. At this point, we have to go beyond mean-field theory and assume that at any given time both solutions are valid in non-overlapping areas of the system. Where those areas border, they are connected by a sharp domain wall (DW) at position xw (see Fig. 2) ̺α for 0 ≤ x ≤ xw , ̺(x) = (9) ̺β for xw ≤ x ≤ 1 . From the dynamics of this domain wall [17] we can deduce important informa-

ρ

1

1 ρβ=1−0.2

β MC

jβ

jα 0.5

1/2

ρα=0.3 0 0

LD

HD xw

x

1

0

1/2

α

1

Fig. 2. (left) Schematic density distribution for α = .3, β = .2: in this situation the particle current jα exceeds the current jβ , thus carrying more particles to the domain wall which is then shifted to the left. (right) Phase diagram of TASEP in α, β phase space shows a low density (LD), high density (HD) and maximal current (MC) phase

tion about the system. The key point of domain wall theory is the identification of particle currents as the cause of domain wall motion. If for example the current jα = ̺α (1 − ̺α ) of the left density solution ̺α is higher than the current in the right part of the system (corresponding to α > β), particles are transported faster to the DW from the left end then they can head on to the right. Thus, the domain wall is shifted to the left. The system is filled up, until finally the whole bulk density has taken the value of ̺β except for a small boundary layer 1 at the system’s left boundary. The opposite happens, if jβ > jα , which is the case for β < α. Hence we can empirically state that the boundary with the smaller rate acts as a bottleneck and imposes its density distribution on the system. To quantify this behavior, one can use a traveling wave solution [17] of the form ̺(x − V t) to obtain the domain wall velocity as V =

jβ − jα =β−α. ̺β − ̺α

(10)

To arrive at a phase diagram in α, β-space we will analyze the bulk density. A positive DW velocity is obtained for α < β and pushes the DW to the right 1

The boundary layer is necessary in order to fulfill the both boundary conditions. Its extend is finite for small systems, but will vanish in the thermodynamic limit N → ∞.

8

Hinsch et al.

side of the system resulting in a bulk density smaller 1/2 called low-density phase (LD). On the contrary, β < α leads to a high-density phase (HD), as illustrated in the phase diagram (Fig. 2). LD and HD phase are connected by a first-order transition. Note, that this phase diagram is equally obtained by the above mentioned exact methods. The transition line α = β < 1/2 in the phase diagram requires special treatment. The velocity of the domain wall yields zero here, but Monte Carlo simulations show that the DW actually will make random steps to either side. This behavior is caused by the stochasticity of the input and output events that are described as rate processes. Hence, the DW makes random steps to either side with equal probability, being nothing else than the famous random walk in a domain with reflecting boundaries. An average over a sufficiently long time will therefore result in a homogeneous probability density over space. The stationary density profile will just be a linear slope connecting the two boundary conditions: ̺(x) = α + (1 − β − α)x .

(11)

Note that mean-field theory neglects fluctuations and hence fails to predict this density distribution, but gives the result of a stationary domain wall. However, coming back to the picture of the domain wall as a random walker, the linear density distribution can readily be made plausible. Interpreting the random walk as free diffusion of the domain wall and keeping in mind that the density constraints ̺(0) = α and ̺(1) = 1 − β hold at all times, the linear density profile is easily derived as solution of the one dimensional diffusion equation with boundary conditions. Finally, consider an increase of e.g. α to values that exceed 1/2, while β > 1/2 (crossing the boundary from the upper left to the upper right quadrant in the phase diagram). In this situation particles are removed sufficiently fast from the system and the supply of particles on the left side acts as a bottleneck limiting density and current. For α < 1/2 any increase in α results in an increased current (think of a highway, where an additional car will result in a higher overall traffic). But above a certain value αC = 1/2 any increase in particle input will not increase the current further 2 but will cause the current to diminish again (as an additional car will further slow down traffic during rush hour). As a result the bulk will keep its current maximum at ̺ = 1/2 and a boundary layer will form at the left side of the system to match the boundary condition. This is of course nothing else than the bulk solution ̺bulk with two boundary layers and was baptized maximal current phase (MC). It is reached via second-order transitions for values α > 1/2 and β > 1/2. A more rigorous treatment of this behavior can be gained by computing the collective velocity of the particles [17]. 2

As j = ̺(1 − ̺) = ̺ − ̺2 has a maximum at ̺ = 1/2. This density-current relation is a convenient tool to characterize traffic processes and its plot is often referred to as fundamental diagram.

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

4

9

Biologically Motivated Generalizations of TASEP

The adoption of lattice gas models for biological systems was followed by a variety of efforts to fit TASEP to different realistic environments. For that purpose some of the simplifying assumptions TASEP is based on had to be questioned. For example, experimental observations have shown that ribosomes typically cover an area on the mRNA that exceeds the lattice spacing by a multiple. To account for this situation the particles in TASEP have to extend over several lattice sites. However it was found that this does not change the phase diagram quantitatively [18]. Another direction of research went back to the original field of MacDonald to elucidate the importance of initiation and prolongation of ribosomes [19] or formation of mRNA loops to facilitate the back transport of ribosomes from the termination site [20]. While particle interactions in TASEP is limited to hard-core potentials, even a small increase in the interaction radius – as it could be caused by charged molecules – leads to qualitative changes in the phase diagram [21]. It was shown that lattice gases with short-range repulsive interaction exhibit a density-current relation with two local maxima in contrast to simple TASEP that leads to one maximum. This behavior results in a qualitative change of the phase diagram that is enriched by four more regions one being a minimal-current phase. Further work was dedicated to the scenario of the interaction of different species of molecular motors that move in opposite directions. This can either happen on the same filament when two different particles are able to surpass each other with a jump rate that differs from the jump rate of either particle to a free site [22] or on two adjacent one-dimensional filaments [23]. In both cases, spontaneous symmetry breaking was observed. Intracellular transport along cytoskeletal filaments has also served as a source of inspiration for driven lattice gas models. While in the TASEP model motors can only bind and unbind on the left and the right boundary respectively, cytoskeletal motors are known to detach from the track to the cytoplasm [2] where they perform Brownian motion and subsequently reattach to the track. The interplay between diffusion in the cytoplasm and directed motion along the filament was studied [24] both in open and closed compartments, focussing on anomalous drift and diffusion behavior, and on maximal current and traffic jams as a function of the motor density. In [25] it has been realized that the on-off kinetics may not only give rise to quantitative changes in the transport efficiency but also to a novel class of driven lattice gas models. It was shown that the interplay between bulk on-off kinetics and driven transport results in a stationary phase exhibiting phase separation. This was achieved by an appropriate scaling of the on-off rates that ensures that particles travel a finite fraction on the lattice even in the limit of large systems. Then, particles spend enough time to “feel” the their mutual interaction and, eventually, produce collective effects. In the following section, we will review the results of these studies [25,26].

10

5

Hinsch et al.

Phase Coexistence

The essential features of cytoskeletal transport are the possibility of bulk attachment and detachment and a finite residence time on the lattice. The latter can be understood as a effect of thermal fluctuations that may overcome the binding energy of the motors that is only of the order of several kB T . Hence, attachment and detachment is a stochastic process whose dynamic rules have to be defined. Parmeggiani et al. [25] chose to use Langmuir kinetics (LK) known as adsorption-desorption kinetics of particles on a one- or two-dimensional lattice coupled to a bulk reservoir [27]. Particles can adsorb at empty sites and desorb from occupied sites and microscopic reversibility demands that the kinetic rates obey detailed balance leading to an evolution towards an equilibrium steady state describable by standard concepts of equilibrium statistical mechanics. In this sense the choice of LK is especially tempting as we are now faced with the competition of two representatives of both equilibrium and non-equilibrium systems. The system – in the following referred to as TASEP/LK – is defined as follows: the well-known TASEP is extended with the possibility of particles to attach to the filament with rate ωA and to detach from an occupied lattice site to the reservoir with rate ωD . According to the type of ensemble (canonical, grand canonical) the reservoir is either finite or infinite. Here, the reservoir is assumed to be infinite and homogeneous throughout space and time. The density on a lattice reached in the equilibrium state of LK is only dependent on the ratio K = ωD /ωA and is completely uncorrelated in both space and time for neglection of any particle interaction except hard-core. This is justified by the assumption that the diffusion in the cytoplasm is fast enough to flatten any deviations of the homogeneous reservoir density. The resulting density profile on the lattice is homogeneous and given as Langmuir isotherm ̺L = K/(K + 1).

Fig. 3. Schematic model of TASEP/LK: the TASEP is extended by possible particle attachment and detachment in the bulk with rate ωA , ωD

If we now consider the combination of TASEP and LK into the model displayed in Fig. 3, attention has to be paid to the different statistical nature of both processes. TASEP evolves into a non-equilibrium state carrying a finite current. Since particles are conserved in the bulk, the system is very sensitive

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

11

to the boundary conditions, whereas LK as a equilibrium process is expected to be robust to any boundary effects especially for large systems. Combining both processes would thus lead to a trivial domination of LK as the bulk rates ωA and ωD that apply to a large number of bulk sites become predominant over the rates α and β that only act on the two boundary sites. To observe any interesting behavior (i.e. real interplay) between the two dynamics, one needs competition. A prerequisite for the two processes to compete are comparable jump rates. To ensure that the rates are of the same order independently of the system size, an appropriate scaling is needed. To this end a N -independent global detachment rate ΩD is introduced, while the local rate per site scales as ωD =

ΩD . N

(12)

For the attachment one proceeds similarly. What does this scaling signify physically? For an explanation it is instructive, to have a look at the time scales involved: a particle on the lattice will perform a certain move on an average time scale, which is the inverse of that moves rate. Therefore, a particle spends an average time τ ≈ 1/ωD on the lattice before it detaches. Bearing in mind that the TASEP jump rate is set to unity, a particles will jump to its adjacent site after a typical time of one unit time step. Therefore the particle will travel a number NT = 1/ωD of sites before leaving the lattice. Compared to the lattice length this corresponds to a fraction of nT = NT /N = 1/(N ωD ). In order to keep this fraction finite in the thermodynamic limit, ωD needs to scale as defined in (12). Only if the fraction nT is finite, a given particle can experience interaction with other particles and give rise to collective phenomena [26]. 5.1

Mean Field Solution of TASEP/LK

To obtain density and current distributions of the TASEP/LK, we use again a mean-field approach and proceed along the lines of Sec. 3.2. First of all, we need to account for the Langmuir kinetics. This is done by adding the following terms to the Heisenberg equation (2) of the simple TASEP to obtain d n ˆ i (t) = n ˆ i−1 (t)(1 − n ˆ i (t)) − n ˆ i (t)(1 − n ˆ i+1 (t)) − ωD n ˆ i (t) + ωA (1 − n ˆ i (t)) . (13) dt where the first added term captures the detachment events and the later the attachment. Neglecting correlations as done before gives for the stationary state 0 = ̺i−1 (1 − ̺i ) − ̺i (1 − ̺i+1 ) − ωD ̺i + ωA (1 − ̺i ) .

(14)

The boundary conditions are not altered by LK. Using again the power series expansion (7) and keeping in mind the scaling of the on and off rates ω as in (12), we obtain the following ODE: ε 2 ∂ ̺ + ∂x ̺(2̺ − 1) − ωD ̺ + ωA (1 − ̺) = 0 . 2 x

(15)

12

Hinsch et al.

The ratio K = ωA /ωD between the attachment and detachment rates will prove an important parameter in the analysis of this differential equation. Since the case K 6= 1 is considerably complicated in its mathematical analysis we refer the reader to reference [26] and restrain our discussion to the case K = 1. In the thermodynamic limit (ε → 0) and with ωD = ωA = Ω the ODE(15) simplifies to first order: (∂x ̺ − Ω)(2̺ − 1) = 0 . (16) Obviously, there are two solutions to this general ODE problem. The homogeneous density ̺L = 1/2 given by the Langmuir isotherm and the linear slope ̺(x) = Ωx + C. The value of C is determined by the boundary conditions and leads to the two solutions ̺α (x) = α + Ωx and ̺β (x) = 1 − β − Ω + Ωx. The complete density profile ̺(x) is the combination of one or several of the three densities above. Depending on how they are matched, we distinguish several phases as explained in the following. 5.2

Phase Diagram and Density Distributions

The only area in the phase diagram that does not change compared to simple TASEP is the upper right quadrant. This is not surprising since the maximal current phase is a bulk controlled regime anyway. Therefore the additional bulk dynamics with the Langmuir isotherm at ̺L = 1/2 do not result in any changes in the density distribution. In this case, non-equilibrium and equilibrium dynamics do not compete but cooperate. As in TASEP different solutions can be matched in various ways, the simplest being the connection by a domain wall between the left solution ̺α and the right solution ̺β . Depending on the current distribution, two possibilities have to be distinguished. As both solutions are non-homogeneous, the corresponding currents jα and jβ will be strictly monotonic (Fig. 4 (left)). If the currents equal each other inside the system at a position xw , the DW is localized at this position, as a displacement to either side would result in a current inequality that drives the DW back to xw (see Fig. 4 (left)). Ergo, the TASEP/LK exhibits multi-phase existence of low and high density regions (LD-HD phase) in the stationary state on all time-scales, opposed to TASEP where this behavior is only observed for short observation times. Recently, this DW localization has been observed experimentally [28]. If the matching of left and right current is not possible inside the system, the known LD and HD phases are found. This is the case for one boundary condition being considerably larger than the other and is evidently depending quantitatively on the slope of the density solutions. This slope is determined by the ratio of the TASEP step rate and the bulk interchange rate Ω. For large Ω any density imposed by the boundary relaxes fast against the Langmuir isotherm of ̺L = 1/2 resulting a steep slope of the density profile. This fact allows for the existence of two other phases with multi-regime coexistence. We could imagine a scenario in which the boundary imposed density

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

13

solutions decay fast enough towards the isotherm to enable a three-regime coexistence of low density, maximal current and high density (LD-MC-HD phase). Furthermore, a combination of a MC phase with a boundary layer on one side and a LD or HD region of finite extend at the other boundary can be imagined. Not all these phases will be realized for every value of Ω. Instead, the phase topology of two-dimensional cuts through the α, β, Ω-phase space changes. An example is shown in Fig. 4 (right).

Fig. 4. (left) The DW connects the two densities ρα and ρβ (both dashed) and is localized at the point where the correspondent currents jα and jβ match. Note the finite extend of the DW (localization length) that is only produced for Monte-Carlo simulations (solid wiggly line) and is not captured by mean-field results (solid line). (right) Topological changes in the phase diagram of TASEP/LK for (a) Ω = 0.3, (b) Ω = 0.5, (c) Ω = 1, from [26]

5.3

Domain Wall Theory

After we have derived an analytical solution for the density profile based on the mean-field differential equation (15), we now have a closer look on the domain wall and its stochastic properties. As mentioned before, the density profile exhibits a discontinuity at xw that is actually a finite size continuous transition between the high and low density for systems of finite size. Only upon increasing the system size N → ∞ a sharp transition between the left and the right solution occurs. However, this is only due to the fact that the lattice spacing decreases to ε → 0+ . So compared to the lattice length L, usually normalized to 1, the domain wall is discontinuous, while on the length scale of lattice sites it will still have a finite extend. This

14

Hinsch et al.

extend is an intrinsic statistical feature and is usually referred to as localization length. The domain wall in its random walk behavior can be either subjected to equal rates (unbiased) as in TASEP or the rates of movement to the left/right can be different (biased). In general, the rates will not only be different, but also depend on the space variable. To begin with, we will show a way how to derive these rates by taking into account fluctuations of particle number [17,29]. Consider a situation where all events that can change the particle number (α, β, Ω) have a typical time scale that is considerably larger than that of jump processes on the lattice. In this case, the time between any entry and exit events is so long, that the system has enough time to “rearrange” (to reach a temporary steady state) in between . Then it is possible to identify jump rates ωl (x) (ωr (x)) for DW movement to the left (right) with the overall rate for entry and exit of any particle at any site. Specifically, if a particle enters the system, the DW is shifted to the left by a distance of ≈ ε/[̺β (xw ) − ̺α (xw )]. Therefore the rate for the DW to move one lattice site to the left is ωl = ωentry /[̺β (xw ) − ̺α (xw )]. If the density distribution ̺(x) is known analytically, it is possible to calculate ωentry , the overall rate of particle entrance, as the sum over all possible entrance events Z 1

ωentry = α(1 − α) +

0

dx ωA (1 − ̺(x)) .

(17)

The first term captures entrance events from the left boundary reservoir and the integral accounts for the Langmuir kinetics. The first multiplicand in both terms is the attempted rate of a jump, whereas the difference in brackets states the probability of the destination site to be vacant. Along the same lines, the exit rate is computed. As we know the analytical density distribution as ̺(x) = α + Ωx + ∆Θ(x − xw ) with the Heavyside function Θ and the DW height 3 ∆ = 1 − α − β − Ω, we can execute the integrals to obtain ωentry = α(1 − α) + Ω(β +

Ω ) + xΩ∆ . 2

(18)

Knowing these rates, one can complete the description of the domain wall as a random walker by calculating the position dependent jump rates to the left and right ωl (x) and ωr (x). The two rates constitute an effective potential displayed in Fig. 5. The DW will always be driven to that point xS in the system where ωexit (xS ) = ωentry (xS ). This position is the center of the DW probability density in a stochastic picture and the stationary DW position in a mean-field picture. It is in agreement with mean field results and yields xS =

Ω−α+β . 2Ω

(19)

The other quantity of interest which - in contrast to the DW position - cannot be determined by mean-field calculations is the localization length. In order to 3

The term −1Ω accounts for the diminished height that is caused by the Langmuir kinetics on the whole length 1 of the system.

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

ω

15

|∆ω| ωr

|∆ω| 0

ωl 0

x

0

xS

1

Fig. 5. (left) The DW will be localized at the point where the jump rate to the left ωl and right ωr intersect. This can be interpreted as an potential |∆ω|. (right) Variance of DW probability distribution over Ω for α = β = 0.01 and α = β = 0.1 in a system of N = 500, predictions according to [30] (solid) and based on particle number fluctuations (dashed) compared to MC simulations (dots)

compute this quantity, fluctuations have to be taken into account as we will show in the following. If we use the notation p(x) for the probability that the domain wall is at position x, then the condition for a stationary DW reads in the continuum limit ωr (x)p(x) = ωl (x + ε)p(x + ε) .

(20)

Introducing now y(x) = ωl (x)p(x) and approximate y ′ (x) = |y(x + ε) − y(x)|/ε we obtain: ωr (x) y ′ (x) + N y(x)(1 − )=0. (21) ωl (x) The solution then is given by p(x) =

1 p˜(x) = exp[−N Z Zωl (x)

Z

x

x0

dx′ (1 −

ωr (x) )] . ωl (x)

(22)

where Z accounts for normalization. In general, Z is not available explicitly, but it has been shown [30] that the unnormalized probability function can be approximated by a Gaussian 2

p˜(x) ∝ e−C(x−xS) ,

(23)

where C is given by the second order derivative of the exponent in (22) as C=

1 d2 [N 2 dx2

Z

x

x0

ωr (x) N (ωl − ωr )′ (xS ) ]= dx′ 1 − . ωl (x) 2ωl (xS )

(24)

p Hence, the variance σ = 1/(2C) of the domain wall can be easily be obtained provided that the jump rates ωl (x) and ωr (x) are available. Evans et al. [30] have assumed those rates to be ωl,r (x) = jα,β /(̺β − ̺α). Kouyos has shown [31] that using the rates (17) derived from the fluctuations of particle number, one

16

Hinsch et al.

arrives at more accurate results compared to Monte Carlo simulations (see Fig. 5). In this case C evaluates to 2N Ω∆ . (25) α(1 − α) + β(1 − β) + Ω √ As the width of the DW is given by σ = 1/ 2C, the localization of the DW scales with N −1/2 . C=

6

Conclusions and Outlook

Much in the same way as MacDonald’s pioneering paper [4] on mRNA translation, recent work on kinesin motors walking along microtubules has spurred progress in nonequilibrium transport phenomena. The ubiquitous exchange of material between the cytoplasm and the molecular track, which originally was thought to only lead to quantitative modifications of the dynamics [24], has recently been identified [25] as the source for qualitatively new phenomena such as phase separation. This introduced a completely new class of non-equilibrium transport models. These lattice gas models are characterized by a scaling of the on-off rates with system size which enables competition of driven motion and equilibrium Langmuir kinetics. Hereby, a finite residence time on the lattice ensures cooperative effects to establish multi-phase coexistence, localized shocks and an enriched phase behavior compared to prior TASEP results. There are now various routes along which one could proceed. The first one is to add more realistic features of the molecular motors, such as the fact that they are dimers [32] or that there is more than one chemo-mechanical state [28]. Such investigations are crucial for a quantitative understanding of intracellular traffic in various ways. One might ask how robust the features of minimal models are with respect to the addition of more molecular details. In the case of dimers the answer is far from obvious since the non-equilibrium dynamics of dimer adsorption shows rich dynamic behavior with anomalously slow relaxation towards the equilibrium state [27]. How this combines with the driven transport along the molecular track was recently analyzed thoroughly [32]. While correlation effects due to the extended nature of dimers invalidate a simple mean-field picture it was found that an extended mean-field scheme can be developed which quantitatively describes the stationary phases. Surprisingly, the topology of the phase diagram and the nature of the phases is similar to the minimal model with monomers. The physical origin of this robustness can be traced back to the form of the current-density relation which exhibits only a single maximum. The second line of research generalizing the minimal model [25] asks for the effect of interactions, more than one molecular traffic lane, “road blocks” such as microtubule associated proteins and various other kinds of “disorder”, bidirectional traffic, coupling of driven and diffusive transport and the like on the stationary density profiles and the dynamics. In almost every instance it is found that this leads to an even richer behavior with new phenomena emerging.

From Intracellular Traffic to a Novel Class of Driven Lattice Gas Models

17

For TASEP it is known that isolated defects (slow sites) depending on their strength may either give rise to local density perturbations for low particle densities or yield macroscopic effects for densities close to the carrying capacity [33]. The interplay between coupling to the motor reservoir in the cytoplasm and the fluctuations in the capacity limit due to disorder along the track gives rise to a number of interesting collective effects [34]. Coupling two lanes by allowing particle exchange at a constant rate along the molecular track also results in novel phenomena. Similar to equilibrium phase transitions described by field theories with two coupled order parameters higher order critical points may emerge [35]. Coupling diffusive and driven transport, the origin of new phenomena is due to a competition of different processes of comparable time scales. The qualitative failure of mean-field theory in some of these systems [36] comes as quite a surprise, since mean-field has proven to predict phase diagrams for large systems with an astonishing accuracy in the lattice gas models mentioned above.

References 1. A.D. Mehta, M. Rief, J.A. Spudich, D.A. Smith, R.M. Simmons: Science 283, 1689 (1999) 2. J. Howard: Mechanics of Motor Proteins and the Cytoskeleton (Sinauer, Sunderland, 2001) 3. M. Schliwa and G. Woehlke: Nature 422, 759 (2003) 4. C.T. MacDonald, J.H. Gibbs, A.C.Pipkin: Biopolymers 6, 1 (1968) 5. N.G. van Kampen: Stochastic Processes in Physics and Chemistry (North Holland, Amsterdam, 1981) 6. A. Parmeggiani, C.F. Schmidt: in Function and Regulation of Cellular Systems (Birkh¨ auser, 2004) 7. F. J¨ ulicher, A. Ajdari, J. Prost: Rev. Mod. Phys. 69, 1269 (1997) 8. N.D. Mermin, H. Wagner: Phys. Rev. Lett. 17, 1133 (1966) 9. B. Derrida, M.R. Evans: J. Physique I 3, 311 (1993) 10. G. Sch¨ utz: Phys. Rev. E 47, 4265 (1993) 11. M. Dudzinski, G.M. Sch¨ utz: J. Phys. A 33, 8351 (2000) 12. H. van Beijeren, R. Kutner, H. Spohn: Phys. Rev. Lett. 54, 2026 (1985) 13. B. Derrida, E. Domany, D. Mukamel: J. Stat. Phys. 69, 667 (1992) 14. G. Sch¨ utz: in Phase Transitions and Critical Phenomena vol 19 edited by C. Domb, J. Lebowitz (Academic, London 2001) 15. N. Goldenfeld: Lectures on phase transitions and the renormalization group (Perseus, Reading, 1992) 16. G.D. Mahan: Many Particle Physics (Kluwer, 2000) 17. A.B. Kolomeisky, G.M. Sch¨ utz, E.B. Kolomeisky, J.P. Straley: J. Phys. A 31, 6911 (1998) 18. L.B. Shaw, P.K.P. Zia, K.H. Lee: Phys. Rev. E 68, 021910 (2003) 19. R. Heinrich, T.A. Rapoport: J. Theor. Bio. 86, 279 (1980) 20. T. Chou: Biophys. J 86, 755 (2003) 21. P. Popkov, G. M. Sch¨ utz: Europhys. Lett. 48, 257 (1999) 22. M.R. Evans, D.P. Foster, C Godreche, D Mukamel: Phys. Rev. Lett. 74, 208 (1995) 23. V. Popkov, I. Peschel: Phys. Rev. E 64, 026126 (2001)

18 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36.

Hinsch et al. R. Lipowsky, S. Klumpp, T.M. Nieuwenhuizen: Phys. Rev. Lett. 87, 108101 (2001) A. Parmeggiani, T. Franosch, E. Frey: Phys. Rev. Lett. 90, 86601 (2003) A. Parmeggiani, T. Franosch, E. Frey: Phys. Rev. E 70, 46101 (2004) E. Frey, A. Vilfan: Chem. Phys. 284, 287 (2002) K. Nishinari, Y. Okada, A. Schadschneider, D. Chowdhury: Phys. Rev. Lett 95, 118101 (2005) L. Santen, C. Appert: J. Stat. Phys 106, 187 (2002) M.R. Evans, R. Juh´ asz, L. Santen: Phys. Rev. E 68, 026117 (2003) R. Kouyos: Disorder and Phase Separation in the Total Asymmetric Exclusion Process, Diploma Thesis, ETH Z¨ urich (2004) P. Pierobon, T. Franosch, E. Frey: Phys. Rev. E 74, 031920 (2006) G. Tripathy and M. Barma: Phys. Rev. E 58, 1911 (1998) P. Pierobon, M. Mobilia, R. Kouyos, and E. Frey: Phys. Rev. E 74, 031906 (2006) T. Reichenbach, T. Franosch, and E. Frey: Phys. Rev. Lett. 97, 050603 (2006) H. Hinsch, E. Frey: Phys. Rev. Lett. 97, 095701 (2006)