A Max-Plus Model of Ribosome Dynamics During mRNA Translation

We examine the dynamics of the translation stage of cellular protein production, in which ribosomes move uni-directionally along mRNA strands building an amino acid chain as they go. We describe the system using a timed event graph - a class of Petri net useful for studying discrete events which take a finite time. We use max-plus algebra to describe a deterministic version of the model, calculating the protein production rate and density of ribosomes on the mRNA. We find exact agreement between these analytical results and numerical simulations of the deterministic case.


Introduction
Messenger RNA (mRNA) translation is one of the steps in protein production in cells [1]. mRNAs are single strands of nucleotides which are transcribed from the DNA. The sequence of nucleotides, grouped in triplets called codons, holds the code for a specific chain of amino acids that makes up a protein.
Translation is performed by molecular machines called ribosomes, which scan along the mRNA adding amino acids to a growing chain which will become the protein. The rates at which different proteins are produced are crucial in determining how a cell grows and functions.
Various statistical models have been used to describe and understand the translation process. In this paper we propose a new model which can be analysed rather completely using algebra on the max-plus semi-ring. The most convenient way to present the model is to write it as a deterministic Petri net [2]. This can then be analysed completely as a linear dynamical system when written in terms of the max-plus algebra [3]. The model gives a more realistic representation of the movement of ribosomes than that of many previous studies, and explicitly considers the exact genetic coding sequence as well as the fact that different codons are translated at different rates. We perform numerical simulations of the Petri net, and find the results to exactly match the behaviour predicted by the max-plus algebra.
Although the general sequence of events in translation is understood, there are many questions which remain unanswered. Translation begins when a ribosome binds to the end of an mRNA strand -a process known as initiation. In this paper we focus on the second stage of translation, which is known as elongation. Here the ribosome moves uni-directionally along the mRNA, pausing at each codon to recruit an amino acid, which is then added to the growing chain. Amino acids are transported within the cytoplasm by carrier molecules called transfer RNAs (tRNAs). Different species of tRNA carry different amino acids, and different codons correspond to different tRNAs. The different tRNAs appear in different concentrations, so the time which a ribosome waits for the required tRNA to arrive differs from codon to codon [4,5]. Since it is often the case that several ribosomes are bound to the mRNA at the same time, the distribution of waiting times can result in traffic jams if the progress of one ribosome is obstructed by another. This picture is complicated by the fact that whilst there are only 20 common amino acids, in the model organism Saccharomyces Cerevisiae there are 41 species of tRNA, i.e., there is redundancy in the genetic code. For some amino acids there are multiple tRNAs; furthermore it is often that case that there is a highly abundant tRNA and a rare tRNA which code for the same amino acid. Sometimes a "slow codon", which will cause a long pause in elongation, is used when a fast codon corresponding to the same amino acid is available. The interest for biologists is therefore to understand how slow codons effect protein production rates and the use of ribosomes [6]. Translation reaches its completion in the termination stage where, through the binding of release factors, the ribosome disassociates from the mRNA and the amino acid chain is released ready for folding or further processing. The quantities of interest which we shall take from the model are the time interval between successive amino acid chains being completed, which we will call the protein production time, and the ribosome occupation density of the mRNA.
In the next section we introduce Petri nets, and detail how these can be used to describe translation. Then in Sec. 3 we introduce the max plus algebra, giving a brief survey of the salient facts required in the rest of the paper. In Sec. 4 we detail the max-plus treatment of the Petri net describing translation. These results are then compared to numerical simulations: we first consider simple "designer mRNAs" where each codon corresponds to the same tRNA, before treating realistic sequences taken from the S. Cerevisiae genome. Finally we consider how introducing stochasticity into the model is likely to alter the results, and how the present work relates to previous models of translation (specifically the totally asymmetric simple exclusion process or TASEP).

