Table Of ContentDriven diffusive systems with mutually interactive Langmuir kinetics
H. D. Vuijk, R. Rens, M. Vahabi, F. C. MacKintosh, and A. Sharma
Department of Physics and Astronomy, VU University, Amsterdam, The Netherlands
(Dated: January 7, 2015)
Weinvestigatethesimpleone-dimensionaldrivenmodel,thetotallyasymmetricexclusionprocess,
coupled to mutually interactive Langmuir kinetics. This model is motivated by recent studies on
clustering of motor proteins on microtubules. In the proposed model, the attachment and detach-
mentratesofaparticlearemodifieddependingupontheoccupancyofneighbouringsites. Wefirst
obtain continuum mean-field equations and in certain limiting cases obtain analytic solutions. We
showhowmutualinteractionsincrease(decrease)theeffectsofboundariesonthephasebehaviorof
the model. We perform Monte Carlo simulations and demonstrate that our analytical approxima-
tions are in good agreement with the numerics over a wide range of model parameters. We present
5
phase diagrams over a selective range of parameters.
1
0
PACSnumbers: 87.16.Uv,05.70.Ln,87.10.Hk
2
n
a I. INTRODUCTION cluded in a study of motor transport is that of mutual
J interactions (MI) between motors. Seitz et al[14] ob-
6 servedthatinpresenceofanobstacle,amolecularmotor
Driven diffusive systems show a very rich behavior.
walking on a microtubule tends to stay attached for a
Evenone-dimensionalsystemsexhibitboundaryinduced
] longer time. Muto et al[15] reported on long-range co-
h phase transitions[1–6] with a complex phase behavior.
p Onesuchmodelisthatofatotallyasymmetricexclusion operative binding of kinesin to a microtubule. The de-
- tachment could depend on the biochemical state of the
process (TASEP)[7–10] coupled to the Langmuir kinet-
o
motor[16, 17] which might itself be determined by the
i ics (LK)[6]. In that model a single species of particles
b presence or absence of neighbouring motors.
performs unidirectional hopping on a 1D lattice. The
.
s particles are assumed to have hard-core repulsion which In a recent study on kinesin-1 motors moving on
c
preventsmorethanoneparticlefromoccupyingthesame microtubule[18], the authors performed numerical simu-
i
s lattice site. Such a system is coupled to Langmuir kinet- lationsofbinding/unbindingdynamicsincorporatingmu-
y
icsbyallowingforadsorption(desorption)ofparticlesat tually attractive interaction between the motor proteins.
h
an empty (filled) lattice site with fixed respective kinetic Their results were in agreement with the experimental
p
[ rates. It was shown in Refs. [5, 6] that the combina- observation, in particular clustering of motors on micro-
tionofTASEPandLKresultsinnonconserveddynamics tubules. Mututal interactions in addition to the hard-
1 with unusual features such as the appearnce of a high- core repulsion introduce additional correlations as in the
v
lowcoexistencephaseseparatedbystablediscontinuities Katz-Lebowitz-Spohn (KLS) model, which is a generic
2
7 in the density profile. The novel phase behavior has its model of interacting driven diffusive systems[19]. By
2 origininthecompetingkineticsofTASEPandLK.How- modifying the hopping rate of particles depending upon
1 ever, in the thermodynamic limit, it is expected that the theoccupancyofnextnearestneighbour,themodelgives
0 bulk effects are predominant with boundaries becoming risetoexoticfeaturessuchaslocalizeddownwardsshocks
1. insignificant. In fact, the competition between bulk and and phase separation into three distinct regimes[20].
0 boundary dynamics can occur only if one rescales the However,inthecaseofmolecularmotors,duetothemu-
5 attachment (detachment) kinetic rates [5, 6] such that tual interactions, the attachment and detachment rates
1 they decrease with increasing system size in a particular of a motor molecule are modified depending upon the
v: fashion. state of the neighbouring sites[18]. Assuming that the
i Besides being fundamentally interesting, understand- hopping rate is unaltered, this corresponds to the ordi-
X
nary TASEP (with no correlations besides the hard-core
ing nonequilibrium physics of driven systems is of par-
r repulsion) with density (local) dependent LK.
a ticular interest in biological systems[6, 11]. One such
particularsystemisthatofmolecularmotorsperforming In this paper, we focus on the TASEP coupled to mu-
directed motion along one-dimensional molecular tracks. tuallyinteractiveLangmuirkinetics. Weinvestigatehow
Typically kinetic rates are such that the fraction of mutual interactions can tilt the balance in favor of pre-
track over which the motor moves before detaching is dominantly bulk effects by enhancing LK. We show that
finite[12]. This allows for the bulk dynamics to com- that this is indeed the case when both the attachment
petewiththeboundary,potentiallygivingrisetounusual and detachment rates are enhanced significantly due to
nonequilibrium stationary states. Recently, exclusion the mutual interactions. In the case of the kinetic rates
processonnetworkshavebeenusedtomodelcytoskeletal being significantly reduced due to the interactions, one
transport[13]. It was shown that active transport pro- suppresses the bulk effects giving rise to rich and com-
cesses spontaneously develop density heterogeneities at plex phase behavior. We also explore the more interest-
various scales. An important aspect that needs to be in- ingscenarioinwhichthekineticratesaremodifiedinan
2
Figure 1: A graphic representation of the TASEP model
with Langmuir dynamics. Particles are injected at the first
sitewithrateαandextractedatthelastsitewithrateβ.The
particles move exclusively to the right, with unitary rate, if Figure 2: The mutual interaction are incorporated in the
the next site is empty. In this model the only interaction TASEP model with LK dynamics by modifying the attach-
between the particles is the hard-core repulsion, this means mentanddetachmentratesifneighbouringsitesareoccupied.
that hopping over particles and multiple occupation are not For each occupied neighbouring site the attachment rate is
allowed. The injection, extraction and hopping to the next multiplied by δ and the detachment rate is multiplied by γ.
site constitute the TASEP model. The Langmuir dynamics
consists of the detachment and attachment from and to the
background. The attachment and detachment rates are ωA Alltheratesaredefinedsuchthatthehoppingrateisuni-
and ωD respectively. For analysis of this model see for ex- tary. Processes a and b are the dynamics of the TASEP
ample Refs. [5, 6]. Mutual interactions are incorporated by model, process c is the interaction with the background.
modifying the attachment and detachment rates; see Fig. 2.
The interaction with the background is called Langmuir
kinetics. We note that the only interaction between the
particles is assumed to be the hard-core repulsion.
asymmetric fashion. The paper is organized in the fol-
For the sake of completion we show the equations for
lowing way. In Sec. II, we present the model composed
site-occupancy below. These equations are the same as
of TASEP coupled to the modified Langmuir Kinetics.
reported in Ref. [5]. The equation for the occupancy of
We first present the modfication to the LK due to the
each site is given by
mutual interactions. We then obtain continuum mean
field equation describing the steady state density profile
d
of particles on the lattice. In Sec. III we study three dif- dtni(t)=ni−1(t)[1−ni(t)]−ni(t)[1−ni+1(t)]
ferentcasesofmodifiedLK.Wefirstconsiderthecasein
+ω [1−n (t)]−ω n (t), (1)
which the unmodified LK rates are assumed to be equal A i D i
andmutualinteractionenhances(suppresses)theattach-
with n the occupancy at site i, which can either be one
ment and detachment rates by the same factor. The sec- i
or zero. The equations for the boundary sites are
ond case corresponds to the unmodified LK rates being
equal, but the mutual interaction enhances one and sup- d
n (t)=α[1−n (t)]−n (t)[1−n (t)], (2)
pressestheotherbythesamefactor. Thelastcaseisthe dt 1 1 1 2
most general one in which all the model parameters are
freely chosen. We do not explore this case in detail. In
d
Sec. IV, we summarize our findings.
n (t)=n (t)[1−n (t)]−βn (t). (3)
dt N N−1 N N
The mutual interactions of the particles are included
II. THE MODEL intheequationbymodifyingtheattachmentanddetach-
mentratesω andω respectively. Theattachment(de-
A D
The model consists of a 1D lattice with sites i = tachment)rateifbothneighbouringsitesareunoccupied
1,2,...,N; see Fig. 1. Each site on the lattice can ei- is ωA (ωD). If either the left or right neighbouring site is
ther be occupied by one particle or no particle. There occupiedtheattachment(detachment)ratebecomesδωA
are three different sub-processes that govern the dynam- (γωD), and if both neighbouring sites are occupied δ2ωA
ics of the system: (γ2ωD);seeFig.2. Inourmodel,thehoppingrateofpar-
ticlesonthelatticeisunalteredinpresenceofmutualin-
(a) Particles are injected at site 1 with rate α and ex- teractions. ThisisincontrastwiththeKLSmodelwhere
tracted at site N with rate β. the hopping rates are modified according to the occu-
pancyofnearestandnext-nearestneighbours[19]whereas
(b) Particles at site i = 1,...N −1 can hop to site i+1 the binding/unbinding kinetics remain unaltered.
if site i+1 is unoccupied. In order to include the mutual interactions of the par-
ticles the following substitutions are needed:
(c) Particles at site i=2,...,N −1 can detach from the
latticewithrateω andattachtositei=2,...,N−1 ω →ω [1+(n +n )(δ−1)+n n (δ2−2δ+1)],
D A A i+1 i−1 i+1 i−1
with rate ω . (4)
A
3
ω →ω [1+(n +n )(γ−1)+n n (γ2−2γ+1)]. other is obtained by matching the currents j(x) [6].
D D i+1 i−1 i+1 i−1
(5)
j(x)=ρ(x)[1−ρ(x)] (9)
At present it is not clear how already bound mo-
tors modify the binding kinetics of motors. It has For a normalized lattice in the continuous limit the
been suggested that presence of a bound motor could crossover region is localized and a discontinuity in the
change the lattice locally somehow leading to a modified density profile appears. Though the crossover region in
binding/unbindingkinetics[21]. However,theunderlying this case is localized it does span a finite number of sites
mechanism remains unkown. implying that in the case of a finite sized lattice the
In order to obtain useful solutions for the distribution crossover region spans a finite fraction of the normalized
of particles on the lattice, the following two steps are lattice.
required. First one goes from the equation for the occu- ThemodelwithoutMIexhibitsaparticle-holesymme-
pation of the sites, where each site can have either value try, in the sense that a particle attaching to the lattice
one or zero, to an equation of the average occupation of means that a vacancy detaches from the lattice and vice
thesites. Secondinthelimitoflargesystemsizesasemi- versa. The same holds for a particle entering (leaving)
continuous variable x instead of the discrete parameter the system on the first (last) site, which can be seen as
i is used for the position on the lattice. This method is a vacancy leaving (entering) the system. And a parti-
the same as used in [6]. The average density at a site is cle hooping to a neighbouring site on the right equals a
definedas(cid:104)ni(t)(cid:105)≡ρi(t). Instationarystatetheaverage vacancy hopping to the left. Due to this symmetry the
(cid:104)ni(cid:105) can either be a time or a sample average. In order following transformations
to take the averages of Eqs. (1, 2, 3) with substitutions
n (t)↔1−n (t), (10a)
(4), (5 the higher order correlations are needed. Instead i N−i
of solving these equations exactly a mean field approach α↔β, (10b)
is used which consists of the approximation: ω ↔ω , (10c)
A D
(cid:104)n (t)n (t)(cid:105)≈(cid:104)n (t)(cid:105)(cid:104)n (t)(cid:105). (6) leave Eqs. (1), (2) and (3) invariant [6]. However this
i i+1 i i+1
particle-hole symmetry is no longer apparent if MI is in-
The lattice constant (cid:15) is defined as (cid:15) ≡ L/N. For cluded.
simplicity the length of the lattice is fixed to one L =
1. For large system sizes, N (cid:29) 1 the quasicontinuous
positionvariablex=i/N isintroduced. Thismeansthat III. MODIFIED LANGMUIR KINETICS
the average density at site i is now defined as (cid:104)n (t)(cid:105) ≡
i
ρ(x,t). The equation for the average density profile in In this section the solutions of Eqn. (7) in the contin-
stationary state, to leading order in (cid:15) becomes uum limit are presented for three different cases. The
first case corresponds to where LK rates are both en-
(cid:15)
0= ∂2ρ+(2ρ−1)∂ ρ+Ω [1+ρ(δ−1)]2(1−ρ) hanced or reduced simultaneously by the same amount,
2 x x A and second case to where Ω is enhanced while Ω is
A D
−ΩD[1+ρ(γ−1)]2ρ. (7) reduced by the same amount. The third case is the most
general one in which the attachment and detachment
The [1+ρ(δ−1)]2 and [1+ρ(γ−1)]2 parts of the equa- rates are independently modified. This case ix explored
tion are due to the mutual interactions. In equation 7 in least details. With these solutions the density profiles
the total detachment/attachment rates are used, defined are constructed and compared with Monte Carlo simu-
as ΩA = ωAN and ΩD = ωDN. The equations for the lations of the model. Besides the density profiles phase
boundarysites(Eqs.(2,3))becometheboundarycondi- diagrams are made which show the characteristics of the
tions. solutions for different values of α and β.
ρ(0)=α , ρ(1)=1−β. (8)
Case 1: mutual interaction with enhanced LK rates
Onecannowtakethecontinuouslimit,(cid:15)→0foranor-
malized lattice this means N → ∞. In order to ensure
The model simplifies significantly in the case that
thattheattachment/detachmentratesperunitlengthdo
Ω = Ω ≡ Ω and δ = γ ≡ 1+η. This means that
not become infinite, the total rates Ω and Ω are kept A D
A D theattachmentanddetachmentratesarebothmultiplied
constant. Inthecontinuouslimit,(cid:15)→0,thesecondorder
by 1+η for each occupied neighbouring site. Values of
differentialequation(7)becomesafirstorderdifferential
η < −1 result in negative rates, therefore η is restricted
equation, but the two boundary conditions remain. This
to values larger than -1. Positive η increases and neg-
means that the problem is overdetermined. However one
ative η decreases the LK dynamics if neighbouring sites
canfindsolutionstotheequationinthecontinuumlimit
areoccupied. InthiscaseEqn.(7)inthecontinuouslimit
that satisfy one of the two boundary conditions. The
becomes
full density profile is is constructed from the possible so-
lutions. The crossover position from one solution to an 0=(2ρ−1)(cid:0)∂ ρ−Ω[1+ρη]2(cid:1). (11)
x
4
The special case where η = 0, the case without mutual phase consists of the ρ solution on the left side and the
α
interactions,correspondstothesymmetriccaseanalysed ρ on the right side of the lattice. The phase diagrams
β
in[6]. Eqn.(11)hasthreesolutions. Aconstantsolution are constructed using the information of the domain of
ρ = 1/2 which is the Langmuir isotherm. The other each of the solutions. The detailed construction of the
l
solutions are ρ and ρ which obey the left and right phase diagram with η =0 is reported in Ref. [6].
α β
boundary condition respectively. The phase diagrams for nonzero η are shown in Fig. 3
The behaviour of the phase diagrams for changes in η
are similar to changing Ω. For increasing Ω the area
α+(1+ηα)Ωx
ρ = (12) occupied by the LD, LD-HD and the HD phase in the
α 1−(1+ηα)ηΩx
phasediagramdecreaseandandeventuallydisappear[6].
The key difference between changing Ω and η is that the
1−β+(η(1−β)+1)Ω(x−1) influence of η is stronger in the phases containing the
ρβ = 1−(η(1−β)+1)ηΩ(x−1) (13) HD phase. As seen in Fig. 3(b) and (d), for increasing
η the HD phase disappears quickly from the phase dia-
The full density profile ρ(x) is a combination of these gram, while the LD phases decreases slowly. This is due
solutions to the high probability of occupied neighbouring sites in
regions of high density. Changes in Ω and η do not ef-
ρα for 0≤x≤xα, fect the area of the MC phase, which is always confined
ρ(x)= ρ for x ≤x≤x , (14) to the upper right quarter of the phase diagram. If one
l α β
ρ for x ≤x≤1. keepsincreasingη eventuallyonlytheLD-MC,MC,MC-
β β
HDandLD-MC-HDphaseremaininthephasediagram.
Where x and x are obtained by equating the currents Further increasing η does not change the phase diagram
α β
of the solutions, j (x ) = j and j (x ) = j [6]. The any more, however the density profiles do change. For
α α l β β l
domainofthesolutionscanbeexplainedasanattraction η → ∞ the LD and HD parts of the density profile oc-
of the density to the Langmuir isotherm as one moves cupy an infinitely small domain on the boundaries of the
away from the boundary. As one moves away from the profile. ThismeansthatduetotheincreaseinLangmuir
boundary the influence of the bulk dynamics i.e. the dynamics the density on the whole lattice is equal to the
Langmuirkineticsbecomesmoredominantandtherefore Langmuir density ρl.
the density tries to reach the Langmuir isotherm ρ .
l
In the case that x ≥ x the constant solution ρ is
α β l
not part of the density profile and the profile becomes Density profiles
(cid:40)
ρ(x)= ρα for 0≤x≤xw, (15) With equations for ρα, ρβ and ρl the density profiles
ρ for x ≤x≤1, areconstructed. InFig.4thedensityprofilesforthefive
β w
different phases with Ω = 0.3 and η = 2.0 are shown
wherex isobtainedbymatchingthecurrentsofthetwo ,thesecorrespondtothephasediagraminFig.3(d). The
w
solutions, j (x )=j (x ) [6]. In this case a discontinu- firstthingtonoticeisthattheρ andρ solutionsarenot
α w β w α β
ity appears at x . straightlines,incontrasttothesolutionsforη =0which
w
are straight lines. For η >0 the ρ and ρ solutions are
α β
concave up with a positive slope. This can be explained
Phase diagrams by an increase in attraction to the Langmuir isotherm
as the density increases. For example in figure 4 (a),
There are seven characteristic density profiles, called if one moves away from the left boundary the density
phases. DependingonΩandη allorsomecanbepresent increases. Thisincreaseindensityleadstoanincreasein
in the phase diagram. We follow the same terminology mutual interactions. In the case of positive η this results
as in Ref. [6]. The simplest three are the high density in an increase in LK dynamics and therefore an increase
(HD), the low density (LD) and the maximum current in attraction to the Langmuir isotherm. This increase
(MC)phase. IntheHD(LD)phasethedensityishigher in attraction causes the slope to increase. The ρ and
α
(lower) than 1/2, and the density profile is given by the ρ solutions with η < 0 are concave down with positive
β
ρ (ρ ) solution. In the MC phase the density profile is slope, this can be explained with the same arguments as
β α
equal to 1/2 over the whole domain. This is called the in the case of η >0.
maximum current phase because the current is maximal The analytical solutions of the density profiles in the
for ρ = 1/2. The density profile in the MC phase does continuum limit are compared with Monte Carlo simula-
not depend on the boundary conditions α and β. tions of the model; see Fig. 4. Due the the finite size of
Due to the interaction with the background it is pos- thelatticeusedinthesimulationsonecanexpectcertain
sible for two or all three solutions to coexist in a density discrepanciesbetweentheanalyticalresultsandthesim-
profile. This happens in the LD-HD, LD-MC, MC-HD ulations. Ifthe ρ ortheρ solutionsarenotpartofthe
α β
and the LD-MC-HD phases. For example the LD-HD density profile, one or both of the boundary conditions
5
Figure 4: Density profiles for Ω = 0.3 and η = 2, this
Figure 3: Five different phase diagrams for Ω = 0.3 and correspondstothephasediagraminFig.3(d). Theboundary
different values of η: (a) η = 0, (b) η = 0.5, (c) η = −0.5, conditionsare(a)α=0.01β =0.01,(b)α=0.2β =0.2,(c)
(d) η = 2.0 and (e) η = −0.9. (a) Corresponds to the case α=0.8β =0.1,(d)α=0.8β =0.8and(e)α=0.1β =0.8.
withoutmutualinteraction. TheinfluenceofΩisanalysedin The solid lines are the analytical solutions for the density
[6]. From(b)and(d)itbecomesclearthatduetotheincrease profiles (Eqs. (14,15)), the constant solution ρ is included
l
in η the HD region disappears more quickly from the phase as a dash-dot line. The solid lines with noise are the result
diagramthanthetheLDphase. Thisisexplainedbythefact of the Monte Carlo simulations with a lattice of 1000 sites,
that the mutual interactions are most apparent in region of averaged over 2000 simulations. The analytical solutions for
high density. In cases (c) and (e) η is decreased. Again it Ω = 0.3 and the same boundary conditions, but without MI
becomesclearfromthesefiguresthatachangeinη hasmore (η =0) are included as a dotted line to emphasize the effect
influenceonthephasescontainingHDregionsthanonphases ofthemutualinteractions. (a)TheMIdochangethedensity
containing LD regions. profile significantly, but do not change the phase. (b) Due to
the MI the LK dynamics increases, and the profile changes
from a LD-HD phase for η = 0 to a LD-MC-HD phase. (c)
TheMIcauseaphasechangefromtheHDphasetotheMC-
HDphase,duetotheincreasedLKdynamics.(d)Thereisno
difference between the solution with or without MI. (e) The
MI induce a phase change from the LD phase to the LD-MC
phase, due to the increased LK dynamics. There are two
types of discrepancies between the analytical results and the
are not met. This happens in the LD, LD-MC, MC and
Monte Carlo simulations. Boundary layers are formed at the
the MC-HD phase. In these phases a so called boundary
leftboundaryin(c)attherightboundaryin(e)andatboth
layer forms where the boundary condition is not met [6].
boundariesin(d). Otherdiscrepanciesbetweentheanalytical
Aboundarylayerisadiscrepancybetweentheanalytical result and the simulations occur where there is a transition
result of the equation in the continuum limit (Eqn. (11)) between the ρ , the ρ and or the ρ solution. Both types of
α β l
andthesimulationresultsofthemodelwithafinitesized discrepanciescanbeexplainedbythefinitesizeofthelattice
lattice. This discrepancy occurs at the boundary where used in the simulations.
the boundary condition is not met, and will occupy an
increasingly small domain for an increasing number of
latticesitesusedinthesimulations. Onecanalsoexpect Case 2: mutual interactions with antisymmetric
a discrepancy between the analytical density profile and modified LK rates
thesimulationresultswherethereisacrossoverfromone
solutiontotheother. Intheanalyticaldensityprofilethe InthepreviouscasetheMIincreasedtheLKdynamics,
crossoverislocalizedonthescaleoftherescaledvariable whichisnottheattractiveinteractionasreportedin[18].
x. However if the lattice has a finite number of sites, In this case the model is analyzed for attractive interac-
thecrossoverwillspanafinitefractionofthenormalized tions. Again Ω and Ω are set equal, Ω = Ω ≡ Ω.
A D A D
lattice. But in this case the mutual attraction are incorporated
6
in an antisymmetric manner, δ is increased and γ is de- of the Langmuir isotherm in the case without MI. The
creased by a factor ψ, δ = 1+ψ and γ = 1−ψ, with density profile is "attracted" to this constant solution,
−1 < ψ < 1. Depending on whether ψ is positive or as explained in the previous section. The solution ρ is
α
negative the interactions between the particles is respec- stable only for α<1/2, and ρ for β ≤1/2 [6].
β
tivelyattractiveorrepulsive. InthiscaseEqn.(7)inthe The full density profile is constructed from the solu-
continuous limit becomes tionsobeyingtheleftandrightboundaryconditions,and
calculating x , the position of the transition from ρ to
0=(2ρ−1)∂ ρ+Ω[1+ρψ]2(1−ρ)−Ω[1−ρψ]2ρ. (16) w α
x ρ , by matching the currents of these solutions.
β
Inordertofindsolutionsforthedensityprofile,theequa-
tion is simplified by neglecting the terms of order ψ2. In
Phase diagrams
the next section the limit of this approximation is dis-
cussed. With this approximation Eqn. (16) simplifies to
Using the same method as in the previous case the
0=(2ρ−1)∂ ρ+Ω−2Ω(1−ψ)ρ. (17) phase diagrams are constructed from the information
x
about the domain of the solutions. There are four pos-
Thisequationhasthesameformastheoneforthemodel
sible phases, these have the same characteristics as the
withoutMIbutwithunequalΩ andΩ (Eqn.(18))[6],
A D phases of the model without MI but with Ω (cid:54)= Ω [6].
A D
Due to the similarity only a short explanation is given
0=(2ρ−1)∂ ρ+Ω −(Ω +Ω )ρ. (18)
x A D A here, for a more elaborate discussion of the phases one
Eqn. (18) exhibits a particle-hole symmetry can consult [6]. In the LD (HD) phase the full density
(Eqn. (10)), therefore this symmetry is also appar- profile is governed by ρα (ρβ); boundary layers appear
ent in Eqn. (17). However this is not a property of the at the right (left) end of the lattice. The condition for
model and is only apparent if the ψ2 terms in Eqn. (16) the LD phase is xw > 1 and α < 1/2, and for the HD
are neglected. The transformation phasetheconditionisxw <0andβ <1/2. TheMphase
occurs for β >1/2 and x <0. This phase is called the
w
ρ(x)→1−ρ(1−x), (19a) "Meissner" phase due to similarities with the Meissner
α→β, (19b) phase in super conducting materials [6]. In the M phase
the density profile is independent of both boundary con-
β →α, (19c)
ditions, therefore boundary layers occur at both ends of
Ω→Ω(1−2ψ), (19d)
the lattice. Because the solution is not stable for these
−ψ valuesofβ theprofileisgivenbyρ (1/2)[6]. Thisphase
ψ → , (19e) β
1−2ψ can be seen as the equivalent of the MC phase in case
without MI [6] or case 1, because a profile in the MC
leaves Eqn. (17) invariant. Using this transformation phase is also independent of the boundary conditions.
density profiles for −1<ψ <0 can be obtained from so- The LD-HD phase, where phase coexistence occurs, is
lutionstoEqn.(17)with0<ψ. Thereforetheanalysisis split in two regions. In the region β < 1/2 the profile
restrictedtopositivevaluesofψ. SolutionstoEqn.(18), obeys both the left and right boundary condition. In
obtained by [6], are Lambert W functions. Using these the region β > 1/2 only the left boundary condition is
solutions one finds that the solutions to Eqn. (17) for obeyed. The right part of the density profile is inde-
ψ >0 are pendent of the right boundary condition and is given by
ρ (1/2),thisphaseisindicatedasLD-HD(M).Profilesin
ψ 1 β
ρ = (W [−y(x)]+1)+ for α<1/2, the LD-HD(M) phase have a boundary layer at the right
α 2(1−ψ) −1 2
end of the lattice. The conditions for the LD-HD phase
(20)
are 0 < x < 1, α < 1/2 and β < 1/2. For the LD-
w
(cid:40) ψ (W [y(x)]+1)+ 1 for 1−β ≥ρ , HD(M) phase the conditions are 0 < xw < 1, α < 1/2
ρ = 2(1−ψ) 0 2 l and β > 1/2. In Fig. 5 three phase diagrams are shown
β ψ (W [−y(x)]+1)+ 1 for 1 ≤1−β ≤ρ ,
2(1−ψ) 0 2 2 l for different values of ψ.
(21) Because Eqs. (20) and (21) are derived using the ap-
where ρα obeys the left and ρβ the right boundary con- proximation ψ2 = 0, the phase diagrams are also an ap-
dition. y(x) Is given by proximation which hold in the limit of small ψ.
(cid:12) (cid:12)
(cid:12)1−ψ (cid:12)
y(x)=(cid:12)(cid:12) ψ (2ρ0−1)−1(cid:12)(cid:12) Density profiles
(cid:20) (1−ψ)2 1−ψ (cid:21)
exp 2Ω (x−x )+ (2ρ −1)−1 , (22)
ψ 0 ψ 0 Using Eqs. (20) and (21) the density profiles can be
constructed. The domain of each of the solutions is
with ρ0 = α, x0 = 0 for ρα and ρ0 = 1−β, x0 = 1 for determined by matching the currents of the solutions,
ρβ. The constant solution ρl = 2(11−ψ) is the equivalent Eqn.(9). InFig.6,fivedensityprofilesaredepicted,one
7
Figure 5: Three phase diagrams obtained with Eqs. (20)
and (21). With Ω = 0.3 and (a) ψ = 0.001, (b) ψ = 0.4, (c)
ψ=0.8. Thephaseswhichcontaintheρ solutiondisappear
α
from the phase diagram for increasing MI, until only the M
andHDphasesareleft. TheareaoftheMphaseinthephase Figure 6: Density profiles for Ω = 0.3 and ψ = 0.4, this
diagramoccupiesanincreasinglylargepartoftheupperhalf corresponds to the phase diagram in Fig. 5 b. The injec-
ofthephasediagramforincreasingψ. Thisstandsincontrast tion/extraction rates are (a) α =0.01, β = 0.7 (LD phase),
withtheMCphaseincase1,wheretheMCphaseisconfined (b)α=0.1,β =0.7(LD-HD(M)phase),(c)α=0.8,β =0.8
totheupperrightquarterofthephasediagramandthearea (M phase), (d) α = 0.05, β = 0.1 (LD-HD phase) and (e)
is independent of any parameter; see Fig. 3. α = 0.3, β = 0.2 (HD phase). The dashed lines are the
Lambert W function solutions to Eqn. (17) obeying the left
and right boundary conditions. The parts of these solutions
thatmakeupthedensityprofilearerepresentedassolidblack
each for a phase in phase diagram 5 (b) (Ω = 0.3,ψ =
lines. The dash-dot line is the constant solution. The ana-
0.4). It is clear from the Fig. 6 that there is good agree-
lytical solutions for the same values but without MI are in-
ment between the simulations and the analytical solu- cludedasdottedlinestoemphasizetheinfluenceoftheMIon
tions(Eqs.(20,21)). Theapparentdiscrepanciesbetween the density profiles. Solid lines with noise are the results of
theanalyticalandnumericalresultsarecausedbythefi- the Monte simulations with a lattice of 1000 sites, averaged
nite size of the lattice used in the simulations. Besides over 2000 simulations. Over all there is good agreement be-
this there is also a small discrepancy between the an- tweenthesimulationsandtheanalyticalresult. Therearetwo
alytical solution and the simulations caused by the ap- causes for the discrepancies, the finite size of the lattice and
proximationψ2 =0. Thiscancauseadiscrepancyinthe the approximation ψ2 ≈ 0. Due to the finite size of the lat-
ticeboundarylayersareformedattherightendof(a),(b),(d)
domainwallposition,ascanbeseeninFig.6(b)and(d).
and at the left end of (c) and (d); and the domain walls in
(b) and (d) are not localized. The discrepancies caused by
theapproximationarevisibleintwoways,thedensityprofile
Approximation limits does not fully coincide with the analytical result and due to
this the position ofthedomain wall is shifted. This is visible
in (b) and (d).
In deriving Eqs. (20) and (21) terms of the order
ψ2 were neglected. The neglected part of Eqn. (16) is
Ωψ2ρ2−2Ωψ2ρ3, which shows that every ψ2 is coupled
to either a ρ2 or a ρ3. This means that the break down Case 3: mutual interactions with arbitrarily
modified LK rates
of the approximation is governed by both ρ and ψ and
theapproximationholdforlowdensitiesregardlessofthe
value of ψ, because MI do not play a significant role in Until now, we have considered enhancement or sup-
low densities due to the low probability of having occu- pression of LK rates in a symmetric or antisymmetric
pied neighbouring sites. Fig. 7 illustrates some of the fashion. The most general case in which all the relevant
limits of the approximation. From Figs. 6(b),(d) and 7 parameters (Ω , Ω , δ, γ, α, β) are assigned randomly
A D
(b)itbecomesclearthatthedomainofthelowandhigh chosen values is extremely difficult to treat analytically.
densitysolutioncandiffersignificantlyeveniftheρ and Due to the large parameter-space (6-dimensional), one
α
ρ differ only little from the simulation result. cannot gain muchinsight by performing numerical simu-
β
8
Figure 7: The same legend as in Fig. 6 is used. In order Figure 8: The density profiles for (a) Ω = 0.5, Ω = 1,
A D
to illustrate the limits of the approximation ψ2 ≈ 0 used in δ = 2, γ = 0.1, α = 1 and β = 1, (b) Ω = 0.5, Ω = 0.5,
A D
deriving the equations for the density profiles, four extreme δ = 2, γ = 0.5, α = 0 and β = 1, (c) Ω = 0.5, Ω = 0.5,
A D
cases are shown. For all figures Ω = 0.3 was used an for (a) δ = 2, γ = 0.5, α = 0 and β = 0, (d) Ω = 0.5, Ω = 1,
A D
ψ=−0.99, α=0.1, β =0.1, (b) ψ=0.8, α=0.1, β =0.9, δ = 3, γ = 0.1, α = 0 and β = 0. Solid line with noise
(c) ψ = 0.5, α = 0.9, β = 0.1 and for (d) ψ = 0.9, α = is the result from the simulation, the dashed line is the so-
0.9, β =0.1. From(a)itisclearthattheapproximationholds lution without MI. In all cases the attachment/detachment
forlowdensities,regardlessofthevaluesofψ. Theanalytical is increased/decreased, which results in a higher density. In
solution in (b) does not fully coincide with the simulations. (b)and(c)thedensityoverlapswithasmallpartoftheana-
There are two causes for the discrepancies. First of all there lytical solution without MI. This can be explained by the in
isaboundarylayerontherightsidecausedbythefinitesize low density regions the effect of the MI is small. Though in
ofthelattice. Secondlythe ρβ solutionishigher thanthere- (d) the solution without MI the density is low, the lattice is
sultsfromthesimulation,thisiscausedbyneglectingtheψ2 almost completely filled due to the increase/decrease in at-
terms. Over all the simulations are in good agreement with tachment/detachment. All parameters in (a) and (d) are the
the analytical solutions, however the domain of the solutions sameexceptfortheboundaryconditions,whichpreventsthe
differssignificantly. TheanalyticalprofileisintheHDphase, lattice in (a) from filling up completely.
x < 0. The simulation result, on the other hand, is in the
w
LD-HD phase, 0 < x < 1. (c) For this value of ψ there
w
is agreement between the analytical solution and the simula-
IV. CONCLUSION
tions, even though density is high. (d) The analytical result
doesnotcoincidewiththesimulationsduetothecombination
ofhighdensityandaψ closeto1. Inthiscasethedensityof Inthiswork,weinvestigatethesimpleone-dimensional
theanalyticalresultexceeds1,whichisphysicallyimpossible. driven model– the totally asymmetric exclusion process,
coupled to a modified form of Langmuir kinetics. This
modelismotivatedbyrecentstudiesonclusteringofmo-
torproteinsonmicrotubules. Withoutaddressingtheun-
derlying mechanism, it is assumed that the attachment
and detachment rates of a particle depend on the occu-
pancy of the nearest neighbours of any given site. Ig-
noring density correlations, we obtain continuum mean-
lations. Therefore,wehavenotexploredthemostgeneral field equation describing the density profile on the lat-
case in any details. However, we consider few represen- tice. Imposing certain conditions, we obtain analytical
tative cases in which the choice of model parameters is solution to the equation and demonstrate using Monte
based on the observation that the resulting density pro- Carlo simulations that our analytical solutions are accu-
files have interesting features when contrasted with the rateoverawiderangeofparameters. Weshowthatwhen
case with no mutual interactions. In Fig. 8 we show bothattachmentanddetachmentratesareenhanceddue
profiles corresponding to four different sets of parame- to mutual interactions, bulk-effects start dominating the
ters. As can be seen the profiles look very different from phase behavior of the model. The two-phase coexistence
thosewithnomutualinteraction. However,asmentioned (low and high density) observed in absence of mutual in-
above, at present our analysis of the most general case is teractions can become three-phase coexistence (low and
very qualitative and highly superficial. It is obvious that highdensitywithmaximumcurrentphase)whenmutual
it requires much further investigation. We leave detailed interactions are attractive. On varying the mutual in-
analysis of the general case to a future study. teraction between particles (attractive or repulsive), we
9
obtain a very rich phase diagram describing the behav- Langmuir kinetics.
ior of driven diffusive system with mutually interactive
[1] J. Krug, Phys. Rev. Lett. 67, 1882 (1991). cytoskeleton (SinauerAssociatesSunderland,MA,2001).
[2] M. Evans, D. Foster, C. Godreche, and D. Mukamel, [13] I. Neri, N. Kern, and A. Parmeggiani, New Journal of
Physical review letters 74, 208 (1995). Physics 15, 085005 (2013).
[3] M.Evans,Y.Kafri,H.Koduvely,andD.Mukamel,Phys- [14] A. Seitz and T. Surrey, The EMBO journal 25, 267
ical review letters 80, 425 (1998). (2006).
[4] Y. Kafri, E. Levine, D. Mukamel, G. Schütz, and [15] E. Muto, H. Sakai, and K. Kaseda, The Journal of cell
J. Török, Physical review letters 89, 035702 (2002). biology 168, 691 (2005).
[5] A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. [16] I. M.-T. Crevel, M. Nyitrai, M. C. Alonso, S. Weiss,
Lett. 90, 086601 (2003). M. A. Geeves, and R. A. Cross, The EMBO journal 23,
[6] A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. E 23 (2004).
70, 046101 (2004). [17] I.A.Telley,P.Bieling,andT.Surrey,Biophysicaljournal
[7] H. Spohn, Large scale dynamics of interacting particles, 96, 3341 (2009).
vol. 825 (Springer, 1991). [18] W. H. Roos, O. Campàs, F. Montel, G. Woehlke, J. P.
[8] V.Privman,Nonequilibrium statistical mechanics in one Spatz, P. Bassereau, and G. Cappello, Physical biology
dimension (Cambridge University Press, 2005). 5, 046004 (2008).
[9] M. E. Cates and M. R. Evans, Soft and Fragile Mat- [19] S. Katz, J. Lebowitz, and H. Spohn, J. Stat. Phys 34,
ter: Nonequilibrium Dynamics, Metastability and Flow 497 (1987).
(PBK) (CRC Press, 2010). [20] V.Popkov,A.Rákos,R.D.Willmann,A.B.Kolomeisky,
[10] C. Domb, Phase transitions and critical phenomena, andG.M.Schütz,PhysicalReviewE67,066117(2003).
vol. 19 (Academic Press, 2000). [21] A. Krebs, K. N. Goldie, and A. Hoenger, Journal of
[11] C.MacDonald, J.Gibbs, andA.Pipkin, Biopolymers6, molecular biology 335, 139 (2004).
1 (1968).
[12] J. Howard et al., Mechanics of motor proteins and the
