A Deterministic Mathematical Model for Bidirectional Excluded Flow with Langmuir Kinetics^{†}^{†}thanks: 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.
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 onedimensional 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 twodirectional, 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 meanfield 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 steadystate 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 steadystate protein production rate.
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).
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 bidirectional, 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 twodirectional 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 kinesin8 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 ribosomeribosome interactions or due to depletion in the concentration of tRNAs [72, 28, 57]. Also, it is known that kinesinfamily 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 slowmoving 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 TASEPlike 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 nonequilibrium 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 meanfield 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 statevariables 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 celldivision process), and proper functioning requires phaselocking 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 meanfield 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 largescale 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, nonfunctional 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 steadystate 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 firstorder 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 nonnegative 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 lefttoright 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 statevariable 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 righthand 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.
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 meanfield 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 statevariables 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
Iiia 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 .
IiiB 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 timevarying 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
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.
IiiC 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 .
At steadystate, i.e. for , the lefthand 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 steadystate output rate. Then , for all , that is, the steadystate 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 steadystate 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 (IIIC) 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 steadystate output rate is .
In general, (IIIC) 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. (IIIC) 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 steadystate 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 steadystate flow from left to right (right to left). The next subsection describes a special case where this question can be answered rigorously.
IiiC1 Bidirectional flow with no Langmuir kinetics
In the case where , , i.e. a system with no internal attachments and detachments, Eq. (IIIC) becomes
(18) 
Proposition 5
Eq. (19) means that in the case of no Langmuir kinetics the steadystate output from the right handside 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 handside 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.
IiiD Entrainment
Assume now that some or all of the rates are timevarying 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 cellcycle 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 statevariables 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 statevariables 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, nonnegative, timevarying 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 phaselocks) 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 .
IiiE 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
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 statevariables by adding to it the equation
that is, (see (II)). Let denote the vector of the statevariables. Clearly, this augmented model admits a first integral . Also, for any initial condition in all the statevariables remain bounded, as the first statevariables 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 nontrivial first integral have many desirable properties (see, e.g. [40]).
IiiF 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 steadystate 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 steadystate density in all the sites. Also, an increase in any of the internal detachment [attachment] rates decreases [increases] the steadystate 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 .
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 dropoff [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, nonfunctional 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 nonstressed 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 nonnegligible. If we consider a dropoff 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 steadystate 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 nonoverlapping pieces. The first piece includes the first codons that are related to various stages of initiation [63]. The other pieces include nonoverlapping 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 riboseq 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 steadystate 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 (permRNA) 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.
V Discussion
In many important processes biological “particles” move along some kind of a onedimensional “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 meanfield 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].^{1}^{1}1Although 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 steadystate production rate is related to the maximal eigenvalue of a certain nonnegative, symmetric tridiagonal matrix with elements that are functions of the RFM rates, i.e. the s [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 linearalgebraic 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 aminoacids 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]).
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 particlehole 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