Timed Petri Nets
Petri nets (attributed to C. A. Petri [7]) are a scheme where a sequence of discrete events is described on a network. We give a brief introduction here, but for a detailed description refer the reader to Ref. [2] and references therein.
A Petri net is a directed graph which contains two types of nodes: places and transitions. These are shown diagrammatically in Fig. 1(a). Directed arcs connect places with transitions, but not places with places or transitions with transitions. Places can contain objects called tokens. A transition is said to be active when each of the upstream places it is connected to contains at least one token. Figure 1(b) shows an example Petri net. When a transition is active it can fire; upon firing one token is removed from all places upstream of the transition, and one token is added to all places downstream of the transition (clearly there is no implicit conservation of tokens since the number of upstream places need not equal the number of downstream places). Events unfold in a discrete manner. Petri nets are often used to model systems where events occur (transitions fire) given that a set of conditions are fulfilled (tokens are present).
A timed Petri net is an extension to this framework in which a waiting time is attached to each place. Whenever a token is put into a place, it will only contribute to the activation of transitions once a time τ associated with that place has elapsed. This can be thought of as a timer on the token which is started as it enters a place. As soon as a transition becomes active it fires. In this way the Petri net describes discrete events occurring in continuous time as conditions are fulfilled. Unlike some other models (for example Monte Carlo based simulations) events can occur simultaneously.
We can describe the movement of ribosomes along the mRNA as a timed Petri net, and this is shown in Fig. 2(a). A sequence of codons of length n is represented by two rows of places, with each pair (upper and lower) representing a codon. A token in the top row indicates that a ribosome is decoding that codon, and a token in the bottom row indicates a vacant codon. Having this two row structure ensures that there can only be one token in each place; i.e., one ribosome translating each codon. This of course requires that the initial state also satisfies this condition. A suitable initial state is having all top row places empty and all bottom row places containing one token; this corresponds to a mRNA initially free of ribosomes.
The firing of a transition T i represents the movement of a ribosome from codon i to codon i+1. There are n + 1 transitions, labelled i = 0, . . . n, with the firing of the zeroth and nth transitions representing initiation and termination events respectively. At the leftmost end of the Petri net we have only an upper row place; this represents a continuous supply of ribosomes ready to begin translation. At the rightmost end of the Petri net there is only a lower row place, representing the cytoplasm which is always ready to accept ribosomes. These places are kept full with a token via an arc to and from the same transition.
The time associated with each upper row place corresponds to the time it takes a ribosome to capture a tRNA of the species corresponding to that codon. Generally the arrival times of tRNAs can be represented by a set of independent Poisson processes, however a deterministic version of the model can be obtained by replacing these random times with their means. In the present paper we concentrate on the deterministic version; the consequences of introducing stochasticity are discussed in Sec. 6. We assume that the mean codon waiting time is proportional to the concentration of the species of tRNA corresponding to that codon. We denote the waiting time for the ith codon counting from the left τ i . We assume that when a ribosome moves forward the codon it was covering is immediately vacant, and the time it takes for the ribosome to physically move is negligible; thus we associate a time zero with each of the lower row of places. In reality a ribosome actually covers more than just the codon it is translating, but for simplicity we do not take this into account here. The first place in the chain -corresponding to initiation -has a time τ in associated with it. Biologically this corresponds to the average time interval between attempts for a ribosome to begin translation, and will depend on the availability of ribosomes and various initiation factors, as well as the presence of any secondary structures in the non-coding leader region of the mRNA. Since our focus here is on elongation, we assume that τ in is constant in time, treating it as a control parameter. The final place in the chain has an associated time τ out which is the time a ribosome waits at the end of the mRNA before releasing the completed protein and detaching from the mRNA. Biologically this could be linked to the availability of a number of release factors, and again we treat it as a control parameter. The model is therefore parametrised by the initiation and termination waiting times τ in and τ out , and the set of internal tRNA capture waiting times {τ i |1 ≤ i ≤ n}. Figure 2(b) shows and example Petri net for a very (unrealistically) short mRNA of length n = 4 codons, with ribosomes occupying codons 2 and 3. Transition T 2 cannot fire because there is no token in its input place in the lower row (representing a vacancy). Transition T 3 will fire as soon as the token in its upper row input place has been there for a time τ 3 . If the token in the place upstream of T 2 (upper row place for codon 2) has been there for a time longer than τ 2 when T 3 fires, then T 2 will also fire at this time. Thus if there is a slow codon (large τ i ), when the corresponding transition fires, it could result in several other transitions simultaneously firing.

