A Deterministic Mathematical Model for Bidirectional Excluded Flow with Langmuir KineticsThe research of MM and TT is partially supported by research grants from the Israeli Ministry of Science, Technology & Space and the Binational Science Foundation. The research of MM is also supported by a research grant form the Israeli Science Foundation.

A Deterministic Mathematical Model for Bidirectional Excluded Flow with Langmuir Kineticsthanks: The research of MM and TT is partially supported by research grants from the Israeli Ministry of Science, Technology & Space and the Binational Science Foundation. The research of MM is also supported by a research grant form the Israeli Science Foundation.

Yoram Zarai, Michael Margaliot, and Tamir Tuller Y. Zarai is with the School of Electrical Engineering, Tel-Aviv University, Tel-Aviv 69978, Israel. E-mail: yoramzar@mail.tau.ac.il M. Margaliot is with the School of Electrical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978, Israel. E-mail: michaelm@eng.tau.ac.il T. Tuller is with the department of Biomedical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978, Israel. E-mail: tamirtul@post.tau.ac.il
Abstract

In many important cellular processes, including mRNA translation, gene transcription, phosphotransfer, and intracellular transport, biological “particles” move along some kind of “tracks”. The motion of these particles can be modeled as a one-dimensional movement along an ordered sequence of sites. The biological particles (e.g., ribosomes, RNAPs, phosphate groups, motor proteins) have volume and cannot surpass one another. In some cases, there is a preferred direction of movement along the track, but in general the movement may be two-directional, and furthermore the particles may attach or detach from various regions along the tracks (e.g. ribosomes may drop off the mRNA molecule before reaching a stop codon).

We derive a new deterministic mathematical model for such transport phenomena that may be interpreted as the dynamic mean-field approximation of an important model from mechanical statistics called the asymmetric simple exclusion process (ASEP) with Langmuir kinetics. Using tools from the theory of monotone dynamical systems and contraction theory we show that the model admits a unique equilibrium, and that every solution converges to this equilibrium. This means that the occupancy in all the sites along the lattice converges to a steady-state value that depends on the parameters but not on the initial conditions. Furthermore, we show that the model entrains (or phase locks) to periodic excitations in any of its forward, backward, attachment, or detachment rates.

We demonstrate an application of this phenomenological transport model for analyzing the effect of ribosome drop off in mRNA translation. One may perhaps expect that drop off from a jammed site may increase the total flow by reducing congestion. Our results show that this is not true. Drop off has a substantial effect on the flow, yet always leads to a reduction in the steady-state protein production rate.

Monotone dynamical systems, systems biology, synthetic biology, mRNA translation, gene transcription, ribosome flow model, ribosome drop off, Langmuir kinetics, bi-directional flow, intracellular transport, contraction theory, contraction after a short transient, entrainment.

I Introduction

Movement is essential for the functioning of cells. Cargoes like organelles and vesicles must be carried between different locations in the cells. The information encoded in DNA and mRNA molecules must be decoded by “biological machines” (RNA polymerases and ribosomes) that move along these molecules in a sequential order.

Many of these important biological transport processes are modeled as the movement of particles along an ordered chain of sites. In the context of intercellular transport, the particles are motor proteins and the chain models actin filaments or microtubules. In transcription, the particles are RNAPs moving along the DNA molecule, and in translation the particles are ribosomes moving along the mRNA molecule (see Figure 1).

Fig. 1: Biological processes that can be studied using the model suggested in this paper.

The movement in such processes may be unidirectional, as in mRNA translation elongation, or bidirectional, as in transcription or translation initiation. Indeed, the normal forward flow of the RNAP may be interrupted, due to transcription errors and various obstacles such as nucleosomes, in which case the RNAP tracks back a few nucleotides and then resumes its normal forward flow [54, 41, 8, 10]. Translation initiation in eukaryotes usually includes diffusion from the 5’end of the transcript towards the start codon [1]. This diffusion process is believed to be bi-directional, but with a preference to the 5’3’ direction. The movement of motor proteins like kinesin and dynein along microtubules is typically unidirectional, but can be two-directional as well [1].

