text:
# Traffic of molecular motors through tube-like compartments

## I Introduction

Molecular motors are proteins that transform the free energy released from chemical reactions into mechanical work. In this article we consider a special class of motor proteins, namely cytoskeletal motors which perform directed walks along cytoskeletal filaments as reviewed in . In the cell, these motors have different functions related, e.g., to vesicle transport and cell division. The best studied examples are kinesins, which walk along microtubules, and certain types of myosins, which walk along actin filaments. After a certain walking time, such a motor unbinds from its filament because its binding energy is finite and can be overcome by thermal activation. For kinesins, this typically happens after 100–150 steps or after a walking time of about 1.2–1.8 seconds, see, e.g., . In many motility assays, the filaments are immobilized on a substrate and are in contact with an aqueous solution. In such a situation, the unbound motors diffuse in the surrounding fluid until they eventually reattach to the same or another filament.

Recently, we have introduced lattice models to study the motors’ random walks, which consist of many diffusional encounters between motors and filaments in open and closed compartments . If many motors are placed in such a compartment, hard core exclusion between the motors has to be taken into account, since the motors are strongly attracted to the binding sites of the filaments, so that the filaments get overcrowded. These models are new variants of driven lattice gas models and exclusion processes, for which the active processes which drive the particles are localized to the filaments.

Lattice models of driven diffusive systems have been studied extensively in the last years, see e.g. . The simplest model is the asymmetric simple exclusion process (ASEP) in one dimension, where particles hop on a one-dimensional lattice with a strong bias towards one direction (in the simplest case there are no backward steps at all) and the only interaction of the particles is hard core exclusion, i.e., steps to occupied lattice sites are forbidden. When coupled to open boundaries, this simple model already exhibits a complex phase diagram, see e.g. , which we will review below in some detail. The first model for the 1–dimensional ASEP was introduced more than 30 years ago by MacDonald et al.  in the context of protein synthesis by ribosomes on messenger RNA (mRNA). At that time it was solved using a mean field approach and used to explain results of radioactive labeling experiments  which showed that protein synthesis gets slower as the ribosome moves on the mRNA template. The model of MacDonald et al. explained this by the steric hindrance between successive ribosomes along the mRNA track. Two years later the same model was discussed by Spitzer as a simple example for interacting particles in probability theory . Since then, the asymmetric simple exclusion process (and variants) has been studied extensively as a generic model for non-equilibrium phase transitions  and interacting stochastic systems  as well as in other applications such as traffic flow . Many properties of the 1–dimensional ASEP are known exactly.

As mentioned, the lattice models for random walks of molecular motors differs from the driven lattice gas models in an important way: Walks of molecular motors are only ’driven’ as long as the motor is attached to a cytoskeletal filament, hence ’driving’ is localized to one or several lines. It can therefore be viewed as an ASEP which has the additional property, that particles (molecular motors) can unbind from the track with a small probability, diffuse in the surroundings and reattach to the same or another filament. In more mathematical terms, the ASEP is coupled to a symmetric exclusion process via adsorption and desorption of particles onto filaments.

Boundary conditions play an important role in driven systems. This becomes apparent, e.g., if one compares a tube–like system with periodic boundary conditions with one with closed boundaries. In the system with closed boundaries, a traffic jam of motors arises at one end of the system and the current of motors bound to the filament is balanced by diffusive currents of unbound motors as first shown in . With periodic boundary conditions, motors arriving at the right end of the system just restart their walk from the left end and a net current through the systems is obtained.

In this article, we study the stationary states of tube-like compartments with open boundaries. These compartments have the shape of a cylinder and contain one filament which is placed along the cylinder axis in order to obtain the simplest possible geometry. The bound motors move along the filament and the unbound motors diffuse within the cylinder, see Fig. 1. At the ends of the tube, motors are inserted and extracted. Such a system is accessible to in vitro experiments using standard motility assays, but it can also be viewed as a strongly simplified model for motor–based transport in an axon , if these motors, which are synthesized in the cell body, are at least partly degraded at the axon terminal , a situation that can be mimicked by insertion and extraction of motors at the ends of a tube. The stationary states depend strongly on the way, in which the motors are inserted and extracted at the boundaries as we will explicitly demonstrate for three different types of boundary conditions, see Fig. 2.

Our article is organized as follows. After introducing the model in section II, we start in section III with periodic boundary conditions. This case can be solved exactly, since it satisfies local balance of currents in the radial direction, see appendix A. In section IV, we discuss the situation in which the density of bound motors on the filament is fixed at the boundaries. Finally, in section V, we consider the case where the filament is shorter than the tube and the motors diffuse into and out of the tube. The main tool to study the open systems are Monte Carlo simulations. These are supplemented by dynamical considerations and self-consistent or mean field calculations. Some details of the latter calculations are presented in appendix B.

## II Theoretical modelling

### A Tube geometry

We consider the motion of molecular motors in a cylindrical tube as shown in Fig. 1. The tube has length $`L`$ and radius $`R`$. The total number of motors within the tube, denoted by $`N_{\mathrm{mo}}`$, defines the overall motor concentration

$$\rho _{\mathrm{mo}}\frac{N_{\mathrm{mo}}}{\pi R^2L}.$$
(1)
Note that this concentration corresponds to the particle number density of the motors and, thus, is insensitive to the size of the motor particles. In general, the latter size depends on the type of motor and on the type of cargo attached to it. In the following, we will implicitly assume that the motor particle has a linear size which is comparable to the basic length scale $`\mathrm{}`$ as defined further below. If one wants to study the dependence on the motor particle size in a systematic way, one should measure the overall motor concentration in terms of the volume fraction of the motor particles.

The cylindrical tube contains one filament located along its symmetry axis which is taken to be the x–axis; the two other Cartesian coordinates are denoted by $`y`$ and $`z`$. Motors bound to the filament undergo directed motion, while unbound motors diffuse freely. Since the motors are strongly attracted by the filament, a large fraction of these motors is in the bound state and mutual exclusion from the binding sites on the filament has to be taken into account even for relatively small overall motor concentrations.

In order to include this mutual exclusion (or hard core repulsion) in the theoretical description, we map the system onto a lattice gas model on a simple cubic lattice. The lattice is oriented in such a way that its three primitive vectors point parallel to the $`x`$–, $`y`$–, and $`z`$–axis, respectively. A rather natural choice for the lattice parameter $`\mathrm{}`$ is the repeat distance of the filament, which is $`8`$nm in the case of kinesin motors moving on microtubule and $`36`$nm for myosin V motors moving on actin filaments.

The discretized tube consists of one line of binding sites, which represents the filament, and $`N_{\mathrm{ch}}`$ unbound ’channels’, i.e. lines of lattice sites parallel to the filament. Thus the cross section $`\varphi `$ of the tube is equal to

$$\varphi =(1+N_{\mathrm{ch}})\mathrm{}^2.$$
(2)
For sufficiently large radii, one has $`\varphi \pi R^2`$, while for small radii, there are corrections due to the underlying lattice.

In the following, we will measure all distances such as $`L`$ and $`R`$ in units of $`\mathrm{}`$. Thus, for the simulations, both $`L`$ and $`R`$ will be quoted as integers. The integer value of $`R`$ corresponds to a certain number of lattice sites along the Cartesian coordinates $`y`$ and $`z`$ which are perpendicular to $`x`$ and run parallel to two basis vectors of the simple cubic lattice. Thus, the interior of the tube contains all channels with $`y^2+z^2R^2`$.

### B Random walks with mutual exclusion

At each time step, a bound motor attempts to make a forward or backward step and to jump to the next lattice site to its right or to its left with probability $`\alpha `$ or $`\beta `$. In addition, the motor attempts to jump to each of the neighboring sites away from the filament with probability $`ϵ/6`$, and does not attempt to jump at all, i.e., to rest at the filament site, with probability $`\gamma `$. Since the sum of these probabilities is equal to one, we have $`\gamma =1\alpha \beta 2ϵ/3`$. The velocity $`v_\mathrm{b}`$ in the bound state is given by $`v_\mathrm{b}=(\alpha \beta )\mathrm{}/\tau `$, where $`\tau `$ is the basic time scale of these random walks. In the following, we will measure all times and rates in units of $`\tau `$ and $`\tau ^1`$, respectively. Since backward steps are rare for cytoskeletal motors, we will focus on the case $`\beta =0`$, i.e., we will ignore backward steps.

If the particle unbinds from the filament, it attempts to jump to all nearest neighbor sites of the simple cubic lattice with equal probability $`1/6`$. By choosing the time scale $`\tau `$ as $`\mathrm{}^2/D_{\mathrm{ub}}`$, this hopping probability can be made to fit the diffusion coefficient of unbound motors, $`D_{\mathrm{ub}}`$. When measured in units of $`\mathrm{}`$ and $`\tau `$, the dimensionless diffusion coefficient $`D_{\mathrm{ub}}=1/6`$. In principle, the resting probability $`\gamma `$ can then be used to account for the ratio $`D_{\mathrm{ub}}/(v_\mathrm{b}\mathrm{})`$, which is quite large in the case of kinesin, and to adapt the velocity $`v_\mathrm{b}`$ to the values obtained from experiments . For simplicity, we will often choose $`\gamma =0`$ in order to eliminate one parameter from the problem. If an unbound motor attempts to hop to a filament site, the motor binds to it with sticking probability $`\pi _{\mathrm{ad}}`$, while the step (and hence binding) is rejected with probability $`1\pi _{\mathrm{ad}}`$.

Both in the bound and in the unbound state, hopping attempts can only be successful, if the target site is not occupied by another motor; otherwise the particle must stay where it is. In the following, we will mainly study overall motor concentrations $`\rho _{\mathrm{mo}}`$ in the range $`0\rho _{\mathrm{mo}}\begin{array}{c}<\\ \end{array}0.05`$. For such small values of $`\rho _{\mathrm{mo}}`$, mutual exclusion of the unbound motors can be safely ignored. However, because the motors are strongly attracted to the filament, the concentration of bound motors is much larger and it is crucial to take their mutual exclusion into account.

## III Periodic boundary conditions

First, we consider a cylindrical tube with periodic boundary conditions in the longitudinal direction. Because of the translational invariance in the direction parallel to the filament, there are no net radial currents in the stationary state. Indeed, a non-zero radial current in a state, which is translationally invariant in the longitudinal direction, would lead to net radial transport of motor particles, which is incompatible with the reflecting radial boundaries. This means that there is a bound current $`j_\mathrm{b}`$ on the filament, but both currents of motors binding to and unbinding from the filament and radial currents of unbound motors  are balanced locally. We call this situation radial detailed balance or radial equilibrium. If there is only one unbound channel, $`N_{\mathrm{ch}}=1`$, (or $`N_{\mathrm{ch}}`$ equivalent channels) radial equilibrium is equivalent to adsorption equilibrium. It is clear, that this will no longer be true, if translational invariance is broken by boundaries or blocked sites on the filament. Another important property of systems with periodic boundary conditions is that, in this case, the number $`N_{\mathrm{mo}}`$ of motors in the system is conserved, which does not apply to open systems.

Because of translational invariance along the x–axis, the bound and unbound motor densities, $`\rho _\mathrm{b}`$ and $`\rho _{\mathrm{ub}}`$, do not depend on $`x`$. In addition, it follows from the absence of radial currents that $`\rho _{\mathrm{ub}}`$ is also independent of the radial coordinate $`r=(y^2+z^2)^{1/2}`$. Hence $`\rho _\mathrm{b}`$ and $`\rho _{\mathrm{ub}}`$ are constant, and radial equilibrium implies the relation

$$\frac{ϵ}{6}\rho _\mathrm{b}(1\rho _{\mathrm{ub}})=\frac{\pi _{\mathrm{ad}}}{6}\rho _{\mathrm{ub}}(1\rho _\mathrm{b})$$
(3)
which leads to