The max-plus semi-ring
If each place in a Petri net has exactly one upstream and one downstream transition (as is the case in those depicted in Fig. 2), then its dynamics can be described using max-plus algebra. Max-plus is an algebra over the semi-ring R max ∈ R ∪ −∞. In this section we give a brief overview of the theorems required in the rest of the paper, and refer the reader to Refs. [3,8] for further details and proofs.
In max-plus algebra the addition and multiplication operators ⊕ and ⊗ are defined as where a, b ∈ R max , and we define ǫ = −∞ and e = 0 which behave as zero and unity respectively. Max-plus algebra refers to the set R max = (R max , ⊗, ⊕, ǫ, e). The elements of R max have much the same properties as in conventional algebra; for example We define Max plus powers in the natural way a ⊗x = a ⊗ a ⊗ · · · ⊗ a x times, and note that in conventional algebra this corresponds to x × a. We note that the ⊗ operator has an inverse which can be expressed as a negative power but there is no inverse of the ⊕ operator. We also use the following notation for sums and products over indices: Vectors and matrices are also defined, and we denote by [A] ij the i-jth component of the matrix A. The sum and product of matrices A ∈ R n×l max and B ∈ R l×m max are then defined as Matrices can have associated eigenvectors and eigenvalues; i.e. a vector u satisfying is an eigenvector of A associated with eigenvalue λ. Note that eigenvectors are not unique, i.e. if u is an eigenvector then so is α ⊗ u with α an arbitrary finite number. In general a matrix can have multiple associated eigenvalues.
The properties of a max-plus matrix can often be determined by considering the associated directed weighted graph. The graph associated with matrix A, denoted G(A), consists of a set of nodes and directed weighted arcs, where if there is a matrix element [A] ij with value a ij = ǫ, then there is an arc from node j → i with weight a ij (note that direction of the arc differs to that in conventional algebra and graph theory). A series of of one or more arcs between two nodes i and j is called a path from i → j, and if there is a path i → i this is called a circuit and is denoted γ. A graph is said to be strongly connected if for any two different nodes there is a path between them, and a matrix A is said to be irreducible if the corresponding graph G(A) is strongly connected. The circuit weight w γ of a circuit γ is defined as the sum of the weights of all arcs in that circuit, and the circuit length l γ as the number of arcs in the circuit. The mean circuit weight is definedw γ = w γ /l γ . If the maximum mean circuit weight in a graph is λ, then a circuit with a mean circuit weight equal to λ is called a critical circuit. The critical graph corresponding to matrix A, denoted G c (A), is defined as the sub-graph of G(A) containing only the nodes and arcs which are in the critical circuits. The cyclicity of a graph σ G is defined as the greatest common divisor of the lengths of all of the circuits in that graph, and the cyclicity of a matrix A, σ A is equal to the cyclicity of the critical graph of A, G c (A).

Theorem 1 Any irreducible matrix A ∈ R n×n max possesses one and only one eigenvalue λ, which is a finite number and is equal to the maximal mean circuit weight of circuits in the graph G(A). (See, for example, Theorem 2.9 in [3] for proof.)
Thus the eigenvalue of an irreducible matrix can be found by considering the corresponding graph. As noted above there is a whole continuum of eigenvectors associated with the eigenvalue of a max-plus matrix. Two eigenvectors x, y are co-linear if there exists a scalar α ∈ R max such that y = α ⊗ x, and such vectors can be projected onto the same object in a projective space [3]. A given eigenvalue of the matrix A can be associated with more that one linearly independent eigenvector, and the number of such eigenvectors can be found by considering the critical graph G c (A).
Theorem 2 If the critical graph G c (A) of an irreducible matrix A has k maximal strongly connected subgraphs (m.s.c.s.), then A has k linearly independent eigenvectors. (For proof see, for example, Theorem 4 in [9].) Linear equations such as where A ∈ R n×n max and x, b ∈ R n max , can often be solved using an object called the Kleene star, defined as The existence of A * can be proven if the graph G(A) has only non-positive circuit weights [3]. We shall see below that the problem studied in this paper involves sequences of vectors x(k), k ∈ N described by a recurrence equation for k ≥ 0 where A ∈ R n×n max is irreducible, and x(0) is some initialisation vector.

Theorem 4 The cyclicity theorem states that there is an integer K 0 such that
where λ and σ A are the eigenvalue and cyclicity of the matrix A respectively. (For proof see [9] or Theorem 3.9 in [3].) This implies some periodicity in the powers of A; the periodic behaviour is characterised by λ and σ A , and K 0 is known as the transient time of A. By using this in Eq. (3) we find If matrix A has cyclicity σ A = 1, then we have i.e., x(k) (or equivalently A ⊗k ⊗ x(0) for any x(0) ∈ R n max ) is an eigenvector of A for k ≥ K 0 . The effect of the initial condition has died out.

