format: text

text:
# The totally asymmetric exclusion process with extended objects, a model for protein synthesis

## Introduction

The process of protein synthesis, important in biological systems, has been the focus of intensive study over the last few decades. We concentrate on protein synthesis in prokaryotes, particularly Escherichia coli, because it is relatively simple and well studied. The mechanism consists of ribosomes “reading” the codons of messenger RNA (mRNA) as the ribosomes move along an mRNA chain, and the recruitment and assembly of amino acids (appropriate to the codons being read) to form a protein. (See, e.g., Stryer , for more details.) Known as “translation,” this process is often described as three steps: initiation, where ribosomes attach themselves, one at a time, at the “start” end of the mRNA; elongation, where the ribosomes move down the chain in a series of steps; and termination, where they detach at the “stop” codon. Since ribosomes cannot overlap, their dynamics is subject to the “excluded volume constraint.” Apart from being impeded by another ribosome (steric hinderance), a ribosome cannot move until the arrival of an appropriate transfer RNA, carrying the appropriate amino acid (a combination known as aminoacyl-tRNA, or aa-tRNA). Thus, the relative abundances of the approximately 60 types NU  of aa-tRNA have significant effects on the elongation rate. Assuming reactant availabilities in a cell are in their steady state, with a time-independent concentration of ribosomes and aa-tRNA, there would be an approximately steady (average) current of ribosomes moving along the mRNA, resulting in a specific production rate of this particular protein. Our goal is the prediction of the protein production rates for various mRNA’s, as a function of the concentration of ribosomes and aa-tRNA’s.

The process of translation is well suited to modeling using a driven lattice gas in one dimension. In most relevant studies of one dimensional driven lattice gases, particles are injected at some rate on one end of a chain of discrete lattice sites, then hop down the chain one site at a time with another rate, and finally exit the chain at the other end with a third rate. These three rates correspond to the rates of initiation, elongation, and termination. The excluded volume constraint is implemented by insuring that each site can accommodate at most one particle. Because the dynamics is stochastic, even this “simple” model, though soluable exactly DDM ; SD , is already not trivial. However, we believe two essential aspects of translation are missing from this model. First, if we model a codon by a lattice site, the ribosome would cover, typically, a dozen sites Heinrich ; Kang . Second, there is non-uniformity in the hopping (elongation) rates along the chain, because a ribosome has to “wait” for the appropriate aa-tRNA before continuing, and the relative abundance of the different aa-tRNA’s is far from unity. Remarkably, the first issue was explored as early as 1968 MGP ; MG , though only at the deterministic, “mean field” level.

In this paper, we present studies which address both of these issues, extending the work on the “simple model” known as TASEP (namely, single-site coverage, Totally Asymmetric Simple Exclusion Processes with open boundaries) GS-DL19 . Our methods involve both Monte Carlo simulations and modern analysis techniques, including domain wall theory KSKS . We have confirmed all key results from earlier studies MGP ; MG , and several new insights have emerged. The paper is organized as follows. The details of the model are delineated first, along with brief summaries of known results. Section II is devoted to a *closed* (i.e., periodic) system with particles of arbitrary size. Though not a direct model for translation, TASEP on a “ring” provides simple solutions as well as useful insights in the form of *exact* relations for relevant parameters, such as current-density relations. In Section III, we turn to the central topic: TASEP with extended objects and open boundaries. Section IV is devoted to non-uniform hopping (elongation) rates. We close with a brief summary and speculate on the relevance of our model as a mechanism for the nonlinear relationship between mRNA and protein levels observed in biological experiments Gygi ; Ideker .

## I The driven lattice gas as a model for protein synthesis

### I.1 Model specifications

We model an mRNA with $`N`$ codons as a chain of sites, each of which is labeled by $`i`$. The first and last sites, $`i=1,N`$, are associated with the start and stop codons respectively. At any time, attached to the mRNA are $`M`$ ribosomes (also referred to as “particles”), which we label by $`\alpha `$ ($`\alpha =1,\mathrm{},M`$). Being a large complex of molecules, each ribosome will cover $`\mathrm{}`$ sites (codons), with $`\mathrm{}=12`$ typically Heinrich ; Kang . By contrast, nearly all studies of the asymmetric simple exclusion processes (ASEP) are devoted to $`\mathrm{}=1`$. Any site may be covered by a single ribosome or none. In case of the latter, we will refer to the site as “empty” or “occupied by a hole.” For convenience, we define $`\stackrel{~}{M}`$ as the number of holes on the chain, so that

$$\stackrel{~}{M}+\mathrm{}M=N.$$
(1)
For open systems, a ribosome at the end can be attached without covering all $`\mathrm{}`$ codons, so that this equality is only approximately true. To locate the ribosome, we arbitrarily choose the *lowest* site covered. For example, if the first $`\mathrm{}`$ sites are empty, a ribosome can bind in an initiation step, and then it is said to be “on site $`i=1.`$” Therefore, a complete specification of the configuration (state, or micro-state) of the mRNA is the set of locations: $`\left\{i_\alpha \right\}`$. One disadvantage of this labeling is that, with each initiation event, the $`\alpha `$ of every ribosome will change (increase by unity). Alternatively, we can use *site occupation* numbers