To increase efficiency, many particles may move simultaneously along the same track thus pipelining the production process. For example, to increase translation efficiency, a number of ribosomes may act simultaneously as polymerases on the same mRNA molecule [66, 4].

The moving biological particles have volume and usually cannot overtake a particle in front of them. This means that a slowly moving particle may lead to the formation of a traffic jam behind it. For example, Leduc et al. [31] have studied Kip3, a yeast kinesin-8 family motor, and demonstrated that motor protein traffic jams can exist, given the right conditions. Other studies have suggested that traffic jams of RNAP [ribosomes] may evolve during transcription [translation] [4, 27, 9].

In some of these biological transport processes, the biological machines may either attach or detach at various sites along the tracks. For example, ribosomes may detach from the mRNA molecule before reaching the stop codon due to traffic “jams” and ribosome-ribosome interactions or due to depletion in the concentration of tRNAs [72, 28, 57]. Also, it is known that kinesin-family motor proteins are more susceptible to dissociation when their pathway is blocked [14, 62].

Defects in these transport processes may lead to severe diseases or may even be lethal. For example, [53] lists the implications of malfunctions of protein motors in disease and developmental defects.

Developing a better understanding of these dynamical biological processes by combining mathematical modeling and biological experiments will have far reaching implications to basic science in fields such as molecular evolution and functional genomics, as well as applications in synthetic biology, biotechnology, human health, and more. Mathematical or computational modeling is especially important in developing approaches for manipulating and controlling these processes, e.g. in order to optimize various goals in biotechnology.

A standard model for such transport processes is the asymmetric simple exclusion process (ASEP) [55, 73]. This is a stochastic model describing particles that hop along an ordered lattice of sites. Each site can be either empty or occupied by a single particle, and a particle can only hop to an empty site. This “simple exclusion principle” represents the fact that the particles have volume and cannot overtake one another. Simple exclusion generates an indirect coupling between the particles. In particular, traffic jams may develop behind a slow-moving particle.

In ASEP, a particle may hop to any of the two neighboring sites (but only if they are free). Typically, a particle can attach the lattice in one of its ends and detach from the other end. When particles can also attach or detach at internal sites along the lattice, the model is referred to as ASEP with Langmuir kinetics. In the special case where the hops are unidirectional, ASEP is sometimes referred to as the totally asymmetric simple exclusion process (TASEP). A TASEP-like system with Langmuir kinetics has been used to model limit order markets in [65], and is often used in modeling molecular motor traffic [42, 43, 32, 33, 16]. More generally, ASEP has become a fundamental model in non-equilibrium statistical mechanics, and has been applied to model numerous natural and artificial processes including traffic and pedestrian flow, the movement of ants, evacuation dynamics, and more [52].

In this paper, we introduce a deterministic mathematical model that may be interpreted as the dynamic mean-field approximation of ASEP with Langmuir kinetics (MFALK). We analyze the MFALK using tools from systems and control theory. In particular, we apply some recent developments in contraction theory to prove that the model is globally asymptotically stable, and that it entrains to periodic excitations in the transition/attachment/detachment rates. In other words, if these rates change periodically in time with some common period  then all the state-variables in the MFALK converge to a periodic solution with period . This is important because many biological processes are excited by periodic signals (e.g. the 24h solar day or the periodic cell-division process), and proper functioning requires phase-locking or entrainment to these excitations.

Our work is motivated by the analysis of a model for mRNA translation called the ribosome flow model (RFM) [48]. This is the mean-field approximation of the unidirectional TASEP without Langmuir kinetics (see, e.g., [52, section 4.9.7] and [6, p. R345]). Recently, the RFM has been studied extensively using tools from systems and control theory [36, 71, 37, 38, 35, 44, 45, 47, 70]. The analysis is motivated by implications to many important biological questions. For example, the sensitivity of the protein production rate to the initiation and elongation rates along the mRNA molecule [45], maximization of protein production rate [44], the effect of ribosome recycling [38, 47], and the consequences of competition for ribosomes on large-scale simultaneous mRNA translation in the cell [46] (see also [19, 2] for some related models).