Analysis of translation with the max-plus algebra
A Petri net can be described using max plus algebra by defining the matrices A 0 , A 1 , . . . , A M which depend on the initial conditions such that [A m ] ij is equal to the maximum of the holding times associated with the places between transitions j and i, which initially contain m tokens. The dynamics are then described by the vector x(k), where the ith component of the vector is equal to the time at which the ith transition fires for the kth time. This vector satisfies the recursion equation For further details see [3,10,11].
In the case of the Petri net which describes mRNA translation (shown in Fig 2(a)), the structure and initial conditions are such that each place can contain either 0 or 1 tokens, i.e., in this case M = 1 and the above equation reduces to In the Petri net we have n + 1 transitions labelled 0, . . . n, so x(k) has components {x i |0 ≤ i ≤ n}. Likewise A 0 and A 1 are (n + 1) × (n + 1) matrices, with elements labelled by indices running from 0 to 1. These matrices are Written in component form Eq. (5) gives If we think of the integer k as labelling the ribosomes (x i (k) gives the time at which the kth ribosome leaves the ith codon) then these equations make sense conceptually since a ribosome will leave codon i either a time τ i after it left codon i − 1, or at the time when the (k − 1)th ribosome leaves codon i (i.e., time x i+1 (k − 1)), whichever is latest.
Although the solution of the recursion equation (5) is not straightforward, since there is no inverse to the ⊕ operator, by identifying the second term as a vector b = A 1 ⊗ x(k − 1), we note that this equation is of the form of Eq. (1), and so via Theorem 3 the solution is From that theorem we know that A * 0 exists since the graph G(A 0 ) contains no circuits of positive weight; this can also be easily demonstrated by considering the first few powers of A 0 [3]. For example for a system with four transitions (n = 3) we see i.e., the number of non-zero entries decreases as the power increases, and actually The Petri net describing translation can therefore be described by the equation where B = A * 0 ⊗ A 1 is given by The eigenvalue and cyclicity of the matrix B can be found from the corresponding graph G(B). In Fig. 3 we show G(B) for the example of n = 3. We note that the graph is strongly connected, so B is irreducible and has exactly one eigenvalue (Theorem 1). Due to non-zero elements on the diagonal of B, there exist circuits of length one, so the cyclicity of the graph is σ G = 1. By inspection of the graph it is clear that one of the circuits of length one will always have a mean circuit weightw γ equal to the maximum mean circuit weight. That is to say, the maximum mean circuit weight, and therefore the eigenvalue of B will be equal to whichever of the waiting times is largest, i.e.,  7)) with n = 3.
It is also clear then that the critical graph of B must always contain a circuit of length one, so σ B = 1.
Due to the cyclicity theorem (Theorem 4), we conclude that for large k the vector x(k) is an eigenvector of B. Regardless of the initial condition x(0), the Petri net describing translation always reaches a steady state described by gives the time at which the ith transition fires for the kth time. The nth component of this vector equation tells us that the transition corresponding to completion of an amino acid chain will fire at time intervals of λ; thus the eigenvalue of A can be interpreted as the protein production time. As noted above, another way of interpreting the eigenvector is that x i (k) is the time at which the kth token (ribosome) leaves the ith place (codon); the token enters the ith place at time x i−1 (k) and the maximum length of time which a ribosome can stay at a given codon is λ. The proportion of time for which there is a token occupying the ith place in the steady state, i.e. the occupation density, is therefore given by We note that ρ i are the components of a vector of length n. It now remains to evaluate the eigenvector. If we define u = x(k − 1), then from Eq. (8), x(k) = λ⊗ u. We can then use this to evaluate the eigenvalue using Eq. (5) which when written as components gives Thus the eigenvectors takes a different form depending on the eigenvalue. We identify the following three cases: We first consider the case where the parameters are such that the waiting time for initiation of translation is longer than that for tRNA capture and termination, i.e., τ in > τ out , τ i for 1 ≤ i ≤ n. We shall denote this the entry limited regime, as it represents the situation where entry of the ribosomes into the system is the rate limiting process. Replacing λ → τ in in Eqs. (10)(11)(12) gives From (15), since τ in = τ out , for consistency it must be the case that τ in ⊗ τ n ⊗ u n−1 > τ out ⊗ u n , and therefore u n = τ n ⊗ u n−1 .
From (14), taking i = n − 1 gives Since τ in = τ n , in order for this to be consistent with Eq. (16) it must be the case that τ in ⊗τ j ⊗u n−2 > u n , meaning Continuing to use Eq. (14) in this way we can find expressions for each component of u in terms of the nth component. The eigenvector is given by Using Eq. (9) we find ρ i = τ i /τ in .
In the case where τ out > τ in , τ i for 1 ≤ i ≤ n, the time it takes a ribosome to leave the mRNA is the limiting process, i.e., this is the exit limited regime. Following the same method as above we make the replacement λ → τ out in Eqs. (12) giving This time we start with the first of the three equations, and note that since τ in = τ out , it must be the case that u 1 > τ in ⊗ u 0 ; therefore We then consider Eq. (18); for i = 1 this gives In order for this to be consistent with Eq. (20), it must be the case that u 2 > τ out ⊗ τ 1 ⊗ u 0 , meaning We continue using (18) to find the other components of the vector; the result is leading to occupation densities ρ i = 1.
The final case is the elongation limited regime, where one or more of the internal waiting times has a value τ , with τ > τ in , τ out , {τ i = τ }. We assume initially that there are two places with waiting times equal to τ deep in the bulk of the mRNA, denoting the first the pth, and the second the qth. That is to say τ p , τ q = τ where 1 < p < q < n. We again take Eqs. (12) and this time make the replacement λ → τ , giving We start with Eq. (21), which since τ in = τ , gives Moving onto Eq. (22) and taking i = 1, if τ 1 = τ it must be the case that We can continue using Eq. (22) in this way until we reach the pth waiting time, since τ p = τ . Taking i = p − 1 and then i = p in Eq. (22) we find respectively Since neither term on the right hand side of the second equation would contradict the first, all this can tell us is that We now consider Eq. (23); since τ out = τ it must be the case that Using Eq. (22) with i = n − 1 gives If τ n = τ then in order not to contradict Eq. (25), it must be the case that τ ⊗ τ n−1 ⊗ u n−2 ≥ u n . Therefore u n−1 = τ n−1 ⊗ u n−2 .
We can continue to iterate backwards using Eq. (22) until we reach the qth waiting time τ q = τ . Taking i = q and then i = q − 1 gives respectively Again the second equation can only give an inequality We are therefore left with two inequalities, Eqs. (24) and (26) which cannot alone tell us about the pthe and the qth components of the eigenvector. We find that if we assume that one of these inequalities is actually an equality, then we can continue using Eq. (22) to find all of the components of the eigenvector consistently with all of the above equations. This leads to two eigenvectors depending on which inequality we set equal where we recall that p < q and the pth and qth are the upstream and downstream most codons with waiting times equal to τ . Any max-plus linear combination either of these two vectors is an eigenvector of B. If we consider the critical graph of B, if there are M non-adjacent codons with waiting times equal to τ , then this graph will contain M m.s.c.s., and therefore (via Theorem 2) there are M linearly independent eigenvectors. We note that the critical graphs in cases (i) and (ii) have only one m.s.c.s.. Each of the two eigenvectors above gives rise to a different occupation density profile using Eq. (9), i.e., which corresponds to queues of ribosomes behind the qth and the pth codons respectively. In mRNAs with more than two isolated codons with waiting times equal to τ , there will be solutions with queues behind each of these "slowest" codons. We will show in the next section that only one of the solutions is realised, and we reserve further discussion until Sec. 5.