$$n_i=\{\begin{array}{cc}1\hfill & \text{if site }i\text{ is covered by any part of a ribosome}\hfill \\ 0\hfill & \text{if site }i\text{ is empty}.\hfill \end{array}$$

With these conventions, we define several density parameters:

* $`\rho _rM/N`$ is the ribosome (or particle) density;
* $`\rho M\mathrm{}/N=_in_i/N`$ is the coverage density;
* $`\rho _h=1\rho `$ is the hole density; and
* $`\rho _s\rho _r+\rho _h`$ is defined for convenience.

All of these quantities are time dependent, because initiation and termination occur *independently*. For mathematical reasons, we will first consider a *closed* system (with periodic boundary conditions, i.e., the ends of the chain tied to form a ring), for which these densities are fixed and Eqn. (1) holds strictly. As will be clear later, it is also convenient to label configurations by specifying the number of holes between successive ribosomes: $`\left\{h_\alpha \right\}`$, where $`h_\alpha `$ is the number of holes *in front of* the $`\alpha ^{th}`$ ribosome. In the terminology of traffic models, $`h_\alpha `$ is also known as the “headway” of this particle. Though not absolutely necessary, we could define $`h_0`$ as the number of holes behind the first ribosome.

Next, we specify the dynamics of our model. An attached ribosome located at site $`i`$ will move to the next site ($`i+1`$) with a rate $`k_i`$, *provided* site $`i+\mathrm{}`$ is empty. For Monte Carlo simulations, it is convenient to update configurations in discrete time units. Then, it is better to use probabilities $`p_i`$ ($`0p_i1`$), so that a ribosome on site $`i`$ will be moved or not with probability $`p_i`$ or $`1p_i`$, respectively. We purposefully associate these hopping probabilities with a site because a site is associated with a particular codon. Thus, the jump rate from that site may depend on the relative abundance of the appropriate aa-tRNA. Apart from these probabilities, another aspect of our stochastics is random sequential updating: i.e., during each Monte Carlo step (MCS), $`M+1`$ particles are chosen at random, in sequence, to attempt moves. They are selected from a pool that includes the $`M`$ particles on the lattice plus another unbound particle that can initiate if there are $`\mathrm{}`$ holes at the beginning of the chain. Let us illustrate with a few examples. First, $`p_0`$ is associated with the start codon and, *if* the first $`\mathrm{}`$ sites are empty, a particle will be placed on the $`i=1`$ site with this probability. Next, a random particle (say at site $`i`$) is chosen and, provided it has a hole in front ($`n_{i+\mathrm{}}=0`$), will be moved with probability $`p_i`$. Naturally, it will not be necessary to check for headway for the “last” ribosome ($`\alpha =M`$ ). Finally, the stop codon will be associated with $`p_N`$. For simulations of the closed system, there will be no “beginning” or “end,” so that there are no special steps for initiation or termination.

In our computational studies, 100 identical systems of $`N`$ sites are simulated in parallel to obtain good statistics. Simulations of closed systems begin with particles evenly distributed around the ring and run for 3600 MCS to ensure that steady state is reached. Open systems begin empty and are run for 12,000 MCS (for $`N<500`$) or $`100N`$ MCS (for $`N500`$) to reach steady state. After steady state is attained, data including the current and density distribution can be collected. Density data is typically collected every 100 MCS. We often use continuous time Monte Carlo BKL  because it runs far more quickly than and provides the same results as standard Monte Carlo.

### I.2 Brief survey of known results

Extensive investigations of the simple totally asymmetric exclusion process (TASEP, defined as point particles hopping with unit rate along a line) with open boundaries can be found in the literature. Simulations have been performed Krug , and exact analytic results for the steady state exist DDM ; SD . Depending on the initiation (or injection) and termination (or depletion) rates, the system will settle into one of three phases. Introduced above as $`p_0`$ and $`p_N`$ respectively, the initiation and termination probabilities are mostly referred to as simply $`\alpha `$ and $`\beta `$ in the literature. From their dominant characteristics, the three phases are known as low density, high density, and maximal current. A phase diagram in this $`\alpha `$-$`\beta `$ plane has been determined, showing second order transitions between the maximal current phase and the others, as well as a first order transition between the high- and low-density regions. Subtle correlations further divide the last two into subregions. When disorder is introduced, i.e., not all the $`p_i`$’s are equal, then methods for exact analytic approaches fail (except in the extremely dilute limit, where only the motion of a single particle is of concern Derrida ). Indeed, even a single slow rate in a *closed* system poses serious difficulties JL1 ; JL2 ; Sep . However, Kolomeisky ABK  obtained approximate steady state solutions and phase diagrams for an *open* system with a single nonuniform rate in the bulk by splitting the system into two smaller systems connected by the nonuniform rate. Tripathy and Barma Barma  considered a closed system, but with a finite fraction of identical slow sites. Based on a combination of Monte Carlo simulations and numerical solutions of mean field equations, they found current-density relations. Though not a model for elongation, a related problem is “particlewise disorder,” in which the unequal hopping rates are associated with the particles rather than the sites Bengrine . Further references on TASEP with disorder may be found in a recent review JK-BJP . Also indirectly related to our one dimensional models are driven lattice gases with quenched disorder in *higher* dimensions LF ; BJ . Finally, we mention that there are many studies on the ASEP in which a particle has a finite probability of stepping backwards ASEP . Back steps are not generally believed to occur in elongation, and we will not consider such processes. All of these studies are restricted to $`\mathrm{}=1`$.

Systems with extended objects ($`\mathrm{}>1`$) have been rarely investigated, despite their introduction over three decades ago as a model for biopolymerization MGP . Using a mean field approach, MacDonald et al. set up mean field equations for the average site occupation $`n_i`$ and considered both closed MGP  and open MG  systems. In the former case, exact solutions were found, leading them to a current vs. density relation. For the latter, the authors resorted to numerical solutions to find the phase diagram for a variety of initiation and termination rates. A phase diagram similar to the simple TASEP, as well as non-trivial density profiles and the associated currents, was obtained. More recently, there is renewed interest in this problem. Naming this system “$`\mathrm{}`$-TASEP,” Sasamoto and Wadati SW  focused on the time dependence of $`M`$ particles in an infinite lattice and, using the Bethe ansatz, found exact results for the conditional probability that the particles are found at certain sites given an initial configuration. The main conclusion is that the dynamics of $`\mathrm{}`$-TASEP lies in the same universality class as ordinary TASEP. This line of inquiry has been further generalized to a system containing a *distribution* of particle sizes AB ; FA . Though these investigations produced interesting results, they are not applicable to our situation, namely, finite systems with open boundaries. In particular, for finite lattices, these studies are restricted to *closed* systems, for which the stationary states are trivial, as we will recapitulate in Section II. Finally, a recent work by Lakatos and Chou LC  considered uniform open systems with extended objects. Using a discrete Tonks gas partition function, they derived the current vs. density relation first presented by MacDonald et al. MGP . Via a refined mean field theory, they extended this result to predict currents and bulk densities for the open system, which they confirmed by Monte Carlo simulations. The phase diagram and its properties are consistent with those initially obtained by MacDonald and Gibbs MG . We are not aware of published results on open systems with both extended objects and quenched disorder.

## II TASEP of extended objects on a ring

For simplicity and mathematical reasons, it is convenient to discuss a uniform closed system with periodic boundary conditions. Here, $`M`$ identical particles of size $`\mathrm{}`$ move on a ring of $`N`$ sites. Due to the excluded volume constraint, there are $`\stackrel{~}{M}=N\mathrm{}M`$ uncovered sites (holes). The system evolves by discrete time steps, with random sequential updates. Specifically, during each MCS, $`M`$ particles are chosen at random in sequence, and each is moved forward one site provided there is headway. Since we have a stochastic dynamics, the complete description of this system is $`P(C,t)`$, the probability that it is found in configuration $`C`$ after $`t`$ steps (starting from some initial $`C_0`$). To label a configuration, we choose an arbitrary particle to be the first ($`\alpha =1`$) and supply $`\left\{h_\alpha \right\}`$, the set of number of holes in front of the $`\alpha ^{th}`$ particle (with $`\alpha =1,\mathrm{},M`$). Clearly, we may also think of a configuration as a series of $`M`$ “gaps” (between the particles) with $`h_\alpha `$ being the number of holes in the $`\alpha ^{th}`$ gap. So, $`P(h_1,h_2,\mathrm{}h_M;t)`$ is an explicit form for $`P(C,t)`$. Note that, since the system is closed, there is a constraint on $`h_\alpha `$, i.e.,

$$\underset{\alpha }{}h_\alpha =\stackrel{~}{M}$$

is a constant in time.

Random sequential updating can be translated into an equation which governs the time evolution of $`P(C,t)`$, namely a master equation. Starting with the initial $`P(C,0)=\delta (C,C_0)`$ (where $`\delta `$ is the Kronecker delta), $`P(C,t)`$ is expected to settle into a unique, *time-independent* distribution, $`P^{}\left(C\right)`$, which we will refer to as the “steady state.” If this system were evolving towards thermal equilibrium, the dynamics would satisfy detailed balance DB , and $`P^{}\left(C\right)`$ would be given by the well-known Boltzmann factor. However, the dynamics of our system definitely violates detailed balance, so that, associated with a *non-equilibrium* steady state, $`P^{}\left(C\right)`$ is not known in general. Fortunately, this closed system belongs to a class for which a simple solution is known Spitzer ; namely, every configuration occurs with equal probability pairwisenote . Thus, $`P^{}\left(C\right)`$ is precisely the reciprocal of the total number of configurations consistent with the given parameters ($`N,M,\mathrm{}`$). With such a simple distribution, we can compute many quantities of interest, such as the probability distribution of the current (and hence the average current) and the headway.

Apart from overall factors, the total number of configurations is just $`Z(\stackrel{~}{M},M),`$ the total number of lists $`\left\{h_\alpha \right\}`$ subject to Eqn. (II). Because this is a well-known combinatorial problem (appearing in, e.g., the Bose gas), we simply quote the result:

$$Z(\stackrel{~}{M},M)=\frac{(\stackrel{~}{M}+M1)!}{\stackrel{~}{M}!(M1)!}\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}+M1}{\stackrel{~}{M}}\right)$$

The actual number of distinct configurations is $`\left(N/M\right)Z(\stackrel{~}{M},M),`$ because there are $`N`$ lattice sites on which the first particle can be placed but the $`M`$ particles are identical. Thus,

$$P^{}\left(C\right)=\frac{\stackrel{~}{M}!M!}{N(\stackrel{~}{M}+M1)!}$$

*independent* of $`C`$.

To compute the probability distribution of the current and its average, we need more detailed information on the above “partition.” Defining the current, $`J`$, as the number of particles which moved in one step (normalized by the system size $`N`$), we see that we will need $`H`$, the number of gaps with one or more holes. The reason is that, for each such gap, the ribosome behind it can move and contributes one “unit” to the current, so that

$$J=\frac{H}{N}.$$

Since all configurations are equally probable, the statistical weight associated with this $`J`$ is just the total number of configurations with a given $`H`$. This quantity, denoted by $`Z(H;\stackrel{~}{M},M)`$ in analog to $`Z(\stackrel{~}{M},M)`$, may be found from its definition:

$`Z(H;\stackrel{~}{M},M)`$ $`=`$ $`{\displaystyle \underset{\left\{h_\alpha \right\}}{}}[\delta (\stackrel{~}{M},{\displaystyle \underset{\alpha }{}}h_\alpha )\times `$
$`\delta (H,{\displaystyle \underset{\alpha }{}}[1\delta (h_\alpha ,0)])],`$
where the sum is over all possible lists $`\left\{h_\alpha \right\}`$ and the Kronecker $`\delta `$’s select only those which satisfy Eqn. (II) and have $`H`$ gaps with $`h_\alpha >0`$. Once $`Z(H;\stackrel{~}{M},M)`$ is known, the full distribution for the current is

$$p(J;\stackrel{~}{M},M)=\frac{Z(H;\stackrel{~}{M},M)}{Z(\stackrel{~}{M},M)}.$$

Now, the explicit form of $`Z(H;\stackrel{~}{M},M))`$ can be obtained either through the generating function

$`W_M(\zeta ,\eta )`$ $``$ $`{\displaystyle \underset{\stackrel{~}{M},H}{}}Z(H;\stackrel{~}{M},M)\zeta ^{\stackrel{~}{M}}\eta ^H`$
$`=`$ $`\left[1+{\displaystyle \frac{\zeta \eta }{1\zeta }}\right]^M`$
or by standard combinatorial techniques:

$$Z(H;\stackrel{~}{M},M)=\left(\genfrac{}{}{0pt}{}{M}{H}\right)\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}1}{H1}\right).$$

Thus, the explicit current distribution is

$$p(J)=\left(\genfrac{}{}{0pt}{}{M}{H}\right)\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}1}{H1}\right)/\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}+M1}{\stackrel{~}{M}}\right),$$
(2)
where $`H`$ on the right stands for $`JN.`$ An alternate form, showing the dependence of this distribution on the control parameters ($`\stackrel{~}{M},M,N=\stackrel{~}{M}+\mathrm{}M`$) is

$`p(J|\stackrel{~}{M},M,N)`$ $`=`$ $`{\displaystyle \frac{1}{M\stackrel{~}{M}\left(\stackrel{~}{M}+M1\right)!}}\times `$
$`{\displaystyle \frac{JN}{(\stackrel{~}{M}JN)!(MJN)!}}\times `$
$`\left[{\displaystyle \frac{M!\stackrel{~}{M}!}{\left(JN\right)!}}\right]^2.`$
To illustrate this distribution, we show in Figure 1 both this prediction and simulation data for the case of $`N=200`$, $`M=15`$, and $`\mathrm{}=12`$. Clearly, there is excellent agreement between theory and simulation.

Containing less information, but easier to grasp, is the average current $`\overline{J}Jp\left(J\right)`$. Its computation is somewhat easier, since it is $`_{H=1}^MHZ(H;\stackrel{~}{M},M)/NZ(\stackrel{~}{M}.M)`$, with the numerator easily gleaned from $`_\eta W_M(\zeta ,\eta )|_{\eta =1}.`$ The result is

$$\overline{J}=\frac{M}{N}\frac{\stackrel{~}{M}}{\stackrel{~}{M}+M1}.$$
(3)
As we will see, the dependence of this average current on the density of particles plays a central role. Expressing this quantity in terms of $`\rho _rM/N`$, we have $`\overline{J}=\rho _r\left(1\mathrm{}\rho _r\right)/\left[1\left(\mathrm{}1\right)\rho _r1/N\right].`$ An appealing form, which displays both the intensive nature of $`J`$ and its underlying particle-hole symmetry, is

$$\overline{J}=\frac{\rho _r\rho _h}{\rho _s1/N}.$$
(4)
A third form, frequently referred to in the literature as the “current-density relationship,” is writing $`\overline{J}`$ as a function of $`\rho `$, the coverage density ($`\rho [0,1]`$):