The MFALK presented here is much more general than the RFM, and can thus be used to model and analyze many transport phenomena, including all the biological processes mentioned above, that cannot be captured using the RFM. We demonstrate this by using the MFALK to model and analyze mRNA translation with ribosome drop off - a feature that cannot be modeled using the RFM.

Ribosome drop off is a fundamental phenomena that has received considerable attention (see, e.g., [57, 25, 24, 68, 7, 61, 18, 28, 22, 20]). In many cases, ribosome drop off is deleterious to the cell since translation is the most energetically consuming process in the cell and, furthermore, drop off yields truncated, non-functional proteins. Thus, transcripts undergo selection to minimize drop off or its energetic cost [67, 63, 61, 18, 28, 20]. There are various hypotheses on the biological advantages of ribosome drop off. For example, Zaher and Green [69] have suggested that ribosome drop off is related to proof reading. One may perhaps expect that another advantage is that drop off from a jammed site may increase the total flow by reducing congestion. Our results using analysis of the MFALK show that this is not true. Drop off has a substantial effect on the flow, yet it always leads to a reduction in the steady-state protein production rate.

The remainder of this paper is organized as follows. The next section describes the new mathematical model. Section III presents our main analysis results. Section IV describes the application of the MFALK to model mRNA translation with ribosome drop off. The final section concludes and describes possible directions for further research. To streamline the presentation, all the proofs are placed in the Appendix.

Ii The model

The MFALK is a set of  first-order nonlinear differential equations, where denotes the number of compartments or sites along the “track”. Each site is associated with a state variable  describing the normalized “level of occupancy” at site  at time , with  [] representing that site  is completely free [full] at time . Since  for all , it may also be interpreted as the probability that site  is occupied at time .

The MFALK contains four sets of non-negative parameters:

  • , , controls the forward transition rate from site to site ,

  • , , controls the backward transition rate from site to site ,

  • , , controls the attachment rate to site ,

  • , , controls the detachment rate from site ,

where we arbitrarily refer to left-to-right flow along the chain as forward flow, and to flow in the other direction as backward flow.

The dynamical equations describing the MFALK are:

(1)

To explain these equations, consider for example the equation for the change in the occupancy in site , namely,

The term  represents the flow from site  to site . This increases with the occupancy in site , and decreases with the occupancy in site . In particular, this term becomes zero when , i.e. when site  is completely full. This is a “soft” version of the hard exclusion principle in ASEP: the effective entry rate into a site decreases as it becomes fuller. Note that the constant  describes the maximal possible transition rate from site  to site . Similarly, the term  represents the flow from site  to site . The term  [] represents the backward flow from site  to site  [site  to site ]. Note that these terms also model soft exclusion. The term  represents attachment of particles from the environment to site , whereas  represents detachment of particles from site  to the environment (see Fig. 2).

The MFALK is a compartmental model [21, 51], as every state-variable describes the occupancy in a compartment (e.g., a site along the the mRNA, gene, microtubule), and the dynamical equations describe the flow between these compartments and the environment. Compartmental models play an important role in pharmacokinetics, enzyme kinetics, basic nutritional processes, cellular growth, and pathological processes, such as tumourigenesis and atherosclerosis (see, e.g., [21, 17] and the references therein). More specifically, the MFALK is a nonlinear tridiagonal compartmental model, as every  directly depends on , and  only.

Note also that

(2)

The term on the right-hand side of the first [second] line here represents the change in  [] due to the flow between the environment and site  [site ], whereas the term on the third line represents the flow between internal sites and the environment.

Output rate
Fig. 2: Topology of the MFALK.

The output rate from site  at time  is the total flow from this site to the environment:

(3)

Note that  may be positive, zero, or negative.