$$\rho _\mathrm{b}=\frac{\rho _{\mathrm{ub}}}{ϵ/\pi _{\mathrm{ad}}+(1ϵ/\pi _{\mathrm{ad}})\rho _{\mathrm{ub}}}.$$
(4)
Here and below, the densities $`\rho _\mathrm{b}`$ and $`\rho _{\mathrm{ub}}`$ are local particle number densities which satisfy

$$0\rho _\mathrm{b}1\mathrm{and}0\rho _{\mathrm{ub}}1.$$
(5)
In dimensionful units, this corresponds to $`0\rho _\mathrm{b}1/\mathrm{}^3`$ and $`0\rho _{\mathrm{ub}}1/\mathrm{}^3`$.

Since the total number $`N_{\mathrm{mo}}`$ of motors is conserved, we can use the normalization condition,

$$\rho _\mathrm{b}+N_{\mathrm{ch}}\rho _{\mathrm{ub}}=\frac{N_{\mathrm{mo}}}{L},$$
(6)
to obtain a quadratic equation for the densities as a function of the system size and the total number of motors. Since one root is always negative, the physically meaningful solution is:

$`\rho _{\mathrm{ub}}={\displaystyle \frac{1}{2N_{\mathrm{ch}}(1\frac{ϵ}{\pi _{\mathrm{ad}}})}}[1N_{\mathrm{ch}}{\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}}+{\displaystyle \frac{N_{\mathrm{mo}}}{L}}(1{\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}})`$ (8)
$`+\sqrt{\left(1N_{\mathrm{ch}}{\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}}+{\displaystyle \frac{N_{\mathrm{mo}}}{L}}(1{\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}})\right)^2+4N_{\mathrm{ch}}(1{\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}}){\displaystyle \frac{ϵ}{\pi _{\mathrm{ad}}}}{\displaystyle \frac{N_{\mathrm{mo}}}{L}}}].`$
From this expression for the unbound density, the bound density $`\rho _\mathrm{b}`$ follows via Eq. (4) and the stationary current is given by $`J=j_\mathrm{b}=v_\mathrm{b}\rho _\mathrm{b}(1\rho _\mathrm{b})`$. The current calculated in this way is shown in Fig. 3 as a function of $`N_{\mathrm{mo}}/L`$. The data points (circles) are the results of Monte Carlo simulations for a system of length $`L=200`$ and radius $`R=25`$ and are in very good agreement with the analytical solution.

It follows from the analytical solution that the current $`J=j_\mathrm{b}`$ vanishes at $`N_{\mathrm{mo}}/L=1+N_{\mathrm{ch}}`$, behaves as

$$\frac{J}{v_\mathrm{b}}\frac{1}{1+\frac{ϵ}{\pi _{\mathrm{ad}}}N_{\mathrm{ch}}}\frac{N_{\mathrm{mo}}}{L}$$
(9)
for small $`N_{\mathrm{mo}}/L`$, and has the maximal value $`\mathrm{max}(J/v_\mathrm{b})=1/4`$ for $`\rho _\mathrm{b}=1/2`$, $`\rho _{\mathrm{ub}}=\frac{ϵ}{\pi _{\mathrm{ad}}}/(1+\frac{ϵ}{\pi _{\mathrm{ad}}})`$ and $`N_{\mathrm{mo}}/L=1/2+N_{\mathrm{ch}}\frac{ϵ}{\pi _{\mathrm{ad}}}/(1+\frac{ϵ}{\pi _{\mathrm{ad}}})`$.

Let us add two remarks:

(i) Eq. (3), as stated here, can be considered as a mean field equation. However, using the quantum Hamiltonian representation of the stochastic process, it can be shown to hold exactly. The calculation is simple, but rather technical and is therefore presented in Appendix A, which shows that stationary states can be constructed as product measures provided the bound and unbound densities satisfy Eq. (3).

(ii) Note that in contrast to a homogeneously driven lattice gas such as the asymmetric simple exclusion process, there is no particle hole symmetry here. Particles attempt to leave the filament with rate $`ϵ/6`$ to a neighboring site while holes do so with rate $`\pi _{\mathrm{ad}}/6`$, i.e. particles are strongly attracted by the filament, while holes are not. However, if one considers only the bound density, the current density relationship $`J=j_\mathrm{b}=v_\mathrm{b}\rho _\mathrm{b}(1\rho _\mathrm{b})`$ is invariant under the exchange of particles and holes. As we will see in the next section, this can lead to an apparent particle hole symmetry for systems with radial equilibrium, because the radial currents vanish and the state of the system can be determined by the bound density alone.

## IV Open boundaries with radial equilibrium

Now, let us consider the more interesting case, where the tube is open and the densities at the left and right boundary are fixed. To be precise, we consider two different sets of boundary conditions, (B) and (C) as shown schematically in Fig. 2. In this section, we study case (B) while case (C) will be considered in the next section V.

For case (B), we add two layers of boundary sites at $`x=0`$ and $`x=L+1`$ with $`y^2+z^2R^2`$. As before, the filament is located at $`y=z=0`$. We then fix the density on the additional filament sites according to

$$\rho _\mathrm{b}(x=0)\rho _{\mathrm{b},\mathrm{in}}\mathrm{and}\rho _\mathrm{b}(x=L+1)\rho _{\mathrm{b},\mathrm{ex}},$$
(10)
see Fig. 2. Furthermore, the densities on the nonfilament boundary sites are chosen in such a way that radial equilibrium as given by (3) holds at both boundaries.

These boundary conditions are implemented by the following choice of random walk probabilities. First, we eliminate two parameters from the problem, namely the jump probability $`\beta `$ to make backward steps on the filament and the resting probability $`\gamma `$ to make no step at all on the filament. Thus, we take $`\beta =\gamma =0`$ throughout this section.

Next, when we choose a site within the left or right boundary layer during the Monte Carlo sweep, we first draw a random number $`\omega `$ which is uniformly distributed over the interval $`0<\omega 1`$. The chosen left or right boundary site is taken to be occupied if $`\omega \rho _{\mathrm{b},\mathrm{in}}`$ or $`\omega \rho _{\mathrm{b},\mathrm{ex}}`$, respectively. If a boundary site with $`y^2+z^2>0`$ is occupied by a motor particle, this particle attempts to jump into the tube with probability $`1/6`$. If the left boundary site with $`y=z=0`$ is occupied, the corresponding particle attempts to jump onto the left end of the filament with probability $`(12ϵ/3)`$. If the right boundary site with $`y=z=0`$ is occupied, this particle cannot enter since $`\beta =0`$. In addition, all particles which jump from a site within the tube, i.e., from a lattice site with $`1xL`$, onto a boundary site, are extracted from the tube.

### A Phase diagram

In the limiting case with $`N_{\mathrm{ch}}=0`$ and $`ϵ=0`$, our system becomes equivalent to the 1–dimensional ASEP, for which the density profiles and the phase diagram are known exactly . Let us therefore summarize some of the known properties of this process, for which we use symbols without the subscript ’$`\mathrm{b}`$’.

For the 1–dimensional ASEP, there are three different phases. If the density $`\rho _{\mathrm{in}}`$ at the left boundary is small and satisfies $`\rho _{\mathrm{in}}<1/2`$ and if the density $`\rho _{\mathrm{ex}}`$ at the right boundary is not too large with $`\rho _{\mathrm{ex}}<1\rho _{\mathrm{in}}`$, the system is in the low density (LD) phase, for which the bulk density $`\rho ^0`$ is equal to the left boundary density, and the current is given by $`J=v\rho _{\mathrm{in}}(1\rho _{\mathrm{in}})`$. Because of the particle hole symmetry, an analogous situation holds for $`\rho _{\mathrm{ex}}>1/2`$ and $`\rho _{\mathrm{in}}>1\rho _{\mathrm{ex}}`$ \[high density (HD) phase\]. Now, the bulk density is given by $`\rho ^0=\rho _{\mathrm{ex}}`$ and the current is $`J=v\rho _{\mathrm{ex}}(1\rho _{\mathrm{ex}})`$. At the line $`\rho _{\mathrm{ex}}=1\rho _{\mathrm{in}}<1/2`$, a discontinuous phase transition takes place, and the bulk density jumps from the left boundary value to the right boundary value. Finally, for $`\rho _{\mathrm{in}}>1/2`$ and $`\rho _{\mathrm{ex}}<1/2`$, the bulk density $`\rho ^0=1/2`$ and the current attains its maximal value $`J=v/4`$. Therefore, this phase is called maximal current (MC) phase. The phase transition towards the maximal current phase is continuous with a diverging correlation length.

The formation of three different phases can be understood in terms of the underlying dynamics of domain walls and density fluctuations . In the low density and high density phases, the selection of the stationary state is governed by domain wall motion. Thus consider a domain wall, which forms between regions of different densities. If its velocity is positive, the domain wall travels to the right boundary and the bulk density is equal to the left boundary density (low density phase). Likewise, the bulk density is given by the right boundary density, if the domain wall velocity is negative (high density phase). At the transition line with $`\rho _{\mathrm{ex}}=1\rho _{\mathrm{in}}<1/2`$ the domain wall velocity is zero, and domain walls diffuse through the system.

The second mechanism is related to density fluctuations. If a small perturbation of the density, corresponding, e.g., to some added particles, enters the system at the left boundary, it can move with positive or negative velocity, i.e., it can spread into the bulk or is driven back towards the boundary. In the maximal current phase the velocity of such density perturbations, is negative and density fluctuations coming from the left boundary are driven back to the boundary. Hence, increasing $`\rho _{\mathrm{in}}`$ does not increase the bulk density, since additional particles cannot enter the system (overfeeding effect). Both the velocities of domain walls and of density fluctuations are governed by the same density current relation which is $`j=v\rho (1\rho )`$ for the 1–dimensional ASEP .

For the filament in a tube as considered here, the velocities of domain walls and density fluctuations are similar to those for the ASEP in one dimension. This can be understood from the fact that the density–current relationship on the filament is the same as for the 1–dimensional ASEP provided one rescales all currents by the factor $`(12ϵ/3)`$ which arises from the possibility to unbind from the filament.

Thus, let us first consider the case, where the behavior on the filament is determined by domain walls, i.e. the low density and high density phases. Drift motion of domain walls is governed by the domain wall velocity on the filament, which is $`v_s=v(1\rho _{\mathrm{ex}}\rho _{\mathrm{in}})`$ for the one-dimensional ASEP . In the tube system considere here, the domain wall velocity is slowed down compared to this value, because the domain wall of the unbound density must follow the bound density domain wall by binding and unbinding of motors to and from the filament. However, the sign of the domain wall velocity is the same for this tube system as for the one–dimensional ASEP, and the domain wall velocity changes sign at the same values of the boundary densities. An explicit expression for $`v_s`$ can be obtained from the general expression for the domain wall velocity as given in Ref.  by integrating the density over the tube cross-section, which leads to

$$v_s=\frac{v_\mathrm{b}\rho _{\mathrm{b},\mathrm{ex}}(1\rho _{\mathrm{b},\mathrm{ex}})v_\mathrm{b}\rho _{\mathrm{b},\mathrm{in}}(1\rho _{\mathrm{b},\mathrm{in}})}{\rho _{\mathrm{b},\mathrm{ex}}\rho _{\mathrm{b},\mathrm{in}}+N_{\mathrm{ch}}\rho _{\mathrm{ub},\mathrm{ex}}N_{\mathrm{ch}}\rho _{\mathrm{ub},\mathrm{in}}}$$
(11)
for the geometry considered here. Remember that the unbound densities are related to the bound densities by radial equilibrium at the boundaries. If a domain wall spreads from the left or right boundary into the system, radial equilibrium will hold approximately in the bulk because of translational invariance, but, in addition, all the way down to the dominating boundary, where we have imposed radial equilibrium via the boundary conditions. Therefore, the current and the bulk density are determined by the bound density at the boundaries.

