Theoretical investigations of asymmetric simple exclusion processes for interacting oligomers

There are many natural and technological processes that can be viewed as non-equilibrium systems where multiple participating particles are driven by constant supply and dissipation of energy [1, 2]. These processes include RNA translation by ribosomes [3, 4], vesicles locomotion [5], intracellular transport [2], transport in ion channels [6], vehicular traffic in highways [7, 8], and many others. One of the most known examples of such systems is motor proteins, also called biological molecular motors [9]. They utilize the free energy of chemical hydrolysis of ATP molecules (or related chemical reactions) to move processively along cytoskeleton protein filaments, such as actin filaments and microtubules, while carrying various cellular cargoes [9]. Although significant number of experimental and theoretical investigations of multi-particle non-equilibrium systems have been performed, in many cases the fundamental aspects of underlying molecular mechanisms remain not well understood.

Various out-of-equilibrium dynamical processes in chemistry, physics and biology have been successfully analyzed using one-dimensional stochastic models called asymmetric simple exclusion processes [10–13]. The simplest version of these models, known as a totally asymmetric simple exclusion process (TASEP), consists of one-dimensional lattice in which particles move unidirectionally following the exclusion principle, i.e. hopping one lattice site forward each time if the next site is available. In this model, the particles also occupy no more than one site, implying that the size of the particles is equal to their forward step size on the lattice. On a ring, TASEP provides useful insights on the mechanisms that lead to the fundamental current-density relations [14]. In an open system, particles are injected into the lattice with a rate α and are ejected out with a rate β, which leads to three distinct stationary phases, such as a low density (LD), a high density (HD) and a maximal current (MC) phase [15, 16]. The simplest TASEP model has been analyzed by many theoretical methods including matrix product ansatz [17], recursion methods [18], Bethe ansatz, and domain wall theory [19]. Moreover, the detailed simple mean-field arguments [18] justify the exact results by pointing out to the absence of correlations in the simplest TASEP model.

Experimental investigations suggest that in many biological transport processes, such as RNA translation and intracellular transport by molecular motors, the participating molecules are larger than their step sizes, and they typically interact with each other beyond the simple exclusion [2, 9, 20–24]. There are several theoretical studies that incorporate the nearest-neighbor particle interactions in extending the original TASEP model with a goal of application of the results to the transport of molecular motors [25–35]. The main idea of these TASEP models for interacting particles is to modify the hopping rates depending on the state of the nearest or next-nearest sites for the given particle. The interactions lead to the appearance of correlations in these systems, which can be analyzed using a time-density functional approach [29–31] and different cluster mean-field theories [34, 35]. In particular, in [34, 35], the interactions were accounted by utilizing a thermodynamically consistent approach, which argued that the process of particles coming together or breaking away from each other can be viewed as effective chemical reactions. It has been suggested that this method might provide a more realistic description of the transport by motor proteins such as kinesins [34, 35]. However, these studies were done only for particles of size 1, i.e. when the step size is the same as their molecular diameter. Similarly, the time-dependent density functional approach was also implemented only for the monomers [29–31].

Furthermore, the TASEP model with extended particles in the absence of interactions has also been well explored [36, 37]. These extended particles might be called oligomers because they occupy several lattice sites, in contrast to monomers in the standard TASEP, which cover only one site. It was found that for open boundary conditions (OBC), the TASEP model for non-interacting oligomers also has three stationary phases in the phase diagram, similar to the simple TASEP. However, the phase boundaries depend on the size of the particles. A simple mean-field description could not explain the dynamic properties of TASEP with extended particle because of the appearance of correlations due to large sizes of the oligomers. These systems were solved analytically using advanced mean-field approaches and an extremal principle arguments based on the domain wall theory [36, 37]. The obtained theoretical predictions results were validated with Monte Carlo computer simulations [36, 37].

Recently, Narasimhan and Baugaetner theoretically investigated the dynamics of interacting particles of arbitrary sizes [38]. They presented a mean-field description of the process based on a so-called discrete Takahashi formalism, which was also tested with computer simulations. Although this elegant theoretical approach was able to describe dynamic properties of the system in some cases (very strong repulsion and very weak attractions), it failed to describe the system in large region of the parameter space, including not very strong repulsions and strong attractions, which also seems to be the most relevant for real biological processes. Also the method was not able to capture a non-monotonic dependence of the MC current as a function of the interaction strength. In addition, particles of size 1–3 were exclusively considered. Furthermore, only the symmetric splitting of interactions with respect to the particle transition rates has been analyzed.

In this paper, we develop a comprehensive theoretical method of analyzing the non-equilibrium dynamics of interacting oligomers moving along the one-dimensional lattice. Our method is based on the cluster mean field approach which explicitly takes into account all dynamic processes inside the cluster and neglects the correlations between different clusters. Specifically, the two-site cluster method is employed in evaluating stationary state properties of the system, and arbitrary sizes of oligomers are considered. It is found that theoretical predictions agree well with extensive Monte Carlo computer simulations for most ranges of the parameters.

2.1. Model

To analyze the one-dimensional dynamics of interacting extended particles, such as the motion of ribosomes of diameter (size)  ∼20 nm on RNA with a step size  ∼1 nm (3 nucleotides), we consider a model presented in figure 1. The RNA molecule is viewed as the one-dimensional lattice with L ($L \gg 1$) sites. The ribosome on the RNA is modeled as a particle of size $l \geqslant 1$ (l-mer or oligomer), which covers l consecutive lattice sites: see figure 1. Each lattice site can be empty or covered by a particle, and we assign occupation variables $\tau_i$ to each site i $(1\leqslant i \leqslant L)$ such that $\tau_i = 0$ when the site is empty and $\tau_i = 1_k$ ($k=1, 2, \cdots, l$) if it is covered by the kth segment of the oligomer, counting from the left end of the extended particle, as shown in figure 1(a). For convenience, the location of the oligomers will be specified by the position of the leftmost segment of the extended particle.