In the particular case where  for all  the MFALK becomes the RFM, i.e. the dynamic mean-field approximation of the unidirectional TASEP with open boundary conditions and without Langmuir kinetics.

Let  denote the solution of (II) at time  for the initial condition . Since the state-variables correspond to normalized occupancy levels, we always assume that  belongs to the closed -dimensional unit cube:

Let  denote the interior of , and let denote the boundary of . The next section analyzes the MFALK defined in (II).

Iii Main Results

Iii-a Invariance and persistence

It is straightforward to show that  is an invariant set for the dynamics of the MFALK, that is, if  then  for all . The following result shows that a stronger property holds. Recall that all the proofs are placed in the Appendix. For notational convenience, let , , and .

Proposition 1

Suppose that at least one of the following two conditions holds:

(4)

or

(5)

Then for any  there exists  such that

(6)

for all , all , and all .

This means in particular that trajectories that emanate from the boundary of  “immediately” enter . This result is useful because as we will see below on the boundary of  the MFALK looses some important properties. For example, the Jacobian matrix of the dynamics (II) is irreducible on , but becomes reducible on some points on the boundary of .

Iii-B Contraction

Differential analysis and in particular contraction theory proved to be a powerful tool for analyzing nonlinear dynamical systems. In a contractive system, trajectories that emanate from different initial conditions contract to each other at an exponential rate [34, 49, 3]. Let  denote the  norm, i.e. for .

Proposition 2

Let

Note that . For any  and any ,

(7)

This implies that the  distance between any two trajectories contracts with the exponential rate . Roughly speaking, this also means that increasing all the sums , , makes the system “more contractive”. Indeed, these parameters have a direct stabilizing effect on the dynamics of site , whereas the other parameters affect the site indirectly via the coupling to the two adjacent sites.

When , (7) only implies that the  distance between trajectories does not increase. This property is not strong enough to prove the asymptotic properties described in the subsections below. Indeed, in this case it is possible that the MFALK will not be contractive with respect to any fixed norm. Fortunately, a certain generalization of contraction turns out to hold in this case.

Consider the time-varying dynamical system

(8)

whose trajectories evolve on a compact and convex set . Let  denote the solution of (8) at time  for the initial condition . System (8) is said to be contractive after a small overshoot (SO) [39] on  w.r.t. a norm  if for any  there exists  such that

for all and all . Intuitively speaking, this means contraction with an exponential rate, but with an arbitrarily small overshoot of .

Proposition 3

Suppose that

(9)

and that at least one of the two conditions (4), (5) holds. Then the MFALK is SO on  w.r.t. the  norm, that is, for any  there exists  such that

(10)

for all and all .

Note that if  for some , that is , then the MFALK decouples into two separate MFALKs: one containing sites , and the other containing sites . Thus, assuming (9) incurs no loss of generality.

There is an important difference between Propositions 2 and 3. If  then Proposition 2 provides an explicit exponential contraction rate. If  then Proposition 3 can be used to deduce SO, but in this result the contraction rate  depends on  and is not given explicitly.

The contraction results above imply that the MFALK satisfies several important asymptotic properties. These are described in the following subsections.

Iii-C Global asymptotic stability

Since the compact and convex set  is an invariant set of the dynamics, it contains an equilibrium point . By Proposition 1, . Applying (10) with  yields the following result.

Corollary 1

Suppose that the conditions in Proposition 3 hold. Then the MFALK admits a unique equilibrium point that is globally asymptotically stable, i.e. , for all .

This means that the rates determine a unique distribution profile along the lattice, and that all trajectories emanating from different initial conditions in  asymptotically converge to this distribution. In addition, perturbations in the occupancy levels along the sites will not change this asymptotic behavior of the dynamics. This also means that various numerical solvers of ODEs will work well for the MFALK (see e.g. [13]).

Example 1

Fig. 3 depicts the trajectories of a MFALK with , , , , , , , , , , , , , for six initial conditions in . It may be seen that all trajectories converge to an equilibrium point 