The drift velocity of fluctuations of the bound density is $`v_c=12\rho _\mathrm{b}`$ . If the behavior on the filament is determined by this velocity, i.e. in the case $`\rho _{\mathrm{b},\mathrm{in}}>1/2`$ and $`\rho _{\mathrm{b},\mathrm{ex}}<1/2`$, radial equilibrium will not hold up to the boundary and unbound motors can enter the tube. However, after a short distance they will bind to the filament and act like an additional particle in the maximal current phase: The density fluctuation generated by the additional particle moves towards the boundary. Therefore in the bulk, the bound density is 1/2 and the current is $`v_\mathrm{b}/4`$.

Summarizing these considerations, we predict to find exactly the same phase diagram as for the 1–dimensional ASEP, see Fig. 4. In fact, we have chosen our boundary conditions (B) in such a way that we only have to replace the density $`\rho `$ of the 1–dimensional ASEP by the bound density $`\rho _\mathrm{b}`$ and the boundary densities $`\rho _{\mathrm{in}}`$ and $`\rho _{\mathrm{ex}}`$ by the boundary densities on the filament, $`\rho _{\mathrm{b},\mathrm{in}}`$ and $`\rho _{\mathrm{b},\mathrm{ex}}`$. The unbound density in the bulk is obtained by radial equilibrium from the bound density. These expectations are confirmed (i) by a detailed analysis of the discrete mean field equations and (ii) by extended Monte Carlo simulations.

A remarkable feature of the phase diagram is that it seems to exhibit particle hole symmetry, while the dynamics does not. The signature of particle hole symmetry in the phase diagram is the symmetry between the high density and low density phase. In both these phases, the bulk density is approximately constant and radial equilibrium holds approximately in the whole system except close to the left or right boundary. Radial equilibrium holds exactly at that boundary which determines the bulk behavior. Therefore radial currents vanish on average and the phase diagram is determined by the bound density alone, resulting in a phase diagram which gives the impression of particle hole symmetry, although this symmetry is broken. Indeed, even though the substitution of the bound density $`\rho _\mathrm{b}`$ by $`1\rho _\mathrm{b}`$ leads to the same bound current $`j_\mathrm{b}`$, it does not lead to the unbound density $`1\rho _{\mathrm{ub}}`$. Two stationary states with bound densities that are related through particle hole symmetry are characterized by two unbound densities which are both smaller than the bound ones and thus break the apparent symmetry.

### B Density profiles

To discuss the concentration profiles of the bound and unbound motors, we use continuum mean field equations and compare the mean field results with simulations. In addition, we make a two-state approximation, i.e., we consider the case of a single unbound channel, so that a motor can be in only two states, bound to the filament or unbound. The two-state model is exact for an arbitrary number $`N_{\mathrm{ch}}`$ of *equivalent* unbound channels, but it can also serve as an approximation for the original tube systems: Since the unbound density depends only weakly on the radial coordinate, we will consider the approximation in which $`\rho _{\mathrm{ub}}`$ is taken to be independent of $`r`$ and the unbound channels are taken to be equivalent. This approximation corresponds to a two–state model in which the bound motors are described by the density $`\rho _\mathrm{b}(x)`$ and the unbound motors by the density $`\rho _{\mathrm{ub}}(x)`$. The effective diffusion coefficient for the unbound motors is given by $`N_{\mathrm{ch}}D_{\mathrm{ub}}`$.

Using again the mean field approximation, the total current $`J`$ parallel to the tube axis is now given by

$$J=v_\mathrm{b}\rho _\mathrm{b}(1\rho _\mathrm{b})D_\mathrm{b}\frac{}{x}\rho _\mathrm{b}N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{}{x}\rho _{\mathrm{ub}}$$
(12)
and the equality between incoming and outgoing currents at any lattice site leads to

$$\frac{}{x}\left[v_\mathrm{b}\rho _\mathrm{b}(1\rho _\mathrm{b})D_\mathrm{b}\frac{}{x}\rho _\mathrm{b}\right]=\stackrel{~}{\pi }_{\mathrm{ad}}\rho _{\mathrm{ub}}(1\rho _\mathrm{b})\stackrel{~}{ϵ}\rho _\mathrm{b}(1\rho _{\mathrm{ub}}),$$
(13)
where $`D_\mathrm{b}`$ and $`D_{\mathrm{ub}}`$ are the diffusion coefficients of motors in the bound and unbound state, respectively. In addition, we have introduced the rescaled binding and unbinding rates, $`\stackrel{~}{\pi }_{\mathrm{ad}}=\frac{2}{3}\pi _{\mathrm{ad}}`$ and $`\stackrel{~}{ϵ}=\frac{2}{3}ϵ`$.

In order to determine the density profiles far from the boundaries, we first calculate the homogeneous, $`x`$–independent solutions, $`\rho _\mathrm{b}^0`$ and $`\rho _{\mathrm{ub}}^0`$, of the two mean field equations (12) and (13). The first equation (12) for the total current $`J`$ then reduces to the current–density relationship

$$J=v_\mathrm{b}\rho _\mathrm{b}^0(1\rho _\mathrm{b}^0)$$
(14)
whereas the second equation (13) becomes

$$\stackrel{~}{\pi }_{\mathrm{ad}}\rho _{\mathrm{ub}}^0(1\rho _\mathrm{b}^0)=\stackrel{~}{ϵ}\rho _\mathrm{b}^0(1\rho _{\mathrm{ub}}^0)$$
(15)
which implies radial equilibrium for the $`x`$–independent solutions. We then decompose the densities according to

$$\rho _\mathrm{b}(x)=\rho _\mathrm{b}^0+\eta _\mathrm{b}(x)\mathrm{and}\rho _{\mathrm{ub}}(x)=\rho _{\mathrm{ub}}^0+\eta _{\mathrm{ub}}(x),$$
(16)
and expand the mean field equations (12) and (13) in powers of the density deviations $`\eta _\mathrm{b}`$ and $`\eta _{\mathrm{ub}}`$. The details of this expansion are described in Appendix B.

#### 1 Low density and high density phases

For the high and low density phases with $`\rho _\mathrm{b}^01/2`$, the expansion of the mean field equations up to first order in $`\eta _\mathrm{b}`$ and $`\eta _{\mathrm{ub}}`$ leads to an exponential approach $`\mathrm{exp}(x/\xi )`$ of the density profiles towards the homogeneous solutions $`\rho _\mathrm{b}^0`$ and $`\rho _{\mathrm{ub}}^0`$, see Appendix B. The corresponding decay length $`\xi `$ satisfies the cubic equation

$$v_\mathrm{b}(12\rho _\mathrm{b}^0)\xi ^3+(D_\mathrm{b}+N_{\mathrm{ch}}D_{\mathrm{ub}}g)\xi ^2+N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{v_\mathrm{b}(12\rho _\mathrm{b}^0)}{A}\xi N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{D_\mathrm{b}}{A}=0$$
(17)
which is solved numerically.

For the ASEP in one dimension, one has $`N_{\mathrm{ch}}=0`$ and (17) reduces to $`v_\mathrm{b}(12\rho _\mathrm{b}^0)\xi ^3+D_\mathrm{b}\xi ^2=\xi ^2(\xi \xi _0)=0`$ with

$$\xi _0D_\mathrm{b}/v_\mathrm{b}(12\rho _\mathrm{b}^0).$$
(18)
Thus, in this limit, mean field theory leads to the correlation length $`\xi =\xi _0`$. For $`N_{\mathrm{ch}}>0`$, we choose the unique solution of (17) which approaches (18) as $`N_{\mathrm{ch}}`$ vanishes. This solution behaves as $`\xi \xi _0(1+gN_{\mathrm{ch}}D_{\mathrm{ub}}/D_\mathrm{b})`$ for small $`N_{\mathrm{ch}}`$.

In addition to the mean field calculation, we again used Monte Carlo simulations in order to determine the density profiles as shown in Figs. 5 and 6. As predicted by the mean field calculation, the constant bulk densities for the bound and unbound states are approached exponentially in the low and high density phases. The corresponding decay length $`\xi `$ is found to be the same for the bound and the unbound density and to diverge as one approaches the maximal current phase.

Some simulation results for the decay length $`\xi `$ are displayed in Fig. 7. In this case, the radius $`R`$ of the tube was varied while the boundary densities were kept fixed. The latter densities were chosen to be $`\rho _{\mathrm{b},\mathrm{in}}=0.38`$ and $`\rho _{\mathrm{b},\mathrm{ex}}=0.6`$ which lies within the low density phase but is close to the phase transition line which separates the low from the high density phase. One surprising feature of the Monte Carlo data for the decay length $`\xi `$ is that they exhibit a pronounced minimum as a function of tube radius $`R`$.

For comparison, we also display in Fig. 7 the $`\xi `$–values as obtained from several mean field approximations corresponding to the dashed, dotted and solid lines. The dashed line is obtained from the solution of equation (17) which we derived from the continuous mean field approximation for the two–state model. We also determined this quantity using a lattice version of this approximation (dotted line) and a more elaborate mean field approximation (solid line) in which we solved the diffusion equation in the cylindrical compartment and matched this solution to the directed transport along the filament.

Inspection of Fig. 7 shows that all three mean field approximations are quite consistent with each other and lead to $`\xi R^2`$ as follows from (17). Such an increase of $`\xi `$ for large $`R`$ is in fair agreement with the Monte Carlo data displayed in Fig. 7. However, in contrast to the Monte Carlo simulations, all three mean field approximations give a monotonic increase of $`\xi `$ with increasing $`R`$.

The largest discrepancy between the mean field results and the Monte Carlo data is found in the limit of small $`R`$ for which one recovers the 1–dimensional ASEP. In this limit, the Monte Carlo data should be quite reliable as one concludes from the value obtained for $`R=0`$ which is in very good agreement with the exact solution for the 1–dimensional ASEP as given by

$$\xi =\left|\frac{1}{\xi _{\mathrm{in}}}+\frac{1}{\xi _{\mathrm{ex}}}\right|^1\mathrm{with}\xi _k\frac{1}{\mathrm{ln}[4\rho _k(1\rho _k)]}$$
(19)
for $`k=`$ in, ex. Thus, we conclude that the decay length $`\xi `$ does indeed exhibit a nonmonotonic dependence on the tube radius $`R`$ and that this behavior is not reproduced by the mean field approximation.

Finally, we note that the different behavior of the decay length $`\xi `$ for large and for small $`R`$ is correlated with a qualitatively different behavior of the corresponding density profiles as observed in the Monte Carlo simulations. As an example, let us consider the density profiles within the low density phase as in Fig. 5. In this case, the bound and unbound densities exhibit plateau regions which are determined by their values $`\rho _{\mathrm{b},\mathrm{in}}`$ and $`\rho _{\mathrm{ub},\mathrm{in}}`$ at the left boundary. As one gets closer to the right boundary where the motors can leave the tube, the densities start to deviate from these constant values, and these deviations grow exponentially as $`\mathrm{exp}(x/\xi )`$. For small $`R`$, the corresponding profiles are convex upwards for all values of $`x`$. For large $`R`$, on the other hand, the profile exhibits an inflection point close to the right boundary. This inflection point moves towards the interior of the tube as $`R`$ is further increased.

#### 2 Maximal current phase

In the maximal current phase, one has $`\rho _\mathrm{b}^0=1/2`$, and one has to consider terms up to second order in the density deviations $`\eta _\mathrm{b}`$ and $`\eta _{\mathrm{ub}}`$, see Appendix B. One then finds that $`\eta _\mathrm{b}`$ satisfies the nonlinear differential equation

$$(D_\mathrm{b}+gN_{\mathrm{ch}}D_{\mathrm{ub}})\frac{}{x}\eta _\mathrm{b}v_\mathrm{b}\eta _\mathrm{b}^2$$
(20)
with

$$g\frac{\stackrel{~}{ϵ}(1\rho _{\mathrm{ub}}^0)+\stackrel{~}{\pi }_{\mathrm{ad}}\rho _{\mathrm{ub}}^0}{\stackrel{~}{\pi }_{\mathrm{ad}}(1\rho _\mathrm{b}^0)+\stackrel{~}{ϵ}\rho _\mathrm{b}^0}.$$
(21)
This differential equation can be solved by separation of variables. The solution behaves as