Zoom In Zoom Out Reset image size

In addition to the exclusion, two neighboring l-mers (located at the sites i and i  +  l, respectively) interact with each other with an energy E (in units of $k_{\rm B}T$) where E  >  0 or E  <  0 represents attractive or repulsive interactions, correspondingly. The l-mer at the site i moves one lattice site to the right provided that the site $(i+l)$ is empty. The hopping rates are specified by the interactions and they depend on the occupancy state of $(i-1)$th and $(i+l+1)$th sites. In the case when the $(i-1)$th site is empty, the hopping rate is $q = e^{\theta E}$ in the presence of an oligomer at the $(i+l+1)$th site; otherwise the rate is 1 (see figure 1(b)(i–ii)). If the site $(i-1)$th is covered, the hopping rate is $r = e^{(\theta - 1)E}$ and 1, respectively for the empty and occupied state of the site $(i+l+1)$ (see figure 1(b)(iii–iv)). Here, the rates q and r, respectively describe the formation and breaking the particle-particle bond. The process of bond breaking and creating is viewed as a reversible chemical reaction, and thus the rates q and r follow the relation $\frac{q}{r} = e^E$, which describes an effective chemical equilibrium for the particles association and disassociation. The parameter θ $(0\leqslant \theta\leqslant 1)$ specifies how the energy E is distributed between the rates q and r.

To understand the fundamental current-density relation in the dynamics of interacting oligomers, we employ the periodic boundary conditions (PBC) in which the total number of extended particles (say M) are conserved. This yields the extended particle density in the steady state as $\rho = M/L$. Since each extended particle of size l covers consecutive l sites, the steady-state coverage density is given as $\rho_c = Ml/L$ and the density of empty sites (holes) is thus identified as $\rho_h = 1-\rho_c=1-l\rho$. In our analysis, we also consider the thermodynamic consistent OBC to obtain the stationary-state dynamic behavior. In the OBC, an l-mer can enter the lattice from the left only if all the first l sites are empty. Furthermore, the entrance rate is α if $(l+1)$th site is empty; otherwise the rate is $q\alpha$ because the entering oligomer will make a bond with already present particle at the $(l+1)$th site. The oligomer at the end of the lattice sitting on the site L  −  l  +  1 exits completely with a rate β if there is no particle on the site L  −  l; otherwise the rate of leaving is $r \beta$ which reflects the breaking of the bond between two oligomers.

2.2. Cluster mean-field theory

For complex dynamic processes that involve correlations simple mean-field methods, which assume that the occupation of neighboring sites are independent from each other, fail to properly describe the systems properties. In this case, more advanced cluster mean-field methods are required [39]. The main idea here is to consider dynamics inside of a cluster of several sites exactly, while the correlations between the clusters can be neglected. This is the approach that we utilize for analyzing TASEP of interacting oligomers.

To explain this method in more detail, we employ a two-site cluster model in which the clusters consisting of n  =  2 consecutive lattice sites are considered. If one takes a sequence of m sites in the configuration $(\tau_{i}, \tau_{i+1}, \cdots, \tau_{i+m-1})$, then the probability of such configuration is given by $P(\tau_{i}, \tau_{i+1}, \cdots, \tau_{i+m-1})$. In the two-cluster method, the probability of the sequence of m−  sites is factorized to the product of two-site cluster probabilities [40] i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(\tau_{i},\tau_{i+1},\cdots,\tau_{i+m-1}) \propto P(\tau_{i},\tau_{i+1})P(\tau_{i+1},\tau_{i+2})\cdots P(\tau_{i+m-2},\tau_{i+m-1}), \label{eq1} \nonumber \end{align} \tag{ 1 }$$

which, after a proper normalization, becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(\tau_{i},\tau_{i+1},\cdots,\tau_{i+m-1}) = \frac{P(\tau_{i},\tau_{i+1})P(\tau_{i+1},\tau_{i+2})\cdots P(\tau_{i+m-2},\tau_{i+m-1})}{P(\tau_{i+1})P(\tau_{i+2})\cdots P(\tau_{i+m-2})}.\label{eq2} \nonumber \end{align} \tag{ 2 }$$

This expression can be understood in the following way. Since the probabilities of the bulk sites $(i+1, i+2, \cdots, i+m-2)$ of the sequence of m sites are accounted twice equation (1), it must be corrected with the division of the probability of the corresponding bulk sites inside the sequence of m sites as given by equation (2).

In the system with oligomers of size l, a cluster of two consecutive sites, $(\tau_{i-1}, \tau_{i})$, can exist in one of $(l+3)$ possible states. This is illustrated in figure 2. There are $(l-1)$ states when both sites are parts of the same extended particle because for the oligomer of size l there exists $(l-1)$ boundaries connecting the segments of the same oligomer. These states can be described as $\{1_{k}, 1_{k+1}\}$ for $k=1, 2, \cdots, l-1$: see figure 2. The other four states are: $\{0, 0\}$ when both sites are empty; {1_l,0} when the right end of the oligomer is followed by the empty site; $\{1_{l}, 1_{1}\}$ when the right end of one particle touches the left end of the second particle; and {0,1_l} when the empty site is before the left end of the oligomer (figure 2). Then the normalization requires that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(0,1_1) + P(0,0) + P(1_l,0) + P(1_1,1_2) + P(1_2,1_3) + \cdots + P(1_l,1_1) = 1.\label{eq3} \nonumber \end{align} \tag{ 3 }$$

Zoom In Zoom Out Reset image size