Fig. 3: Trajectories of the MFALK in Example 1 for six initial conditions in .

The MFALK (II) can be written as

(11)

where

(12)

At steady-state, i.e. for , the left-hand side of all the equations in (11) is zero, so

(13)

Let  denote the parameters of the MFALK. It follows from (13) that if we multiply all these parameters by  then  will not change, that is, . Let

(14)

denote the steady-state output rate. Then , for all , that is, the steady-state production rate is homogeneous of order one w.r.t. the parameters. By (13),

(15)

This yields the following set of recursive equations relating the steady-state occupancy levels and the output rate in the MFALK:

(16)
and also

For a given , this is a set of  equations in the  unknowns: .

Example 2

Consider the MFALK with dimension . Then (III-C) becomes

(17)
and also

This yields the polynomial equation , where

with and .

Note that the polynomial equation admits several solutions , but only one solution corresponds to the unique equilibrium point . For example, for , , , and  for all  the polynomial equation becomes . This admits two solutions  and , with . Substituting  in (2) yields , with , so this is not a feasible solution. Substituting  in (2) yields (all numbers are to four digit accuracy) , which is the unique feasible solution. Thus, the steady-state output rate is 

In general, (III-C) can be transformed into a polynomial equation for . The next result shows that the degree of this polynomial equation grows quickly with .

Proposition 4

Consider the MFALK with dimension  and with , , , for all . Then generically Eq. (III-C) may be written as , where is a polynomial equation in  of degree , and with coefficients that are algebraic functions of the rates.

We note that this is exponential increase in the degree of the polynomial equation is a feature of the MFALK that does not take place in the RFM. Indeed, in the RFM the degree of the polynomial equation for the steady-state production rate grows linearly with .

Let  denote the sign function, i.e.

An interesting question is what is . Indeed, if this is positive (negative) then this means that there is a net steady-state flow from left to right (right to left). The next subsection describes a special case where this question can be answered rigorously.

Iii-C1 Bidirectional flow with no Langmuir kinetics

In the case where , , i.e. a system with no internal attachments and detachments, Eq. (III-C) becomes

(18)
Proposition 5

Consider the case where , , and suppose that (9) holds. Then

(19)

In particular, if  then , and

(20)

Eq. (19) means that in the case of no Langmuir kinetics the steady-state output from the right hand-side of the chain will be positive [negative] if the the product of the forward rates is larger [smaller] than the product of the backward rates. In transcription and translation the steady state flow from the right hand-side of the chain should always be positive, but in other cases, e.g. transport along microtubules, the steady state flow may be either positive or negative.

Iii-D Entrainment

Assume now that some or all of the rates are time-varying periodic functions with the same period . This may be interpreted as a periodic excitation of the system. Many biological processes are affected by such excitations due for example to the periodic 24h solar day or the periodic cell-cycle division process. For example, translation elongation factors, tRNAs, translation and transcription initiation factors, ATP levels, and more may change in a periodic manner and affect various rates that appear in the MFALK.

A natural question is will the state-variables of the MFALK converge to a periodic pattern with period ? We will show that this is indeed so, i.e. the MFALK entrains to a periodic excitation in the rates. In order to understand what this means, consider a different setting, namely, using the MFALK to model traffic flow. Then the rates may correspond to traffic lights, changing in a periodic manner, and the state-variables are the density of the moving particles (cars) along different sections of the road, so entrainment corresponds to what is known as the “green wave” (see e.g. [26] and the references therein).

We say that a function is -periodic if for all . Assume that the s, s, s and s are uniformly bounded, non-negative, time-varying functions satisfying:

  • there exists a (minimal) such that all the s, s, s, and s are -periodic.

  • there exist  such that at least one of the following two conditions holds for all time 

    (21)
    (22)
  • there exists  such that

    (23)

We refer to this model as the Periodic MFALK (PMFALK).

Theorem 1

Consider the PMFALK with dimension . There exists a unique function , that is -periodic, and for any  the trajectory  converges to  as .