$$\overline{J}(\rho )=\frac{\rho }{\mathrm{}}\frac{1\rho }{1\rho +\rho /\mathrm{}1/N}.$$
(5)
In the limit $`N\mathrm{}`$, this result was first presented by MacDonald et al. MGP . A generalization of the well-known expression for $`\mathrm{}=1`$ (i.e., $`\overline{J}=\rho (1\rho )`$), this $`\overline{J}(\rho )`$ is no longer symmetric about $`\rho =1/2`$. Instead, the optimal density increases from $`1/2`$ to $`\sqrt{\mathrm{}}/\left(1+\sqrt{\mathrm{}}\right)`$, while the maximum current is lowered from $`1/4`$ to $`\left(1+\sqrt{\mathrm{}}\right)^2`$. As these quantities will appear frequently, we will denote them by

$$\widehat{\rho }\frac{\sqrt{\mathrm{}}}{1+\sqrt{\mathrm{}}}\text{and}\widehat{J}\left(1+\sqrt{\mathrm{}}\right)^2.$$
(6)
To appreciate this shift graphically, we present, in Figure 2, both analytic and simulation results. In addition, to show the effects of the finite size corrections (due to the $`1/N`$ term), we include a curve of the limiting form MGP ; LC

$$\overline{J}\frac{\rho _r\rho _h}{\rho _s}$$
(7)
for the $`N=40`$ case.

In connection with these expressions, a conclusion for  *conditional* probabilities can be drawn. Since $`\overline{J}`$ is precisely the joint probability of finding a “covered”-hole pair, we see that $`\rho _r/\rho _s`$ is the probability that site $`i`$ is covered, given that site $`i+1`$ is empty and, similarly, $`\rho _h/\rho _s`$ is the probability that site $`i+1`$ is empty, given that site $`i`$ is covered. These conditional probabilities will play a role in our understanding of the behavior of open systems.

In addition to the average current, we can also compute its fluctuations exactly.

$$\mathrm{\Delta }J^2J^2p\left(J\right)\overline{J}^2=\frac{\overline{J}^2}{\rho _sN}.$$

Typical of non-critical thermodynamic systems, in which the (fractional) deviations are $`O\left(N^{1/2}\right)`$, the full distribution $`p(J)`$ approaches the standard Gaussian form: $`\mathrm{exp}\left[\frac{N}{2}\rho _s\left(J/\overline{J}1\right)^2\right]`$. Indeed, this is the form we see in the example shown in Figure 1 above.

Finally, we turn to another quantity of interest: the statistics of “headway.” Since there are $`M`$ gaps in each configuration and there are $`Z(\stackrel{~}{M},M)`$ configurations, we have a total of $`MZ(\stackrel{~}{M},M)`$ gaps. Out of these, we wish to compute the number of gaps which contain precisely $`m`$ holes, which we denote by $`Z(m;\stackrel{~}{M},M)`$. From its definition

$$Z(m;\stackrel{~}{M},M)\underset{\left\{h_\alpha \right\}}{}\delta (\stackrel{~}{M},\underset{\alpha }{}h_\alpha )\left(\underset{\alpha }{}\delta (m,h_\alpha )\right),$$

we find the associated generating function

$`W(\zeta ,\eta )`$ $``$ $`{\displaystyle \underset{\stackrel{~}{M},m}{}}Z(m;\stackrel{~}{M},M)\zeta ^{\stackrel{~}{M}}\eta ^m`$
$`=`$ $`M\left({\displaystyle \frac{1}{1\zeta \eta }}\right)\left({\displaystyle \frac{1}{1\zeta }}\right)^{M1},`$
leading to the distribution of head spacings

$`p(m;\stackrel{~}{M},M)`$ $``$ $`{\displaystyle \frac{Z(m;\stackrel{~}{M},M)}{MZ(\stackrel{~}{M},M)}}`$ (8)
$`=`$ $`{\displaystyle \frac{\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}+Mm2}{\stackrel{~}{M}m}\right)}{\left(\genfrac{}{}{0pt}{}{\stackrel{~}{M}+M1}{\stackrel{~}{M}}\right)}}.`$
In the limit of large $`N`$, this expression simplifies to

$$p(m)\frac{\rho _r}{\rho _s}\left(\frac{\rho _h}{\rho _s}\right)^m.$$

This distribution is reproduced faithfully in Monte Carlo simulations, as shown by an example in Figure 3. It is easy to understand this result intuitively if we regard $`\rho _h/\rho _s`$ as the probability of having a single hole in the headway. With independent hole statistics, we have $`p(m)\left(\rho _h/\rho _s\right)^m`$. The average number of holes in the headway and the associated standard deviation can be computed easily: $`\rho _h/\rho _r`$ and $`\sqrt{\left(\rho _h/\rho _r\right)^2+\rho _h/\rho _r}`$, respectively.

Though the system considered in this section, particles traveling on a ring, bears little resemblance to the translation process, it is sufficiently simple for us to derive a number of exact results. Apart from their own interest, these results provide crucial insights, such as the current-density relationship, for postulating appropriate equations in a coarse grained, mean field approach to the physical problem at hand.

## III Extended objects in open systems

In this section, we attempt the next step towards a realistic model for protein synthesis by considering systems with open boundaries. The section is organized as follows. We describe the problem and introduce some terminology. Next, we study the effects of the open boundaries alone, by analyzing the continuum limit of a symmetric exclusion process. Finally, we consider the asymmetric exclusion process with open boundaries. We present the phase diagram found from simulations and an extremal principle analysis, and we show the ability of the continuum limit to predict steady state density profiles.

In the open system, the first site ($`i=1`$) is no longer “in front of” the $`i=N`$ site. Instead, a particle (of extent $`\mathrm{}`$) will be placed at $`i=1`$ with probability $`\alpha `$ (previously labeled by $`p_0`$), *provided* allthe first $`\mathrm{}`$ sites are empty. This models the initiation process. As for elongation, we continue to restrict ourselves to uniform rates here. In the language of TASEP, every (randomly) chosen particle *will move* by one site (increasing $`i`$ by unity) if it has some headway. For particles on the last $`\mathrm{}`$ sites, there is no hinderance, so they will always move if chosen. Finally, to simulate termination, a particle on the $`N^{th}`$ site will be removed from the system with rate $`\beta `$ (previously labeled by $`p_N`$). To repeat, in a Monte Carlo step (MCS), $`M+1`$ particles are chosen randomly (in sequence) to attempt a move. They are selected from a pool including the $`M`$ particles on the lattice plus an unbound particle that may initiate.

Since this system is no longer closed, $`M`$ and $`\stackrel{~}{M}`$ are fluctuating quantities. Of course, their average values will be controlled by the rates $`\alpha `$ and $`\beta `$. Our goal is to find the average densities and the average current of such a model, as functions of $`(\alpha ,\beta )`$. In the $`\mathrm{}=1`$ case, it is known that the system exhibits three different “phases” as $`\alpha `$ and $`\beta `$ are varied, only one of which resembles the closed system above (in the sense that the current approaches $`\widehat{J}`$ for large $`N`$). Beyond overall averages, we seek the density profile, which will not only be non-trivial, but which also displays drastically different properties as we move about in the $`\alpha `$-$`\beta `$ plane. The main goal of this paper is to study the effects of $`\mathrm{}>1`$ on both the phase diagram and these density profiles. Unfortunately, $`P^{}\left(C\right)`$ for a system with open boundaries is not known in general. Even for the $`\mathrm{}=1`$ case, only a limited set of quantities may be computed exactly. To make progress, we resort here to a more phenomenological approach, in the spirit of hydrodynamics or Lanudau-Ginzburg free-energy functionals. Considering coarse-grained densities and the continuum limit, we postulate equations of motion, based on some of the properties of the closed system. In principle, such equations can be “derived” from the master equation for $`P(C,t)`$, using the mean-field approximation (i.e., ignoring all correlations). In a sense, this is also the approach of MacDonald et al. MGP , except that they focused on $`n_i^{\left(L\right)}\left(t\right)`$ and $`n_i^{\left(0\right)}\left(t\right),`$ the probabilities for site $`i`$ to be occupied by a ribosome and a hole, respectively, at time $`t`$. By keeping the spatial co-ordinate discrete, they faced difference equations and succeeded in finding solutions only numerically. In contrast, by using a continuous spatial co-ordinate, we have ordinary differential equations instead, giving us *analytic* solutions. However, lost in the notion of “coarse-graining” are the period $`\mathrm{}`$ structures which feature prominently behind “blockages.” Nevertheless, our approach appears to capture the essential (gross) features of these systems.

In our heuristic approach, we imagine $`N`$ to be large enough to justify taking the continuum limit, i.e., replacing the discrete site label $`i`$ by a continuous co-ordinate: $`x`$. For simplicity, define

$$xi/N$$

so that $`x`$ lies within the unit interval. (If physical units of length are desired, we may introduce $`a`$ as the lattice spacing, corresponding to the length of a codon, i.e., three bases of the mRNA. Then $`L_{mRNA}Na`$ would be the length of the mRNA in question.) Similarly, continuous *local* densities will take the place of the discrete occupation variables. For example, the coverage density $`\rho \left(x\right)`$ will be used instead of $`n_i`$. Since the maximum occupancy is unity, the hole density is just

$$\rho _h\left(x\right)=1\rho \left(x\right).$$

Following our considerations above, we also define the ribosome density by

$$\rho _r\left(x\right)=\rho \left(x\right)/\mathrm{}.$$

Note that, despite the continuum limit, these equations display the meaning of $`\mathrm{}`$, which serves as a measure of the “size” (or “extent”) of ribosome. Of course, we are also interested in their time dependence, so that we must consider $`\rho (x,t)`$ in general. Now, ribosomes rarely detach from the mRNA during elongation, so that we are justified in regarding these densities as *conserved* fields. Thus, the appropriate equation of motion is the continuity equation, i.e.,

