## Abstract

A single mammalian cell includes an order of 10^{4}–10^{5} mRNA molecules and as many as 10^{5}–10^{6} ribosomes. Large-scale simultaneous mRNA translation induces correlations between the mRNA molecules, as they all compete for the finite pool of available ribosomes. This has important implications for the cell's functioning and evolution. Developing a better understanding of the intricate correlations between these simultaneous processes, rather than focusing on the translation of a single isolated transcript, should help in gaining a better understanding of mRNA translation regulation and the way elongation rates affect organismal fitness. A model of simultaneous translation is specifically important when dealing with highly expressed genes, as these consume more resources. In addition, such a model can lead to more accurate predictions that are needed in the interconnection of translational modules in synthetic biology. We develop and analyse a general dynamical model for large-scale simultaneous mRNA translation and competition for ribosomes. This is based on combining several ribosome flow models (RFMs) interconnected via a pool of free ribosomes. We use this model to explore the interactions between the various mRNA molecules and ribosomes at steady state. We show that the compound system always converges to a steady state and that it always entrains or phase locks to periodically time-varying transition rates in any of the mRNA molecules. We then study the effect of changing the transition rates in one mRNA molecule on the steady-state translation rates of the other mRNAs that results from the competition for ribosomes. We show that increasing any of the codon translation rates in a specific mRNA molecule yields a local effect, an increase in the translation rate of this mRNA, and also a global effect, the translation rates in the other mRNA molecules all increase or all decrease. These results suggest that the effect of codon decoding rates of endogenous and heterologous mRNAs on protein production is more complicated than previously thought. In addition, we show that increasing the length of an mRNA molecule decreases the production rate of all the mRNAs.

## 1. Introduction

Various processes in the cell use the same finite pool of available resources. This means that the processes actually compete for these resources, leading to an indirect coupling between the processes. This is particularly relevant when many identical intracellular processes, all using the same resources, take place in parallel. The indirect coupling induced by competition may affect the cell's structure and functioning in ways that cannot be understood when studying each process as a single, isolated process.

Biological evidence suggests that the competition for RNA polymerase (RNAP) and ribosomes, and various transcription and translation factors, is a key factor in the cellular economy of gene expression. The limited availability of these resources is one of the reasons why the levels of genes, mRNA and proteins produced in the cell do not necessarily correlate [1–7].

It was estimated that in a yeast cell there are approximately 60 000 mRNA molecules. These can be translated in parallel [8,9], with possibly many ribosomes scanning the same transcript concurrently. The number of ribosomes is limited (in a yeast cell it is approx. 240 000) and this leads to a competition for ribosomes. For example, if more ribosomes bind to a certain mRNA molecule then the pool of free ribosomes in the cell is depleted, and this may lead to lower initiation rates in the other mRNAs (figure 1).

There is a growing interest in computational or mathematical models that take into account the competition for available resources in translation and/or transcription (see e.g. [10–15]). One such model that explicitly considers the movement of the ribosomes (RNAP) along the mRNA (DNA) is based on a set of *asymmetric simple exclusion processes* (ASEPs) interconnected to a pool of free ribosomes. ASEP is an important model from non-equilibrium statistical physics describing particles that hop randomly from one site to the next along an ordered lattice of sites, but only if the next site is empty. This form of ‘rough exclusion’ models the fact that the particles cannot overtake one another. ASEP has been used to model and analyse numerous multiagent systems with local interactions including the flow of ribosomes along the mRNA molecule [16,17]. In this context, the lattice represents the mRNA molecule, and the particles are the ribosomes. For more on mathematical and computational models of translation, see the survey paper [18].

Ha & den Nijs [19] considered a closed system composed of a single ASEP connected to a pool (or reservoir) of ‘free’ particles. The total number of particles is conserved. This is sometimes referred to as the *parking garage problem*, with the lattice/particles/pool modelling a road/cars/parking garage, respectively. Cook & Zia [20] studied a similar system using domain wall theory. Cook & Zia [21] (see also [22]) considered a network composed of two ASEPs connected to a finite pool of particles. The analysis in these papers focuses on the phase diagram of the compound system with respect to certain parameters, and on how the phase of one ASEP affects the phase of the other ASEPs. These studies rely on the phase diagram of a single ASEP that is well understood only in the case where all the transition rates inside the chain (the elongation rates) are equal. Thus, the network is typically composed of *homogeneous* ASEPs. Another model [15] combines non-homogeneous ASEPs in order to study competition between multiple species of mRNA molecules for a pool of tRNA molecules. This study was based on the *Saccharomyces cerevisiae* genome. However, in this case (and similar models, such as [14]) analysis seems intractable and one must resort to simulations only.

Our approach is based on the *ribosome flow model* (RFM) [23]. This is a deterministic, continuous-time, synchronous model for translation that can be derived via the mean-field approximation of ASEP [24]. The RFM includes *n* state variables describing the ribosomal density in *n* consecutive sites along the mRNA molecule, and positive parameters: the translation initiation rate , and the translation elongation rate from site *i* to site , for . The RFM has a unique equilibrium point , and any trajectory emanating from a feasible initial condition converges to *e* [25] (see also [26]). This means that the system always converges to a steady-state ribosomal density that depends on the rates, but not on the initial ribosomal density profile. In particular, the production rate converges to a steady-state value denoted by *R*. The mapping from the rates to *R* is a concave function, so maximizing *R* subject to a suitable constraint on the rates is a convex optimization problem [27,28]. Sensitivity analysis of the RFM with respect to the rates has been studied in Poker *et al.* [29]. These results are important in the context of optimizing the protein production rate in synthetic biology. Margaliot *et al.* [26] have shown that when the rates are time-periodic functions, with a common minimal period *T*, then every state variable converges to a periodic solution with period *T*. In other words, the ribosomal densities *entrain* to periodic excitations in the rates (due e.g. to periodically varying abundances of tRNA molecules).

In ASEP with *periodic boundary conditions*, a particle that hops from the last site returns to the first one. The mean-field approximation of this model is called the *ribosome flow model on a ring* (RFMR). The periodic boundary conditions mean that the total number of ribosomes is conserved. Raveh *et al.* [30] analysed the RFMR using the theory of monotone dynamical systems that admit a first integral.