Thus, the PMFALK entrains (or phase-locks) to the periodic excitation in the parameters. In particular, this means that the output rate in (3) converges to the unique -periodic function:

Note that since a constant function is a periodic function for all , Theorem 1 implies that entrainment holds also in the particular case where a single parameter is oscillating (with period ), while all other parameters are constant. Note also that Corollary 1 follows from Theorem 1.

Example 3

Consider the MFALK with dimension , parameters: , , , , , , , , , , , , , , and initial condition . Note that all the rates here are periodic, with a minimal common period . Fig. 4 depicts , , as a function of . It may be seen that each state variable converges to a periodic function with period 

Fig. 4: State variables [solid line]; [dashed line]; and [dotted line] as a function of in Example 3. Note that each state variable converges to a periodic function with a period .

Iii-E Strong Monotonicity

Recall that a proper cone  defines a partial ordering in  as follows. For two vectors , we write  if ; if  and ; and  if . The system  is called monotone if  implies that  for all . In other words, the flow preserves the partial ordering [60]. It is called strongly monotone if  implies that  for all .

From here on we consider the particular case where the cone is . Then  if  for all , and  if  for all . A system that is monotone with respect to this partial ordering is called cooperative.

Proposition 6

For any , with , the solutions of the MFALK satisfy

(24)

Furthermore, if (9) holds then

(25)

To explain this, consider two initial densities and  with  for all , that is, corresponds to a larger or equal density at each site. Then the trajectories  and  emanating from these initial conditions continue to satisfy the same relationship between the densities, namely, , for all  and for all time .

The MFALK is thus a strongly cooperative tridiagonal system (SCTS) on . Some of the properties deduced above using contraction theory can also be deduced using this property [59].

Remark 1

Suppose that we augment the MFALK into a model of  ODEs in  state-variables by adding to it the equation

that is,  (see (II)). Let  denote the vector of the  state-variables. Clearly, this augmented model admits a first integral . Also, for any initial condition in  all the state-variables remain bounded, as the first  state-variables remain in  and  for all . It is straightforward to verify that the augmented system is a cooperative system, and that if (9) holds then it is a SCTS. SCTS systems that admit a non-trivial first integral have many desirable properties (see, e.g. [40]).

Iii-F Effect of attachment and detachment

One may perhaps expect that detachment from a jammed site may increase the total flow by reducing congestion. The next result shows that this is not so. Detachment always decreases the steady-state production rate . Similarly, attachment always increases .

Proposition 7

Consider a MFALK with dimension . Suppose that the conditions in Proposition 3 hold. Then , and , for all . Also, , and  for all .

This means that an increase in any of the detachment [attachment] rates decreases [increases] the steady-state density in all the sites. Also, an increase in any of the internal detachment [attachment] rates decreases [increases] the steady-state production rate. The next example demonstrates this.

Example 4

Consider the MFALK with , , , , , . Fig. 5 depicts as a function of and . It may be seen that decreases with both  and 

Fig. 5: as a function of and for the MFALK in Example 4.

We note that the analytical results in Proposition 7 agree well with the simulation results obtained using a TASEP model for translation that included alternative initiation along the mRNA and ribosome drop-off [74].

The next section describes an application of the MFALK to a biological process.

Iv An application: modeling mRNA translation with ribosome drop off

It is believed that during mRNA translation ribosome movement is unidirectional from the 5’ end to the 3’ end, and that ribosomes do not enter in the middle of the coding regions. However, ribosomes can detach from various sites along the mRNA molecule due for example to collisions between ribosomes. This is known as ribosome drop off.

As mentioned in the introduction, ribosome drop off has been the topic of numerous studies [57, 25, 24, 68, 7, 61, 18, 28, 22, 20, 29]. It was suggested that in some cases ribosome drop off is important for proof reading [69], and also that ribosome stalling/abortion plays a role in translational regulation (e.g. see [56, 74]).