The average occupancy of any site is given as $P(1_k) = \rho$, where ρ is the density of the oligomers in the system. Since the particle covers l consecutive sites simultaneously, the hole density is given as $P(0) = 1 - l\rho$. Using the definition of density, symmetry arguments and Kolmogorov consistency conditions the following relations between probabilities of the two-cluster sites can be obtained:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(0,1_1)+ P(1_l,1_1) = \rho,\label{eq4} \nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(1_1,1_2) = P(1_2,1_3)= \cdots P(1_{l-1},1_l)=\rho,\label{eq5} \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(1_l,0) + P(0,0) = 1- l\rho,\label{eq6} \nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P(1_l,0) = P(0,1_1).\label{eq7} \nonumber \end{align} \tag{ 7 }$$

To simplify calculations, we denote the three unknown probabilities $P(1_l, 1_1)$, P(0,1₁) and $P(0, 0)$ by x, y and z, respectively, which reduces the above equations to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x + y = \rho, \label{eq8} \nonumber \end{align} \tag{ 8 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y+z = (1-l\rho). \label{eq9} \nonumber \end{align} \tag{ 9 }$$

Furthermore, the ratio $\frac{q}{r}$ can be viewed as an effective equilibrium constant for the process of creating and breaking the bonds between two neighboring oligomers, yielding

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{q}{r} = \frac{P(1_1,1_2,\cdots, {1_l},1_1,1_2, \cdots, 1_l)P(0,0)}{P(1_1,1_2,\cdots ,1_l,0)P(1_1,1_2,\cdots, 1_l,0)}, \nonumber \end{align} \tag{ 10 }$$

which using the correlation relation equation (2) reduces to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y^2 = \frac{rxz}{q}.\label{eq16} \nonumber \end{align} \tag{ 11 }$$

Equations (8), (9) and (11) can be easily solved,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle x = \left\{\begin{array}{@{}l@{}}\rho + \frac{1}{2(q-r)}\Big[ r(1 - (l-1)\rho)-\sqrt{r^2(1- (l-1)\rho)^2+4(q-r)r \rho(1-l\rho)} \Big]; \nonumber \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q,r \neq 1 \nonumber \\ \nonumber \\ \frac{\rho^2}{1-(l-1)\rho};~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q,r = 1, \end{array} \right.\label{eq5a} \nonumber \end{align} \tag{ 12 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y = \left\{\begin{array}{@{}l@{}} -\frac{1}{2(q-r)}\Big[ r(1 - (l-1)\rho)-\sqrt{r^2(1- (l-1)\rho)^2+4(q-r)r \rho(1-l\rho)} \Big]; \nonumber \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q,r \neq 1 \nonumber \\ \nonumber \\ \frac{\rho(1-l\rho)}{(1-(l-1)\rho)};~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q,r = 1, \end{array} \right.\label{eq5b} \nonumber \end{align} \tag{ 13 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle z = \left\{\begin{array}{@{}l@{}} (1-l\rho)+ \frac{1}{2(q-r)}\Big[ r(1 - (l-1)\rho)-\sqrt{r^2(1- (l-1)\rho)^2+4(q-r)r \rho(1-l\rho)} \Big]; \nonumber \\~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q,r \neq 1 \nonumber \\ \nonumber \\ \frac{(1-l \rho)^2}{1-(l-1)\rho};~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{\it if}~ q, r = 1. \end{array} \right. \label{eq5c} \nonumber \end{align} \tag{ 14 }$$

The above expressions for x, y and z provide a full description of the stationary-state probabilities for $(l+3)$ two-site cluster states in terms of the particle density and transitions rates q and r. For the simple case of monomers, i.e. for l  =  1, the above expressions reduce, as expected, to the respective probabilities for the TASEP model of interacting particles of size 1 [41]. For strong repulsions between l-mers ($E \rightarrow -\infty$), these equations simplify into $y=\rho$, x  =  0 and $z=1-(l+1)\rho$, which is identical to corresponding expressions for the case of non-interacting $(l+1)$-mers. In the opposite case of strong attractions ($E \rightarrow +\infty$) we have y  =  0, $x=\rho$ and $z=1-l\rho$.

It is important to mention that by specifying every hole and every oligomer with the size l on a lattice with L sites, respectively, as a hole and a new effective monomer on another lattice with $L'$ sites, we obtain a mapping $L'/L = 1-(l-1)\rho_l$, where $L' = L-(l-1)M$ and $\rho_l$ denotes the density of the system with M l-mers and L sites. In the light of the mapping, $\rho_l = (L'/L)\rho_1$, where $\rho_1$ denotes the density of the system with M unit-size particles and $(L-lM)$ holes. Using this mapping, the two-cluster probabilities x, y and z for l-mers given by equations (12)–(14) can be alternatively obtained from relations $x_l = (L'/L)x_1$, $y_l = (L'/L)y_1$ and $z_l = (L'/L)z_1$, respectively. Note that the mapping works only for the PBC with the fixed density of particles.

Now let us investigate the fundamental current-density relation for the TASEP model with interacting oligomers. There are four different contributions to the particle flux through the system as indicated in figure 1(b). The bulk current arising from each of the configurations (i)–(iv) in figure 1(b) can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm i} = P(0,1_1,1_2,\cdots,1_l,0,0), \nonumber \end{align} \tag{ 15 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm ii} = q P(0,1_1,1_2,\cdots,1_l,0,1_1), \nonumber \end{align} \tag{ 16 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm iii} = r P(1_l,1_1,1_2,\cdots,1_l,0,0), \nonumber \end{align} \tag{ 17 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm iv} = P(1_l,1_1,1_2,\cdots,1_l,0,1_1). \nonumber \end{align} \tag{ 18 }$$

Using the two-cluster mean field theory (equation (2)) and equations (4)–(9), the total particle current $J_{\rm bulk} = J_{\rm i} + J_{\rm ii} + J_{\rm iii} + J_{\rm iv}$ can be written as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm bulk}= \frac{y^2(x+z) + y(qy^2 + rxz)}{\rho (1- l \rho)}. \nonumber \end{align} \tag{ 19 }$$

After substituting the values of x and z from equations (12) and (14), we obtain the following fundamental current-density relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm bulk} = ry + \frac{(1-r)(1-(l-1)\rho)}{\rho (1-l\rho)}y^2 + \frac{(q+r-2)}{\rho (1-l \rho)}y^3,\label{eq22} \nonumber \end{align} \tag{ 20 }$$

where y is a function of ρ as given by equation (13). For l  =  1, the above current-density relation reduces to the relation obtained in [41] using the same approach. For the case of no interactions (q  =  r  =  1), $J_{\rm bulk} = y = \frac{\rho(1-l\rho)}{(1- (l-1)\rho)}$, which matches with the current-density relation for the non-interacting l-mers [36]. For strong repulsive interactions, $E \rightarrow - \infty$, any two consecutive l-mers are always separated by at least a single hole because it is energetically unfavorable to make a bond between two neighboring particles. As a result, the interacting l- mers behave as non-interacting $(l+1)$-mers at this limit. This is justified as for $E \rightarrow - \infty$ we obtain $x \rightarrow 0$, $y \rightarrow \rho$ and $z \rightarrow 1-(l+1)\rho$, and under these conditions the bulk current for the l-mers approaches to $J_{\rm bulk} = \frac{\rho(1-(l+1)\rho)}{(1-l\rho)}$, which is exactly the same as the bulk current for the non-interacting $(l+1)$-mers [42]. For strong attractions, $E \rightarrow +\infty$, the current goes to zero, $J_{\rm bulk} \rightarrow 0$, for any value of l. This is easy to understand since in this case the oligomers group together to form large jamming clusters of particles which hinder their movement.

The current-density relation obtained in equation (20) depends on the particle size l as well as on the interaction energy E. To investigate the variation of the current-density relation with respect to both energy of interactions and the size of the oligomers, in figure 3(a) we present a 3D plot of these functions for the symmetric splitting of the interaction on transition rates ($\theta=0.5$). One can see that, similarly to the case of interacting monomers [41], the fundamental diagram changes its behavior when the interaction strength is varied. For attractive interactions and weak repulsions the curve has only one maximum, while for stronger repulsions there are two maxima and one minimum in the current-density curves. It is also clear from figure 3(b) where the fundamental diagrams are presented for E  =  −5 $k_{\rm B}T$ and for various sizes of the oligomers.

Zoom In Zoom Out Reset image size

The existence of the single MC for some intermediate density is easy to understand. The particle flux, $J \sim \rho v$, is maximal for the density that does not decrease significantly the velocity of each particle. However, the current-density curve with several extremal points is more unexpected. To explain this observation, let us consider the case of strong repulsions ($E \rightarrow -\infty$) when densities and fluxes at the extremal points can be obtained analytically. At the minimum position the flux goes to zero ($J_{\rm min} \rightarrow 0$) because for strong repulsions the system of interacting l-mers is identical to the system of non-interacting $(l+1)$-mers. This argument suggest that the density at the minimal point is $\rho_{\rm min} \rightarrow\frac{1}{(l+1)}$. The whole lattice will be covered by the particles of size l  +  1, which should lead to zero current. The density corresponding to the first maximum in the current-density curve can be obtained by solving the equation $\frac{\partial J_{\rm bulk}}{\partial \rho} = 0$. This leads to to

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle \rho_1 = \frac{1}{l}\bigg(1 - \frac{1}{\sqrt{(1+l)}}\bigg)=\frac{1}{(l+1)+\sqrt{l+1}}. \nonumber \end{align} \tag{ 21 }$$

The second maximum can be obtained by noting that for densities larger than $1/(l+1)$ it is more convenient to look at the particle flux in one direction as the flux of interacting 'holes' in the opposite direction. Because of the strong repulsion, the dynamics in the system at these conditions can be viewed as the motion of the non-interacting 'holes' of size l  +  1, for which the current is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{h} = \frac{\rho_{h}(1-(l+1)\rho_{h})}{(1-l\rho_{h})}. \nonumber \end{align} \tag{ 22 }$$

This flux has the maximum for $\rho_{h}^{*}=1/(l+1+\sqrt{l+1})$, which gives the location of the second maximum in the fundamental current-density relation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{2}=1-l\rho_{h}^{*}=\frac{l+\sqrt{l+1}}{l(l+1+\sqrt{l+1})}. \nonumber \end{align} \tag{ 23 }$$

Of the two humps in the curve, the particle current is larger at the density $\rho_{\rm max} \rightarrow \rho_1$, and it is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm max} \simeq \frac{1}{(\sqrt{l+1}+1)^{2}}. \nonumber \end{align} \tag{ 24 }$$

The density and the current at the first maximal point can be also obtained by considering the non-interacting oligomers of size l  +  1 as explained above. The reason for the second maximum and smaller particle flux at this point can be understood as follows. As we already noticed above, the motion of interacting oligomers of size l moving in one direction can be viewed as a motion of interacting holes of size 1 in the opposite direction. Thus, for densities larger than $1/(l+1)$, the number of holes is small and the corresponding current at the second maximum, which represents the maximum flow of holes, becomes smaller than the MC at the first maximum.

To evaluate the coordinates of the extremal points in the current-density curve for general conditions, we again solve the equation $\frac{\partial J_{\rm bulk}}{\partial \rho} = 0$. It is found that the number of extremal points depends on the values of the interaction energy E, on the splitting parameters θ and on the size of the oligomer l. By decreasing the interaction energy from strong attractions to strong repulsions, it can be shown that there is a critical energy $E_c(\theta, l)$ for $0 < \theta \leqslant 1$ below which the current-density relation always has three extrema, and above $E_c(\theta, l)$, the number of extremal points reduces to one. Table 1 displays the numerical values for critical interaction strength $E_c(\theta)$ for different sizes of oligomers and for different splitting parameters θ. We observe that as the size of the oligomers increases, the critical interaction strength becomes more repulsive for any value of θ ($0 < \theta \leqslant 1$). This can be explained by analyzing equation (20). The first term in this equation is always positive, the second term is negative for E  <  0, and the third term could be negative or positive depending upon the sign of $(q+r-2)$. The parameter y, which is the probability to have a configuration with the particle being followed by the hole, is a small number. Then in most cases the first and the second term in equation (20) dominate the behavior. It can be argued that the critical interaction approximately corresponds to the situation when the first term is roughly compensated by the second term. Since the parameter y decreases with the increase in the size of the oligomers then to compensate the transition rate r must be large at the critical condition. This corresponds to more negative interactions. The results presented in table 1 also suggest that for the fixed size of the oligomer lowering the value of the parameter θ makes the critical interaction more repulsive. For the special case of $\theta = 0$, when interactions only affect the bond-breaking rate r, with q  =  1 at all conditions, we find that equation (20) always has one extremal point corresponding to the MC for any size l (see figure 4(a)).

Zoom In Zoom Out Reset image size

It is interesting to note here that, in contrast to the case of interacting monomers (l  =  1), the density ($\rho_{\rm max}$) corresponding to the largest MC for interacting l-mers (l  >  1) is influenced by the strength of interactions and by varying the splitting parameter θ. Moreover, $\rho_{\rm max}$ also depends on the size of the oligomers. When attractions affect the bond-breaking rate more than the bond-making rate, i.e. for $0 \leqslant \theta \leqslant 0.5$, $\rho_{\rm max} \rightarrow 0$ as $E \rightarrow \infty$ (see figures 4(a) and (b)). This is justified as for very strong attractions, we have $r \ll 1$ that causes the oligomers to group together, and in this situation the particle flux can be maximum only if there are less number of particles in the system. This obviously should decrease the maximal particle current. In the case when the bond-making rate is influenced more by the interactions in comparison to the bond-breaking rate (for $\theta > 0.5$) in the limit of strong attractions the density of the MC state approaches to $\rho_{\rm max} \rightarrow 1/2l$ (see figure 4(c)). This can be explained by considering the case of $\theta=1$ using the following arguments. At these conditions we have r  =  1, and the bulk current can be rewritten as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm bulk} = y + \frac{(q-1)y^{3}}{\rho (1-l \rho)}.\label{eq_special} \nonumber \end{align} \tag{ 25 }$$

For strong attractions, $q \gg 1$ and from equation (13) the parameter y takes the following form

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle y \simeq \frac{\sqrt{\rho(1-l\rho)}}{\sqrt{q}}. \nonumber \end{align} \tag{ 26 }$$

Then the particle flux (equation (25) is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J \simeq \frac{2\sqrt{\rho(1-l\rho)}}{\sqrt{q}}, \nonumber \end{align} \tag{ 27 }$$

which reaches a maximum at $\rho=1/2l$.

Now let us consider a system of interacting oligomers under OBC. The particle flux in the bulk of the system is given by equation (20). The entrance current for the oligomers of size l can be expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm entr} = \alpha P\underbrace{(0,0,\cdots,0,0)}_{{(l+1)}} + q\alpha P\underbrace{(0,0,\cdots,0,1)}_{{(l+1)}}, \nonumber \end{align} \tag{ 28 }$$

which in the two-site cluster mean-field analysis gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{j-entr} J_{\rm entr} = \alpha z + q\alpha y. \nonumber \end{align} \tag{ 29 }$$

For a special case of no interactions (E  =  0), we have $J_{\rm entr} = \alpha(1-l\rho)$, which is expected for the case of non-interacting oligomers of size l. When the repulsion is very strong ($E \rightarrow -\infty$), $J_{\rm entr} \rightarrow \alpha (1-(l+1)\rho)$, which agrees with the known results of the TASEP with non-interacting $(l+1)$-mers [36, 37]. For strong attractions ($E \rightarrow \infty$), it can be shown that $J_{\rm entr} \rightarrow 0$, as expected. Similar arguments for the current of exiting particles lead to the following expression,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{j-exit} J_{\rm exit} = \frac{(\beta y + r \beta x)}{(1-(l-1)\rho)}. \nonumber \end{align} \tag{ 30 }$$

When there are no interactions between particles (E  =  0) and for strong repulsions ($E \rightarrow -\infty$), we obtain $J_{\rm exit}= \beta \rho/(1-(l-1)\rho)$, in agreement with the previous analysis for TASEP with the non-interacting oligomers [36] and for the interacting monomers (l  =  1) [35].

4.1. Phase diagrams

Now we investigate the effect of interactions on the stationary phase diagrams for oligomers of arbitrary size. It is to be noted that the thermodynamically consistent boundary conditions in our model for interacting oligomers do not allow us to apply the maximum and minimum current principles to predict all possible stationary phases [41]. To evaluate what dynamic phases are realized, we apply the following method. At the phase boundaries the particle fluxes corresponding to different phases (entrance, bulk or exit), which can be obtained from our two-site cluster calculations as explained above, must be equal. Thus we identify three possible stationary phases in the system of interacting oligomers, namely, LD, HD and MC. This is very similar to all other TASEP models, and interactions are influencing the positions of the phase boundaries. Let us discuss the phase diagram in detail.

4.1.1. MC phase.

In the MC phase, the particle current is maximal and it is independent of dynamics at the system's boundaries. We differentiate J_bulk in equation (20) with respect to the particle density ρ and equate it to zero in order to obtain the stationary density profile $\rho_{\rm MC}$ in the MC phase. Then the specific expressions for the particle current is given by $J_{\rm bulk}(\rho=\rho_{\rm MC})$. For E  =  0 (no interactions) it can be easily shown that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{\rm MC}=\frac{1}{l+\sqrt{l}}, \quad J_{\rm MC}=\frac{1}{(\sqrt{l}+1)^{2}} \nonumber \end{align} \tag{ 31 }$$

while for strong repulsions we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{\rm MC}=\frac{1}{(l+1)+\sqrt{l+1}}, \quad J_{\rm MC}=\frac{1}{(\sqrt{l+1}+1)^{2}}. \nonumber \end{align} \tag{ 32 }$$

Again, one can clearly see that the dynamics of interacting l-mers at strong repulsions is identical to the dynamics of non-interacting $(l+1)$-mers. For strong attractive interactions, the particle density depends on the value of the splitting parameter θ. As we argued above, $\rho_{\rm MC} \rightarrow 0$ for $0 \leqslant \theta \leqslant 0.5$, while $\rho_{\rm MC} \rightarrow 1/2l$ for $0.5 < \theta \leqslant 1$. But in both cases the particle fluxes are disappearing at this limit ($J_{\rm MC} \rightarrow 0$).

The results of our calculations for the particle fluxes as a function of the interaction strength for different sizes of oligomers are presented in figure 5 where they are also compared with Monte Carlo computer simulations. Excellent agreement is observed for all values of the splitting parameter θ except for the limiting cases of $\theta=0$ and $\theta=1$. It is found that the particle current decreases with the size of oligomers. There is a maximal particle flux observed at weak repulsions, while at strong repulsions the current approaches a constant value, and the current disappears for strong attractions. It is also found that the position of the most optimal current does not depend on the size of oligomers. This picture is generally valid for oligomers of any size and for $0 < \theta < 1$.

Zoom In Zoom Out Reset image size

For $\theta=0$, computer simulations predict a monotonic decrease in the particle current as a function of the interaction strength. While our theoretical calculations qualitatively agree with this behavior, the particles fluxes for repulsive interactions are overestimated. These deviations between the theory and computer simulations can be explained using the following arguments. In this regime, the formation of bonds between oligomers are not affected by the interactions (q  =  1), but the bond breaking is affected much stronger (r  =  e^−E), and the effect is especially strong for repulsive interactions. For $\theta=1$, computer simulations suggest that the particle flux is an increasing function of the interaction strength. Our theory correctly describes the behavior for the repulsive interactions, but incorrectly predicts the maximum in the particle current and it fails for the attractions. In this regime, the formation of bonds between oligomers strongly depends on interaction energy (q  =  e^E), and is strongest for large attractions, while the breaking of bonds between the oligomers is not affected by E at all (r  =  1). These arguments suggest that our theoretical method cannot capture correctly all the aspects of the dynamics at these limiting regimes, while it describes perfectly the large fraction of the parameters space outside of the limiting cases.

4.1.2. LD phase.

In the LD phase, the dynamics are governed by processes of particles entering into the system. The relation between the oligomers density and the entrance rate can be obtained by equating the expressions for J_bulk and J_entr, yielding

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha = &~\Bigg[\frac{A+B}{\Big(2(q-r)^2\rho (l\rho - 1)\Big)\Big(r((l+1)\rho -1) - B + q(2(1-l\rho) + A+B)\Big)}\Bigg]\nonumber \\ & \times \Big(4q^2r\rho(l\rho - 1) + r^2(1-\rho(l+1))^2 + AB -2qr^2((1-l\rho)^2 +\rho^2)+q(1+(1-l)\rho)B\nonumber \\& + qr(1+(1-6l+l^2)\rho^2 + 2B - 2\rho(l-3+(l-1)B))\Big),\label{eqac} \nonumber \end{align} \tag{ 33 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle A = r((l-1)\rho -1), \nonumber \end{align} \tag{ 34 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle B = \sqrt{4qr\rho(1-l\rho) + r^2(1-\rho(l+1)^2)}. \nonumber \end{align} \tag{ 35 }$$

The equation (33) can be solved to obtain the density in the LD phase, $\rho_{\rm LD}(\alpha)$, for any set of values for l, θ and E. For the special case of strong repulsions $(E \rightarrow - \infty)$, the above equation yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{\rm LD} = \frac{\alpha}{(1 + l\alpha)}, \nonumber \end{align} \tag{ 36 }$$

which is the expected outcome because it corresponds to the case of non-interacting oligomers of size $(l+1)$ [36]. For strong attractions, we obtain $\alpha \rightarrow 0$, which implies that the LD phase does not exist at these conditions for any size of the oligomers. For the case of no interaction (E  =  0), the condition of $J_{\rm bulk} = J_{\rm entr}$ yields $\alpha = \frac{\rho}{1-\rho(l-1)}$ and then $\rho_{\rm LD} = \frac{\alpha}{1+\alpha(l-1)}$, as expected for the case of non-interacting oligomers of size l [36]. The particle current in the LD phase can be evaluated from equation (29) with $\rho=\rho_{\rm LD}$. For example, for strong repulsions it gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm LD}=\frac{\alpha(1-\alpha)}{1+l\alpha}. \nonumber \end{align} \tag{ 37 }$$

4.1.3. HD phase.

The HD phase is defined by the particles exiting processes. Equating J_bulk and J_exit, we obtain the following relation:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \beta =&~ \Bigg[\frac{A(A+B)}{r(2(q-r)^2\rho (l\rho - 1))(r^2((l+1)\rho -1) - B + r(1+\rho(1-l-2q)+B))}\Bigg]\nonumber \\& \times \Big(4q^2r\rho(l\rho - 1) + r^2(1-\rho(l+1))^2 + AB -2qr^2((1-l\rho)^2 +\rho^2)+q(1+(1-l)\rho)B\nonumber \\& + qr(1+(1-6l+l^2)\rho^2 + 2B - 2\rho(l-3+(l-1)B))\Big).\label{eqbc} \nonumber \end{align} \tag{ 38 }$$

This equation can be solved to evaluate the density in the HD phase as a function of the exit rate β for any values of the parameters l, θ and E. When oligomers strongly repel each other ($E \rightarrow -\infty$), it can be shown that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \rho_{\rm HD}=\frac{l(2-\beta)-\sqrt{\beta^{2}l^{2}+4(1-\beta)}}{2(l^{2}-1)}. \nonumber \end{align} \tag{ 39 }$$

For monomers (l  =  1), this reduces to $\rho_{\rm HD} = (1-\beta)/(2-\beta)$, which has been found before for TASEP of interacting monomers [35]. For the case of strong attractions ($E \rightarrow \infty$), as expected, the particle current approaches to zero and the HD region spans the complete phase diagram. In the absence of interactions between the neighboring oligomers, the condition for the HD phase yields $\rho_{\rm HD} = (1-\beta)/l$, which is fully consistent with known results for the non-interacting particles of size l [37]. Finally, the particle current in the HD phase is found by using equation (30) with $\rho=\rho_{\rm HD}$. For example, in the case of strong repulsions, we have

$$\begin{align} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \displaystyle J_{\rm HD}= \bigg(\frac{\beta}{l-1} \bigg) \Bigg(\frac{l(2-\beta)-\sqrt{\beta^{2}l^{2}+4(1-\beta)}}{2+l\beta +\sqrt{\beta^{2}l^{2}+4(1-\beta)}} \Bigg). \nonumber \end{align} \tag{ 40 }$$

We, now, estimate the positions of boundaries between different stationary phases. It can be shown that for entrance rates $\alpha < \alpha_{c}$, the system is in the LD phase with the particle density $\rho_{\rm LD} (\alpha)$, while for $\alpha > \alpha_{c}$, the system crosses into the MC phase with the particle density $\rho_{\rm MC}$, which is independent of the entrance rate. Thus, the line $\alpha = \alpha_c$ provides a second-order continuous phase transition line between the LD and the MC phases. The value of $\alpha_{c}$ can be obtained by substituting $\rho_{\rm MC}$ in equation (33). For the case of strong repulsions ($E \rightarrow -\infty$), equation (33) yields

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{alpha_c} \alpha_c = \frac{\sqrt{l+1}-1}{l}=\frac{1}{\sqrt{l+1}+1}, \nonumber \end{align} \tag{ 41 }$$

which for l  =  1, gives $\alpha_c = \sqrt{2}-1$, in agreement with the results for non-interacting dimers [36]. In the case of no interactions (E  =  0), our calculations produce $\alpha_c = 1/(1+\sqrt{l})$, which again reproduces the earlier results on TASEP for non-interacting particles [36].

Similar calculations can be done to show that the line $\beta = \beta_c$ gives the second-order continuous phase transition line between the LD and the MC phases. For strong repulsions, this leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{beta_c} \beta_c = \frac{(l^2 - l -4) + (4-l)\sqrt{l+1}}{l^2(\sqrt{l+1}-1)}=\frac{\sqrt{l+1}+2}{(\sqrt{l+1}+1)^{2}}. \nonumber \end{align} \tag{ 42 }$$

For E  =  0, we have $\beta_c = 1/(1+\sqrt{l})$, as found earlier for the non-interacting l-mers—see [36]. Furthermore, analyzing equations (41) and (42), one might conclude that increasing the size of oligomers should decrease the fraction of the phase diagram occupied by LD and HD phases in favor of the MC phase because for $l \gg 1$, we have $\alpha_{c}=\beta_{c} \simeq 1/\sqrt{l}\rightarrow 0$ (see figure 6).

Zoom In Zoom Out Reset image size

The LD–HD phase boundary line describes a first-order phase transition at which the density changes abruptly from $\rho_{\rm LD}$ to $\rho_{\rm HD}$. The phase transition line can be obtained by equating J_HD and J_LD. All three phase transition lines meet together at a special point with the coordinates $(\alpha_c, \beta_c)$, which is known as a triple point. Stationary phase diagrams for interacting oligomers of different sizes are presented in figure 6. One can see that increasing the size of the interacting particles shrinks the LD and the HD phases and increases the MC phase, as we already predicted above. This figure also shows that our theory agrees quite well with computer simulations, and the agreement is better for LD-MC phase boundary than for the HD-MC phase boundary. It is possible that the last observation is the result of weaker correlations in LD phase in comparison with HD phase because of the smaller particle density.

The extended size of the particles and the presence of interactions between them induce correlations in the system, which strongly influence the dynamics and the stationary behavior of the system. Our theoretical method allows us to explicitly evaluate the strength and the sign of such correlations. For this purpose, we introduce a correlation function C,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{correlations} C=P(\tau_{i}=1_1,\cdots,\tau_{i+l}=1_1)-P(\tau_i = 1_1)P(\tau_{i+l}= {1_1 })= x - \rho^2. \nonumber \end{align} \tag{ 43 }$$

The physical meaning of the correlation function C is the following. If the left end of the oligomer of size l is found at the site i, it affects the probability of finding the left end of another oligomer at the site i  +  l, and the correlation function C measures this effect. If there are no correlations, we have C  =  0, i.e. the probability of finding the oligomer at the neighboring site is independent of the occupancy of the given site. The case of C  <  0 means that the presence of the particle at the given site decreases the probability of finding another oligomer next to the original oligomer, i.e. the negative correlations. The case of C  >  0 describes the positive correlations when the presence of the particle at the site i increases the probability of finding the oligomer at the site i  +  l.

Figure 7 shows the results of explicit calculations for the correlation functions of oligomers of different sizes l for the different splitting parameters θ as a function of the interaction energy E. The calculations are done at the conditions that correspond to the MC phase, and we also present Monte Carlo computer simulations to compare with our theoretical predictions.

Zoom In Zoom Out Reset image size

Our theory suggests that for repulsive interactions the correlations are negative (C  <  0), they are almost independent of the value of the splitting parameter θ and the amplitude of correlations decreases with the size of the particle. One can see it more clearly when $E \rightarrow -\infty$, the equations (12) and (43) yield

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C=-\rho^{2}=-\frac{1}{(l+1+\sqrt{l+1})^{2}}. \nonumber \end{align} \tag{ 44 }$$

For attractive interactions, we predict that the correlations are positive (C  >  0), they depend on θ and generally they are stronger (although also decreasing with l). For $E \rightarrow \infty$, it can be shown that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C=\rho(1-\rho). \nonumber \end{align} \tag{ 45 }$$

The calculations indicate that the correlations disappear for strong attractive interactions when $0 \leqslant \theta \leqslant 0.5$, and C  =  (2l  −  1)/(4l²) for $0.5 < \theta \leqslant 1$. Thus, our theory predicts a non-monotonic behavior for C as a function of the interaction strength. Surprising results are found for very weak interactions between the oligomers. For E  =  0, it is found that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle C=\frac{(l-1)\rho^{3}}{(1-(l-1)\rho)}=\frac{\sqrt{l}-1}{l(\sqrt{l}+1)^2}. \nonumber \end{align} \tag{ 46 }$$

This suggests that, in contrast to interacting monomers (l  =  1), correlations are not zero for l-mers when there are no interactions: see figure 7(b) inset. The correlations for interacting oligomers disappear for weak repulsive interactions. These observations can be explained using the following arguments. In the MC phase the fraction of the empty sites decreases with the size of the particle as  ∼$1/\sqrt{l}$ even when there is no interaction. This means that the probability of finding two oligomers next to each increases with l. This corresponds to the appearance of effective interactions, which is canceled by weak repulsions. Intriguingly, the maximal possible current (see figure 5(b)) is found close to these conditions.

Comparing our theoretical predictions with the results from Monte Carlo simulations, we notice that a very good agreement is found for the repulsive interactions in all regimes, while for attractive interactions the agreement is only semi-quantitative but it strongly improves for larger l. As we explained before, the least successful performance of our theory is found at the limiting cases for $\theta=0$ (for repulsions) and for $\theta=1$ (for attractions). These observations can be explained by noting that for repulsions the density correlations are weaker than for the attractive interactions. In addition, the correlations decrease with the size of oligomers l. Because our theory is based on taking into account short correlations (inside the cluster of two neighboring lattice sites) it works well for the situations with relatively weak correlations, as found for the repulsions and generally for large l particles.

To summarize, stimulated by biological transport phenomena we developed a new theoretical description to analyze the dynamics of interacting oligomers moving on one-dimensional lattices. Our approach is based on extending the TASEP models for extended particles to include the interactions between them, which are taken into account using thermodynamically consistent arguments. Because interactions and the extended sizes of the particle induce significant correlations in the system, our method employs cluster mean-field calculations that partially account for them. Analytical calculations are also supported by extensive Monte Carlo computer simulations.

Our theoretical calculations indicate that the dynamics of interacting oligomers differ significantly from the dynamics of interacting monomers. First, the TASEP model for interacting oligomers is considered for PBC. It is found that the fundamental current-density relation depends on the strength of interactions. A single-maximum curve for attractions and for weak repulsions is modified into a two-maximum curve for stronger repulsions at the critical interaction energy. This critical interaction is more negative for larger oligomers, and it also depends on the value of the splitting parameter θ, which specifies how interactions affect the transition rates in the system. Next, our analysis shifts to the dynamics of interacting oligomers under the OBC. Similar to many TASEP models, three possible stationary phases (MC, LD and HD) are identified and fully described for thermodynamically consistent OBC. It is found that phase boundaries, particle fluxes and bulk densities depend on the size of the particle, on the strength of interactions and how these interactions affect transition rates in the system. The analysis suggests that the largest particle flux can be achieved in the MC phase for weak repulsions, and it is independent of the size of oligomers. To explain the dynamic properties of interacting oligomers, correlations in the system are analyzed in detail. We determined that correlations are negative for repulsive interactions, and positive for attractive interactions. But the amplitude of correlations decreases with the size of extended particles, and it is generally stronger for attractive interactions. It is also shown that in the case of no interactions, the correlations are positive for l-mers (although decreasing with the size), in contrast to the monomers where the correlations are zero for E  =  0. Our calculations indicate that for oligomers the correlations disappear in the case of weak repulsions, which might be related with the observation of the MC at these conditions. The presented theoretical method based on the cluster two-site mean field analysis shows a very good agreement with computer simulations for most ranges of the parameters, with the exception of some limiting cases where the qualitatively correct behavior is still predicted. We also argue that our theory works better for larger sizes of the oligomers.

Although the presented theoretical approach seems to capture the main physical features of the system of interacting extended-size particles, there are several issues that still needs to be addressed. Our theory does not work well in the limiting cases when the interaction affects only the formation of the bonds or only the breaking the bonds between the oligomers. But possibly these problems can be resolved by extending our method into three-site or more sites cluster mean-field calculations. In addition, real biological transport phenomena include dissociations and associations at any locations on the lattice, and the transition rates might be also position-dependent. It will be interesting to extend the presented theoretical method to take these realistic features into account.

ABK acknowledges support from the Welch Foundation (Grant C-1559), from the NSF (Grant CHE-1664218) and from the Center for Theoretical and Biological Physics sponsored by the NSF (Grant PHY-1427654). We also thank Dr. S L Narasimhan for sending us his preprint and for useful comments and discussions. We also thank one of the anonymous reviewers for giving an insight into the mapping approach. The research work initiated during a visit for participating in the program—Collective Dynamics of-, on- and around Filaments in Living Cells: Motors, MAPs, TIPs and Tracks (Code: ICTS/CDFLC/2017/10) held at International Center for Theoretical Sciences (ICTS), Bangalore. It was also carried out with the support of CNPq, Conselho Nacional de Desenvolvimento Científico e Tecnológico–Brasil.