Both the RFM and the RFMR model mRNA translation on a single mRNA molecule. In this paper, we introduce a new model, called the *RFM network with a pool* (RFMNP), that includes a network of RFMs, interconnected through a dynamical pool of free ribosomes, to model and analyse simultaneous translation and competition for ribosomes in the cell. To the best of our knowledge, this is the first study of a network of RFMs. The total number of ribosomes in the RFMNP is conserved, leading to a first integral of the dynamics. Applying the theory of monotone dynamical systems that admit a first integral, we prove several mathematical properties of the RFMNP: it admits a continuum of equilibrium points, every trajectory converges to an equilibrium point, and any two solutions emanating from initial conditions corresponding to an equal total number of ribosomes in the system converge to the same equilibrium point. These results hold for any set of rates and in particular when the RFMs in the network are not necessarily homogeneous. These stability results are important because they provide a rigorous framework for studying questions such as how does a change in one RFM affects the behaviour of all the other RFMs in the network? Indeed, since a steady-state exists, this can be reduced to asking: how does a change in one RFM in the network affects the *steady-state* behaviour of the network?

To analyse competition for ribosomes, we consider the effect of increasing one of the rates in one RFM, say RFM #1. This means that the ribosomes traverse RFM #1 more quickly. We show that this always leads to an increase in the production rate of RFM #1. All the other RFMs are always affected in the same manner, that is, either all the other production rates increase or they all decrease. Our analysis shows that this can be explained as follows. Increasing the rate in RFM #1 tends to increase the steady-state density in sites , and decrease the density in site *i* of this RFM. The *total* density (i.e. the sum of all the densities on the different sites along RFM #1) can either decrease or increase. In the first case, more ribosomes are freed to the pool, and this increases the initiation rates in all the other RFMs leading to higher production rates. The second case leads to the opposite result. The exact outcome of increasing one of the rates thus depends on the many parameter values defining the pool and the set of RFMs in the network.

Our model takes into account the dynamics of the translation elongation stage, yet is still amenable to analysis. This allows to develop a rigorous understanding of the effect of competition for ribosomes. Previous studies on this topic were either based on simulations (see, for example, [14,15]) or did not include a dynamical model of translation elongation (e.g. [11,13]). For example, in an interesting paper, combining mathematical modelling and biological experiments, Gyorgy *et al.* [11] study the expression levels of two adjacent reporter genes on a plasmid in *Escherichia coli* based on measurements of fluorescence levels, that is, protein levels. These are of course the result of all the gene expression steps (transcription, translation, mRNA degradation, protein degradation) making it difficult to separately study the effect of competition for ribosomes or to study specifically the translation elongation step. Their analysis yields that the attainable output of the two proteins satisfies the formula1.1where *Y* is related to the total number of ribosomes (but also other translation factors and possibly additional gene expression factors), and are constants that depend on parameters such as the plasmid copy number, dissociation constants of the ribosomes binding to the ribosomal binding site (RBS), etc. This equation implies that increasing the production of one protein always leads to a decrease in the production of the other protein (although more subtle correlations may take place via the effects on the constants *α* and *β*). A similar conclusion also has been derived for other models ([31], ch. 7).

In our model, improving the translation rate of a codon in one mRNA may either increase or decrease the translation rates of *all* other mRNAs in the cell. Indeed, the effect on the other genes depends on the change in the *total density* of ribosomes on the modified mRNA molecule, highlighting the importance of modelling the dynamics of the translation elongation step. We show, however, that when increasing the initiation rate in an RFM in the network, the total density in this RFM always increases and, consequently, the production rate in all the other RFMs decreases. This special case agrees with the results in [11].

Another recent study [32] showed that the hidden layer of interactions among genes arising from competition for shared resources can dramatically change network behaviour. For example, a cascade of activators can behave like an effective repressor, and a repression cascade can become bistable. This agrees with several previous studies in the field (e.g. [7,33]).

The remainder of this paper is organized as follows. Section 2 summarizes the biological implications of our theoretical results. This includes no equations, and was written in order to make this paper more accessible to a wider audience. Section 3 describes the new model, and demonstrates using several examples how it can be used to study translation at the *cell level*. Section 4 describes our main theoretical results, and details their biological implications. To streamline the presentation, the proofs of the results are in the electronic supplementary material. The final section concludes and describes several directions for further research, including suggestions for biological experiments inspired from the theoretical results.

## 2. Biological rationale and summary of the main results

The goal of our study is to understand the intricate coupling between the translation rates of different mRNA molecules induced by the competition for the finite pool of available resources and, in particular, ribosomes. We develop a new mathematical model that includes a finite pool of ‘free’ ribosomes, and an arbitrary number of mRNA molecules. This implies that the different mRNA molecules ‘compete’ for the finite pool of ribosomes (figure 1).

Our model encapsulates the fundamental dynamical aspects of the translation process: the codon decoding rates and initiation rates may vary among different genes, and within a coding region; the movement of ribosomes is unidirectional from the 5′ end to the 3′ of the coding region; when there are multiple ribosomes on the mRNA, a ribosome cannot overtake another ribosome, etc.

The initiation and codon translation rates in the model are related to various biophysical mechanisms and transcript features including adaptation of codons to the tRNA pool (codons that are recognized by a tRNA with a higher intracellular abundance tend to be translated more quickly [34]); local folding of the mRNA (stronger folding tends to decrease elongation ratex [35]); interaction/hybridization between the ribosomal RNA and the mRNA [36] (there are nucleotide sub-sequences that tend to interact with the ribosomal RNA); the interaction between the nascent peptide and the exit tunnel of the ribosome [37,38]; and more. The initiation rates may also be related to various mechanisms and transcript features including ribosome binding sites in prokaryote [39]; mRNA folding near the start codon [2,40] and nucleotide composition surrounding the start codon in eukaryotes [41,42]. However, the theoretical results proved here hold for any set of feasible parameter values covering any possible physiological condition. For example, they hold when the initiation rate is very low (and thus rate limiting), but also when the initiation rate is high, causing ribosomal ‘traffic jams’ along the mRNA. Similarly, the results hold if the codon decoding rates are constant along different codons and when they significantly vary along the codons. Thus, the reported results are relevant for various organisms (e.g. bacteria, fungi, mammals) in various regimes and conditions (e.g. normal conditions, starvation, viral infection and diseases such as cancer, etc.). In addition, they are relevant both for endogenous and heterologous/engineered genes.

We show that the compound system always converges to a steady state. In other words, the ribosome density profile on every mRNA molecule and the number of ribosomes in the pool converge, as time goes to infinity, to some fixed value. We then analyse how this steady state changes as we change a transition rate in one of the mRNA molecules. Note that the change in one mRNA induces a change throughout the system because of the mutual coupling via the competition for the available ribosomes.