$$\eta _\mathrm{b}\frac{D_\mathrm{b}+gN_{\mathrm{ch}}D_{\mathrm{ub}}}{v_\mathrm{b}x}\mathrm{for}\mathrm{large}x$$
(22)
i.e., the deviation of the bound density from its asymptotic bulk value $`\rho _\mathrm{b}^0=1/2`$ decays as $`1/x`$ within the mean field approximation.

For $`N_{\mathrm{ch}}=0`$, this becomes identical with the mean field solution for the ASEP in one dimension as discussed in . In this latter case, an exact solution is available , which shows that the density profile decays as $`1/x^{1/2}`$, i.e., with a different exponent. This is related to the fact that density fluctuations spread superdiffusively in the 1–dimensional ASEP . The dispersion $`(\mathrm{\Delta }x)^2`$ of an ensemble of particles in such a system behaves as $`(\mathrm{\Delta }x)^2t^{4/3}`$ for large times $`t`$. This superdiffusive spreading of density fluctuations in one dimension can be taken into account within a mean field approximation if one considers a scale-dependent diffusion coefficient as shown in Ref. . In the following, we will extend this approach to the tube geometry considered here.

Thus, we now replace $`D_\mathrm{b}`$ in the mean field equation (20) by a scale–dependent diffusion coefficient $`D_\mathrm{b}(x)`$ and consider the modified mean field equation

$$[D_\mathrm{b}(x)+gN_{\mathrm{ch}}D_{\mathrm{ub}}]\frac{}{x}\eta _\mathrm{b}v_\mathrm{b}\eta _\mathrm{b}^2.$$
(23)
A convenient choice for $`D_\mathrm{b}(x)`$ which embodies the correct superdiffusive behavior for the ASEP in one dimension is

$$D_\mathrm{b}(x)D_{sc}\left(1+\sqrt{\frac{xx_0}{x_{sc}}}\right).$$
(24)
The left boundary is located at $`x=x_0`$, and $`D_{sc}`$ and $`x_{sc}`$ represent two scale parameters. With this choice, the modified mean field equation (23) can again be solved by separation of variables. As a result, we obtain

$$\eta _\mathrm{b}(x)=a\left[\frac{a}{\eta _\mathrm{b}(x_0)}+\sqrt{y}b\mathrm{ln}(1+\sqrt{y}/b)\right]^1$$
(25)
where we have introduced the abbreviations

$$y(xx_0)/x_{sc},$$
(26)
$$aD_{sc}/2v_bx_{sc},$$
(27)
and

$$b1+gN_{\mathrm{ch}}D_{\mathrm{ub}}/D_{sc}.$$
(28)
The ’initial’ value $`\eta _\mathrm{b}(x_0)=\rho _{\mathrm{b},\mathrm{in}}\frac{1}{2}`$ denotes the density deviation at the left boundary.

Far from the left boundary, i.e., for large values of $`yx`$, the expression (25) leads to the asymptotic behavior

$$\eta _\mathrm{b}(x)\frac{a}{\sqrt{y}}\left[1+b\frac{\mathrm{ln}(\sqrt{y})}{\sqrt{y}}\right].$$
(29)
Thus, the deviation of the bound density from its asymptotic value now decays as $`1/\sqrt{y}1/\sqrt{x}`$ as for the ASEP in one dimension.

In addition, we also obtain a correction term in (29) which depends on $`b=1+gN_{\mathrm{ch}}D_{\mathrm{ub}}/D_{sc}`$. For large tube radius $`R`$, one has $`N_{\mathrm{ch}}R^2`$ and $`bN_{\mathrm{ch}}R^2`$. Therefore, the correction term becomes large for large $`R`$. We can now define a crossover length $`y=y_{}`$ at which the correction term has the same size as the leading term. This leads to the implicit equation

$$\sqrt{y_{}}/\mathrm{ln}(\sqrt{y_{}})=b$$
(30)
and, thus, to

$$y_{}[b\mathrm{ln}(b)]^2[R^2\mathrm{ln}(R^2)]^2\mathrm{for}\mathrm{large}bR^2.$$
(31)
Close to the left boundary, i.e., for small values of $`yxx_0`$, the expression (25) for the deviation of the bound density from its asymptotic value $`\rho _\mathrm{b}^0=1/2`$ leads to

$$\eta _\mathrm{b}(x)\eta _\mathrm{b}(x_0)[1\frac{1}{2ab}\eta _\mathrm{b}(x_0)y].$$
(32)
We may now define a second crossover length or extrapolation length $`y=y_{}`$ at which the two terms in (32) cancel. This extrapolation length is given by

$$y_{}=2ab/\eta _\mathrm{b}(x_0)=2ab/(\rho _{\mathrm{b},\mathrm{in}}\frac{1}{2}).$$
(33)
For large tube radius $`R`$, the latter length scale grows as $`y_{}bR^2`$. In general, the extrapolation length can have both signs but, for the maximal current phase, one always has $`\rho _{\mathrm{b},\mathrm{in}}>1/2`$ and, thus, $`y_{}>0`$.

There is also an intermediate range of $`y`$–values defined by $`y_{}yb^2y_{}`$. For these $`y`$–values, the deviation $`\eta _\mathrm{b}`$ of the bound density from its asymptotic value as given by (25) simplifies and becomes

$$\eta _\mathrm{b}\frac{2ab}{y}=\frac{D_{sc}+gN_{\mathrm{ch}}D_{\mathrm{ub}}}{v_\mathrm{b}x_{sc}}\frac{1}{y}.$$
(34)
Thus, for these intermediate $`y`$–values the density deviation decays again as $`1/y`$.

In summary, the theory described here indicates (i) that the ’initial’ value at the left boundary is felt up to an extrapolation length $`y_{}bR^2`$, (ii) that the true asymptotic behavior of the density deviation is obtained for $`y>y_{}R^4`$ as follows from (31), and (iii) that the density deviation decays as $`1/y`$ on intermediate length scales with $`y_{}yy_{}`$.

These conclusions agree with the results of Monte Carlo simulations. In these simulations, the tube length $`L`$ is necessarily finite. This implies that the profiles observed in the simulations may not reach the asymptotic behavior present in a tube of infinite length. Since both the crossover length $`y_{}`$ and the extrapolation length $`y_{}`$ increase quickly with increasing tube radius $`R`$, we expect to find the true asymptotic behavior only for sufficiently small values of $`R`$. This expectation is confirmed by the simulation results. For sufficiently small $`R`$, the density deviation $`\eta _\mathrm{b}`$ is found to decay as $`1/x^{1/2}`$ for large $`x`$, as shown in Fig. 9 for $`R=3`$. For sufficiently large values of $`R`$, on the other hand, the observed profiles decay as $`1/x`$, see the Monte Carlo data in Fig. 9 for $`R=5`$ and $`R=7`$. In the latter cases, the true asymptotic behavior is not accessible and is cut off by the finite value of $`L`$.

In addition, a more detailed analysis of the Monte Carlo results for small values of $`R`$ but large values of $`L`$ explicitly shows that the decay of the bound density deviation behaves as $`1/\sqrt{x}`$ for large $`x`$ but decays faster than $`1/\sqrt{x}`$ for smaller values of $`x`$. This crossover behavior is shown in Fig. 10 where the simulation data are compared with the thin dashed line corresponding to the decay law $`1/\sqrt{x}`$. Furthermore, the data can be well fitted with a density profile as given by (25) if one makes an appropriate choice for the scale parameters $`D_{sc}`$ and $`x_{sc}`$.

The Monte Carlo data shown in Fig. 10 correspond to a tube of radius $`R=3`$ and length $`L=6000`$ with boundary densities $`\rho _{\mathrm{b},\mathrm{in}}=0.8`$ and $`\rho _{\mathrm{b},\mathrm{ex}}=0.5`$. In this case, a least–squares fit of the data for $`x>20`$ leads to the parameter values $`D_{sc}0.81`$ and $`x_{sc}8.57`$. The corresponding fitting curve corresponds to the thick dashed line in Fig. 10. The fit becomes less reliable close to the left boundary, where the assumption that $`\eta _\mathrm{b}`$ is small is no longer fulfilled. For fixed boundary densities, both fitting parameters $`D_{sc}`$ and $`x_{sc}`$ are found to depend on $`R`$ and to increase with increasing $`R`$. In addition, for $`R4`$, it becomes rather difficult to determine these parameters since one would have to simulate rather long systems with $`L6000`$ in order to observe the true asymptotic behavior.

For the ASEP in one dimension, it has been argued that the power–law decay of the density profile as given by the asymptotic form $`c/\sqrt{x/\mathrm{}}`$ for large $`x`$ is characterized by the universal amplitude $`c=1/2\sqrt{\pi }0.282`$.  Inspection of the relations (26), (27), and (29) shows that, in the present situation, $`c=D_{sc}/2v_b\sqrt{x_{sc}}`$ where the parameters $`D_{sc}`$ and $`x_{sc}`$ depend on the tube radius $`R`$. This implies that the amplitude $`c`$ will also depend on $`R`$. Using the numerically determined values for $`D_{sc}`$ and $`x_{sc}`$ and the bound state velocity $`v_b0.99`$, we obtain the estimates $`c0.27`$ for $`R=0`$, which corresponds to the 1–dimensional ASEP and should be compared with the exact value $`c=1/2\sqrt{\pi }0.282`$, and $`c0.14`$ for $`R=3`$. Thus, for small values of $`R`$, the amplitude $`c`$ is found to decrease with increasing $`R`$.

## V Diffusive injection and extraction of motors

Finally, we consider the second type of open boundary conditions corresponding to case (C) in Fig. 2. The length of the tube is again denoted by $`L`$. This tube is now longer than the filament which has length $`L_\mathrm{F}<L`$. The left end of the tube is located at $`x=0`$ as before but the left end of the filament is at $`x=\mathrm{\Delta }L(LL_\mathrm{F})/2`$. Likewise, the right end of the filament is at $`x=\mathrm{\Delta }L+L_\mathrm{F}`$ whereas the right end of the tube is at $`x=2\mathrm{\Delta }L+L_\mathrm{F}=L`$. Thus, a motor particle which enters the tube on the left must diffuse over a distance $`\mathrm{\Delta }L`$ before it can come into contact with the filament, and a motor particle which leaves the filament at its right end must also diffuse over a distance $`\mathrm{\Delta }L`$ before it can leave the tube.

At the left and right end of the tube, we now prescribe constant boundary densities as given by

$$\rho (x=0,y,z)=\rho _{\mathrm{ub},\mathrm{in}}\mathrm{and}\rho (x=L+1,y,z)=\rho _{\mathrm{ub},\mathrm{ex}}$$
(35)
for all values of $`y`$ and $`z`$ with $`y^2+z^2R^2`$.

As before, the jump probability $`\beta `$ to make backward steps on the filament is taken to be zero. Likewise, the resting probability $`\gamma `$ to make no step at all on the filament is also zero with a possible exception at the ’last’ filament site with $`(x,y,z)=(\mathrm{\Delta }L+L_\mathrm{F},0,0)`$. Indeed, in order to define the system in a unique way, we still have to specify the probability to make a forward step at this ’last’ filament site. Two possible choices appear rather natural: (i) Active unbinding from the ’last’ filament site, i.e., the motor particle attempts to step forward with probability $`\alpha `$ and makes a step if the adjacent nonfilament site is unoccupied. In this case, the forward step at the last filament site is governed by the same probabilities as all other forward steps along the filament; and (ii) Thermal unbinding in which the motor particle unbinds with probability $`ϵ/6`$ both in the forward direction and in the four orthogonal directions. In the latter case, one has to choose the resting probability $`\gamma `$ to be nonzero and to be given by $`\gamma =15ϵ/6`$.

The choice (i) is suggested by the results of recent experiments on microtubules and kinesin motors  which indicate that these motors unbind quickly at the filament ends. In the following subsections, we will first consider this choice (i) corresponding to active unbinding from the last filament site. In the last subsection V D, we will show that the choice (ii) leads to rather similar behavior.