$$\frac{}{t}\rho _r(x,t)=\stackrel{}{}\stackrel{}{J}_r=\frac{}{x}J_r(x,t),$$
(9)
where $`J_r(x,t)`$ is the local (ribosome) current. Our first task is to find how this current depends on the (local) density $`\rho _r(x,t)`$, i.e., to find the functional form for $`J_r\left[\rho _r\right]`$. Then we will arrive at an acceptable equation for $`\rho _r`$. To find the steady state profile $`\rho _r^{}\left(x\right)`$, which satisfies $`_t\rho _r^{}=0`$, we see that this state is associated with a constant ($`x,t`$-independent) current $`J_r^{}`$. Therefore, our problem consists of finding the solution to

$$J_r\left[\rho _r^{}\right]=constant,$$

subject to the appropriate boundary conditions. To keep the notation simple, from here on we will drop the subscript $`r`$ and write $`J`$ or $`J(x,t)`$ for the local current (as well as $`J^{}`$ for the constant current in the stationary state).

### III.1 The case of free diffusion and an effective density variable for extended objects

To show how we build an appropriate equation for the TASEP, let us begin with a study of the effects of open boundaries *alone*. In other words, let us consider a *symmetric* exclusion process for extended objects, i.e., a system of large particles diffusing freely on a line, subjected to the excluded volume constraint and the controls $`\alpha ,\beta `$. For simulations, a randomly chosen particle is moved one site forward or backward with equal probability ($`0.5`$), provided it does not run into its neighbor. The only exception is the first particle, which is prohibited from jumping backwards into the “source.”

For the $`\mathrm{}=1`$ case ($`\rho =\rho _r`$), the steady state profile is trivially linear:

$$\rho ^{}\left(x\right)=\left(1x\right)\rho ^{}\left(0\right)+\left(x\right)\rho ^{}\left(1\right)$$

with current

$$J^{}=\frac{1}{N}\left[\rho ^{}\left(0\right)\rho ^{}\left(1\right)\right].$$

(Note that, since the current is controlled by the gradient of the local density only, it vanishes necessarily in the $`N\mathrm{}`$ limit. So, we must keep the $`N`$ explicitly here.) Meanwhile, the boundary densities are fixed by matching the injection/depletion rates to the internal current:

$$\alpha \rho _h^{}\left(0\right)=\alpha \left[1\rho ^{}\left(0\right)\right]=J^{}=\beta \rho ^{}\left(1\right),$$

and our problem is completely solved. These well known results can be traced to the fact that Eqn. (9) assumes the form of the simple diffusion equation $`_t\rho _x^2\rho `$, since $`J\left[\rho \right]\rho `$. Though there are many ways to arrive at this result, it is less clear how to generalize it to the case of extended objects. In particular, as displayed in Figure 4, the profile from simulations (for $`N=200,\mathrm{}=12,\alpha =\beta =1`$) is far from linear! Here, we will show that there is a natural generalization for the (particle) density in the case of extended objects and that its profile is again linear.

We proceed by returning to the basics, starting with the current $`J_r`$ being proportional to both the conductivity (mobility) and the drive. For the former, we take the result from the previous section. There, the drive is constant, so that the right hand side of Eqn. (7), $`\rho _r\rho _h/\rho _s`$, can be interpreted as the conductivity or mobility (neglecting the finite-size effect term from the closed system). Meanwhile, the drive for free diffusion should be the gradient of a pressure $`𝒫`$, so that

$$J_r\left[\rho \right]=D\left[\frac{\rho _r\rho _h}{\rho _s}\right]\left[𝒫\right],$$
(10)
where $`D`$ is a constant, to be fitted to data. Since we have scaled the system size to unity, we should keep in mind that $`D`$ is a quantity of $`O\left(1/N\right)`$, to be consistent with the continuum limit. For the pressure $`𝒫`$, we follow the standard route of statistical mechanics and write

$$𝒫=\frac{\delta }{\delta \rho _r},$$

where $``$ is a free energy functional (e.g., the Landau-Ginzburg “Hamiltonian” in case of the ordinary lattice gas). Here, we have a non-interacting system, so that a reasonable form for $``$ is just the entropy entropynote :

$$=𝑑x\left[\rho _r\mathrm{ln}\rho _r+\rho _h\mathrm{ln}\rho _h\rho _s\mathrm{ln}\rho _s\right].$$

Using

$$\rho _h=1\mathrm{}\rho _r,\rho _s=1(\mathrm{}1)\rho _r$$

and carrying out the steps, we arrive at a *modified* diffusion equation:

$$_t\rho _r=D_x\left[\rho _s^2_x\rho _r\right].$$
(11)
Note that, for $`\mathrm{}>1`$, the effective “diffusion constant,” $`D/\rho _s^2`$, is *density-dependent*. Thus, the steady state profile will not be linear in $`x`$. Instead, it satisfies

$$\frac{D}{\left(\rho _s^{}\right)^2}\frac{\rho _r^{}}{x}=J^{},$$
(12)
which will definitely lead to non-linear profiles.

Although we have derived this diffusion equation via phenomenological techniques, we speculate that it could be obtained through a more rigorous derivation. Specifically, if a matrix product description of the $`\mathrm{}>1`$ system becomes available, the methods of Derrida et al. DLS1  may extend to calculate a large deviation functional, which would play the role of a nonequilibrium free energy for this system. Should such a functional exist, we expect that our diffusion equation would result from a variation of this functional.

To find $`\rho _r^{}`$ explicitly, notice that, due to $`\rho _s=1\left(\mathrm{}1\right)\rho _r`$, the left hand side is just a simple derivative of $`1/\rho _s^{}`$. Thus, the profile $`1/\rho _s^{}\left(x\right)`$ will be *linear.* However, $`\rho _s`$ does not reduce to a sensible variable as $`\mathrm{}1`$, leading us to define

$$\chi \frac{1}{\mathrm{}1}\left[\frac{1}{\rho _s^{}}1\right]=\frac{\rho ^{}}{\mathrm{}\left(\mathrm{}1\right)\rho ^{}}.$$

Not only does $`\chi `$ lie in the unit interval; it reduces to the usual particle density $`\rho ^{}`$ in the limit $`\mathrm{}1`$. Most importantly, it is uniquely related to $`\rho ^{}`$ for all $`\mathrm{}`$. First introduced in FA , $`\chi `$ is the “natural” generalization of the particle density for extended objects. We will refer to it as the effective particle density (EPD). Meanwhile, a hole remains a single-site entity, so its density needs no modification. In terms of the EPD, the various densities are

$`\rho _s^{}={\displaystyle \frac{1}{1+\left(\mathrm{}1\right)\chi }};`$
$`\rho _r^{}={\displaystyle \frac{\chi }{1+\left(\mathrm{}1\right)\chi }};`$
$`\rho _h^{}={\displaystyle \frac{1\chi }{1+\left(\mathrm{}1\right)\chi }}`$
so that $`\chi `$ is just

$$\chi \rho _r^{}/\rho _s^{}.$$
(13)
The current-density relationship (Eqn. 7) for the steady state in the *ring* is again a simple product:

$$\overline{J}_{ring}=\chi \rho _h.$$
(14)
Note that this product will serve as the *mobility* in our mean field approach here and can be expressed in terms of $`\chi `$ alone:

$$\chi \rho _h=\frac{\chi \left(1\chi \right)}{1+\left(\mathrm{}1\right)\chi }.$$

The main advantage of using the EPD is the re-emergence of the familiar combination $`\chi \left(1\chi \right)`$ in the numerator. Finally, corresponding to the optimal density (Eqn. 6), we have

$$\widehat{\chi }=\frac{1}{1+\sqrt{\mathrm{}}},$$

a quantity which plays a significant role in systems with open boundaries. Note that the corresponding hole density ($`1\widehat{\rho }`$) also assumes the same value, so that $`\widehat{J}`$ is just $`\widehat{\chi }^2`$. The underlying “particle-hole” symmetry generalizes to one under exchange of $`\chi `$ $`\rho _h`$, or

$$\chi \frac{\left(1\chi \right)}{1+\left(\mathrm{}1\right)\chi },$$

or, in terms of the more physical $`\rho `$,

$$\rho \frac{1\rho }{1\rho +\rho /\mathrm{}}.$$

We now find that Eqn. (12) simplifies to the familiar equation from ordinary diffusion:

$$D\frac{\chi }{x}=J^{}.$$

As a result, the profile in terms of $`\chi `$ is again linear in $`x`$ linearnote . For completeness, we write the solution:

$$\chi \left(x\right)=\left(1x\right)\chi \left(0\right)+\left(x\right)\chi \left(1\right),$$

with

$$\chi \left(0\right)\chi \left(1\right)=J^{}/D.$$

The explicit solution ensues once the constraints of injection/depletion rates are imposed. Note that, in general, we would have $`J^{}D`$, so that the current inherits the $`O\left(1/N\right)`$ from $`D`$, in contrast to the $`O\left(1\right)`$ behavior for systems with non-zero drive.

To close, we illustrate how well the theory agrees with data by showing the analytic, non-linear profile of $`\rho ^{}\left(x\right)`$ with boundary conditions $`\rho \left(0\right)=1`$, $`\rho \left(1\right)=0`$:

$$\rho ^{}\left(x\right)=\frac{1x}{1\left(11/\mathrm{}\right)x},$$

as a curve in Figure 4.

### III.2 TASEP in open systems

We now turn to the other extreme case, where backward jumps are completely excluded (TASEP), generalizing the model of DDM  to extended objects. With open boundaries, this model incorporates another essential feature of translation. Our interest is again the average current and density profile, as a function of $`\alpha `$ and $`\beta `$, the rates of initiation and termination.

#### III.2.1 Phase diagram from extremal principle and simulations

Although we are unable to generalize the methods of DDM ; SD  to arrive at exact solutions for an open system with particles covering $`\mathrm{}>1`$ sites, we find that the phase diagram determined by the mean field techniques of MG ; LC  can be understood also by extending the extremal principle, an hypothesis first proposed by Popkov and Schütz PS  and later exploited successfully by others (e.g., Hager1 ; Hager2 ). In this approach, the open boundaries are regarded as connections to reservoirs with appropriate densities, so that by keeping the same jump rates as in the bulk, $`\alpha `$ and $`\beta `$ are realized. Defining $`\rho _{}`$ and $`\rho _+`$ as the reservoir densities at the initiation and termination boundaries, respectively, the extremal principle relates the current in the open system to the $`J[\rho ]`$ for a closed, periodic system with the same bulk dynamics inflimitnote :

