Geometric stability considerations of the ribosome flow model with pool
aa r X i v : . [ q - b i o . S C ] O c t Geometric stability considerations of the ribosome flow model with pool ∗ Wolfgang Halter , Jan Maximilian Montenbruck and Frank Allg¨ower Abstract — In order to better understand the process of genetranslation, the ribosome flow model (RFM) with pool wasintroduced recently. This model describes the movement ofseveral ribosomes along an mRNA template and simultaneouslycaptures the dynamics of the finite pool of ribosomes. Studyingthis system with respect to the number and stability of itsequilibria was so far based on monotone systems theory [1].We extend the results obtained therein by using a geometricapproach, showing that the equilibria of the system constitutea normally hyperbolic invariant submanifold. Subsequently, weanalyze the Jacobi linearization of the system evaluated at theequilibria in order to show that the equilibria are asymptoticallystable relative to certain affine subspaces. As this approach doesnot require any monotonicity features of the system, it mayalso be applied for more complex systems of the same kindsuch as bi-directional ribosome flows or time-varying templatenumbers.
I. I
NTRODUCTION
The ribosome flow model (RFM) as presented in [2] de-scribes the movement of ribosomes along an mRNA templateand can be used to study the dynamics of mRNA translation,an important step in the process of protein synthesis. Inorder to better understand and eventually design geneticregulatory networks (GRNs), the RFM with pool can be usedto provide a simulation framework which not only simulatesthe production of certain proteins but also takes into accountthe allocation of the ribosomes. Thus, a better evaluation ofthe performance of artificial GRNs, an important topic in thefield of synthetic biology ([3]), will be possible.Several versions of the RFM have been proposed recently, allof them built upon the works [4] and [5] where a probabilisticmodel of a growth center moving along a nucleic acidtemplate is considered. This model, also known as the totallyasymmetric exclusion process (TASEP) was then simplifiedand adapted in [2] to obtain the RFM, a deterministicdescription of the movement of several ribosomes alonga strand of mRNA. For this classical RFM, where theamount of ribosomes is assumed to be abundant and only asingle mRNA is studied, several results on model propertiessuch as uniqueness of the steady state and convergence tothis equilibrium point were presented in [1] wherein theauthors make use of the theory of monotone systems andthe contraction principle. Recently, [6] studied a network ofseveral mRNA templates in interaction with a finite poolof ribosomes, following up the argumentation of [1], theauthors show that the solution of the RFM network with *Supported by the research cluster BW Authors are with Institute for Systems Theory and AutomaticControl, 70569 Stuttgart, Germany.
E-Mail of correspondingauthor: [email protected] . pool monotonically converges to a certain equilibrium pointwhich is determined solely by the initial condition.In contrast to the monotone systems theory approach, wepropose a geometric approach to study the stability propertiesof the model in [6] in order to provide a concept to studynon-monotonic variants of the RFM in future.In particular, after introducing the detailed system descriptionof the RFM with pool, we show that the equilibria of theRFM with pool constitute a normally hyperbolic invariantsubmanifold and therefore the restriction of the linearizationof this system to the normal spaces of this submanifold canbe used to study the stability properties of the submanifold.Finally, we conclude that these equilibria are asymptoticallystable relative to certain affine subspaces.II. T HE RIBOSOME FLOW MODEL WITH POOL
As described in [2], the mechanism of translation canbe approximated as an initiation event followed by severalelongation steps. To be more precise, after binding to themRNA, the ribosomes perform a unidirectional movementalong the mRNA until they reach its end and subsequentlyunbind. In general, the speed of this motion is not constant,but to focus on the system theoretic analysis and for thesake of simplicity however, a constant elongation speed isassumed in the remainder.The RFM with pool can be modeled as the differentialequation ˙ R = − λR (1 − x ) + λ c x n (1) ˙ x ˙ x ... ˙ x n = λR (1 − x ) − λ c x (1 − x ) λ c x (1 − x ) − λ c x (1 − x ) ... λ c x n − (1 − x n ) − λ c x n , (2)with initial conditions R ( t = 0) = R tot ∈ [0 , ∞ ) (3) x i ( t = 0) = 0 , i = 1 , . . . , n. (4)The model states and parameters are explained in Table(I). The number of discretization points is determined suchthat each state can take up exactly one ribosome, therefore n = mlrl spec with ml the length of the mRNA template and rl spec the specific length a single ribosome occupies on themRNA template. Typically, the elongation rate λ c is cell typedependent and therefore might be determined from biologicaldata. The remaining free variables are ml , λ as well as R tot ,the total amount of available ribosomes. ABLE IS
TATES AND PARAMETERS OF THE
RFM
WITH POOL . R ∈ [0 , ∞ ) : molecular amount of free ribosomes x i ∈ [0 , : avg. ribosome density at mRNA location i = 1 , . . . , nλ ∈ R + : initiation rate λ c ∈ R + : elongation rate n ∈ N : number of discretization points on mRNA template For simpler notation, we collect all states such that x = (cid:2) R x x . . . x n (cid:3) ⊤ (5) ˙ x = f ( x ) (6)with x = R and x ∈ R n +1 . Further, as x , . . . , x n aredensities, it only makes sense to consider solutions t x ( t ) for which ∀ t ∈ [0 , ∞ ) , x ( t ) ∈ Ω := [0 , ∞ ) × [0 , × . . . × [0 , and as shown in [6], all solutions starting in Ω will not onlystay in this set but also are separated from the boundary ∂ Ω in finite time. III. M ODEL PROPERTIES
Similarly to the work of [6] and [1], we are interestedin the number of equilibra of System (6) as well as theirstability properties. In order to analyze similar models withslightly different conditions such as bidirectional flow onthe template ([7]) or even combinations of such systems,a geometric point of view for the stability analysis will turnout to be beneficial.
A. Equilibrium points
In order to find all equilibra of the RFM with pool, webring System (6) into the form ˙ x = A ( x ) x (7)with A ( x ) = − λ (1 − x ) 0 · · · λ c λ (1 − x ) − λ c (1 − x ) 0 · · · λ c (1 − x ) − λ c (1 − x ) 0 ...... ... ... · · · λ c (1 − x n ) − λ c (8)and rewrite the nullspace of A as ker A ( x ) = span( v ( x )) (9)with v ( x ) = h λ c λ (1 − x ) 11 − x · · · − x n i ⊤ . (10)This curve now represents the continuum of all equilibria of(6) through f ( x ) = 0 ⇔ x ∈ span( v ( x )) . However, notation(9) is quite unhandy and for that reason we will first give analternative representation of (9).Instead of using the span of a state dependent vector wecan define a parameterization γ : s γ ( s ) of the nullspace with s a scalar independent variable. Specifically, we define γ such that ∀ x ∈ int Ω : x ∈ ker A ( x ) ∃ s > x = γ ( s ) . (11)Calculating γ is straight forward and can be achieved bymultiplying v with the independent variable s and thenrecursively solving γ ( s ) = sv ( x ) and substituting all x i ,which then results in γ : s (cid:2) γ ( s ) . . . γ n ( s ) (cid:3) ⊤ (12)with the components γ i ( s ) given recursively as a series ofcontinued fractions with γ i ( s ) = λ c sλ (1 − γ ( s )) i = 0 s − γ i +1 ( s ) i = 1 . . . ( n − s i = n. (13)In the remainder, we restrict our attention to solutionsinitialized in int Ω , the interior of Ω . This is justified assolutions initialized on ∂ Ω attain values in int Ω in finitetime. Under this assumption we notice that ∃ ¯ s > s ∈ (0 , ¯ s ) ⇒ γ ( s ) ∈ int Ω . (14)We henceforth restrict the domain of γ to (0 , ¯ s ) . This isa rather technical assumption to ensure that we study theRFM with pool in a domain which makes sense biologically.It is further possible to show that for s ∗ > ¯ s , γ ( s ∗ ) / ∈ Ω and for this case, the model ceases to have any biologicalmeaning. The value of ¯ s is only dependent on the numberof discretization points n and further is a solution of thepolynomial equation j n − k X j =0 ( − ¯ s ) j +1 (cid:18) n − jj + 1 (cid:19) = 0 . (15)The derivation of this result is given in the appendix and wenote on the side that this ¯ s can be used to calculate an upperbound on the protein production rate κ = λ c ¯ s of the studiedmRNA template.With this representation of the equilibria of (6) at hand, wecan proceed with studying their stability properties. B. Stability of equilibria
In this section, we consider the Jacobian of f in order tostudy the stability properties of the equilibria: J f ( x ) = ∂f ∂x ( x ) · · · ∂f ∂x n ( x ) ... . . . ... ∂f n ∂x ( x ) · · · ∂f n ∂x n ( x ) (16) = A ( x ) + λx · · · − λx λ c x .. . ... . . . − λ c x .. . ... .. . λ c x n − · · · − λ c x n − (17)nd evaluate this Jacobian at γ ( s ) , i.e. J f ( γ ) = − λ (1 − γ ) λγ · · · λ c λ (1 − γ ) − λ c (1 − γ ) − λγ λ c γ λ c (1 − γ ) − λ c (1 − γ + γ ) ...... λ c (1 − γ ) ... λ c γ n − · · · λ c (1 − γ n ) − λ c (1 + γ n − ) (18)where we omitted the argument s for the sake of readability. Theorem 1:
For all s > the Jacobian matrix of f ,evaluated at γ ( s ) , has an eigenvalue equal to . Further, allremaining eigenvalues have real parts strictly smaller thanzero. Proof:
In (18) one can see that all diagonal elementsof J f ( γ ) are strictly negative as γ i < for i = 1 , . . . , n andfurther, with J k,if ( γ ) being the element of J f ( γ ) in the k -throw and i -th column, X k = i | J k,if ( γ ) | = − J i,if ( γ ) . (19)This means that, using Gerschgorin circles ([8]) with theircenter coordinates given by the value of the diagonal ele-ments J i,if ( γ ) and their radius given by P k = i | J k,if ( γ ) | , alleigenvalues have a real part smaller or equal to zero.It remains to show that J f ( γ ) has precisely one eigenvalueequal to zero which then also implies that the remainingeigenvalues have a real part strictly smaller than zero. There-fore we use the theorem on the reduced row echelon formand bring J f ( γ ) into upper triangular form, such that J f ( γ ) = LU, with (20) L = · · · − . . . . . . ...... − · · · − (21) U = − λ (1 − γ ) λγ · · · λ c − λ c (1 − γ ) λ c γ λ c ... − λ c (1 − γ ) ... ...... ... λ c γ n − λ c − λ c (1 − γ n ) λ c (1 + γ n − )0 · · · . (22)Now that the rank of U is n − , this concludes the proof.One might be tempted to use Lyapunov’s indirect method([9]), that asymptotic stability of a hyperbolic equilibriumis determined by the Jacobi linearization, or the theoremof Hartman-Grobman ([10]), that a vector field and itslinearization are conjugate in a neighborhood of a hyperbolicequilibrium, in order to draw conclusions about the stabilityproperties of the equilibria. However, the existence of thezero eigenvalue means that the equilibria are non-hyperbolicand these methods are not applicable. Yet, [11] offers anextension to normally hyperbolic invariant manifolds in thefollowing sense: a vector field and the restriction of its linearization to the normal spaces of a given normally hy-perbolic invariant manifold are conjugate in a neighborhoodof the normally hyperbolic invariant manifold. Therefore,following the notation of [12], we will show in the nextsection that f is normally hyperbolic at γ ((0 , ¯ s )) in order tocontinue with the stability analysis. C. Normal hyperbolicity
In this section, we proceed as follows: we already notedthat γ ((0 , ¯ s )) is a manifold of equilibria (and thus invari-ant) and that J f ( γ ) , the Jacobian of f evaluated on γ ,has precisely one zero eigenvalue. Next, we show that theeigenvector associated with the zero eigenvalue lies in thetangent space T γ ( s ) γ ((0 , ¯ s )) = span( { ˙ γ ( s ) } ) of γ ((0 , ¯ s )) at γ ( s ) , no matter at which s ∈ (0 , ¯ s ) we evaluate γ . Further, asshown in Theorem (1), all other eigenvalues of the Jacobianare negative if evaluated on γ ((0 , ¯ s )) and we will show thattheir eigenvectors span N γ ( s ) γ ((0 , ¯ s )) , the normal space of γ ((0 , ¯ s )) at γ ( s ) in R n +1 . In conclusion, γ ((0 , ¯ s )) is notonly normally hyperbolic, but further asymptotically stablein a sense we detail further below. Lemma 1:
The eigenvector associated with the zero eigen-value of J f ( γ ( s )) is linearly dependent on ˙ γ ( s ) = dds γ ( s ) . Proof:
It is sufficient to show that J f ( γ ) ˙ γ = 0 . (23)We thus consecutively show that J f ( γ ) ˙ γ = 0 (24) J if ( γ ) ˙ γ = 0 i = 1 , . . . , n − (25) J nf ( γ ) ˙ γ = 0 (26)with J if the i + 1 -th row of the Jacobian J f . We start withshowing (26): J nf ( γ ) ˙ γ = λ c (1 − γ n ) ˙ γ n − − λ c (1 + γ n − ) ˙ γ n (27) = λ c (1 − s ) 1( s − − λ c (1 + s − s ) (28) = λ c (cid:18) − s − s − (cid:19) = 0 . (29)Showing (25) can now be achieved in a general form sincefor all i = 1 , . . . , n − the structure of J if is identical, namely J if ( γ ) ˙ γ = λ c (1 − γ i ) ˙ γ i − − λ c (1 − γ i +1 + γ i − ) ˙ γ i + λ c γ i ˙ γ i +1 . (30)We rearrange the last equation to get J if ( γ ) ˙ γ = λ c ( ˙ γ i − − γ i ˙ γ i − − γ i − ˙ γ i ) − λ c ( ˙ γ i − γ i +1 ˙ γ i − γ i ˙ γ i +1 ) (31) = λ c (cid:16) ˙ γ i − − ˙( γ i − γ i ) (cid:17) − λ c (cid:16) ˙ γ i − ˙( γ i γ i +1 ) (cid:17) . (32)ow, we can utilize the generating equation (13) again torealize that ˙ γ i − ˙( γ i γ i +1 ) = ( λ c λ i = 01 i = 1 , . . . , n − (33)and using (33) for i = 1 , . . . , n − to arrive at the equality J if ( γ ) ˙ γ = 0 i = 1 , . . . , n − . (34)Finally, we merely need to verify whether this is also truefor (24): J f ( γ ) ˙ γ = − λ (1 − γ ) ˙ γ + λγ ˙ γ + λ c ˙ γ n (35) = λ ( − ˙ γ + ˙ γ γ + ˙ γ γ ) + λ c ˙ γ n . (36)Now with ˙ γ n = 1 , and using equation (33) for i = 0 , J f ( γ ) ˙ γ = λ (cid:18) − λ c λ (cid:19) + λ c = 0 . (37)This concludes the proof.One says that f is normally hyperbolic at γ ((0 , ¯ s )) ifthe derivative of f evaluated at γ ( s ) leaves the continuoussplitting R n +1 = N u ⊕ T γ ( s ) γ ((0 , ¯ s )) ⊕ N s (38)invariant and if the normal behavior dominates the tangentone. N u and N s in that respect are the subspaces spannedby the normal eigenvectors of the Jacobian associated withpositive and negative eigenvalues respectively. Due to The-orem (1) we know that N u = ∅ . Lemma (1) shows thatthe dynamics of f on T γ ( s ) γ ((0 , ¯ s )) is determined by thezero eigenvalue and it remains to show that the (generalized)eigenvectors of the remaining eigenvalues span the normalspace of γ ((0 , ¯ s )) at any γ ( s ) . In order to do so, we definethe affine subspaces S p := ( x (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) n X i =0 x i = p ) (39)and note the following proposition: Proposition 1:
Any solution of System (6) initialized ona certain S p will stay on S p for all times and p = R tot . Proof:
Assume the initial state of System (6) lies in S p , i.e. p = R tot , then choose V ( x ) = n X i =0 x i (40)as a Lyapunov function such that S p is just the level set V − ( { p } ) of V . Now consider the Lie derivative L f V of V along f , L f V = ∂V∂ x · f (41) = (cid:2) · · · (cid:3) f (42) = 0 . (43)This shows that any level set of V is invariant under (6).With Proposition (1) at hand, it remains to show that theintersection of γ with all S p is always transversal and nevertangential to conclude normal hyperbolicity of f at γ ((0 , ¯ s )) . Lemma 2:
For all p > , the curve γ intersects S p uniquely and transversely.In other words, this means that the continuum of equilibriaof system (6) represented by the curve γ never intersects the n -dimensional affine subspace S p tangentially. Proof:
We start with showing the transversality of theintersection of γ and S p as the uniqueness of this intersectionthen follows as we will show at the end of the proof. In orderto do so, we note that the affine subspace S p is always thesame subspace M translated in direction of x S p = { e p } + M (44)with e the canonical unit vector and the subspace M definedas the image of Λ = (cid:2) µ µ . . . µ N (cid:3) = − − · · · −
11 0 · · ·
00 1 ...... · · · . (45)This means that the only vectors perpendicular to M aremultiples of the all ones vector and it is thus sufficient toshow that the velocity vector of γ never is perpendicular to , in other words h ˙ γ, i 6 = 0 . (46)Showing this will be achieved by noting that ˙ γ i is positivefor all i = 0 , . . . , n . This can be achieved by the followinginductive argument starting with the highest appearing index n and then reducing it step wise.First, we consider γ n ( s ) = s and see that ˙ γ n ( s ) = 1 > . (47)Next, we assume that the statement that ˙ γ i > (48)is true. Now further reducing the index, we need to showthat (48) still holds for ˙ γ i − . Therefore we differentiate thegenerating equation (13) for index i − to arrive at ˙ γ i − ( s ) = (1 − γ i ) + s ˙ γ i (1 − γ i ) . (49)As shown in (14) we already know that γ i ∈ (0 , ∀ i > .Together with (48), this means that ˙ γ i − ( s ) > (50)which concludes the proof for the transversality of theintersection. Now that we further know that all derivativesof γ with respect to s are larger than zero for s ∈ (0 , ¯ s ) ,the uniqueness of the intersection follows directly from thecombination of these arguments.Fig. (1) illustrates two affine subspaces S (red) and S (green) and the continuum of equilibria γ ((0 , )) (blue) inthe first three coordinates. One can see that the subspaces arejust shifted by the difference in total amount of ribosomes R tot in R direction and that the curve intersects with each . . x x R S S γ ( s ) Fig. 1. Two affine subspaces S and S and the nullspace of A ( x ) depictedfor the first three dimensions. subspace uniquely and not tangentially. We are now able toformulate our main result. Theorem 2:
The invariant set γ ((0 , ¯ s )) of (6) is asymptot-ically stable. Proof:
With Lemmata (1) and (2) at hand we canconclude that f is normally hyperbolic at γ ((0 , ¯ s )) andtherefore f and the restriction of its linearization to the nor-mal spaces of γ ((0 , ¯ s )) are conjugate in a neighborhood of γ ((0 , ¯ s )) . As shown in Theorem (1), the eigenvalues of thislinearization evaluated at γ ( s ) are all strictly negative exceptfor one, which is exactly zero. Now that the eigenvectorassociated with this zero eigenvalue lies in the tangent space T γ ( s ) γ ((0 , ¯ s )) and we are only interested in the eigenvectorslying in the normal spaces of γ ((0 , ¯ s )) , which thus are allassociated with strictly negative eigenvalues, this concludesthe proof.We further showed in Proposition (1) that the affine sub-spaces S p are also invariant under the system dynamics.This further means that any solution initialized on a certain S p with p > will converge to the unique equilibriumgiven by the intersection of γ ( s ) with S p . In other words,this intersection point is asymptotically stable relative to S p ,following the terminology of [13].IV. C ONCLUSION
We considered the RFM with pool, a model describing themovement of ribosomes along a single mRNA template aswell as the dynamics of a pool of available ribosomes. Wefound that the equilibria of the system can be characterizedby a curve γ and that there exist affine subspaces S p whichare invariant under the dynamics of the system. In orderto characterize the stability of the equilibria we studiedthe Jacobi linearization of the system evaluated on γ and found that all eigenvalues are smaller than zero except forprecisely one which is equal to zero, therefore the equilibriaare non-hyperbolic. In order to draw any conclusions fromthe linearization we showed that the system under studyis normally hyperbolic at γ . This was achieved in twosteps, first, showing that the eigenvector associated with thezero eigenvalue of the Jacobian evaluated on γ is linearlydependent on the velocity vector of γ and second, showingthat γ intersects all S p transversely. This insight then enabledus to apply the results of [11] in order to conclude that theequilibria of the system are asymptotically stable relative tothe affine subspaces S p .In previous works on the RFM ([6]) monotone systemstheory was used to show that every equilibrium point issemistable in a sense that any solution initialized on a certain S p monotonically converges to a unique equilibrium pointwhich is dependent on the initial condition. Our results nowoffer a more detailed characterization of the stability of theRFM with pool and further introduce a geometric approachfor studying similar systems with higher complexity as thisapproach does not require any monotonicity features of thesystem. Such more complex systems may for instance bemodels where the copy number of templates is varying overtime or the flow of ribosomes is allowed to be bi-directional.R EFERENCES[1] M. Margaliot and T. Tuller, “Stability Analysis of the RibosomeFlow Model,”
IEEE/ACM Transactions on Computational Biologyand Bioinformatics , vol. 9, no. 5, pp. 1545–1552, sep 2012.[2] S. Reuveni, I. Meilijson, M. Kupiec, E. Ruppin, and T. Tuller,“Genome-scale analysis of translation elongation with a ribosome flowmodel,”
PLoS Computational Biology , vol. 7, no. 9, 2011.[3] F. Ceroni, R. Algar, G.-B. Stan, and T. Ellis, “Quantifying cellularcapacity identifies gene expression designs with reduced burden,”
Nature Methods , vol. 12, no. 5, pp. 415–418, 2015.[4] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, “Kinetics ofbiopolymerization on nucleic acid templates,”
Biopolymers , vol. 6,no. 1, pp. 1–25, jan 1968.[5] R. Heinrich and T. A. Rapoport, “Mathematical modelling oftranslation of mRNA in eucaryotes; steady states, time-dependentprocesses and application to reticulocytest,”
Journal of TheoreticalBiology , vol. 86, no. 2, pp. 279–313, sep 1980.[6] A. Raveh, M. Margaliot, E. D. Sontag, and T. Tuller, “A Modelfor Competition for Ribosomes in the Cell,” arXiv:1508.02408[q-bio.GN] , aug 2015.[7] S. Edri, E. Gazit, E. Cohen, and T. Tuller, “The RNA Polymerase FlowModel of Gene Transcription,”
IEEE Transactions on BiomedicalCircuits and Systems , vol. 8, no. 1, pp. 54–64, feb 2014.[8] S. Gerschgorin, “ ¨Uber die Abgrenzung der Eigenwerte einer Matrix,”
Izv. Akad. Nauk. USSR Otd. Fiz.-Mat. Nauk , no. 6, pp. 749–754, 1931.[9] H. K. Khalil,
Nonlinear Systems . Prentice Hall, 2002.[10] P. Hartman, “A lemma in the theory of structural stability ofdifferential equations,”
Proceedings of the American MathematicalSociety , vol. 11, no. 4, pp. 610–610, apr 1960.[11] C. Pugh and M. Shub, “Linearization of normally hyperbolicdiffeomorphisms and flows,”
Inventiones Mathematicae , vol. 10,no. 3, pp. 187–198, sep 1970.[12] M. W. Hirsch, C. C. Pugh, and M. Shub,
Invariant Manifolds , ser.Lecture Notes in Mathematics. Springer Berlin Heidelberg, 1977,vol. 583.[13] N. P. Bhatia and G. P. Szeg¨o,
Stability Theory of Dynamical Systems ,1970.[14] R. M. Dudley, “Some Inequalities for Continued Fractions,”
Mathe-matics of Computation , vol. 49, no. 180, pp. 585–593, 1987.
PPENDIX