### A Diffusive bottlenecks

In order to understand the behavior found for boundary condition (C), it is instructive to partition the tube into three compartments which are defined as follows: (i) A left compartment with $`1x<\mathrm{\Delta }L`$ where transport is purely diffusive; (ii) A middle compartment with $`\mathrm{\Delta }Lx\mathrm{\Delta }L+L_\mathrm{F}`$ where all directed (or active) transport occurs; and (iii) A right compartment with $`\mathrm{\Delta }L+L_\mathrm{F}<xL`$ where the transport is again purely diffusive.

For a stationary state, the total current through the tube must be constant and, thus, must be the same in all three compartments. The current through the middle compartment is given by the bound current $`j_\mathrm{b}=v_\mathrm{b}\rho _\mathrm{b}^0(1\rho _\mathrm{b}^0)`$. Thus, the diffusive currents $`J_{\mathrm{dif},\mathrm{L}}`$ and $`J_{\mathrm{dif},\mathrm{R}}`$ in the left and right tube segment must be equal and must satisfy the simple relation

$$J_{\mathrm{dif},\mathrm{L}}=J_{\mathrm{dif},\mathrm{R}}=v_\mathrm{b}\rho _\mathrm{b}^0(1\rho _\mathrm{b}^0).$$
(36)
The relation as given by (36) is easily checked in the simulations since the density profile is found to be approximately linear in the left and in the right compartments provided $`\mathrm{\Delta }L`$ is sufficiently large. For such a linear density profile in the right compartment, the diffusive current $`J_{\mathrm{dif},\mathrm{R}}`$ can be estimated as

$$J_{\mathrm{dif},\mathrm{R}}(1+N_{\mathrm{ch}})D_{\mathrm{ub}}\frac{|\rho _{\mathrm{ub}}(x=\mathrm{\Delta }L+L_\mathrm{F})\rho _{\mathrm{ub},\mathrm{ex}}|}{\mathrm{\Delta }L}.$$
(37)
Since the maximal density difference is one (in the units used here), the maximal diffusive current behaves as

$$\mathrm{max}(J_{\mathrm{dif},\mathrm{R}})(1+N_{\mathrm{ch}})D_{\mathrm{ub}}\frac{1}{\mathrm{\Delta }L}R^2D_{\mathrm{ub}}\frac{1}{\mathrm{\Delta }L}.$$
(38)
This shows that the maximal diffusive current depends on the tube radius $`R`$, the lateral size $`\mathrm{\Delta }L`$ of the diffusive segments, and the diffusion coefficient $`D_{\mathrm{ub}}`$ of the unbound motors. Since these parameters can be chosen independently from the bound motor velocity $`v_\mathrm{b}`$, the diffusive current $`J_{\mathrm{dif},\mathrm{L}}`$ can be made smaller than the maximal bound current $`v_\mathrm{b}/4`$ on the filament. In the latter case, the diffusive compartments act as diffusive bottlenecks and the maximal current phase characterized by the current $`v_\mathrm{b}/4`$ is expected to be absent from the phase diagram. This expectation is indeed confirmed by the simulations as discussed next.

### B Phase diagram without maximal current phase

One geometry, for which no maximal current phase has been observed in the simulations, is provided by a tube with radius $`R=5`$, length $`L=600`$, and filament length $`L_\mathrm{F}=590`$. The corresponding phase diagram as determined by the Monte Carlo simulations is shown in Fig. 11. The largest part of the phase diagram is covered by the high density phase; in addition, a low density phase is found for small values of the boundary densities $`\rho _{\mathrm{ub},\mathrm{in}}`$ and $`\rho _{\mathrm{ub},\mathrm{ex}}`$.

The transition line displayed in Fig. 11 has been determined from the functional dependence of the bound density in the bulk, $`\rho _\mathrm{b}^0`$, on the left boundary density $`\rho _{\mathrm{ub},\mathrm{in}}`$ as shown in Fig. 12. Inspection of this latter figure shows that the bound bulk density $`\rho _\mathrm{b}^0`$ jumps at a certain value of the left boundary density $`\rho _{\mathrm{ub},\mathrm{in}}`$. Since $`L`$ is finite, this jump occurs over a small but finite interval of $`\rho _{\mathrm{ub},\mathrm{in}}`$. Thus, the jump can be characterized by two $`\rho _{\mathrm{ub},\mathrm{in}}`$–values, say $`\rho _{\mathrm{ub},\mathrm{in}}^<`$ and $`\rho _{\mathrm{ub},\mathrm{in}}^>`$, which represent the left and the right ’corner’ of the numerically determined jump. Both $`\rho _{\mathrm{ub},\mathrm{in}}`$–values have been included in the phase diagram of Fig. 11.

For $`\rho _{\mathrm{ub},\mathrm{in}}^<<\rho _{\mathrm{ub},\mathrm{in}}<\rho _{\mathrm{ub},\mathrm{in}}^>`$, the simulations do not reach a stationary state within two days of computation. Simulations also become very slow in the high density phase especially when the overall motor concentration gets so large that the bound density is close to one and the unbound density is no longer small compared to the bound density.

### C Presence of the maximal current phase

In order to estimate the set of parameters, for which the phase diagram exhibits a maximal current phase, we return to the estimate (37) for the diffusive current $`J_{\mathrm{dif},\mathrm{R}}`$ and make the simplifying assumption that the bound and unbound densities in the middle compartment of the tube are essentially constant. Thus, we replace the true unbound density $`\rho _{\mathrm{ub}}(x=\mathrm{\Delta }L+L_\mathrm{F})`$ at the right end of the filament by its bulk value $`\rho _{\mathrm{ub}}^0`$.

The bulk density $`\rho _{\mathrm{ub}}^0`$ of the unbound motors is related to the bulk density $`\rho _\mathrm{b}^0`$ of the bound motors via the radial equilibrium relation (3). For the maximal current phase, one has $`\rho _\mathrm{b}^0=1/2`$ and (3) leads to $`\rho _{\mathrm{ub}}^0=\frac{ϵ/\pi _{\mathrm{ad}}}{1+ϵ/\pi _{\mathrm{ad}}}ϵ/\pi _{\mathrm{ad}}`$ for small $`ϵ`$. In this way, we arrive at the estimate

$$J_{\mathrm{dif},\mathrm{R}}(1+N_{\mathrm{ch}})D_{\mathrm{ub}}\frac{ϵ}{\pi _{\mathrm{ad}}\mathrm{\Delta }L}\mathrm{for}\rho _{\mathrm{ub},\mathrm{ex}}=0.$$
(39)
The maximal current phase should be present in the phase diagram if the diffusive current $`J_{\mathrm{dif},\mathrm{R}}`$ exceeds the maximal current $`v_\mathrm{b}/4`$ on the filament. It then follows from (39) that the maximal current phase should be present for $`\rho _{\mathrm{ub},\mathrm{ex}}=0`$ provided the tube radius $`R`$ satisfies

$$R^2>R_{}^2\frac{v_\mathrm{b}}{4\pi }\frac{\pi _{\mathrm{ad}}\mathrm{\Delta }L}{D_{\mathrm{ub}}ϵ}$$
(40)
where $`1+N_{\mathrm{ch}}\pi R^2`$ has been used. Using the same line of reasoning, a second, less restrictive condition can be obtained from an estimate for the diffusive current within the left compartment.

The threshold value $`R_{}`$ for the tube radius as given by (40) has been confirmed by Monte Carlo simulations for the jump probabilities $`\alpha =12ϵ/3`$ with $`ϵ=1/100`$, the sticking probability $`\pi _{\mathrm{ad}}=1`$, the compartment size $`\mathrm{\Delta }L=5`$, and the boundary densities $`\rho _{\mathrm{ub},\mathrm{in}}=0.2`$ and $`\rho _{\mathrm{ub},\mathrm{ex}}=0`$. The jump probabilities imply the bound motor velocity $`v_\mathrm{b}=12/300`$; the diffusion coefficient $`D_{\mathrm{ub}}`$ of the unbound motors has the value $`D_{\mathrm{ub}}=1/6`$ (since the resting probability $`\gamma =0`$ as mentioned above). When these parameter values are inserted into (40), one obtains the estimate $`R_{}15.4`$.

The corresponding Monte Carlo data are displayed in Fig. 13. Inspection of this figure shows that the current does indeed attain its maximal value $`v_\mathrm{b}/4`$ for $`RR_{}`$ with $`R_{}16`$. The same Fig. 13 also shows the transition from the low density to the high density phase which occurs for a tube radius $`R_{}`$ which satisfies $`4<R_{}<5`$.

A complete phase diagram for a tube with radius $`R=17`$ is shown in Fig. 14. Again, most of the phase diagram is covered by the high density phase (HD), while a low density phase (LD) is found only for very small values of $`\rho _{\mathrm{ub},\mathrm{in}}`$. The maximal current phase (MC) is present now but only for very small values of $`\rho _{\mathrm{ub},\mathrm{ex}}`$.

It is interesting to note that similar effects also occur in the purely one-dimensional system if one considers a driven system which is bounded by two segments which exhibit only diffusive transport. In this 1–dimensional case, quantitative predictions can be made using a mean field approximation as will be discussed elsewhere .

### D Active versus thermal unbinding from the ’last’ filament site

The Monte Carlo data displayed so far have been obtained for active unbinding of a motor particle which is bound to the ’last’ filament site. As mentioned above, another possibility is that the motor particle gets stuck at the ’last’ filament site and unbinds only by thermal excitations, i.e., with unbinding probability $`ϵ/6`$. For these two different unbinding mechanisms, one will, in general, obtain different density profiles. However, this difference is not dramatic as one can see from Fig. 15 which exhibits density profiles for both cases. Although the probability for a forward step at the ’last’ filament site differs by two orders of magnitude for the two cases, the bulk density exhibits a relatively small difference. The bound density, on the other hand increases and decreases close to the right end of the filament for thermal and active unbinding from the ’last’ site, respectively.

## VI Summary and conclusions

Let us summarize the main results and add a few comments. We have studied a lattice model for the motion of many molecular motors in an open tube which contains a single filament. When bound to the filament, the motor particles undergo an asymmetric simple exclusion process (ASEP). In addition, motors can unbind from the filament and then diffuse freely in the tube. As for the ASEP in one dimension, the motor traffic in open tubes can exhibit three different phases: high density and low density phases which are characterized by an exponential decay of the density deviations from their bulk values and maximal current phases characterized by an algebraic decay. Therefore, the molecular motor traffic in open tubes are promising candidates for the experimental observation of boundary–induced non–equilibrium phase transitions.

In general, the location of the transition lines is found to depend on the precise choice of the boundary conditions. Apart from periodic boundary conditions, case (A), we studied two different boundary conditions (B) and (C) for an open tube. In case (B), the bound and unbound densities are kept fixed at the boundaries and satisfy radial equilibrium. For this case, the location of the transition lines is independent of the model parameters, and the phase diagram of the ASEP in one dimension is recovered. In case (C), the active compartment of the tube is bounded by two compartments where the transport is purely diffusive. In this latter case, the phase diagram depends on the geometry of the tube and on the transport properties in the bound and unbound motor states. In many cases, the maximal current phase is completely suppressed by the coupling to the diffusive compartments which act as bottle necks for the transport.

The theoretical results described here should be accessible to experiments on cytoskeletal filaments and motors. In particular, the motor traffic through open tubes as discussed here provides new opportunities to study the transport properties of ASEPs by systematic experiments.

## Acknowledgments

We thank Theo M. Nieuwenhuizen for comments and discussions during the initial stage of this work.

## A Radial equilibrium for periodic boundary conditions

In the following, we show that Eq. (3), the condition for radial equilibrium, holds exactly in the case of periodic boundary conditions using the quantum Hamiltonian formalism . The exact stationary master equation can be written in the form