Simulation Results
In this section we compare the results from the max-plus algebra detailed above with direct numerical simulation of the Petri net shown in Fig. 2(a). We conduct the simulations by considering an initial state where all bottom row places contain a token and all top row places are empty (the first and last places in the chain both contain one token). We then move forward in time from one event (firing of a transition) to the next, updating token positions at each step. We allow the system to reach a steady state, and then examine the protein production time P , and the codon occupation density of each place ρ i , as well as the mean density, defined as We first consider a very simple mRNA where each codon is assumed to be identical, i.e., it codes for the same tRNA, and the waiting time for each is the same. The internal hopping rates are chosen to be τ i =τ for 1 ≤ i ≤ n, andτ is used as the unit of time. We then measure P , ρ i and ρ for different values of τ in and τ out .
In Sec. 5.2 we consider real mRNA sequences from the S. Cerevisiae genome. The waiting time for a particular tRNA species is assumed to be inversely proportional to the concentration of those molecules found in a typical cell. We estimate this from the gene copy number of each tRNA [12]. We choose the internal hopping rates such that the average of the waiting times is equal toτ , and again use this as the unit of time. Figure 4 shows colour maps for the protein production time and mean density as a function of τ in and τ out for a uniform mRNA of length n = 500, as generated from numerical simulations. There are three regimes, with the dynamics depending on the relative magnitudes ofτ , τ in , and τ out . These correspond to the three cases found from the max-plus algebra, i.e., the initiation, termination and tRNA capture limited regimes. We summarise the results in Table 1. Figure 5 shows how P and ρ vary as a function of (a) τ in for fixed τ out , and (b) τ out for fixed τ in . The numerical and analytic results match exactly.   Table 1: Details of each regime as given by the max-plus algebra, for a uniform mRNA, where the tRNA capture waiting times are τ i =τ for 1 ≤ i ≤ n. These are found to exactly match the numerical simulation results.