It is clear that ribosome abortion has drawbacks. Indeed, translation is the most energetically consuming process in the cell, and abortion results in truncated, non-functional and possibly deleterious proteins. It is believed that transcripts undergo evolutionary selection to minimize abortion and/or its energetic cost [67, 63, 61, 18, 28, 20]. Nevertheless, there seems to be a certain minimal abortion rate even in non-stressed conditions [57, 29]. This basal value was estimated (see more details below) to be of the order or abortion events per codon in E. coli. In other words, in every codon one out of decoding ribosomes aborts. This value is non-negligible. If we consider a drop-off rate of per codon along a coding region of  codons (approximately the average coding region length for E. coli) then on average, around out of every ribosomes will fail to complete the translation of the mRNA.

To model translation with ribosome drop off, we use the MFALK with  (i.e. no backwards motion) and  (i.e. no attachment to internal sites along the chain) for all . Changing the values of the s allows to model and analyze the effect of ribosome drop off at different sites along the mRNA molecule. We assume that

(26)

as otherwise the chain decouples into two smaller, disconnected chains. Note that (26) implies that the conditions in Proposition 3 hold, so the model is SO on  w.r.t. the  norm, and thus admits a unique globally asymptotically stable equilibrium point .

We study the effect of ribosome drop off on the steady-state protein production rate and ribosome density using real biological data. To this end, we considered S. cerevisiae genes (see Figures 6 and 7) with various mRNA levels (all genes were sorted according to their mRNA levels and genes were uniformly sampled from the list). Similarly to the approach used in [48], we divided the mRNAs related to these genes to non-overlapping pieces. The first piece includes the first codons that are related to various stages of initiation [63]. The other pieces include  non-overlapping codons each, except for the last one that includes between  and codons.

To model the translation dynamics in these mRNAs using MFALK, we model every piece of mRNA as a site. We estimated the elongation rates  at each site using ribo-seq data for the codon decoding rates [12], normalized so that the median elongation rate of all S. cerevisiae mRNAs becomes  codons per second [23]. The site rate is , where site time is the sum over the decoding times of all the codons in the piece of mRNA corresponding to this site. These rates thus depend on various factors including availability of tRNA molecules, amino acids, Aminoacyl tRNA synthetase activity and concentration, and local mRNA folding [12, 1, 63].

The initiation rate  (that corresponds to the first piece) was estimated based on the ribosome density per mRNA levels, as this value is expected to be approximately proportional to the initiation rate when initiation is rate limiting [48, 36]. Again we applied a normalization that brings the median initiation rate of all S. cerevisiae mRNAs to be [9].

We analyzed the effect of uniform ribosome drop off with a rate in the range of to per codon. This corresponds to , i.e., all the s are equal, and denote their common value. Since we assumed codons per site, values range from to (ten times the rate associated with a single codon). This makes sense as in the MFALK the level of occupancy in a site is related to the probability to see a ribosome in this site.

Let

denote the steady-state mean ribosomal density. Figures 6 and 7 depict and  in our model as a function of . In these figures the genes in the legends are sorted according to their expression levels: the gene at the top (YGR192C) has the highest mRNA levels while the gene at the bottom (YER106W) has the lowest levels. It may be seen that as the drop off (detachment) rate  increases from to , decreases by about , and decreases by about . This demonstrate the significant ramifications that ribosomal drop off is expected to have on translation and the importance of modeling drop off.

Note also that there is a strong variability in the effect of drop off on the different genes: for mRNAs with higher expression levels (i.e. mRNAs with higher copy number in the cell) the drop off effect is weaker. It is possible that this is related to stronger evolutionary selection for lower drop off rate in genes with higher mRNA levels. Indeed, highly expressed genes “consume” more ribosomes (due to higher mRNA levels), so a given (per-mRNA) drop off rate is expected to be more deleterious to the cell, and a mutation which decreases the drop of rate in such genes has a higher probability of fixation.

Fig. 6: Reduction in steady-state mean density in percent as a function of for S. cerevisiae genes.
Fig. 7: Reduction in steady-state output rate (production rate) in percent as a function of for S. cerevisiae genes.