$$H|\rho =0,$$
(A1)
where $`H`$ is the ’quantum Hamiltonian’ of the stochastic process and $`|\rho `$ is a vector in a product Hilbert space; each lattice site is represented by a two-state system with the orthogonal vectors $`|1`$ for an occupied lattice site (”spin down”) and $`|0`$ for a vacancy (”spin up”). Because of translational invariance in the direction parallel to the filament, there cannot be any radial currents in our system, and the unbound density is independent of the radial coordinate; thus we can restrict the analysis to the case of one unbound channel. We denote by $`|\rho _{k,\mathrm{b}}`$ and $`|\rho _{k,\mathrm{ub}}`$ the state of site $`k`$ of the bound and unbound channel, respectively. Using the general recipe given in chapter 2 of Ref.  we construct the ’quantum Hamiltonian’ $`H`$ of our system:

$$H=H_1+H_2+H_3,$$
(A2)
where $`H_1`$ represents the dynamics of the asymmetric exclusion process in the bound channel, $`H_2`$ the symmetric exclusion process in the unbound channel and $`H_3`$ the coupling of the two channels. Each term can be written as a sum $`H_i=_kh_k^{(i)}`$ with

$`h_k^{(1)}`$ $`=`$ $`v_\mathrm{b}\left(n_{k,\mathrm{b}}v_{k+1,\mathrm{b}}s_{k,\mathrm{b}}^+s_{k+1,\mathrm{b}}^{}\right)`$ (A3)
$`h_k^{(2)}`$ $`=`$ $`D_{\mathrm{ub}}\left(n_{k,\mathrm{ub}}v_{k+1,\mathrm{ub}}s_{k,\mathrm{ub}}^+s_{k+1,\mathrm{ub}}^{}+v_{k,\mathrm{ub}}n_{k+1,\mathrm{ub}}s_{k,\mathrm{ub}}^{}s_{k+1,\mathrm{ub}}^+\right)`$ (A4)
$`h_k^{(3)}`$ $`=`$ $`{\displaystyle \frac{ϵ}{6}}\left(n_{k,\mathrm{b}}v_{k,\mathrm{ub}}s_{k,\mathrm{b}}^+s_{k,\mathrm{ub}}^{}\right)+{\displaystyle \frac{\pi _{\mathrm{ad}}}{6}}\left(v_{k,\mathrm{b}}n_{k,\mathrm{ub}}s_{k,\mathrm{b}}^{}s_{k,\mathrm{ub}}^+\right).`$ (A5)
$`s_{k,\mathrm{b}}^{}`$ is a creation operator for a particle at site $`k`$ of the bound channel, and $`s_{k,\mathrm{b}}^+`$ is the corresponding annihilation operator. $`n_{k,\mathrm{b}}`$ is the particle number operator at site $`k`$ of the bound channel and $`v_{k,\mathrm{b}}=1n_{k,\mathrm{b}}`$. Operators for the unbound channel are defined in an analogous manner. The off-diagonal parts of the operators $`h_k^{(i)}`$ represent hopping of particles, while the diagonal parts are determined by conservation of probability, see chapter 2 of Ref. .

We now show that the product measure

$$|\rho =\underset{k}{}\left(\left[(1\rho _\mathrm{b})|0_{k,\mathrm{b}}+\rho _\mathrm{b}|1_{k,\mathrm{b}}\right]\left[(1\rho _{\mathrm{ub}})|0_{k,\mathrm{ub}}+\rho _{\mathrm{ub}}|1_{k,\mathrm{ub}}\right]\right),$$
(A6)
is a stationary state, if Eq. (3) holds, i.e. that the radial equilibrium condition implies $`H|\rho =0`$. The product measure $`|\rho `$ defines a state where the density is $`\rho _\mathrm{b}`$ at each lattice site of the bound channel and $`\rho _{\mathrm{ub}}`$ at each site of the unbound channel and spatial correlations vanish.

$`H_1`$ and $`H_2`$ do not couple the bound and unbound channel, therefore we can consider them separately and refer to the result, that for the symmetric as well as for the asymmetric exclusion process, the product measure $`|\rho `$ is stationary in the case of periodic boundary conditions . For a proof using the quantum Hamiltonian formalism, see chapter 7.1.2 of Ref. : In the case of the ASEP in one dimension, for example, one can easily check, that $`h_k^{(1)}|\rho =v_\mathrm{b}(n_{k,\mathrm{b}}n_{k+1,\mathrm{b}})|\rho `$, therefore the summation gives zero for periodic boundary conditions. Hence $`H_1|\rho =0`$ and $`H_2|\rho =0`$.

Concerning $`H_3`$, it is sufficient to consider a single site $`k`$ in both channels. Doing some steps, we find

$$\left(n_{k,\mathrm{b}}v_{k,\mathrm{ub}}s_{k,\mathrm{b}}^+s_{k,\mathrm{ub}}^{}\right)|\rho _{k,\mathrm{b}}|\rho _{k,\mathrm{ub}}=\rho _\mathrm{b}(1\rho _{\mathrm{ub}})\left(|1_{k,\mathrm{b}}|0_{k,\mathrm{ub}}|0_{k,\mathrm{b}}|1_{k,\mathrm{ub}}\right)$$
(A7)
and

$$\left(v_{k,\mathrm{b}}n_{k,\mathrm{ub}}s_{k,\mathrm{b}}^{}s_{k,\mathrm{ub}}^+\right)|\rho _{k,\mathrm{b}}|\rho _{k,\mathrm{ub}}=\rho _{\mathrm{ub}}(1\rho _\mathrm{b})\left(|0_{k,\mathrm{b}}|1_{k,\mathrm{ub}}|1_{k,\mathrm{b}}|0_{k,\mathrm{ub}}\right)$$
(A8)
for the product measure (A6). Hence from

$$h_k^{(3)}|\rho _{k,\mathrm{b}}|\rho _{k,\mathrm{ub}}=\left(\frac{ϵ}{6}\rho _\mathrm{b}(1\rho _{\mathrm{ub}})\frac{\pi _{\mathrm{ad}}}{6}\rho _{\mathrm{ub}}(1\rho _\mathrm{b})\right)\left(|1_{k,\mathrm{b}}|0_{k,\mathrm{ub}}|0_{k,\mathrm{b}}|1_{k,\mathrm{ub}}\right)=0,$$
(A9)
we obtain Eq. (3), which is the condition for radial equilibrium.

## B Continuum two-state mean field equations for open boundaries

In this appendix, we give the details of the continuum mean field approximation for the two–state model introduced in section IV B. We insert the decomposition $`\rho _\mathrm{b}=\rho _\mathrm{b}^0+\eta _\mathrm{b}`$ and $`\rho _{\mathrm{ub}}=\rho _{\mathrm{ub}}^0+\eta _{\mathrm{ub}}`$ as introduced in (16) into the mean field equations (12) and (13) and expand these equations up to second order in the density deviations $`\eta _\mathrm{b}`$ and $`\eta _{\mathrm{ub}}`$. As a result, we obtain

$$D_\mathrm{b}\frac{}{x}\eta _\mathrm{b}+N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{}{x}\eta _{\mathrm{ub}}=v_\mathrm{b}(12\rho _\mathrm{b}^0)\eta _\mathrm{b}v_\mathrm{b}\eta _\mathrm{b}^2$$
(B1)
and

$$v_\mathrm{b}(12\rho _\mathrm{b}^0)\frac{}{x}\eta _\mathrm{b}2v_\mathrm{b}\eta _\mathrm{b}\frac{}{x}\eta _\mathrm{b}D_\mathrm{b}\frac{^2}{x^2}\eta _\mathrm{b}=A\eta _{\mathrm{ub}}B\eta _\mathrm{b}+(\stackrel{~}{ϵ}\stackrel{~}{\pi }_{\mathrm{ad}})\eta _{\mathrm{ub}}\eta _\mathrm{b},$$
(B2)
with

$$A\stackrel{~}{\pi }_{\mathrm{ad}}(1\rho _\mathrm{b}^0)+\stackrel{~}{ϵ}\rho _\mathrm{b}^0\mathrm{and}B\stackrel{~}{ϵ}(1\rho _{\mathrm{ub}}^0)+\stackrel{~}{\pi }_{\mathrm{ad}}\rho _{\mathrm{ub}}^0.$$
(B3)
Note that radial equilibrium for the densities $`\rho _\mathrm{b}`$ and $`\rho _{\mathrm{ub}}`$ would imply that the right hand side of (B2) vanishes.

Let us first consider the case $`\rho _\mathrm{b}^01/2`$ which applies to the high and low density phases. In this case, we can neglect the second order terms and obtain two equations which are linear in the density deviations. The second equation (B2) can then be solved for $`\eta _{\mathrm{ub}}`$ which leads to

$$\eta _{\mathrm{ub}}=g\eta _\mathrm{b}+\frac{v_\mathrm{b}(12\rho _\mathrm{b}^0)}{A}\frac{}{x}\eta _\mathrm{b}\frac{D_\mathrm{b}}{A}\frac{^2}{x^2}\eta _\mathrm{b}$$
(B4)
with

$$g\frac{B}{A}=\frac{\stackrel{~}{ϵ}(1\rho _{\mathrm{ub}}^0)+\stackrel{~}{\pi }_{\mathrm{ad}}\rho _{\mathrm{ub}}^0}{\stackrel{~}{\pi }_{\mathrm{ad}}(1\rho _\mathrm{b}^0)+\stackrel{~}{ϵ}\rho _\mathrm{b}^0}.$$
(B5)
When this expression is inserted into (B1), we obtain the first order relation

$$(D_\mathrm{b}+N_{\mathrm{ch}}D_{\mathrm{ub}}g)\frac{}{x}\eta _\mathrm{b}+N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{v_\mathrm{b}(12\rho _\mathrm{b}^0)}{A}\frac{^2}{x^2}\eta _\mathrm{b}N_{\mathrm{ch}}D_{\mathrm{ub}}\frac{D_\mathrm{b}}{A}\frac{^3}{x^3}\eta _\mathrm{b}=v_\mathrm{b}(12\rho _\mathrm{b}^0)\eta _\mathrm{b}.$$
(B6)
for $`\eta _\mathrm{b}`$. We now make the Ansatz $`\eta _\mathrm{b}\mathrm{exp}(x/\xi )`$ which leads to the cubic equation (17) for the decay length $`\xi `$.

Now we consider the maximal current phase, i.e. the case $`\rho _\mathrm{b}^0=1/2`$. In this case the linear terms are zero. Furthermore, we neglect terms of order $`\eta _\mathrm{b}\frac{\eta _\mathrm{b}}{x}`$ and $`\frac{^2\eta _\mathrm{b}}{x^2}`$ as can be justified a posteriori since $`\eta _\mathrm{b}`$ is found to decay as an inverse power of $`x`$. Thus, up to leading order, we can ignore the left hand side of (B2), and the right hand side of (B2) vanishes. This implies radial equilibrium for the asymptotic decay to the homogeneous solution. Up to this order, we find

$$\eta _{\mathrm{ub}}=\frac{B\eta _\mathrm{b}}{A+(\stackrel{~}{ϵ}\stackrel{~}{\pi }_{\mathrm{ad}})\eta _\mathrm{b}}g\eta _\mathrm{b}\frac{B(\stackrel{~}{ϵ}\stackrel{~}{\pi }_{\mathrm{ad}})}{A^2}\eta _\mathrm{b}^2$$
(B7)
for small $`\eta _\mathrm{b}`$ and therefore

$$\frac{}{x}\eta _{\mathrm{ub}}=g\frac{}{x}\eta _\mathrm{b}+O(\eta _\mathrm{b}\frac{\eta _\mathrm{b}}{x}).$$
(B8)
If this latter expression is inserted into (B1), we obtain

$$(D_\mathrm{b}+gN_{\mathrm{ch}}D_{\mathrm{ub}})\frac{}{x}\eta _\mathrm{b}v_\mathrm{b}\eta _\mathrm{b}^2$$
(B9)
which leads to

