Driven diffusive systems with mutually interactive Langmuir kinetics
Abstract
We investigate the simple onedimensional driven model, the totally asymmetric exclusion process, 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 detachment rates of a particle are modified depending upon the occupancy of neighbouring sites. We first obtain continuum meanfield equations and in certain limiting cases obtain analytic solutions. We show how mutual interactions increase (decrease) the effects of boundaries on the phase behavior of the model. We perform Monte Carlo simulations and demonstrate that our analytical approximations are in good agreement with the numerics over a wide range of model parameters. We present phase diagrams over a selective range of parameters.
pacs:
87.16.Uv, 05.70.Ln, 87.10.HkI introduction
Driven diffusive systems show a very rich behavior. Even onedimensional systems exhibit boundary induced phase transitionsKrug (1991); Evans et al. (1995, 1998); Kafri et al. (2002); Parmeggiani et al. (2003, 2004) with a complex phase behavior. One such model is that of a totally asymmetric exclusion process (TASEP)Spohn (1991); Privman (2005); Cates and Evans (2010); Domb (2000) coupled to the Langmuir kinetics (LK)Parmeggiani et al. (2004). In that model a single species of particles performs unidirectional hopping on a 1D lattice. The particles are assumed to have hardcore repulsion which prevents more than one particle from occupying the same lattice site. Such a system is coupled to Langmuir kinetics by allowing for adsorption (desorption) of particles at an empty (filled) lattice site with fixed respective kinetic rates. It was shown in Refs. Parmeggiani et al. (2003, 2004) that the combination of TASEP and LK results in nonconserved dynamics with unusual features such as the appearnce of a highlow coexistence phase separated by stable discontinuities in the density profile. The novel phase behavior has its origin in the competing kinetics of TASEP and LK. However, in the thermodynamic limit, it is expected that the bulk effects are predominant with boundaries becoming insignificant. In fact, the competition between bulk and boundary dynamics can occur only if one rescales the attachment (detachment) kinetic rates Parmeggiani et al. (2003, 2004) such that they decrease with increasing system size in a particular fashion.
Besides being fundamentally interesting, understanding nonequilibrium physics of driven systems is of particular interest in biological systemsMacDonald et al. (1968); Parmeggiani et al. (2004). One such particular system is that of molecular motors performing directed motion along onedimensional molecular tracks. Typically kinetic rates are such that the fraction of track over which the motor moves before detaching is finiteHoward et al. (2001). This allows for the bulk dynamics to compete with the boundary, potentially giving rise to unusual nonequilibrium stationary states. Recently, exclusion process on networks have been used to model cytoskeletal transportNeri et al. (2013). It was shown that active transport processes spontaneously develop density heterogeneities at various scales. An important aspect that needs to be included in a study of motor transport is that of mutual interactions (MI) between motors. Seitz et alSeitz and Surrey (2006) observed that in presence of an obstacle, a molecular motor walking on a microtubule tends to stay attached for a longer time. Muto et alMuto et al. (2005) reported on longrange cooperative binding of kinesin to a microtubule. The detachment could depend on the biochemical state of the motorCrevel et al. (2004); Telley et al. (2009) which might itself be determined by the presence or absence of neighbouring motors.
In a recent study on kinesin1 motors moving on microtubuleRoos et al. (2008), the authors performed numerical simulations of binding/unbinding dynamics incorporating mutually attractive interaction between the motor proteins. Their results were in agreement with the experimental observation, in particular clustering of motors on microtubules. Mututal interactions in addition to the hardcore repulsion introduce additional correlations as in the KatzLebowitzSpohn (KLS) model, which is a generic model of interacting driven diffusive systemsKatz et al. (1987). By modifying the hopping rate of particles depending upon the occupancy of next nearest neighbour, the model gives rise to exotic features such as localized downwards shocks and phase separation into three distinct regimesPopkov et al. (2003). However, in the case of molecular motors, due to the mutual interactions, the attachment and detachment rates of a motor molecule are modified depending upon the state of the neighbouring sitesRoos et al. (2008). Assuming that the hopping rate is unaltered, this corresponds to the ordinary TASEP (with no correlations besides the hardcore repulsion) with density (local) dependent LK.
In this paper, we focus on the TASEP coupled to mutually interactive Langmuir kinetics. We investigate how mutual interactions can tilt the balance in favor of predominantly bulk effects by enhancing LK. We show that that this is indeed the case when both the attachment and detachment rates are enhanced significantly due to the mutual interactions. In the case of the kinetic rates being significantly reduced due to the interactions, one suppresses the bulk effects giving rise to rich and complex phase behavior. We also explore the more interesting scenario in which the kinetic rates are modified in an asymmetric fashion. The paper is organized in the following way. In Sec. II, we present the model composed of TASEP coupled to the modified Langmuir Kinetics. We first present the modfication to the LK due to the mutual interactions. We then obtain continuum mean field equation describing the steady state density profile of particles on the lattice. In Sec. III we study three different cases of modified LK. We first consider the case in which the unmodified LK rates are assumed to be equal and mutual interaction enhances (suppresses) the attachment and detachment rates by the same factor. The second case corresponds to the unmodified LK rates being equal, but the mutual interaction enhances one and suppresses the other by the same factor. The last case is the most general one in which all the model parameters are freely chosen. We do not explore this case in detail. In Sec. IV, we summarize our findings.
Ii The Model
The model consists of a 1D lattice with sites ; see Fig. 1.
Each site on the lattice can either be occupied by one particle or no particle. There are three different subprocesses that govern the dynamics of the system:

Particles are injected at site with rate and extracted at site with rate .

Particles at site can hop to site if site is unoccupied.

Particles at site can detach from the lattice with rate and attach to site with rate .
All the rates are defined such that the hopping rate is unitary. Processes and are the dynamics of the TASEP model, process is the interaction with the background. The interaction with the background is called Langmuir kinetics. We note that the only interaction between the particles is assumed to be the hardcore repulsion.
For the sake of completion we show the equations for siteoccupancy below. These equations are the same as reported in Ref. Parmeggiani et al. (2003). The equation for the occupancy of each site is given by
(1) 
with the occupancy at site , which can either be one or zero. The equations for the boundary sites are
(2) 
(3) 
The mutual interactions of the particles are included in the equation by modifying the attachment and detachment rates and respectively. The attachment (detachment) rate if both neighbouring sites are unoccupied is (). If either the left or right neighbouring site is occupied the attachment (detachment) rate becomes (), and if both neighbouring sites are occupied (); see Fig. 2. In our model, the hopping rate of particles on the lattice is unaltered in presence of mutual interactions. This is in contrast with the KLS model where the hopping rates are modified according to the occupancy of nearest and nextnearest neighboursKatz et al. (1987) whereas the binding/unbinding kinetics remain unaltered.
In order to include the mutual interactions of the particles the following substitutions are needed:
(4) 
(5) 
At present it is not clear how already bound motors modify the binding kinetics of motors. It has been suggested that presence of a bound motor could change the lattice locally somehow leading to a modified binding/unbinding kineticsKrebs et al. (2004). However, the underlying mechanism remains unkown.
In order to obtain useful solutions for the distribution of particles on the lattice, the following two steps are required. First one goes from the equation for the occupation of the sites, where each site can have either value one or zero, to an equation of the average occupation of the sites. Second in the limit of large system sizes a semicontinuous variable instead of the discrete parameter is used for the position on the lattice. This method is the same as used in (Parmeggiani et al., 2004). The average density at a site is defined as . In stationary state the average can either be a time or a sample average. In order to take the averages of Eqs. (1, 2, 3) with substitutions (4), (5 the higher order correlations are needed. Instead of solving these equations exactly a mean field approach is used which consists of the approximation:
(6) 
The lattice constant is defined as . For simplicity the length of the lattice is fixed to one . For large system sizes, the quasicontinuous position variable is introduced. This means that the average density at site is now defined as . The equation for the average density profile in stationary state, to leading order in becomes
(7) 
The and parts of the equation are due to the mutual interactions. In equation 7 the total detachment/attachment rates are used, defined as and . The equations for the boundary sites (Eqs. (2, 3)) become the boundary conditions.
(8) 
One can now take the continuous limit, for a normalized lattice this means . In order to ensure that the attachment/detachment rates per unit length do not become infinite, the total rates and are kept constant. In the continuous limit, , the second order differential equation (7) becomes a first order differential equation, but the two boundary conditions remain. This means that the problem is overdetermined. However one can find solutions to the equation in the continuum limit that satisfy one of the two boundary conditions. The full density profile is is constructed from the possible solutions. The crossover position from one solution to an other is obtained by matching the currents (Parmeggiani et al., 2004).
(9) 
For a normalized lattice in the continuous limit the crossover region is localized and a discontinuity in the density profile appears. Though the crossover region in this case is localized it does span a finite number of sites implying that in the case of a finite sized lattice the crossover region spans a finite fraction of the normalized lattice.
The model without MI exhibits a particlehole symmetry, in the sense that a particle attaching to the lattice means that a vacancy detaches from the lattice and vice versa. The same holds for a particle entering (leaving) the system on the first (last) site, which can be seen as a vacancy leaving (entering) the system. And a particle hooping to a neighbouring site on the right equals a vacancy hopping to the left. Due to this symmetry the following transformations
(10a)  
(10b)  
(10c) 
leave Eqs. (1), (2) and (3) invariant (Parmeggiani et al., 2004). However this particlehole symmetry is no longer apparent if MI is included.
Iii Modified Langmuir Kinetics
In this section the solutions of Eqn. (7) in the continuum limit are presented for three different cases. The first case corresponds to where LK rates are both enhanced or reduced simultaneously by the same amount, and second case to where is enhanced while is reduced by the same amount. The third case is the most general one in which the attachment and detachment rates are independently modified. This case ix explored in least details. With these solutions the density profiles are constructed and compared with Monte Carlo simulations of the model. Besides the density profiles phase diagrams are made which show the characteristics of the solutions for different values of and .
Case 1: mutual interaction with enhanced LK rates
The model simplifies significantly in the case that and . This means that the attachment and detachment rates are both multiplied by for each occupied neighbouring site. Values of result in negative rates, therefore is restricted to values larger than 1. Positive increases and negative decreases the LK dynamics if neighbouring sites are occupied. In this case Eqn. (7) in the continuous limit becomes
(11) 
The special case where , the case without mutual interactions, corresponds to the symmetric case analysed in (Parmeggiani et al., 2004). Eqn. (11) has three solutions. A constant solution which is the Langmuir isotherm. The other solutions are and which obey the left and right boundary condition respectively.
(12) 
(13) 
The full density profile is a combination of these solutions
(14) 
Where and are obtained by equating the currents of the solutions, and (Parmeggiani et al., 2004). The domain of the solutions can be explained as an attraction of the density to the Langmuir isotherm as one moves away from the boundary. As one moves away from the boundary the influence of the bulk dynamics i.e. the Langmuir kinetics becomes more dominant and therefore the density tries to reach the Langmuir isotherm .
In the case that the constant solution is not part of the density profile and the profile becomes
(15) 
where is obtained by matching the currents of the two solutions, (Parmeggiani et al., 2004). In this case a discontinuity appears at .
Phase diagrams
There are seven characteristic density profiles, called phases. Depending on and all or some can be present in the phase diagram. We follow the same terminology as in Ref. (Parmeggiani et al., 2004). The simplest three are the high density (HD), the low density (LD) and the maximum current (MC) phase. In the HD (LD) phase the density is higher (lower) than , and the density profile is given by the () solution. In the MC phase the density profile is equal to over the whole domain. This is called the maximum current phase because the current is maximal for . The density profile in the MC phase does not depend on the boundary conditions and .
Due to the interaction with the background it is possible for two or all three solutions to coexist in a density profile. This happens in the LDHD, LDMC, MCHD and the LDMCHD phases. For example the LDHD phase consists of the solution on the left side and the on the right side of the lattice. The phase diagrams are constructed using the information of the domain of each of the solutions. The detailed construction of the phase diagram with is reported in Ref. (Parmeggiani et al., 2004).
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 occupied by the LD, LDHD and the HD phase in the phase diagram decrease and and eventually disappear (Parmeggiani et al., 2004). The key difference between changing and is that the influence of is stronger in the phases containing the HD phase. As seen in Fig. 3(b) and (d), for increasing the HD phase disappears quickly from the phase diagram, while the LD phases decreases slowly. This is due to the high probability of occupied neighbouring sites in regions of high density. Changes in and do not effect the area of the MC phase, which is always confined to the upper right quarter of the phase diagram. If one keeps increasing eventually only the LDMC, MC, MCHD and LDMCHD phase remain in the phase diagram. Further increasing does not change the phase diagram any more, however the density profiles do change. For the LD and HD parts of the density profile occupy an infinitely small domain on the boundaries of the profile. This means that due to the increase in Langmuir dynamics the density on the whole lattice is equal to the Langmuir density .
Density profiles
With equations for , and the density profiles are constructed. In Fig. 4 the density profiles for the five different phases with and are shown ,these correspond to the phase diagram in Fig. 3(d). The first thing to notice is that the and solutions are not straight lines, in contrast to the solutions for which are straight lines. For the and solutions are concave up with a positive slope. This can be explained by an increase in attraction to the Langmuir isotherm as the density increases. For example in figure 4 (a), if one moves away from the left boundary the density increases. This increase in density leads to an increase in mutual interactions. In the case of positive this results in an increase in LK dynamics and therefore an increase in attraction to the Langmuir isotherm. This increase in attraction causes the slope to increase. The and solutions with are concave down with positive slope, this can be explained with the same arguments as in the case of .
The analytical solutions of the density profiles in the continuum limit are compared with Monte Carlo simulations of the model; see Fig. 4. Due the the finite size of the lattice used in the simulations one can expect certain discrepancies between the analytical results and the simulations. If the or the solutions are not part of the density profile, one or both of the boundary conditions are not met. This happens in the LD, LDMC, MC and the MCHD phase. In these phases a so called boundary layer forms where the boundary condition is not met (Parmeggiani et al., 2004). A boundary layer is a discrepancy between the analytical result of the equation in the continuum limit (Eqn. (11)) and the simulation results of the model with a finite sized lattice. This discrepancy occurs at the boundary where the boundary condition is not met, and will occupy an increasingly small domain for an increasing number of lattice sites used in the simulations. One can also expect a discrepancy between the analytical density profile and the simulation results where there is a crossover from one solution to the other. In the analytical density profile the crossover is localized on the scale of the rescaled variable . However if the lattice has a finite number of sites, the crossover will span a finite fraction of the normalized lattice.
Case 2: mutual interactions with antisymmetric modified LK rates
In the previous case the MI increased the LK dynamics, which is not the attractive interaction as reported in (Roos et al., 2008). In this case the model is analyzed for attractive interactions. Again and are set equal, . But in this case the mutual attraction are incorporated in an antisymmetric manner, is increased and is decreased by a factor , and , with . Depending on whether is positive or negative the interactions between the particles is respectively attractive or repulsive. In this case Eqn. (7) in the continuous limit becomes
(16) 
In order to find solutions for the density profile, the equation is simplified by neglecting the terms of order . In the next section the limit of this approximation is discussed. With this approximation Eqn. (16) simplifies to
(17) 
This equation has the same form as the one for the model without MI but with unequal and (Eqn. (18)) (Parmeggiani et al., 2004),
(18) 
Eqn. (18) exhibits a particlehole symmetry (Eqn. (10)), therefore this symmetry is also apparent in Eqn. (17). However this is not a property of the model and is only apparent if the terms in Eqn. (16) are neglected. The transformation
(19a)  
(19b)  
(19c)  
(19d)  
(19e) 
leaves Eqn. (17) invariant. Using this transformation density profiles for can be obtained from solutions to Eqn. (17) with . Therefore the analysis is restricted to positive values of . Solutions to Eqn. (18), obtained by (Parmeggiani et al., 2004), are Lambert W functions. Using these solutions one finds that the solutions to Eqn. (17) for are
(20) 
(21) 
where obeys the left and the right boundary condition. Is given by
(22) 
with , for and , for . The constant solution is the equivalent of the Langmuir isotherm in the case without MI. The density profile is "attracted" to this constant solution, as explained in the previous section. The solution is stable only for , and for (Parmeggiani et al., 2004).
The full density profile is constructed from the solutions obeying the left and right boundary conditions, and calculating , the position of the transition from to , by matching the currents of these solutions.
Phase diagrams
Using the same method as in the previous case the phase diagrams are constructed from the information about the domain of the solutions. There are four possible phases, these have the same characteristics as the phases of the model without MI but with (Parmeggiani et al., 2004). Due to the similarity only a short explanation is given here, for a more elaborate discussion of the phases one can consult (Parmeggiani et al., 2004). In the LD (HD) phase the full density profile is governed by (); boundary layers appear at the right (left) end of the lattice. The condition for the LD phase is and , and for the HD phase the condition is and . The M phase occurs for and . This phase is called the "Meissner" phase due to similarities with the Meissner phase in super conducting materials (Parmeggiani et al., 2004). In the M phase the density profile is independent of both boundary conditions, therefore boundary layers occur at both ends of the lattice. Because the solution is not stable for these values of the profile is given by (Parmeggiani et al., 2004). This phase can be seen as the equivalent of the MC phase in case without MI (Parmeggiani et al., 2004) or case 1, because a profile in the MC phase is also independent of the boundary conditions. The LDHD phase, where phase coexistence occurs, is split in two regions. In the region the profile obeys both the left and right boundary condition. In the region only the left boundary condition is obeyed. The right part of the density profile is independent of the right boundary condition and is given by , this phase is indicated as LDHD(M). Profiles in the LDHD(M) phase have a boundary layer at the right end of the lattice. The conditions for the LDHD phase are , and . For the LDHD(M) phase the conditions are , and . In Fig. 5 three phase diagrams are shown for different values of .
Density profiles
Using Eqs. (20) and (21) the density profiles can be constructed. The domain of each of the solutions is determined by matching the currents of the solutions, Eqn. (9). In Fig. 6, five density profiles are depicted, one each for a phase in phase diagram 5 (b) (). It is clear from the Fig. 6 that there is good agreement between the simulations and the analytical solutions (Eqs. (20,21)). The apparent discrepancies between the analytical and numerical results are caused by the finite size of the lattice used in the simulations. Besides this there is also a small discrepancy between the analytical solution and the simulations caused by the approximation . This can cause a discrepancy in the domain wall position, as can be seen in Fig. 6(b) and (d).
Approximation limits
In deriving Eqs. (20) and (21) terms of the order were neglected. The neglected part of Eqn. (16) is , which shows that every is coupled to either a or a . This means that the break down of the approximation is governed by both and and the approximation hold for low densities regardless of the value of , because MI do not play a significant role in low densities due to the low probability of having occupied neighbouring sites. Fig. 7 illustrates some of the limits of the approximation. From Figs. 6(b),(d) and 7 (b) it becomes clear that the domain of the low and high density solution can differ significantly even if the and differ only little from the simulation result.
Case 3: mutual interactions with arbitrarily modified LK rates
Until now, we have considered enhancement or suppression of LK rates in a symmetric or antisymmetric fashion. The most general case in which all the relevant parameters (, , , , , ) are assigned randomly chosen values is extremely difficult to treat analytically. Due to the large parameterspace (6dimensional), one cannot gain much insight by performing numerical simulations. Therefore, we have not explored the most general case in any details. However, we consider few representative cases in which the choice of model parameters is based on the observation that the resulting density profiles have interesting features when contrasted with the case with no mutual interactions. In Fig. 8 we show profiles corresponding to four different sets of parameters. As can be seen the profiles look very different from those with no mutual interaction. However, as mentioned above, at present our analysis of the most general case is very qualitative and highly superficial. It is obvious that it requires much further investigation. We leave detailed analysis of the general case to a future study.
Iv conclusion
In this work, we investigate the simple onedimensional driven model– the totally asymmetric exclusion process, coupled to a modified form of Langmuir kinetics. This model is motivated by recent studies on clustering of motor proteins on microtubules. Without addressing the underlying mechanism, it is assumed that the attachment and detachment rates of a particle depend on the occupancy of the nearest neighbours of any given site. Ignoring density correlations, we obtain continuum meanfield equation describing the density profile on the lattice. Imposing certain conditions, we obtain analytical solution to the equation and demonstrate using Monte Carlo simulations that our analytical solutions are accurate over a wide range of parameters. We show that when both attachment and detachment rates are enhanced due to mutual interactions, bulkeffects start dominating the phase behavior of the model. The twophase coexistence (low and high density) observed in absence of mutual interactions can become threephase coexistence (low and high density with maximum current phase) when mutual interactions are attractive. On varying the mutual interaction between particles (attractive or repulsive), we obtain a very rich phase diagram describing the behavior of driven diffusive system with mutually interactive Langmuir kinetics.
References
 Krug (1991) J. Krug, Phys. Rev. Lett. 67, 1882 (1991).
 Evans et al. (1995) M. Evans, D. Foster, C. Godreche, and D. Mukamel, Physical review letters 74, 208 (1995).
 Evans et al. (1998) M. Evans, Y. Kafri, H. Koduvely, and D. Mukamel, Physical review letters 80, 425 (1998).
 Kafri et al. (2002) Y. Kafri, E. Levine, D. Mukamel, G. Schütz, and J. Török, Physical review letters 89, 035702 (2002).
 Parmeggiani et al. (2003) A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. Lett. 90, 086601 (2003).
 Parmeggiani et al. (2004) A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. E 70, 046101 (2004).
 Spohn (1991) H. Spohn, Large scale dynamics of interacting particles, vol. 825 (Springer, 1991).
 Privman (2005) V. Privman, Nonequilibrium statistical mechanics in one dimension (Cambridge University Press, 2005).
 Cates and Evans (2010) M. E. Cates and M. R. Evans, Soft and Fragile Matter: Nonequilibrium Dynamics, Metastability and Flow (PBK) (CRC Press, 2010).
 Domb (2000) C. Domb, Phase transitions and critical phenomena, vol. 19 (Academic Press, 2000).
 MacDonald et al. (1968) C. MacDonald, J. Gibbs, and A. Pipkin, Biopolymers 6, 1 (1968).
 Howard et al. (2001) J. Howard et al., Mechanics of motor proteins and the cytoskeleton (Sinauer Associates Sunderland, MA, 2001).
 Neri et al. (2013) I. Neri, N. Kern, and A. Parmeggiani, New Journal of Physics 15, 085005 (2013).
 Seitz and Surrey (2006) A. Seitz and T. Surrey, The EMBO journal 25, 267 (2006).
 Muto et al. (2005) E. Muto, H. Sakai, and K. Kaseda, The Journal of cell biology 168, 691 (2005).
 Crevel et al. (2004) I. M.T. Crevel, M. Nyitrai, M. C. Alonso, S. Weiss, M. A. Geeves, and R. A. Cross, The EMBO journal 23, 23 (2004).
 Telley et al. (2009) I. A. Telley, P. Bieling, and T. Surrey, Biophysical journal 96, 3341 (2009).
 Roos et al. (2008) W. H. Roos, O. Campàs, F. Montel, G. Woehlke, J. P. Spatz, P. Bassereau, and G. Cappello, Physical biology 5, 046004 (2008).
 Katz et al. (1987) S. Katz, J. Lebowitz, and H. Spohn, J. Stat. Phys 34, 497 (1987).
 Popkov et al. (2003) V. Popkov, A. Rákos, R. D. Willmann, A. B. Kolomeisky, and G. M. Schütz, Physical Review E 67, 066117 (2003).
 Krebs et al. (2004) A. Krebs, K. N. Goldie, and A. Hoenger, Journal of molecular biology 335, 139 (2004).