Heterogeneous mRNAs
We now consider mRNA sequences from genes YJL136C and YDR382W in S. Cerevisiae, which we here on denote mRNA A and mRNA B, respectively. For a given mRNA of length n codons, each of the waiting times {τ i |1 ≤ i ≤ n} takes a value chosen from the set of 41 times {s i |1 ≤ i ≤ 41} corresponding to each of the 41 tRNA/codon species, and (1/41) 41 i=1 s i =τ . We takeτ as the unit of time; biological experiments estimate thatτ ≈ 0.1 s −1 [1]. Figure 6 shows the waiting times at each codon for each of the sequences. Figures 7 and 8 show simulation results for P and ρ for each mRNA respectively. We again observe three regimes depending on the relative magnitude of τ in , τ out , and the largest of the τ i which we denote τ . Again the simulation results exactly match those from the max-plus algebra, which are summarised in Table 2.
In two of the three regimes the occupation density differs from codon to codon. In Fig. 9 we show the density profile for each mRNA when parameters are chosen such that the system is in the entry limited phase. We note that for most codons i, ρ i < 0.5; in other models (i.e., the TASEP which we discuss in Sec. 7) the equivalent regime is usually called the low density phase. In Fig. 10 we show the density profiles for the tRNA capture limited phase. For codons located upstream of the most upstream codon with the longest waiting time (τ p = τ ) the density ρ i ≈ 1 for i < p; i.e., the leftmost instance of the slowest codon causes queueing of ribosomes. In the case of mRNA B the first "slowest codon" is at i = 1, so no queue is observed. Although from the max-plus algebra there are other possible solutions, only one is realised in the simulations. Instances of the slowest codon appearing further downstream do not give rise to queues, since the rate at which ribosomes arrive at those codons is too low. In terms of the max plus algebra, the different possible eigenvectors are reachable by different initial conditions [13].
(iii) tRNA capture limited (τ > τin, τout) Protein production time P τin τout τ Occupation density ρi τi/τin Table 2: Details of each phase as given by the max-plus algebra, for a non-uniform mRNA. Here τ is the longest internal waiting time (τ = max{τ i |1 ≤ i ≤ n}), and p is the position of the leftmost codon with waiting time equal to τ .   and x i (0) = ǫ for i = 0; i.e. we start counting the elapsed time from 0 at the first initiation event. Only the eigenvectors corresponding to a queue behind the most upstream slowest codon are reachable from this initial condition.
The other solution for the density profile of mRNA B can be accessed by making a change to the waiting times during the simulation. If we first allow the system to reach steady state, and then increase the waiting time on the slow codon at position i = 104 by some small amount δτ , then a queue will begin to form behind this codon. If we reverse the change after the system reaches steady state, then the queue behind codon 104 will persist (data not shown).