$$\eta _\mathrm{b}\frac{D_\mathrm{b}+gN_{\mathrm{ch}}D_{\mathrm{ub}}}{v_\mathrm{b}x}\mathrm{for}\mathrm{large}x.$$
(B10)
For $`N_{\mathrm{ch}}=0`$, this asymptotic behavior is identical to the one found from the mean field approximation for the ASEP in one dimension.

## C List of symbols

| $`a`$ | abbreviation for $`D_{sc}/2v_bx_{sc}`$ used in subsection IV B 2 |
| --- | --- |
| (A) | periodic boundary conditions, see Fig. 2 |
| $`A`$ | abbreviation used in Appendix B, defined in Eq. (B3) |
| $`\alpha `$ | probability for a forward step of the bound motor |
| $`b`$ | abbreviation for $`1+gN_{\mathrm{ch}}D_{\mathrm{ub}}/D_{sc}`$ used in subsection IV B 2 |
| (B) | open tube, radial equilibrium at the boundaries, see Fig. 2 |
| $`B`$ | abbreviation used in Appendix B, defined in Eq. (B3) |
| $`\beta `$ | probability for a backward step of the bound motor |
| (C) | open tube, diffusive injection and extraction of motors, see Fig. 2 |
| $`D`$ | diffusion coefficient for the 1–dimensional ASEP |
| $`D_\mathrm{b}`$ | diffusion coefficient of the bound motor particle |
| $`\mathrm{\Delta }L`$ | distance between filament end and tube end for boundary conditions (C) |
| $`D_{\mathrm{ub}}`$ | diffusion coefficient of the unbound motor particle |
| $`D_{sc}`$ | scale parameter for the scale-dependent diffusion coefficient |
| $`ϵ`$ | unbinding probability of a bound motor |
| $`\stackrel{~}{ϵ}`$ | rescaled unbinding probability $`2ϵ/3`$ |
| $`\eta _\mathrm{b}(x)`$ | deviation of the bound density from the constant bulk value |
| $`\eta _{\mathrm{ub}}(x)`$ | deviation of the unbound density from the constant bulk value |
| $`g`$ | parameter defined by Eq. (21) |
| $`\gamma `$ | resting probability of a bound motor |
| $`J`$ | (global) current through the tube |
| $`j_\mathrm{b}`$ | local bound current |
| $`j_{\mathrm{ub}}`$ | local unbound current |
| $`J_{\mathrm{dif},\mathrm{L}}`$ | diffusive current in the left diffusive compartment for boundary conditions (C) |
| $`J_{\mathrm{dif},\mathrm{R}}`$ | diffusive current in the right diffusive compartment for boundary conditions (C) |
| $`L`$ | length of the tube |
| $`\mathrm{}`$ | lattice constant on filament |
| $`L_F`$ | length of the filament |
| $`N_{\mathrm{ch}}`$ | number of unbound channels |
| $`N_{\mathrm{mo}}`$ | number of motors in the tube which is fixed for periodic boundary conditions (A) |
| $`\varphi `$ | cross section of the tube |
| $`\pi _{\mathrm{ad}}`$ | sticking probability for a motor hopping to the filament |
| $`\stackrel{~}{\pi }_{\mathrm{ad}}`$ | rescaled sticking probability $`2\pi _{\mathrm{ad}}/3`$ |
| $`r`$ | radial spatial coordinate |
| $`R`$ | radius of the tube |
| $`R_{}`$ | minimal tube radius for which a maximal current phase is found in case (C) |
| $`\rho (x)`$ | density in the case of the 1–dimensional ASEP |
| $`\rho _\mathrm{b}(x)`$ | density of motors bound to the filament |
| $`\rho _{\mathrm{b},\mathrm{ex}}`$ | right boundary density on the filament |
| $`\rho _{\mathrm{b},\mathrm{in}}`$ | left boundary density on the filament |
| $`\rho _\mathrm{b}^0`$ | constant bulk density on the filament |
| $`\rho _{\mathrm{ex}}`$ | right boundary density for the 1–dimensional ASEP |
| $`\rho _{\mathrm{in}}`$ | left boundary density for the 1–dimensional ASEP |
| $`\rho _{\mathrm{mo}}`$ | overall motor concentration |
| $`\rho ^0`$ | constant bulk density for the 1–dimensional ASEP |
| $`\rho _{\mathrm{ub}}(x,r)`$ | density of unbound motors |
| $`\rho _{\mathrm{ub},\mathrm{ex}}`$ | density of unbound motors at the right boundary |
| $`\rho _{\mathrm{ub},\mathrm{in}}`$ | density of unbound motors at the left boundary |
| $`\rho _{\mathrm{ub}}^0`$ | unbound constant bulk density |
| $`\tau `$ | basic time unit |
| $`v`$ | velocity for the 1–dimensional ASEP |
| $`v_\mathrm{b}`$ | velocity of bound motor |
| $`v_c`$ | velocity of density fluctuations |
| $`v_s`$ | domain wall velocity |
| $`x`$ | spatial coordinate parallel to the filament |
| $`x_{sc}`$ | scale parameter for the scale-dependent diffusion coefficient, see Eq. (24) |
| $`x_{},x_{}`$ | crossover lengths |
| $`\xi `$ | decay length for the density deviations |
| $`\xi _0`$ | localization length for the 1-dimensional ASEP |
| $`x_0`$ | spatial location of the left tube end |
| $`y`$ | rescaled and shifted spatial coordinate, $`y(xx_0)/x_{sc}`$ |
| $`y_{},y_{}`$ | rescaled crossover lengths |


format: text

text:
# 1 Introduction

## 1 Introduction

The wave function renormalization constant $`Z`$ satisfies the unitarity bound $`0Z1`$, because it is interpreted as the probability of finding a bare particle in the physical one-particle state . In perturbation theory, this bound implies the positivity of the coefficient $`c`$ in the one-loop contribution to the wave function renormalization constant,

$$Z=1cg^2\mathrm{log}\frac{\mathrm{\Lambda }^2}{\mu ^2}1,(c>0)$$
(1.1)
where $`g,\mathrm{\Lambda }`$ and $`\mu `$ denote the coupling constant, the ultraviolet cut-off and the renormalization mass scale, respectively. This inequality for $`Z`$ implies the positivity of the anomalous dimension in the one-loop approximation:

$$\gamma =\frac{1}{2}\mu \frac{}{\mu }\mathrm{log}Z=cg^20.$$
(1.2)
The basic assumption of the unitarity bound is the positivity of the Hilbert space. However, often this assumption does not hold in the covariant quantization of gauge theories, where the time component of the gauge field has a negative signature relative to the transverse components. The anomalous dimension of the electron in the one-loop approximation, for example,

$$\gamma _e=\frac{e^2\alpha }{(4\pi )^2},$$
(1.3)
may be negative, depending on the choice of the gauge parameter $`\alpha `$. Although the total Hilbert space has an indefinite metric in gauge theories, the physical Hilbert space has a positive semidefinite norm and has a consistent probability interpretation . The anomalous dimensions of the physical operators, of course, satisfy positivity and are gauge independent.

We sometimes encounter negative anomalous dimensions in nonlinear sigma models (NL$`\sigma `$Ms), which are renormalizable in two dimensions . For example, the $`O(N+1)`$ invariant NL$`\sigma `$M, defined on the sphere $`S^N`$, is described by the Lagrangian

$$=\frac{1}{2}\frac{\left(_\mu \stackrel{}{\phi }\right)^2}{\left(1+\frac{g^2}{4}\stackrel{}{\phi }^2\right)^2}\left(\stackrel{}{\phi }=(\phi ^1,\mathrm{},\phi ^N)\right),$$
(1.4)
and has the negative anomalous dimension

$$\gamma _\phi =\frac{N}{8\pi }g^2<0.$$
(1.5)
There are neither higher-order derivative terms, which often introduce a ghost with indefinite metric, nor gauge fields in this model. What is wrong in the proof of the unitarity bound of the wave function renormalization constant? This is the question we study in this report.

This paper is organized as follows. In §2, we briefly review the unitarity bound of the wave function renormalization constant. In §3, we discuss the wave function renormalization constant in NL$`\sigma `$M in two dimensions. In §4, we describe the coordinate dependence of the anomalous dimension. In §5, we summarize our results.

## 2 Unitarity bound of the wave function renormalization constant

Here we give a brief derivation of the unitarity bound of the wave function renormalization constant. In the Källén-Lehmann representation , the commutation relation of the Heisenberg operators $`\varphi (x)`$ and $`\varphi (y)`$ is expressed as a superposition of the Jordan-Pauli function $`i\mathrm{\Delta }(xy;m^2)`$ for a particle with mass $`m`$:

$`0\left|[\varphi (x),\varphi (y)]\right|0={\displaystyle _0^{\mathrm{}}}\rho (m^2)i\mathrm{\Delta }(xy;m^2)𝑑m^2.`$ (2.1)
Here, the spectral function

$$\rho (k^2)2\pi \underset{n}{}\delta ^{(2)}(p_nk)|0\left|\varphi (0)\right|n|^2$$
(2.2)
depends only on the invariant mass squared, $`m^2p_n^2`$, of the complete set of states $`|n`$ and has support for non-negative $`m^2`$. The positivity of the metric of the Hilbert space implies the positivity of the spectral function:

$$\rho (m^2)0.$$
(2.3)
The Jordan-Pauli invariant function $`\mathrm{\Delta }(xy;m^2)`$ is normalized according to

$$\mathrm{\Delta }(x,m^2)|_{x^0=0}=0,\frac{\mathrm{\Delta }(x,m^2)}{x^0}|_{x^0=0}=\delta (x^1).$$
(2.4)
Now, the canonical commutation relations for the unrenormalized Heisenberg operators,

$$[\varphi (t,x^1),\varphi (t,y^1)]=0,[\varphi (t,x^1),\dot{\varphi }(t,y^1)]=i\delta (x^1y^1),$$
(2.5)
imply a sumrule for the spectral function. In fact, differentiating (2.1) with respect to $`y^0`$ and taking the equal-time limit, we obtain

$$_0^{\mathrm{}}\rho (m^2)𝑑m^2=1.$$
(2.6)
The spectral function $`\rho (m^2)`$ consists of an isolated delta function contribution from the one-particle state and a continuous spectrum contribution from scattering states of multiple particles:

$$\rho (m^2)=Z\delta (m^2M^2)+\sigma (m^2).$$
(2.7)
Here, the wave function renormalization constant $`Z`$ represents the probability of finding a bare particle in the dressed particle state, while $`\sigma (m^2)`$ represents the continuous spectrum of multiple particles. The sum rule (2.6) for $`\rho (m^2)`$ reads

$$Z+\sigma (m^2)𝑑m^2=1.$$
(2.8)
Since the positivity (2.3) of $`\rho (m^2)`$ implies the positivity of both $`Z0`$ and $`\sigma (m^2)0`$, we find

$$1Z=\sigma (m^2)𝑑m^20.$$
(2.9)
The desired unitarity bound for the renormalization constant follows from this inequality and the positivity of $`Z`$:

$$0Z1.$$
(2.10)
From this unitarity bound, we usually conclude that the anomalous dimension of a field $`\varphi `$ is positive in the one-loop approximation, as shown in the introduction.

## 3 Wave function renormalization constant in NL$`\sigma `$Ms in two dimensions

Let us mow discuss the unitarity bound for the wave function renormalization constant in NL$`\sigma `$Ms described by the Lagrangian

$$=\frac{1}{2}g_{ij}(\phi )_\mu \phi ^i^\mu \phi ^j.$$
(3.1)
Here, we assume that the metric of the target space is positive definite to avoid the appearance of negative norm states. Furthermore, we do not incorporate higher-order derivative terms in our Lagrangian, since higher derivative terms inevitably introduce negative norm states. In the example discussed in the introduction, the target space is the coset space $`S^N=SO(N+1)/SO(N)`$, parameterized by $`\stackrel{}{\phi }=(\phi ^1,\mathrm{},\phi ^N)`$, and the metric of the target space is given by