The first fundamental conclusion of our study describes the local and global effects of modifying the *initiation rate* in one mRNA molecule on the total protein production in the cell (figure 2*a*). We show that increasing [decreasing] the initiation rate in a certain mRNA molecule always increases [decreases] the ribosome density profile and the translation rate of that mRNA, and also decreases [increases] the ribosome density and translation rate in *all* the other mRNA molecules. In other words, the local and global effects are opposed. This conclusion is important in the context of biotechnology and synthetic biology, as it implies that increasing the initiation rate of an engineered (e.g. heterologous) mRNA molecule will lead to a decrease in the translation rates of all the other mRNA molecules. This is expected to yield a decrease in the growth rate of the host, and may eventually decrease the production rate of the heterologous gene. This conclusion is also relevant in the context of molecular evolution, as it shows that a mutation near the 5′ end of the coding region in one mRNA affects the organismal fitness not only locally but also globally.

The second fundamental conclusion of our study considers increasing [decreasing] a codon elongation rate (but not the initiation rate) in a specific mRNA molecule. The local effect is *always* an increase [decrease] in the translation rate of the modified mRNA. The global effect is homogeneous: either all the translation rates of the other mRNAs increase or they all decrease (figures 2*b* and 3*b*). Which of these two outcomes will indeed take place depends on all the rates in the system, and can be predicted using our model. We can provide an intuitive explanation of the two possible outcomes as follows. Suppose that the modified rate is a bottleneck rate, i.e. it leads to ribosomal traffic jams. Then increasing this rate reduces the traffic jams, thus freeing more ribosomes to the pool. This increases the initiation rates in *all* the other mRNAs leading to a global increase in translation rates (figure 2*b*).

On the other hand, suppose that the modified codon is located close to the beginning of the coding region and that there is a bottleneck codon, generating a traffic jam, further away along the coding region. Then an increase in the transition rate of the modified codon actually worsens the traffic jam. This depletes the pool, and thus yields a global decrease in the translation rates of all the other mRNAs (figure 3*b*).

This second conclusion is important in the context of biotechnology and synthetic biology, as it demonstrates that evaluating the effect of (heterologous or endogenous) gene engineering on its production rate should take into account the change in the ribosomal pool in order to estimate the global, and not only local, effect of the modification. Our model allows to predict this effect.

Specifically, consider the effect of introducing a ‘ramp’, i.e. a region with relatively slow codons at the beginning of the coding region [1,4]. This is expected to reduce the traffic jams in the modified mRNA because fewer ribosomes can enter and pile up in the modified mRNA. This then should increase the number of free ribosomes and thus the global effect is expected to be an increase in the initiation and translation rates in all the other mRNAs, leading to an increase in organismal growth rate and fitness (figure 3*a*).

This conclusion is also relevant in the context of molecular evolution, as it allows a better understanding of the *trans* effect of mutations in the coding region on the organismal fitness via the global, and not only local, effect of the modification.

Our model and its analysis also suggest several novel biological experiments. Some of these are described in §5.

## 3. The model and some examples

Since our model is based on a network of interconnected RFMs, we begin with a brief review of the RFM.

### 3.1. Ribosome flow model

The RFM models the traffic flow of ribosomes along the mRNA (figure 4). The mRNA chain is divided into a set of *n* compartments or sites, where each site may correspond to a codon or a group of codons. The state variable , , describes the ribosome occupancy at site *i* at time *t*, normalized such that [] implies that site *i* is completely empty [full] at time *t*. Since takes values in [0, 1], one may also view as the probability that site *i* is occupied at time *t*. The dynamical equations of the RFM are3.1

These equations describe the movement of ribosomes along the mRNA chain. The *transition rates* are all positive numbers (units = 1/time). To explain this model, consider the equation . The term represents the flow of particles from site 1 to site 2. This is proportional to the occupancy at site 1 and also to , i.e. the flow decreases as site 2 becomes fuller. In particular, if , i.e. site 2 is completely full, the flow from site 1 to site 2 is zero. This is a ‘soft’ version of the rough exclusion principle in ASEP. Note that the maximal possible flow rate from site 1 to site 2 is the transition rate . The term represents the flow of particles from site 2 to site 3.

The dynamical equations for the other state variables are similar. Note that controls the initiation rate into the chain, and that
is the rate of flow of ribosomes out of the chain, that is the translation (or protein production) rate at time *t*. The RFM topology is depicted in figure 4.

The RFM encapsulates simple exclusion and unidirectional movement along the lattice just as in ASEP. This is not surprising, as the RFM can be derived via a mean-field approximation of ASEP (e.g. [24], p. R345, and [43], p. 1919). However, the analysis of these two models is quite different, as the RFM is a deterministic, continuous-time, synchronous model, whereas ASEP is a stochastic, discrete-type, asynchronous one.

In order to study a network of interconnected RFMs, it is useful to first extend the RFM into a single-input single-output control system:3.2

Here the translation rate becomes the output *y* of the system, and the flow into site 1 is multiplied by a time-varying control , representing the flow of ribosomes into the mRNA from the ‘outside world’. The time-varying function is thus related to the rate ribosomes diffuse to the 5′ end (in eukaryotes) or the RBS (in prokaryotes) of the mRNA at time *t*. Of course, mathematically one can absorb into *u*, but we do not do this because we think of as representing some intrinsic local mRNA-specific features (e.g. the strength of the Kozak sequences in eukaryotes or the RBS in prokaryote).

The set of admissible controls is the set of bounded and measurable functions taking values in for all time . Equation (3.2), referred to as the RFM with input and output (RFMIO) [44], facilitates the study of RFMs with feedback connections. We note in passing that (3.2) is a *monotone control system* as defined in [45]. From now on we write (3.2) as3.3

Let denote the closed unit cube in . Since the state variables in the RFM represent normalized occupancy levels, we always consider initial conditions . It is straightforward to verify that is an invariant set of (3.3), i.e. for any and any , the trajectory satisfies for all . This means that if the initial density is well defined then it remains so for all time.

### 3.2. Ribosome flow model network with a pool

To model competition for ribosomes in the cell, we consider a set of RFMIOs, representing *m* different mRNA molecules in the cell (figure 5). The *i*th RFMIO has length , input function , output function and rates . The dynamics of these RFMIOs is thus given by3.4