Introduction of stochasticity
In the above work we have examined a deterministic version of the model, where rather than choosing waiting times from a Poisson distribution, we instead use the mean of the distribution. Simulation of a Petri net with waiting times chosen from a distribution is straightforward. The description using max-plus algebra is less so. At its essence the problem consists of a product of irreducible max-plus matrices such that x(k + 1) = A(k) ⊗ x(k) for k ≥ 0, where {A(k)|k ∈ N} is an independent and identically distributed (i.d.d.) sequence of random matrices. Some work on such sequences can be found in the literature [3,14], and it can be shown that for an i.d.d. sequence of matrices which have fixed support and are irreducible, a max-plus Lyapunov exponent exists, and this is equal to the asymptotic growth rate. However, there is currently no method for calculating this quantity, which in the present model would be related to the mean protein production time. For this reason, here we only briefly discuss how simulation results for a stochastic Petri net compare with the deterministic case.
We consider here a uniform mRNA of length n = 500 codons, where waiting times for initiation, termination and tRNA capture are chosen from exponential distributions with means τ in , τ out andτ respectively. We perform simulations where we allow the system to reach steady state before measuring the time averaged protein production time and codon occupation density, which we again denote P and ρ respectively. As shown in Figs. 11(a) and (d) we again see three regimes depending on which waiting time is the largest (τ in and τ out compared toτ ). We note by comparing with Fig. 4 that in the stochastic case the boundary for the tRNA capture limited regime is at τ in , τ out ≈ 2τ , compared to τ in , τ out =τ in the deterministic case. Also in that regime, the density is significantly reduced compared to the deterministic case (ρ ≈ 0.75 compared to ρ = 1) and the protein production time is increased (P ≈ 2.7τ compared to P =τ ). This is as expected since allowing the internal waiting times to vary will lead to gaps between ribosomes. The fact that the difference between the stochastic and deterministic models in this regime is so large can be attributed to the fact that the rate limiting process involves choosing waiting times for each of the n = 500 places representing the codons. In the other regimes the effect of the stochasticity is less severe, as the limiting process involves only one place. Deep within the initiation or termination limited regimes ρ and P are almost the same as in the deterministic model. The general effect of stochastic waiting times in these regimes is that both density and protein production time increase. If we consider ribosomes occupying two consecutive codons, and if the waiting time drawn for the leading ribosome is long, it could hold up the ribosome behind it; if the time drawn for the leading ribosome is short it is unlikely to effect the ribosome behind it. Thus we expect ribosomes will in general move more slowly along the mRNA; as we shall see below, often the maximum of two exponentially distributed random numbers is what determines the overall behaviour.
For non-uniform mRNAs the effect of stochasticity in the initiation and termination limited regimes are similar: both ρ and P are increased compared to the deterministic model. In the tRNA capture regime the situation is more complicated, and features such as multiple queues are observed. Discussion of such phenomena is beyond the scope of the present paper. Here we only briefly consider one consequence of stochasticity which has played an important role in other translation models [15,16]. It has been observed in a widely used alternative model for translation (the TASEP, as discussed in section 7) that several slow codons in close proximity have a more dramatic effect on the protein production rate from   an mRNA than slow codons in isolation. In the deterministic Petri net the token release time interval from a place corresponding to a slow codon (a slow place) which is in isolation is the same as that from the second of a pair of adjacent identical slow places. Once steady state has been reached the "timer" on the token in the first slow place will end at the same time as on the second, so a token will always enter the second place as soon as the previous token leaves -the time interval between tokens leaving the second place is exactly the waiting time for that place. The deterministic Petri net does not predict any effects due to clusters of slow codons. For a simple toy mRNA sequence it is possible to quantitatively predict how a cluster affects the protein production rate in the stochastic version of the Petri net model.
We consider an mRNA with very fast initiation and termination containing only very fast codons, except for a single slow codon somewhere in the bulk. We focus on the release time interval of tokens from the place corresponding to the slow codon. If the slow place has mean waiting time τ slow , and we assume that all of the other waiting times are so fast that this place is refilled with a new token practically immediately after the previous one leaves, then the mean release time interval for tokens from this place is equal to τ slow . Now consider the same system, but with a pair of slow codons labelled A and B -see Fig. 12. In the Petri net we denote the places corresponding to the slow codons place A and B. We denote the waiting times drawn from the (identical) distributions for each place τ A 1 , τ A 2 , . . . , and τ B 1 , τ B 2 , . . . .
We set t = 0 when the system is in the steady state and a token enters place B (since the other codons are very fast, a token also enters place A at this time). The first token leaves place B at time t = τ B 1 ; the next token leaves B at t = τ B 2 + max{τ A 1 , τ B 1 }, and the next at t = τ B 3 + max{τ A 2 , τ B 2 } + max{τ A 1 , τ B 1 }. That is, the ith token leaves place B at The ith interval ∆ i between two consecutive tokens leaving B is therefore Taking the average, and noting that the distribution for each slow place is the same, gives i.e., the mean release time interval is the mean of the maximum of two numbers drawn consecutively from the same exponential distribution. This can easily be shown 1 to be 3τ slow /2, where τ slow is the mean waiting time of the slow codons. This agrees with our earlier assertion that, in general, stochastic waiting times leads to slower movement of ribosomes along the mRNA.

