Behavior of Totally Positive Differential Systems Near a Periodic Solution
Chengshuai Wu, Lars Gruene, Thomas Kriecherbauer, Michael Margaliot
aa r X i v : . [ m a t h . D S ] J a n Behavior of Totally Positive Differential Systems Near a Periodic Solution
Chengshuai Wu, Lars Gr¨une, Thomas Kriecherbauer, and Michael Margaliot
Abstract — A time-varying nonlinear dynamical system iscalled a totally positive differential system (TPDS) if its Ja-cobian admits a special sign pattern: it is tri-diagonal withpositive entries on the super- and sub-diagonals. If the vectorfield of a TPDS is T -periodic then every bounded trajectoryconverges to a T -periodic solution. In particular, when thevector field is time-invariant every bounded trajectory of aTPDS converges to an equlbrium. Here, we use the spectraltheory of oscillatory matrices to analyze the behavior near aperiodic solution of a TPDS. This yields information on theperturbation directions that lead to the fastest and slowestconvergence to or divergence from the periodic solution. Wedemonstrate the theoretical results using a model from systemsbiology called the ribosome flow model. I. I
NTRODUCTION
Consider the time-varying nonlinear system: ˙ x ( t ) = f ( t, x ( t )) , (1)where f : R + × Ω → R n is continuously differentiable and Ω is a convex subset of R n . For t ≥ t ≥ and a ∈ Ω ,let x ( t, t , x ) denote the solution of (1) at time t withthe initial condition x ( t ) = x . From here on we alwaystake t = 0 and write x ( t, x ) for x ( t, , x ) . We assumethroughout that for any x ∈ Ω the solution x ( t, x ) ∈ Ω forall t ≥ .Let J ( t, x ) := ∂∂x f ( t, x ) denote the Jacobian of thevector field f . Let M + ⊂ R n × n denote the set of n × n tri-diagonal matrices with positive entries of the super-and sub-diagonals, i.e. the subset of Jacobi matrices. Sys-tem (1) is called a totally positive differential system (TPDS)if J ( t, x ) ∈ M + for all ( t, x ) . In particular, J ( t, x ) isa Metzler matrix (i.e. a matrix whose off-diagonal entriesare all non-negative), and also an irreducible matrix. Thisimplies that (1) is strongly cooperative [27], i.e. if a, b ∈ Ω ,with a = b , and a − b ∈ R n + := { x ∈ R n : x i ≥ , i = 1 , . . . , n } then x ( t, a ) − x ( t, b ) ∈ R n ++ := { x ∈ R n : x i > , i = 1 , . . . , n } for all t > . Cooperative systems play an importantrole in models from biology and neuroscience, where it isknown that the state-variables represent quantities that cannever attain negative values, e.g. concentration of molecules,average number of spikes in a neuron, etc. Cooperative CW and MM are with the School of Electrical Engineering, Tel AvivUniversity, Israel. LG and TK are with the Mathematical Institute, Universityof Bayreuth, Germany. This research is partially supported by a researchgrant from the ISF. Correspondence: [email protected] systems enjoy a well-ordered behavior. By Hirsch’s quasi-convergence theorem [27], almost every bounded trajectoryof a time-invariant (i.e. f ( t, x ) = f ( x ) ) cooperative systemconverges to the set of equilibria.TPDSs require a stronger condition on the Jacobian andconsequently enjoy stronger properties. Smillie [25] provedthat if the vector field is time-invariant then every boundedtrajectory converges to an equilibrium. Smith [26] general-ized this result and showed that if the vector field f of aTPDS is T -periodic, that is, f ( t, z ) = f ( t + T, z ) for all t ≥ , z ∈ R n then every bounded trajectory converges to a T -periodicsolution of the system. These results found many applicationsin fields such as systems biology and neuroscience [13], [4],[6], as well as several generalizations (see, e.g. [9], [28]).The proofs of Smillie and Smith are based on using thenumber of sign changes in the vector of derivatives ˙ x ( t ) asan integer-valued Lyapunov function. Recently, it was shownthat these results are closely related to the sign-variationdiminishing property of totally positive matrices [12]. Recallthat a matrix A ∈ R n × n is called totally non-negative [TN]if all its minors are non-negative, and totally positive [TP]if they are all positive. Such matrices have a rich andbeautiful structure [18], [8]. For example, if A ∈ R n × n is TPand v ∈ R n \ { } , then the number of sign variations in Av is smaller or equal to the number of sign variations in v .Schwarz [23] introduced and analyzed linear TPDSs. Heshowed in particular that the following two conditions areequivalent: 1) A ∈ M + ; 2) the transition matrix exp( At ) of the linear time-invariant (LTI) system ˙ x = Ax is TP forany t > . Unfortunately, his results were almost forgotten.In many dynamical systems it is important to under-stand how a small perturbation affects the behavior near aperiodic solution (and, in particular, near an equilibrium).In general, this is a difficult problem, as there is littleexplicit information on the periodic solution. Here, we usethe spectral theory of TP matrices to analyze trajectoriesof the T -periodic TPDS (1) in the vicinity of a periodicsolution (and in the special case of a time-invariant vectorfield, near an equilibrium). Our main result characterizesthe sign pattern of the “most stable” and “most unstable”perturbation directions near a periodic solution. We alsoprovide an intuitive explanation for the structure of these signpatterns. We demonstrate an application of these theoreticalresults to an important model from systems biology calledthe ribosome flow model (RFM) [32].We note that when n = 2 and (1) is a strictly cooperativesystem then it is also a TPDS (as a × matrix is alwaysri-diagonal), so our results hold in this special case as well.The remainder of this note is organized as follows. Thenext section briefly reviews several known definitions andresults that are used later on. Section III presents our mainresults, and describes an application of these results tothe RFM. The final section concludes. We use standardnotation. Vectors [matrices] are denoted by small [capital]letters. For A ∈ R n × m , A T is the transpose of A .II. P RELIMINARIES
We begin by briefly reviewing known results that will beused later on.
A. Sign variation diminishing property
An important property of TP matrices is that multiplyinga vector by a TP matrix can never increase the number ofsign variations. To explain this, we recall two definitionsfor the number of sign changes in a vector. For z ∈ R n \ { } , let s − ( z ) denote the number of sign changes inthe vector z , after deleting all the zero entries (with s − (0) defined as zero). For example, s − ( (cid:2) − (cid:3) T ) =2 . For z ∈ R n , let s + ( z ) denote the maximal possiblenumber of sign changes in the vector z , after replacingevery zero entry by either plus one or minus one. Forexample, s + ( (cid:2) − (cid:3) T ) = 4 . Note that thesedefinitions imply that ≤ s − ( z ) ≤ s + ( z ) ≤ n − , for any z ∈ R n . (2)The next result describes the sign variation diminishingproperty of TP matrices. Proposition 1: [8] Let A ∈ R n × n be TP. Then s + ( Ax ) ≤ s − ( x ) for all x ∈ R n \ { } . B. Oscillatory matrices
A matrix A ∈ R n × n is called oscillatory if it is TN andthere exists an integer k > such that A k is TP [10]. Thesmallest such k is called the exponent of the oscillatorymatrix A [30]. A TN matrix is oscillatory if and only if (iff)it is non-singular and irreducible [10]. Oscillatory matriceshave special spectral properties: their eigenvalues are allreal, simple, and positive, and the corresponding eigenvectorshave a special sign pattern [18]. The next result summarizesthis spectral structure. Proposition 2: [18] If A ∈ R n × n is oscillatory then allits eigenvalues are real, positive, and simple. Denote theseeigenvalues as α > α > · · · > α n > , and let v k ∈ R n \ { } denote an eigenvector corresponding to α k . Forany ≤ i ≤ j ≤ n , let V ij := span { v i , v i +1 , . . . , v j } . Then i − ≤ s − ( z ) ≤ s + ( z ) ≤ j − for any z ∈ V ij \ { } . Remark 3:
Note that this implies in particular that s − ( v i ) = s + ( v i ) = i − , for any i ∈ { , . . . , n } . (3)In general, the eigenvectors may include zero entries. How-ever, the zeros must be located such that (3) holds. In partic-ular, v cannot include any zero entries because s + ( v ) = 0 , and v n cannot include any zero entries because s − ( v n ) = n − . Example 1:
Let A = . . All the × mi-nors of A (i.e., its entries) are non-negative. The × minors are , , , . , , , . , . , , and the single × minor is det( A ) = 25 . . Thus, A is TN. Calculating allthe minors of A shows that they are all positive, so A is TP. Thus, A is oscillatory. Its eigenvalues are α =5 . , α = 3 . , α = 1 . , (all numericalvalues in this paper are to 5-digit accuracy) with correspond-ing eigenvectors: v = (cid:2) . . . (cid:3) T , v = (cid:2) . . − . (cid:3) T , and v = (cid:2) . − . . (cid:3) T . Note that s − ( v i ) = s + ( v i ) = i − , i ∈ { , , } .For our purposes, it is important to know when a tri-diagonal matrix is an oscillatory matrix. C. Tri-diagonal oscillatory matrices
Consider the n × n tri-diagonal matrix T ( a i , b i , c i ) := a b . . . c a b . . . c a . . . ... ... ... . . . ... ... ... . . . c n − a n − b n − . . . c n − a n . (4)Assume that b i , c i ≥ , and that the entries satisfy the dominance condition : a i ≥ b i + c i − , i = 1 , . . . , n, (5)with b n := 0 , and c := 0 . Then all the minors of T arenon-negative, i.e. T is TN [8, Chapter 0]. If, furthermore, T is non-singular and irreducible then it is oscillatory [10]. D. Linear TPDSs
Schwarz [23] considered the system ˙ x ( t ) = A ( t ) x ( t ) , (6)with t → A ( t ) continuous. The system is called a linearTPDS if the corresponding transition matrix Φ( t, t ) (i.e, thematrix such that x ( t ) = Φ( t, t ) x ( t ) for all t ≥ t ≥ and x ( t ) ∈ R n ) is TP for any t > t ≥ . Schwarz showedthat this holds iff: 1) a ij ( t ) ≡ for all | i − j | > , 2) a ij ( t ) ≥ for all | i − j | = 1 , and 3) any of the functions in 2) doesnot vanish identically on any interval of positive length. Inparticular, if A ( t ) ≡ A is constant then the system is a TPDSiff A ∈ M + . One important property of a linear TPDS, thatfollows from Prop. 1, is that for any non-trivial solution x ( t ) , s + ( x ( t )) ≤ s − ( x ( t )) , for all t > t ≥ . (7)Thus, the number of sign variations along a solution ofthe system is non-increasing. Furthermore, it can be shownthat s − ( x ( t )) = s + ( x ( t )) except perhaps at up to n − isolated points [12].ne implication of this is that if γ ( t ) is a T -periodic solu-tion of a linear TPDS (that is not the trivial solution γ ( t ) ≡ )then s − ( γ ( t )) = s + ( γ ( t )) = s − ( γ (0)) for all t ≥ . The next section describes our main results. It is usefulto begin with case of a linear TPDS, and then consider thenonlinear case. III. M
AIN R ESULTS A. T -periodic linear TPDS Consider the linear-time varying (LTV) system (6)with A ( t ) a continuous and T -periodic matrix. We alsoassume that A ( t ) ∈ M + for all t ∈ [0 , T ) , so (6) is alinear TPDS. Then any trajectory with an initial conditionin Ω converges to a T -periodic solution of (6) (that is notnecessarily unique).The transition matrix satisfies Φ(0) = I and by Floquettheory [3], Φ( t + T ) = Φ( t )Φ( T ) , for all t ≥ . Let B :=Φ( T ) . The eigenvalues of B are called the characteristicmultipliers of (6). Since A ( t ) ∈ M + , Φ( t ) is TP for any t > and in particular B is TP. Let λ > λ > · · · > λ n > , (8)denote the eigenvalues of B , and let v i ∈ R n denote theeigenvector corresponding to λ i .Let γ ( t ) denote a T -periodic solution of (6) with γ (0) = 0 .Then γ ( T ) = Bγ (0) and since γ ( T ) = γ (0) , this impliesthat there exists p ∈ { , . . . , n } such that λ > λ > · · · > λ p = 1 > · · · > λ n > , (9)Fix a “perturbation direction” w ∈ R n and consider thedifference z ( t ) := x ( t, γ (0) + w ) − x ( t, γ (0)) . (10)Then for any integer k > , z ( kT ) = Φ( kT ) z (0)= B k z (0)= B k w. Let c i ∈ R be such that w = P ni =1 c i v i . Then z ( kT ) = n X i =1 c i λ ki v i . (11)Combining this with (9) implies the following result. Proposition 4:
Consider the linear TPDS (6), and assumethat it admits a T -periodic solution γ ( t ) with γ (0) = 0 .There are two possible cases.Case 1. If λ > then the periodic solution is not orbitallystable, v is the “worst” perturbation direction in the sensethat this perturbation takes the state away from γ as quicklyas possible. Also, λ n ≤ , and if λ n < then v n is the “best”perturbation direction in the sense that this perturbation takesthe state back to γ as quickly as possible. Case 2. If λ = 1 then (9) implies that the periodic solutionis stable and orbitally asymptotically stable, and that v [ v n ]is the “worst” [“best”] perturbation direction in the sense thatthis perturbation takes the state back to γ as slowly [quickly]as possible.Note that for a TPDS, even if we do not have an explicitexpression for γ ( t ) , we always know that the eigenvalues λ i of B = Φ( T ) are real, positive, and distinct, and also thespecial sign pattern of each of the eigenvectors v i .The next simple example demonstrates Prop. 4. Example 2:
Consider (6) with n = 2 and A ( t ) = (cid:20) − t )2 + sin( t ) − (cid:21) . Note that A ( t ) ∈ M + for all t ≥ , and that A ( t ) is T -periodic for T = 2 π . The transition matrix is Φ( t ) = exp( − t ) (cid:20) cosh( c ( t )) sinh( c ( t ))sinh( c ( t )) cosh( c ( t )) (cid:21) , with c ( t ) := 2 t − cos( t ) + 1 . Note that every entry of Φ( t ) ispositive and that det(Φ( t )) > (i.e. Φ( t ) is TP) for all t > .Here, B = Φ(2 π ) = exp( − π ) (cid:20) cosh(4 π ) sinh(4 π )sinh(4 π ) cosh(4 π ) (cid:21) . The eigenvalues of B are λ = 1 , λ = exp( − π ) , with corresponding eigenvectors v = (cid:2) (cid:3) T , v = (cid:2) − (cid:3) T . Hence, a π -periodic solution is γ ( t ) := x ( t, v )= Φ( t ) v = exp( − t )(cosh( c ( t )) + sinh( c ( t ))) v = exp(1 − cos( t )) v . As expected, v [ v ] has zero [one] sign changes. To demon-strate (11) in this case, define z as in (10). Any w ∈ R canbe written as w = w + w v + w − w v , so z ( kT ) = B k w = w + w v + w − w − kπ ) v . B. T -periodic nonlinear TPDS Consider the nonlinear system (1) with a T -periodicvector field. Assume that its trajectories evolve on the state-space Ω ⊆ R n , and that J ( t, z ) ∈ M + for all t ∈ [0 , T ) , z ∈ Ω , (12)that is, (1) is a TPDS.We first consider the case where f is time-invariant (andhence T -periodic for any T ). The next result describes thespecial spectral structure of J ( z ) at any point z ∈ Ω . Proposition 5:
Fix z ∈ Ω . The eigenvalues of J ( z ) areeal and simple. Denote these eigenvalues by α ( z ) > α ( z ) > · · · > α n ( z ) , and let v k ( z ) ∈ R n denote the eigenvector correspondingto α k ( z ) . For any ≤ i ≤ j ≤ n , let V ij ( z ) :=span { v i ( z ) , v i +1 ( z ) , . . . , v j ( z ) } . Then i − ≤ s − ( y ) ≤ s + ( y ) ≤ j − , for any y ∈ V ij ( z ) \ { } . (13) Proof:
Fix z ∈ Ω . Then J ( z ) is a Jacobi matrix.Pick s > large enough so that sI + J ( z ) satisfies thedominance condition (5), and is nonsingular. Then sI + J ( z ) is TN, nonsingular and irreducible, so it is oscillatory.Applying Prop. 2 completes the proof.In particular, if e is an equilibrium then the set ofeigenvalues and eigenvectors of J ( e ) provides a rathercomplete picture of the dynamical behavior near e . Indeed,the Hartman–Grobman theorem [3] asserts that the phaseportrait near e (assuming that J ( e ) has no eigenvalues witha zero real part) is the same as the phase portrait of theLTI system ˙ y = J ( e ) y , up to a continuous change ofcoordinates x = p ( y ) , with p : R n → R n satisfying p (0) = e . Suppose for example that all the eigenvalues are negative.Then the eigenvector v ( e ) [ v n ( e ) ] corresponds to the “slow-est” [“fastest”] eigenvalue α ( e ) [ α n ( e ) ]. A perturbation ofthe steady-state in the form e → p ( εv ( e )) [ e → p ( εv n ( e )) ],with | ε | small, will be “compensated” at the slowest [fastest]possible rate. Example 3:
Consider the system: ˙ x = − x + tanh( x ) + 2 tanh( x ) , ˙ x = − x + (1 /
2) tanh( x ) + tanh( x ) . (14)Such systems appear in models of neural networkswith tanh( · ) as the activation function. It is straightforwardto verify that there are three equilibrium points in R : e = (cid:2) (cid:3) T , e = (cid:2) . . (cid:3) T , and e = − e .The Jacobian of (14) is J ( x ) = (cid:20) − − ( x ) 2 cosh − ( x )(1 /
2) cosh − ( x ) − − ( x ) (cid:21) , (15)so J ( x ) ∈ M + for all x ∈ R . The eigenvalues of J ( e ) are ( √ − / , ( −√ − / , with corresponding eigenvectors (cid:2) √ − (cid:3) T , (cid:2) −√ − (cid:3) T . The eigenvalues of J ( e ) (and of J ( e ) ) are − . , − . , with corresponding eigenvectors (cid:2) . . (cid:3) T , (cid:2) . − . (cid:3) T . As expected, the eiegnvalues are real and distinct, and theeigenvectors have the specified sign pattern. (cid:3)
We now turn to consider the case where the nonlinear system is time-varying, and T -periodic with a minimalperiod T > . For x ∈ Ω , let H ( t, x ) := ∂∂x x ( t, x ) , that is, H maps a change in the initial condition at time to the change in the solution at time t . Then H (0 , x ) = I, and ˙ H ( t, x ) := ddt H ( t, x )= ∂∂x f ( t, x ( t, x ))= J ( t, x ( t, x )) H ( t, x ) . (16)Since J ∈ M + , this is a linear time-varying TPDS. Let a := γ ′ (0) , where γ ( t ) is a T -periodic solution of (1). Thenfor x = a the linear TPDS (16) is also T -periodic,so Φ( t ) := H ( t, a ) satisfies Φ( t + T ) = Φ( t ) B, with B :=Φ( T ) . We conclude that ∂∂x x ( kT, a ) = B k . For ε > and ω ∈ R n , let z ( t ) := x ( t, a + εω ) − x ( t, a ) . Then for a fixed time t ≥ , z ( t ) = ε Φ( t ) ω + o ( ε ) , so for afixed k , z ( kT ) = ε Φ( kT ) ω + o ( ε ) = εB k ω + o ( ε ) . Thus, as in the cases described above, B always has as an eigenvalue and a corresponding eigenvector a . Theeigenvalues and eigenvectors of the TP matrix B determinethe response to a small perturbation a → a + εω .The special sign pattern of the eigenvectors v and v n of B can be explained intuitively using the cooperativeand tridiagonal structure of the Jacobian of a TPDS. Wedemonstrate this using a model from systems biology calledthe ribosome flow model (RFM). C. Application: the RFM
The RFM is a phenomenological model for the flow of“material” along a 1D chain composed of n consecutivesites. The RFM is the dynamic mean field approximation ofa fundamental model from statistical mechanics called the totally asymmetric simple exclusion process (TASEP) [2],[11]. This model describes the stochastic motion of particleshopping along a 1D chain of sites while restricted by thesimple exclusion principle, namely, two particles cannotoccupy the same site at the same time.The RFM, and more generally networks of intercon-nected RFMs, has been extensively used to model the flowof ribosomes along the mRNA during translation [19], [14],[13], [20], [34], [21], [17], [31], [33], [15], and also otherimportant processes such as the transfer of the phosphorylgroup from the sensor kinases to the ultimate target duringphosphorelay [1].The RFM is a set of n ODEs: ˙ x i ( t ) = λ i − x i − ( t )(1 − x i ( t )) − λ i x i ( t )(1 − x i +1 ( t )) , (17)ith i = 1 , . . . , n , x ( t ) ≡ , and x n +1 ( t ) ≡ . Thestate-variable x i ( t ) represents the density at site i at time t ,normalized such that x i ( t ) = 0 [ x i ( t ) = 1 ] corresponds tosite i being empty [full] at time t . The positive parameter λ i represents the transition rate from site i to site i + 1 . Inparticular, λ [ λ n ] is called the initiation [exit] rate, and theother λ i s are the elongation rates.To explain (17), consider the equation for ˙ x , namely, ˙ x = λ x (1 − x ) − λ x (1 − x ) . This asserts that the change in the density at site is theflow from site to site minus the flow from site tosite . The flow from site to site , given by λ x (1 − x ) ,is proportional to the density at site , the amount of “freespace” (1 − x ) at site , and the transition rate λ from site to site . In particular, as the density in site increases, theflow from site to site decreases. This is a “soft version”of the simple exclusion principle in TASEP.An important property of the RFM, “inherited”from TASEP, is that allows to study the generation of“traffic jams” along the chain. Indeed, if some λ i is verysmall w.r.t. the other rates then the exit rate from site i will be small. Then the density at site i will increase, andconsequently, the transition rate from site i − will decrease.In this way, a traffic jam of high density sites will evolve“behind” site i . The dynamic evolution and implications oftraffic jams of “biological particles” like ants, ribosomes,and molecular motors are attracting considerable interest(see, e.g. [22], [5], [24], [7]).The output flow from the last site is R ( t ) := λ n x n ( t ) .This represents the rate at which ribosomes exit the mRNA,i.e. the protein production rate.The state space of RFM is the unit cube [0 , n . For a ∈ [0 , n , let x ( t, a ) denote the solution of the RFM at time t with x (0) = a . The flow possesses a unique globallystable (GAS) equilibrium e ∈ (0 , n [13], [16], that is, lim t →∞ x ( t, a ) = e , for all a ∈ [0 , n . This steady-staterepresents a set of densities e , . . . , e n for which the flowinto each site is equal to the flow out of this site. In physics,this is sometimes referred to as a nonequilibrium stationarystate.Let R := λ n e n . This is the flow of ribosomes out of thechain, and thus the protein production rate, at equilibrium.It follows from (17) that λ i e i (1 − e i +1 ) = R, i = 0 , . . . , n, (18)where e := 1 and e n +1 := 0 .Let J ( x ) denote the Jacobian of the vector field inthe RFM. Then J ( x ) admits the tridiagonal structure J ( x ) = T ( a i ( x ) , b i ( x ) , c i ( x )) in (4) with a i ( x ) = − λ i − x i − − λ i (1 − x i +1 ) , b i ( x ) = λ i x i , and c i ( x ) = λ i (1 − x i +1 ) ,where x := 1 , and x n +1 := 0 . Thus, the RFM is a TPDSon the invariant set (0 , n , and Theorem 5 implies that theeigenvalues of J ( e ) are real and simple and the correspond-ing eigenvectors satisfy the sign pattern (13). Furthermore,since the RFM is contractive [16], the eigenvalues of J ( e ) are all real, simple, and negative.The sign pattern of the eigenvectors corresponding to theslowest [fastest] eigenvalue α ( e ) [ α n ( e ) ] can be explainedas follows. Recall that v ( e ) has zero sign variations. Wemay assume that v i ( e ) > for all i . Then e → e + εv ( e ) isdifficult to compensate because it increases all the densitiesin the chain. This corresponds a to very congested situation,and returning to the steady-state e takes more time.On the other-hand, a perturbation in the form e → e + εv n ( e ) is easier to compensate because v n ( e ) has analternating sign pattern, that is, we may assume that v ni ( e ) > for i odd, and v ni ( e ) < for i even. This represents apositive addition to e , a negative addition to e , a positiveaddition to e , and so on. The RFM dynamics compensatesfor this variation quite naturally: the additional density in anodd site generates an increased flow to the consecutive siteand in this way also “automatically” takes care of the lowerdensities in the even sites.The next example considers specific values for the ratesfor which e (and thus J ( e ) ) is known explicitly for any n . Example 4:
Consider the RFM with λ = λ n = 1 / ,and λ i = 1 for i = 1 , . . . , n − . It can be shown using thespectral representation of e [19] that in this case e i = 1 / , i = 1 , . . . , n . This implies that J ( e ) = T ( a i ( e ) , b i ( e ) , c i ( e )) is an n × n Toeplitz tridiagonal matrix with − on the maindiagonal, and / on the sub- and super-diagonals. It is well-known [29] that the eigenvalues of such a matrix are α k = − kπ/ ( n + 1)) , k = 1 , . . . , n, and the eigenvector corresponding to α k is v k = (cid:2) sin( kπn +1 ) sin( kπn +1 ) . . . sin( nkπn +1 ) (cid:3) T . In particular, v := (cid:2) sin( πn +1 ) sin( πn +1 ) . . . sin( nπn +1 ) (cid:3) T . and v n := h sin( nπn +1 ) sin( nπn +1 ) . . . sin( n πn +1 ) i T . As expected, all the eigenvalues are real, negative andsimple, all the entries of v are positive, and the entriesof v n have the sign pattern (+ , − , + , − , . . . ) . The maximaleigenvalue is α = − π/ ( n + 1)) . Thus, lim n →∞ log( − α ( n ))log( n ) = lim n →∞ log( π n +1) − π n +1) + . . . )log( n )= − . (19)Thus, the relaxation time of the system behaves like n forlarge n .The minimal eigenvalue is α n = − πn/ ( n + 1)) , so lim n →∞ α n = − Thus, the fastest convergence rate to e , corresponding to aperturbation in the form e + εv n , is independent of n forlarge n . (cid:3) Remark 6:
The asymptotic behavior described in (19)seems to hold for other rates as well. For example, Fig. 1
Fig. 1. log( − α ( n )) as a function of log( n ) . depicts log( − α ( n )) as a function of log( n ) for the casewhere λ i = 1 for i ∈ { , . . . , n } . It may be seen that forlarge values of n the graph behaves like a line with slope − .IV. C ONCLUSION
Any bounded solution of a T -periodic TPDS converges toa T -periodic solution of the TPDS. If the T -periodic vectorfield represents a T -periodic excitation then this implies thatthe dynamical system entrains to the excitation. In particular,any bounded solution of a time-invariant TPDS converges toan equilibrium [12].It is important to understand how a perturbation froma periodic solution affects the dynamics and, in particular,what are the convergence or divergence rates associated withdifferent perturbation directions. In general, Floquet theorydoes not provide explicit answers to theses questions. Weused the spectral theory of TP matrices to analyze thisquestion for TPDSs and showed that the “best” and “worst”perturbation directions have special sign patterns that admitan intuitive interpretation.An interesting topic for future research is to use this infor-mation to design appropriate control algorithms for TPDSs.R EFERENCES[1] E. Bar-Shalom, A. Ovseevich, and M. Margaliot, “Ribosome flowmodel with different site sizes,”
SIAM J. Applied Dynamical Systems ,vol. 19, pp. 541–576, 2020.[2] R. A. Blythe and M. R. Evans, “Nonequilibrium steady states ofmatrix-product form: a solver’s guide,”
J. Phys. A: Math. Gen. , vol. 40,no. 46, pp. R333–R441, 2007.[3] C. Chicone,
Ordinary Differential Equations with Applications , 2nd ed.New York: Springer, 2006.[4] L. O. Chua and T. Roska, “Stability of a class of nonreciprocal cellularneural networks,”
IEEE Trans. Circuits and Systems , vol. 37, no. 12,pp. 1520–1527, 1990.[5] A. Diament, A. Feldman, E. Schochet, M. Kupiec, Y. Arava, andT. Tuller, “The extent of ribosome queuing in budding yeast,”
PLOSComputational Biology , vol. 14, pp. 1–21, 2018.[6] P. Donnell, S. A. Baigent, and M. Banaji, “Monotone dynamics oftwo cells dynamically coupled by a voltage-dependent gap junction,”
J. Theoretical Biology , vol. 261, no. 1, pp. 120–125, 2009.[7] A. Dussutour, J.-L. Deneubourg, and V. Fourcassie, “Temporal orga-nization of bi-directional traffic in the ant Lasius niger (L.),”
J. Exp.Biol. , vol. 208, pp. 2903–2912, 2005.[8] S. M. Fallat and C. R. Johnson,
Totally Nonnegative Matrices . Prince-ton, NJ: Princeton University Press, 2011.[9] C. Fang, M. Gyllenberg, and Y. Wang, “Floquet bundles for tridiagonalcompetitive-cooperative systems and the dynamics of time-recurrentsystems,”
SIAM J. Math. Anal. , vol. 45, no. 4, pp. 2477–2498, 2013. [10] F. R. Gantmacher and M. G. Krein,
Oscillation Matrices and Kernelsand Small Vibrations of Mechanical Systems . Providence, RI:American Mathematical Society, 2002, translation based on the 1941Russian original.[11] T. Kriecherbauer and J. Krug, “A pedestrian‘s view on interactingparticle systems, KPZ universality and random matrices,”
Journal ofPhysics A: Mathematical and Theoretical , vol. 43, no. 40, p. 403001,2010.[12] M. Margaliot and E. D. Sontag, “Revisiting totally positive differentialsystems: A tutorial and new results,”
Automatica , vol. 101, pp. 1–14,2019.[13] M. Margaliot and T. Tuller, “Stability analysis of the ribosome flowmodel,”
IEEE/ACM Trans. Comput. Biol. Bioinf. , vol. 9, pp. 1545–1552, 2012.[14] ——, “Ribosome flow model with positive feedback,”
J. Royal SocietyInterface , vol. 10, p. 20130267, 2013.[15] M. Margaliot, W. Huleihel, and T. Tuller, “Variability in mRNAtranslation: A random matrix theory approach,”
Scientific Reports ,2020, to appear.[16] M. Margaliot, T. Tuller, and E. D. Sontag, “Checkable conditions forcontraction after small transients in time and amplitude,” in
FeedbackStabilization of Controlled Dynamical Systems: In Honor of LaurentPraly , N. Petit, Ed. Cham, Switzerland: Springer InternationalPublishing, 2017, pp. 279–305.[17] I. Nanikashvili, Y. Zarai, A. Ovseevich, T. Tuller, and M. Margaliot,“Networks of ribosome flow models for modeling and analyzingintracellular traffic,”
Scientific Reports , vol. 9, p. 1703, 2019.[18] A. Pinkus,
Totally Positive Matrices . Cambridge, UK: CambridgeUniversity Press, 2010.[19] G. Poker, Y. Zarai, M. Margaliot, and T. Tuller, “Maximizing proteintranslation rate in the nonhomogeneous ribosome flow model: a convexoptimization approach,”
J. Royal Society Interface , vol. 11, no. 100,2014.[20] A. Raveh, Y. Zarai, M. Margaliot, and T. Tuller, “Ribosome flow modelon a ring,”
IEEE/ACM Trans. Comput. Biol. Bioinf. , vol. 12, no. 6, pp.1429–1439, 2015.[21] A. Raveh, M. Margaliot, E. D. Sontag, and T. Tuller, “A model forcompetition for ribosomes in the cell,”
J. Royal Society Interface ,vol. 13, no. 116, 2016.[22] J. L. Ross, “The impacts of molecular motor traffic jams,”
Proceedingsof the National Academy of Sciences , vol. 109, no. 16, pp. 5911–5912,2012.[23] B. Schwarz, “Totally positive differential systems,”
Pacific J. Math. ,vol. 32, no. 1, pp. 203–229, 1970.[24] S. A. Small, S. Simoes-Spassov, R. Mayeux, and G. A. Petsko,“Endosomal traffic jams represent a pathogenic hub and therapeutictarget in Alzheimer’s disease,”
Trends Neurosci. , vol. 40, pp. 592–602,2017.[25] J. Smillie, “Competitive and cooperative tridiagonal systems of differ-ential equations,”
SIAM J. Math. Anal. , vol. 15, pp. 530–534, 1984.[26] H. L. Smith, “Periodic tridiagonal competitive and cooperative systemsof differential equations,”
SIAM J. Math. Anal. , vol. 22, no. 4, pp.1102–1109, 1991.[27] ——,
Monotone Dynamical Systems: An Introduction to the Theoryof Competitive and Cooperative Systems , ser. Mathematical Surveysand Monographs. Providence, RI: Amer. Math. Soc., 1995, vol. 41.[28] E. Weiss and M. Margaliot, “A generalization of Smillie’s theorem onstrongly cooperative tridiagonal systems,” in
Proc. 57th IEEE Conf.on Decision and Control , Miami Beach, FL, 2018, to appear.[29] W.-C. Yueh, “Eigenvalues of several tridiagonal matrices,”
AppliedMathematics E-Notes , vol. 5, pp. 66–74, 2005.[30] Y. Zarai and M. Margaliot, “On the exponent of several classes ofoscillatory matrices,”
Linear Algebra Appl. , vol. 608, pp. 363–386,2021.[31] Y. Zarai, M. Margaliot, and T. Tuller, “Optimal down regulation ofmRNA translation,”
Scientific Reports , vol. 7, p. 41243, 2017.[32] ——, “Modeling and analyzing the flow of molecular machines ingene expression,” in
Systems Biology , N. Rajewsky, S. Jurga, andJ. Barciszewski, Eds. Springer, 2018, pp. 275–300.[33] ——, “On the ribosomal density that maximizes protein translationrate,”
PLOS ONE , vol. 11, no. 11, pp. 1–26, 2016. [Online].Available: https://doi.org/10.1371/journal.pone.0166481[34] ——, “The ribosome flow model with extended objects,”