$$J=\{\begin{array}{cc}\mathrm{max}J[\rho ]\hfill & \text{for }\rho _+<\rho <\rho _{}\hfill \\ \mathrm{min}J[\rho ]\hfill & \text{for }\rho _{}<\rho <\rho _+\hfill \end{array}.$$

Unfortunately, there is no prescription for finding $`\rho _{},\rho _+`$ from $`\alpha ,\beta `$ in general. Exploiting results from the exactly soluable $`\mathrm{}=1`$ case (where $`\rho _{}=\alpha `$ and $`\rho _+=1\beta `$ ) and from previous studies of the $`\mathrm{}>1`$ case MG ; LC , we argue in favor of

$$\rho _{}(\alpha )=\frac{\mathrm{}\alpha }{1+\alpha \left(\mathrm{}1\right)}\text{and}\rho _+(\beta )=1\beta .$$
(15)
These reservoir densities can be understood as follows. Recall that, for a closed system, the probability for site $`i1`$ to be filled, *given* that site $`i`$ is empty, is $`\rho _r/\rho _s=\left(\rho /\mathrm{}\right)/\left(1\rho +\rho /\mathrm{}\right)`$. We now argue that when the first site of an open system is empty and will be filled with probability $`\alpha `$, it can be thought of as being coupled to a reservoir of the appropriate density, i.e., $`\rho _{}`$ such that $`\alpha =\left(\rho _{}/\mathrm{}\right)/\left(1\rho _{}+\rho _{}/\mathrm{}\right)`$. Solving for $`\rho _{}`$ leads to the expression above for $`\rho _{}(\alpha )`$. The expression for $`\rho _+\left(\beta \right)`$ is not readily explained by similar arguments, so discussion of its origin will be deferred until the following section. An important feature associated with this choice is that the current is symmetric under $`\alpha \beta `$, as observed in simulations. Most importantly, these choices lead to phase diagrams in good agreement with Monte Carlo data.

By combining these choices with the extremal principle, we find that, although the result is qualitatively similar to the $`\mathrm{}=1`$ system, there are quantitative changes for the $`\mathrm{}>1`$ case. First, there is a shift in the location of transition lines, from $`1/2`$ to $`\widehat{\chi }=1/\left(1+\sqrt{\mathrm{}}\right)`$ (as shown in Figure 5). Then, the current and bulk densities are also modified:

$$\overline{J}(\alpha ,\beta )=\{\begin{array}{cc}\frac{\alpha (1\alpha )}{1+\alpha (\mathrm{}1)}\hfill & \text{for }\alpha <\beta <\widehat{\chi }\text{ (low density)}\hfill \\ \frac{\beta (1\beta )}{1+\beta (\mathrm{}1)}\hfill & \text{for }\beta \alpha <\widehat{\chi }\text{ (high density)}\hfill \\ \widehat{J}\hfill & \text{for }\alpha ,\beta \widehat{\chi }\text{ (max current)}\text{ ,}\hfill \end{array}$$
(16)
and

$$\overline{\rho }(\alpha ,\beta )=\{\begin{array}{cc}\rho _{}\hfill & \text{for }\alpha <\beta <\widehat{\chi }\text{ (low density)}\hfill \\ \rho _+\hfill & \text{for }\beta <\alpha <\widehat{\chi }\text{ (high density)}\hfill \\ \widehat{\rho }\hfill & \text{for }\alpha ,\beta \widehat{\chi }\text{ (max current)}\text{ .}\hfill \end{array}$$
(17)
As examples of how well these predictions fit simulation data, we show plots of $`\overline{J}(1,\beta )`$ (Figure 6) and $`\overline{\rho }(\alpha ,0.1)`$ (Figure 7). As expected, when $`\alpha `$ was rate-limiting (low density phase), a bulk coverage density of $`\rho _{}(\alpha )`$ was induced, and when $`\beta `$ was rate-limiting, a bulk coverage density of $`\rho _+(\beta )`$ was induced. Given the success of the extremal principle analysis, it would be appropriate to search for more fundamental theories from which the extremal principle could be derived.

The transition between low and high density phases is clearly first order, displayed as a jump in the bulk density plot (Figure 7). As in the $`\mathrm{}=1`$ case KSKS , domain wall theory can be used to understand our results. In particular, for $`\alpha =\beta <\widehat{\chi }`$, we have observed a shock front (between the low and high density regions) diffusing freely along the lattice. As a result, the *average* profile is linear, as shown in Figure 8. To appreciate the shock, we display in the inset a typical configuration (showing a shock) where each ribosome is represented by a dot. In addition, we exhibit another aspect of this co-existence by sampling the average (coverage) density in the central 10% of a large ($`N=1000`$) lattice. Compiling a histogram (closed diamonds in Figure 9), we find a bimodal distribution typical of systems at a first order transition. Peaks are expected to correspond to the densities $`\rho _{}\left(\alpha \right)`$ and $`\rho _+\left(\beta \right)`$. To connect with local fluctuations in the “pure” systems, we compile similar histograms for a closed system (i.e., 10% of a ring with $`1000`$ sites) with overall density set at $`\rho _{}(\alpha )`$ and $`\rho _+(\beta )`$. For comparison, we show a simple linear combination of such distributions in Figure 9 (open squares). We expect correspondence between densities in the window in the open system, which alternate between approximately $`\rho _{}\left(\alpha \right)`$ and $`\rho _+\left(\beta \right)`$, and an appropriately weighted average of window densities in the two closed systems. The deviations can be understood as the contribution of configurations with the shock in the sampling window. Roughly, this may occur about 10% of the time, which is also the order of magnitude of the deviation from the simple average of pure systems.

While most of the features presented here were known to MacDonald et al. MGP ; MG  and have also been obtained by Lakatos and Chou LC , our efforts are to go beyond mean field approaches, showing the different perspective offered by domain wall theory and the extremal principle. Shock effects on the co-existence line are particularly well described by this approach.

#### III.2.2 Steady state profiles from the continuum approach

Following an understanding of the phase diagram, we turn to more details of the system, namely, steady state density profiles. Our continuum approach again leads us to differential equations. As considerable insight can be gained by studying simple differential equations, we believe it is worthwhile to devote a section to the continuum approach.

Starting with the case of free diffusion, we extend Eqn. (11) to the driven system by simply adding an “Ohmic” term to the current: $`J=\sigma E`$, where the drive is associated with strength $`E`$ and the conductivity will be the density dependent factor found before (first used in Eqn. 10). Thus, our starting equation is

$$\frac{\rho _r}{t}=\frac{}{x}J\left(x\right)=\frac{}{x}\left[\frac{D}{\rho _s^2}\frac{\rho _r}{x}\frac{\rho _r\rho _h}{\rho _s}E\right].$$
(18)
As before, we suggest that this driven diffusion equation might be obtained from a variation of an appropriate large deviation functional, the $`\mathrm{}=1`$ version of which has been derived previously DLS2 ; DLS3 .

As a reminder, our $`x`$ is actually the *fractional* length along the mRNA, so that it is dimensionless. Similarly, the densities are also unitless (e.g., $`\rho _h[0,1]`$), so that $`D/E`$ carries the information of the real length of the chain. An estimate of this ratio can be obtained by taking the naïve continuum limit of a discrete hopping model, with the result $`1/2N`$. In the Appendix, we will show how this arises from taking such a limit of the MacDonald et al. MGP  current equation. For here, we may regard $`D/E`$ as a phenomenological parameter. Further, since we are focusing only on *totally* asymmetric processes, we could just as well set $`E`$ to unity (as in Eqn. 7), and measure $`J`$ and $`D`$ “in units of $`E`$.”

Taking these considerations into account, we seek the steady state profile $`\rho _r^{}\left(x\right)`$, as in Eqn. (12), by setting the square bracketed terms in Eqn. (18) to the steady state current:

$$D\left(\rho _s^{}\right)^2\frac{\rho _r^{}}{x}\rho _r^{}\rho _h^{}\left(\rho _s^{}\right)^1=J^{}.$$

In terms of our “natural” variable $`\chi `$ (Eqn. 13), this equation reduces to

$`\chi ^{}`$ $``$ $`{\displaystyle \frac{\chi }{x}}`$ (19)
$`=`$ $`{\displaystyle \frac{1}{D\left[1+\left(\mathrm{}1\right)\chi \right]}}\left[\stackrel{~}{J}\chi +J^{}\chi \left(1\chi \right)\right],`$
where $`\stackrel{~}{J}J^{}\left(\mathrm{}1\right)`$. Thus, we see that this $`\mathrm{}>1`$ generalization is quite similar to the $`\mathrm{}=1`$ equation: $`\rho ^{}=\left[\rho \left(1\rho \right)J^{}\right]/D`$.

Now, the zeros of $`\chi ^{}`$ will play an important role, occurring at

$$\chi _>\frac{1\stackrel{~}{J}+R}{2}\text{ and }\chi _<\frac{1\stackrel{~}{J}R}{2},$$
(20)
where

$$R\sqrt{\left(1\stackrel{~}{J}\right)^24J^{}}.$$

It is not surprising that the maximal current $`\widehat{J}`$ (i.e., $`\left(1+\sqrt{\mathrm{}}\right)^2`$ in our system of units) is a key player, so that $`R`$ can be written in the following form:

$$R=\left(\mathrm{}1\right)\sqrt{\left(J^{}\widehat{J}\right)\left(J^{}\stackrel{ˇ}{J}\right)}.$$

Here, $`\stackrel{ˇ}{J}`$ $`\left(1\sqrt{\mathrm{}}\right)^2`$ appears as a natural counterpart to $`\widehat{J}`$. Another advantage of this form is its relationship to the extremal current principle. As we will see, for finite systems, $`J^{}`$ will be slightly ($`O\left(1/N^2\right)`$) larger than $`\widehat{J}`$ for the maximal current phase, giving us complex roots and profiles with inflections. On the other hand, for the other two phases, $`J^{}<\widehat{J}`$, leading to fixed points in the $`N\mathrm{}`$ limit.

Eqn. (19) can be integrated to find an implicit function for the density profile in each phase. For explicit solutions, we must specify the boundary conditions:

$`J^{}`$ $`=`$ $`\alpha \left[1\rho \left(0\right)\right]`$
$`J^{}`$ $`=`$ $`{\displaystyle \frac{\beta }{1+\beta \left(\mathrm{}1\right)}}\rho \left(N\right)`$ (21)
so that

$$\chi \left(0\right)=\frac{\rho \left(0\right)/\mathrm{}}{1\rho \left(0\right)+\rho \left(0\right)/\mathrm{}}=\frac{\alpha J^{}}{\alpha +\left(\mathrm{}1\right)J^{}}$$
(22)
and

$`\chi \left(1\right)`$ $`=`$ $`{\displaystyle \frac{\rho \left(N\right)/\mathrm{}}{1\rho \left(N\right)+\rho \left(N\right)/\mathrm{}}}`$ (23)
$`=`$ $`{\displaystyle \frac{\left[1+\beta \left(\mathrm{}1\right)\right]J^{}}{\mathrm{}\beta \left(\mathrm{}1\right)\left[1+\beta \left(\mathrm{}1\right)\right]J^{}}}.`$
The boundary conditions can be understood as follows. From $`\rho _{}\left(\alpha \right)`$ in Eqn. (15), we can show that the effective reservoir at the initiation boundary has an EPD of $`\chi _{}\left(\alpha \right)=\alpha `$. Eqn. (14) then implies the initiation boundary condition in Eqn. (21). Particles on the final $`\mathrm{}`$ lattice sites experience no steric hindrance, so the current through the final $`\mathrm{}`$ sites is just

$`J^{}`$ $`=`$ $`\rho _r\left(i\right)\text{ for }i=N\mathrm{}+1,\mathrm{},N1`$
$`J^{}`$ $`=`$ $`\beta \rho _r\left(N\right).`$
These relations can be used to express $`\rho \left(N\right)=_{i=N\mathrm{}+1}^N\rho _r\left(i\right)`$ in terms of $`J^{}`$, leading to the termination boundary condition in Eqn. (21). The boundary conditions for $`\chi `$ will serve to fix the constant of integration as well as $`J^{}`$, which is still an unknown in Eqn. (19). (The current thus determined is expected to match closely the current predicted by Eqn. 16.) We discuss the various phases separately.

For smaller values of $`J^{}`$ (i.e., $`J^{}<\widehat{J}`$), $`\chi \left(x\right)`$ has fixed points at $`\chi _>`$ and $`\chi _<`$. Boundary conditions determine whether the steady state profile approaches $`\chi _<`$ (low density) or $`\chi _>`$(high density). Each solution corresponds to part of the phase diagram shown in Figure 5. There is a “kink” solution when $`\chi \left(0\right),\chi \left(1\right)(\chi _<,\chi _>)`$ with $`\chi \left(0\right)\chi _<`$ and $`\chi \left(1\right)\chi _>`$ that corresponds to the first order transition line.

In the high density phase ($`\beta <\alpha `$) where $`\chi >\chi _>`$, the density profile is given implicitly by

$$\left[1+\left(\mathrm{}1\right)\chi _>\right]\mathrm{ln}\left(\frac{\chi \chi _>}{\chi _0\chi _>}\right)\left[1+\left(\mathrm{}1\right)\chi _<\right]\mathrm{ln}\left(\frac{\chi \chi _<}{\chi _0\chi _<}\right)=\frac{Rx}{D}$$

(where $`\chi _0\chi \left(0\right)`$ is the EPD at initiation). It may be noted that as $`N\mathrm{}`$, the bulk density and termination boundary density both approach $`\chi _>`$. Setting $`\chi \left(1\right)=\chi _>`$ and eliminating $`J^{}`$ from the termination boundary condition (Eqn. 23) and the definition of $`\chi _>`$ (Eqn. 20) yields a bulk $`\chi `$ of $`\left(1\beta \right)/\left[1+\beta \left(\mathrm{}1\right)\right]`$. We expect the bulk density in the high density phase to match the termination reservoir density. Indeed, the bulk $`\chi `$ value calculated here is consistent with the reservoir density $`\rho _+=1\beta `$ given above (Eqn. 15).

Though it is impossible to write a closed form for $`\chi \left(x\right)`$ in general, a convenient expression is

$$\chi =\chi _>+\left(\chi _0\chi _>\right)\mathrm{exp}\left[\mu _+x\right]\left(\frac{\chi \chi _<}{\chi _0\chi _<}\right)^\gamma $$
(24)
where

$$\mu _+\frac{R}{D\left[1+\left(\mathrm{}1\right)\chi _>\right]}$$

and

$$\gamma =\frac{1+\left(\mathrm{}1\right)\chi _<}{1+\left(\mathrm{}1\right)\chi _>}.$$

Note that $`R`$, $`\chi _<`$, and $`\chi _>`$ can be conveniently estimated using the $`J^{}`$ value given by Eqn. (16). The advantage here is that the unknown factor $`\left(\frac{\chi \chi _<}{\chi _0\chi _<}\right)^\gamma `$ displays only limited variation, since $`\chi `$ lies outside the interval $`[\chi _<,\chi _>]`$. Thus, we expect to find better and better approximations by exploiting an *iterative* scheme. Starting with the first iteration

$$\chi ^{(1)}=\chi _>+\left(\chi _0\chi _>\right)\mathrm{exp}\left[\mu _+x\right],$$

this procedure involves substituting repeatedly for the $`\chi `$ in the unknown factor of Eqn. (24), namely,

$$\chi ^{(k)}=\chi _>+\left(\chi _0\chi _>\right)\mathrm{exp}\left[\mu _+x\right]\left(\frac{\chi ^{(k1)}\chi _<}{\chi _0\chi _<}\right)^\gamma .$$

In Figure 10, we show an example of three such iterations, converging rapidly to the real $`\chi \left(x\right)`$, and in the inset a comparison with the $`\chi \left(x\right)`$ observed in simulations. Near the termination end of the system, the EPD from simulations shows a depletion below the bulk density. This reduction is characteristic of the high density phase and was originally observed in the numerical results of MG . The continuum limit does not capture this feature.

A similar analysis may be carried out for the low density phase when $`\chi <\chi _<`$. In this case, the boundary layer occurs at the termination end, so it is convenient to apply the $`\chi \left(1\right)`$ boundary condition explicitly and use the $`\chi \left(0\right)`$ boundary condition to determine $`J^{}`$. This method also works for high and low density phase profiles when $`\chi \left(0\right),\chi \left(1\right)(\chi _<,\chi _>)`$, as long as $`\alpha `$ and $`\beta `$ are not too similar. The analysis fails close to the first order transition line, as we would expect since Eqn. (19) predicts a “kink” rather than a linear density profile on the transition line.

For the maximal current phase, we know that $`J^{}>\widehat{J}`$, so that $`R`$ is purely imaginary. Since $`\chi ^{}`$ is negative definite, we have a downward sloping profile, which is the generalization of the $`\mathrm{tan}\left(x/\xi \right)`$ (or $`\mathrm{cot}\left(x/\xi \right)`$) profiles in the ordinary driven lattice gas BSZ . In the large $`N`$ limit, a macroscopically large region of the profile will be almost flat (corresponding to $`\chi ^{}0`$), where the density assumes the optimal value $`\widehat{\rho }`$ (or $`\chi \widehat{\chi }`$). Meanwhile, $`J^{}`$ will approach the maximal current $`\widehat{J}`$ from above with $`O\left(1/N^2\right)`$ terms. The details are somewhat involved, so that we will show the solution and discuss its properties only for large $`N`$.

Defining the real quantity $`\stackrel{~}{R}=iR`$, Eqn. (19) can be integrated to show that the steady state density profile is given by

$$\mathrm{ln}\left[\frac{\chi ^2\left(1\stackrel{~}{J}\right)\chi +J^{}}{\chi _0^2\left(1\stackrel{~}{J}\right)\chi _0+J^{}}\right]+\frac{4+2\left(\mathrm{}1\right)\left(1\stackrel{~}{J}\right)}{\stackrel{~}{R}\left(\mathrm{}1\right)}\left(\theta \theta _0\right)=\frac{2x}{D\left(\mathrm{}1\right)},$$
(25)
where

$$\theta \left(x\right)\mathrm{arctan}\frac{\stackrel{~}{R}}{1\stackrel{~}{J}2\chi }$$

and $`\chi _0=\chi \left(0\right)`$ is given by Eqn. (22). Evaluating Eqn. (25) at the termination end ($`x=1,\chi =\chi \left(1\right)`$) and using Eqn. (23), we arrive at an equation for determining $`J^{}`$ in terms of the control parameters $`\alpha ,\beta `$. Needless to say, this equation is too complex to solve analytically. Nevertheless, we can gain some insight by considering the large $`N`$ limit. Defining $`ϵ1/N`$ for convenience, recall that $`D`$ is $`O\left(ϵ\right)`$. Next, let us assume that $`J^{}=\widehat{J}+O\left(ϵ^2\right)`$ and show that it is justified later. This leads us to $`\stackrel{~}{R}`$ being $`O\left(ϵ\right)`$ and $`\left(1\stackrel{~}{J}\right)/2=\widehat{\chi }+O\left(ϵ^2\right)`$, so that $`\chi _{>,<}=\widehat{\chi }+O\left(ϵ^2\right)\pm iO\left(ϵ\right)`$. Thus, the quadratic form in the argument of the $`\mathrm{ln}`$ in Eqn. (25) never becomes smaller than $`\stackrel{~}{R}^2/4=O\left(ϵ^2\right)`$, so that this term never exceeds $`O\left(\mathrm{ln}N\right)`$. Since the other terms are generally $`O\left(N\right)`$ (from $`D`$ and $`\stackrel{~}{R}`$), we can write an approximate equation by neglecting the $`\mathrm{ln}`$ term, leaving an expression for $`\theta ,`$ and taking the co-tangent of both sides:

$$\chi \widehat{\chi }=\frac{1}{2}\stackrel{~}{R}\mathrm{cot}[x/\xi +const.]+O\left(ϵ^2\right),$$
(26)
where $`\xi 2D\sqrt{\mathrm{}}/\stackrel{~}{R}`$ is of $`O\left(1\right)`$. Note that the constant here can be written as $`\mathrm{arccot}\left[2\left(\chi _0\widehat{\chi }\right)/\stackrel{~}{R}\right]`$, which vanishes as $`\stackrel{~}{R}/\left(\chi _0\widehat{\chi }\right)`$ in the limit $`ϵ0`$. Thus, for typical values of $`x`$, we see that $`\chi =\widehat{\chi }+O\left(ϵ\right)`$. For the termination end, $`J^{}`$ must be carefully fixed so that the argument of the co-tangent in Eqn. (26) approaches $`\pi `$ in an appropriate way for the right hand side to equal $`\chi \left(1\right)\widehat{\chi }`$. This requirement is consistent with the original assumption that $`J^{}=\widehat{J}+O\left(ϵ^2\right)`$.

Returning to the finite $`N`$ case, although $`\chi \left(x\right)`$ cannot be determined explicitly, we find that an iterative method is again helpful and converges rapidly. This time, successive iterations of $`\chi `$ are substituted into the term $`\mathrm{ln}\left[\frac{\chi ^2\left(1\stackrel{~}{J}\right)\chi +J^{}}{\chi _0^2\left(1\stackrel{~}{J}\right)\chi _0+J^{}}\right]`$ of Eqn. (25). As Figure 11 shows, there is good agreement between the steady state profiles from simulation data and this mean field theory.

## IV Disordered hopping rates in the open system

We finally turn to the difficult problem of the TASEP with extended objects and quenched disorder. At each codon of the mRNA, a ribosome translating the mRNA must wait for an appropriate aa-tRNA to decode that codon. aa-tRNA availability ranges over about an order of magnitude Solom , and it is thought that the elongation (hopping) rates may also range over an order of magnitude, with the slowest ones comparable to the initiation and termination rates Bergmann .

Since previous studies of the TASEP with quenched disorder involve $`\mathrm{}=1`$, their relevance to the process of translation may be questioned. Indeed, we believe that the system will display essential differences when these point particles are replaced by extended objects. First, as we have seen in the uniform elongation studies above, generalizing to the $`\mathrm{}>1`$ case led to significant changes. Now, with non-uniform rates, steric hinderance should play an even larger role, since extended particles can block multiple sites with potentially different rates. Thus, we will devote this section to a limited study of the effects of $`\mathrm{}>1`$ on open systems with non-uniform elongation rates: (a) systems with a single slow site, (b) bounds on the current in general open systems, and (c) as an illustration of systems with full disorder, simulation results for a real mRNA sequence.

### IV.1 Simple model of a single internal blockage

We begin with the simplest form of quenched disorder, namely, having a single *internal* site $`i`$ with a reduced elongation rate: $`r<1`$. If this blockage were moved to the termination site, then $`r`$ would carry the label $`\beta `$ to conform with the notation above. In most of the simulation studies for this section, we focus on $`\alpha =\beta =1`$ with the blockage at the center of the mRNA (i.e., $`i=100`$ with $`N=200`$). Though the naïve expectation is that the blockage should limit the current in the same way, regardless of its location, we observe that the current is noticeably reduced if the slow site is internal. This reduction can be understood as follows. Once a particle moves past the slow site, the particle behind it is not necessarily free to move to the slow site. Instead, because of the random sequential updating, the particle behind has a nonzero chance of being blocked by the particle ahead. With just two particles in the system, this effect can be accounted for exactly. Beyond the scope of this paper, the details will be published elsewhere. Here, let us present some simple heuristic arguments which lead us to results that agree well with simulation data.

Extrapolating from the current in the uniform system (Eqn. 16), we speculate that

$$J=\rho _{source}\frac{P_{free}}{T},$$
(27)
where $`\rho _{source}`$ is the average coverage density behind the slow site and $`P_{free}`$ is the approximate probability for a particle just beyond the slow site to advance. Lastly, $`T`$ is the average time to travel through the slow site and the $`\mathrm{}1`$ sites that precede it. For example, in the termination-limited case, $`\rho _{source}=1\beta `$, $`T=\mathrm{}1+1/\beta `$, and $`P_{free}=1`$ (no steric hindrance beyond the termination site), giving $`J=\beta \left(1\beta \right)/\left[1+\beta \left(\mathrm{}1\right)\right]`$. We estimate $`\rho _{source}`$ from the bulk density induced if the lattice were truncated immediately after the slow site. $`P_{free}`$ can be estimated from the density-dependent head spacing in a closed system (Eqn. 8), using the bulk density that would be induced if the lattice began with the slow site.

To determine $`T`$ approximately, we consider the behavior of two particles in an infinite system with a single slow site (located at the origin). Denoting the position of the leading and following particles by $`\zeta `$ and $`\eta `$, respectively, the evolution of this system can be regarded as a random walk in the $`\zeta `$-$`\eta `$ plane, confined to $`\zeta \eta \mathrm{}`$. During each time step, the walker moves, with equal probability, either “upwards” ($`\eta \eta +1`$) or “to the right” ($`\zeta \zeta +1`$). When it arrives at the $`\zeta =\eta +\mathrm{}`$ line, the walker remains stationary with probability $`\frac{1}{2}`$ (and moves “to the right” otherwise). Assuming that there are $`m+1`$ holes between the particles immediately after the leading one leaves the slow site, so that the walker is “initially” located at $`(\zeta ,\eta )=(1,\mathrm{}m)`$, we computed numerically $`\tau _m`$, the average time for the walker to arrive at the $`\eta =0`$ line (i.e., for the second particle to reach the slow site). Note that these $`\tau `$’s take into account the steric hindrance due to the leading particle, which we assume is always free to move. To extract $`T`$, we make the assumption that the left particle advances its first $`m+1`$ steps in time $`m+1`$, leaving $`\mathrm{}1`$ more steps to reach the slow site. Now, the probability $`p\left(m\right)`$ for finding a gap of $`m`$ holes before the leading particle passes the slow site can be estimated by inserting the bulk density before the slow site ($`\rho _{source}`$) in Eqn. (8). Taking the average over these initial starting positions and accounting for the time to move over the slow site ($`1/r`$), we obtain $`T=_m\left[\tau _m\left(m+1\right)\right]p\left(m\right)+1/r`$. Predictions for the current from Eqn. (27), compared with the actual current from Monte Carlo simulations, are shown in Table 1 for several values of the slow rate $`r`$. Despite the approximations involved in this approach, the agreement is surprisingly good.

### IV.2 Bounds for the current

To model a real mRNA, we must allow for arbitrary translation rates associated with each codon. Let us denote the rate at codon $`i`$ by $`k_i`$. Due to the excluded volume constraint, however, it is meaningful to consider also $`K_{\mathrm{},i}`$, the *maximum* rate for a ribosome to translate a stretch of $`\mathrm{}`$ sites beginning with site $`i`$:

$$K_{\mathrm{},i}\left(\underset{q=i}{\overset{i+\mathrm{}1}{}}\frac{1}{k_q}\right)^1.$$

Now, consider a “window” of any stretch of $`\mathrm{}`$ consecutive sites in the lattice. If one particle is moving within this window, the following particle must wait until the first one passes through entirely before it can begin translating these $`\mathrm{}`$ sites. Thus the characteristic time associated with the current must be at least the time required to translate the slowest stretch of $`\mathrm{}`$ codons. In this way, we find an upper bound to the current, i.e.,

$$J\underset{i\{1,\mathrm{},N\mathrm{}+1\}}{\mathrm{min}}K_{\mathrm{},i}.$$
(28)
One might imagine the slowest segment of $`\mathrm{}`$ codons acting as a “gatekeeper” and preventing the current from exceeding the value in Eqn. (28).

To arrive at a lower bound, we need only replace the elongation rates at each site by the slowest elongation rate. From Eqn. (16) above, we have the current of a system with uniform rate unity. Thus, the minimum current for the disordered system is simply

$$J\left(\underset{i\{1,\mathrm{},N1\}}{\mathrm{min}}k_i\right)\overline{J}(\alpha ,\beta ),$$
(29)
where $`\alpha `$ is the ratio of the initiation rate to the slowest elongation rate, and similarly for $`\beta `$. Though the gap between these bounds for a real system may be too large to be of significant predictive value, they can provide some guide to our understanding of the current.

### IV.3 Simulation of a real gene sequence

To illustrate the full problem of disorder, we have simulated translation of several real mRNA sequences from Escherichia coli strain MG1655, obtained from MG1655 . Elongation rates at each codon were estimated using commonly accepted values for the availability of tRNA in E. coli Solom . The rate at each codon was assumed proportional (with an arbitrary proportionality constant) to the availability of the tRNA decoding that codon, as in Lesnik . Corresponding data were not available for estimating initiation and termination rates, so a range of initiation and termination rates was studied. We assumed that ribosomes cover $`\mathrm{}=12`$ codons Kang ; Heinrich .

Figure 12 shows the steady state coverage density profile for the reasonably well studied gene ompA when the elongation rates are limiting. The lower curve in the figure shows the maximum translation rate $`K_{12}`$ for each window of $`\mathrm{}=12`$ codons, calculated from

$$K_{12}=\left(\underset{q=i}{\overset{i+11}{}}\frac{1}{k_q}\right)^1$$

for each site $`i`$. The minimum of the $`K_{12}`$ values is 4.51 in arbitrary units, thus giving an upper bound for the current (Eqn. 28). This value is significantly higher than the actual current of 3.52. (For comparison, Eqn. 29 gives a lower bound of 0.99 for the current in this system.) It should be noted that the minimum $`K_{12}`$ occurs at codon 71, which is approximately the location behind which the ribosome density is very high, due to ribosomes “piled up” behind the slow region. In general, lower values for the rate $`K_{12}`$ correspond to higher ribosome densities, and higher $`K_{12}`$ to lower ribosome densities, leading to an approximate symmetry between $`K_{12}`$ and $`\rho `$. Thus the $`K_{12}`$ values are useful in understanding the ribosome density profiles observed.

## V Conclusions

This work generalizes the well-studied $`\mathrm{}=1`$ TASEP model to particles with extended sizes. Since there is a difference between particle density $`\rho _r`$ and density of occupied sites $`\rho `$, the familiar particle-hole symmetry ($`\rho 1\rho `$) takes the form $`\rho _r1\rho `$ here. Exact results for TASEP on a uniform ring, including probability distributions for the current and for particle headway, were found. Particularly useful in the latter part of this study is the new current-density relationship, which is interpreted as a novel density dependent mobility factor (or “diffusion constant”).

An extremal principle PS  based on domain wall theory allowed the closed system current-density relation to predict currents and bulk densities in the uniform open system as functions of the initiation and termination rates $`\alpha `$ and $`\beta `$. The phase diagram for the open system was thus determined using the extremal principle. Domain wall theory also provided an explanation for the linear density profiles and other unique characteristics observed at the first order transition between the high and low density phases. Given the ability of this theory to describe particles with length $`\mathrm{}>1`$, it might be exploited further to determine fluctuations in number of bound particles in the steady state and behavior in the pre-steady state regime, as has been done for the $`\mathrm{}=1`$ system SA .

Based on the new mobility factor, a simple continuum limit led to a differential equation for the density profile in open systems. Though non-linear in general, this equation can be transformed, in the case of a symmetric exclusion process, to the familiar *linear* diffusion equation for an effective particle density $`\chi `$. As a result, the stationary density profile, expressed in terms of $`\chi `$, is again linear. For the asymmetric exclusion process, the differential equation is more complex, due to an extra term to account for the drive. Analytic expressions, albeit implicit, for the stationary profiles were obtained and solved numerically through a rapidly converging, iterative procedure. The results matched simulation data closely. However, the predictive power of this approach is limited, since the “coarse-grained” parameter $`D`$ cannot be derived from the microscopic jump rates. Naïve estimates for it are of the right order of magnitude but quantitatively inaccurate, reflecting the importance of particle correlations. Regarding $`D`$ as a phenomenological quantity, we simply fit it to data. At a more detailed level, open questions about density profiles remain. In particular, our continuum limit cannot explain the particle depletion and period-$`\mathrm{}`$ structure near the termination end of high density systems. Profiles of the particle density $`\rho _r`$ are also of interest, though beyond the scope of this simple continuum theory. Absent from the $`\mathrm{}=1`$ TASEP, peaks with spacing $`\mathrm{}`$ extend far behind a blockage, with a decay length well beyond typical microscopic scales (data not shown). Evidently, when extended objects are included, even a simple uniform TASEP tantalizes us with a rich variety of novel behavior.

Finally, systems with quenched disorder in the particle hopping rates were briefly considered. Effects of a single internal slow site on the current were estimated with fair accuracy by considering the average delay a particle near the slow site experiences due to a particle ahead of it. Also, bounds on the current were determined for general disordered systems. A real gene sequence was simulated, leading to a complicated density profile. The parameter $`K_{\mathrm{},i}`$, the maximum rate to translate a stretch of $`\mathrm{}`$ sites beginning with site $`i`$, proved helpful in understanding the shape of the density profile, but the disordered system remains far from solved.

We close with some speculation about the relevance of this work to an understanding of translation. The advent of functional genomics technologies to measure simultaneously mRNA and protein expression profiles from many thousands of genes provides special opportunities to begin to understand gene expression regulation. Our results for currents (i.e., protein production rates) and ribosome densities for uniform systems cannot be directly compared with such experimental data from typical bacterial cells because translation is not approximated well enough by a uniform system. However, it is possible to use the uniform system results to interpret data from an mRNA artificially constructed to be uniform. Although there are no reports in the literature of systems that are approximately uniform, the use of an in vitro translation system Cen  provides an opportunity to make such experimental observations.

It is known that the relationship between mRNA and protein levels in typical cells is nonlinear Gygi ; Ideker . Specifically, when cells are grown under two different sets of conditions, the amount of protein corresponding to a particular gene may be down-regulated while its corresponding mRNA is up-regulated, and the opposite may be true for other genes measured from the same samples. We expect that in biological systems, the initiation rate $`\alpha `$ should be an increasing function of the availability of ribosomes within the cell. The protein production rate (current $`J`$, Eqn. 16) is, in turn, a non-decreasing function of $`\alpha `$. This analysis thus suggests that the observed nonlinear relationship can arise from changes in the availability of ribosomes given the nonlinear relationship between $`J`$ and $`\alpha `$. However, this situation would cause all protein production rates to change in the same direction, that in which the ribosome availability changes. It would not permit some proteins to be up-regulated while others are down-regulated—which is the situation observed experimentally. Further nonlinearity would arise if mRNA’s were to compete for available aa-tRNA as well as for ribosomes. Detailed modeling of this effect will require a better understanding of systems with quenched disorder, in which the elongation rates result from aa-tRNA availability.

###### Acknowledgements.

We thank Johannes Hager, James Sethna, Beate Schmittmann, Vassily Hatzimanikatis, Amit Mehra, Tom Chou, and Anatoly Kolomeisky for their helpful suggestions. RKPZ acknowledges support from the National Science Foundation through grants DMR-9727574 and 0088451. KHL acknowledges support from the NSF through grants BES-0120315 and 9874938. LBS was supported by an NSF Graduate Research Fellowship and a Corning Foundation Fellowship. This research was conducted using the resources of the Cornell Theory Center, which receives funding from Cornell University, New York State, federal agencies, foundations, and corporate partners.

*

## Appendix A Derivation of modified diffusion equation from discrete version

A modified diffusion equation qualitatively equivalent to Eqn. (19) can be obtained by taking the naïve continuum limit of the discrete mean field equations of MacDonald et al. We begin with Eqn. (7) of MGP  for the current from lattice site $`j`$ to $`j+1`$:

$$q_j=\frac{n_j^{\left(L\right)}n_{j+1}^{\left(0\right)}}{n_{j+1}^{\left(0\right)}+n_{j+L}^{\left(L\right)}}$$

and make the following correspondences with our continuum notation:

$`n_j^{(0)}`$ $``$ $`\rho _h\left(x\right)`$
$`n_j^{(L)}`$ $``$ $`\rho _r\left(x\right)`$
$`L`$ $``$ $`\mathrm{}`$
$`q`$ $``$ $`J`$
where $`x`$ lies in $`[0,N]`$ here. Thus we have

$$J\left(jj+1\right)=\frac{\rho _r\left(x\right)\rho _h\left(x+1\right)}{\rho _h\left(x+1\right)+\rho _r\left(x+\mathrm{}\right)}$$
(30)
Performing a series expansion of Eqn. (30) and keeping terms up to second order, we find that

$`J\left(jj+1\right)`$ $`=`$ $`{\displaystyle \frac{\rho _r\left(x\right)\left[\rho _h\left(x\right)+\left(1\right)\frac{\rho _h\left(x\right)}{x}+\frac{1}{2!}\frac{^2\rho _h\left(x\right)}{x^2}+\mathrm{}\right]}{\left[\rho _h\left(x\right)+\left(1\right)\frac{\rho _h\left(x\right)}{x}+\frac{1}{2!}\frac{^2\rho _h\left(x\right)}{x^2}+\mathrm{}\right]+\left[\rho _r\left(x\right)+\mathrm{}\frac{\rho _r\left(x\right)}{x}+\frac{\mathrm{}^2}{2!}\frac{^2\rho _r\left(x\right)}{x^2}+\mathrm{}\right]}}`$
$`=`$ $`{\displaystyle \frac{\rho _r\rho _h}{\rho _s}}+{\displaystyle \frac{\rho _r\rho _h^{}}{\rho _s}}+{\displaystyle \frac{1}{2}}{\displaystyle \frac{\rho _r}{\rho _s^2}}\left[\rho _r\rho _h^{\prime \prime }\mathrm{}^2\rho _h\rho _r^{\prime \prime }\right]`$
by using $`\rho _h\left(x\right)=1\mathrm{}\rho _r\left(x\right)`$. Similarly,

$$J\left(j1j\right)=\frac{\rho _r\rho _h}{\rho _s}\frac{\rho _r^{}\rho _h}{\rho _s}+\frac{\rho _r\rho _h}{\rho _s^2}\rho _s^{}+\frac{\rho _r\rho _h}{2\rho _s}\left[2\left(\frac{\rho _s^{}}{\rho _s}\right)^22\frac{\rho _r^{}}{\rho _r}\frac{\rho _s^{}}{\rho _s}+\frac{\rho _r^{\prime \prime }}{\rho _r}\left(\mathrm{}1\right)^2\frac{\rho _r^{\prime \prime }}{\rho _s}\right]$$

It can then be shown that

$`{\displaystyle \frac{J}{x}}`$ $`=`$ $`J\left(jj+1\right)J\left(j1j\right)`$
$`=`$ $`{\displaystyle \frac{}{x}}\left[{\displaystyle \frac{1}{2\rho _s^2}}\left[\mathrm{}\left(1\mathrm{}\right)\rho _r^2+1\right]{\displaystyle \frac{\rho _r}{x}}\right]+{\displaystyle \frac{}{x}}\left({\displaystyle \frac{\rho _r\rho _h}{\rho _s}}\right).`$
This leads to the modified diffusion equation

$$J=D\left(\frac{1+\mathrm{}\left(1\mathrm{}\right)\rho _r^2}{\rho _s^2}\right)\left(\frac{\rho _r}{x}\right)+E\frac{\rho _r\rho _h}{\rho _s},$$
(31)
where we have defined $`D=1/2`$ and $`E=1`$.

Finally, we rewrite Eqn. (31) in terms of the effective particle density $`\chi `$. Thus we find that the mean field equations of MGP  predict

$$\frac{\chi }{x}=\frac{1+\left(\mathrm{}1\right)\chi }{1+\left(\mathrm{}1\right)\chi \left(\mathrm{}1\right)\chi ^2}\left[\stackrel{~}{J}\chi +J\chi \left(1\chi \right)\right].$$
(32)
We expect Eqn. (32) to be comparable to Eqn. (19), the steady state density profile equation derived previously from simple arguments. Indeed, both equations give the same fixed points for $`\chi `$, thus producing qualitatively identical density profiles. Further, their quantitative differences have little effect on the shape of the density profile (data not shown).