The *m* RFMIOs are interconnected through a pool of free ribosomes (i.e. ribosomes that are not attached to any mRNA molecule). We use to denote the pool occupancy at time *t* (this may also be interpreted as the average number of free ribosomes in the pool at time *t*). The pool feeds the initiation location in all the mRNAs (figure 5). To model this we can simply take for all *i* (recall that describes the rate at which ribosomes from the ‘outside world’ feed RFMIO #i), but we allow a somewhat more general setting modelled by3.5

We assume that each satisfies the following properties: (i) ; (ii) is continuously differentiable and for all (so is strictly increasing on ) and (iii) there exists such that for all sufficiently small. The first property implies that if the pool is empty then no ribosomes can bind to the mRNA molecules, and the second means that as the pool becomes fuller the initiation rates to the RFMIOs increase. The third property is a technical assumption that is needed for one of our proofs below.

Typical examples for functions satisfying these properties include the linear function, say, , and also the function , with . In the first case, the flow of ribosomes into the first site of RFM is given by , and the product here can be justified via mass-action kinetics. The use of tanh may be suitable for modelling a saturating function. This is in fact a standard function in ASEP models with a pool [21,46], because it is zero when *z* is zero, uniformly bounded and strictly increasing for . Also, for the function takes values in so it can also be interpreted as a probability function [47].

In the context of a shared pool, it is natural to consider the special case where for all . This assumption is reasonable since there is no *a priori* reason to expect a ‘preference’ between different mRNA species. However, the differences between the initiation sites in the mRNA molecules may still be modelled by taking different s.

The output of each RFMIO is fed into the pool, and the pool feeds the initiation locations in the mRNAs (figure 5). Thus, the pool dynamics is described by3.6

The above equation means that the flow into the pool is the sum of all output rates of the RFMIOs minus the total flow of ribosomes that bind to an mRNA molecule . Recall that the term represents the exclusion, i.e. as the first site in RFMIO #*j* becomes fuller, less ribosomes can bind to it. Summarizing, the RFMNP is given by equations (3.4), (3.5) and (3.6). This is a dynamical system with state variables.

Note that combining the properties of the functions with (3.6) implies that if then for all . Thus, the pool occupancy is always non-negative.

### Example 3.1.

Consider a network with RFMIOs, the first [second] with dimension []. Then the RFMNP is given by3.7

Note that this system has state variables. ▪

An important property of the RFMNP is, that being a closed system, the total occupancy3.8is conserved (i.e. the total number of ribosomes is conserved), that is,3.9

In other words, *H* is a first integral of the dynamics. In particular, this means that for all , i.e. the pool occupancy is uniformly bounded by the constant .

The RFMNP models mRNAs that compete for ribosomes because the total number of ribosomes is conserved. As more ribosomes bind to the RFMIOs, the pool depletes, decreases and the effective initiation rate to all the RFMIOs decreases (figure 5). This allows to systematically address important biological questions on large-scale simultaneous translation under competition for ribosomes. The following examples demonstrate this. We prove in §4 that all the state variables in the RFMNP converge to a steady-state. Let denote the steady-state occupancy in site *j* in RFMIO #*i*, and let denote the steady-state occupancy in the pool. In the examples below we always consider these steady-state values (obtained numerically by simulating the differential equations).

### Example 3.2.

Although we are mainly interested in modelling large-scale simultaneous translation, it is natural to first consider a model with a single mRNA molecule connected to a pool of ribosomes. From a biological perspective, this models the case where there is one gene that is highly expressed with respect to all other genes (e.g. an extremely highly expressed heterologous gene).

Consider an RFMNP that includes a single RFMIO (i.e. ), with dimension , rates , , and a pool with output function . We simulated this system for the initial condition for all *i*, and for various values of *c*. Note that . Figure 6 depicts the steady-state values of the state-space variables in the RFMIO, and the steady-state pool occupancy . It may be seen that for small values of *c* the steady-state ribosomal densities and thus the production rates are very low. This is simply because there are not enough ribosomes in the network. The ribosomal densities increase with *c*. For large values of *c*, the output function of the pool saturates, as , and so does the initiation rate in the RFMIO. Thus, the densities in the RFMIO saturate to the values corresponding to the initiation rate , and then all the remaining ribosomes accumulate in the pool. Using a different pool output function, for example , leads to the same qualitative behaviour, but with higher saturation values for the ribosomal densities in the RFMIO. (Note that the ribosomal densities in an RFM are finite even when [48].) ▪

This simple example already demonstrates the coupling between the ribosomal pool, initiation rate and elongation rates. When the ribosomal pool is small the initiation rate is low. Thus, the ribosomal densities on the mRNA are low and there are no interactions between ribosomes (i.e. no ‘traffic jams’) along the mRNA. The initiation rate becomes the rate limiting step of translation. On the other hand, when there are many ribosomes in the pool the initiation rate increases, the elongation rates become rate limiting and ‘traffic jams’ along the mRNA evolve. At some point, a further increase in the number of ribosomes in the pool will have a negligible effect on the production rate.

It is known that there can be very large changes in the number of ribosomes in the cell during, e.g. exponential growth. For example, Bremer & Dennis [49] report changes in the range 6800–72 000. The example above demonstrates how these large changes in the number of ribosomes are expected to affect the translational regimes; specifically, it may cause a switch between the different regimes mentioned above.

The next example describes an RFMNP with several mRNA chains. Let denote the vector of *n* ones.

### Example 3.3.

Consider an RFMNP with RFMIOs of dimensions , and rates

In other words, every RFMIO has homogeneous rates. Suppose also that , for . We simulated this RFMNP for different values of *c* with the initial condition , , and . Thus, in all the simulations. For each value of *c*, every state variable in the RFMNP converges to a steady state. Figure 7 depicts the steady-state value and the steady-state output in each RFMIO. It may be seen that increasing *c*, i.e. increasing all the elongation rates in RFMIO #1 leads to an increase in the steady-state translation rates in *all* the RFMIOs in the network. Also, it leads to an increase in the steady-state occupancy of the pool. It may seem that this contradicts (3.9) but this is not so. Increasing *c* indeed increases all the steady-state translation rates, but it decreases the steady-state occupancies inside each RFMIO so that the total is conserved.

Let , i.e. the averaged steady-state occupancy (ASSO) on RFMIO . Figure 8 depicts the ASSO in each RFMIO as a function of *c*. It may be seen as *c* increases the ASSO in RFMIO #1 decreases quickly, yet the ASSOs in the other two RFMIOs slowly increase. Indeed, since the ribosomes spend less time on RFMIO #1 (due to increased *c*) they are now available for translating the other RFMIOs, leading to the increased ASSO in the other mRNAs. ▪

From a biological point of view this example corresponds to a situation where accelerating *one* of the mRNA chains increases the protein production rates in *all* the mRNAs and also increases the number of free ribosomes. Surprisingly, perhaps, it also suggests that a relatively larger number of free ribosomes in the cell corresponds to higher protein production rates. This agrees with evolutionary, biological and synthetic biology studies that have suggested that (specifically) highly expressed genes (that are transcribed into many mRNA molecules) undergo selection to include codons with improved elongation rates [1,2,50]. Specifically, two mechanisms by which improved codons affect translation efficiency and the organismal fitness are the following [50]: (i) *global mechanism*: selection for improved codons contributes toward improved ribosomal recycling and global allocation; the increased number of free ribosomes improves the effective translation initiation rate of all genes, and thus improves global translation efficiency, and (ii) *local mechanism*: the improved translation elongation rate of an mRNA contributes directly to its protein production rate.

The example above demonstrates both mechanisms, as improvement of the translation elongation rates of one RFM increases the translation rate of this mRNA (local translation efficiency), and also of the other RFMs (global translation efficiency). In addition, as can be seen, the decrease in ASSO in RFMIO #1 is significantly higher than the increase in ASSO in the other RFMIO. Thus, the simulation also demonstrates that increasing the translation rate *c* may contribute to decreasing ribosomal collision (and possibly ribosomal abortion).

We prove in §4 that when one of the rates in one of the RFMIOs increases two outcomes are possible: either all the production rates in the other RFMs increase (as in this example) or they all decrease. As discussed below, we believe that this second case is less likely to occur in endogenous genes, but may occur in heterologous gene expression.

The next example describes the effect of changing the length of one RFMIO in the network. Let denote a vector of *n* zeros.

### Example 3.4.

Consider an RFMNP with RFMIOs of dimensions and , rates
and , . In other words, both RFMIOs have the same homogeneous rates. We simulated this RFMNP for different values of with the initial condition , and . Thus, in all the simulations. For each value of , every state variable in the RFMNP converges to a steady state. Figure 9 depicts the steady-state values of *z*, and the steady-state output in each RFMIO. It may be seen that increasing , i.e. increasing the length of RFMIO #1 leads to a decrease in the steady-state production rates and in the steady-state pool occupancy. This is reasonable, as increasing means that ribosomes that bind to the first chain remain on it for a longer period of time. This decreases the production rate and, by the competition for ribosomes, also decreases the pool occupancy and thus decreases . ▪

From a biological point of view this suggests that decreasing the length of mRNA molecules contributes locally and globally to improving translation efficiency. A shorter coding sequence improves the translation rate of the mRNA and, by competition, may also improve the translation rates in all other mRNAs. Thus, we should expect to see selection for shorter coding sequences, specifically in highly expressed genes and in organisms with large population size. Indeed, previous studies have reported that in some organisms the coding regions of highly expressed genes tend to be shorter [14]; other studies have shown that other (non-coding) parts of highly expressed genes tend to be shorter [51–54]. Decreasing the length of different parts of the gene should contribute to organismal fitness via improving the energetic cost of various gene expression steps. For example, shorter genes should improve the metabolic cost of synthesizing mRNA and proteins; it can also reduce the energy spent for splicing and processing of RNA and proteins. However, there are of course various functional and regulatory constraints that also contribute to shaping the gene length (e.g. [55]). Our results and these previous studies suggest that in some cases genes are expected to undergo selection also for short *coding* regions, as this reduces the required number of translating ribosomes.

The next section describes various fundamental properties related to the RFMNP. All the proofs are placed in the Appendix (see the electronic supplementary material).

## 4. Properties of the ribosome flow model network with a pool

Let
denote the state space of the RFMNP (recall that every takes values in and ). For an initial condition , let denote the solution of the RFMNP at time *t*. It is straightforward to show that the solution remains in *Ω* for all .

### 4.1. Stability

For , let denote a *level set* of the first integral *H*:

In other words, means that *y* is a condition corresponding to a total occupancy of *s* ribosomes in the system.

### Theorem 4.1.

*Every level set* , , *contains a unique equilibrium point* *of the RFMNP*, *and for any initial condition* , *the solution of the RFMNP converges to* *. Furthermore*, *for any* ,4.1

In particular, this means that every trajectory converges to an equilibrium point, representing steady-state ribosomal densities in the RFMIOs and the pool. Equation (4.1) means that the continuum of equilibrium points, namely, , are linearly ordered. In other word, given any two equilibrium points, and , with , then every entry in is strictly smaller than the corresponding entry in . In other words, if the RFMNP is initiated with two initial conditions, with the first one corresponding to a higher number of ribosomes in the network than the second then the steady-state corresponding to the first condition will have a higher density than the steady state corresponding to the second in each site in each mRNA and in the pool.

It follows from proposition 4.8, stated in §4.3 below, that for any , . In other words, the steady state densities will never include a density that is either zero or one, and the steady-state pool occupancy is always strictly positive.

### Example 4.2.

Consider an RFMNP with RFMIOs with dimensions , and , , i.e.4.2

Note that even in this simple case the RFMNP is a nonlinear system. Assume that , and that , and denote this value simply by . Pick an initial condition in *Ω*, and let , so that the trajectory belong to for all . Any equilibrium point satisfies

This yields two solutions4.3and

It is straightforward to verify that in the latter solution , so this is not a feasible solution. The solution (4.3) does belong to , so the system admits a unique equilibrium in . Figure 10 depicts trajectories of (4.2) for three initial conditions in , namely , and , and the equilibrium point (4.3) for . It may be seen that every one of these trajectories converges to *e*. ▪

Various possible intracellular mechanisms may affect any of the transition rates. For example, synonymous mutations/changes (in endogenous or heterologous) genes inside the coding region may modify the adaptation of codons to the tRNA pool (codons that are recognized by tRNAs with higher intracellular abundance usually tend to be translated more quickly [34]); the local folding of the mRNA (stronger folding tends to decrease elongation rates [35]); the interaction/hybridization between the ribosomal RNA and the mRNA [36] (there are nucleotide sub-sequence that tend to interact with the ribosomal RNA, causing transient pausing of the ribosome and thus reducing the translation elongation rate). Non-synonymous mutations or changes inside the coding region may also affect the elongation rates, for example via the interaction between the nascent peptide and the exit tunnel of the ribosome [37,38]. In addition, intracellular changes in various translation factors (e.g. tRNA levels, translation elongation factors, concentrations of amino acids, concentrations of aminoacyl tRNA synthetase) and, as explained above, the mRNA levels can also affect elongation rates. Furthermore, various recent studies have demonstrated that manipulating the codons of a heterologous gene tend to result in significant changes in the translation rates and protein levels of the gene [2,40,56]. Because of the competition for the finite pool of ribosomes a change in one molecule will also affect the decoding of other mRNA molecules. The next subsection describes the main result in the paper, namely, the local and global effect of mutations in one mRNA on the others resulting from the competition for the finite resources.

### 4.2. Competition

A natural question is how will a change in the parameters (that is, the transition rates) of the RFMNP affect the equilibrium point of the network (and, in particular, the steady-state production rates). For example, if we increase some transition rate in RFMIO , how will this affect the steady-state production rate in the other RFMIOs? Of course, this effect is due to the competition for the finite resources. Without loss of generality, we assume that the change is in a transition rate of RFMIO #1.

### Theorem 4.3.

*Consider an RFMNP with m RFMIOs with dimensions* . *Let* *denote the set of all parameters of the RFMNP, and let*
*denote the equilibrium point of the RFMNP on some fixed level set of H. Pick* *. Consider the RFMNP obtained by modifying* *to* , *with* . *Let* *denote the equilibrium point in the new RFMNP and let* *. Then*4.44.54.6(In the case , condition (4.4) is vacuous.)

Increasing means that ribosomes flow ‘more easily’ from site *i* to site in RFMIO #1. Equation (4.4) means that the effect on the density in this RFMIO is that the number of ribosomes in site *i* decreases, whereas (4.5) implies that the number of ribosomes in all the sites to the right of site *i* increases. Equation (4.6) describes the effect on the steady-state densities in all the other RFMIOs and the pool: either *all* these steady-state values increase or they all decrease. The first case agrees with the results in example 3.3 above.

Note that the theorem does not provide any information on the change in , . Our simulations show that any of these values may either increase or decrease, with the outcome depending on the various parameter values. Thus, the amount of information provided by (4.5) depends on *i*. In particular, when is changed to then the information provided by (4.5) is only that

Much more information is available when .

### Corollary 4.4.

*Suppose that* *is changed to* *. Then*4.7*and*4.8

Indeed, for , (4.5) yields (4.7). Also, we know that the changes in the densities in all other RFMIOs and the pool have the same sign. This sign cannot be positive, as combining this with (4.7) contradicts the conservation of ribosomes, so (4.8) follows.

In other words, increasing , the intrinsic initiation rate in RFMIO #1, yields an increase in *all* the densities in RFMIO #1, and a decrease in all the densities in all the other RFMIOs. This makes sense, as increasing means that it is easier for ribosomes to bind to the mRNA molecule. This increases the total number of ribosomes along this molecule and, by competition, decreases all the densities in the other molecules and the pool. Note that this special case agrees well with the results described in [11] (see (1.1)). Thus, our model allows studying fundamental biological phenomena that are not covered by models that ignore the elongation dynamics.

### Example 4.5.

Consider the RFMNP in (3.7) with , , , and initial condition . We consider a range of values for . For each fixed value, we simulated the dynamics until steady state for two cases: and . Figure 11 depicts for the various fixed values of . It may be seen that we always have and . Figure 12 depicts , , and for the various fixed values of . It may be seen that for a small value of all the 's and are negative, whereas for large values of they all become positive. ▪

Intuitively, this can be explained as follows. When is small it is the bottleneck rate in RFMIO #1, and increasing only generates more ‘traffic jams’ along RFMIO #1. This depletes the pool, and thus decreases the production rate in the second RFMIO. On the other hand, when is large becomes the bottleneck rate, and increasing it to allows ribosomes to traverse RFMIO #1 more quickly. This decreases the number of ribosomes on this RFMIO, increases the pool occupancy, and this leads to a higher initiation rate in the second RFMIO.

Theorem 4.3 implies that when the codons of a gene are modified into ‘faster’ codons (either via synthetic engineering or during evolution) then the translation rates of the other genes either all increase or all decrease. However, theorem 4.3 does not provide information on when each of these two cases actually happens. In order to address this, we need to calculate derivatives of the equilibrium point coordinates with respect to the rates. The next result shows that these derivatives are well defined. Denote the mapping from the parameters to the unique equilibrium point in by , that is, , , . (Recall that is the vector of all the parameters in the RFMIOs.)

### Proposition 4.6.

*The derivative* *exists for all* .

The next example uses these derivatives to obtain information on the two cases that can take place as we change one of the rates.

### Example 4.7.

Consider an RFMNP with RFMIOs with lengths *n* and . To simplify the notation, let denote the equilibrium point of RFMIO #1 [RFMIO #2], and let , , denote the rates along RFMIO #1. Suppose that is changed to . Differentiating the steady-state equations

w.r.t. yields where we use the notation . These two equations yield

Recall that , , for all *j*, and for all *p*. Also, by theorem 4.3, , , for all , and , for all *k*. Thus,4.9

This means that the sign of the change in the densities in all the other RFMIOs and the pool depends on several steady-state quantities including terms related to the initiation rate and exit rate in RFMIO #1, and also the change in the total density in this RFMIO.

In the particular case (i.e. a very short RFM), equation (4.9) becomes4.10

Note that [] is the steady-state initiation [exit] rate in RFMIO #1. Thus, means that it is ‘easier’ for ribosomes to exit than to enter RFMIO #1, and in this case (4.10) means that when is increased the change in all other densities will be positive. This is intuitive, as more ribosomes will exit the modified molecule and this will improve the production rates in the other molecules. On the other hand, if , then it is ‘easier’ for ribosomes to enter than to exit RFMIO #1, so increasing will lead to an increased number of ribosomes in RFMIO #1 and, by competition, to a decrease in the production rate in all the other RFMIOs. ▪

Note that in the example above, increasing always increases the steady-state production rate in RFM #1 (recall that ). One may expect that this will always lead to an increase in the production rate in the second RFMIO as well. However, the behaviour in the RFMNP is more complicated because the shared pool generates a feedback connection between the RFMIOs in the network. In particular, the effect on the other RFMIOs depends not only on the modified production rate of RFMIO #1, but also on other factors including the change in the *total* ribosome density in RFMIO #1 (see (4.9)).

This analysis of a very short RFM suggests that the steady-state initiation rate of the mRNA with the modified codon plays an important role in determining the effect of modifications in the network. If this initiation rate is relatively low (so it becomes the rate limiting factor), as believed to be the case in most endogenous genes [57], then the increase in the rate of one codon of the mRNA increases the translation rate in all the other mRNAs, whereas when this initiation rate is high then the opposite effect is obtained. This latter case may occur, for example, when a heterologous gene is highly expressed and thus ‘consumes’ some of the available elongation/termination factors making the elongation rates the rate limiting factors.

To recap, figure 13 describes the biological implications of our analysis. Making a codon ‘faster’ has two possible global outcomes. Figure 13*a* depicts the case of *increasing* the initiation rate of an mRNA molecule: the translation rate in the modified mRNA and also the number of ribosomes on the mRNA will increase; due to the increased ribosome density on the mRNA there are less free ribosomes in the pool, and thus the initiation and translation rates of all other mRNA molecules *decreases*. Figure 13*b* describes the second possible outcome. Improving a slow codon, that is a bottleneck rate in an mRNA, alleviates the ‘traffic jam’ along this mRNA, so the number of free ribosomes in the pool increases. This improves the translation rate of all other the mRNAs.

We now describe other mathematical properties of the RFMNP.

### 4.3. Persistence

The next result shows that for any initial condition, the following property holds. After an arbitrarily short time all the densities are larger than zero and smaller than one, and the pool occupancy is strictly positive.

### Proposition 4.8.

*For any* *there exists* , *with* *when* , *such that for all* , *all* , *all* *and all* ,
*and*

In other words, after any time the solution is -separated from the boundary of *Ω* . This result is useful because on the boundary of *Ω*, denoted , the RFMNP looses some desirable properties. For example, its Jacobian matrix may become reducible on . Proposition 4.8 allows us to overcome this technical difficulty, as it implies that any trajectory is separated from the boundary after an arbitrarily short time.

### 4.4. Strong monotonicity

Recall that a cone defines a partial order in as follows. For two vectors , we write if ; if and ; and if . A dynamical system is called *monotone* if implies that for all . In other words, monotonicity means that the flow preserves the partial ordering [58]. 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 *i*, and if for all *i*. A system that is monotone with respect to this partial order is called *cooperative*.

The next result analyses the cooperativity of the RFMNP. Let denote the dimension of the RFMNP.

### Proposition 4.9.

*For any* *with* *,*4.11

*Furthermore, if* *then*4.12

This means the following. Consider the RFMNP initiated with two initial conditions such that the ribosomal densities in every site and the pool corresponding to the first initial condition are smaller or equal to the densities in the second initial condition. Then this correspondence between the densities remains true for all time .

### 4.5. Contraction

Contraction theory is a powerful tool for analysing nonlinear dynamical systems (e.g. [59]), with applications to many models from systems biology [60–62]. In a contractive system, the distance between any two trajectories decreases at an exponential rate. It is clear that the RFMNP is not a contractive system on *Ω*, with respect to any norm, as it admits more than a single equilibrium point. Nevertheless, the next result shows that the RFMNP is *non-expanding* with respect to the norm defined by .

### Proposition 4.10.

*For any* ,4.13

In other words, the distance between trajectories can never increase.

Pick , and let . Substituting in (4.13) yields4.14

This means that the convergence to the equilibrium point is monotone in the sense that the distance to can never increase.

### 4.6. Entrainment

Many important biological processes are periodic. Examples include circadian clocks and the cell-cycle division process. Proper functioning requires certain biological systems to follow these periodic patterns, i.e. to *entrain* to the periodic excitation.

In the context of translation, it has been shown that both the RFM [26] and the RFMR [30] entrain to periodic translation rates, i.e. if all the transition rates are periodic time-varying functions, with a common (minimal) period then each state variable converges to a periodic trajectory, with a period *T*. Here we show that the same property holds for the RFMNP.

We say that a function *f* is *T*-periodic if for all *t*. Assume that the 's in the RFMNP are time-varying functions satisfying

— there exist such that for all and all , .

— there exists a (minimal) such that all the 's are

*T*-periodic.

We refer to the model in this case as the *periodic ribosome flow model network with a pool* (PRFMNP).

### Theorem 4.11.

*Consider the PRFMNP. Fix an arbitrary* *. There exists a unique function* , *that is T-periodic, and for any* *the solution of the PRFMNP converges to* .

In other words, every level set of *H* contains a unique periodic solution, and every solution of the PRFMNP emanating from converges to this solution. Thus, the PRFMNP entrains (or phase locks) to the periodic excitation in 's. This implies in particular that all the protein production rates converge to a periodic pattern with period *T*.

Note that since a constant function is a periodic function for any *T*, theorem 4.11 implies entrainment to a periodic trajectory in the particular case where one of 's oscillates, and all the other are constant. Note also that the stability result in theorem 4.1 follows from theorem 4.11.

### Example 4.12.

Consider the RFMNP (3.7) with , and all rates equal to one except for . In other words, there is a single time-varying periodic rate in RFMIO #2. Note that all these rates are periodic with a common minimal period . Figure 14 depicts the solution of this PRFMNP as a function of time *t* for . The initial condition is for all . It may be seen that all the state variables converge to a periodic solution. In particular, all state variables converge to a periodic solution with (minimal) period , and so does the pool occupancy . 's also converge to a periodic solution, but it is not possible to tell from the figure whether there are small oscillations with period or the convergence is to a constant (of course, in both cases this is a periodic solution with period ). However, it can be shown using the first two equations in (3.7) that if converges to a periodic solution then so do and . Note that the peaks in are correlated with dips in , this is because when is high on some time interval, i.e. the transition rate from site 2 to site 3 is high, there is a high flow of ribosomes from site 2 to site 3 during this interval. ▪

From the biophysical point of view, this means that the competition for free ribosomes induces a coupling between the mRNA molecules. This can induce periodic oscillations in *all* the protein production rates even when all the transition rates in the molecules are constant, except for a single rate in a single molecule that oscillates periodically. The translation rate of codons is affected among others by the tRNA supply (i.e. the intracellular abundance of the different tRNA species) and demand (i.e. total number of codons from each type on all the mRNA molecules) [63]. Thus, the translation rate of a codon(s) is affected by changes in the demand (e.g. oscillations in mRNA levels) or by changes in the supply (e.g. oscillations in tRNA levels). The results reported here may suggest that oscillations in the mRNA levels of some genes or in the concentration of some tRNA species (that occur, for example, during the cell cycle [64,65]), can induce oscillations in the translation rates of the rest of the genes.

## 5. Discussion

We introduced a new model, the RFMNP, for large-scale simultaneous translation and competition for ribosomes that combines several RFMIOs interconnected via a dynamic pool of free ribosomes. To the best of our knowledge, this is the first model of a network composed of interconnected RFMIOs. The RFMNP is amenable to analysis because it is a monotone dynamical system that admits a non-trivial first integral. The fact that the total number of ribosomes in the network is conserved means that local properties of any mRNA molecule (e.g. the abundance of the corresponding tRNA molecules) affects its own translation rate, and via competition, also *globally* affects the translation rates of all the other mRNAs in the network. We provide several examples demonstrating the behaviour of the RFMNP, and rigorously analyse some of its mathematical properties. In particular, the RFMNP is an irreducible cooperative dynamical system admitting a continuum of linearly ordered equilibrium points, and every trajectory converges to an equilibrium point. The RFMNP is also on the ‘verge of contraction’ with respect to the norm, and it entrains to periodic transition rates with a common period.

An important implication of our analysis and simulation results is that there are regimes and parameter values where there is a strong coupling between the different ‘translation components’ (ribosomes and mRNAs) in the cell. Such regimes cannot be studied using models for translation of a single isolated mRNA molecule. The RFMNP is specifically important when studying highly expressed genes with many mRNA molecules and ribosomes translating them because the dynamics of such genes strongly affects the ribosomal pool. For example, changes in the translation dynamics of a heterologous gene which is expressed with a very strong promoter, resulting in very high mRNA copy number should affect the entire tRNA pool, and thus the translation of other endogenous genes. Highly expressed endogenous genes ‘consume’ many ribosomes. Thus, a mutation that affects their (local) translation rate is expected to affect also the translation dynamics of other mRNA molecules. Studying the evolution of such genes should be based on understanding the global effect of such mutations using a computational model such as the RFMNP.

On the other hand, we can approximate the dynamics of genes that are not highly expressed (e.g. a gene with mRNA levels that are 0.01% of the mRNA levels in the cell) using a single RFM. In this case, the relative effect of the mRNA on all other mRNAs is expected to be limited.

Our analysis shows that increasing the translation initiation rate of a heterologous gene will always have a negative effect on the translation rate of the other genes (i.e. their translation rates decrease) and vice versa. The effect of increasing [decreasing] the translation rate of a codon of the heterologous gene on the translation rate of other genes is more complicated: while it always increases [decreases] the translation rate of the heterologous gene it may either increase or decrease the translation rate of all other genes. The specific outcome of such a manipulation can be predicted using the RFMNP with parameter values that are based on the biophysical properties of the heterologous genes and the host genome.

Some of our analytical results may be tested experimentally. This can be done by designing and expressing a library of heterologous reporter genes [2,40,56]. The initiation rates can be manipulated based on the engineering of the nucleotides surrounding the START codon (e.g. [40–42]). The elongation rates of codons (and thus the transition rates) can be estimated using ribosome profiling [66] in *in vivo* experiments [34,40,67]. We can simulate a model such as the RFM (or TASEP) based on the inferred initiation and elongation rates (e.g. [23]), and use the models described here to design and demonstrate examples of heterologous genes where increasing [decreasing] the rate of a codon increases [decreases] the initiation (and translation) rate of all other mRNAs.

The effect of the manipulation of a codon (i.e. increasing or decreasing its rate) of the heterologous reporter gene on the ribosomal densities and translation rates of all the mRNAs (endogenous and heterologous) can be performed via ribosome profiling [66] in addition to measurements of mRNA levels, translation rates and protein levels [68]. One can also measure the fluorescence level of the protein related to an *additional* reporter gene following the expression of our library, as was recently done in [69]. Variants in the library that decrease/increase the free ribosomal pool are expected to decrease/increase the fluoresce level of the protein related to the additional reporter gene. Furthermore, we can measure the growth rate of the host and examine whether indeed a global increase [decrease] in the initiation rates yields a positive [negative] effect on the host growth rate.

Our analysis suggests that the effect of improving the transition rate of a codon in an mRNA molecule on the production rate of other genes and the pool of ribosomes depends on the initiation rate in the modified mRNA. When the initiation rate is very low the effect is expected to be positive (all other production rates increase). However, if the initiation rate is high the effect may be negative. This may partially explain the selection for slower codons in highly expressed genes that practically decrease the initiation rate [1,4]. This relation may also suggest a new factor that contributes to the evolution of highly expressed genes towards higher elongation and termination rates (i.e. the tendency of highly expressed genes to include ‘fast’ codons). Indeed, lower elongation rates (and thus a relatively high initiation rate) may decrease the production rates of other mRNAs that are needed for proper functioning of the organism.

The RFM, and thus also the model described here, does not capture certain aspects of mRNA translation. For example, eukaryotic ribosomes may translate mRNAs in multiple cycles before entering the free ribosomal pool [18,44,70]. This phenomenon may perhaps be modelled by adding positive feedback [44] in the RFMNP. In addition, different genes are transcribed at different rates, resulting in a different number of (identical) mRNA copies for different genes. This can be modelled using a set of identical RFMs for each gene. Such a model can help in understanding how changes in mRNA levels of one gene affect the translation rates of all the mRNAs. The analysis here suggests that modifying the mRNA levels of a gene will affect the translation rates of all other genes in the same way. These and other aspects of biological translation may be integrated in our model in future studies. Hamadeh & Del Vecchio [71] develop the notion of the realizable region for steady-state gene expression under resource limitations, and methods for mitigating the effects of ribosome competition. Another interesting research direction is studying these topics in the context of the RFMNP.

We believe that networks of interconnected RFMIOs may also prove to be powerful modelling and analysis tools for other natural and artificial systems. These include communication networks, intracellular trafficking in the cell, coordination of large groups of organisms (e.g. ants), traffic control and more.

## Competing interests

We declare we have no competing interests.

## Funding

The research of M.M. and T.T. is partially supported by a research grant from the Israeli Ministry of Science, Technology, and Space. The work of E.D.S. is supported in part by grants AFOSR FA9550-14-1-0060 and ONR 5710003367.

## Acknowledgement

We thank Yoram Zarai for helpful comments. We are grateful to the anonymous reviewers for comments that helped us to improve this paper.

- Received December 12, 2015.
- Accepted February 15, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.