Relationship to other translation models: TASEP
An alternative framework for describing mRNA translation, is the totally asymmetric simple exclusion process, or TASEP [17,18,15]. In that model ribosomes are represented by particles which hop in one direction along a 1D lattice of sites, which represent the codons of the mRNA. The current model shows many similarities to the TASEP, but the dynamics unfold in an essentially different way. The dynamics of the TASEP are most often simulated via Monte Carlo methods using a randomsequential update rule, or using continuous time Monte Carlo [19]. In the former case, time is discretised and in each time step lattice sites are chosen at random with uniform probability. If the chosen site contains a particle and the next site is vacant the particle is advanced with some probability p; particles are injected from the leftmost site and removed from the rightmost site with probabilities α and β, respectively. For a lattice with n sites this is repeated n times in each time step. In the continuous time version, exactly one particle is moved in each time step but the particle is chosen based on its probability of hopping, and the duration of the time step drawn from the corresponding Poisson distribution. Both methods are equivalent in that they are the realisation of the usual master equation in continuous time; this last point means that multiplying all probabilities by a common factor leads only to a rescaling of time. Under some conditions the random sequential TASEP can be solved exactly [18], and also a mean field treatment which ignore spatial correlations is often used.
The major difference between this and our timed Petri net picture is that in the TASEP only one particle can move at a certain instant in time. In the Petri net, as soon as any token move becomes possible, it is executed; i.e., events can happen simultaneously. Since in translation ribosomes actually take a finite time to move from one codon to the next, i.e., they are able to move at the same time, we propose that the Petri net more closely models the microscopic dynamics of this system. Results from the Petri net model can most easily be compared with the TASEP by identifying αδt = 1/τ in , βδt = 1/τ out and, in the case of a uniform mRNA, pδt = 1/τ, where δt is the time step. Although the  Table 3: Comparison of the results from the Petri net and TASEP models for a uniform mRNA. The TASEP results are from a mean field treatment which has been shown to very closely match randomsequential Monte Carlo simulations for large systems [18].
uniform TASEP produces a very similar phase diagram to that of Fig 4, the subtle differences in the dynamics do lead to some important differences in the macroscopic behaviour. The analytic results from the max-plus treatment for the Petri net model and from a mean field approximation for the TASEP are summarised in 3. As well as the density and protein production time differing in the two models, the boundaries of the regimes are also different. Specifically the boundaries for the tRNA capture limited regime (usually called the maximal current phase in the TASEP literature) are at τ in , τ out = 2τ instead ofτ in the deterministic Petri net. We note that the phase boundaries in the TASEP are the same as those in the stochastic Petri net model discussed in the previous section; we could therefore ascribe this to the stochasticity. However in the tRNA capture limited regime we have ρ = 1/2 and P = 4τ in the TASEP, compared to ρ = 1 and P =τ in the deterministic Petri net and ρ ≈ 3/4 and P ≈ 2.7τ in the stochastic Petri net. The dynamics of the TASEP are still different to those of the stochastic Petri net. A TASEP with an ordered-sequential update has been studied by some authors [20,21,22], and is particularly applicable to, for example, traffic flow models, movement of molecular motors, and -as we study here -mRNA translation. This update rule corresponds more closely to the Petri net model. Sites are taken in turn from right to left and updated according to the hopping probabilities. However, this has remained a less favourable update rule since analytic approaches run into several difficulties, particularly in the physical interpretation of hopping probabilities. Unlike the random-sequential case, for the ordered update time is not continuous, i.e., scaling the probabilities does not lead to a simple scaling of time (for further discussion see [21]). The advantage of the Petri net picture is therefore clear: we can include the fact that multiple events can happen simultaneously, whilst retaining (at least in the deterministic case) an analytically soluble mathematical framework, i.e., algebra on the max-plus semi-ring.
Another difference between the current model and the TASEP is in the biological interpretation of the waiting time of the ribosome on each codon. In the random-sequential TASEP the hopping probability from a site is chosen based on the abundance of tRNAs which corresponds to that codon. A recent study [23] shows that this is in fact the wrong interpretation of the hopping probability. Ciandrini et al. have developed an extension to the TASEP where ribosomes take two internal states. They identify two times: the waiting time for a ribosome to capture the correct tRNA τ capture , and the time it takes for the ribosome to physically move from one site to the next τ move . The former of these depends on the availability of tRNAs whilst the latter does not, and they argue that the capture of the tRNA can occur independently of whether or not there is a vacancy to the right of the ribosome. A commonly held misconception in applying the TASEP to translation is that the model represents the limit τ move → 0; Ciandrini's model shows that it is the opposite τ capture → 0 limit which recovers the original TASEP model; importantly, this limit does not describe the biologically relevant regime (τ move ≪ τ capture ). In the timed Petri net the place waiting times correspond to τ capture , hence we explicitly operate in the biologically relevant τ move → 0 limit.

Concluding Remarks
We have presented here a new model of ribosome dynamics during mRNA translation. Algebra on the max-plus semi ring lends its self to describing systems in which discrete events occur depending on the fulfilment of conditions, and has previously be used to study, for example, distributed software systems, automated manufacturing or industrial control systems. We have shown here that biological systems represent another area where such methods can be applied. The Petri net is a useful tool for visually representing discrete event systems, and the wealth of previous work on describing Petri nets using maxplus has allowed us to quickly develop a framework for predicting the dynamical behaviour of elongating ribosomes given only the codon sequence of an mRNA.
The analytic treatment we have presented is applicable to a deterministic version of the model, and we have discussed the modifications to the system when ribosome waiting times are chosen from a distribution. The present work therefore represents the first step in a new direction for protein synthesis modelling. Whilst numerical simulation of a stochastic Petri net is straightforward, there are currently few methods or algorithms available for the study of sequences of i.d.d. max-plus matrices. Therefore this work also represents a new source of motivation in this area.
There is also scope for introducing features which bring the model closer to the biology. For example including a finite time for the physical movement of the ribosomes [23], taking into account the fact that ribosomes cover more than one codon at a time [24,25], or allowing for a finite pool of ribosomes or other resources [26,27,28].
In summary we have presented a new model of translation which uses max-plus algebra to solve a Petri net description of the system. The microscopic dynamics are more realistic than in other translation models, and calculation of quantities such as density profiles and protein production rates is straightforward. Unlike other translation models the max-plus algebra also allows for exact analytic solutions for inhomogeneous mRNA sequences.