$$g_{ij}=\frac{\delta _{ij}}{\left(1+\frac{g^2}{4}\stackrel{}{\phi }^2\right)^2}.(i,j=1,2,\mathrm{},N)$$
(3.2)
In general, a NL$`\sigma `$M is defined on a coset space $`G/H`$. It describes the interaction of massless Nambu-Goldstone bosons that appears when the symmetry group $`G`$ breaks down to its subgroup $`HG`$. While the global symmetry $`G`$ is realized nonlinearly, its subgroup $`H`$ is realized linearly. We further assume that the fields $`\stackrel{}{\phi }=(\phi ^1,\mathrm{},\phi ^N)`$ belong to some irreducible representation of the subgroup $`H`$.

With these assumptions, the Källén-Lehmann representation can be generalized to 333The vacuum is invariant under the unbroken symmetry group $`H`$.

$$0\left|[\phi ^i(x),\phi ^j(y)]\right|0=\delta ^{ij}_0^{\mathrm{}}\rho (m^2)i\mathrm{\Delta }(xy;m^2)𝑑m^2,$$
(3.3)
where the right-hand side is proportional to the $`H`$-invariant symmetric second-rank tensor $`\delta ^{ij}`$, because the fields $`\stackrel{}{\phi }`$ belong to an irreducible representation of $`H`$. If this were not the case, we would have several independent spectral functions for various representations.

The essential ingredient for the unitarity bound of the wave function renormalization constant is the positivity of the Hilbert space, as well as the canonical commutation relation (2.5). In NL$`\sigma `$Ms, the derivative interaction term changes the definition of the canonical momentum, so that the equal-time commutation relation (2.5) is no longer valid. The canonical commutation relation in this case reads

$$[\phi ^i(t,x^1),g_{jk}\dot{\phi }^k(t,y^1)]=i\delta _j^i\delta (x^1y^1),$$
(3.4)
or

$$[\phi ^i(t,x^1),\dot{\phi }^j(t,y^1)]=ig^{ij}\delta (x^1y^1),$$
(3.5)
where we have used the identity

$$g^{ik}g_{kj}=\delta _j^i.$$
(3.6)
Taking the vacuum expectation value of (3.5), we find

$$0\left|[\phi ^i(t,x^1),\dot{\phi }^j(t,y^1)]\right|0=i0\left|g^{ij}\right|0\delta (x^1y^1).$$
(3.7)
Because of the $`H`$-invariance of the vacuum, the right-hand side of this relation must be proportional to $`\delta ^{ij}`$, and thus we have

$$0\left|g^{ij}(\phi )\right|0=\delta ^{ij}(1+B)$$
(3.8)
for some quantity $`B`$. In the example of $`S^N`$, defined by Eq.(3.2), $`B`$ is given by

$$B(S^N)=0\left|\left(1+\frac{g^2}{4}\stackrel{}{\phi }^2\right)^2\right|010.$$
(3.9)
Differentiating (3.3) with respect to $`y^0`$ and taking the equal-time limit, we obtain

$$0\left|[\phi ^i(t,x^1),\dot{\phi }^j(t,y^1)]\right|0=i\delta ^{ij}\delta (x^1y^1)\rho (m^2)𝑑m^2.$$
(3.10)
Comparing this result with Eqs.(3.7) and (3.8), we find the sum-rule of the spectral function,

$$\rho (m^2)𝑑m^2=1+B.$$
(3.11)
The expression (2.7) of the spectral function as a sum of contributions from a single-particle and multi-particle states leads to

$`Z`$ $`=`$ $`1+BC1+B,`$ (3.12)
$`C`$ $`=`$ $`{\displaystyle \sigma (m^2)𝑑m^2}0.`$ (3.13)
We find from (3.12) that when there are no derivative interactions, $`B=0`$ implies the unitarity bound $`0Z1`$. In NL$`\sigma `$Ms, however, the wave function renormalization constant $`Z`$ can be greater than $`1`$, since $`B`$ may be positive, as shown in the example (3.9). Therefore, more generally, we have

$$0Z1+B.$$
(3.14)
The quantity $`B`$ comes from the matrix elements of local operators, and thus it represents the contributions of bubble diagrams, while $`\sigma (m^2)`$ gives the positive contributions of multi-particle states.

As an example, the lowest order contributions to $`B`$ and $`C`$ are shown in Figs.1 and 2. In the one-loop approximation, only $`B`$ in (3.9) receives a nonvanishing contribution, 444We use the Euclidean metric after the Wick rotation.

$$B=\frac{g^2}{2}\frac{d^2k}{(2\pi )^2}\frac{N}{k^2}=\frac{Ng^2}{8\pi }\mathrm{log}\frac{\mathrm{\Lambda }^2}{\mu ^2},$$
(3.15)
from the bubble diagram in Fig.1, while $`C`$ vanishes in this approximation. This value of $`B`$ and Eqs.(1.2) and (3.12) reproduce the negative one-loop anomalous dimension (1.5) given in the introduction.

It is worthwhile mentioning that the origin of the bound (3.14) is the invalid normalization (3.10) of field variables. In fact, if we define a rescaled field $`\stackrel{~}{\phi }`$ by

$$\phi ^i(x)\sqrt{1+B}\stackrel{~}{\phi }^i(x),$$
(3.16)
$`\stackrel{~}{\phi }(x)`$ satisfies the equal-time commutation relation

$$0\left|[\stackrel{~}{\phi }^i(t,x^1),\dot{\stackrel{~}{\phi }}^j(t,y^1)]\right|0=i\delta ^{ij}\delta (x^1y^1).$$
(3.17)
The spectral function $`\stackrel{~}{\rho }(m^2)`$ for the rescaled field $`\stackrel{~}{\phi }^i(x)`$ satisfies the normal sumrule (2.6). Therefore, we have the usual unitarity bound,

$$0\stackrel{~}{Z}1,$$
(3.18)
for the rescaled wave function renormalization constant $`\stackrel{~}{Z}`$, defined by

$$\stackrel{~}{\rho }(m^2)=\stackrel{~}{Z}\delta (m^2M^2)+\stackrel{~}{\sigma }(m^2).$$
(3.19)
The rescaling (3.16) has to be carried out order by order in perturbation theory, because $`B`$ depends on the coupling constant. The total renormalization constant $`Z`$ is proportional to $`\stackrel{~}{Z}`$,

$$Z=\stackrel{~}{Z}(1+B),$$
(3.20)
and satisfies the bound (3.14).

## 4 Coordinate dependence of the anomalous dimension

The NL$`\sigma `$M is a field theory in which field variables take values on a curved manifold, like $`S^N`$. To describe the curved manifold, we have the freedom of choosing any coordinates we wish. In fact we used a stereographic projection to express $`S^N`$ in the coordinates of $`N`$-dimensional Euclidean space. Instead, we could use other parameterizations, and in this section, we use the simplest coordinates ($`\varphi ^1,\varphi ^2,\mathrm{},\varphi ^{N+1}`$) to express $`S^N`$, defined by

$$(\varphi ^1)^2+(\varphi ^2)^2+\mathrm{}+(\varphi ^{N+1})^2=\frac{1}{g^2}.$$
(4.1)
If we express $`\varphi ^{N+1}`$ in terms of the other variables, $`\stackrel{}{\varphi }=(\varphi ^1,\varphi ^2,\mathrm{},\varphi ^N)`$, we have

$$\varphi ^{N+1}=\pm \sqrt{\frac{1}{g^2}\stackrel{}{\varphi }^2},$$
(4.2)
and the line element is given by

$$ds^2=\underset{i=1}{\overset{N+1}{}}(d\varphi ^i)^2=\underset{i,j=1}{\overset{N}{}}g_{ij}d\varphi ^id\varphi ^j,$$
(4.3)
where the metric in these coordinates reads

$$g_{ij}=\delta _{ij}+\frac{\varphi ^i\varphi ^j}{(\varphi ^{N+1})^2}.$$
(4.4)
Because the inverse of this metric is

$$g^{ij}=\delta ^{ij}g^2\varphi ^i\varphi ^j,$$
(4.5)
the one-loop contribution to $`B`$ defined by (3.8) is negative:

$$B=g^2\frac{d^2k}{(2\pi )^2}\frac{1}{k^2}=\frac{g^2}{4\pi }\mathrm{log}\frac{\mathrm{\Lambda }^2}{\mu ^2}.$$
(4.6)
This gives a positive anomalous dimension,

$$\gamma _\varphi =\frac{g^2}{4\pi },$$
(4.7)
in contrast to the negative anomalous dimension given in (1.5).

The sign of the anomalous dimension of the NL$`\sigma `$M depends on the parameterization of the manifold. How should we understand this fact? The anomalous dimension is an off-shell quantity describing the asymptotic behavior of the propagator in the ultraviolet region. Various parameterizations of a manifold give different metrics defining different field theories. There must be common quantities in these theories, since they all describe the same target manifold. In fact, under the coordinate transformation from $`\stackrel{}{\phi }`$ with the metric (3.2) to $`\stackrel{}{\varphi }`$ with the metric (4.4),

$$\stackrel{}{\varphi }=\frac{\stackrel{}{\phi }}{1+\frac{g^2}{4}\stackrel{}{\phi }^2},$$
(4.8)
physical quantities defined on the mass-shell are invariant. This is a special case of the theorem due to Kamefuchi, O’Raifeartaigh and Salam, which states that S-matrix elements are invariant under transformations of the type

$$\varphi ^i(x)=\phi ^i(x)+\underset{jk}{}a_{jk}^i\phi ^j(x)\phi ^k(x)+\underset{jkl}{}a_{jkl}^i\phi ^j(x)\phi ^k(x)\phi ^l(x)+\mathrm{}.$$
(4.9)
S-matrix elements are defined as pole residues of Green’s functions. Only the linear term in the transformation (4.9) contributes to the pole residues. Quadratic and higher-order terms do not contribute to poles unless there are bound states. Therefore, all S-matrix elements, that is, all physical quantities, are invariant under the transformation (4.9).

In a previous work , we derived the Wilsonian renormalization group equation in two-dimensional NL$`\sigma `$Ms with $`𝒩=2`$ supersymmetry. We used the so-called Kähler normal coordinate  to derive the anomalous dimensions. Because the anomalous dimensions depend on the choice of the coordinate, we should interpret the results obtained in Ref as the anomalous dimensions in the Kähler normal coordinate system.

## 5 Summary

The wave function renormalization constant $`Z`$ represents the probability of finding the bare particle in the physical particle, in the case that there are no derivative interactions. It satisfies the unitarity bound $`0Z1`$ in field theories without indefinite metric states. This unitarity bound implies the positivity of the anomalous dimension of the field in the one-loop approximation.

In NL$`\sigma `$Ms, however, this bound does not hold, due to the presence of the derivative interactions. The field dependence of the canonical momentum modifies the canonical commutation relation to

$$0\left|[\phi ^i(t,x^1),\dot{\phi }^j(t,y^1)]\right|0=i\delta ^{ij}\delta (x^1y^1)(1+B).$$
(5.1)
In this case, the actual unitarity bound of the wave function renormalization constant is

$$0Z1+B.$$

The contribution of the bubble diagrams to the self-energy diagrams, $`B`$, can be either positive or negative, while the contributions from multi-particle states, $`C`$, satisfies positivity, as expected. When $`B`$ is positive, the anomalous dimension in the one-loop approximation becomes negative. The sign of the contribution from the bubble diagrams, $`B`$, depends on the choice of the coordinates parameterizing the target manifold.

It is convenient to introduce the properly normalized field $`\stackrel{~}{\phi }`$ by (3.16). The renormalization constant $`\stackrel{~}{Z}`$ of the field $`\stackrel{~}{\phi }`$ satisfies the usual unitarity bound (3.18) and has a probabilistic interpretation.

## 6 Acknowledgements

This work was supported in part by Grants-in-Aid for Scientific Research (#13640283 and #13135215) and a Sasakawa Scientific Research Grant from The Japan Science Society. We would like to thank Kazuhiko Nishijima, Thomas E. Clark and Muneto Nitta for useful communications.