V Discussion

In many important processes biological “particles” move along some kind of a one-dimensional “track”. Examples include gene transcription and translation, cellular transport, and more. The flow can be either bidirectional (as in the case of transcription) or unidirectional (as in the case of translation), with the possibility of both attachment and detachment of particles at different sites along the track. For example, motor proteins like kinesin and dynein that move along a certain microtubule may detach and attach to an overlapping microtubule.

To rigorously model and analyze such processes, we introduced a new deterministic mathematical model that can be derived as the dynamic mean-field approximation of ASEP with Langmuir kinetics, called the MFALK. Our main results show that the MFALK is a monotone and contractive dynamical system. This implies that it admits a globally asymptotically unique equilibrium point, and that it entrains to periodic excitations (with a common period ) in any of its rates, i.e. the densities along the chain, as well as the output rate, converge to unique period solutions with period .

It is important to note that several known models are special cases of the MFALK. These include for example the RFM [48], the model used in [15] for DNA transcription, and the model of phosphorelays in [11].111Although in this model the occupancy levels are normalized differently.

Topics for further research include the following. In the RFM, it has been shown that the steady-state production rate is related to the maximal eigenvalue of a certain non-negative, symmetric tridiagonal matrix with elements that are functions of the RFM rates, i.e. the [44]. This implies that the mapping is strictly concave, and that sensitivity analysis of is an eigenvalue sensitivity problem [45]. An interesting research topic is whether in the MFALK can also be described using such a linear-algebraic approach.

The application of the MFALK to model ribosome drop off suggests an interesting direction for further study, namely, how to design genes that minimize the drop off rate.

Another research direction is motivated by the fact that many of the transport phenomena that can be modeled using the MFALK do not take place in isolation. For example, many mRNA molecules are translated in parallel in the cell. Thus, a natural next step is to study networks of interconnected MFALKs. Graph theory can be used to describe the interconnections between the various MFALKs in the network. In this context, ribosome drop off may perhaps increase the total production rate in the entire system, as it allows ribosomes to detach from slow sites, enter the pool of free ribosomes, and then attach to the initiation sites of other, less crowded, mRNA molecules. However, drop off still incurs the biological “cost” associated to the synthesis of a chain of amino-acids that is only a part of the desired protein. The fact that the MFALK is contractive may prove useful in analyzing networks of MFALKs, as there exist interesting results proving the overall contractivity of a network based on contractivity of the subsystems and their couplings (see, e.g. [5, 50]).

Another interesting topic for further research is studying the effect of controlled detachment rates on the formation of traffic jams. Indeed, it is known that kinesin-family motor proteins are more susceptible to dissociation when their pathway is blocked [14, 62].

Appendix: Proofs

We begin by discussing some symmetry properties of the MFALK, as these will be useful in the proofs later on.

Symmetry

The MFALK enjoys two symmetries that will be useful later on. First, let , . In other words,  is the amount of “free space” at site  at time . Then using (II) yields

(27)

This is just the MFALK (II), but with the parameters permuted as follows: , , , and  for all . The symmetry here follows from the fact that we can replace the roles of the forward and backward flows in the MFALK.

Next, let , . In other words,  is the amount of “free space” at site  at time . Then using (II) yields

(28)

This is just the MFALK (II), but with the parameters permuted as follows: , , , and  for all . Note that (Symmetry) is simply (Symmetry) with the variable renaming , .

Both symmetries are reminiscent of the particle-hole symmetry in ASEP [6, 30]: the basic idea is that the progression of a particle from left to right is also the progression of a hole from right to left.

Proof of Proposition 1. If (4) holds then the MFALK satisfies property (BR) in [35], and [35, Lemma 1] implies (6). If (5) holds then (Symmetry) satisfies property (BR) in [35], and this implies (6). ∎

Proof of Proposition 2. Write the MFALK as . A calculation shows that the Jacobian matrix  satisfies